double agm(const double a, const double b)
{
//if (a <= 0 || b <= 0) throw std::invalid_argument("agm(a,b) is defined for a>0, b>0");
if (a > b) return agm(b, a);
double aLast = 0;
double bLast = 0;
double aNew = a;
double bNew = b; // temporary variables for storing a(n), b(n)
while (aLast != aNew)
{
aLast = aNew;
bLast = bNew;
aNew = std::sqrt(aLast * bLast);
bNew = (aLast + bLast) / 2;
}
//#ifndef NDEBUG // Assertion module
// const double eps = 1e-9;
// double a_recursive = calc_a(10, a, b);
// assert(std::abs(1 - a_recursive / aLast) < eps);
//#endif
return aLast;
}