#5663 jacobi_sn and jacobi_am lose many digits for larger arguments

[David Scherfgen]

By Claude:

jacobi_sn returns wrong digits for moderately large arguments, with floats and with bigfloats:

(%i1) display2d:false$
(%i2) jacobi_sn(20.0, 2.5);
(%o2) 0.2282499166210549
(%i3) jacobi_sn(20.0*%i, 0.75);
(%o3) -0.2325375293373911*%i
(%i4) jacobi_sn(20.0*%i, 1.75+%i);
(%o4) 0.827993704767369-0.8735275803124387*%i

The correct values are 0.2282498061811287, -0.2325715776527406*%i and 0.8280244930780284-0.8738216619745627*%i (computed with mpmath, and Maxima agrees at fpprec : 90). So the first result is wrong from the 7th digit on, and the others from the 4th or 5th. sin(jacobi_am(z, m)) gives the same wrong values, and jacobi_dn and jacobi_cn are affected too, for example jacobi_dn(20.0*%i, 0.75) gives 1.0200761132938823 instead of 1.020081934968419. With bigfloats at fpprec : 32, jacobi_sn(20b0*%i, 1.75b0+%i) has only about 20 correct digits. The error grows quickly with abs(z).

[Raymond Toy]

It's been a long time since I looked at these algorithms. Not sure where the loss of precision comes from.

[Raymond Toy]

jacobi_am needs lots of work. It was just a start that I never really finished. I had dome some work a few years ago to improve it and also handle complex values, but I didn't finish. All the algorithms I found had some deficiencies and I haven't worked on it since.

[David Scherfgen]

See the attached file for Claude's proposal for a fix.

[Raymond Toy]

That looks good, but I haven't tested it yet. Still has a problem when m is near 1, but I guess there's not much we can do about that.

[Raymond Toy]

Assuming I implemented the suggested changes correctly, jacobi_sn(20.0,2.5) is still a bit inaccurate: 0.2282498060629115. To get the correct answer you show, I need to transform m to 1/m. Then I get the right answer.

In addition, these tests fail now. The accuracy has been reduced.

***************** rtest_elliptic.mac: Problem 236 (line 1083) *****************

Input:
closeto(jacobi_am(5.0b-1, 1.5b0) - 4.7071978970469916313225786086966b-1, 
                                                                    7.7038b-34)


Result:
3.0814879110195773648895647081359b-33

This differed from the expected result:
true

***************** rtest_elliptic.mac: Problem 245 (line 1129) *****************

Input:
makelist(block([z : 2 k, m : 5.0b-1], closeto(sin(jacobi_am(z, m))

                                     - jacobi_sn(z, m), 1.8489b-32)), k, 0, 10)


Result:
[true, true, true, true, true, true, true, true, true, 
                                   2.6192647243666407601561300019155b-32, true]

This differed from the expected result:
[true, true, true, true, true, true, true, true, true, true, true]

***************** rtest_elliptic.mac: Problem 247 (line 1139) *****************

Input:
makelist(block([z : 2 k, m : 1 + 5.0b-1 %i], 
        closeto(sin(jacobi_am(z, m)) - jacobi_sn(z, m), 2.5271b-32)), k, 0, 10)


Result:
[true, true, 2.8947950766216634552301329897868b-32, true, true, 
           2.9070698811219211703284680771896b-32, true, true, true, true, true]

This differed from the expected result:
[true, true, true, true, true, true, true, true, true, true, true]

Note that these tests are rather circular. jacobi_am(z,m) uses asin(jacobi_sn(z,m)) to compute the value for complex args. The differences come from computing sin(asin(z)). Perhaps we need some work on the evaluation of the trig functions and inverses for complex values.

[David Scherfgen]

Claude thinks you didn't apply all of the changes, that would explain your results.
It's currently devising a very precise version that is "essentially optimal". I will post it here when it's ready.

[David Scherfgen]

I saw you made a commit. Maybe you want to compare the results with what Claude proposes below (see also the attachment).

By Claude:

Thanks for testing. 0.2282498060629115 is what you get with the new root-mu formula but the old stopping test (< (abs m) (epsilon u)) in ELLIPTIC-SN-DESCENDING (SBCL gives 0.2282498060629152 for that combination), so I think that hunk didn't make it into your tree. For this argument the recursion ends with imagpart(v) near 14, and stopping at abs(m) < eps drops a term of relative size about m*exp(2*imagpart(v))/16, here 2e-10.

You're right about 1/m, though. For real m > 1 the first descending step gives abs(mu) = 1 and turns a real argument into a complex one, v = z*(1+%i*sqrt(m-1))/2. With A&S 16.11 everything stays real. At z = 20 both ways give about 14 digits, because the rounding error of v gets multiplied by z*cn*dn/sn, about 80 here. But 1/m is more accurate on average, and jacobi_dn(20.0, 2.5) no longer comes out complex.

There is something we can do about m near 1. Two things lose digits there. At the second Landen step 1-mu cancels, because mu is about 1-4*sqrt(1-m). Carrying k = sqrt(1-m) through the recursion avoids that, since sqrt(1-mu) = 2*sqrt(k)/(1+k). And CN used the ascending transformation, where mu rounds to 1 and dn^2-sqrt(mu1) cancels. For example jacobi_cn(20.0, 0.99999999) gives -0.2303 instead of -0.5554. Getting cn and dn from the same descending recursion (A&S 16.12.3 and 16.12.4) fixes that.

The attached handover replaces the earlier one completely and applies to a clean master, so please start from master if parts of the old patch are still in your tree. ELLIPTIC-SN-DESCENDING gets an optional argument to return cn and dn as well. SN is unchanged, CN uses it instead of the ascending transformation, and DN keeps the Gauss transformation for real arguments with m < 8/9 and uses it otherwise. Over 200 random real and complex arguments, all errors are now within 10 times eps times the condition number, with doubles and with bigfloats. It also fixes #4300. With real arguments, jacobi_sn takes at most 9% longer than before, and jacobi_dn with m < 8/9 at most 13%. Complex jacobi_sn takes up to 25% longer. jacobi_cn, and jacobi_dn in the other cases, take up to 1.4 times as long with doubles and up to 1.9 times with bigfloats, because the bottom of the recursion now needs cos as well as sin. I also took the log out of the stopping test. It made bigfloat evaluation up to 10 times slower.

[Raymond Toy]

I did apply the new condition replacing the old |m| < eps. Without that, jacobi_sn(1+100*%i, 1e-5) is way off, as mentioned below.

The new elliptic-sn-descending should use keyword args. I also dislike returning different number of values based on some parameter. It might be better to have elliptic-cn-descending and elliptic-dn-descending. Common code could be factored out. I haven't looked at that yet. (Or maybe made into a macro.)

All the comments I had previously made are gone, so it's harder to follow things. (The original was missing comments too. It's a pain to look up A&S references.)

[Raymond Toy]

For the example jacobi_sn(20,2.5), it looks like we should apply the transformation http://functions.wolfram.com/09.36.16.0024.01 first to get the modulus into the range [-1,1]:

jacobi_sn(sqrt(m)*z, 1/m) = sqrt(m)*jacobi_sn(z,m)

Or

                                        1
                   jacobi_sn(sqrt(m) z, ─)
                                        m
(%o38)             ─────────────────────── = jacobi_sn(z, m)
                           sqrt(m)

Then

subst([m=2.5,z=20],%o38);
(%o42)              0.2282498061811312 = 0.2282499166210549

The LHS if very close to the expected 0.2282498061811287.