diff options
Diffstat (limited to 'boost/geometry/util/math.hpp')
-rw-r--r-- | boost/geometry/util/math.hpp | 224 |
1 files changed, 174 insertions, 50 deletions
diff --git a/boost/geometry/util/math.hpp b/boost/geometry/util/math.hpp index 22c02168ad..4042f4e4cd 100644 --- a/boost/geometry/util/math.hpp +++ b/boost/geometry/util/math.hpp @@ -1,11 +1,11 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -// Copyright (c) 2008-2014 Bruno Lalande, Paris, France. -// Copyright (c) 2009-2014 Mateusz Loskot, London, UK. +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2015. +// Modifications copyright (c) 2014-2015, Oracle and/or its affiliates. // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -23,7 +23,12 @@ #include <cmath> #include <limits> +#include <boost/core/ignore_unused.hpp> + #include <boost/math/constants/constants.hpp> +#ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS +#include <boost/math/special_functions/fpclassify.hpp> +#endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS #include <boost/math/special_functions/round.hpp> #include <boost/numeric/conversion/cast.hpp> #include <boost/type_traits/is_fundamental.hpp> @@ -40,11 +45,104 @@ namespace math namespace detail { +template <typename T> +inline T const& greatest(T const& v1, T const& v2) +{ + return (std::max)(v1, v2); +} + +template <typename T> +inline T const& greatest(T const& v1, T const& v2, T const& v3) +{ + return (std::max)(greatest(v1, v2), v3); +} + +template <typename T> +inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4) +{ + return (std::max)(greatest(v1, v2, v3), v4); +} + +template <typename T> +inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5) +{ + return (std::max)(greatest(v1, v2, v3, v4), v5); +} + + +template <typename T, + bool IsFloatingPoint = boost::is_floating_point<T>::value> +struct abs +{ + static inline T apply(T const& value) + { + T const zero = T(); + return value < zero ? -value : value; + } +}; + +template <typename T> +struct abs<T, true> +{ + static inline T apply(T const& value) + { + return fabs(value); + } +}; + + +struct equals_default_policy +{ + template <typename T> + static inline T apply(T const& a, T const& b) + { + // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17 + return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1)); + } +}; + +template <typename T, + bool IsFloatingPoint = boost::is_floating_point<T>::value> +struct equals_factor_policy +{ + equals_factor_policy() + : factor(1) {} + explicit equals_factor_policy(T const& v) + : factor(greatest(abs<T>::apply(v), T(1))) + {} + equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3) + : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1), + abs<T>::apply(v2), abs<T>::apply(v3), + T(1))) + {} + + T const& apply(T const&, T const&) const + { + return factor; + } + + T factor; +}; + +template <typename T> +struct equals_factor_policy<T, false> +{ + equals_factor_policy() {} + explicit equals_factor_policy(T const&) {} + equals_factor_policy(T const& , T const& , T const& , T const& ) {} + + static inline T apply(T const&, T const&) + { + return T(1); + } +}; -template <typename Type, bool IsFloatingPoint> +template <typename Type, + bool IsFloatingPoint = boost::is_floating_point<Type>::value> struct equals { - static inline bool apply(Type const& a, Type const& b) + template <typename Policy> + static inline bool apply(Type const& a, Type const& b, Policy const&) { return a == b; } @@ -53,25 +151,31 @@ struct equals template <typename Type> struct equals<Type, true> { - static inline Type get_max(Type const& a, Type const& b, Type const& c) + template <typename Policy> + static inline bool apply(Type const& a, Type const& b, Policy const& policy) { - return (std::max)((std::max)(a, b), c); - } + boost::ignore_unused(policy); - static inline bool apply(Type const& a, Type const& b) - { if (a == b) { return true; } - // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17, - // FUTURE: replace by some boost tool or boost::test::close_at_tolerance - return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * get_max(std::abs(a), std::abs(b), 1.0); + return abs<Type>::apply(a - b) <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b); } }; -template <typename Type, bool IsFloatingPoint> +template <typename T1, typename T2, typename Policy> +inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy) +{ + return detail::equals + < + typename select_most_precise<T1, T2>::type + >::apply(a, b, policy); +} + +template <typename Type, + bool IsFloatingPoint = boost::is_floating_point<Type>::value> struct smaller { static inline bool apply(Type const& a, Type const& b) @@ -85,7 +189,7 @@ struct smaller<Type, true> { static inline bool apply(Type const& a, Type const& b) { - if (equals<Type, true>::apply(a, b)) + if (equals<Type, true>::apply(a, b, equals_default_policy())) { return false; } @@ -94,8 +198,11 @@ struct smaller<Type, true> }; -template <typename Type, bool IsFloatingPoint> -struct equals_with_epsilon : public equals<Type, IsFloatingPoint> {}; +template <typename Type, + bool IsFloatingPoint = boost::is_floating_point<Type>::value> +struct equals_with_epsilon + : public equals<Type, IsFloatingPoint> +{}; template < @@ -120,28 +227,52 @@ struct square_root } }; -template <> -struct square_root<float, true> +template <typename FundamentalFP> +struct square_root_for_fundamental_fp { - typedef float return_type; + typedef FundamentalFP return_type; - static inline float apply(float const& value) + static inline FundamentalFP apply(FundamentalFP const& value) { - // for float use std::sqrt +#ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS + // This is a workaround for some 32-bit platforms. + // For some of those platforms it has been reported that + // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value. + // For those platforms we need to define the macro + // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument + // to std::sqrt is checked appropriately before passed to std::sqrt + if (boost::math::isfinite(value)) + { + return std::sqrt(value); + } + else if (boost::math::isinf(value) && value < 0) + { + return -std::numeric_limits<FundamentalFP>::quiet_NaN(); + } + return value; +#else + // for fundamental floating point numbers use std::sqrt return std::sqrt(value); +#endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS } }; template <> -struct square_root<long double, true> +struct square_root<float, true> + : square_root_for_fundamental_fp<float> { - typedef long double return_type; +}; - static inline long double apply(long double const& value) - { - // for long double use std::sqrt - return std::sqrt(value); - } +template <> +struct square_root<double, true> + : square_root_for_fundamental_fp<double> +{ +}; + +template <> +struct square_root<long double, true> + : square_root_for_fundamental_fp<long double> +{ }; template <typename T> @@ -156,7 +287,10 @@ struct square_root<T, true> // Note: in C++98 the only other possibility is double; // in C++11 there are also overloads for integral types; // this specialization works for those as well. - return std::sqrt(boost::numeric_cast<double>(value)); + return square_root_for_fundamental_fp + < + double + >::apply(boost::numeric_cast<double>(value)); } }; @@ -240,44 +374,36 @@ inline T relaxed_epsilon(T const& factor) template <typename T1, typename T2> inline bool equals(T1 const& a, T2 const& b) { - typedef typename select_most_precise<T1, T2>::type select_type; return detail::equals < - select_type, - boost::is_floating_point<select_type>::type::value - >::apply(a, b); + typename select_most_precise<T1, T2>::type + >::apply(a, b, detail::equals_default_policy()); } template <typename T1, typename T2> inline bool equals_with_epsilon(T1 const& a, T2 const& b) { - typedef typename select_most_precise<T1, T2>::type select_type; return detail::equals_with_epsilon < - select_type, - boost::is_floating_point<select_type>::type::value - >::apply(a, b); + typename select_most_precise<T1, T2>::type + >::apply(a, b, detail::equals_default_policy()); } template <typename T1, typename T2> inline bool smaller(T1 const& a, T2 const& b) { - typedef typename select_most_precise<T1, T2>::type select_type; return detail::smaller < - select_type, - boost::is_floating_point<select_type>::type::value + typename select_most_precise<T1, T2>::type >::apply(a, b); } template <typename T1, typename T2> inline bool larger(T1 const& a, T2 const& b) { - typedef typename select_most_precise<T1, T2>::type select_type; return detail::smaller < - select_type, - boost::is_floating_point<select_type>::type::value + typename select_most_precise<T1, T2>::type >::apply(b, a); } @@ -336,8 +462,7 @@ sqrt(T const& value) template<typename T> inline T abs(T const& value) { - T const zero = T(); - return value < zero ? -value : value; + return detail::abs<T>::apply(value); } /*! @@ -356,12 +481,11 @@ static inline int sign(T const& value) \ingroup utility \note If the source T is NOT an integral type and Result is an integral type the value is rounded towards the closest integral value. Otherwise it's - just casted. + casted. */ template <typename Result, typename T> inline Result round(T const& v) { - // NOTE: boost::round() could be used instead but it throws in some situations return detail::round<Result, T>::apply(v); } |