48 lines · plain
1// (C) Copyright Nick Thompson 2020.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_AGM_HPP7#define BOOST_MATH_TOOLS_AGM_HPP8#include <limits>9#include <cmath>10 11namespace boost { namespace math { namespace tools {12 13template<typename Real>14Real 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#endif48