agm.hpp 1.1 KB

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