Chris@16: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@16: 
Chris@101: // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
Chris@101: // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
Chris@101: // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
Chris@101: 
Chris@101: // This file was modified by Oracle on 2014, 2015.
Chris@101: // Modifications copyright (c) 2014-2015, Oracle and/or its affiliates.
Chris@101: 
Chris@101: // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
Chris@101: // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
Chris@16: 
Chris@16: // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
Chris@16: // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
Chris@16: 
Chris@16: // Use, modification and distribution is subject to the Boost Software License,
Chris@16: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
Chris@16: #define BOOST_GEOMETRY_UTIL_MATH_HPP
Chris@16: 
Chris@16: #include <cmath>
Chris@16: #include <limits>
Chris@16: 
Chris@101: #include <boost/core/ignore_unused.hpp>
Chris@101: 
Chris@16: #include <boost/math/constants/constants.hpp>
Chris@101: #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
Chris@101: #include <boost/math/special_functions/fpclassify.hpp>
Chris@101: #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
Chris@101: #include <boost/math/special_functions/round.hpp>
Chris@101: #include <boost/numeric/conversion/cast.hpp>
Chris@101: #include <boost/type_traits/is_fundamental.hpp>
Chris@16: 
Chris@16: #include <boost/geometry/util/select_most_precise.hpp>
Chris@16: 
Chris@16: namespace boost { namespace geometry
Chris@16: {
Chris@16: 
Chris@16: namespace math
Chris@16: {
Chris@16: 
Chris@16: #ifndef DOXYGEN_NO_DETAIL
Chris@16: namespace detail
Chris@16: {
Chris@16: 
Chris@101: template <typename T>
Chris@101: inline T const& greatest(T const& v1, T const& v2)
Chris@101: {
Chris@101:     return (std::max)(v1, v2);
Chris@101: }
Chris@16: 
Chris@101: template <typename T>
Chris@101: inline T const& greatest(T const& v1, T const& v2, T const& v3)
Chris@101: {
Chris@101:     return (std::max)(greatest(v1, v2), v3);
Chris@101: }
Chris@101: 
Chris@101: template <typename T>
Chris@101: inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4)
Chris@101: {
Chris@101:     return (std::max)(greatest(v1, v2, v3), v4);
Chris@101: }
Chris@101: 
Chris@101: template <typename T>
Chris@101: inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5)
Chris@101: {
Chris@101:     return (std::max)(greatest(v1, v2, v3, v4), v5);
Chris@101: }
Chris@101: 
Chris@101: 
Chris@101: template <typename T,
Chris@101:           bool IsFloatingPoint = boost::is_floating_point<T>::value>
Chris@101: struct abs
Chris@101: {
Chris@101:     static inline T apply(T const& value)
Chris@101:     {
Chris@101:         T const zero = T();
Chris@101:         return value < zero ? -value : value;
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: template <typename T>
Chris@101: struct abs<T, true>
Chris@101: {
Chris@101:     static inline T apply(T const& value)
Chris@101:     {
Chris@101:         return fabs(value);
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: 
Chris@101: struct equals_default_policy
Chris@101: {
Chris@101:     template <typename T>
Chris@101:     static inline T apply(T const& a, T const& b)
Chris@101:     {
Chris@101:         // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17
Chris@101:         return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1));
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: template <typename T,
Chris@101:           bool IsFloatingPoint = boost::is_floating_point<T>::value>
Chris@101: struct equals_factor_policy
Chris@101: {
Chris@101:     equals_factor_policy()
Chris@101:         : factor(1) {}
Chris@101:     explicit equals_factor_policy(T const& v)
Chris@101:         : factor(greatest(abs<T>::apply(v), T(1)))
Chris@101:     {}
Chris@101:     equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3)
Chris@101:         : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1),
Chris@101:                           abs<T>::apply(v2), abs<T>::apply(v3),
Chris@101:                           T(1)))
Chris@101:     {}
Chris@101: 
Chris@101:     T const& apply(T const&, T const&) const
Chris@101:     {
Chris@101:         return factor;
Chris@101:     }
Chris@101: 
Chris@101:     T factor;
Chris@101: };
Chris@101: 
Chris@101: template <typename T>
Chris@101: struct equals_factor_policy<T, false>
Chris@101: {
Chris@101:     equals_factor_policy() {}
Chris@101:     explicit equals_factor_policy(T const&) {}
Chris@101:     equals_factor_policy(T const& , T const& , T const& , T const& ) {}
Chris@101: 
Chris@101:     static inline T apply(T const&, T const&)
Chris@101:     {
Chris@101:         return T(1);
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: template <typename Type,
Chris@101:           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
Chris@16: struct equals
Chris@16: {
Chris@101:     template <typename Policy>
Chris@101:     static inline bool apply(Type const& a, Type const& b, Policy const&)
Chris@16:     {
Chris@16:         return a == b;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <typename Type>
Chris@16: struct equals<Type, true>
Chris@16: {
Chris@101:     template <typename Policy>
Chris@101:     static inline bool apply(Type const& a, Type const& b, Policy const& policy)
Chris@16:     {
Chris@101:         boost::ignore_unused(policy);
Chris@16: 
Chris@16:         if (a == b)
Chris@16:         {
Chris@16:             return true;
Chris@16:         }
Chris@16: 
Chris@101:         return abs<Type>::apply(a - b) <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b);
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@101: template <typename T1, typename T2, typename Policy>
Chris@101: inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy)
Chris@101: {
Chris@101:     return detail::equals
Chris@101:         <
Chris@101:             typename select_most_precise<T1, T2>::type
Chris@101:         >::apply(a, b, policy);
Chris@101: }
Chris@101: 
Chris@101: template <typename Type,
Chris@101:           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
Chris@16: struct smaller
Chris@16: {
Chris@16:     static inline bool apply(Type const& a, Type const& b)
Chris@16:     {
Chris@16:         return a < b;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <typename Type>
Chris@16: struct smaller<Type, true>
Chris@16: {
Chris@16:     static inline bool apply(Type const& a, Type const& b)
Chris@16:     {
Chris@101:         if (equals<Type, true>::apply(a, b, equals_default_policy()))
Chris@16:         {
Chris@16:             return false;
Chris@16:         }
Chris@16:         return a < b;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: 
Chris@101: template <typename Type,
Chris@101:           bool IsFloatingPoint = boost::is_floating_point<Type>::value>
Chris@101: struct equals_with_epsilon
Chris@101:     : public equals<Type, IsFloatingPoint>
Chris@101: {};
Chris@101: 
Chris@101: template
Chris@101: <
Chris@101:     typename T,
Chris@101:     bool IsFundemantal = boost::is_fundamental<T>::value /* false */
Chris@101: >
Chris@101: struct square_root
Chris@101: {
Chris@101:     typedef T return_type;
Chris@101: 
Chris@101:     static inline T apply(T const& value)
Chris@101:     {
Chris@101:         // for non-fundamental number types assume that sqrt is
Chris@101:         // defined either:
Chris@101:         // 1) at T's scope, or
Chris@101:         // 2) at global scope, or
Chris@101:         // 3) in namespace std
Chris@101:         using ::sqrt;
Chris@101:         using std::sqrt;
Chris@101: 
Chris@101:         return sqrt(value);
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: template <typename FundamentalFP>
Chris@101: struct square_root_for_fundamental_fp
Chris@101: {
Chris@101:     typedef FundamentalFP return_type;
Chris@101: 
Chris@101:     static inline FundamentalFP apply(FundamentalFP const& value)
Chris@101:     {
Chris@101: #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
Chris@101:         // This is a workaround for some 32-bit platforms.
Chris@101:         // For some of those platforms it has been reported that
Chris@101:         // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value.
Chris@101:         // For those platforms we need to define the macro
Chris@101:         // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument
Chris@101:         // to std::sqrt is checked appropriately before passed to std::sqrt
Chris@101:         if (boost::math::isfinite(value))
Chris@101:         {
Chris@101:             return std::sqrt(value);
Chris@101:         }
Chris@101:         else if (boost::math::isinf(value) && value < 0)
Chris@101:         {
Chris@101:             return -std::numeric_limits<FundamentalFP>::quiet_NaN();
Chris@101:         }
Chris@101:         return value;
Chris@101: #else
Chris@101:         // for fundamental floating point numbers use std::sqrt
Chris@101:         return std::sqrt(value);
Chris@101: #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: template <>
Chris@101: struct square_root<float, true>
Chris@101:     : square_root_for_fundamental_fp<float>
Chris@101: {
Chris@101: };
Chris@101: 
Chris@101: template <>
Chris@101: struct square_root<double, true>
Chris@101:     : square_root_for_fundamental_fp<double>
Chris@101: {
Chris@101: };
Chris@101: 
Chris@101: template <>
Chris@101: struct square_root<long double, true>
Chris@101:     : square_root_for_fundamental_fp<long double>
Chris@101: {
Chris@101: };
Chris@101: 
Chris@101: template <typename T>
Chris@101: struct square_root<T, true>
Chris@101: {
Chris@101:     typedef double return_type;
Chris@101: 
Chris@101:     static inline double apply(T const& value)
Chris@101:     {
Chris@101:         // for all other fundamental number types use also std::sqrt
Chris@101:         //
Chris@101:         // Note: in C++98 the only other possibility is double;
Chris@101:         //       in C++11 there are also overloads for integral types;
Chris@101:         //       this specialization works for those as well.
Chris@101:         return square_root_for_fundamental_fp
Chris@101:             <
Chris@101:                 double
Chris@101:             >::apply(boost::numeric_cast<double>(value));
Chris@101:     }
Chris@101: };
Chris@16: 
Chris@16: 
Chris@16: /*!
Chris@16: \brief Short construct to enable partial specialization for PI, currently not possible in Math.
Chris@16: */
Chris@16: template <typename T>
Chris@16: struct define_pi
Chris@16: {
Chris@16:     static inline T apply()
Chris@16:     {
Chris@16:         // Default calls Boost.Math
Chris@16:         return boost::math::constants::pi<T>();
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <typename T>
Chris@16: struct relaxed_epsilon
Chris@16: {
Chris@16:     static inline T apply(const T& factor)
Chris@16:     {
Chris@16:         return factor * std::numeric_limits<T>::epsilon();
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@101: // ItoF ItoI FtoF
Chris@101: template <typename Result, typename Source,
Chris@101:           bool ResultIsInteger = std::numeric_limits<Result>::is_integer,
Chris@101:           bool SourceIsInteger = std::numeric_limits<Source>::is_integer>
Chris@101: struct round
Chris@101: {
Chris@101:     static inline Result apply(Source const& v)
Chris@101:     {
Chris@101:         return boost::numeric_cast<Result>(v);
Chris@101:     }
Chris@101: };
Chris@101: 
Chris@101: // FtoI
Chris@101: template <typename Result, typename Source>
Chris@101: struct round<Result, Source, true, false>
Chris@101: {
Chris@101:     static inline Result apply(Source const& v)
Chris@101:     {
Chris@101:         namespace bmp = boost::math::policies;
Chris@101:         // ignore rounding errors for backward compatibility
Chris@101:         typedef bmp::policy< bmp::rounding_error<bmp::ignore_error> > policy;
Chris@101:         return boost::numeric_cast<Result>(boost::math::round(v, policy()));
Chris@101:     }
Chris@101: };
Chris@16: 
Chris@16: } // namespace detail
Chris@16: #endif
Chris@16: 
Chris@16: 
Chris@16: template <typename T>
Chris@16: inline T pi() { return detail::define_pi<T>::apply(); }
Chris@16: 
Chris@16: template <typename T>
Chris@16: inline T relaxed_epsilon(T const& factor)
Chris@16: {
Chris@16:     return detail::relaxed_epsilon<T>::apply(factor);
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: // Maybe replace this by boost equals or boost ublas numeric equals or so
Chris@16: 
Chris@16: /*!
Chris@16:     \brief returns true if both arguments are equal.
Chris@16:     \ingroup utility
Chris@16:     \param a first argument
Chris@16:     \param b second argument
Chris@16:     \return true if a == b
Chris@16:     \note If both a and b are of an integral type, comparison is done by ==.
Chris@16:     If one of the types is floating point, comparison is done by abs and
Chris@16:     comparing with epsilon. If one of the types is non-fundamental, it might
Chris@16:     be a high-precision number and comparison is done using the == operator
Chris@16:     of that class.
Chris@16: */
Chris@16: 
Chris@16: template <typename T1, typename T2>
Chris@16: inline bool equals(T1 const& a, T2 const& b)
Chris@16: {
Chris@16:     return detail::equals
Chris@16:         <
Chris@101:             typename select_most_precise<T1, T2>::type
Chris@101:         >::apply(a, b, detail::equals_default_policy());
Chris@16: }
Chris@16: 
Chris@16: template <typename T1, typename T2>
Chris@16: inline bool equals_with_epsilon(T1 const& a, T2 const& b)
Chris@16: {
Chris@16:     return detail::equals_with_epsilon
Chris@16:         <
Chris@101:             typename select_most_precise<T1, T2>::type
Chris@101:         >::apply(a, b, detail::equals_default_policy());
Chris@16: }
Chris@16: 
Chris@16: template <typename T1, typename T2>
Chris@16: inline bool smaller(T1 const& a, T2 const& b)
Chris@16: {
Chris@16:     return detail::smaller
Chris@16:         <
Chris@101:             typename select_most_precise<T1, T2>::type
Chris@16:         >::apply(a, b);
Chris@16: }
Chris@16: 
Chris@16: template <typename T1, typename T2>
Chris@16: inline bool larger(T1 const& a, T2 const& b)
Chris@16: {
Chris@16:     return detail::smaller
Chris@16:         <
Chris@101:             typename select_most_precise<T1, T2>::type
Chris@16:         >::apply(b, a);
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: double const d2r = geometry::math::pi<double>() / 180.0;
Chris@16: double const r2d = 1.0 / d2r;
Chris@16: 
Chris@16: /*!
Chris@16:     \brief Calculates the haversine of an angle
Chris@16:     \ingroup utility
Chris@16:     \note See http://en.wikipedia.org/wiki/Haversine_formula
Chris@16:     haversin(alpha) = sin2(alpha/2)
Chris@16: */
Chris@16: template <typename T>
Chris@16: inline T hav(T const& theta)
Chris@16: {
Chris@16:     T const half = T(0.5);
Chris@16:     T const sn = sin(half * theta);
Chris@16:     return sn * sn;
Chris@16: }
Chris@16: 
Chris@16: /*!
Chris@16: \brief Short utility to return the square
Chris@16: \ingroup utility
Chris@16: \param value Value to calculate the square from
Chris@16: \return The squared value
Chris@16: */
Chris@16: template <typename T>
Chris@16: inline T sqr(T const& value)
Chris@16: {
Chris@16:     return value * value;
Chris@16: }
Chris@16: 
Chris@16: /*!
Chris@101: \brief Short utility to return the square root
Chris@101: \ingroup utility
Chris@101: \param value Value to calculate the square root from
Chris@101: \return The square root value
Chris@101: */
Chris@101: template <typename T>
Chris@101: inline typename detail::square_root<T>::return_type
Chris@101: sqrt(T const& value)
Chris@101: {
Chris@101:     return detail::square_root
Chris@101:         <
Chris@101:             T, boost::is_fundamental<T>::value
Chris@101:         >::apply(value);
Chris@101: }
Chris@101: 
Chris@101: /*!
Chris@16: \brief Short utility to workaround gcc/clang problem that abs is converting to integer
Chris@16:        and that older versions of MSVC does not support abs of long long...
Chris@16: \ingroup utility
Chris@16: */
Chris@16: template<typename T>
Chris@16: inline T abs(T const& value)
Chris@16: {
Chris@101:     return detail::abs<T>::apply(value);
Chris@16: }
Chris@16: 
Chris@16: /*!
Chris@16: \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
Chris@16: \ingroup utility
Chris@16: */
Chris@16: template <typename T>
Chris@101: static inline int sign(T const& value)
Chris@16: {
Chris@16:     T const zero = T();
Chris@16:     return value > zero ? 1 : value < zero ? -1 : 0;
Chris@16: }
Chris@16: 
Chris@101: /*!
Chris@101: \brief Short utility to calculate the rounded value of a number.
Chris@101: \ingroup utility
Chris@101: \note If the source T is NOT an integral type and Result is an integral type
Chris@101:       the value is rounded towards the closest integral value. Otherwise it's
Chris@101:       casted.
Chris@101: */
Chris@101: template <typename Result, typename T>
Chris@101: inline Result round(T const& v)
Chris@101: {
Chris@101:     return detail::round<Result, T>::apply(v);
Chris@101: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: 
Chris@16: 
Chris@16: }} // namespace boost::geometry
Chris@16: 
Chris@16: #endif // BOOST_GEOMETRY_UTIL_MATH_HPP