math.hpp 5.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
  3. // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
  4. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
  5. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  6. // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  7. // Use, modification and distribution is subject to the Boost Software License,
  8. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  9. // http://www.boost.org/LICENSE_1_0.txt)
  10. #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
  11. #define BOOST_GEOMETRY_UTIL_MATH_HPP
  12. #include <cmath>
  13. #include <limits>
  14. #include <boost/math/constants/constants.hpp>
  15. #include <boost/geometry/util/select_most_precise.hpp>
  16. namespace boost { namespace geometry
  17. {
  18. namespace math
  19. {
  20. #ifndef DOXYGEN_NO_DETAIL
  21. namespace detail
  22. {
  23. template <typename Type, bool IsFloatingPoint>
  24. struct equals
  25. {
  26. static inline bool apply(Type const& a, Type const& b)
  27. {
  28. return a == b;
  29. }
  30. };
  31. template <typename Type>
  32. struct equals<Type, true>
  33. {
  34. static inline Type get_max(Type const& a, Type const& b, Type const& c)
  35. {
  36. return (std::max)((std::max)(a, b), c);
  37. }
  38. static inline bool apply(Type const& a, Type const& b)
  39. {
  40. if (a == b)
  41. {
  42. return true;
  43. }
  44. // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17,
  45. // FUTURE: replace by some boost tool or boost::test::close_at_tolerance
  46. return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * get_max(std::abs(a), std::abs(b), 1.0);
  47. }
  48. };
  49. template <typename Type, bool IsFloatingPoint>
  50. struct smaller
  51. {
  52. static inline bool apply(Type const& a, Type const& b)
  53. {
  54. return a < b;
  55. }
  56. };
  57. template <typename Type>
  58. struct smaller<Type, true>
  59. {
  60. static inline bool apply(Type const& a, Type const& b)
  61. {
  62. if (equals<Type, true>::apply(a, b))
  63. {
  64. return false;
  65. }
  66. return a < b;
  67. }
  68. };
  69. template <typename Type, bool IsFloatingPoint>
  70. struct equals_with_epsilon : public equals<Type, IsFloatingPoint> {};
  71. /*!
  72. \brief Short construct to enable partial specialization for PI, currently not possible in Math.
  73. */
  74. template <typename T>
  75. struct define_pi
  76. {
  77. static inline T apply()
  78. {
  79. // Default calls Boost.Math
  80. return boost::math::constants::pi<T>();
  81. }
  82. };
  83. template <typename T>
  84. struct relaxed_epsilon
  85. {
  86. static inline T apply(const T& factor)
  87. {
  88. return factor * std::numeric_limits<T>::epsilon();
  89. }
  90. };
  91. } // namespace detail
  92. #endif
  93. template <typename T>
  94. inline T pi() { return detail::define_pi<T>::apply(); }
  95. template <typename T>
  96. inline T relaxed_epsilon(T const& factor)
  97. {
  98. return detail::relaxed_epsilon<T>::apply(factor);
  99. }
  100. // Maybe replace this by boost equals or boost ublas numeric equals or so
  101. /*!
  102. \brief returns true if both arguments are equal.
  103. \ingroup utility
  104. \param a first argument
  105. \param b second argument
  106. \return true if a == b
  107. \note If both a and b are of an integral type, comparison is done by ==.
  108. If one of the types is floating point, comparison is done by abs and
  109. comparing with epsilon. If one of the types is non-fundamental, it might
  110. be a high-precision number and comparison is done using the == operator
  111. of that class.
  112. */
  113. template <typename T1, typename T2>
  114. inline bool equals(T1 const& a, T2 const& b)
  115. {
  116. typedef typename select_most_precise<T1, T2>::type select_type;
  117. return detail::equals
  118. <
  119. select_type,
  120. boost::is_floating_point<select_type>::type::value
  121. >::apply(a, b);
  122. }
  123. template <typename T1, typename T2>
  124. inline bool equals_with_epsilon(T1 const& a, T2 const& b)
  125. {
  126. typedef typename select_most_precise<T1, T2>::type select_type;
  127. return detail::equals_with_epsilon
  128. <
  129. select_type,
  130. boost::is_floating_point<select_type>::type::value
  131. >::apply(a, b);
  132. }
  133. template <typename T1, typename T2>
  134. inline bool smaller(T1 const& a, T2 const& b)
  135. {
  136. typedef typename select_most_precise<T1, T2>::type select_type;
  137. return detail::smaller
  138. <
  139. select_type,
  140. boost::is_floating_point<select_type>::type::value
  141. >::apply(a, b);
  142. }
  143. template <typename T1, typename T2>
  144. inline bool larger(T1 const& a, T2 const& b)
  145. {
  146. typedef typename select_most_precise<T1, T2>::type select_type;
  147. return detail::smaller
  148. <
  149. select_type,
  150. boost::is_floating_point<select_type>::type::value
  151. >::apply(b, a);
  152. }
  153. double const d2r = geometry::math::pi<double>() / 180.0;
  154. double const r2d = 1.0 / d2r;
  155. /*!
  156. \brief Calculates the haversine of an angle
  157. \ingroup utility
  158. \note See http://en.wikipedia.org/wiki/Haversine_formula
  159. haversin(alpha) = sin2(alpha/2)
  160. */
  161. template <typename T>
  162. inline T hav(T const& theta)
  163. {
  164. T const half = T(0.5);
  165. T const sn = sin(half * theta);
  166. return sn * sn;
  167. }
  168. /*!
  169. \brief Short utility to return the square
  170. \ingroup utility
  171. \param value Value to calculate the square from
  172. \return The squared value
  173. */
  174. template <typename T>
  175. inline T sqr(T const& value)
  176. {
  177. return value * value;
  178. }
  179. /*!
  180. \brief Short utility to workaround gcc/clang problem that abs is converting to integer
  181. and that older versions of MSVC does not support abs of long long...
  182. \ingroup utility
  183. */
  184. template<typename T>
  185. inline T abs(T const& value)
  186. {
  187. T const zero = T();
  188. return value < zero ? -value : value;
  189. }
  190. /*!
  191. \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
  192. \ingroup utility
  193. */
  194. template <typename T>
  195. static inline int sign(T const& value)
  196. {
  197. T const zero = T();
  198. return value > zero ? 1 : value < zero ? -1 : 0;
  199. }
  200. } // namespace math
  201. }} // namespace boost::geometry
  202. #endif // BOOST_GEOMETRY_UTIL_MATH_HPP