ceph/src/boost/boost/geometry/util/math.hpp

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