#5334 inverse_jacobi_sn/cn/dn with bigfloat m and a rational argument lose precision

[David Scherfgen]

This bug was discovered by Claude Opus 5.5 and reproduced on Maxima Git commit c62924c (2026-09-26) with SBCL 2.2.9.

Where: src/ellipt.lisp, SIMP-%INVERSE_JACOBI_SN, SIMP-%INVERSE_JACOBI_CN, SIMP-%INVERSE_JACOBI_DN (def-simplifier inverse_jacobi_sn, _cn, _dn)

The bigfloat branches pass u and m through BIGFLOAT:TO without $bfloat, so a rational u reaches BIGFLOAT:ACOS or the unary SQRT inside BF-RF as a CL rational and is computed in single float. Converting with $bfloat first fixes it.

Steps

1. Verify that the existing behavior is incorrect: run the reproducer below; it gives the result under "Before the fix". Judge for yourself that this is wrong instead of relying on the ticket: where the bug allows it, try some related inputs of your own, or check the result another way, for example numerically or against the documentation.
2. Edit the file src/ellipt.lisp as described under "Fix".
3. Rebuild Maxima.
4. Verify that the behavior is now correct: the reproducer gives the result under "Correct result". Judge for yourself that this is right, not only that it matches: where the bug allows it, try some related inputs of your own, including some that the fix should leave unchanged.
5. Add the test case(s) under "Test cases" to tests/rtest_elliptic.mac, at the end of the file.
6. Run the complete test suite, share tests included (make check, or run_testsuite(share_tests=true); in ./maxima-local), and check that the change breaks no test. Failures that also occur without the change (a run on the unchanged build shows which) are pre-existing and can be ignored.
7. Add the line * \#5334 inverse_jacobi_sn/cn/dn with bigfloat m and a rational argument lose precision to ChangeLog, under "Bug fixes for numbered bugs:".
8. Commit the fix; a proposed Git commit message is under "Commit message".
9. Push the commit. Then set this ticket's status to "closed" and add a comment naming the commit, for example Fixed by commit [1a2b3c]., with the first six characters of the commit's SHA in the brackets (git rev-parse HEAD | cut -c1-6 prints them). SourceForge links this to the commit.

Reproducer

fpprec: 40$ inverse_jacobi_cn(1/2, 0.5b0);

Before the fix:

1.142429144434984062103964282253040692846b0

Correct result:

1.142429058045777255501395526626045782232b0

Fix

In src/ellipt.lisp (SIMP-%INVERSE_JACOBI_SN, bigfloat-numerical-eval-p branch after "(def-simplifier inverse_jacobi_sn (u m)"), replace

          ((bigfloat-numerical-eval-p u m)
           (let ((uu (bigfloat:to u))
                 (mm (bigfloat:to m)))
             (to (bigfloat:* uu

with

          ((bigfloat-numerical-eval-p u m)
           (let ((uu (bigfloat:to ($bfloat u)))
                 (mm (bigfloat:to ($bfloat m))))
             (to (bigfloat:* uu

In src/ellipt.lisp (SIMP-%INVERSE_JACOBI_SN, complex-bigfloat-numerical-eval-p branch after "(def-simplifier inverse_jacobi_sn (u m)"), replace

               args
             (let ((uu (bigfloat:to u))
                   (mm (bigfloat:to m)))
             (to (bigfloat:* uu

with

               args
             (let ((uu (bigfloat:to ($bfloat u)))
                   (mm (bigfloat:to ($bfloat m))))
             (to (bigfloat:* uu

In src/ellipt.lisp (SIMP-%INVERSE_JACOBI_CN, bigfloat-numerical-eval-p branch after "(def-simplifier inverse_jacobi_cn (u m)"), replace

           (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to u))
                                        (bigfloat:to m))))

with

           (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to ($bfloat u)))
                                        (bigfloat:to ($bfloat m)))))

In src/ellipt.lisp (SIMP-%INVERSE_JACOBI_CN, complex-bigfloat-numerical-eval-p branch after "(def-simplifier inverse_jacobi_cn (u m)"), replace

             (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to u))
                                          (bigfloat:to m)))))

with

             (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to ($bfloat u)))
                                          (bigfloat:to ($bfloat m))))))

In src/ellipt.lisp (SIMP-%INVERSE_JACOBI_DN, bigfloat-numerical-eval-p branch after "(def-simplifier inverse_jacobi_dn (u m)"), replace

          ((bigfloat-numerical-eval-p u m)
           (let ((uu (bigfloat:to u))
                 (mm (bigfloat:to m)))
             (to (bigfloat::bf-inverse-jacobi-dn uu mm))))

with

          ((bigfloat-numerical-eval-p u m)
           (let ((uu (bigfloat:to ($bfloat u)))
                 (mm (bigfloat:to ($bfloat m))))
             (to (bigfloat::bf-inverse-jacobi-dn uu mm))))

