math.hpp 18 KB

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