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