1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 
3 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
4 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
5 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
6 
7 // This file was modified by Oracle on 2014, 2015.
8 // Modifications copyright (c) 2014-2015, Oracle and/or its affiliates.
9 
10 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
12 
13 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
14 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
15 
16 // Use, modification and distribution is subject to the Boost Software License,
17 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
18 // http://www.boost.org/LICENSE_1_0.txt)
19 
20 #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
21 #define BOOST_GEOMETRY_UTIL_MATH_HPP
22 
23 #include <cmath>
24 #include <limits>
25 
26 #include <boost/core/ignore_unused.hpp>
27 
28 #include <boost/math/constants/constants.hpp>
29 #include <boost/math/special_functions/fpclassify.hpp>
30 //#include <boost/math/special_functions/round.hpp>
31 #include <boost/numeric/conversion/cast.hpp>
32 #include <boost/type_traits/is_fundamental.hpp>
33 #include <boost/type_traits/is_integral.hpp>
34 
35 #include <boost/geometry/core/cs.hpp>
36 
37 #include <boost/geometry/util/select_most_precise.hpp>
38 
39 namespace boost { namespace geometry
40 {
41 
42 namespace math
43 {
44 
45 #ifndef DOXYGEN_NO_DETAIL
46 namespace detail
47 {
48 
49 template <typename T>
greatest(T const & v1,T const & v2)50 inline T const& greatest(T const& v1, T const& v2)
51 {
52     return (std::max)(v1, v2);
53 }
54 
55 template <typename T>
greatest(T const & v1,T const & v2,T const & v3)56 inline T const& greatest(T const& v1, T const& v2, T const& v3)
57 {
58     return (std::max)(greatest(v1, v2), v3);
59 }
60 
61 template <typename T>
greatest(T const & v1,T const & v2,T const & v3,T const & v4)62 inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4)
63 {
64     return (std::max)(greatest(v1, v2, v3), v4);
65 }
66 
67 template <typename T>
greatest(T const & v1,T const & v2,T const & v3,T const & v4,T const & v5)68 inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5)
69 {
70     return (std::max)(greatest(v1, v2, v3, v4), v5);
71 }
72 
73 
74 template <typename T>
bounded(T const & v,T const & lower,T const & upper)75 inline T bounded(T const& v, T const& lower, T const& upper)
76 {
77     return (std::min)((std::max)(v, lower), upper);
78 }
79 
80 template <typename T>
bounded(T const & v,T const & lower)81 inline T bounded(T const& v, T const& lower)
82 {
83     return (std::max)(v, lower);
84 }
85 
86 
87 template <typename T,
88           bool IsFloatingPoint = boost::is_floating_point<T>::value>
89 struct abs
90 {
applyboost::geometry::math::detail::abs91     static inline T apply(T const& value)
92     {
93         T const zero = T();
94         return value < zero ? -value : value;
95     }
96 };
97 
98 template <typename T>
99 struct abs<T, true>
100 {
applyboost::geometry::math::detail::abs101     static inline T apply(T const& value)
102     {
103         using ::fabs;
104         using std::fabs; // for long double
105 
106         return fabs(value);
107     }
108 };
109 
110 
111 struct equals_default_policy
112 {
113     template <typename T>
applyboost::geometry::math::detail::equals_default_policy114     static inline T apply(T const& a, T const& b)
115     {
116         // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17
117         return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1));
118     }
119 };
120 
121 template <typename T,
122           bool IsFloatingPoint = boost::is_floating_point<T>::value>
123 struct equals_factor_policy
124 {
equals_factor_policyboost::geometry::math::detail::equals_factor_policy125     equals_factor_policy()
126         : factor(1) {}
equals_factor_policyboost::geometry::math::detail::equals_factor_policy127     explicit equals_factor_policy(T const& v)
128         : factor(greatest(abs<T>::apply(v), T(1)))
129     {}
equals_factor_policyboost::geometry::math::detail::equals_factor_policy130     equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3)
131         : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1),
132                           abs<T>::apply(v2), abs<T>::apply(v3),
133                           T(1)))
134     {}
135 
applyboost::geometry::math::detail::equals_factor_policy136     T const& apply(T const&, T const&) const
137     {
138         return factor;
139     }
140 
141     T factor;
142 };
143 
144 template <typename T>
145 struct equals_factor_policy<T, false>
146 {
equals_factor_policyboost::geometry::math::detail::equals_factor_policy147     equals_factor_policy() {}
equals_factor_policyboost::geometry::math::detail::equals_factor_policy148     explicit equals_factor_policy(T const&) {}
equals_factor_policyboost::geometry::math::detail::equals_factor_policy149     equals_factor_policy(T const& , T const& , T const& , T const& ) {}
150 
applyboost::geometry::math::detail::equals_factor_policy151     static inline T apply(T const&, T const&)
152     {
153         return T(1);
154     }
155 };
156 
157 template <typename Type,
158           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
159 struct equals
160 {
161     template <typename Policy>
applyboost::geometry::math::detail::equals162     static inline bool apply(Type const& a, Type const& b, Policy const&)
163     {
164         return a == b;
165     }
166 };
167 
168 template <typename Type>
169 struct equals<Type, true>
170 {
171     template <typename Policy>
applyboost::geometry::math::detail::equals172     static inline bool apply(Type const& a, Type const& b, Policy const& policy)
173     {
174         boost::ignore_unused(policy);
175 
176         if (a == b)
177         {
178             return true;
179         }
180 
181         if (boost::math::isfinite(a) && boost::math::isfinite(b))
182         {
183             // If a is INF and b is e.g. 0, the expression below returns true
184             // but the values are obviously not equal, hence the condition
185             return abs<Type>::apply(a - b)
186                 <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b);
187         }
188         else
189         {
190             return a == b;
191         }
192     }
193 };
194 
195 template <typename T1, typename T2, typename Policy>
equals_by_policy(T1 const & a,T2 const & b,Policy const & policy)196 inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy)
197 {
198     return detail::equals
199         <
200             typename select_most_precise<T1, T2>::type
201         >::apply(a, b, policy);
202 }
203 
204 template <typename Type,
205           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
206 struct smaller
207 {
applyboost::geometry::math::detail::smaller208     static inline bool apply(Type const& a, Type const& b)
209     {
210         return a < b;
211     }
212 };
213 
214 template <typename Type>
215 struct smaller<Type, true>
216 {
applyboost::geometry::math::detail::smaller217     static inline bool apply(Type const& a, Type const& b)
218     {
219         if (!(a < b)) // a >= b
220         {
221             return false;
222         }
223 
224         return ! equals<Type, true>::apply(b, a, equals_default_policy());
225     }
226 };
227 
228 template <typename Type,
229           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
230 struct smaller_or_equals
231 {
applyboost::geometry::math::detail::smaller_or_equals232     static inline bool apply(Type const& a, Type const& b)
233     {
234         return a <= b;
235     }
236 };
237 
238 template <typename Type>
239 struct smaller_or_equals<Type, true>
240 {
applyboost::geometry::math::detail::smaller_or_equals241     static inline bool apply(Type const& a, Type const& b)
242     {
243         if (a <= b)
244         {
245             return true;
246         }
247 
248         return equals<Type, true>::apply(a, b, equals_default_policy());
249     }
250 };
251 
252 
253 template <typename Type,
254           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
255 struct equals_with_epsilon
256     : public equals<Type, IsFloatingPoint>
257 {};
258 
259 template
260 <
261     typename T,
262     bool IsFundemantal = boost::is_fundamental<T>::value /* false */
263 >
264 struct square_root
265 {
266     typedef T return_type;
267 
applyboost::geometry::math::detail::square_root268     static inline T apply(T const& value)
269     {
270         // for non-fundamental number types assume that sqrt is
271         // defined either:
272         // 1) at T's scope, or
273         // 2) at global scope, or
274         // 3) in namespace std
275         using ::sqrt;
276         using std::sqrt;
277 
278         return sqrt(value);
279     }
280 };
281 
282 template <typename FundamentalFP>
283 struct square_root_for_fundamental_fp
284 {
285     typedef FundamentalFP return_type;
286 
applyboost::geometry::math::detail::square_root_for_fundamental_fp287     static inline FundamentalFP apply(FundamentalFP const& value)
288     {
289 #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
290         // This is a workaround for some 32-bit platforms.
291         // For some of those platforms it has been reported that
292         // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value.
293         // For those platforms we need to define the macro
294         // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument
295         // to std::sqrt is checked appropriately before passed to std::sqrt
296         if (boost::math::isfinite(value))
297         {
298             return std::sqrt(value);
299         }
300         else if (boost::math::isinf(value) && value < 0)
301         {
302             return -std::numeric_limits<FundamentalFP>::quiet_NaN();
303         }
304         return value;
305 #else
306         // for fundamental floating point numbers use std::sqrt
307         return std::sqrt(value);
308 #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
309     }
310 };
311 
312 template <>
313 struct square_root<float, true>
314     : square_root_for_fundamental_fp<float>
315 {
316 };
317 
318 template <>
319 struct square_root<double, true>
320     : square_root_for_fundamental_fp<double>
321 {
322 };
323 
324 template <>
325 struct square_root<long double, true>
326     : square_root_for_fundamental_fp<long double>
327 {
328 };
329 
330 template <typename T>
331 struct square_root<T, true>
332 {
333     typedef double return_type;
334 
applyboost::geometry::math::detail::square_root335     static inline double apply(T const& value)
336     {
337         // for all other fundamental number types use also std::sqrt
338         //
339         // Note: in C++98 the only other possibility is double;
340         //       in C++11 there are also overloads for integral types;
341         //       this specialization works for those as well.
342         return square_root_for_fundamental_fp
343             <
344                 double
345             >::apply(boost::numeric_cast<double>(value));
346     }
347 };
348 
349 
350 
351 template
352 <
353     typename T,
354     bool IsFundemantal = boost::is_fundamental<T>::value /* false */
355 >
356 struct modulo
357 {
358     typedef T return_type;
359 
applyboost::geometry::math::detail::modulo360     static inline T apply(T const& value1, T const& value2)
361     {
362         // for non-fundamental number types assume that a free
363         // function mod() is defined either:
364         // 1) at T's scope, or
365         // 2) at global scope
366         return mod(value1, value2);
367     }
368 };
369 
370 template
371 <
372     typename Fundamental,
373     bool IsIntegral = boost::is_integral<Fundamental>::value
374 >
375 struct modulo_for_fundamental
376 {
377     typedef Fundamental return_type;
378 
applyboost::geometry::math::detail::modulo_for_fundamental379     static inline Fundamental apply(Fundamental const& value1,
380                                     Fundamental const& value2)
381     {
382         return value1 % value2;
383     }
384 };
385 
386 // specialization for floating-point numbers
387 template <typename Fundamental>
388 struct modulo_for_fundamental<Fundamental, false>
389 {
390     typedef Fundamental return_type;
391 
applyboost::geometry::math::detail::modulo_for_fundamental392     static inline Fundamental apply(Fundamental const& value1,
393                                     Fundamental const& value2)
394     {
395         return std::fmod(value1, value2);
396     }
397 };
398 
399 // specialization for fundamental number type
400 template <typename Fundamental>
401 struct modulo<Fundamental, true>
402     : modulo_for_fundamental<Fundamental>
403 {};
404 
405 
406 
407 /*!
408 \brief Short constructs to enable partial specialization for PI, 2*PI
409        and PI/2, currently not possible in Math.
410 */
411 template <typename T>
412 struct define_pi
413 {
applyboost::geometry::math::detail::define_pi414     static inline T apply()
415     {
416         // Default calls Boost.Math
417         return boost::math::constants::pi<T>();
418     }
419 };
420 
421 template <typename T>
422 struct define_two_pi
423 {
applyboost::geometry::math::detail::define_two_pi424     static inline T apply()
425     {
426         // Default calls Boost.Math
427         return boost::math::constants::two_pi<T>();
428     }
429 };
430 
431 template <typename T>
432 struct define_half_pi
433 {
applyboost::geometry::math::detail::define_half_pi434     static inline T apply()
435     {
436         // Default calls Boost.Math
437         return boost::math::constants::half_pi<T>();
438     }
439 };
440 
441 template <typename T>
442 struct relaxed_epsilon
443 {
applyboost::geometry::math::detail::relaxed_epsilon444     static inline T apply(const T& factor)
445     {
446         return factor * std::numeric_limits<T>::epsilon();
447     }
448 };
449 
450 // This must be consistent with math::equals.
451 // By default math::equals() scales the error by epsilon using the greater of
452 // compared values but here is only one value, though it should work the same way.
453 // (a-a) <= max(a, a) * EPS       -> 0 <= a*EPS
454 // (a+da-a) <= max(a+da, a) * EPS -> da <= (a+da)*EPS
455 template <typename T, bool IsFloat = boost::is_floating_point<T>::value>
456 struct scaled_epsilon
457 {
applyboost::geometry::math::detail::scaled_epsilon458     static inline T apply(T const& val)
459     {
460         return (std::max)(abs<T>::apply(val), T(1))
461                     * std::numeric_limits<T>::epsilon();
462     }
463 };
464 
465 template <typename T>
466 struct scaled_epsilon<T, false>
467 {
applyboost::geometry::math::detail::scaled_epsilon468     static inline T apply(T const&)
469     {
470         return T(0);
471     }
472 };
473 
474 // ItoF ItoI FtoF
475 template <typename Result, typename Source,
476           bool ResultIsInteger = std::numeric_limits<Result>::is_integer,
477           bool SourceIsInteger = std::numeric_limits<Source>::is_integer>
478 struct rounding_cast
479 {
applyboost::geometry::math::detail::rounding_cast480     static inline Result apply(Source const& v)
481     {
482         return boost::numeric_cast<Result>(v);
483     }
484 };
485 
486 // TtoT
487 template <typename Source, bool ResultIsInteger, bool SourceIsInteger>
488 struct rounding_cast<Source, Source, ResultIsInteger, SourceIsInteger>
489 {
applyboost::geometry::math::detail::rounding_cast490     static inline Source apply(Source const& v)
491     {
492         return v;
493     }
494 };
495 
496 // FtoI
497 template <typename Result, typename Source>
498 struct rounding_cast<Result, Source, true, false>
499 {
applyboost::geometry::math::detail::rounding_cast500     static inline Result apply(Source const& v)
501     {
502         return boost::numeric_cast<Result>(v < Source(0) ?
503                                             v - Source(0.5) :
504                                             v + Source(0.5));
505     }
506 };
507 
508 } // namespace detail
509 #endif
510 
511 
512 template <typename T>
pi()513 inline T pi() { return detail::define_pi<T>::apply(); }
514 
515 template <typename T>
two_pi()516 inline T two_pi() { return detail::define_two_pi<T>::apply(); }
517 
518 template <typename T>
half_pi()519 inline T half_pi() { return detail::define_half_pi<T>::apply(); }
520 
521 template <typename T>
relaxed_epsilon(T const & factor)522 inline T relaxed_epsilon(T const& factor)
523 {
524     return detail::relaxed_epsilon<T>::apply(factor);
525 }
526 
527 template <typename T>
scaled_epsilon(T const & value)528 inline T scaled_epsilon(T const& value)
529 {
530     return detail::scaled_epsilon<T>::apply(value);
531 }
532 
533 
534 // Maybe replace this by boost equals or so
535 
536 /*!
537     \brief returns true if both arguments are equal.
538     \ingroup utility
539     \param a first argument
540     \param b second argument
541     \return true if a == b
542     \note If both a and b are of an integral type, comparison is done by ==.
543     If one of the types is floating point, comparison is done by abs and
544     comparing with epsilon. If one of the types is non-fundamental, it might
545     be a high-precision number and comparison is done using the == operator
546     of that class.
547 */
548 
549 template <typename T1, typename T2>
equals(T1 const & a,T2 const & b)550 inline bool equals(T1 const& a, T2 const& b)
551 {
552     return detail::equals
553         <
554             typename select_most_precise<T1, T2>::type
555         >::apply(a, b, detail::equals_default_policy());
556 }
557 
558 template <typename T1, typename T2>
equals_with_epsilon(T1 const & a,T2 const & b)559 inline bool equals_with_epsilon(T1 const& a, T2 const& b)
560 {
561     return detail::equals_with_epsilon
562         <
563             typename select_most_precise<T1, T2>::type
564         >::apply(a, b, detail::equals_default_policy());
565 }
566 
567 template <typename T1, typename T2>
smaller(T1 const & a,T2 const & b)568 inline bool smaller(T1 const& a, T2 const& b)
569 {
570     return detail::smaller
571         <
572             typename select_most_precise<T1, T2>::type
573         >::apply(a, b);
574 }
575 
576 template <typename T1, typename T2>
larger(T1 const & a,T2 const & b)577 inline bool larger(T1 const& a, T2 const& b)
578 {
579     return detail::smaller
580         <
581             typename select_most_precise<T1, T2>::type
582         >::apply(b, a);
583 }
584 
585 template <typename T1, typename T2>
smaller_or_equals(T1 const & a,T2 const & b)586 inline bool smaller_or_equals(T1 const& a, T2 const& b)
587 {
588     return detail::smaller_or_equals
589         <
590             typename select_most_precise<T1, T2>::type
591         >::apply(a, b);
592 }
593 
594 template <typename T1, typename T2>
larger_or_equals(T1 const & a,T2 const & b)595 inline bool larger_or_equals(T1 const& a, T2 const& b)
596 {
597     return detail::smaller_or_equals
598         <
599             typename select_most_precise<T1, T2>::type
600         >::apply(b, a);
601 }
602 
603 
604 template <typename T>
d2r()605 inline T d2r()
606 {
607     static T const conversion_coefficient = geometry::math::pi<T>() / T(180.0);
608     return conversion_coefficient;
609 }
610 
611 template <typename T>
r2d()612 inline T r2d()
613 {
614     static T const conversion_coefficient = T(180.0) / geometry::math::pi<T>();
615     return conversion_coefficient;
616 }
617 
618 
619 #ifndef DOXYGEN_NO_DETAIL
620 namespace detail {
621 
622 template <typename DegreeOrRadian>
623 struct as_radian
624 {
625     template <typename T>
applyboost::geometry::math::detail::as_radian626     static inline T apply(T const& value)
627     {
628         return value;
629     }
630 };
631 
632 template <>
633 struct as_radian<degree>
634 {
635     template <typename T>
applyboost::geometry::math::detail::as_radian636     static inline T apply(T const& value)
637     {
638         return value * d2r<T>();
639     }
640 };
641 
642 template <typename DegreeOrRadian>
643 struct from_radian
644 {
645     template <typename T>
applyboost::geometry::math::detail::from_radian646     static inline T apply(T const& value)
647     {
648         return value;
649     }
650 };
651 
652 template <>
653 struct from_radian<degree>
654 {
655     template <typename T>
applyboost::geometry::math::detail::from_radian656     static inline T apply(T const& value)
657     {
658         return value * r2d<T>();
659     }
660 };
661 
662 } // namespace detail
663 #endif
664 
665 template <typename DegreeOrRadian, typename T>
as_radian(T const & value)666 inline T as_radian(T const& value)
667 {
668     return detail::as_radian<DegreeOrRadian>::apply(value);
669 }
670 
671 template <typename DegreeOrRadian, typename T>
from_radian(T const & value)672 inline T from_radian(T const& value)
673 {
674     return detail::from_radian<DegreeOrRadian>::apply(value);
675 }
676 
677 
678 /*!
679     \brief Calculates the haversine of an angle
680     \ingroup utility
681     \note See http://en.wikipedia.org/wiki/Haversine_formula
682     haversin(alpha) = sin2(alpha/2)
683 */
684 template <typename T>
hav(T const & theta)685 inline T hav(T const& theta)
686 {
687     T const half = T(0.5);
688     T const sn = sin(half * theta);
689     return sn * sn;
690 }
691 
692 /*!
693 \brief Short utility to return the square
694 \ingroup utility
695 \param value Value to calculate the square from
696 \return The squared value
697 */
698 template <typename T>
sqr(T const & value)699 inline T sqr(T const& value)
700 {
701     return value * value;
702 }
703 
704 /*!
705 \brief Short utility to return the square root
706 \ingroup utility
707 \param value Value to calculate the square root from
708 \return The square root value
709 */
710 template <typename T>
711 inline typename detail::square_root<T>::return_type
sqrt(T const & value)712 sqrt(T const& value)
713 {
714     return detail::square_root
715         <
716             T, boost::is_fundamental<T>::value
717         >::apply(value);
718 }
719 
720 /*!
721 \brief Short utility to return the modulo of two values
722 \ingroup utility
723 \param value1 First value
724 \param value2 Second value
725 \return The result of the modulo operation on the (ordered) pair
726 (value1, value2)
727 */
728 template <typename T>
729 inline typename detail::modulo<T>::return_type
mod(T const & value1,T const & value2)730 mod(T const& value1, T const& value2)
731 {
732     return detail::modulo
733         <
734             T, boost::is_fundamental<T>::value
735         >::apply(value1, value2);
736 }
737 
738 /*!
739 \brief Short utility to workaround gcc/clang problem that abs is converting to integer
740        and that older versions of MSVC does not support abs of long long...
741 \ingroup utility
742 */
743 template<typename T>
abs(T const & value)744 inline T abs(T const& value)
745 {
746     return detail::abs<T>::apply(value);
747 }
748 
749 /*!
750 \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
751 \ingroup utility
752 */
753 template <typename T>
sign(T const & value)754 inline int sign(T const& value)
755 {
756     T const zero = T();
757     return value > zero ? 1 : value < zero ? -1 : 0;
758 }
759 
760 /*!
761 \brief Short utility to cast a value possibly rounding it to the nearest
762        integral value.
763 \ingroup utility
764 \note If the source T is NOT an integral type and Result is an integral type
765       the value is rounded towards the closest integral value. Otherwise it's
766       casted without rounding.
767 */
768 template <typename Result, typename T>
rounding_cast(T const & v)769 inline Result rounding_cast(T const& v)
770 {
771     return detail::rounding_cast<Result, T>::apply(v);
772 }
773 
774 } // namespace math
775 
776 
777 }} // namespace boost::geometry
778 
779 #endif // BOOST_GEOMETRY_UTIL_MATH_HPP
780