math.hpp 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
  3. // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
  4. // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
  5. // This file was modified by Oracle on 2014, 2015.
  6. // Modifications copyright (c) 2014-2015, Oracle and/or its affiliates.
  7. // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
  8. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  9. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  10. // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  11. // Use, modification and distribution is subject to the Boost Software License,
  12. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  13. // http://www.boost.org/LICENSE_1_0.txt)
  14. #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
  15. #define BOOST_GEOMETRY_UTIL_MATH_HPP
  16. #include <cmath>
  17. #include <limits>
  18. #include <boost/core/ignore_unused.hpp>
  19. #include <boost/math/constants/constants.hpp>
  20. #include <boost/math/special_functions/fpclassify.hpp>
  21. //#include <boost/math/special_functions/round.hpp>
  22. #include <boost/numeric/conversion/cast.hpp>
  23. #include <boost/type_traits/is_fundamental.hpp>
  24. #include <boost/type_traits/is_integral.hpp>
  25. #include <boost/geometry/core/cs.hpp>
  26. #include <boost/geometry/util/select_most_precise.hpp>
  27. namespace boost { namespace geometry
  28. {
  29. namespace math
  30. {
  31. #ifndef DOXYGEN_NO_DETAIL
  32. namespace detail
  33. {
  34. template <typename T>
  35. inline T const& greatest(T const& v1, T const& v2)
  36. {
  37. return (std::max)(v1, v2);
  38. }
  39. template <typename T>
  40. inline T const& greatest(T const& v1, T const& v2, T const& v3)
  41. {
  42. return (std::max)(greatest(v1, v2), v3);
  43. }
  44. template <typename T>
  45. inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4)
  46. {
  47. return (std::max)(greatest(v1, v2, v3), v4);
  48. }
  49. template <typename T>
  50. inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T const& v5)
  51. {
  52. return (std::max)(greatest(v1, v2, v3, v4), v5);
  53. }
  54. template <typename T,
  55. bool IsFloatingPoint = boost::is_floating_point<T>::value>
  56. struct abs
  57. {
  58. static inline T apply(T const& value)
  59. {
  60. T const zero = T();
  61. return value < zero ? -value : value;
  62. }
  63. };
  64. template <typename T>
  65. struct abs<T, true>
  66. {
  67. static inline T apply(T const& value)
  68. {
  69. using ::fabs;
  70. using std::fabs; // for long double
  71. return fabs(value);
  72. }
  73. };
  74. struct equals_default_policy
  75. {
  76. template <typename T>
  77. static inline T apply(T const& a, T const& b)
  78. {
  79. // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17
  80. return greatest(abs<T>::apply(a), abs<T>::apply(b), T(1));
  81. }
  82. };
  83. template <typename T,
  84. bool IsFloatingPoint = boost::is_floating_point<T>::value>
  85. struct equals_factor_policy
  86. {
  87. equals_factor_policy()
  88. : factor(1) {}
  89. explicit equals_factor_policy(T const& v)
  90. : factor(greatest(abs<T>::apply(v), T(1)))
  91. {}
  92. equals_factor_policy(T const& v0, T const& v1, T const& v2, T const& v3)
  93. : factor(greatest(abs<T>::apply(v0), abs<T>::apply(v1),
  94. abs<T>::apply(v2), abs<T>::apply(v3),
  95. T(1)))
  96. {}
  97. T const& apply(T const&, T const&) const
  98. {
  99. return factor;
  100. }
  101. T factor;
  102. };
  103. template <typename T>
  104. struct equals_factor_policy<T, false>
  105. {
  106. equals_factor_policy() {}
  107. explicit equals_factor_policy(T const&) {}
  108. equals_factor_policy(T const& , T const& , T const& , T const& ) {}
  109. static inline T apply(T const&, T const&)
  110. {
  111. return T(1);
  112. }
  113. };
  114. template <typename Type,
  115. bool IsFloatingPoint = boost::is_floating_point<Type>::value>
  116. struct equals
  117. {
  118. template <typename Policy>
  119. static inline bool apply(Type const& a, Type const& b, Policy const&)
  120. {
  121. return a == b;
  122. }
  123. };
  124. template <typename Type>
  125. struct equals<Type, true>
  126. {
  127. template <typename Policy>
  128. static inline bool apply(Type const& a, Type const& b, Policy const& policy)
  129. {
  130. boost::ignore_unused(policy);
  131. if (a == b)
  132. {
  133. return true;
  134. }
  135. if (boost::math::isfinite(a) && boost::math::isfinite(b))
  136. {
  137. // If a is INF and b is e.g. 0, the expression below returns true
  138. // but the values are obviously not equal, hence the condition
  139. return abs<Type>::apply(a - b)
  140. <= std::numeric_limits<Type>::epsilon() * policy.apply(a, b);
  141. }
  142. else
  143. {
  144. return a == b;
  145. }
  146. }
  147. };
  148. template <typename T1, typename T2, typename Policy>
  149. inline bool equals_by_policy(T1 const& a, T2 const& b, Policy const& policy)
  150. {
  151. return detail::equals
  152. <
  153. typename select_most_precise<T1, T2>::type
  154. >::apply(a, b, policy);
  155. }
  156. template <typename Type,
  157. bool IsFloatingPoint = boost::is_floating_point<Type>::value>
  158. struct smaller
  159. {
  160. static inline bool apply(Type const& a, Type const& b)
  161. {
  162. return a < b;
  163. }
  164. };
  165. template <typename Type>
  166. struct smaller<Type, true>
  167. {
  168. static inline bool apply(Type const& a, Type const& b)
  169. {
  170. if (!(a < b)) // a >= b
  171. {
  172. return false;
  173. }
  174. return ! equals<Type, true>::apply(b, a, equals_default_policy());
  175. }
  176. };
  177. template <typename Type,
  178. bool IsFloatingPoint = boost::is_floating_point<Type>::value>
  179. struct smaller_or_equals
  180. {
  181. static inline bool apply(Type const& a, Type const& b)
  182. {
  183. return a <= b;
  184. }
  185. };
  186. template <typename Type>
  187. struct smaller_or_equals<Type, true>
  188. {
  189. static inline bool apply(Type const& a, Type const& b)
  190. {
  191. if (a <= b)
  192. {
  193. return true;
  194. }
  195. return equals<Type, true>::apply(a, b, equals_default_policy());
  196. }
  197. };
  198. template <typename Type,
  199. bool IsFloatingPoint = boost::is_floating_point<Type>::value>
  200. struct equals_with_epsilon
  201. : public equals<Type, IsFloatingPoint>
  202. {};
  203. template
  204. <
  205. typename T,
  206. bool IsFundemantal = boost::is_fundamental<T>::value /* false */
  207. >
  208. struct square_root
  209. {
  210. typedef T return_type;
  211. static inline T apply(T const& value)
  212. {
  213. // for non-fundamental number types assume that sqrt is
  214. // defined either:
  215. // 1) at T's scope, or
  216. // 2) at global scope, or
  217. // 3) in namespace std
  218. using ::sqrt;
  219. using std::sqrt;
  220. return sqrt(value);
  221. }
  222. };
  223. template <typename FundamentalFP>
  224. struct square_root_for_fundamental_fp
  225. {
  226. typedef FundamentalFP return_type;
  227. static inline FundamentalFP apply(FundamentalFP const& value)
  228. {
  229. #ifdef BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
  230. // This is a workaround for some 32-bit platforms.
  231. // For some of those platforms it has been reported that
  232. // std::sqrt(nan) and/or std::sqrt(-nan) returns a finite value.
  233. // For those platforms we need to define the macro
  234. // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS so that the argument
  235. // to std::sqrt is checked appropriately before passed to std::sqrt
  236. if (boost::math::isfinite(value))
  237. {
  238. return std::sqrt(value);
  239. }
  240. else if (boost::math::isinf(value) && value < 0)
  241. {
  242. return -std::numeric_limits<FundamentalFP>::quiet_NaN();
  243. }
  244. return value;
  245. #else
  246. // for fundamental floating point numbers use std::sqrt
  247. return std::sqrt(value);
  248. #endif // BOOST_GEOMETRY_SQRT_CHECK_FINITENESS
  249. }
  250. };
  251. template <>
  252. struct square_root<float, true>
  253. : square_root_for_fundamental_fp<float>
  254. {
  255. };
  256. template <>
  257. struct square_root<double, true>
  258. : square_root_for_fundamental_fp<double>
  259. {
  260. };
  261. template <>
  262. struct square_root<long double, true>
  263. : square_root_for_fundamental_fp<long double>
  264. {
  265. };
  266. template <typename T>
  267. struct square_root<T, true>
  268. {
  269. typedef double return_type;
  270. static inline double apply(T const& value)
  271. {
  272. // for all other fundamental number types use also std::sqrt
  273. //
  274. // Note: in C++98 the only other possibility is double;
  275. // in C++11 there are also overloads for integral types;
  276. // this specialization works for those as well.
  277. return square_root_for_fundamental_fp
  278. <
  279. double
  280. >::apply(boost::numeric_cast<double>(value));
  281. }
  282. };
  283. template
  284. <
  285. typename T,
  286. bool IsFundemantal = boost::is_fundamental<T>::value /* false */
  287. >
  288. struct modulo
  289. {
  290. typedef T return_type;
  291. static inline T apply(T const& value1, T const& value2)
  292. {
  293. // for non-fundamental number types assume that a free
  294. // function mod() is defined either:
  295. // 1) at T's scope, or
  296. // 2) at global scope
  297. return mod(value1, value2);
  298. }
  299. };
  300. template
  301. <
  302. typename Fundamental,
  303. bool IsIntegral = boost::is_integral<Fundamental>::value
  304. >
  305. struct modulo_for_fundamental
  306. {
  307. typedef Fundamental return_type;
  308. static inline Fundamental apply(Fundamental const& value1,
  309. Fundamental const& value2)
  310. {
  311. return value1 % value2;
  312. }
  313. };
  314. // specialization for floating-point numbers
  315. template <typename Fundamental>
  316. struct modulo_for_fundamental<Fundamental, false>
  317. {
  318. typedef Fundamental return_type;
  319. static inline Fundamental apply(Fundamental const& value1,
  320. Fundamental const& value2)
  321. {
  322. return std::fmod(value1, value2);
  323. }
  324. };
  325. // specialization for fundamental number type
  326. template <typename Fundamental>
  327. struct modulo<Fundamental, true>
  328. : modulo_for_fundamental<Fundamental>
  329. {};
  330. /*!
  331. \brief Short constructs to enable partial specialization for PI, 2*PI
  332. and PI/2, currently not possible in Math.
  333. */
  334. template <typename T>
  335. struct define_pi
  336. {
  337. static inline T apply()
  338. {
  339. // Default calls Boost.Math
  340. return boost::math::constants::pi<T>();
  341. }
  342. };
  343. template <typename T>
  344. struct define_two_pi
  345. {
  346. static inline T apply()
  347. {
  348. // Default calls Boost.Math
  349. return boost::math::constants::two_pi<T>();
  350. }
  351. };
  352. template <typename T>
  353. struct define_half_pi
  354. {
  355. static inline T apply()
  356. {
  357. // Default calls Boost.Math
  358. return boost::math::constants::half_pi<T>();
  359. }
  360. };
  361. template <typename T>
  362. struct relaxed_epsilon
  363. {
  364. static inline T apply(const T& factor)
  365. {
  366. return factor * std::numeric_limits<T>::epsilon();
  367. }
  368. };
  369. // This must be consistent with math::equals.
  370. // By default math::equals() scales the error by epsilon using the greater of
  371. // compared values but here is only one value, though it should work the same way.
  372. // (a-a) <= max(a, a) * EPS -> 0 <= a*EPS
  373. // (a+da-a) <= max(a+da, a) * EPS -> da <= (a+da)*EPS
  374. template <typename T, bool IsFloat = boost::is_floating_point<T>::value>
  375. struct scaled_epsilon
  376. {
  377. static inline T apply(T const& val)
  378. {
  379. return (std::max)(abs<T>::apply(val), T(1))
  380. * std::numeric_limits<T>::epsilon();
  381. }
  382. };
  383. template <typename T>
  384. struct scaled_epsilon<T, false>
  385. {
  386. static inline T apply(T const&)
  387. {
  388. return T(0);
  389. }
  390. };
  391. // ItoF ItoI FtoF
  392. template <typename Result, typename Source,
  393. bool ResultIsInteger = std::numeric_limits<Result>::is_integer,
  394. bool SourceIsInteger = std::numeric_limits<Source>::is_integer>
  395. struct rounding_cast
  396. {
  397. static inline Result apply(Source const& v)
  398. {
  399. return boost::numeric_cast<Result>(v);
  400. }
  401. };
  402. // TtoT
  403. template <typename Source, bool ResultIsInteger, bool SourceIsInteger>
  404. struct rounding_cast<Source, Source, ResultIsInteger, SourceIsInteger>
  405. {
  406. static inline Source apply(Source const& v)
  407. {
  408. return v;
  409. }
  410. };
  411. // FtoI
  412. template <typename Result, typename Source>
  413. struct rounding_cast<Result, Source, true, false>
  414. {
  415. static inline Result apply(Source const& v)
  416. {
  417. return boost::numeric_cast<Result>(v < Source(0) ?
  418. v - Source(0.5) :
  419. v + Source(0.5));
  420. }
  421. };
  422. } // namespace detail
  423. #endif
  424. template <typename T>
  425. inline T pi() { return detail::define_pi<T>::apply(); }
  426. template <typename T>
  427. inline T two_pi() { return detail::define_two_pi<T>::apply(); }
  428. template <typename T>
  429. inline T half_pi() { return detail::define_half_pi<T>::apply(); }
  430. template <typename T>
  431. inline T relaxed_epsilon(T const& factor)
  432. {
  433. return detail::relaxed_epsilon<T>::apply(factor);
  434. }
  435. template <typename T>
  436. inline T scaled_epsilon(T const& value)
  437. {
  438. return detail::scaled_epsilon<T>::apply(value);
  439. }
  440. // Maybe replace this by boost equals or boost ublas numeric equals or so
  441. /*!
  442. \brief returns true if both arguments are equal.
  443. \ingroup utility
  444. \param a first argument
  445. \param b second argument
  446. \return true if a == b
  447. \note If both a and b are of an integral type, comparison is done by ==.
  448. If one of the types is floating point, comparison is done by abs and
  449. comparing with epsilon. If one of the types is non-fundamental, it might
  450. be a high-precision number and comparison is done using the == operator
  451. of that class.
  452. */
  453. template <typename T1, typename T2>
  454. inline bool equals(T1 const& a, T2 const& b)
  455. {
  456. return detail::equals
  457. <
  458. typename select_most_precise<T1, T2>::type
  459. >::apply(a, b, detail::equals_default_policy());
  460. }
  461. template <typename T1, typename T2>
  462. inline bool equals_with_epsilon(T1 const& a, T2 const& b)
  463. {
  464. return detail::equals_with_epsilon
  465. <
  466. typename select_most_precise<T1, T2>::type
  467. >::apply(a, b, detail::equals_default_policy());
  468. }
  469. template <typename T1, typename T2>
  470. inline bool smaller(T1 const& a, T2 const& b)
  471. {
  472. return detail::smaller
  473. <
  474. typename select_most_precise<T1, T2>::type
  475. >::apply(a, b);
  476. }
  477. template <typename T1, typename T2>
  478. inline bool larger(T1 const& a, T2 const& b)
  479. {
  480. return detail::smaller
  481. <
  482. typename select_most_precise<T1, T2>::type
  483. >::apply(b, a);
  484. }
  485. template <typename T1, typename T2>
  486. inline bool smaller_or_equals(T1 const& a, T2 const& b)
  487. {
  488. return detail::smaller_or_equals
  489. <
  490. typename select_most_precise<T1, T2>::type
  491. >::apply(a, b);
  492. }
  493. template <typename T1, typename T2>
  494. inline bool larger_or_equals(T1 const& a, T2 const& b)
  495. {
  496. return detail::smaller_or_equals
  497. <
  498. typename select_most_precise<T1, T2>::type
  499. >::apply(b, a);
  500. }
  501. template <typename T>
  502. inline T d2r()
  503. {
  504. static T const conversion_coefficient = geometry::math::pi<T>() / T(180.0);
  505. return conversion_coefficient;
  506. }
  507. template <typename T>
  508. inline T r2d()
  509. {
  510. static T const conversion_coefficient = T(180.0) / geometry::math::pi<T>();
  511. return conversion_coefficient;
  512. }
  513. #ifndef DOXYGEN_NO_DETAIL
  514. namespace detail {
  515. template <typename DegreeOrRadian>
  516. struct as_radian
  517. {
  518. template <typename T>
  519. static inline T apply(T const& value)
  520. {
  521. return value;
  522. }
  523. };
  524. template <>
  525. struct as_radian<degree>
  526. {
  527. template <typename T>
  528. static inline T apply(T const& value)
  529. {
  530. return value * d2r<T>();
  531. }
  532. };
  533. template <typename DegreeOrRadian>
  534. struct from_radian
  535. {
  536. template <typename T>
  537. static inline T apply(T const& value)
  538. {
  539. return value;
  540. }
  541. };
  542. template <>
  543. struct from_radian<degree>
  544. {
  545. template <typename T>
  546. static inline T apply(T const& value)
  547. {
  548. return value * r2d<T>();
  549. }
  550. };
  551. } // namespace detail
  552. #endif
  553. template <typename DegreeOrRadian, typename T>
  554. inline T as_radian(T const& value)
  555. {
  556. return detail::as_radian<DegreeOrRadian>::apply(value);
  557. }
  558. template <typename DegreeOrRadian, typename T>
  559. inline T from_radian(T const& value)
  560. {
  561. return detail::from_radian<DegreeOrRadian>::apply(value);
  562. }
  563. /*!
  564. \brief Calculates the haversine of an angle
  565. \ingroup utility
  566. \note See http://en.wikipedia.org/wiki/Haversine_formula
  567. haversin(alpha) = sin2(alpha/2)
  568. */
  569. template <typename T>
  570. inline T hav(T const& theta)
  571. {
  572. T const half = T(0.5);
  573. T const sn = sin(half * theta);
  574. return sn * sn;
  575. }
  576. /*!
  577. \brief Short utility to return the square
  578. \ingroup utility
  579. \param value Value to calculate the square from
  580. \return The squared value
  581. */
  582. template <typename T>
  583. inline T sqr(T const& value)
  584. {
  585. return value * value;
  586. }
  587. /*!
  588. \brief Short utility to return the square root
  589. \ingroup utility
  590. \param value Value to calculate the square root from
  591. \return The square root value
  592. */
  593. template <typename T>
  594. inline typename detail::square_root<T>::return_type
  595. sqrt(T const& value)
  596. {
  597. return detail::square_root
  598. <
  599. T, boost::is_fundamental<T>::value
  600. >::apply(value);
  601. }
  602. /*!
  603. \brief Short utility to return the modulo of two values
  604. \ingroup utility
  605. \param value1 First value
  606. \param value2 Second value
  607. \return The result of the modulo operation on the (ordered) pair
  608. (value1, value2)
  609. */
  610. template <typename T>
  611. inline typename detail::modulo<T>::return_type
  612. mod(T const& value1, T const& value2)
  613. {
  614. return detail::modulo
  615. <
  616. T, boost::is_fundamental<T>::value
  617. >::apply(value1, value2);
  618. }
  619. /*!
  620. \brief Short utility to workaround gcc/clang problem that abs is converting to integer
  621. and that older versions of MSVC does not support abs of long long...
  622. \ingroup utility
  623. */
  624. template<typename T>
  625. inline T abs(T const& value)
  626. {
  627. return detail::abs<T>::apply(value);
  628. }
  629. /*!
  630. \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
  631. \ingroup utility
  632. */
  633. template <typename T>
  634. inline int sign(T const& value)
  635. {
  636. T const zero = T();
  637. return value > zero ? 1 : value < zero ? -1 : 0;
  638. }
  639. /*!
  640. \brief Short utility to cast a value possibly rounding it to the nearest
  641. integral value.
  642. \ingroup utility
  643. \note If the source T is NOT an integral type and Result is an integral type
  644. the value is rounded towards the closest integral value. Otherwise it's
  645. casted without rounding.
  646. */
  647. template <typename Result, typename T>
  648. inline Result rounding_cast(T const& v)
  649. {
  650. return detail::rounding_cast<Result, T>::apply(v);
  651. }
  652. } // namespace math
  653. }} // namespace boost::geometry
  654. #endif // BOOST_GEOMETRY_UTIL_MATH_HPP