1 // Boost.Geometry (aka GGL, Generic Geometry Library)
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.
7 // This file was modified by Oracle on 2014, 2015, 2018, 2019.
8 // Modifications copyright (c) 2014-2019, Oracle and/or its affiliates.
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 // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program
14 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
15 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
17 // Use, modification and distribution is subject to the Boost Software License,
18 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
19 // http://www.boost.org/LICENSE_1_0.txt)
21 #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
22 #define BOOST_GEOMETRY_UTIL_MATH_HPP
27 #include <boost/core/ignore_unused.hpp>
29 #include <boost/math/constants/constants.hpp>
30 #include <boost/math/special_functions/fpclassify.hpp>
31 //#include <boost/math/special_functions/round.hpp>
32 #include <boost/numeric/conversion/cast.hpp>
33 #include <boost/type_traits/is_fundamental.hpp>
34 #include <boost/type_traits/is_integral.hpp>
36 #include <boost/geometry/core/cs.hpp>
38 #include <boost/geometry/util/select_most_precise.hpp>
40 namespace boost { namespace geometry
46 #ifndef DOXYGEN_NO_DETAIL
51 inline T const& greatest(T const& v1, T const& v2)
53 return (std::max)(v1, v2);
57 inline T const& greatest(T const& v1, T const& v2, T const& v3)
59 return (std::max)(greatest(v1, v2), v3);
63 inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4)
65 return (std::max)(greatest(v1, v2, v3), v4);
69 inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5)
71 return (std::max)(greatest(v1, v2, v3, v4), v5);
76 inline T bounded(T const& v, T const& lower, T const& upper)
78 return (std::min)((std::max)(v, lower), upper);
82 inline T bounded(T const& v, T const& lower)
84 return (std::max)(v, lower);
89 bool IsFloatingPoint = boost::is_floating_point<T>::value>
92 static inline T apply(T const& value)
95 return value < zero ? -value : value;
102 static inline T apply(T const& value)
105 using std::fabs; // for long double
112 struct equals_default_policy
114 template <typename T>
115 static inline T apply(T const& a, T const& b)
117 // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17
118 return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1));
122 template <typename T,
123 bool IsFloatingPoint = boost::is_floating_point<T>::value>
124 struct equals_factor_policy
126 equals_factor_policy()
128 explicit equals_factor_policy(T const& v)
129 : factor(greatest(abs<T>::apply(v), T(1)))
131 equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3)
132 : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1),
133 abs<T>::apply(v2), abs<T>::apply(v3),
137 T const& apply(T const&, T const&) const
145 template <typename T>
146 struct equals_factor_policy<T, false>
148 equals_factor_policy() {}
149 explicit equals_factor_policy(T const&) {}
150 equals_factor_policy(T const& , T const& , T const& , T const& ) {}
152 static inline T apply(T const&, T const&)
158 template <typename Type,
159 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
162 template <typename Policy>
163 static inline bool apply(Type const& a, Type const& b, Policy const&)
169 template <typename Type>
170 struct equals<Type, true>
172 template <typename Policy>
173 static inline bool apply(Type const& a, Type const& b, Policy const& policy)
175 boost::ignore_unused(policy);
182 if (boost::math::isfinite(a) && boost::math::isfinite(b))
184 // If a is INF and b is e.g. 0, the expression below returns true
185 // but the values are obviously not equal, hence the condition
186 return abs<Type>::apply(a - b)
187 <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b);
196 template <typename T1, typename T2, typename Policy>
197 inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy)
199 return detail::equals
201 typename select_most_precise<T1, T2>::type
202 >::apply(a, b, policy);
205 template <typename Type,
206 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
209 static inline bool apply(Type const& a, Type const& b)
215 template <typename Type>
216 struct smaller<Type, true>
218 static inline bool apply(Type const& a, Type const& b)
220 if (!(a < b)) // a >= b
225 return ! equals<Type, true>::apply(b, a, equals_default_policy());
229 template <typename Type,
230 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
231 struct smaller_or_equals
233 static inline bool apply(Type const& a, Type const& b)
239 template <typename Type>
240 struct smaller_or_equals<Type, true>
242 static inline bool apply(Type const& a, Type const& b)
249 return equals<Type, true>::apply(a, b, equals_default_policy());
254 template <typename Type,
255 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
256 struct equals_with_epsilon
257 : public equals<Type, IsFloatingPoint>
263 bool IsFundemantal = boost::is_fundamental<T>::value /* false */
267 typedef T return_type;
269 static inline T apply(T const& value)
271 // for non-fundamental number types assume that sqrt is
273 // 1) at T's scope, or
274 // 2) at global scope, or
275 // 3) in namespace std
283 template <typename FundamentalFP>
284 struct square_root_for_fundamental_fp
286 typedef FundamentalFP return_type;
288 static inline FundamentalFP apply(FundamentalFP const& value)
290 #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
291 // This is a workaround for some 32-bit platforms.
292 // For some of those platforms it has been reported that
293 // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value.
294 // For those platforms we need to define the macro
295 // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument
296 // to std::sqrt is checked appropriately before passed to std::sqrt
297 if (boost::math::isfinite(value))
299 return std::sqrt(value);
301 else if (boost::math::isinf(value) && value < 0)
303 return -std::numeric_limits<FundamentalFP>::quiet_NaN();
307 // for fundamental floating point numbers use std::sqrt
308 return std::sqrt(value);
309 #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
314 struct square_root<float, true>
315 : square_root_for_fundamental_fp<float>
320 struct square_root<double, true>
321 : square_root_for_fundamental_fp<double>
326 struct square_root<long double, true>
327 : square_root_for_fundamental_fp<long double>
331 template <typename T>
332 struct square_root<T, true>
334 typedef double return_type;
336 static inline double apply(T const& value)
338 // for all other fundamental number types use also std::sqrt
340 // Note: in C++98 the only other possibility is double;
341 // in C++11 there are also overloads for integral types;
342 // this specialization works for those as well.
343 return square_root_for_fundamental_fp
346 >::apply(boost::numeric_cast<double>(value));
355 bool IsFundemantal = boost::is_fundamental<T>::value /* false */
359 typedef T return_type;
361 static inline T apply(T const& value1, T const& value2)
363 // for non-fundamental number types assume that a free
364 // function mod() is defined either:
365 // 1) at T's scope, or
366 // 2) at global scope
367 return mod(value1, value2);
373 typename Fundamental,
374 bool IsIntegral = boost::is_integral<Fundamental>::value
376 struct modulo_for_fundamental
378 typedef Fundamental return_type;
380 static inline Fundamental apply(Fundamental const& value1,
381 Fundamental const& value2)
383 return value1 % value2;
387 // specialization for floating-point numbers
388 template <typename Fundamental>
389 struct modulo_for_fundamental<Fundamental, false>
391 typedef Fundamental return_type;
393 static inline Fundamental apply(Fundamental const& value1,
394 Fundamental const& value2)
396 return std::fmod(value1, value2);
400 // specialization for fundamental number type
401 template <typename Fundamental>
402 struct modulo<Fundamental, true>
403 : modulo_for_fundamental<Fundamental>
409 \brief Short constructs to enable partial specialization for PI, 2*PI
410 and PI/2, currently not possible in Math.
412 template <typename T>
415 static inline T apply()
417 // Default calls Boost.Math
418 return boost::math::constants::pi<T>();
422 template <typename T>
425 static inline T apply()
427 // Default calls Boost.Math
428 return boost::math::constants::two_pi<T>();
432 template <typename T>
433 struct define_half_pi
435 static inline T apply()
437 // Default calls Boost.Math
438 return boost::math::constants::half_pi<T>();
442 template <typename T>
443 struct relaxed_epsilon
445 static inline T apply(const T& factor)
447 return factor * std::numeric_limits<T>::epsilon();
451 // This must be consistent with math::equals.
452 // By default math::equals() scales the error by epsilon using the greater of
453 // compared values but here is only one value, though it should work the same way.
454 // (a-a) <= max(a, a) * EPS -> 0 <= a*EPS
455 // (a+da-a) <= max(a+da, a) * EPS -> da <= (a+da)*EPS
456 template <typename T, bool IsIntegral = boost::is_integral<T>::value>
457 struct scaled_epsilon
459 static inline T apply(T const& val)
461 return (std::max)(abs<T>::apply(val), T(1))
462 * std::numeric_limits<T>::epsilon();
465 static inline T apply(T const& val, T const& eps)
467 return (std::max)(abs<T>::apply(val), T(1))
472 template <typename T>
473 struct scaled_epsilon<T, true>
475 static inline T apply(T const&)
480 static inline T apply(T const&, T const&)
487 template <typename Result, typename Source,
488 bool ResultIsInteger = std::numeric_limits<Result>::is_integer,
489 bool SourceIsInteger = std::numeric_limits<Source>::is_integer>
492 static inline Result apply(Source const& v)
494 return boost::numeric_cast<Result>(v);
499 template <typename Source, bool ResultIsInteger, bool SourceIsInteger>
500 struct rounding_cast<Source, Source, ResultIsInteger, SourceIsInteger>
502 static inline Source apply(Source const& v)
509 template <typename Result, typename Source>
510 struct rounding_cast<Result, Source, true, false>
512 static inline Result apply(Source const& v)
514 return boost::numeric_cast<Result>(v < Source(0) ?
520 } // namespace detail
524 template <typename T>
525 inline T pi() { return detail::define_pi<T>::apply(); }
527 template <typename T>
528 inline T two_pi() { return detail::define_two_pi<T>::apply(); }
530 template <typename T>
531 inline T half_pi() { return detail::define_half_pi<T>::apply(); }
533 template <typename T>
534 inline T relaxed_epsilon(T const& factor)
536 return detail::relaxed_epsilon<T>::apply(factor);
539 template <typename T>
540 inline T scaled_epsilon(T const& value)
542 return detail::scaled_epsilon<T>::apply(value);
545 template <typename T>
546 inline T scaled_epsilon(T const& value, T const& eps)
548 return detail::scaled_epsilon<T>::apply(value, eps);
551 // Maybe replace this by boost equals or so
554 \brief returns true if both arguments are equal.
556 \param a first argument
557 \param b second argument
558 \return true if a == b
559 \note If both a and b are of an integral type, comparison is done by ==.
560 If one of the types is floating point, comparison is done by abs and
561 comparing with epsilon. If one of the types is non-fundamental, it might
562 be a high-precision number and comparison is done using the == operator
566 template <typename T1, typename T2>
567 inline bool equals(T1 const& a, T2 const& b)
569 return detail::equals
571 typename select_most_precise<T1, T2>::type
572 >::apply(a, b, detail::equals_default_policy());
575 template <typename T1, typename T2>
576 inline bool equals_with_epsilon(T1 const& a, T2 const& b)
578 return detail::equals_with_epsilon
580 typename select_most_precise<T1, T2>::type
581 >::apply(a, b, detail::equals_default_policy());
584 template <typename T1, typename T2>
585 inline bool smaller(T1 const& a, T2 const& b)
587 return detail::smaller
589 typename select_most_precise<T1, T2>::type
593 template <typename T1, typename T2>
594 inline bool larger(T1 const& a, T2 const& b)
596 return detail::smaller
598 typename select_most_precise<T1, T2>::type
602 template <typename T1, typename T2>
603 inline bool smaller_or_equals(T1 const& a, T2 const& b)
605 return detail::smaller_or_equals
607 typename select_most_precise<T1, T2>::type
611 template <typename T1, typename T2>
612 inline bool larger_or_equals(T1 const& a, T2 const& b)
614 return detail::smaller_or_equals
616 typename select_most_precise<T1, T2>::type
621 template <typename T>
624 static T const conversion_coefficient = geometry::math::pi<T>() / T(180.0);
625 return conversion_coefficient;
628 template <typename T>
631 static T const conversion_coefficient = T(180.0) / geometry::math::pi<T>();
632 return conversion_coefficient;
636 #ifndef DOXYGEN_NO_DETAIL
639 template <typename DegreeOrRadian>
642 template <typename T>
643 static inline T apply(T const& value)
650 struct as_radian<degree>
652 template <typename T>
653 static inline T apply(T const& value)
655 return value * d2r<T>();
659 template <typename DegreeOrRadian>
662 template <typename T>
663 static inline T apply(T const& value)
670 struct from_radian<degree>
672 template <typename T>
673 static inline T apply(T const& value)
675 return value * r2d<T>();
679 } // namespace detail
682 template <typename DegreeOrRadian, typename T>
683 inline T as_radian(T const& value)
685 return detail::as_radian<DegreeOrRadian>::apply(value);
688 template <typename DegreeOrRadian, typename T>
689 inline T from_radian(T const& value)
691 return detail::from_radian<DegreeOrRadian>::apply(value);
696 \brief Calculates the haversine of an angle
698 \note See http://en.wikipedia.org/wiki/Haversine_formula
699 haversin(alpha) = sin2(alpha/2)
701 template <typename T>
702 inline T hav(T const& theta)
704 T const half = T(0.5);
705 T const sn = sin(half * theta);
710 \brief Short utility to return the square
712 \param value Value to calculate the square from
713 \return The squared value
715 template <typename T>
716 inline T sqr(T const& value)
718 return value * value;
722 \brief Short utility to return the square root
724 \param value Value to calculate the square root from
725 \return The square root value
727 template <typename T>
728 inline typename detail::square_root<T>::return_type
731 return detail::square_root
733 T, boost::is_fundamental<T>::value
738 \brief Short utility to return the modulo of two values
740 \param value1 First value
741 \param value2 Second value
742 \return The result of the modulo operation on the (ordered) pair
745 template <typename T>
746 inline typename detail::modulo<T>::return_type
747 mod(T const& value1, T const& value2)
749 return detail::modulo
751 T, boost::is_fundamental<T>::value
752 >::apply(value1, value2);
756 \brief Short utility to workaround gcc/clang problem that abs is converting to integer
757 and that older versions of MSVC does not support abs of long long...
761 inline T abs(T const& value)
763 return detail::abs<T>::apply(value);
767 \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
770 template <typename T>
771 inline int sign(T const& value)
774 return value > zero ? 1 : value < zero ? -1 : 0;
778 \brief Short utility to cast a value possibly rounding it to the nearest
781 \note If the source T is NOT an integral type and Result is an integral type
782 the value is rounded towards the closest integral value. Otherwise it's
783 casted without rounding.
785 template <typename Result, typename T>
786 inline Result rounding_cast(T const& v)
788 return detail::rounding_cast<Result, T>::apply(v);
792 \brief Evaluate the sine and cosine function with the argument in degrees
793 \note The results obey exactly the elementary properties of the trigonometric
794 functions, e.g., sin 9° = cos 81° = − sin 123456789°.
795 If x = −0, then \e sinx = −0; this is the only case where
796 −0 is returned.
799 inline void sin_cos_degrees(T const& x,
803 // In order to minimize round-off errors, this function exactly reduces
804 // the argument to the range [-45, 45] before converting it to radians.
805 T remainder; int quotient;
807 remainder = math::mod(x, T(360));
808 quotient = floor(remainder / 90 + T(0.5));
809 remainder -= 90 * quotient;
811 // Convert to radians.
812 remainder *= d2r<T>();
814 T s = sin(remainder), c = cos(remainder);
816 switch (unsigned(quotient) & 3U)
818 case 0U: sinx = s; cosx = c; break;
819 case 1U: sinx = c; cosx = -s; break;
820 case 2U: sinx = -s; cosx = -c; break;
821 default: sinx = -c; cosx = s; break; // case 3U
824 // Set sign of 0 results. -0 only produced for sin(-0).
827 sinx += T(0); cosx += T(0);
832 \brief Round off a given angle
835 inline T round_angle(T const& x) {
836 static const T z = 1/T(16);
845 // z - (z - y) must not be simplified to y.
846 y = y < z ? z - (z - y) : y;
848 return x < 0 ? -y : y;
853 \brief The error-free sum of two numbers.
856 inline T sum_error(T const& u, T const& v, T& t)
858 volatile T s = u + v;
859 volatile T up = s - v;
860 volatile T vpp = s - up;
870 \brief Evaluate the polynomial in x using Horner's method.
872 // TODO: adl1995 - Merge these functions with formulas/area_formulas.hpp
873 // i.e. place them in one file.
874 template <typename NT, typename IteratorType>
875 inline NT horner_evaluate(NT const& x,
880 IteratorType it = end;
883 result = result * x + *--it;
890 \brief Evaluate the polynomial.
892 template<typename IteratorType, typename CT>
893 inline CT polyval(IteratorType first,
897 int N = std::distance(first, last) - 1;
900 CT y = N < 0 ? 0 : *(first + (index++));
904 y = y * eps + *(first + (index++));
911 \brief Short utility to calculate the power
914 template <typename T1, typename T2>
915 inline T1 pow(T1 const& a, T2 const& b)
924 }} // namespace boost::geometry
926 #endif // BOOST_GEOMETRY_UTIL_MATH_HPP