On 6/11/10 4:54 PM, Je Suis wrote:
> 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:
>
>
[snip]
> 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
>
>
Think about what you're calculating here. x=1 and y is small.
((x-y)/(x+y))^2 is very nearly one. Note also that elliptic_kc has a
pole at 1, so any tiny error in the argument of elliptic_kc is greatly
magnified.
> 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.
>
We can work around some of the problems and produce an answer close to
what Mathematica produces by using bfloats.
fpprec:32;
x=1;
y=1b-15;
elliptic_kc(((x-y)/(x+y))^2);
17.9265353780153059...
Ray