14 Jul
2009
14 Jul
'09
7:15 p.m.
AGM16(a,b)= \\ return AGM(a,b) { local(P,M,R); a = sqrt(a); b = sqrt(b); while ( abs(a-b) > eps, \\ eps == 10^(-precision) P=(a+b)/2; M=(a-b)/2; R=(P^4-M^4)^(1/4); a=(P+R)/2; b=(P*R*(P^2+R^2)/2)^(1/4); ); return( a^2 ); } The rate of convergence is _very_ impressive: ? t = AGM16(1.0, 0.5) [1.00000000000000, 0.292893218813452] \\ == [a, a-b] [0.853460904507906, 2.94235252084151 E-17] [0.853460904507906, 3.16563558120117 E-273] [0.853460904507906, 1.02028230345421 E-4368] [0.853460904507906, 1.38315649004407 E-69896] [0.853460904507906, 0.E-100009] \\ using a precision of 100,000 digits 0.728395515523453 \\ returned value ? t-agm(a,b) -4.08056715461103 E-100010 \\ check