]>
Commit | Line | Data |
---|---|---|
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 <cmath> | |
9 | ||
10 | namespace boost { namespace math { namespace tools { | |
11 | ||
12 | template<typename Real> | |
13 | Real agm(Real a, Real g) | |
14 | { | |
15 | if (a < g) | |
16 | { | |
17 | // Mathematica, mpfr, and mpmath are all symmetric functions: | |
18 | return agm(g, a); | |
19 | } | |
20 | // Use: M(rx, ry) = rM(x,y) | |
21 | if (a <= 0 || g <= 0) { | |
22 | if (a < 0 || g < 0) { | |
23 | return std::numeric_limits<Real>::quiet_NaN(); | |
24 | } | |
25 | return Real(0); | |
26 | } | |
27 | ||
28 | // The number of correct digits doubles on each iteration. | |
29 | // Divide by 512 for some leeway: | |
30 | const Real scale = sqrt(std::numeric_limits<Real>::epsilon())/512; | |
31 | while (a-g > scale*g) | |
32 | { | |
33 | Real anp1 = (a + g)/2; | |
34 | g = sqrt(a*g); | |
35 | a = anp1; | |
36 | } | |
37 | ||
38 | // Final cleanup iteration recovers down to ~2ULPs: | |
39 | return (a + g)/2; | |
40 | } | |
41 | ||
42 | ||
43 | }}} | |
44 | #endif |