an+1 = (an + bn) / 2
bn+1 = √(an bn)

a0 > b0 とする。

a1 = (a0 + b0) / 2 < (a0 + a0) / 2 = a0
b1 = √(a0 b0) > √(b0 b0) = b0 

a1 - b1 = (a0 ＋b0) / 2 - √(a0 b0) = 1 / 2 (√a0 - √b0) > 0 

帰納的に 

a0 > … > an > an+1 > bn+1 > bn > … > b0 

より数列 {an}, {bn} は単調数列となり同じ値に収束する。 

#include 

#ifndef DBL_EPSILON
#define DBL_EPSILON     2.2204460492503131E-16
#endif

double
_agm(a, b)
    double a;
    double b;
{
    double t;
    while(fabs(a / b - 1.0) > DBL_EPSILON){
        t = (a + b) * 0.5;
        b = sqrt(a * b);
        a = t;
    }
    return(a);
}

• 数学選書 特殊関数入門, 一松　信(ひとつまつ・しん), 森北出版, 1999年

• An implementation of Complete Elliptic Integral of the First Kind for C