1234567891011121314151617181920212223242526272829303132333435363738394041424344454647 |
- // (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
|