ceph/src/boost/boost/math/tools/agm.hpp

   1 //  (C) Copyright Nick Thompson 2020.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_AGM_HPP
   7 #define BOOST_MATH_TOOLS_AGM_HPP
   8 #include <limits>
   9 #include <cmath>
  10
  11 namespace boost { namespace math { namespace tools {
  12
  13 template<typename Real>
  14 Real agm(Real a, Real g)
  15 {
  16     using std::sqrt;
  17
  18     if (a < g)
  19     {
  20         // Mathematica, mpfr, and mpmath are all symmetric functions:
  21         return agm(g, a);
  22     }
  23     // Use: M(rx, ry) = rM(x,y)
  24     if (a <= 0 || g <= 0) {
  25         if (a < 0 || g < 0) {
  26             return std::numeric_limits<Real>::quiet_NaN();
  27         }
  28         return Real(0);
  29     }
  30
  31     // The number of correct digits doubles on each iteration.
  32     // Divide by 512 for some leeway:
  33     const Real scale = sqrt(std::numeric_limits<Real>::epsilon())/512;
  34     while (a-g > scale*g)
  35     {
  36         Real anp1 = (a + g)/2;
  37         g = sqrt(a*g);
  38         a = anp1;
  39     }
  40
  41     // Final cleanup iteration recovers down to ~2ULPs:
  42     return (a + g)/2;
  43 }
  44
  45
  46 }}}
  47 #endif