Numerical issues with jacobi_sn



Raymond Toy wrote:
> I did some further experiments.  We can compute dn from sn using dn =
> sqrt(1-m*sn) (A&S 16.9.5).  For real arguments, there's no ambiguity
> with sqrt since dn >= 0.  Not sure about complex values, but simple
> experiments indicate that it also works for complex values.
>   
Some more interesting examples.

jacobi_sn(1.785063352082689082581887b0 *%i +
8.890858759528509578661528b-1, 9.434463451695984398149033b-1 -
1.476052973708684178844821b-1 * %i)
-> 1.016169837568035196962361b0 %i + 3.327138718909161319167269b0

This differs from the descending Landen version and from what
Mathematica says the answer should be:

1.345233534539675700312281b0 - 7.599023926615176284214056b-2 %i

And we can't use dn = sqrt(1-m*sn^2).  Use the same parameters above,
sqrt(1-m*sn^2) is
8.617730683333292717095686b-1 %i + 2.663978258141280808361839b-1

But the true answer is the negative of this.

At this point, I'm stumped.  Neither AGM nor Landen/Gauss
transformations produce the correct answers.

Ray