math/tools/agm.hpp

0001 //  (C) Copyright Nick Thompson 2020.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_AGM_HPP
0007 #define BOOST_MATH_TOOLS_AGM_HPP
0008 #include <limits>
0009 #include <cmath>
0010
0011 namespace boost { namespace math { namespace tools {
0012
0013 template<typename Real>
0014 Real agm(Real a, Real g)
0015 {
0016     using std::sqrt;
0017
0018     if (a < g)
0019     {
0020         // Mathematica, mpfr, and mpmath are all symmetric functions:
0021         return agm(g, a);
0022     }
0023     // Use: M(rx, ry) = rM(x,y)
0024     if (a <= 0 || g <= 0) {
0025         if (a < 0 || g < 0) {
0026             return std::numeric_limits<Real>::quiet_NaN();
0027         }
0028         return Real(0);
0029     }
0030
0031     // The number of correct digits doubles on each iteration.
0032     // Divide by 512 for some leeway:
0033     const Real scale = sqrt(std::numeric_limits<Real>::epsilon())/512;
0034     while (a-g > scale*g)
0035     {
0036         Real anp1 = (a + g)/2;
0037         g = sqrt(a*g);
0038         a = anp1;
0039     }
0040
0041     // Final cleanup iteration recovers down to ~2ULPs:
0042     return (a + g)/2;
0043 }
0044
0045
0046 }}}
0047 #endif