In src/ellipt.lisp (SIMP-%INVERSE_JACOBI_DN, complex-bigfloat-numerical-eval-p branch after "(def-simplifier inverse_jacobi_dn (u m)"), replace

             (to (bigfloat::bf-inverse-jacobi-dn (bigfloat:to u) (bigfloat:to m)))))

with

             (to (bigfloat::bf-inverse-jacobi-dn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m))))))

The source file indents with tabs (8 columns); the code above shows them as spaces.

Test cases

Add these lines to tests/rtest_elliptic.mac, separated from the tests before them by two empty lines:

/* Bug #5334: "inverse_jacobi_sn/cn/dn with bigfloat m and a rational argument lose precision" */

block([fpprec : 40],
  [is(abs(inverse_jacobi_sn(1/2, bfloat(1/2))

          - inverse_jacobi_sn(bfloat(1/2), bfloat(1/2))) < 1b-38),
   is(abs(inverse_jacobi_cn(1/2, bfloat(1/2))
          - inverse_jacobi_cn(bfloat(1/2), bfloat(1/2))) < 1b-38),
   is(abs(inverse_jacobi_dn(1/2, bfloat(1/2))
          - inverse_jacobi_dn(bfloat(1/2), bfloat(1/2))) < 1b-38)]);
[true, true, true]$

Commit message

ellipt.lisp: fix inverse_jacobi_sn/cn/dn bigfloat precision

The bigfloat branches of SIMP-%INVERSE_JACOBI_SN,
SIMP-%INVERSE_JACOBI_CN and SIMP-%INVERSE_JACOBI_DN passed u and m to
BIGFLOAT:TO without $BFLOAT, so a rational u reached BIGFLOAT:ACOS or
BF-RF as a CL rational and was computed in single float. The arguments
now go through $BFLOAT first.

This fixes bug #5334.

Found and fixed by Claude Opus 5.5.

[Raymond Toy]

There are other issues with these functions. I'm looking into it and Claude is helping out.

[Raymond Toy]

I don't think it's trivial anymore. Removing that label.

[Raymond Toy]

Some of the issues:

1. Cancellation in the $R_F$ arguments. sn formed $1-u^2$ and
   $1-mu^2$ directly. cn went through $F(\arccos u \mid m)$, which
   rebuilds $1-m\sin^2\varphi$ and loses the $mu^2$ term relative to
   $1-m$. dn went through acn(w,1/m)/sqrt(m), which forms $1/m$.
   Errors reached $10^{11}$ to $10^{12}$ ulps for $u$ near $\pm 1$, for
   $m$ near 1, and for small $u$ with $m$ near 1. This held for real and
   complex arguments, in float and in bigfloat.

2. elliptic_f cancellation for $m < 0$. The A&S 17.4.17 branch
   cancels, so cn lost accuracy for $u$ near 1. The patch bypasses
   this for sn, cn and dn, but elliptic_f itself still has it.

3. Overflow for large arguments. Double float signalled
   FLOATING-POINT-OVERFLOW once $\lvert u\rvert \gtrsim 10^{154}$, or
   once $\lvert m u^2\rvert$ was very large.

4. Very small $m$ in dn. It overflowed inside bf-rf for
   $\lvert m\rvert \sim 10^{-300}$, and gave a wrong bigfloat value for
   small imaginary $w$.

5. Lost imaginary part of $1-u^2$. For tiny complex $u$, both
   factors of $(1-u)(1+u)$ round to have real part exactly 1, so the
   product loses its tiny imaginary part. That part decides which side
   of an $R_F$ cut the result is on.

6. Inconsistent sides of the cut in double float. SBCL's
   real + complex addition turns an imaginary part of $-0$ into $+0$, so
   $R_F$ arguments that are equal in exact arithmetic could land on
   opposite sides of the cut. Float and bigfloat could then disagree.

7. Float and bigfloat used different code, so they could pick
   different branches at points on a cut. dn on
   $\mathrm{Re}\thinspace w = 0$ is one example.

8. Wrong-sign results in cd, dc, sd and ds.

   • cd and dc drop the sign of $u$.
   • sd and ds are wrong when $1+mu^2 < 0$.
   • Other conjugated results come from the rectform pipeline.

[Raymond Toy]

I think these are, for the most part, fixed except for elliptic_f mentioned above. The other inverse functions need to be fixed still. Probably best in a different bug and/or branch.

[Raymond Toy]

Fixed in [9724f8] and related commits by merging in the branch bug-5334-inverse-jacobi-inaccurate

This only handles jacobi sn, cn, and dn. The other issues will be handled in a separate bug. This is a big enough change.

Also added doc/info/impl-notes/inverse-jacobi-accuracy.texi.m4 that describes the algorithms in more detail.