Hello,
I have to use the arithmetic-geometric mean and ...I don't know what to think anymore.
Instead of using the a[i], b[i] approach, which gets seriously slower as [i] goes higher, I used the elliptic integral formula:
agm(x,y):=%pi/4*(x+y)/elliptic_kc( ((x-y)/(x+y))^2 );
which, for numbers not too different in terms of order, it is accurate. However, when x=1 and y=1e-15 were involved (number signifying the largest difference needed), there seemed to be a diference in the results. I tried the first method and it works. Thankfully, the fast convergence only needs 7 steps for the following precision.
Here is what Maxima shows:
M(x,y):=%pi/4*((x+y)/(elliptic_kc( ((x-y)/(x+y))^2 )))$
N:7$
a[0](x,y):=(x+y)/2$ g[0](x,y):=sqrt(x*y)$
a[i](x,y):=(a[i-1](x,y)+g[i-1](x,y))/2$ g[i](x,y):=sqrt(a[i-1](x,y)*g[i-1](x,y))$
for i:2 thru N do a[i](x,y):=(a[i-1](x,y)+g[i-1](x,y))/2$
AGM(x,y):=a[N](x,y)$
m:1$ n:1e-15$
float(M(m,n)); float(AGM(m,n));
(%o138) .04378916557404113
(%o139) .04372423752376933
float(K( ((1-1e-15)/(1+1e-15))^2 ));
(%o149) 17.93590156609527
For confirmation, I tried the following in Scilab:
-->x=1; y=1e-15;
-->%pi*(x+y)/4/%k( ((x-y)/(x+y))^2)
ans =
0.0437892
-->%k( ((x-y)/(x+y))^2)
ans =
17.935902
and in Mathematica (online):
http://www.wolframalpha.com/input/?i=elliptick(+((1-1e-15)/(1%2B1e-15))^2)
http://www.wolframalpha.com/input/?i=pi/4*(1%2B1e-15)/(elliptick(((1-1e-15)/(1%2B1e-15))^2)
I am confused: which one is accurate? This is not a hunt for bugs, but an accidental, unfortunate discovery.
Regards,
Vlad.