48 lines
1.1 KiB
C++
48 lines
1.1 KiB
C++
|
// (C) Copyright Nick Thompson 2020.
|
||
|
|
// Use, modification and distribution are subject to the
|
||
|
|
// Boost Software License, Version 1.0. (See accompanying file
|
||
|
|
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||
|
|
|
||
|
|
#ifndef BOOST_MATH_TOOLS_AGM_HPP
|
||
|
|
#define BOOST_MATH_TOOLS_AGM_HPP
|
||
|
|
#include <limits>
|
||
|
|
#include <cmath>
|
||
|
|
|
||
|
|
namespace boost { namespace math { namespace tools {
|
||
|
|
|
||
|
|
template<typename Real>
|
||
|
|
Real agm(Real a, Real g)
|
||
|
|
{
|
||
|
|
using std::sqrt;
|
||
|
|
|
||
|
|
if (a < g)
|
||
|
|
{
|
||
|
|
// Mathematica, mpfr, and mpmath are all symmetric functions:
|
||
|
|
return agm(g, a);
|
||
|
|
}
|
||
|
|
// Use: M(rx, ry) = rM(x,y)
|
||
|
|
if (a <= 0 || g <= 0) {
|
||
|
|
if (a < 0 || g < 0) {
|
||
|
|
return std::numeric_limits<Real>::quiet_NaN();
|
||
|
|
}
|
||
|
|
return Real(0);
|
||
|
|
}
|
||
|
|
|
||
|
|
// The number of correct digits doubles on each iteration.
|
||
|
|
// Divide by 512 for some leeway:
|
||
|
|
const Real scale = sqrt(std::numeric_limits<Real>::epsilon())/512;
|
||
|
|
while (a-g > scale*g)
|
||
|
|
{
|
||
|
|
Real anp1 = (a + g)/2;
|
||
|
|
g = sqrt(a*g);
|
||
|
|
a = anp1;
|
||
|
|
}
|
||
|
|
|
||
|
|
// Final cleanup iteration recovers down to ~2ULPs:
|
||
|
|
return (a + g)/2;
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
}}}
|
||
|
|
#endif
|