Numerical issues with jacobi_sn



>>>>> "Raymond" == Raymond Toy <toy.raymond at gmail.com> writes:

    Raymond> Just a note that jacobi_sn has some numerical issues for complex
    Raymond> arguments.  I'm not sure why.

    Raymond> For example, with fpprec = 25,

    Raymond> jacobi_sn(1b0+%i*1b0, .7b0) ->
    Raymond>  3.522523469224944528936666b-1 %i + 1.134045971912365284394387b0

    Raymond> This uses the new implementation using AGM.  The old version using
    Raymond> Landen's descending transformation produces:

Some further experiments indicate that AGM works fine for real
arguments.  It just loses accuracy for complex arguments, and quite a
few digits are lost.  The descending Landen transform for sn appears
to work fine for both real and complex arguments.

But the ascending and descending Landen transformations (A&S 16.14,
16.12) don't work very well for cn and dn when the argument gets
larger than 10 or 100, so we can't use that.

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.

Finally, we need cn.  This can be obtained from the ascending Landen
transformation since it only involves dn, which we've already
computed.

I think I need to do some more experiments, but this seems to be a
reasonable approach to get accurate values for the Jacobian elliptic
functions.

I'm kind of curious why AGM doesn't seem to produce very accurate
results for complex arguments.  Is it a bug in Maxima's implementation
of the trig functions?   Something inherent in AGM?  Something else?

Ray