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.
8 // Modifications copyright (c) 2014-2015, 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
13 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
14 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
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)
20 #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
21 #define BOOST_GEOMETRY_UTIL_MATH_HPP
26 #include <boost/core/ignore_unused.hpp>
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>
35 #include <boost/geometry/core/cs.hpp>
37 #include <boost/geometry/util/select_most_precise.hpp>
39 namespace boost { namespace geometry
45 #ifndef DOXYGEN_NO_DETAIL
50 inline T const& greatest(T const& v1, T const& v2)
52 return (std::max)(v1, v2);
56 inline T const& greatest(T const& v1, T const& v2, T const& v3)
58 return (std::max)(greatest(v1, v2), v3);
62 inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4)
64 return (std::max)(greatest(v1, v2, v3), v4);
68 inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5)
70 return (std::max)(greatest(v1, v2, v3, v4), v5);
75 bool IsFloatingPoint = boost::is_floating_point<T>::value>
78 static inline T apply(T const& value)
81 return value < zero ? -value : value;
88 static inline T apply(T const& value)
91 using std::fabs; // for long double
98 struct equals_default_policy
100 template <typename T>
101 static inline T apply(T const& a, T const& b)
103 // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17
104 return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1));
108 template <typename T,
109 bool IsFloatingPoint = boost::is_floating_point<T>::value>
110 struct equals_factor_policy
112 equals_factor_policy()
114 explicit equals_factor_policy(T const& v)
115 : factor(greatest(abs<T>::apply(v), T(1)))
117 equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3)
118 : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1),
119 abs<T>::apply(v2), abs<T>::apply(v3),
123 T const& apply(T const&, T const&) const
131 template <typename T>
132 struct equals_factor_policy<T, false>
134 equals_factor_policy() {}
135 explicit equals_factor_policy(T const&) {}
136 equals_factor_policy(T const& , T const& , T const& , T const& ) {}
138 static inline T apply(T const&, T const&)
144 template <typename Type,
145 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
148 template <typename Policy>
149 static inline bool apply(Type const& a, Type const& b, Policy const&)
155 template <typename Type>
156 struct equals<Type, true>
158 template <typename Policy>
159 static inline bool apply(Type const& a, Type const& b, Policy const& policy)
161 boost::ignore_unused(policy);
168 if (boost::math::isfinite(a) && boost::math::isfinite(b))
170 // If a is INF and b is e.g. 0, the expression below returns true
171 // but the values are obviously not equal, hence the condition
172 return abs<Type>::apply(a - b)
173 <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b);
182 template <typename T1, typename T2, typename Policy>
183 inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy)
185 return detail::equals
187 typename select_most_precise<T1, T2>::type
188 >::apply(a, b, policy);
191 template <typename Type,
192 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
195 static inline bool apply(Type const& a, Type const& b)
201 template <typename Type>
202 struct smaller<Type, true>
204 static inline bool apply(Type const& a, Type const& b)
206 if (!(a < b)) // a >= b
211 return ! equals<Type, true>::apply(b, a, equals_default_policy());
215 template <typename Type,
216 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
217 struct smaller_or_equals
219 static inline bool apply(Type const& a, Type const& b)
225 template <typename Type>
226 struct smaller_or_equals<Type, true>
228 static inline bool apply(Type const& a, Type const& b)
235 return equals<Type, true>::apply(a, b, equals_default_policy());
240 template <typename Type,
241 bool IsFloatingPoint = boost::is_floating_point<Type>::value>
242 struct equals_with_epsilon
243 : public equals<Type, IsFloatingPoint>
249 bool IsFundemantal = boost::is_fundamental<T>::value /* false */
253 typedef T return_type;
255 static inline T apply(T const& value)
257 // for non-fundamental number types assume that sqrt is
259 // 1) at T's scope, or
260 // 2) at global scope, or
261 // 3) in namespace std
269 template <typename FundamentalFP>
270 struct square_root_for_fundamental_fp
272 typedef FundamentalFP return_type;
274 static inline FundamentalFP apply(FundamentalFP const& value)
276 #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
277 // This is a workaround for some 32-bit platforms.
278 // For some of those platforms it has been reported that
279 // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value.
280 // For those platforms we need to define the macro
281 // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument
282 // to std::sqrt is checked appropriately before passed to std::sqrt
283 if (boost::math::isfinite(value))
285 return std::sqrt(value);
287 else if (boost::math::isinf(value) && value < 0)
289 return -std::numeric_limits<FundamentalFP>::quiet_NaN();
293 // for fundamental floating point numbers use std::sqrt
294 return std::sqrt(value);
295 #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
300 struct square_root<float, true>
301 : square_root_for_fundamental_fp<float>
306 struct square_root<double, true>
307 : square_root_for_fundamental_fp<double>
312 struct square_root<long double, true>
313 : square_root_for_fundamental_fp<long double>
317 template <typename T>
318 struct square_root<T, true>
320 typedef double return_type;
322 static inline double apply(T const& value)
324 // for all other fundamental number types use also std::sqrt
326 // Note: in C++98 the only other possibility is double;
327 // in C++11 there are also overloads for integral types;
328 // this specialization works for those as well.
329 return square_root_for_fundamental_fp
332 >::apply(boost::numeric_cast<double>(value));
341 bool IsFundemantal = boost::is_fundamental<T>::value /* false */
345 typedef T return_type;
347 static inline T apply(T const& value1, T const& value2)
349 // for non-fundamental number types assume that a free
350 // function mod() is defined either:
351 // 1) at T's scope, or
352 // 2) at global scope
353 return mod(value1, value2);
359 typename Fundamental,
360 bool IsIntegral = boost::is_integral<Fundamental>::value
362 struct modulo_for_fundamental
364 typedef Fundamental return_type;
366 static inline Fundamental apply(Fundamental const& value1,
367 Fundamental const& value2)
369 return value1 % value2;
373 // specialization for floating-point numbers
374 template <typename Fundamental>
375 struct modulo_for_fundamental<Fundamental, false>
377 typedef Fundamental return_type;
379 static inline Fundamental apply(Fundamental const& value1,
380 Fundamental const& value2)
382 return std::fmod(value1, value2);
386 // specialization for fundamental number type
387 template <typename Fundamental>
388 struct modulo<Fundamental, true>
389 : modulo_for_fundamental<Fundamental>
395 \brief Short constructs to enable partial specialization for PI, 2*PI
396 and PI/2, currently not possible in Math.
398 template <typename T>
401 static inline T apply()
403 // Default calls Boost.Math
404 return boost::math::constants::pi<T>();
408 template <typename T>
411 static inline T apply()
413 // Default calls Boost.Math
414 return boost::math::constants::two_pi<T>();
418 template <typename T>
419 struct define_half_pi
421 static inline T apply()
423 // Default calls Boost.Math
424 return boost::math::constants::half_pi<T>();
428 template <typename T>
429 struct relaxed_epsilon
431 static inline T apply(const T& factor)
433 return factor * std::numeric_limits<T>::epsilon();
437 // This must be consistent with math::equals.
438 // By default math::equals() scales the error by epsilon using the greater of
439 // compared values but here is only one value, though it should work the same way.
440 // (a-a) <= max(a, a) * EPS -> 0 <= a*EPS
441 // (a+da-a) <= max(a+da, a) * EPS -> da <= (a+da)*EPS
442 template <typename T, bool IsFloat = boost::is_floating_point<T>::value>
443 struct scaled_epsilon
445 static inline T apply(T const& val)
447 return (std::max)(abs<T>::apply(val), T(1))
448 * std::numeric_limits<T>::epsilon();
452 template <typename T>
453 struct scaled_epsilon<T, false>
455 static inline T apply(T const&)
462 template <typename Result, typename Source,
463 bool ResultIsInteger = std::numeric_limits<Result>::is_integer,
464 bool SourceIsInteger = std::numeric_limits<Source>::is_integer>
467 static inline Result apply(Source const& v)
469 return boost::numeric_cast<Result>(v);
474 template <typename Source, bool ResultIsInteger, bool SourceIsInteger>
475 struct rounding_cast<Source, Source, ResultIsInteger, SourceIsInteger>
477 static inline Source apply(Source const& v)
484 template <typename Result, typename Source>
485 struct rounding_cast<Result, Source, true, false>
487 static inline Result apply(Source const& v)
489 return boost::numeric_cast<Result>(v < Source(0) ?
495 } // namespace detail
499 template <typename T>
500 inline T pi() { return detail::define_pi<T>::apply(); }
502 template <typename T>
503 inline T two_pi() { return detail::define_two_pi<T>::apply(); }
505 template <typename T>
506 inline T half_pi() { return detail::define_half_pi<T>::apply(); }
508 template <typename T>
509 inline T relaxed_epsilon(T const& factor)
511 return detail::relaxed_epsilon<T>::apply(factor);
514 template <typename T>
515 inline T scaled_epsilon(T const& value)
517 return detail::scaled_epsilon<T>::apply(value);
521 // Maybe replace this by boost equals or boost ublas numeric equals or so
524 \brief returns true if both arguments are equal.
526 \param a first argument
527 \param b second argument
528 \return true if a == b
529 \note If both a and b are of an integral type, comparison is done by ==.
530 If one of the types is floating point, comparison is done by abs and
531 comparing with epsilon. If one of the types is non-fundamental, it might
532 be a high-precision number and comparison is done using the == operator
536 template <typename T1, typename T2>
537 inline bool equals(T1 const& a, T2 const& b)
539 return detail::equals
541 typename select_most_precise<T1, T2>::type
542 >::apply(a, b, detail::equals_default_policy());
545 template <typename T1, typename T2>
546 inline bool equals_with_epsilon(T1 const& a, T2 const& b)
548 return detail::equals_with_epsilon
550 typename select_most_precise<T1, T2>::type
551 >::apply(a, b, detail::equals_default_policy());
554 template <typename T1, typename T2>
555 inline bool smaller(T1 const& a, T2 const& b)
557 return detail::smaller
559 typename select_most_precise<T1, T2>::type
563 template <typename T1, typename T2>
564 inline bool larger(T1 const& a, T2 const& b)
566 return detail::smaller
568 typename select_most_precise<T1, T2>::type
572 template <typename T1, typename T2>
573 inline bool smaller_or_equals(T1 const& a, T2 const& b)
575 return detail::smaller_or_equals
577 typename select_most_precise<T1, T2>::type
581 template <typename T1, typename T2>
582 inline bool larger_or_equals(T1 const& a, T2 const& b)
584 return detail::smaller_or_equals
586 typename select_most_precise<T1, T2>::type
591 template <typename T>
594 static T const conversion_coefficient = geometry::math::pi<T>() / T(180.0);
595 return conversion_coefficient;
598 template <typename T>
601 static T const conversion_coefficient = T(180.0) / geometry::math::pi<T>();
602 return conversion_coefficient;
606 #ifndef DOXYGEN_NO_DETAIL
609 template <typename DegreeOrRadian>
612 template <typename T>
613 static inline T apply(T const& value)
620 struct as_radian<degree>
622 template <typename T>
623 static inline T apply(T const& value)
625 return value * d2r<T>();
629 template <typename DegreeOrRadian>
632 template <typename T>
633 static inline T apply(T const& value)
640 struct from_radian<degree>
642 template <typename T>
643 static inline T apply(T const& value)
645 return value * r2d<T>();
649 } // namespace detail
652 template <typename DegreeOrRadian, typename T>
653 inline T as_radian(T const& value)
655 return detail::as_radian<DegreeOrRadian>::apply(value);
658 template <typename DegreeOrRadian, typename T>
659 inline T from_radian(T const& value)
661 return detail::from_radian<DegreeOrRadian>::apply(value);
666 \brief Calculates the haversine of an angle
668 \note See http://en.wikipedia.org/wiki/Haversine_formula
669 haversin(alpha) = sin2(alpha/2)
671 template <typename T>
672 inline T hav(T const& theta)
674 T const half = T(0.5);
675 T const sn = sin(half * theta);
680 \brief Short utility to return the square
682 \param value Value to calculate the square from
683 \return The squared value
685 template <typename T>
686 inline T sqr(T const& value)
688 return value * value;
692 \brief Short utility to return the square root
694 \param value Value to calculate the square root from
695 \return The square root value
697 template <typename T>
698 inline typename detail::square_root<T>::return_type
701 return detail::square_root
703 T, boost::is_fundamental<T>::value
708 \brief Short utility to return the modulo of two values
710 \param value1 First value
711 \param value2 Second value
712 \return The result of the modulo operation on the (ordered) pair
715 template <typename T>
716 inline typename detail::modulo<T>::return_type
717 mod(T const& value1, T const& value2)
719 return detail::modulo
721 T, boost::is_fundamental<T>::value
722 >::apply(value1, value2);
726 \brief Short utility to workaround gcc/clang problem that abs is converting to integer
727 and that older versions of MSVC does not support abs of long long...
731 inline T abs(T const& value)
733 return detail::abs<T>::apply(value);
737 \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
740 template <typename T>
741 inline int sign(T const& value)
744 return value > zero ? 1 : value < zero ? -1 : 0;
748 \brief Short utility to cast a value possibly rounding it to the nearest
751 \note If the source T is NOT an integral type and Result is an integral type
752 the value is rounded towards the closest integral value. Otherwise it's
753 casted without rounding.
755 template <typename Result, typename T>
756 inline Result rounding_cast(T const& v)
758 return detail::rounding_cast<Result, T>::apply(v);
764 }} // namespace boost::geometry
766 #endif // BOOST_GEOMETRY_UTIL_MATH_HPP