maxima-commits Mailing List for Maxima -- GPL CAS based on DOE-MACSYMA

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, bug-5663-1-jacobi_sn-inaccurate has been updated
       via  a61d7d256834a872567a2f50d042635a071d58da (commit)
      from  67c22aafcf0afb82e2ab90ec78e874d0f5f0aab8 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit a61d7d256834a872567a2f50d042635a071d58da
Author: Raymond Toy <toy...@gm...>
Date:   Sat Oct 17 18:18:52 2026 -0700

    Fix a few more minor comment typos

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index a35df6005..95a0a93a7 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -93,7 +93,7 @@
   ;; m/4*cos(u)^2 and m/4*sin(u)^2, which grow like
   ;; exp(2*abs(imagpart(u))) relative to sin(u) or cos(u).  The terms
   ;; with u shift the phase by about m*u/4, less than the rounding
-  ;; error in u when abs(m) < eps.  For dn, the error is m/2*sin(u)^2
+  ;; error in u when abs(m) < eps.  For dn, the error is m/2*sin(u)^2,
   ;; which grows like m*exp(2*|y|)/8, which is still covered by the
   ;; test.
   (and (< (abs m) (epsilon u))
@@ -144,7 +144,7 @@
      ;;
      ;;   dn(u|m) = (1-p)/(1+p)
      ;;
-     ;; where p = 1-sqrt(mu)*sn^2(v|mu).
+     ;; where p = sqrt(mu)*sn^2(v|mu).
      (multiple-value-bind (v mu root-mu root-mu1)
 	  (descending-transform u m root-m1)
 	(multiple-value-bind (s c d)

-----------------------------------------------------------------------

Summary of changes:
 src/ellipt.lisp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, bug-5663-1-jacobi_sn-inaccurate has been updated
       via  67c22aafcf0afb82e2ab90ec78e874d0f5f0aab8 (commit)
      from  20399fc616b901333b9885ad9e1dc5fbb2130c4c (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 67c22aafcf0afb82e2ab90ec78e874d0f5f0aab8
Author: Raymond Toy <toy...@gm...>
Date:   Sat Oct 17 18:15:33 2026 -0700

    Clean up code and comments
    
    Fix some typos and redundancies in the comments.
    
    Rename some variables in the last cond sexp in dn so it matches what
    descending-transform returns.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index 6ef3989e6..a35df6005 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -54,28 +54,6 @@
 (declaim (inline descending-transform))
 
 ;; Landen descending transform (Gauss' transform).
-;;
-;; Given u and m, compute:
-;;
-;;   mu = ((1-sqrt(m1))/(1+sqrt(m1)))^2
-;;   v = u/sqrt(mu)
-;;
-;; where m1 = 1 - m.  However, compute sqrt(mu) =
-;; (1-sqrt(m1))/(1+sqrt(m1) and square this to get mu.  Doing it in
-;; the order above loses information when m is complex.
-;;
-;; Also (1-sqrt(m1))/(1+sqrt(m1)) loses accuracy when m is small.
-;; Instead compute the equivalent m/(1+sqrt(m1))^2.  (Multiply the
-;; former by (1+sqrt(m1))/(1+sqrt(m1))).
-;;
-;; Returns 3 values:
-;;
-;;   v
-;;   mu
-;;   sqrt(1-mu)
-;;
-;; However sqrt(1-mu) is computed as 2*sqrt(m1)/(1+sqrt(m1)) to get
-;; rid of the round-off in 1-mu when m is near 1.
 (defun descending-transform (u m root-m1)
   ;; A&S 16.12.1 (https://dlmf.nist.gov/22.7.E1).  ROOT-M1 is
   ;; sqrt(1-m).  Return v = u/(1+sqrt(mu)), mu, sqrt(mu), and
@@ -115,9 +93,9 @@
   ;; m/4*cos(u)^2 and m/4*sin(u)^2, which grow like
   ;; exp(2*abs(imagpart(u))) relative to sin(u) or cos(u).  The terms
   ;; with u shift the phase by about m*u/4, less than the rounding
-  ;; error in u when abs(m) < eps.  The base case for dn, includes the
-  ;; first-order term, so its error is O(m^2*sin(u)^4) and this test
-  ;; is more than enough for it.
+  ;; error in u when abs(m) < eps.  For dn, the error is m/2*sin(u)^2
+  ;; which grows like m*exp(2*|y|)/8, which is still covered by the
+  ;; test.
   (and (< (abs m) (epsilon u))
        (or (realp u)
 	   (zerop m)
@@ -161,6 +139,12 @@
      ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4) for dn:
      ;;
      ;;   dn(u|m) = (dn^2(v|mu) - (1-sqrt(mu)))/((1+sqrt(mu))-dn^2(v|mu))
+     ;;
+     ;; Using dn^2(v|mu) = 1-mu*sn^2(v|mu), we get
+     ;;
+     ;;   dn(u|m) = (1-p)/(1+p)
+     ;;
+     ;; where p = 1-sqrt(mu)*sn^2(v|mu).
      (multiple-value-bind (v mu root-mu root-mu1)
 	  (descending-transform u m root-m1)
 	(multiple-value-bind (s c d)
@@ -211,7 +195,7 @@
 	   (multiple-value-bind (s c d)
 	       (elliptic-sn-cn-dn-descending (* root u) (/ (- m) (- 1 m)) (/ root))
 	     (declare (ignore c))
-	     (/ s (* root d)))))        
+	     (/ s (* root d)))))
 	(t
 	 (elliptic-sn-descending u m (sqrt (- 1 m))))))
 
@@ -230,7 +214,7 @@
 	 ;;   dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
 	 (/ (cosh u)))
 	((and (realp m) (> m 1))
-	 ;; A&S 16.11.3 (https://dlmf.nist.gov/22.17.E3):
+	 ;; A&S 16.11.4 (https://dlmf.nist.gov/22.17.E4):
          ;;
          ;;   dn(u|m) = cn(sqrt(m)*u|1/m)
          (multiple-value-bind (s c d)
@@ -253,15 +237,15 @@
 	((not (and (realp u) (realp m) (< m 8/9)))
 	 ;; dn is computed from one step of the descending transformation,
 	 ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4),
- 	 ;;
+	 ;;
 	 ;;   dn(u|m) = (1-p)/(1+p),  p = sqrt(mu)*sn(v|mu)^2,
-         ;;
+	 ;;
 	 ;; with v = u/(1+sqrt(mu)) and sqrt(mu) =
 	 ;; (1-sqrt(1-m))/(1+sqrt(1-m)).  For real u and 0 <= m < 8/9,
 	 ;; sqrt(mu) < 1/2, so 0 <= p < 1/2 and 1-p cannot cancel; the
 	 ;; branch below then needs only sn.  Otherwise 1-p can cancel,
 	 ;; for example when m is near 1 and u is near an odd multiple of
-	 ;; K(m)
+	 ;; K(m).
          (multiple-value-bind (s c d)
 	     (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))
            (declare (ignore s c))
@@ -270,10 +254,10 @@
 	 ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4), as above.
 	 ;; This is also the Gauss transformation of
 	 ;; http://functions.wolfram.com/09.29.16.0013.01.
-	 (multiple-value-bind (z root-sq root root-1-root-sq)
+	 (multiple-value-bind (v mu root-mu root-mu1)
 	     (descending-transform u m (sqrt (- 1 m)))
-	   (let* ((s (elliptic-sn-descending z root-sq root-1-root-sq))
-		  (p (* root s s)))
+	   (let* ((s (elliptic-sn-descending v mu root-mu1))
+		  (p (* root-mu s s)))
 	     (/ (- 1 p)
 		(+ 1 p)))))))
 
@@ -285,8 +269,8 @@
 	 ;; cn(u,m) = cos(u) + 1/4*m*(u-sin(u)*cos(u))*sin(u)
 	 (cos u))
 	((= m 1)
-	 ;; jacobi_cn(u,1) = sech(u).  Should we use A&S 16.15.3
-	 ;; (https://dlmf.nist.gov/22.10.E9) if m is close enough to
+	 ;; jacobi_cn(u,1) = sech(u).  Should we use A&S 16.15.2
+	 ;; (https://dlmf.nist.gov/22.10.E8) if m is close enough to
 	 ;; 1?
 	 ;;
 	 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)

-----------------------------------------------------------------------

Summary of changes:
 src/ellipt.lisp | 54 +++++++++++++++++++-----------------------------------
 1 file changed, 19 insertions(+), 35 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, bug-5663-1-jacobi_sn-inaccurate has been updated
       via  20399fc616b901333b9885ad9e1dc5fbb2130c4c (commit)
       via  a0998dd9b33e6f95c2c6bbbc7c95bd58b62cb0ab (commit)
      from  277371813ae8635cc124edd1b98a215691cb3de0 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 20399fc616b901333b9885ad9e1dc5fbb2130c4c
Author: Raymond Toy <toy...@gm...>
Date:   Sat Oct 17 15:44:36 2026 -0700

    Add additional tests
    
    These were generated by Claude, including the reference values
    computed by Claude using mpmath.
    
    I updated the `closeto` thresholds to be much more restrictive.
    
    The tests from m = -1e10 shows that A&S 16.10 improves accuracy.
    Without it, we lose 4 or so digits of accuracy.

diff --git a/tests/rtest_elliptic.mac b/tests/rtest_elliptic.mac
index d2498e873..dc8046180 100755
--- a/tests/rtest_elliptic.mac
+++ b/tests/rtest_elliptic.mac
@@ -1225,6 +1225,60 @@ true;
 closeto(jacobi_dn(20.0*%i, 0.75) - 1.020081934968419, 6.6614e-16);
 true;
 
+block([z : 20.0*%i],
+  [closeto(jacobi_dn(z, 0.75) - 1.020081934968419, 6.6614e-16),
+   closeto(jacobi_cn(z, 0.75) - 1.026688627935405, 1.33323e-15),
+   closeto(jacobi_dn(z, 1.75+%i) - (1.0765867842487509 + 1.2123267066947074*%i), 1.0663e-14),
+   closeto(jacobi_cn(z, 1.75+%i) - (-1.2004970668366945 - 0.6027051283045847*%i), 7.1504e-15)]);
+[true, true, true, true]$
+
+block([u : 20.0, m : 2.5],
+  [closeto(jacobi_sn(u, m) - 0.2282498061811287, 4.7740e-15),
+   closeto(jacobi_cn(u, m) - 0.973602601669838, 1.1103e-15),
+   closeto(jacobi_dn(u, m) + 0.9326065970953095, 3.2197e-15),
+   imagpart(jacobi_dn(u, m))]);
+[true, true, true, 0]$
+
+block([u : 20.0, m : 0.99999999],
+  [closeto(jacobi_sn(u, m) - 0.8315901214124292, 5.5512e-16),
+   closeto(jacobi_cn(u, m) + 0.5553898360334489, 7.7716e-16),
+   closeto(jacobi_dn(u, m) - 0.5553898422591851, 1.1103e-15)]);
+[true, true, true]$
+
+block([fpprec : 32, m],
+  m : 1 - 2b0^-60,
+  [closeto(jacobi_cn(20b0, m) - 4069705528989350766860108776551655/10^42, 4.8214e-40),
+   closeto(jacobi_dn(20b0, m) - 4174908960764880545715816189280978/10^42, 4.5918e-40)]);
+[true, true]$
+
+/* m = 1-2^-53, the largest float below 1.  sqrt(mu) must be computed
+ * from sqrt(1-m) alone, or cn and dn lose about 4 more digits.
+ */
+block([u : 20.0, m : 1 - 2.0^-53],
+  [closeto(jacobi_cn(u, m) + 2.6107123887567846e-9, 4.7977e-23),
+   closeto(jacobi_dn(u, m) - 1.0855326878511021e-8, 1.3235e-23)]);
+[true, true]$
+
+/* Large negative m.  Without A&S 16.10, 1+sqrt(mu) cancels and about
+ * log10(sqrt(-m)) digits are lost.  The condition number with respect
+ * to u is about 2e5 here, so the tolerances are relatively large.
+ */
+block([u : 3.0, m : -1e10],
+  [closeto(jacobi_sn(u, m) - 0.11593314977968198, 2.5943e-12),
+   closeto(jacobi_cn(u, m) + 0.99325701849126736, 3.0254e-13),
+   closeto(jacobi_dn(u, m) - 11593.315021096501, 2.5943e-7)]);
+[true, true, true]$
+
+/* Bug #4300: "jacobi_sn(99.0, 1.5) incorrect" */
+
+block([u : 99.0, m : 1.5],
+  [closeto(jacobi_sn(u, m) + 0.37384023685931953, 2.1650e-15),
+   closeto(jacobi_cn(u, m) - 0.9274931144245589, 9.9921e-16),
+   closeto(jacobi_dn(u, m) - 0.8890248680197038, 1.3323e-15),
+   imagpart(jacobi_dn(u, m))]);
+[true, true, true, 0]$
+
+
 (fpprec:oldfpprec,kill(l2,l3,test_deriv,oldfpprec),done);
 done;
 

commit a0998dd9b33e6f95c2c6bbbc7c95bd58b62cb0ab
Author: Raymond Toy <toy...@gm...>
Date:   Sat Oct 17 15:41:27 2026 -0700

    Add code to handle real negative m
    
    We accidentally left this out, but we need to have code to handle
    large real negative m because we lose up to 4 digits of accuracy
    without it.  This is basically A&S 16.10.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index a69865a37..6ef3989e6 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -200,6 +200,18 @@
 	 (let ((root-m (sqrt m)))
 	   (/ (elliptic-sn-descending (* root-m u) (/ m) (sqrt (/ (- m 1) m)))
 	      root-m)))
+	((and (realp m) (minusp m))
+	 ;; A&S 16.10.2 (https://dlmf.nist.gov/22.17.E6):
+         ;;
+         ;;   sn(u|m) = sn(v|mu)/(sqrt(1-m)*dn(v|mu)),
+         ;;
+         ;; where v = sqrt(1-m)*u and mu = -m/(1-m).  sqrt(1-mu) =
+         ;; 1/sqrt(1-m).
+	 (let ((root (sqrt (- 1 m))))
+	   (multiple-value-bind (s c d)
+	       (elliptic-sn-cn-dn-descending (* root u) (/ (- m) (- 1 m)) (/ root))
+	     (declare (ignore c))
+	     (/ s (* root d)))))        
 	(t
 	 (elliptic-sn-descending u m (sqrt (- 1 m))))))
 
@@ -227,6 +239,17 @@
 					   (sqrt (/ (- m 1) m)))
            (declare (ignore s d))
            c))
+	((and (realp m) (minusp m))
+	 ;; A&S 16.10.4 (https://dlmf.nist.gov/22.17.E8):
+         ;;
+         ;;   dn(u|m) = 1/dn(v|mu),
+         ;;
+         ;; where v = sqrt(1-m)*u and mu = -m/(1-m).  sqrt(1-mu) =
+         ;; 1/sqrt(1-m).
+	 (let ((root (sqrt (- 1 m))))
+	   (/ (nth-value 2 (elliptic-sn-cn-dn-descending (* root u)
+							 (/ (- m) (- 1 m))
+							 (/ root))))))
 	((not (and (realp u) (realp m) (< m 8/9)))
 	 ;; dn is computed from one step of the descending transformation,
 	 ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4),
@@ -278,6 +301,18 @@
 					   (sqrt (/ (- m 1) m)))
            (declare (ignore s c))
            d))
+	((and (realp m) (minusp m))
+	 ;; A&S 16.10.3 (https://dlmf.nist.gov/22.17.E7):
+         ;;
+         ;;   cn(u|m) = cn(v|mu)/dn(v|mu),
+         ;;
+         ;; where v = sqrt(1-m)*u and mu = -m/(1-m).  sqrt(1-mu) =
+         ;; 1/sqrt(1-m).
+	 (let ((root (sqrt (- 1 m))))
+	   (multiple-value-bind (s c d)
+	       (elliptic-sn-cn-dn-descending (* root u) (/ (- m) (- 1 m)) (/ root))
+	     (declare (ignore s))
+	     (/ c d))))
 	(t
          (multiple-value-bind (s c d)
              (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))

-----------------------------------------------------------------------

Summary of changes:
 src/ellipt.lisp          | 35 +++++++++++++++++++++++++++++++
 tests/rtest_elliptic.mac | 54 ++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 89 insertions(+)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, bug-5663-1-jacobi_sn-inaccurate has been updated
       via  277371813ae8635cc124edd1b98a215691cb3de0 (commit)
      from  6e996e774a13183ff4b5b4792f116188e22c2560 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 277371813ae8635cc124edd1b98a215691cb3de0
Author: Raymond Toy <toy...@gm...>
Date:   Sat Oct 17 06:58:07 2026 -0700

    Add tests given in bug #5663

diff --git a/tests/rtest_elliptic.mac b/tests/rtest_elliptic.mac
index 94e5d92a4..d2498e873 100755
--- a/tests/rtest_elliptic.mac
+++ b/tests/rtest_elliptic.mac
@@ -1212,6 +1212,18 @@ true;
 closeto('elliptic_e(1.0, 0.5) - 0.9273298836244401, 1e-15);
 true;
 
+/* Bug 5663 */
+closeto(jacobi_sn(20.0, 2.5) - 0.2282498061811287, 4.7740e-15);
+true;
+
+closeto(jacobi_sn(20.0*%i, 0.75) - (-0.2325715776527406*%i), 3.8858e-15);
+true;
+
+closeto(jacobi_sn(20.0*%i, 1.75 + %i) - (0.8280244930780284-0.8738216619745627*%i), 8.0428e-15);
+true;
+
+closeto(jacobi_dn(20.0*%i, 0.75) - 1.020081934968419, 6.6614e-16);
+true;
 
 (fpprec:oldfpprec,kill(l2,l3,test_deriv,oldfpprec),done);
 done;

-----------------------------------------------------------------------

Summary of changes:
 tests/rtest_elliptic.mac | 12 ++++++++++++
 1 file changed, 12 insertions(+)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, bug-5663-1-jacobi_sn-inaccurate has been updated
       via  6e996e774a13183ff4b5b4792f116188e22c2560 (commit)
       via  0b8a62d727287802ec6672e8982a054edec3f739 (commit)
      from  9ce370cf45bc8304e8739a46c0d02a434b15f2b2 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 6e996e774a13183ff4b5b4792f116188e22c2560
Author: Raymond Toy <toy...@gm...>
Date:   Fri Oct 16 18:27:51 2026 -0700

    More comment clean-ups and DLMF links
    
    Also replace `nth-value` calls with `multiple-value-bind`.  This makes
    it a bit easier to see quickly which value we really want from
    elliptic-sn-cn-dn-descending.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index a0b1cd565..a69865a37 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -77,8 +77,9 @@
 ;; However sqrt(1-mu) is computed as 2*sqrt(m1)/(1+sqrt(m1)) to get
 ;; rid of the round-off in 1-mu when m is near 1.
 (defun descending-transform (u m root-m1)
-  ;; A&S 16.12.1.  ROOT-M1 is sqrt(1-m).  Return v = u/(1+sqrt(mu)),
-  ;; mu, sqrt(mu), and sqrt(1-mu).
+  ;; A&S 16.12.1 (https://dlmf.nist.gov/22.7.E1).  ROOT-M1 is
+  ;; sqrt(1-m).  Return v = u/(1+sqrt(mu)), mu, sqrt(mu), and
+  ;; sqrt(1-mu).
   ;;
   ;; Don't calculate mu first, as given in 16.12.1.  We should
   ;; calculate sqrt(mu) = (1-sqrt(m1))/(1+sqrt(m1)), and then compute
@@ -102,7 +103,9 @@
 
 (defun descending-base-p (u m)
   ;; True if m is small enough to stop the descending transformation.
-  ;; By A&S 16.13.1, 16.13.2, and 16.13.3
+  ;; By A&S 16.13.1 (https://dlmf.nist.gov/22.10.E4), 16.13.2
+  ;; (https://dlmf.nist.gov/22.10.E5), and 16.13.3
+  ;; (https://dlmf.nist.gov/22.10.E6):
   ;;
   ;;   sn(u|m) = sin(u) - m/4*(u-sin(u)*cos(u))*cos(u) + O(m^2)
   ;;   cn(u|m) = cos(u) + m/4*(u-sin(u)*cos(u))*sin(u) + O(m^2)
@@ -126,36 +129,36 @@
 (defun elliptic-sn-descending (u m root-m1)
   (cond
     ((descending-base-p u m)
-     ;; A&S 16.13.1
+     ;; A&S 16.13.1 (https://dlmf.nist.gov/22.10.E4)
      (sin u))
     (t
-     ;; A&S 16.12.2:
+     ;; A&S 16.12.2 (https://dlmf.nist.gov/22.7.E2):
      ;;
      ;;   sn(u|m) = (1+sqrt(mu))*sn(v|mu)/(1+sqrt(mu)*sn^2(v|mu)).
      (multiple-value-bind (v mu root-mu root-mu1)
 	 (descending-transform u m root-m1)
        (let ((s (elliptic-sn-descending v mu root-mu1)))
-	 ;; A&S 16.12.2
 	 (/ (* (1+ root-mu) s)
-	    (1+ (* root-mu s s)))))))
+	    (1+ (* root-mu s s))))))))
 
 ;; sn(u|m), cn(u|m), and dn(u|m) for m /= 1.  ROOT-M1 is sqrt(1-m).
 (defun elliptic-sn-cn-dn-descending (u m root-m1)
   (cond
     ((descending-base-p u m)
-     ;; A&S 16.13.1 to 16.13.3, dropping the linear (and higher) terms
-     ;; in m.
+     ;; A&S 16.13.1 (https://dlmf.nist.gov/22.10.E4) to 16.13.3
+     ;; (https://dlmf.nist.gov/22.10.E6), dropping the linear (and
+     ;; higher) terms in m.
      (values (sin u) (cos u) 1))
     (t
-     ;; A&S 16.12.2:
+     ;; A&S 16.12.2 (https://dlmf.nist.gov/22.7.E2):
      ;;
      ;;   sn(u|m) = (1+sqrt(mu))*sn(v|mu)/(1+sqrt(mu)*sn^2(v|mu)).
      ;;
-     ;; A&S 16.12.3 for cn:
+     ;; A&S 16.12.3 (https://dlmf.nist.gov/22.7.E3) for cn:
      ;;
      ;;   cn(u|m) = cn(v|mu)*dn(v|mu)/(1+sqrt(mu)*sn^2(v|mu))
      ;;
-     ;; A&S 16.12.4 for dn:
+     ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4) for dn:
      ;;
      ;;   dn(u|m) = (dn^2(v|mu) - (1-sqrt(mu)))/((1+sqrt(mu))-dn^2(v|mu))
      (multiple-value-bind (v mu root-mu root-mu1)
@@ -202,26 +205,48 @@
 
 (defun dn (u m)
   (cond ((zerop m)
-	 ;; jacobi_dn(u,0) = 1.  Should we use A&S 16.13.3 for small m?
+	 ;; jacobi_dn(u,0) = 1.  Should we use A&S 16.13.3
+	 ;; (https://dlmf.nist.gov/22.10.E6) for small m?
 	 ;;
-	 ;; dn(u,m) = 1 - 1/2*m*sin(u)^2
+	 ;;   dn(u,m) = 1 - 1/2*m*sin(u)^2
 	 1)
 	((= m 1)
-	 ;; jacobi_dn(u,1) = sech(u).  Should we use A&S 16.15.3 if m
-	 ;; is close enough to 1?
+	 ;; jacobi_dn(u,1) = sech(u).  Should we use A&S 16.15.3
+	 ;; (https://dlmf.nist.gov/22.10.E9) if m is close enough to
+	 ;; 1?
 	 ;;
-	 ;; dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
+	 ;;   dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
 	 (/ (cosh u)))
 	((and (realp m) (> m 1))
-	 ;; A&S 16.11.3: dn(u|m) = cn(sqrt(m)*u|1/m)
-	 (nth-value 1 (elliptic-sn-cn-dn-descending (* (sqrt m) u) (/ m)
-						    (sqrt (/ (- m 1) m)))))
+	 ;; A&S 16.11.3 (https://dlmf.nist.gov/22.17.E3):
+         ;;
+         ;;   dn(u|m) = cn(sqrt(m)*u|1/m)
+         (multiple-value-bind (s c d)
+             (elliptic-sn-cn-dn-descending (* (sqrt m) u)
+                                           (/ m)
+					   (sqrt (/ (- m 1) m)))
+           (declare (ignore s d))
+           c))
 	((not (and (realp u) (realp m) (< m 8/9)))
-	 ;; The Gauss transformation below can lose digits through the
-	 ;; cancellation in 1-p unless u and m are real and m < 8/9.
-	 (nth-value 2 (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))))
+	 ;; dn is computed from one step of the descending transformation,
+	 ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4),
+ 	 ;;
+	 ;;   dn(u|m) = (1-p)/(1+p),  p = sqrt(mu)*sn(v|mu)^2,
+         ;;
+	 ;; with v = u/(1+sqrt(mu)) and sqrt(mu) =
+	 ;; (1-sqrt(1-m))/(1+sqrt(1-m)).  For real u and 0 <= m < 8/9,
+	 ;; sqrt(mu) < 1/2, so 0 <= p < 1/2 and 1-p cannot cancel; the
+	 ;; branch below then needs only sn.  Otherwise 1-p can cancel,
+	 ;; for example when m is near 1 and u is near an odd multiple of
+	 ;; K(m)
+         (multiple-value-bind (s c d)
+	     (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))
+           (declare (ignore s c))
+           d))
 	(t
-         ;; 
+	 ;; A&S 16.12.4 (https://dlmf.nist.gov/22.7.E4), as above.
+	 ;; This is also the Gauss transformation of
+	 ;; http://functions.wolfram.com/09.29.16.0013.01.
 	 (multiple-value-bind (z root-sq root root-1-root-sq)
 	     (descending-transform u m (sqrt (- 1 m)))
 	   (let* ((s (elliptic-sn-descending z root-sq root-1-root-sq))
@@ -231,23 +256,33 @@
 
 (defun cn (u m)
   (cond ((zerop m)
-	 ;; jacobi_cn(u,0) = cos(u).  Should we use A&S 16.13.2 for
-	 ;; small m?
+	 ;; jacobi_cn(u,0) = cos(u).  Should we use A&S 16.13.2
+	 ;; (https://dlmf.nist.gov/22.10.E5) for small m?
 	 ;;
 	 ;; cn(u,m) = cos(u) + 1/4*m*(u-sin(u)*cos(u))*sin(u)
 	 (cos u))
 	((= m 1)
-	 ;; jacobi_cn(u,1) = sech(u).  Should we use A&S 16.15.3 if m
-	 ;; is close enough to 1?
+	 ;; jacobi_cn(u,1) = sech(u).  Should we use A&S 16.15.3
+	 ;; (https://dlmf.nist.gov/22.10.E9) if m is close enough to
+	 ;; 1?
 	 ;;
 	 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)
 	 (/ (cosh u)))
 	((and (realp m) (> m 1))
-	 ;; A&S 16.11.2: cn(u|m) = dn(sqrt(m)*u|1/m)
-	 (nth-value 2 (elliptic-sn-cn-dn-descending (* (sqrt m) u) (/ m)
-						    (sqrt (/ (- m 1) m)))))
+	 ;; A&S 16.11.3 (https://dlmf.nist.gov/22.17.E3):
+         ;;
+         ;;   cn(u|m) = dn(sqrt(m)*u|1/m)
+         (multiple-value-bind (s c d)
+             (elliptic-sn-cn-dn-descending (* (sqrt m) u)
+                                           (/ m)
+					   (sqrt (/ (- m 1) m)))
+           (declare (ignore s c))
+           d))
 	(t
-	 (nth-value 1 (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))))))
+         (multiple-value-bind (s c d)
+             (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))
+           (declare (ignore s d))
+           c))))
 
 
 ;; Arithmetic-Geometric Mean algorithm for real or complex numbers.

commit 0b8a62d727287802ec6672e8982a054edec3f739
Author: Raymond Toy <toy...@gm...>
Date:   Fri Oct 16 18:25:56 2026 -0700

    Adjust thresholds for rtest_elliptic tests 237 and 247
    
    Increase the thresholds for these tests just slightly so that they
    pass.  The failures were:
    
    ```
    Running tests in rtest_elliptic:
    ***************** rtest_elliptic.mac: Problem 237 (line 1086) *****************
    
    Input:
    closeto(jacobi_am(1.5b0, 1.5b0 + %i) - (9.3405421687007830327435083030859b-1
                            - 3.7239604521460716553950026072423b-1 %i), 2.3112b-32)
    
    Result:
    2.3162460526365255507658857966452b-32
    
    This differed from the expected result:
    true
    
    ***************** rtest_elliptic.mac: Problem 247 (line 1152) *****************
    
    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, true, true, true, true, true, true, true, true,
                                             2.6328244252901964229102084718221b-32]
    
    This differed from the expected result:
    [true, true, true, true, true, true, true, true, true, true, true]
    
    ```

diff --git a/tests/rtest_elliptic.mac b/tests/rtest_elliptic.mac
index 62ede1b4d..94e5d92a4 100755
--- a/tests/rtest_elliptic.mac
+++ b/tests/rtest_elliptic.mac
@@ -1083,7 +1083,7 @@ true$
 closeto(jacobi_am(0.5b0, 1.5b0) - 0.470719789704699163132257860869655024b0, 7.7038b-34);
 true$
 
-closeto(jacobi_am(1.5b0, 1.5b0+%i) - (0.934054216870078303274350830308588452b0 - 0.372396045214607165539500260724230487b0*%i), 2.3112b-32);
+closeto(jacobi_am(1.5b0, 1.5b0+%i) - (0.934054216870078303274350830308588452b0 - 0.372396045214607165539500260724230487b0*%i), 2.3163b-32);
 true$
 
 closeto(jacobi_am(2.0b0-.5b0*%i, .5b0+.5b0*%i) - (1.46937336079225563334301589326200829b0 - 0.72457112357459480134743514393673164b0*%i), 1.2422b-32);
@@ -1148,7 +1148,7 @@ makelist(block([z : 2*k, m : 2.5b0],
 [true, true, true, true, true, true, true, true, true, true, true];
 
 makelist(block([z : 2*k, m : 1+.5b0*%i],
-    closeto(sin(jacobi_am(z, m))-jacobi_sn(z, m), 2.5271b-32)),
+    closeto(sin(jacobi_am(z, m))-jacobi_sn(z, m), 2.6329b-32)),
   k, 0, 10);
 [true, true, true, true, true, true, true, true, true, true, true];
 

-----------------------------------------------------------------------

Summary of changes:
 src/ellipt.lisp          | 99 ++++++++++++++++++++++++++++++++----------------
 tests/rtest_elliptic.mac |  4 +-
 2 files changed, 69 insertions(+), 34 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  f9ea46acc39004ba8394d2a0eeaed53b04d28a52 (commit)
      from  a7f328122815ef2e700626619dcbe2a4016e1a51 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit f9ea46acc39004ba8394d2a0eeaed53b04d28a52
Author: Dan Gildea <dgildea>
Date:   Fri Oct 16 17:36:30 2026 -0400

    rtest_nfloat.mac
    
    increase accuracy of expected result
    test 25 now passes in cmucl, sbcl, ecl, gcl

diff --git a/src/testsuite.lisp b/src/testsuite.lisp
index 5069a9360..2c83570b8 100644
--- a/src/testsuite.lisp
+++ b/src/testsuite.lisp
@@ -158,8 +158,7 @@
 	;; The tests that failed with abcl 1.5.0
 	((mlist simp) "rtest_hg"
 	 #+abcl ((mlist simp) 120))
-	((mlist simp) "rtest_nfloat"
-	 #-gcl((mlist simp) 25))
+	((mlist simp) "rtest_nfloat")
 	((mlist simp) "rtest_ilt")
 	((mlist simp) "ulp_tests"
 	 ;; Clisp doesn't have denormals
diff --git a/tests/rtest_nfloat.mac b/tests/rtest_nfloat.mac
index c483e108c..8db1dad1f 100644
--- a/tests/rtest_nfloat.mac
+++ b/tests/rtest_nfloat.mac
@@ -87,7 +87,7 @@ nfloat([a = b, a < b, a <= b, a > b, a >- b],[a = 2, b = 4],10);
 [2.0 = 4.0, 2.0 < 4.0, 2.0 <= 4.0, 2.0 > 4.0, 2.0 > - 4.0]$
 
 nfloat(map('rhs, solve(x^3 + x + 1,x)),[], 10);
-[0.341163901914009 - 1.161541399997251 * %i, 1.161541399997251 * %i + 0.341163901914009, - 0.682327803828018]$
+[0.3411639019140098 - 1.161541399997252 * %i, 1.161541399997252 * %i + 0.3411639019140098, - 0.6823278038280196]$
 
 nfloat([acos(x), asin(x), atan(x), acosh(x), asinh(x), atanh(x)], [x=0.12],10);
 [1.4505064444001086, 0.12028988239478806, 0.11942892601833845, 1.4505064444001083 *%i,

-----------------------------------------------------------------------

Summary of changes:
 src/testsuite.lisp     | 3 +--
 tests/rtest_nfloat.mac | 2 +-
 2 files changed, 2 insertions(+), 3 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, bug-5663-1-jacobi_sn-inaccurate has been updated
       via  9ce370cf45bc8304e8739a46c0d02a434b15f2b2 (commit)
       via  da8fa6804f1bcd82abad1e85c4c6612b72c28649 (commit)
       via  880b9be8889c507af2b5ae04effd57012acac976 (commit)
       via  fd2dd75c88ca50e30d7b256e0cc78b5da7384207 (commit)
       via  92979cf7b51ae21222022e5f6fb8d515d50d7619 (commit)
       via  3966a3f14d527b668c36b2d93ec1913a2783141b (commit)
       via  e0e8329ddd727ad8aca560cdccdf78956714b05a (commit)
       via  ab431db5bc95475aa05776daeae4ab344d16c215 (commit)
       via  038ac188a0addb320d2c3f0e7db1982db61769e6 (commit)
       via  7999b7670542ad0b6b56764eabf9e160c00442f5 (commit)
       via  6f25daa244dc0c665f81b5004d71acc9aab6348e (commit)
       via  c497c1cf55b076499c53f3ef5a44e608473b659f (commit)
       via  b0b8d67a05bca2ac348a9dfcdb782693b5567f64 (commit)
       via  3602e648167882b995f15ef824ca392fa8257a92 (commit)
       via  c634df5c157e7fd2f9576655ffae316f3b538afd (commit)
       via  c05cb1b1f728d0c18a617bcefbbe5d881f72c28b (commit)
       via  f01a0209053ab2ebf1f59c4d73e08cd078ccc23d (commit)
       via  bd0f63ebe5a2bc2863f3a6b7aba22a17d61c2a6d (commit)
       via  5e015cf5c094021c1ec5a99854ef3ed34eb9844f (commit)
       via  8606f3175bd7bef36f93a622457594b0e714b35d (commit)
       via  eec2fa78019a654daca57434ca9c6428ea5ee35f (commit)
       via  4fa98b94a5ff3f7753b12260e6d88a0df3f1054a (commit)
       via  3d4c7d2af299a0caac29c68aaa479064fe203aa6 (commit)
       via  317711cba47bcd47ca52990a54ec43558f2a3b52 (commit)
      from  28e0b35b2a35a60991a5a845893b4ea536530432 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 9ce370cf45bc8304e8739a46c0d02a434b15f2b2
Author: Raymond Toy <toy...@gm...>
Date:   Fri Oct 16 08:05:06 2026 -0700

    Clean up comments and add DLMF links
    
    Clean up comments and hopefully make it clearer what is actually being
    done.  And add DLMF links that correspond to the A&S formulas.
    
    It's unfortunate that the DLMF uses k instead of m.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index 920d9b33d..a0b1cd565 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -102,32 +102,63 @@
 
 (defun descending-base-p (u m)
   ;; True if m is small enough to stop the descending transformation.
-  ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
-  ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
+  ;; By A&S 16.13.1, 16.13.2, and 16.13.3
+  ;;
+  ;;   sn(u|m) = sin(u) - m/4*(u-sin(u)*cos(u))*cos(u) + O(m^2)
+  ;;   cn(u|m) = cos(u) + m/4*(u-sin(u)*cos(u))*sin(u) + O(m^2)
+  ;;   dn(u|m) = 1 - 1/2*m*sin^2(u) + O(m^2)
+  ;;
+  ;; Relative to sin(u) and cos(u), the terms with sin(u)*cos(u) are
+  ;; m/4*cos(u)^2 and m/4*sin(u)^2, which grow like
+  ;; exp(2*abs(imagpart(u))) relative to sin(u) or cos(u).  The terms
+  ;; with u shift the phase by about m*u/4, less than the rounding
+  ;; error in u when abs(m) < eps.  The base case for dn, includes the
+  ;; first-order term, so its error is O(m^2*sin(u)^4) and this test
+  ;; is more than enough for it.
   (and (< (abs m) (epsilon u))
        (or (realp u)
 	   (zerop m)
-	   (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u)))))))))
+	   (< (abs m)
+              (* (epsilon u)
+                 (exp (* -2 (abs (imagpart u)))))))))
 
 ;; sn(u|m) for m /= 1.  ROOT-M1 is sqrt(1-m).
 (defun elliptic-sn-descending (u m root-m1)
-  (if (descending-base-p u m)
-      ;; A&S 16.13.1
-      (sin u)
-      (multiple-value-bind (v mu root-mu root-mu1)
-	  (descending-transform u m root-m1)
-	(let ((s (elliptic-sn-descending v mu root-mu1)))
-	  ;; A&S 16.12.2
-	  (/ (* (1+ root-mu) s)
-	     (1+ (* root-mu s s)))))))
+  (cond
+    ((descending-base-p u m)
+     ;; A&S 16.13.1
+     (sin u))
+    (t
+     ;; A&S 16.12.2:
+     ;;
+     ;;   sn(u|m) = (1+sqrt(mu))*sn(v|mu)/(1+sqrt(mu)*sn^2(v|mu)).
+     (multiple-value-bind (v mu root-mu root-mu1)
+	 (descending-transform u m root-m1)
+       (let ((s (elliptic-sn-descending v mu root-mu1)))
+	 ;; A&S 16.12.2
+	 (/ (* (1+ root-mu) s)
+	    (1+ (* root-mu s s)))))))
 
 ;; sn(u|m), cn(u|m), and dn(u|m) for m /= 1.  ROOT-M1 is sqrt(1-m).
 (defun elliptic-sn-cn-dn-descending (u m root-m1)
-  (if (descending-base-p u m)
-      ;; A&S 16.13.1 to 16.13.3
-      (let ((s (sin u)))
-	(values s (cos u) (- 1 (* m s s 1/2))))
-      (multiple-value-bind (v mu root-mu root-mu1)
+  (cond
+    ((descending-base-p u m)
+     ;; A&S 16.13.1 to 16.13.3, dropping the linear (and higher) terms
+     ;; in m.
+     (values (sin u) (cos u) 1))
+    (t
+     ;; A&S 16.12.2:
+     ;;
+     ;;   sn(u|m) = (1+sqrt(mu))*sn(v|mu)/(1+sqrt(mu)*sn^2(v|mu)).
+     ;;
+     ;; A&S 16.12.3 for cn:
+     ;;
+     ;;   cn(u|m) = cn(v|mu)*dn(v|mu)/(1+sqrt(mu)*sn^2(v|mu))
+     ;;
+     ;; A&S 16.12.4 for dn:
+     ;;
+     ;;   dn(u|m) = (dn^2(v|mu) - (1-sqrt(mu)))/((1+sqrt(mu))-dn^2(v|mu))
+     (multiple-value-bind (v mu root-mu root-mu1)
 	  (descending-transform u m root-m1)
 	(multiple-value-bind (s c d)
 	    (elliptic-sn-cn-dn-descending v mu root-mu1)
@@ -136,14 +167,17 @@
 		 (den (1+ p)))
 	    (values (/ (* (1+ root-mu) s) den)
 		    (/ (* c d) den)
-		    ;; 1-p cancels when p is near 1.  Then use
-		    ;; 1-p = (1-sqrt(mu)) + sqrt(mu)*cn^2, where
-		    ;; 1-sqrt(mu) = 2*sqrt(m1)/(1+sqrt(m1)).
+		    ;; 1-p cancels when p is near 1, which happens
+		    ;; when m is near 1 and u is near an odd multiple
+		    ;; of K(m), where sn(v|mu) is near 1.  Then use
+		    ;; 1-p = 1-sqrt(mu)*sn^2 = 1-sqrt(mu)*(1-cn^2) =
+		    ;; (1-sqrt(mu)) + sqrt(mu)*cn^2, where 1-sqrt(mu)
+		    ;; = 2*sqrt(m1)/(1+sqrt(m1)).
 		    (/ (if (< (* 2 (abs p)) 1)
 			   (- 1 p)
 			   (+ (/ (* 2 root-m1) (+ 1 root-m1))
 			      (* root-mu c c)))
-		       den)))))))
+		       den))))))))
 
 (defun sn (u m)
   (cond ((zerop m)
@@ -187,31 +221,7 @@
 	 ;; cancellation in 1-p unless u and m are real and m < 8/9.
 	 (nth-value 2 (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))))
 	(t
-	 ;; Use the Gauss transformation from
-	 ;; http://functions.wolfram.com/09.29.16.0013.01:
-	 ;;
-	 ;;
-	 ;; dn((1+sqrt(m))*z, 4*sqrt(m)/(1+sqrt(m))^2)
-	 ;;   =  (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
-	 ;;
-	 ;; So
-	 ;;
-	 ;; dn(y, mu) = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
-	 ;;
-	 ;; where z = y/(1+sqrt(m)) and mu=4*sqrt(m)/(1+sqrt(m))^2.
-	 ;;
-	 ;; Solve for m, and we get
-	 ;;
-	 ;; sqrt(m) = -(mu+2*sqrt(1-mu)-2)/mu or (-mu+2*sqrt(1-mu)+2)/mu.
-	 ;;
-	 ;; I don't think it matters which sqrt we use, so I (rtoy)
-	 ;; arbitrarily choose the first one above.
-	 ;;
-	 ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as
-	 ;; -(mu+2*sqrt(1-mu)-2)/mu and as mu/(1+sqrt(1-mu))^2.  The
-	 ;; last one is the most accurate for small mu.  The parameter
-	 ;; of sn is root^2, whose complement has the square root
-	 ;; 2*sqrt(sqrt(1-mu))/(1+sqrt(1-mu)).
+         ;; 
 	 (multiple-value-bind (z root-sq root root-1-root-sq)
 	     (descending-transform u m (sqrt (- 1 m)))
 	   (let* ((s (elliptic-sn-descending z root-sq root-1-root-sq))

commit da8fa6804f1bcd82abad1e85c4c6612b72c28649
Author: Raymond Toy <toy...@gm...>
Date:   Thu Oct 15 18:39:03 2026 -0700

    Split elliptic-sn-descending into 2 routines
    
    elliptic-sn-descending now only handles sn.  elliptic-cn-dn-descending
    handles cn and dn because these transforms require computing sn, cn,
    and dn to get the final result.
    
    I think this is easier to read as well.  This also simplifies the
    calling convention too so that there no optional args anymore.
    
    This code was generated by Claude.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index a4d6ada78..920d9b33d 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -51,34 +51,7 @@
 ;;; functions and integrals.
 (in-package #:bigfloat)
 
-(declaim (inline descending-transform ascending-transform))
-
-(defun ascending-transform (u m)
-  ;; A&S 16.14.1
-  ;;
-  ;; Take care in computing this transform.  For the case where
-  ;; m is complex, we should compute sqrt(mu1) first as
-  ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
-  ;; If not, we may choose the wrong branch when computing
-  ;; sqrt(mu1).
-  (let* ((root-m (sqrt m))
-	 (mu (/ (* 4 root-m)
-		(expt (1+ root-m) 2)))
-	 (root-mu1 (/ (- 1 root-m) (+ 1 root-m)))
-	 (v (/ u (1+ root-mu1))))
-    (values v mu root-mu1)))
-
-#+nil
-(defun descending-transform (u m)
-  ;; Note: Don't calculate mu first, as given in 16.12.1.  We
-  ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
-  ;; then compute mu = sqrt(mu)^2.  If we calculate mu first,
-  ;; sqrt(mu) loses information when m or m1 is complex.
-  (let* ((root-m1 (sqrt (- 1 m)))
-	 (root-mu (/ (- 1 root-m1) (+ 1 root-m1)))
-	 (mu (* root-mu root-mu))
-	 (v (/ u (1+ root-mu))))
-    (values v mu root-mu)))
+(declaim (inline descending-transform))
 
 ;; Landen descending transform (Gauss' transform).
 ;;
@@ -104,216 +77,73 @@
 ;; However sqrt(1-mu) is computed as 2*sqrt(m1)/(1+sqrt(m1)) to get
 ;; rid of the round-off in 1-mu when m is near 1.
 (defun descending-transform (u m root-m1)
+  ;; A&S 16.12.1.  ROOT-M1 is sqrt(1-m).  Return v = u/(1+sqrt(mu)),
+  ;; mu, sqrt(mu), and sqrt(1-mu).
+  ;;
+  ;; Don't calculate mu first, as given in 16.12.1.  We should
+  ;; calculate sqrt(mu) = (1-sqrt(m1))/(1+sqrt(m1)), and then compute
+  ;; mu = sqrt(mu)^2.  If we calculate mu first, sqrt(mu) loses
+  ;; information when m or m1 is complex.
+  ;;
+  ;; When |sqrt(m1)| > 1/2, sqrt(mu) is computed as m/(1+sqrt(m1))^2
+  ;; instead, which avoids the cancellation in 1-sqrt(m1) for small m.
+  ;; Otherwise m is near 1 and has an absolute error that is large
+  ;; relative to m1, so sqrt(mu) is computed from ROOT-M1 alone.
+  ;; sqrt(1-mu) is computed as 2*sqrt(sqrt(m1))/(1+sqrt(m1)), without
+  ;; the cancellation in 1-mu when m is near 1.
   (let* ((den (+ 1 root-m1))
-	 (root-mu (/ m (* den den)))
+	 (root-mu (if (<= (* 2 (abs root-m1)) 1)
+		      (/ (- 1 root-m1) den)
+		      (/ m (* den den))))
 	 (mu (* root-mu root-mu))
 	 (v (/ u (1+ root-mu))))
     (values v mu root-mu (/ (* 2 (sqrt root-m1)) den))))
 
-;;
-;; This appears to work quite well for both real and complex values
-;; of u.
-#+nil
-(defun elliptic-sn-descending (u m)
-  (cond ((= m 1)
-	 ;; A&S 16.6.1
-	 (tanh u))
-	((< (abs m) (epsilon u))
-	 ;; A&S 16.6.1
-	 (sin u))
-	(t
-	 (multiple-value-bind (v mu root-mu)
-	     (descending-transform u m)
-	   (let* ((new-sn (elliptic-sn-descending v mu)))
-	     (/ (* (1+ root-mu) new-sn)
-		(1+ (* root-mu new-sn new-sn))))))))
-
-#+nil
-(defun elliptic-sn-descending (u m &optional cn-dn-p
-				      (root-m1 (sqrt (- 1 m)) root-m1-p))
-
-   (cond ((= m 1)
-	 ;; A&S 16.6.1 to 16.6.3
-	 (values (tanh u) (/ (cosh u)) (/ (cosh u))))
-	((and (not root-m1-p) (realp m) (> m 1))
-	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
-	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
-	 ;; descending transformation goes through complex values even
-	 ;; for real u, and loses digits.  This can only happen at the
-	 ;; top, where ROOT-M1 is not passed.
-	 (let ((root-m (sqrt m)))
-	   (multiple-value-bind (s c d)
-	       (elliptic-sn-descending (* root-m u) (/ m) cn-dn-p
-				       (sqrt (/ (- m 1) m)))
-	     (values (/ s root-m) d c))))
-	((and (< (abs m) (epsilon u))
-	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
-	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
-	      (or (realp u)
-		  (zerop m)
-		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
-	 ;; A&S 16.13.1 to 16.13.3
-	 (if cn-dn-p
-	     (let ((s (sin u)))
-	       (values s (cos u) (- 1 (* m s s 1/2))))
-	     (sin u)))
-	(t
-	 (multiple-value-bind (v mu root-mu root-mu1)
-	     (descending-transform u m root-m1)
-	   (multiple-value-bind (new-sn new-cn new-dn)
-	       (elliptic-sn-descending v mu cn-dn-p root-mu1)
-	     ;; A&S 16.12.2 to 16.12.4
-	     (let* ((p (* root-mu new-sn new-sn))
-		    (den (1+ p))
-		    (sn (/ (* (1+ root-mu) new-sn) den)))
-	       (if cn-dn-p
-		   (values sn
-			   (/ (* new-cn new-dn) den)
-			   ;; 1-p cancels when p is near 1.  Then use
-			   ;; 1-p = (1-root-mu) + root-mu*new-cn^2.
-			   (/ (if (< (abs p) 1/2)
-				  (- 1 p)
-				  (+ (/ (* 2 root-m1) (+ 1 root-m1))
-				     (* root-mu new-cn new-cn)))
-			      den))
-		   sn)))))))
-
-(defun impl-elliptic-sn-descending (u m root-m1)
-   (cond ((= m 1)
-	 ;; A&S 16.6.1 to 16.6.3
-	 (values (tanh u)))
-	((and (< (abs m) (epsilon u))
-	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
-	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
-	      (or (realp u)
-		  (zerop m)
-		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
-	 ;; A&S 16.13.1 to 16.13.3
-	 (sin u))
-	(t
-	 (multiple-value-bind (v mu root-mu root-mu1)
-	     (descending-transform u m root-m1)
-	   (multiple-value-bind (new-sn)
-	       (impl-elliptic-sn-descending v mu root-mu1)
-	     ;; A&S 16.12.2 to 16.12.4
-	     (let* ((p (* root-mu new-sn new-sn))
-		    (den (1+ p))
-		    (sn (/ (* (1+ root-mu) new-sn) den)))
-	       sn))))))
-
-(defun elliptic-sn-descending (u m)
-  (cond ((= m 1)
-         (tanh u))
-        ((and (realp m) (> m 1))
-	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
-	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
-	 ;; descending transformation goes through complex values even
-	 ;; for real u, and loses digits.  This can only happen at the
-	 ;; top, where ROOT-M1 is not passed.
-	 (let ((root-m (sqrt m)))
-	   (multiple-value-bind (s)
-	       (impl-elliptic-sn-descending (* root-m u) (/ m) (sqrt (/ (- m 1) m)))
-	     (/ s root-m))))
-	(t
-         (impl-elliptic-sn-descending u m (sqrt (- 1 m))))))
-        
-
-(defun elliptic-cn-descending (u m &optional cn-dn-p
-				      (root-m1 (sqrt (- 1 m)) root-m1-p))
-   (cond ((= m 1)
-	 ;; A&S 16.6.1 to 16.6.3
-	 (values (tanh u) (/ (cosh u)) (/ (cosh u))))
-	((and (not root-m1-p) (realp m) (> m 1))
-	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
-	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
-	 ;; descending transformation goes through complex values even
-	 ;; for real u, and loses digits.  This can only happen at the
-	 ;; top, where ROOT-M1 is not passed.
-	 (let ((root-m (sqrt m)))
-	   (multiple-value-bind (s c d)
-	       (elliptic-cn-descending (* root-m u) (/ m) cn-dn-p
-				       (sqrt (/ (- m 1) m)))
-	     (values (/ s root-m) d c))))
-	((and (< (abs m) (epsilon u))
-	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
-	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
-	      (or (realp u)
-		  (zerop m)
-		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
-	 ;; A&S 16.13.1 to 16.13.3
-	 (if cn-dn-p
-	     (let ((s (sin u)))
-	       (values s (cos u) (- 1 (* m s s 1/2))))
-	     (sin u)))
-	(t
-	 (multiple-value-bind (v mu root-mu root-mu1)
-	     (descending-transform u m root-m1)
-	   (multiple-value-bind (new-sn new-cn new-dn)
-	       (elliptic-cn-descending v mu cn-dn-p root-mu1)
-	     ;; A&S 16.12.2 to 16.12.4
-	     (let* ((p (* root-mu new-sn new-sn))
-		    (den (1+ p))
-		    (sn (/ (* (1+ root-mu) new-sn) den)))
-	       (if cn-dn-p
-		   (values sn
-			   (/ (* new-cn new-dn) den)
-			   ;; 1-p cancels when p is near 1.  Then use
-			   ;; 1-p = (1-root-mu) + root-mu*new-cn^2.
-			   (/ (if (< (abs p) 1/2)
-				  (- 1 p)
-				  (+ (/ (* 2 root-m1) (+ 1 root-m1))
-				     (* root-mu new-cn new-cn)))
-			      den))
-		   sn)))))))
-
-(defun elliptic-dn-descending (u m &optional cn-dn-p
-				      (root-m1 (sqrt (- 1 m)) root-m1-p))
-
-   (cond ((= m 1)
-	 ;; A&S 16.6.1 to 16.6.3
-	 (values (tanh u) (/ (cosh u)) (/ (cosh u))))
-	((and (not root-m1-p) (realp m) (> m 1))
-	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
-	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
-	 ;; descending transformation goes through complex values even
-	 ;; for real u, and loses digits.  This can only happen at the
-	 ;; top, where ROOT-M1 is not passed.
-	 (let ((root-m (sqrt m)))
-	   (multiple-value-bind (s c d)
-	       (elliptic-dn-descending (* root-m u) (/ m) cn-dn-p
-				       (sqrt (/ (- m 1) m)))
-	     (values (/ s root-m) d c))))
-	((and (< (abs m) (epsilon u))
-	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
-	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
-	      (or (realp u)
-		  (zerop m)
-		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
-	 ;; A&S 16.13.1 to 16.13.3
-	 (if cn-dn-p
-	     (let ((s (sin u)))
-	       (values s (cos u) (- 1 (* m s s 1/2))))
-	     (sin u)))
-	(t
-	 (multiple-value-bind (v mu root-mu root-mu1)
-	     (descending-transform u m root-m1)
-	   (multiple-value-bind (new-sn new-cn new-dn)
-	       (elliptic-dn-descending v mu cn-dn-p root-mu1)
-	     ;; A&S 16.12.2 to 16.12.4
-	     (let* ((p (* root-mu new-sn new-sn))
-		    (den (1+ p))
-		    (sn (/ (* (1+ root-mu) new-sn) den)))
-	       (if cn-dn-p
-		   (values sn
-			   (/ (* new-cn new-dn) den)
-			   ;; 1-p cancels when p is near 1.  Then use
-			   ;; 1-p = (1-root-mu) + root-mu*new-cn^2.
-			   (/ (if (< (abs p) 1/2)
-				  (- 1 p)
-				  (+ (/ (* 2 root-m1) (+ 1 root-m1))
-				     (* root-mu new-cn new-cn)))
-			      den))
-		   sn)))))))
+
+(defun descending-base-p (u m)
+  ;; True if m is small enough to stop the descending transformation.
+  ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
+  ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
+  (and (< (abs m) (epsilon u))
+       (or (realp u)
+	   (zerop m)
+	   (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u)))))))))
+
+;; sn(u|m) for m /= 1.  ROOT-M1 is sqrt(1-m).
+(defun elliptic-sn-descending (u m root-m1)
+  (if (descending-base-p u m)
+      ;; A&S 16.13.1
+      (sin u)
+      (multiple-value-bind (v mu root-mu root-mu1)
+	  (descending-transform u m root-m1)
+	(let ((s (elliptic-sn-descending v mu root-mu1)))
+	  ;; A&S 16.12.2
+	  (/ (* (1+ root-mu) s)
+	     (1+ (* root-mu s s)))))))
+
+;; sn(u|m), cn(u|m), and dn(u|m) for m /= 1.  ROOT-M1 is sqrt(1-m).
+(defun elliptic-sn-cn-dn-descending (u m root-m1)
+  (if (descending-base-p u m)
+      ;; A&S 16.13.1 to 16.13.3
+      (let ((s (sin u)))
+	(values s (cos u) (- 1 (* m s s 1/2))))
+      (multiple-value-bind (v mu root-mu root-mu1)
+	  (descending-transform u m root-m1)
+	(multiple-value-bind (s c d)
+	    (elliptic-sn-cn-dn-descending v mu root-mu1)
+	  ;; A&S 16.12.2 to 16.12.4
+	  (let* ((p (* root-mu s s))
+		 (den (1+ p)))
+	    (values (/ (* (1+ root-mu) s) den)
+		    (/ (* c d) den)
+		    ;; 1-p cancels when p is near 1.  Then use
+		    ;; 1-p = (1-sqrt(mu)) + sqrt(mu)*cn^2, where
+		    ;; 1-sqrt(mu) = 2*sqrt(m1)/(1+sqrt(m1)).
+		    (/ (if (< (* 2 (abs p)) 1)
+			   (- 1 p)
+			   (+ (/ (* 2 root-m1) (+ 1 root-m1))
+			      (* root-mu c c)))
+		       den)))))))
 
 (defun sn (u m)
   (cond ((zerop m)
@@ -328,12 +158,13 @@
 	 ;;
 	 ;; sn(u,m) = tanh(u) + 1/4*(1-m)*(sinh(u)*cosh(u)-u)*sech(u)^2
 	 (tanh u))
+	((and (realp m) (> m 1))
+	 ;; A&S 16.11.1: sn(u|m) = sn(sqrt(m)*u|1/m)/sqrt(m)
+	 (let ((root-m (sqrt m)))
+	   (/ (elliptic-sn-descending (* root-m u) (/ m) (sqrt (/ (- m 1) m)))
+	      root-m)))
 	(t
-	 ;; Use the ascending Landen transformation to compute sn.
-	 (let ((s (elliptic-sn-descending u m)))
-	   (if (and (realp u) (realp m))
-	       (realpart s)
-	       s)))))
+	 (elliptic-sn-descending u m (sqrt (- 1 m))))))
 
 (defun dn (u m)
   (cond ((zerop m)
@@ -347,12 +178,14 @@
 	 ;;
 	 ;; dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
 	 (/ (cosh u)))
+	((and (realp m) (> m 1))
+	 ;; A&S 16.11.3: dn(u|m) = cn(sqrt(m)*u|1/m)
+	 (nth-value 1 (elliptic-sn-cn-dn-descending (* (sqrt m) u) (/ m)
+						    (sqrt (/ (- m 1) m)))))
 	((not (and (realp u) (realp m) (< m 8/9)))
 	 ;; The Gauss transformation below can lose digits through the
-	 ;; cancellation in 1-p unless u and m are real and m < 8/9, and
-	 ;; m > 1 needs A&S 16.11.  ELLIPTIC-SN-DESCENDING handles that,
-	 ;; but has to compute cn as well.
-	 (nth-value 2 (elliptic-dn-descending u m t)))
+	 ;; cancellation in 1-p unless u and m are real and m < 8/9.
+	 (nth-value 2 (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))))
 	(t
 	 ;; Use the Gauss transformation from
 	 ;; http://functions.wolfram.com/09.29.16.0013.01:
@@ -375,19 +208,16 @@
 	 ;; arbitrarily choose the first one above.
 	 ;;
 	 ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as
-	 ;; -(mu+2*sqrt(1-mu)-2)/mu.  Also, the former is more
-	 ;; accurate for small mu.
-	 (let* (#+nil
-                (root (let ((root-1-m (sqrt (- 1 m))))
-			(/ (- 1 root-1-m)
-			   (+ 1 root-1-m))))
-                (root (let ((root-1-m (sqrt (- 1 m))))
-                        (/ m (expt (+ 1 root-1-m) 2))))
-		(z (/ u (+ 1 root)))
-		(s (elliptic-dn-descending z (* root root)))
-		(p (* root s s )))
-	   (/ (- 1 p)
-	      (+ 1 p))))))
+	 ;; -(mu+2*sqrt(1-mu)-2)/mu and as mu/(1+sqrt(1-mu))^2.  The
+	 ;; last one is the most accurate for small mu.  The parameter
+	 ;; of sn is root^2, whose complement has the square root
+	 ;; 2*sqrt(sqrt(1-mu))/(1+sqrt(1-mu)).
+	 (multiple-value-bind (z root-sq root root-1-root-sq)
+	     (descending-transform u m (sqrt (- 1 m)))
+	   (let* ((s (elliptic-sn-descending z root-sq root-1-root-sq))
+		  (p (* root s s)))
+	     (/ (- 1 p)
+		(+ 1 p)))))))
 
 (defun cn (u m)
   (cond ((zerop m)
@@ -402,17 +232,13 @@
 	 ;;
 	 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)
 	 (/ (cosh u)))
-	#+nil
-        (t
-	 ;; Use the ascending Landen transformation, A&S 16.14.3.
-	 (multiple-value-bind (v mu root-mu1)
-	     (ascending-transform u m)
-	   (let ((d (dn v mu)))
-	     (* (/ (+ 1 root-mu1) mu)
-		(/ (- (* d d) root-mu1)
-		   d)))))
-        (t
-	 (nth-value 1 (elliptic-cn-descending u m t)))))
+	((and (realp m) (> m 1))
+	 ;; A&S 16.11.2: cn(u|m) = dn(sqrt(m)*u|1/m)
+	 (nth-value 2 (elliptic-sn-cn-dn-descending (* (sqrt m) u) (/ m)
+						    (sqrt (/ (- m 1) m)))))
+	(t
+	 (nth-value 1 (elliptic-sn-cn-dn-descending u m (sqrt (- 1 m)))))))
+
 
 ;; Arithmetic-Geometric Mean algorithm for real or complex numbers.
 ;; See https://dlmf.nist.gov/22.20.ii.

commit 880b9be8889c507af2b5ae04effd57012acac976
Author: Raymond Toy <toy...@gm...>
Date:   Mon Oct 12 14:47:34 2026 -0700

    Refactor elliptic-sn-descending into multiple functions.
    
    First split elliptic-sn-descending into 3 functions, one for sn, cn,
    and dn.
    
    Second, split elliptic-sn-descending into a two functions: one that
    implements the recursion and the main entry that computes sn (named
    elliptic-sn-descending) that takes just two args.  The
    impl (impl-elliptic-sn-descending) takes 3 args and is used by
    elliptic-sn-descending.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index db0a556df..a4d6ada78 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -128,6 +128,7 @@
 	     (/ (* (1+ root-mu) new-sn)
 		(1+ (* root-mu new-sn new-sn))))))))
 
+#+nil
 (defun elliptic-sn-descending (u m &optional cn-dn-p
 				      (root-m1 (sqrt (- 1 m)) root-m1-p))
 
@@ -175,7 +176,144 @@
 				  (+ (/ (* 2 root-m1) (+ 1 root-m1))
 				     (* root-mu new-cn new-cn)))
 			      den))
-		   sn))))))) 
+		   sn)))))))
+
+(defun impl-elliptic-sn-descending (u m root-m1)
+   (cond ((= m 1)
+	 ;; A&S 16.6.1 to 16.6.3
+	 (values (tanh u)))
+	((and (< (abs m) (epsilon u))
+	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
+	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
+	      (or (realp u)
+		  (zerop m)
+		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
+	 ;; A&S 16.13.1 to 16.13.3
+	 (sin u))
+	(t
+	 (multiple-value-bind (v mu root-mu root-mu1)
+	     (descending-transform u m root-m1)
+	   (multiple-value-bind (new-sn)
+	       (impl-elliptic-sn-descending v mu root-mu1)
+	     ;; A&S 16.12.2 to 16.12.4
+	     (let* ((p (* root-mu new-sn new-sn))
+		    (den (1+ p))
+		    (sn (/ (* (1+ root-mu) new-sn) den)))
+	       sn))))))
+
+(defun elliptic-sn-descending (u m)
+  (cond ((= m 1)
+         (tanh u))
+        ((and (realp m) (> m 1))
+	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
+	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
+	 ;; descending transformation goes through complex values even
+	 ;; for real u, and loses digits.  This can only happen at the
+	 ;; top, where ROOT-M1 is not passed.
+	 (let ((root-m (sqrt m)))
+	   (multiple-value-bind (s)
+	       (impl-elliptic-sn-descending (* root-m u) (/ m) (sqrt (/ (- m 1) m)))
+	     (/ s root-m))))
+	(t
+         (impl-elliptic-sn-descending u m (sqrt (- 1 m))))))
+        
+
+(defun elliptic-cn-descending (u m &optional cn-dn-p
+				      (root-m1 (sqrt (- 1 m)) root-m1-p))
+   (cond ((= m 1)
+	 ;; A&S 16.6.1 to 16.6.3
+	 (values (tanh u) (/ (cosh u)) (/ (cosh u))))
+	((and (not root-m1-p) (realp m) (> m 1))
+	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
+	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
+	 ;; descending transformation goes through complex values even
+	 ;; for real u, and loses digits.  This can only happen at the
+	 ;; top, where ROOT-M1 is not passed.
+	 (let ((root-m (sqrt m)))
+	   (multiple-value-bind (s c d)
+	       (elliptic-cn-descending (* root-m u) (/ m) cn-dn-p
+				       (sqrt (/ (- m 1) m)))
+	     (values (/ s root-m) d c))))
+	((and (< (abs m) (epsilon u))
+	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
+	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
+	      (or (realp u)
+		  (zerop m)
+		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
+	 ;; A&S 16.13.1 to 16.13.3
+	 (if cn-dn-p
+	     (let ((s (sin u)))
+	       (values s (cos u) (- 1 (* m s s 1/2))))
+	     (sin u)))
+	(t
+	 (multiple-value-bind (v mu root-mu root-mu1)
+	     (descending-transform u m root-m1)
+	   (multiple-value-bind (new-sn new-cn new-dn)
+	       (elliptic-cn-descending v mu cn-dn-p root-mu1)
+	     ;; A&S 16.12.2 to 16.12.4
+	     (let* ((p (* root-mu new-sn new-sn))
+		    (den (1+ p))
+		    (sn (/ (* (1+ root-mu) new-sn) den)))
+	       (if cn-dn-p
+		   (values sn
+			   (/ (* new-cn new-dn) den)
+			   ;; 1-p cancels when p is near 1.  Then use
+			   ;; 1-p = (1-root-mu) + root-mu*new-cn^2.
+			   (/ (if (< (abs p) 1/2)
+				  (- 1 p)
+				  (+ (/ (* 2 root-m1) (+ 1 root-m1))
+				     (* root-mu new-cn new-cn)))
+			      den))
+		   sn)))))))
+
+(defun elliptic-dn-descending (u m &optional cn-dn-p
+				      (root-m1 (sqrt (- 1 m)) root-m1-p))
+
+   (cond ((= m 1)
+	 ;; A&S 16.6.1 to 16.6.3
+	 (values (tanh u) (/ (cosh u)) (/ (cosh u))))
+	((and (not root-m1-p) (realp m) (> m 1))
+	 ;; A&S 16.11: sn(u|m) = sn(v|1/m)/sqrt(m), cn(u|m) = dn(v|1/m)
+	 ;; and dn(u|m) = cn(v|1/m), where v = sqrt(m)*u.  For m > 1, the
+	 ;; descending transformation goes through complex values even
+	 ;; for real u, and loses digits.  This can only happen at the
+	 ;; top, where ROOT-M1 is not passed.
+	 (let ((root-m (sqrt m)))
+	   (multiple-value-bind (s c d)
+	       (elliptic-dn-descending (* root-m u) (/ m) cn-dn-p
+				       (sqrt (/ (- m 1) m)))
+	     (values (/ s root-m) d c))))
+	((and (< (abs m) (epsilon u))
+	      ;; sn(u|m) = sin(u) + O(m*sin(u)*cos(u)^2), A&S 16.13.1, and
+	      ;; cos(u)^2 grows like exp(2*abs(imagpart(u))).
+	      (or (realp u)
+		  (zerop m)
+		  (< (abs m) (* (epsilon u) (exp (* -2 (abs (imagpart u))))))))
+	 ;; A&S 16.13.1 to 16.13.3
+	 (if cn-dn-p
+	     (let ((s (sin u)))
+	       (values s (cos u) (- 1 (* m s s 1/2))))
+	     (sin u)))
+	(t
+	 (multiple-value-bind (v mu root-mu root-mu1)
+	     (descending-transform u m root-m1)
+	   (multiple-value-bind (new-sn new-cn new-dn)
+	       (elliptic-dn-descending v mu cn-dn-p root-mu1)
+	     ;; A&S 16.12.2 to 16.12.4
+	     (let* ((p (* root-mu new-sn new-sn))
+		    (den (1+ p))
+		    (sn (/ (* (1+ root-mu) new-sn) den)))
+	       (if cn-dn-p
+		   (values sn
+			   (/ (* new-cn new-dn) den)
+			   ;; 1-p cancels when p is near 1.  Then use
+			   ;; 1-p = (1-root-mu) + root-mu*new-cn^2.
+			   (/ (if (< (abs p) 1/2)
+				  (- 1 p)
+				  (+ (/ (* 2 root-m1) (+ 1 root-m1))
+				     (* root-mu new-cn new-cn)))
+			      den))
+		   sn)))))))
 
 (defun sn (u m)
   (cond ((zerop m)
@@ -214,7 +352,7 @@
 	 ;; cancellation in 1-p unless u and m are real and m < 8/9, and
 	 ;; m > 1 needs A&S 16.11.  ELLIPTIC-SN-DESCENDING handles that,
 	 ;; but has to compute cn as well.
-	 (nth-value 2 (elliptic-sn-descending u m t)))
+	 (nth-value 2 (elliptic-dn-descending u m t)))
 	(t
 	 ;; Use the Gauss transformation from
 	 ;; http://functions.wolfram.com/09.29.16.0013.01:
@@ -246,7 +384,7 @@
                 (root (let ((root-1-m (sqrt (- 1 m))))
                         (/ m (expt (+ 1 root-1-m) 2))))
 		(z (/ u (+ 1 root)))
-		(s (elliptic-sn-descending z (* root root)))
+		(s (elliptic-dn-descending z (* root root)))
 		(p (* root s s )))
 	   (/ (- 1 p)
 	      (+ 1 p))))))
@@ -274,7 +412,7 @@
 		(/ (- (* d d) root-mu1)
 		   d)))))
         (t
-	 (nth-value 1 (elliptic-sn-descending u m t)))))
+	 (nth-value 1 (elliptic-cn-descending u m t)))))
 
 ;; Arithmetic-Geometric Mean algorithm for real or complex numbers.
 ;; See https://dlmf.nist.gov/22.20.ii.

commit fd2dd75c88ca50e30d7b256e0cc78b5da7384207
Merge: 28e0b35b2 92979cf7b
Author: Raymond Toy <toy...@gm...>
Date:   Mon Oct 12 14:17:30 2026 -0700

    Merge branch 'master' into bug-5663-1-jacobi_sn-inaccurate


-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                   |   8 ++
 src/compar.lisp             | 100 +++++++++++-----
 src/defint.lisp             |  29 ++++-
 src/defmfun-check.lisp      |   8 +-
 src/ellipt.lisp             | 272 ++++++++++++++++++++------------------------
 src/float.lisp              |  39 +++++++
 src/mdot.lisp               |   2 +-
 src/mlisp.lisp              |  16 ++-
 src/numth.lisp              |   2 +-
 src/rat3c.lisp              |  22 +++-
 src/simp.lisp               |  31 ++---
 src/testsuite.lisp          |   6 +-
 src/trigi.lisp              |   4 +-
 tests/rtest16.mac           |   4 +-
 tests/rtest_elliptic.mac    |  53 ++++++---
 tests/rtest_gamma.mac       |  20 ++--
 tests/rtest_integrate.mac   |  17 +++
 tests/rtest_limit_extra.mac |   2 +-
 tests/rtest_rules.mac       |   9 ++
 tests/rtest_scalarp.mac     |  24 ++--
 tests/rtest_sign.mac        |  35 ++++++
 tests/rtest_trig.mac        |   9 ++
 22 files changed, 466 insertions(+), 246 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  a7f328122815ef2e700626619dcbe2a4016e1a51 (commit)
      from  739b9b02e0b2a13391cdeccd5b7b449c7ca31a60 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit a7f328122815ef2e700626619dcbe2a4016e1a51
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 13:52:48 2026 +0200

    Better fix for bug #5681
    
    Substitute the gensym back to '$ZEROA or '$ZEROB, otherwise a noun form
    would still contain the gensym and could reach the user.

diff --git a/src/compar.lisp b/src/compar.lisp
index a107464ad..30314cca9 100644
--- a/src/compar.lisp
+++ b/src/compar.lisp
@@ -912,16 +912,17 @@
 	(t (or (has-int-symbols (car e))
 	       (has-int-symbols (cdr e))))))
 
+;; The limit of A as EPS, which is '$ZEROA or '$ZEROB, goes to 0 from direction
+;; DIR. limit() rejects '$ZEROA and '$ZEROB as its variable, since they're
+;; constants, so use a fresh variable.
+(defun infinitesimal-limit (a eps dir)
+  (let ((v (gensym)))
+    (maxima-substitute eps v ($limit (maxima-substitute v eps a) v 0 dir))))
+
 ;;; Do substitutions for special symbols.
 (defun nmr (a)
-  ;; zeroa and zerob are constants, which limit() rejects as its variable,
-  ;; so take these limits in a fresh variable.
-  (unless (free a '$zeroa)
-    (let ((v (gensym)))
-      (setq a ($limit (maxima-substitute v '$zeroa a) v 0 '$plus))))
-  (unless (free a '$zerob)
-    (let ((v (gensym)))
-      (setq a ($limit (maxima-substitute v '$zerob a) v 0 '$minus))))
+  (unless (free a '$zeroa) (setq a (infinitesimal-limit a '$zeroa '$plus)))
+  (unless (free a '$zerob) (setq a (infinitesimal-limit a '$zerob '$minus)))
   (unless (free a 'z**) (setq a ($limit a 'z** 0 '$plus)))
   (unless (free a '*z*) (setq a ($limit a '*z* 0 '$plus)))
   (unless (free a 'epsilon) (setq a ($limit a 'epsilon 0 '$plus)))

-----------------------------------------------------------------------

Summary of changes:
 src/compar.lisp | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  739b9b02e0b2a13391cdeccd5b7b449c7ca31a60 (commit)
      from  5faac5ad1f98399a846127405930a1d4fa004ff5 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 739b9b02e0b2a13391cdeccd5b7b449c7ca31a60
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 10:39:22 2026 +0200

    Fix nested products from merging powers of %e
    
    %e^(log(x)/2+y)*%e^(log(x)/2) simplifies to x*%e^y. Inside a longer
    product, as in %e^(log(x)/2+y)*%e^(log(x)/2)*x, that result was kept as
    one factor, giving the nested product x*(%e^y*x), and
    is(%e^(log(x)/2+y)*%e^(log(x)/2)*x = x^2*%e^y) returned false. Now its
    factors are multiplied in one by one, as is already done for merged
    powers of abs().
    
    This fixes bug #5683.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 43753be26..e698b54c9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -227,6 +227,7 @@ Bug fixes for numbered bugs:
 * \#5679 Bigfloat atan(x) is wrong if x has a precision other than fpprec
 * \#5680 integrate(x*sec(x)^2,x,0,1) and some limits return results with gensyms
 * \#5681 integrate(tan(x),x,0,%pi/4) and limit(log(cos(x)),x,%pi/4) signal an error
+* \#5683 Merging powers of %e in a product can leave a nested product
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/simp.lisp b/src/simp.lisp
index 1fcd7a91c..d271e9409 100644
--- a/src/simp.lisp
+++ b/src/simp.lisp
@@ -2501,7 +2501,12 @@
                                           (list x w))))
             (if (and (not *timesinp*) (not (eq x '$%i)))
                 (let ((*timesinp* t))
-                  (setq temp (simplifya temp t))))))
+                  (setq temp (simplifya temp t))))
+            (when (mtimesp temp)
+              ;; x^w simplified to a product, as %e^(log(u)+v) = u*%e^v does.
+              ;; Multiply in its factors, not the product as one.
+              (setq x temp)
+              (go times))))
      (setq x (if (mexptp temp)
                  (cdr temp)
                  (list temp 1)))
diff --git a/tests/rtest16.mac b/tests/rtest16.mac
index 3754e4483..2922b3014 100644
--- a/tests/rtest16.mac
+++ b/tests/rtest16.mac
@@ -4481,6 +4481,14 @@ block([browser : browser, url_base : url_base, b : browser, u : url_base],
 [true, true]$
 
 
+/* Bug #5683: "Merging powers of %e in a product can leave a nested product" */
+
+[%e^(log(x)/2 + y)*%e^(log(x)/2)*x,
+ is(%e^(log(x)/2 + y)*%e^(log(x)/2)*x = x^2*%e^y),
+ %e^(%i*%pi/8 + y)*%e^(%i*%pi/8)*x];
+[%e^y*x^2, true, %e^y*(%i/sqrt(2) + 1/sqrt(2))*x]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog         | 1 +
 src/simp.lisp     | 7 ++++++-
 tests/rtest16.mac | 8 ++++++++
 3 files changed, 15 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  5faac5ad1f98399a846127405930a1d4fa004ff5 (commit)
      from  84d534aa83a6ed368535a9cdd85e241d535d1efe (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 5faac5ad1f98399a846127405930a1d4fa004ff5
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 10:30:49 2026 +0200

    Fix limit() and integrate() errors about zeroa
    
    limit(log(cos(x)),x,%pi/4,minus) and integrate(tan(x),x,0,%pi/4)
    failed with "limit: second argument must be a variable, not a
    constant; found: zeroa". Other limits of log() at points like %pi/3
    or %pi/4 failed the same way. Maxima 5.50.0 got them right.
    
    When the limit code asks for the sign of an expression that contains
    zeroa or zerob, sign() first takes a limit with zeroa or zerob as the
    variable. Since 5.50, zeroa and zerob are constants, which limit()
    does not accept as a variable. This went unnoticed until the limit of
    log() started to ask sign() this way. Now that inner limit uses an
    ordinary variable.
    
    This fixes bug #5681.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 0f2343383..43753be26 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -226,6 +226,7 @@ Bug fixes for numbered bugs:
 * \#5678 ratsimp() overflows the stack under some equal() facts
 * \#5679 Bigfloat atan(x) is wrong if x has a precision other than fpprec
 * \#5680 integrate(x*sec(x)^2,x,0,1) and some limits return results with gensyms
+* \#5681 integrate(tan(x),x,0,%pi/4) and limit(log(cos(x)),x,%pi/4) signal an error
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/compar.lisp b/src/compar.lisp
index 89035cba7..a107464ad 100644
--- a/src/compar.lisp
+++ b/src/compar.lisp
@@ -914,8 +914,14 @@
 
 ;;; Do substitutions for special symbols.
 (defun nmr (a)
-  (unless (free a '$zeroa) (setq a ($limit a '$zeroa 0 '$plus)))
-  (unless (free a '$zerob) (setq a ($limit a '$zerob 0 '$minus)))
+  ;; zeroa and zerob are constants, which limit() rejects as its variable,
+  ;; so take these limits in a fresh variable.
+  (unless (free a '$zeroa)
+    (let ((v (gensym)))
+      (setq a ($limit (maxima-substitute v '$zeroa a) v 0 '$plus))))
+  (unless (free a '$zerob)
+    (let ((v (gensym)))
+      (setq a ($limit (maxima-substitute v '$zerob a) v 0 '$minus))))
   (unless (free a 'z**) (setq a ($limit a 'z** 0 '$plus)))
   (unless (free a '*z*) (setq a ($limit a '*z* 0 '$plus)))
   (unless (free a 'epsilon) (setq a ($limit a 'epsilon 0 '$plus)))
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index c2572ea3d..deb6e7a36 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7259,6 +7259,15 @@ block([r : integrate(x*sec(x)^2, x, 0, %pi/4)],
 [[], 0]$
 
 
+/* Bug #5681: "integrate(tan(x),x,0,%pi/4) and limit(log(cos(x)),x,%pi/4) signal an error" */
+
+[integrate(tan(x), x, 0, %pi/4),
+ integrate(tan(x), x, 0, %pi/3),
+ integrate(sin(x)/cos(x), x, 0, %pi/4),
+ integrate(tan(x), x, -%pi/4, %pi/3)];
+[log(2)/2, log(2), log(2)/2, log(2)/2]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/
diff --git a/tests/rtest_limit.mac b/tests/rtest_limit.mac
index 60085de8a..18b82332c 100644
--- a/tests/rtest_limit.mac
+++ b/tests/rtest_limit.mac
@@ -1737,6 +1737,15 @@ inf;
 [y, y, y, y]$
 
 
+/* Bug #5681: "integrate(tan(x),x,0,%pi/4) and limit(log(cos(x)),x,%pi/4) signal an error" */
+
+[limit(log(cos(x)), x, %pi/4, minus),
+ limit(log(cos(x)), x, %pi/4),
+ limit(log(cos(x)), x, %pi/3, plus),
+ limit(log(sin(x)), x, %pi/4, plus)];
+[-log(2)/2, -log(2)/2, -log(2), -log(2)/2]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/compar.lisp           | 10 ++++++++--
 tests/rtest_integrate.mac |  9 +++++++++
 tests/rtest_limit.mac     |  9 +++++++++
 4 files changed, 27 insertions(+), 2 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  84d534aa83a6ed368535a9cdd85e241d535d1efe (commit)
      from  fa96ec969f968b18f298bc0c4c94d8c484658e5b (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 84d534aa83a6ed368535a9cdd85e241d535d1efe
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 10:25:20 2026 +0200

    Fix gensyms in results of limit() and integrate()
    
    limit() could return results containing internal symbols, such as
    y*g1587+y*(1-g1587) for limit((sin(x)^2+cos(x)^2)*y,x,0,plus). The
    symbols cancel, so the value was right, but not even float() could
    get rid of them. Many definite integrals are evaluated with limit(),
    so integrate() was affected too, e.g. integrate(x*sec(x)^2,x,0,1).
    
    limit() replaces sin(x)^2+cos(x)^2 by 1 with the help of an internal
    symbol. In an expanded sum such as y*sin(x)^2+y*cos(x)^2, the symbol
    cancels only after rational simplification. limit() checked that it
    cancels, but then kept the expression from before that simplification.
    Now it keeps the simplified one.
    
    This fixes bug #5680.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 17d248928..0f2343383 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -225,6 +225,7 @@ Bug fixes for numbered bugs:
 * \#5677 risch() returns wrong polylogarithms for some log integrands
 * \#5678 ratsimp() overflows the stack under some equal() facts
 * \#5679 Bigfloat atan(x) is wrong if x has a precision other than fpprec
+* \#5680 integrate(x*sec(x)^2,x,0,1) and some limits return results with gensyms
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/limit.lisp b/src/limit.lisp
index 47b5a76a8..41038fbd5 100644
--- a/src/limit.lisp
+++ b/src/limit.lisp
@@ -805,7 +805,9 @@ ignoring dummy variables and array indices."
                     (setq z (gensym))
                     (setq ee (maxima-substitute z (power (ftake '%cos g) 2) e))
                     (setq ee (maxima-substitute (sub 1 z) (power (ftake '%sin g) 2) ee))
-                    (when (freeof z (sratsimp ee))
+                    (unless (freeof z ee)
+                      (setq ee (sratsimp ee)))
+                    (when (freeof z ee)
                         (setq e ee)))
                 e)
             ;; maybe this isn't needed, but it's not wrong.
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 0229402ba..c2572ea3d 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7247,6 +7247,18 @@ block([g, G, a, e, q],
 [0, 0]$
 
 
+/* Bug #5680: "integrate(x*sec(x)^2,x,0,1) and some limits return results with gensyms" */
+
+block([r, v : tan(1) + log(cos(1))],
+  r : [integrate(x*sec(x)^2, x, 0, 1), integrate(x*tan(x)^2, x, 0, 1)],
+  [listofvars(r), is(lmax(abs(float(r - [v, v - 1/2]))) < 1e-12)]);
+[[], true]$
+
+block([r : integrate(x*sec(x)^2, x, 0, %pi/4)],
+  [listofvars(r), radcan(r - (%pi/4 - log(2)/2))]);
+[[], 0]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/
diff --git a/tests/rtest_limit.mac b/tests/rtest_limit.mac
index bedd0a2b4..60085de8a 100644
--- a/tests/rtest_limit.mac
+++ b/tests/rtest_limit.mac
@@ -1728,6 +1728,15 @@ limit(1/x + a/2^x, x, 0, plus);
 inf;
 
 
+/* Bug #5680: "integrate(x*sec(x)^2,x,0,1) and some limits return results with gensyms" */
+
+[limit((sin(x)^2 + cos(x)^2)*y, x, 0),
+ limit((sin(x)^2 + cos(x)^2)*y, x, 0, plus),
+ limit((sin(x)^2 + cos(x)^2)*y, x, inf),
+ limit((sin(2*x)^2 + cos(2*x)^2)*y + x, x, 0, plus)];
+[y, y, y, y]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/limit.lisp            |  4 +++-
 tests/rtest_integrate.mac | 12 ++++++++++++
 tests/rtest_limit.mac     |  9 +++++++++
 4 files changed, 25 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  fa96ec969f968b18f298bc0c4c94d8c484658e5b (commit)
      from  bcc4ffc889d6057e6191ef01cbf290d8bd127936 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit fa96ec969f968b18f298bc0c4c94d8c484658e5b
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 10:20:58 2026 +0200

    Fix atan() of a bigfloat made at another fpprec
    
    atan() of a real bigfloat returned a wrong value if the bigfloat had
    been made at another fpprec than the current one: atan() of 0.5b0
    made at the default fpprec returned 5.55b-17 at fpprec : 32 instead
    of 4.64b-1, and atan() of 2.0b0 made at fpprec : 32 returned nearly
    %pi/2 at fpprec : 16. Now atan() rounds the argument to the current
    precision first, as the other bigfloat functions already do.
    
    This fixes bug #5679.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 800ef24a7..17d248928 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -224,6 +224,7 @@ Bug fixes for numbered bugs:
 * \#5676 factor() returns a multiple of some polynomials with %i
 * \#5677 risch() returns wrong polylogarithms for some log integrands
 * \#5678 ratsimp() overflows the stack under some equal() facts
+* \#5679 Bigfloat atan(x) is wrong if x has a precision other than fpprec
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/float.lisp b/src/float.lisp
index 9e4694e88..afa20dc23 100644
--- a/src/float.lisp
+++ b/src/float.lisp
@@ -983,7 +983,7 @@
            (sub v
                 (mul '$%i u))))
          (t
-         (bcons (fpatan (cdr x))))))
+         (bcons (fpatan (cdr (bigfloatp x)))))))
       
 ;; atan(y/x) taking into account the quadrant.  (Also equal to
 ;; arg(x+%i*y).)
diff --git a/tests/rtest_trig.mac b/tests/rtest_trig.mac
index 3bc3a92fb..9d19d4f99 100644
--- a/tests/rtest_trig.mac
+++ b/tests/rtest_trig.mac
@@ -950,3 +950,13 @@ block([m],
   declare(m, integer),
   [cos(%pi*(4*m+3)^3), cos(%pi*(6*m+5)^3), sin(%pi*(4*m+3)^3+%pi/2)]);
 [-1, -1, -1]$
+
+
+/* Bug #5679: "Bigfloat atan(x) is wrong if x has a precision other than fpprec" */
+
+block([fpprec : 32, x : [1/2, 3/4, 2, -2], lo, hi],
+  lo : block([fpprec : 16], bfloat(x)),
+  hi : block([fpprec : 64], bfloat(x)),
+  x : bfloat(x),
+  [is(map(atan, lo) = map(atan, x)), is(map(atan, hi) = map(atan, x))]);
+[true, true]$

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog            |  1 +
 src/float.lisp       |  2 +-
 tests/rtest_trig.mac | 10 ++++++++++
 3 files changed, 12 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  bcc4ffc889d6057e6191ef01cbf290d8bd127936 (commit)
      from  8ebdaff39a963c38d4973b655c9ae40eb6eef5ce (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit bcc4ffc889d6057e6191ef01cbf290d8bd127936
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 00:18:15 2026 +0200

    sign(): fix stack overflow with equal() facts
    
    Under an equal() fact about a square root, such as
    equal(y,sqrt(e*x-(b*e)/(2*a)+(a*e)/(2*b))), ratsimp() of that root with
    algebraic : true recursed until the stack overflowed.
    
    To find abs(b) positive, sign() checks whether b is nonzero. That check
    applies the equal() facts, and simplifying the result can need the sign
    of abs(b) again. Now sign() does not start the check again while it is
    running.
    
    This fixes bug #5678.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 193d34a3b..800ef24a7 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -223,6 +223,7 @@ Bug fixes for numbered bugs:
 * \#5675 factor() and integrate() hang when %i is a common factor
 * \#5676 factor() returns a multiple of some polynomials with %i
 * \#5677 risch() returns wrong polylogarithms for some log integrands
+* \#5678 ratsimp() overflows the stack under some equal() facts
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/compar.lisp b/src/compar.lisp
index 296e8d7c8..89035cba7 100644
--- a/src/compar.lisp
+++ b/src/compar.lisp
@@ -2234,6 +2234,9 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
         odds (unless (member sign '($pos $neg $zero))
                (ncons (if (member arg-sign '($pos $pz)) (sub arg 1) x))))))
 
+(defvar *sign-mabs-args* nil
+  "The arguments of abs() that SIGN-MABS is testing with MNQP.")
+
 (defun sign-mabs (x)
   (let ((*complexsign* t))
     (sign (cadr x))
@@ -2241,7 +2244,12 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
 	  ((member sign '($neg $pn) :test #'eq) (setq sign '$pos))
 	  ;; abs(nonzero) > 0. If SIGN found nothing ('$PNZ), pass NIL as the third
 	  ;; argument to MNQP, allowing it to skip a redundant second $CSIGN call.
-	  ((eq t (mnqp 0 (cadr x) (not (eq sign '$pnz)))) (setq sign '$pos))
+	  ;; MNQP can need the sign of abs() of the same argument again, as when
+	  ;; it applies an equal() fact, so do not test it again then.
+	  ((and (not (member (cadr x) *sign-mabs-args* :test #'alike1))
+		(eq t (let ((*sign-mabs-args* (cons (cadr x) *sign-mabs-args*)))
+			(mnqp 0 (cadr x) (not (eq sign '$pnz))))))
+	   (setq sign '$pos))
 	  (t (setq sign '$pz minus nil evens (nconc odds evens) odds nil)))))
 
 (defun sign-asin/acos/atanh (x)
diff --git a/tests/rtest_sign.mac b/tests/rtest_sign.mac
index 7f1595b23..2fc944bb4 100644
--- a/tests/rtest_sign.mac
+++ b/tests/rtest_sign.mac
@@ -3126,6 +3126,16 @@ sqrt(1-2*x^2)/x$
 [pnz, sqrt(%e^2-1)]$
 
 
+/* Bug #5678: "ratsimp() overflows the stack under some equal() facts" */
+
+block([r],
+  assume(equal(y, sqrt(e*x-(b*e)/(2*a)+(a*e)/(2*b)))),
+  r : block([algebraic : true], ratsimp(sqrt(e*x-(b*e)/(2*a)+(a*e)/(2*b)))),
+  forget(equal(y, sqrt(e*x-(b*e)/(2*a)+(a*e)/(2*b)))),
+  r);
+sqrt((2*a*b*e*x+(a^2-b^2)*e)/(a*b))/sqrt(2)$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog            |  1 +
 src/compar.lisp      | 10 +++++++++-
 tests/rtest_sign.mac | 10 ++++++++++
 3 files changed, 20 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  8ebdaff39a963c38d4973b655c9ae40eb6eef5ce (commit)
      from  d0ec62c8077363cca291b57878094393ffd1eec0 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 8ebdaff39a963c38d4973b655c9ae40eb6eef5ce
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 00:08:17 2026 +0200

    risch(): fix polylogs with mixed-up variables
    
    risch((log(1-%i*x/sqrt(q-e^2))-log(1+%i*x/sqrt(q-e^2)))/(x+a),x)
    returns a wrong antiderivative: one of its polylogarithms is
    li[2](-((x-e)/(e+a))), where e stands in for %i*sqrt(q-e^2).
    DILOG writes the arguments of two logarithms as polynomials and divides
    one by the other, but the two polynomials can have different main
    variables, and then the quotient is wrong. Now DILOG gives up in that
    case, as it does when the quotient is not constant, and that part of the
    integral stays unevaluated.
    
    This fixes bug #5677.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 6e6f0c69b..193d34a3b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -222,6 +222,7 @@ Bug fixes for numbered bugs:
 * \#5674 factor() returns a wrong sign or hangs for some quotients with %i and square roots
 * \#5675 factor() and integrate() hang when %i is a common factor
 * \#5676 factor() returns a multiple of some polynomials with %i
+* \#5677 risch() returns wrong polylogarithms for some log integrands
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/risch.lisp b/src/risch.lisp
index 0a88b0dfe..3d0d70a0a 100644
--- a/src/risch.lisp
+++ b/src/risch.lisp
@@ -637,7 +637,11 @@
     (desetq (rc . rd) (ratdivide rn ro))
     (let
 	((result
-	  (cond ((and (freeof risch-intvar (rdis rc)) ;; can't use risch-constp because varlist
+	  ;; RN must be a polynomial in RISCH-VAR as RO is, as RATDIVIDE divides by
+	  ;; their main variables.
+	  (cond ((and (not (pcoefp (car rn)))
+		      (eq (p-var (car rn)) risch-var)
+		      (freeof risch-intvar (rdis rc)) ;; can't use risch-constp because varlist
 		      (freeof risch-intvar (rdis rd))) ;; is not set up with vars in correct order.
 		 (setq narg ($ratsimp (sub 1 (div narg (rdis rd)))))
 		 (mul* coef (power -1 (1+ risch-degree))
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 7c8c45b69..0229402ba 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7237,6 +7237,16 @@ block([x],
 (2*%e^x^(3/2)*(x^(3/2)-1))/3$
 
 
+/* Bug #5677: "risch() returns wrong polylogarithms for some log integrands" */
+
+block([g, G, a, e, q],
+  g : (log(1-%i*x/sqrt(q-e^2))-log(1+%i*x/sqrt(q-e^2)))/(x+a),
+  G : risch(g, x),
+  map(lambda([p], radcan(rectform(subst(p, diff(G, x) - g)))),
+      [[a = 2, e = 1, q = 5, x = 1/3], [a = -1, e = 2, q = 13, x = 3]]));
+[0, 0]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/risch.lisp            |  6 +++++-
 tests/rtest_integrate.mac | 10 ++++++++++
 3 files changed, 16 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  d0ec62c8077363cca291b57878094393ffd1eec0 (commit)
      from  f1d122e8f1e9a9b8db8a41b29799b01b0580e364 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit d0ec62c8077363cca291b57878094393ffd1eec0
Author: David Scherfgen <d.s...@go...>
Date:   Thu Oct 15 00:01:28 2026 +0200

    Fix wrong factor() results with %i
    
    factor() returned a constant multiple of many polynomials with %i in
    their coefficients, such as (4*%i+3)*(25*x-4*%i+3)*(2*y+x) for the
    expanded form of ((3+4*%i)*x+1)*(2*y+x). Of 400 random products of such
    polynomials, 221 came back wrong. ratsimp() could also fail with
    "CQUOTIENT: quotient is not exact" for some quotients with %i, and
    integrate() failed the same way for some integrands with %i.
    
    With %i or another algebraic number in the coefficients, the gcd of
    the coefficients of a polynomial can be a multiple of its content.
    OLDCONTENT then divided the polynomial by it and dropped the
    denominator of the quotient. Now it divides the gcd by that
    denominator if that is exact, and otherwise leaves the polynomial
    whole.
    
    This fixes bug #5676.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 7cfe1b6df..6e6f0c69b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -221,6 +221,7 @@ Bug fixes for numbered bugs:
 * \#5669 bfloat(%e^-1)^0.5b0 hangs for some values of fpprec
 * \#5674 factor() returns a wrong sign or hangs for some quotients with %i and square roots
 * \#5675 factor() and integrate() hang when %i is a common factor
+* \#5676 factor() returns a multiple of some polynomials with %i
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/rat3c.lisp b/src/rat3c.lisp
index 96a954450..6cc9d8426 100644
--- a/src/rat3c.lisp
+++ b/src/rat3c.lisp
@@ -327,10 +327,17 @@
   (cond ((pcoefp x) (list x 1))
 	((null (p-red x))
 	 (list (p-lc x) (make-poly (p-var x) (p-le x) 1)))
-	(t (let ((u (contsort (cdr x))) v)
+	(t (let ((u (contsort (cdr x))) v w)
 	     (setq u (oldcontent2 (cdr u) (car u))
-		   v (cond ($algebraic (car (rquotient x u)))
-			   (t (pcquotient x u))))
+		   v (cond ($algebraic (rquotient x u))
+			   (t (cons (pcquotient x u) 1))))
+	     ;; With algebraic numbers, U can be a multiple of the content
+	     ;; and X/U a fraction. Then U is divided by the denominator if
+	     ;; that is exact, else X is kept whole.
+	     (cond ((equal (cdr v) 1) (setq v (car v)))
+		   ((setq w (pquotientchk-safe u (cdr v)))
+		    (setq u w v (car v)))
+		   (t (setq u 1 v x)))
 	     (cond ((pminusp v) (list (pminus u) (pminus v)))
 		   (t (list u v)))))))
 
diff --git a/tests/rtest_gcd.mac b/tests/rtest_gcd.mac
index 3cc6b4a25..d779ef59a 100644
--- a/tests/rtest_gcd.mac
+++ b/tests/rtest_gcd.mac
@@ -203,3 +203,13 @@ block([e : (%i*a*x+1)^(3/2)/(a^2*c^(3/2)*x^2*(1-%i*a*x)^(3/2)*sqrt(a^2*x^2+1)
            +c^(3/2)*(1-%i*a*x)^(3/2)*sqrt(a^2*x^2+1))],
   ratsimp(block([algebraic : true, gcd : subres], factor(e)) - e));
 0$
+
+
+/* Bug #5676: "factor() returns a multiple of some polynomials with %i" */
+
+[factor(expand(((3+4*%i)*x+1)*(2*y+x))),
+ factor(expand((x-%i)*((2*x+%i)*y+2*x+a))),
+ block([algebraic : true, gcd : subres],
+   ratsimp((x-%i)*((2*x+%i)*y+2*x+a)/(2*y+2*%i)))];
+[(4*%i*x+3*x+1)*(2*y+x), (x-%i)*(2*x*y+%i*y+2*x+a),
+ ((2*x^2-%i*x+1)*y+2*x^2+(a-2*%i)*x-%i*a)/(2*y+2*%i)]$

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog           |  1 +
 src/rat3c.lisp      | 13 ++++++++++---
 tests/rtest_gcd.mac | 10 ++++++++++
 3 files changed, 21 insertions(+), 3 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  f1d122e8f1e9a9b8db8a41b29799b01b0580e364 (commit)
      from  f072b36b97de2983a3e24a11aa6b038910289ef5 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit f1d122e8f1e9a9b8db8a41b29799b01b0580e364
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 23:30:05 2026 +0200

    sign(): fix stack overflow when comparing powers
    
    sign() could recurse until the stack overflowed, which often crashed
    SBCL, as for sign(%e^x-%e^2) under assume(x>0), for
    limit(sqrt(%e^x-1),x,2), and for integrate(%e^(x^(3/2))*x^2,x) under
    assume(x>1).
    
    To compare two powers, sign() checks whether their bases are equal.
    That check can need the sign of a difference of powers again, and then
    it starts over. Now the rules for such differences are skipped while
    the check is running.
    
    This fixes bug #5668.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 59868d5ae..7cfe1b6df 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -217,6 +217,7 @@ Bug fixes for numbered bugs:
 * \#5665 integrate(exp(p(x)),x,a,b) with a cubic or quartic p takes minutes or exhausts the heap
 * \#5666 is(abs(z) >= 0) returns false when z contains %i
 * \#5667 constantp(log(2)) returns false once log has tellsimp rules
+* \#5668 limit(sqrt(%e^x-1), x, 2) causes stack overflow
 * \#5669 bfloat(%e^-1)^0.5b0 hangs for some values of fpprec
 * \#5674 factor() returns a wrong sign or hangs for some quotients with %i and square roots
 * \#5675 factor() and integrate() hang when %i is a common factor
diff --git a/src/compar.lisp b/src/compar.lisp
index 69f76ee6a..296e8d7c8 100644
--- a/src/compar.lisp
+++ b/src/compar.lisp
@@ -1712,11 +1712,14 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
       (setq odds (ncons x) evens nil minus nil))
     retval))
 
+(defvar *in-expt-of-base* nil
+  "T inside EXPT-OF-BASE, whose MEQP can recurse into SIGNDIFF-SPECIAL.")
+
 (defun expt-of-base (expr q)
   "If EXPR is a power of base Q (or EXPR = Q, then the exponent is implicitly 1),
   returns the exponent (1 or k). Otherwise returns NIL.
   Equality of EXPR's base and Q or EXPR and Q is checked via MEQP."
-  (let ($ratprint)
+  (let ($ratprint (*in-expt-of-base* t))
     (if (mexptp expr)
       (if (eq t (meqp (cadr expr) q))
         (caddr expr))
@@ -1743,6 +1746,7 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
 			;; Do NOT apply when Spos is itself a power of Qpos (e.g. S = Q):
 			;; The reduction would just toggle the exponent R <-> 1/R and recurse forever.
 			;; That same-base case is handled by the exponent-comparison rule further below.
+			(not *in-expt-of-base*) ; prevent endless recursion
 			(not (expt-of-base xrhs (cadr xlhs)))
 			(eq (sign* (sub (cadr xlhs)
 					(power xrhs (div 1 (caddr xlhs)))))
@@ -1825,7 +1829,8 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
     ;;     Q > 1 => sign(m - n),
     ;; 0 < Q < 1 => sign(n - m).
     ;; Resolves e.g. sign(x^a - x) = pos for x > 1, a > 1, and always terminates.
-    (when (null sgn)
+    (when (and (null sgn)
+               (not *in-expt-of-base*)) ; prevent endless recursion
       (let ((q (cond ((mexptp xlhs) (cadr xlhs))
                      ((mexptp xrhs) (cadr xrhs)))))
         (when q
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index df10c7ba4..7c8c45b69 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7221,6 +7221,22 @@ block([f, F],
 0$
 
 
+/* Bug #5668: "limit(sqrt(%e^x-1), x, 2) causes stack overflow"
+   (not actually a limit() bug) */
+
+block([x],
+  local(x),
+  assume(x > 1),
+  [integrate(%e^(x^(3/2))*x^2, x), integrate(%e^(x^(5/2))*x^4, x)]);
+[(2*%e^x^(3/2)*(x^(3/2)-1))/3, (2*%e^x^(5/2)*(x^(5/2)-1))/5]$
+
+block([x],
+  local(x),
+  assume(x > 0, x < 1),
+  integrate(%e^(x^(3/2))*x^2, x));
+(2*%e^x^(3/2)*(x^(3/2)-1))/3$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/
diff --git a/tests/rtest_sign.mac b/tests/rtest_sign.mac
index f21951a4d..7f1595b23 100644
--- a/tests/rtest_sign.mac
+++ b/tests/rtest_sign.mac
@@ -3107,6 +3107,25 @@ block([a, b, c, d, e, f, g],
 [[], [], []]$
 
 
+/* Bug #5668: "limit(sqrt(%e^x-1), x, 2) causes stack overflow"
+   (not actually a limit() bug) */
+
+block([r, x, y],
+  local(x, y),
+  assume(x > 0, equal(y, sqrt(1/x^2-2))),
+  [radcan(sqrt(1/x^2-2)), sign(x^2-1/2)]);
+[sqrt(1-2*x^2)/x, pnz]$
+
+block([r, x, y],
+  local(x, y),
+  assume(x > 0, x < 1, equal(y, sqrt(1/x^2-2))),
+  radcan(sqrt(1/x^2-2)));
+sqrt(1-2*x^2)/x$
+
+[block(local(x), assume(x > 0), sign(%e^x-%e^2)), limit(sqrt(%e^x-1), x, 2)];
+[pnz, sqrt(%e^2-1)]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/compar.lisp           |  9 +++++++--
 tests/rtest_integrate.mac | 16 ++++++++++++++++
 tests/rtest_sign.mac      | 19 +++++++++++++++++++
 4 files changed, 43 insertions(+), 2 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  f072b36b97de2983a3e24a11aa6b038910289ef5 (commit)
      from  6c20745f843a721413a46a12c46641d736875c15 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit f072b36b97de2983a3e24a11aa6b038910289ef5
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 23:14:52 2026 +0200

    factor(): fix endless loop on a unit such as %i
    
    To find how often a common factor divides the numerator and the
    denominator, factor() divided by it as long as the division was exact.
    With algebraic:true, a unit such as %i divides everything, so this
    never ended. integrate(%e^(3*%i*atan(a*x))/(a^2*c*x^2+c)^(3/2),x) hung
    this way.
    
    Now factor() does not count the powers of such algebraic factors.
    Further powers of one that is not a unit, such as 1+2^(1/3), then stay
    in the other factors, so the result is still correct.
    
    This needs the previous commit. Without it, some quotients for which
    factor() looped give a wrong sign instead.
    
    This fixes bug #5675.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index c521e6a30..59868d5ae 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -219,6 +219,7 @@ Bug fixes for numbered bugs:
 * \#5667 constantp(log(2)) returns false once log has tellsimp rules
 * \#5669 bfloat(%e^-1)^0.5b0 hangs for some values of fpprec
 * \#5674 factor() returns a wrong sign or hangs for some quotients with %i and square roots
+* \#5675 factor() and integrate() hang when %i is a common factor
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/lesfac.lisp b/src/lesfac.lisp
index 7b89fc67d..e5aaaf355 100644
--- a/src/lesfac.lisp
+++ b/src/lesfac.lisp
@@ -139,13 +139,16 @@
 	 (cond ((equal c 1) (go nextb))
 	       ((equal ai 1) (go bloop)))
 	aloop
-	 (when (setq t1 (testdivide ai c))
+	 ;; Do not count how often an algebraic C divides AI and BJ: a unit
+	 ;; such as %i divides them any number of times.
+	 (when (and (not (alg c)) (setq t1 (testdivide ai c)))
            (setq ai t1)
            (incf d1)
            (go aloop))
 	bloop
 	 (and (= d1 1)
 	      (not (equal bj 1))
+	      (not (alg c))
 	      (do ((t1
 		    (testdivide bj c)
 		    (testdivide bj c)))
diff --git a/tests/rtest_gcd.mac b/tests/rtest_gcd.mac
index df640b220..3cc6b4a25 100644
--- a/tests/rtest_gcd.mac
+++ b/tests/rtest_gcd.mac
@@ -195,3 +195,11 @@ block([N : x^2*sqrt(1-x)*sqrt(x+1)+%i*x^3-%i*x,
        D : %i*x*(x^2-1)*sqrt(1-x)*sqrt(x+1)-(x^2-1)^2],
   ratsimp(block([gcd : subres], factor(1/N+1/D))*N*D/(N+D)));
 1$
+
+
+/* Bug #5675: "factor() and integrate() hang when %i is a common factor" */
+
+block([e : (%i*a*x+1)^(3/2)/(a^2*c^(3/2)*x^2*(1-%i*a*x)^(3/2)*sqrt(a^2*x^2+1)
+           +c^(3/2)*(1-%i*a*x)^(3/2)*sqrt(a^2*x^2+1))],
+  ratsimp(block([algebraic : true, gcd : subres], factor(e)) - e));
+0$
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index a7fbcd8a2..df10c7ba4 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7212,6 +7212,15 @@ block([f, F],
 0$
 
 
+/* Bug #5675: "factor() and integrate() hang when %i is a common factor" */
+
+block([f, F],
+  f : %e^(3*%i*atan(a*x))/(a^2*c*x^2+c)^(3/2),
+  F : integrate(f, x),
+  if freeof('integrate, F) then radcan(logarc(exponentialize(diff(F, x) - f))) else 'fail);
+0$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 | 1 +
 src/lesfac.lisp           | 5 ++++-
 tests/rtest_gcd.mac       | 8 ++++++++
 tests/rtest_integrate.mac | 9 +++++++++
 4 files changed, 22 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  6c20745f843a721413a46a12c46641d736875c15 (commit)
      from  1080555c39e10850c6ccf10bbf67ab9e6fad5785 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 6c20745f843a721413a46a12c46641d736875c15
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 21:49:52 2026 +0200

    factor(): keep the coefficient of a common factor
    
    factor() could drop the coefficient of a factor that numerator and
    denominator have in common, such as -%i, and return a wrong result.
    With P = %i*x*sqrt(1-x)*sqrt(x+1)-x^2+1,
    block([gcd:subres], factor(expand(-%i*P)/P)) returned %i instead of
    -%i. For some similar quotients factor() looped endlessly instead.
    
    MAKPRODG, which records such a factor, kept only its main variable
    when the factor was a single term. Now it keeps the whole term.
    
    This fixes bug #5674.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 5e605357e..c521e6a30 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -218,6 +218,7 @@ Bug fixes for numbered bugs:
 * \#5666 is(abs(z) >= 0) returns false when z contains %i
 * \#5667 constantp(log(2)) returns false once log has tellsimp rules
 * \#5669 bfloat(%e^-1)^0.5b0 hangs for some values of fpprec
+* \#5674 factor() returns a wrong sign or hangs for some quotients with %i and square roots
 
 Bug fixes for unnumbered bugs:
 ------------------------------
diff --git a/src/lesfac.lisp b/src/lesfac.lisp
index 7d496f0b5..7b89fc67d 100644
--- a/src/lesfac.lisp
+++ b/src/lesfac.lisp
@@ -170,9 +170,11 @@
      (return (list g a b))))
 
 (defun makprodg (p sw)
-  (if (pcoefp p)
-      p
-      (car (makprod p sw))))
+  (cond ((pcoefp p) p)
+        ;; A monomial such as -%i needs a kernel of its own, since its main
+        ;; variable alone, which MAKPROD would give, drops the coefficient.
+        ((and (pureprod p) (not (equal p (pget (car p))))) (maksym p))
+        (t (car (makprod p sw)))))
 
 (defun dopgcdcofacts (x y)
   (let (($gcd $gcd)
diff --git a/tests/rtest_gcd.mac b/tests/rtest_gcd.mac
index 082adf577..df640b220 100644
--- a/tests/rtest_gcd.mac
+++ b/tests/rtest_gcd.mac
@@ -180,3 +180,18 @@ ratsimp(expand((y^2+1)*(x+2))/expand((y^2+x)*(x+2)));
 
 block([gcd : 'mod, x, y], gcd(expand((y + 1) * (y*x + 2)), y^2 - 1));
 y+1$
+
+
+/* Bug #5674: "factor() returns a wrong sign or hangs for some quotients with %i and square roots" */
+
+block([P : %i*x*sqrt(1-x)*sqrt(x+1)-x^2+1, gcd : subres],
+  factor(expand(-%i*P)/P));
+-%i$
+
+factor((x^2*sqrt(1-x)*sqrt(x+1)+%i*x^3-%i*x)/(%i*x*(x^2-1)*sqrt(1-x)*sqrt(x+1)-(x^2-1)^2)), gcd : subres;
+-%i*x/((x-1)*(x+1))$
+
+block([N : x^2*sqrt(1-x)*sqrt(x+1)+%i*x^3-%i*x,
+       D : %i*x*(x^2-1)*sqrt(1-x)*sqrt(x+1)-(x^2-1)^2],
+  ratsimp(block([gcd : subres], factor(1/N+1/D))*N*D/(N+D)));
+1$

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog           |  1 +
 src/lesfac.lisp     |  8 +++++---
 tests/rtest_gcd.mac | 15 +++++++++++++++
 3 files changed, 21 insertions(+), 3 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  1080555c39e10850c6ccf10bbf67ab9e6fad5785 (commit)
      from  17452499740c87a052ba5b523f82d3537de30527 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 1080555c39e10850c6ccf10bbf67ab9e6fad5785
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 19:20:02 2026 +0200

    Fix wrong antiderivatives from risch()
    
    risch() and integrate() could return an antiderivative that does not
    differentiate back to the integrand, for example
    -%e^(x/log(x))*x*log(x)^2 for %e^(x/log(x))*(x-(x+2)*log(x)). A
    linear system inside the Risch integrator was solved without checking
    that its solution has the form the system was set up for. Now that is
    checked, and when it fails another method is tried.
    
    As a side effect, risch() now leaves integrals like log(x)^n
    unevaluated instead of integrating them by parts into ones that are
    no easier, and finds x*log(x)^n for log(x)^n+n*log(x)^(n-1).
    
    This fixes bug #5273.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index fabc960a2..5e605357e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -106,6 +106,7 @@ Bug fixes for numbered bugs:
 * \#5270 risch stops when the coefficient system is short of an equation
 * \#5271 risch stops when the inner monomial's coefficients are not constants
 * \#5272 RISCHEXPVAR answers with the pair its walk reached, not the one before it
+* \#5273 risch(%e^(x/log(x))*(x-(x+2)*log(x)), x) returns a wrong antiderivative
 * \#5274 solve(tan(x) = sqrt(2)+1) -> []
 * \#5277 newcontext(foo)$ newcontext(bar)$ is(y>0)$ killcontext(foo); -> Lisp error
 * \#5278 killcontext doesn't always refuse to kill an active context
diff --git a/src/risch.lisp b/src/risch.lisp
index cabf6b3fb..0a88b0dfe 100644
--- a/src/risch.lisp
+++ b/src/risch.lisp
@@ -1267,6 +1267,16 @@
      (incf rmu)
      (unless (> rmu rrmu)
        (go l3))
+     ;; LSA solves these equations over a field that contains the variable of
+     ;; integration, but they hold only for coefficients that are constants.
+     ;; A solution that depends on that variable does not solve the
+     ;; differential equation, so walk for coefficients that are functions.
+     (unless (rzerop (r- tt
+			 (r* s risch-y)
+			 (r* r (spderivative risch-y risch-mainvar))))
+       (setq mu rrmu)
+       (setq fromlsa t)
+       (go recurse))
      (return (if (null flag)
 		 (ratqu risch-y p)
 		 (cons (r* (ratexpt (cons (list expg 1 1) 1) n) (ratqu risch-y p)) '(0))))))
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 12219e57c..a7fbcd8a2 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7204,6 +7204,14 @@ block([L],
 [0, 0, 0, 0]$
 
 
+/* Bug #5273: "risch(%e^(x/log(x))*(x-(x+2)*log(x)), x) returns a wrong antiderivative" */
+
+block([f, F],
+  f : %e^(x/log(x))*(x-(x+2)*log(x)),
+  ratsimp(diff(risch(f, x), x) - f));
+0$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/risch.lisp            | 10 ++++++++++
 tests/rtest_integrate.mac |  8 ++++++++
 3 files changed, 19 insertions(+)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  17452499740c87a052ba5b523f82d3537de30527 (commit)
      from  129972e293540043facd367bee871e8c2da822d6 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 17452499740c87a052ba5b523f82d3537de30527
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 18:57:17 2026 +0200

    Make risch() leave simpler integrals behind
    
    When risch() and integrate() cannot find a complete antiderivative,
    they may return one plus a remaining integral. For integrands such as
    z^2*%e^(c*z)/(b+a*z), the Risch integrator stopped one step short of
    what it had already computed and left a harder integral than
    necessary. Now it takes that step, unless it is the only one, which
    would amount to a plain integration by parts.
    
    Nine expected results in rtest_integrate change accordingly.
    
    This fixes bug #5272.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 6582f2f8f..fabc960a2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -105,6 +105,7 @@ Bug fixes for numbered bugs:
 * \#5269 risch loses the constant of integration in an exponential tower
 * \#5270 risch stops when the coefficient system is short of an equation
 * \#5271 risch stops when the inner monomial's coefficients are not constants
+* \#5272 RISCHEXPVAR answers with the pair its walk reached, not the one before it
 * \#5274 solve(tan(x) = sqrt(2)+1) -> []
 * \#5277 newcontext(foo)$ newcontext(bar)$ is(y>0)$ killcontext(foo); -> Lisp error
 * \#5278 killcontext doesn't always refuse to kill an active context
diff --git a/src/risch.lisp b/src/risch.lisp
index b1307f4ff..cabf6b3fb 100644
--- a/src/risch.lisp
+++ b/src/risch.lisp
@@ -833,7 +833,7 @@
 		  #'(lambda ()
 		      (ratpl (ratti risch-alphar (polcoef r risch-beta risch-var) t)
 			     (polcoef s risch-gamma risch-var))))))
-     (setq risch-y 0)
+     (setq risch-y (rzero))
    loop
      (setq yn (polcoef (ratnumerator tt) (if (zerop risch-alphar) risch-gamma k) risch-var)
       yd (r* (ratdenominator tt) ;DENOM MAY BE 0
@@ -884,13 +884,24 @@
 	       (return (list (ratqu (r* ytemp (ratexpt (cons (list expg 1 1) 1) n)) p)
 			     (erfarg2 (r* n (get expg 'rischarg))
 				      ttemp risch-ratform risch-intvar risch-mainvar))))
+	      ;; With nothing accumulated in RISCH-Y, the last step alone is a
+	      ;; single integration by parts, which usually leaves an integral no
+	      ;; easier than the one it started from. Hand that one back unevaluated
+	      ;; instead of splitting it.
+	      ((rzerop risch-y)
+	       (return (cxerfarg (rzero) expg n a risch-ratform risch-intvar
+				 risch-mainvar)))
 	      (t
+	       ;; YTEMP and TTEMP are the pair the last step reached, one
+	       ;; term further along than RISCH-Y and TT and just as
+	       ;; consistent. TTEMP is already divided by LCM: the
+	       ;; ERFARG2 clause above divides it in its own test.
 	       (return
 		 (cxerfarg
-		  (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
+		  (ratqu (r* ytemp (ratexpt (cons (list expg 1 1) 1) n)) p)
 		  expg
 		  n
-		  (ratqu tt lcm)
+		  ttemp
 		  risch-ratform
 		  risch-intvar
 		  risch-mainvar))))))
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 5a95b9df8..12219e57c 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -559,32 +559,34 @@ z*%e^(c*z)/(a*c*z+b*c)-b*'integrate(%e^(c*z)/(a^2*c*z^2+2*a*b*c*z+b^2*c),z);
 
 /* Problem 120: NOUN */
 integrate(z^2*%e^(c*z)/(b+a*z),z);
-(c*z^2-z)*%e^(c*z)/(a*c^2*z+b*c^2)-'integrate(
-                                   (b*c*z-b)*%e^(c*z)
-                                    /(a^2*c^2*z^2+2*a*b*c^2*z+b^2*c^2),z);
+b^2*'integrate(%e^(c*z)/(a^3*c*z^2+2*a^2*b*c*z+a*b^2*c),z)
+ +(a*c*z^2-a*z-b)*%e^(c*z)/(a^2*c^2*z+a*b*c^2);
 
 /* Problem 121: NOUN */
 integrate(z*%e^(c*z)/(b+a*z)^2,z);
-'integrate((a*z-b)*%e^(c*z)/(a^3*c*z^3+3*a^2*b*c*z^2+3*a*b^2*c*z+b^3*c),z)
- +z*%e^(c*z)/(a^2*c*z^2+2*a*b*c*z+b^2*c);
+(2*a-2*b*c)*'integrate(%e^(c*z)/(a^3*c^2*z^3+3*a^2*b*c^2*z^2+3*a*b^2*c^2*z
+                                 +b^3*c^2),z)
+ +(c*z+1)*%e^(c*z)/(a^2*c^2*z^2+2*a*b*c^2*z+b^2*c^2);
 
 /* Problem 122: NOUN */
 integrate(z^2*%e^(c*z)/(b+a*z)^2,z);
-z^2*%e^(c*z)/(a^2*c*z^2+2*a*b*c*z+b^2*c)
- -2*b*'integrate(z*%e^(c*z)/(a^3*c*z^3+3*a^2*b*c*z^2+3*a*b^2*c*z+b^3*c),z);
+(a*c*z^2-2*b)*%e^(c*z)/(a^3*c^2*z^2+2*a^2*b*c^2*z+a*b^2*c^2)
+ -(4*a*b-2*b^2*c)*'integrate(%e^(c*z)/(a^4*c^2*z^3+3*a^3*b*c^2*z^2
+                                       +3*a^2*b^2*c^2*z+a*b^3*c^2),z);
 
 /* Problem 123: NOUN */
 integrate(z^3*%e^(c*z)/(b+a*z)^2,z);
-'integrate((2*b^2*c*z+2*b^2)*%e^(c*z)/(a^4*c^2*z^3+3*a^3*b*c^2*z^2
-                                       +3*a^2*b^2*c^2*z+a*b^3*c^2),z) 
-  +(a*c*z^3-a*z^2-2*b*z)*%e^(c*z)/(a^3*c^2*z^2+2*a^2*b*c^2*z+a*b^2*c^2);
+(6*a*b^2-2*b^3*c)*'integrate(%e^(c*z)/(a^5*c^2*z^3+3*a^4*b*c^2*z^2
+                                       +3*a^3*b^2*c^2*z+a^2*b^3*c^2),z)
+ +(a^2*c*z^3-a^2*z^2-2*a*b*z+2*b^2)*%e^(c*z)
+  /(a^4*c^2*z^2+2*a^3*b*c^2*z+a^2*b^2*c^2);
 
 /* Problem 124: NOUN */
 integrate(z^4*%e^(c*z)/(b+a*z)^2,z);
-'integrate(((2*a*b^2*c-2*b^3*c^2)*z-2*b^3*c-4*a*b^2)*%e^(c*z)
-            /(a^5*c^3*z^3+3*a^4*b*c^3*z^2+3*a^3*b^2*c^3*z+a^2*b^3*c^3),z)
- +(a^2*c^2*z^4-2*a^2*c*z^3+(2*a^2-2*a*b*c)*z^2+(2*b^2*c+4*a*b)*z)*%e^(c*z)
-  /(a^4*c^3*z^2+2*a^3*b*c^3*z+a^2*b^2*c^3);
+(a^3*c^2*z^4-2*a^3*c*z^3+(2*a^3-2*a^2*b*c)*z^2+(2*a*b^2*c+4*a^2*b)*z-2*b^3*c
+  +2*a*b^2)*%e^(c*z)/(a^5*c^3*z^2+2*a^4*b*c^3*z+a^3*b^2*c^3)
+ -(8*a*b^3-2*b^4*c)*'integrate(%e^(c*z)/(a^6*c^2*z^3+3*a^5*b*c^2*z^2
+                                         +3*a^4*b^2*c^2*z+a^3*b^3*c^2),z);
 
 /* ----- Involving (a*z^2+b)^-n and base %e --------------------------------- */
 
@@ -594,8 +596,8 @@ integrate(%e^(c*z)/(b+a*z^2),z);
 
 /* Problem 126: NOUN */
 integrate(z*%e^(c*z)/(b+a*z^2),z);
-'integrate((a*z^2-b)*%e^(c*z)/(a^2*c*z^4+2*a*b*c*z^2+b^2*c),z)
- +z*%e^(c*z)/(a*c*z^2+b*c);
+'integrate((2*a*z-2*b*c)*%e^(c*z)/(a^2*c^2*z^4+2*a*b*c^2*z^2+b^2*c^2),z)
+ +(c*z+1)*%e^(c*z)/(a*c^2*z^2+b*c^2);
 
 /* Problem 127: NOUN */
 integrate(%e^(c*z)/(b+a*z^2)^2,z);
@@ -603,8 +605,9 @@ integrate(%e^(c*z)/(b+a*z^2)^2,z);
 
 /* Problem 128: NOUN */
 integrate(z*%e^(c*z)/(b+a*z^2)^2,z);
-'integrate((3*a*z^2-b)*%e^(c*z)/(a^3*c*z^6+3*a^2*b*c*z^4+3*a*b^2*c*z^2+b^3*c),z)
- +z*%e^(c*z)/(a^2*c*z^4+2*a*b*c*z^2+b^2*c);
+'integrate((12*a*z-4*b*c)*%e^(c*z)/(a^3*c^2*z^6+3*a^2*b*c^2*z^4
+                                    +3*a*b^2*c^2*z^2+b^3*c^2),z)
+ +(c*z+3)*%e^(c*z)/(a^2*c^2*z^4+2*a*b*c^2*z^2+b^2*c^2);
 
 /* ----- Involving (a*z^2+b*z+c)^-n ----------------------------------------- */
 
@@ -614,9 +617,10 @@ integrate(%e^(d*z)/(c+b*z+a*z^2),z);
 
 /* Problem 130: NOUN */
 integrate(z*%e^(d*z)/(c+b*z+a*z^2),z);
-'integrate((a*z^2-c)*%e^(d*z)/(a^2*d*z^4+2*a*b*d*z^3+(2*a*c*d+b^2*d)*z^2
-                                        +2*b*c*d*z+c^2*d),z)
- +z*%e^(d*z)/(a*d*z^2+b*d*z+c*d);
+'integrate(((2*a-b*d)*z-2*c*d+b)*%e^(d*z)/(a^2*d^2*z^4+2*a*b*d^2*z^3
+                                           +(2*a*c*d^2+b^2*d^2)*z^2
+                                           +2*b*c*d^2*z+c^2*d^2),z)
+ +(d*z+1)*%e^(d*z)/(a*d^2*z^2+b*d^2*z+c*d^2);
 
 /* Problem 131: NOUN */
 integrate(%e^(d*z)/(c+b*z+a*z^2)^2,z);
@@ -624,12 +628,13 @@ integrate(%e^(d*z)/(c+b*z+a*z^2)^2,z);
 
 /* Problem 132: NOUN */
 integrate(z*%e^(d*z)/(c+b*z+a*z^2)^2,z);
-'integrate((3*a*z^2+b*z-c)*%e^(d*z)/(a^3*d*z^6+3*a^2*b*d*z^5
-                                              +(3*a^2*c*d+3*a*b^2*d)*z^4
-                                              +(6*a*b*c*d+b^3*d)*z^3
-                                              +(3*a*c^2*d+3*b^2*c*d)*z^2
-                                              +3*b*c^2*d*z+c^3*d),z)
- +z*%e^(d*z)/(a^2*d*z^4+2*a*b*d*z^3+(2*a*c*d+b^2*d)*z^2+2*b*c*d*z+c^2*d);
+'integrate(((12*a-2*b*d)*z-4*c*d+6*b)*%e^(d*z)/(a^3*d^2*z^6+3*a^2*b*d^2*z^5
+                                                +(3*a^2*c*d^2+3*a*b^2*d^2)*z^4
+                                                +(6*a*b*c*d^2+b^3*d^2)*z^3
+                                                +(3*a*c^2*d^2+3*b^2*c*d^2)*z^2
+                                                +3*b*c^2*d^2*z+c^3*d^2),z)
+ +(d*z+3)*%e^(d*z)/(a^2*d^2*z^4+2*a*b*d^2*z^3+(2*a*c*d^2+b^2*d^2)*z^2
+                    +2*b*c*d^2*z+c^2*d^2);
 
 /* ----- Involving (a*z*b)^w and base %e ------------------------------------ */
 
@@ -7188,6 +7193,17 @@ risch(diff((x+1)*%e^(1/x)*%e^(x*%e^(-1/x))/x, x), x);
 (%e^(%e^-(1/x)*x+1/x)*(x+1))/x$
 
 
+/* Bug #5272: "RISCHEXPVAR answers with the pair its walk reached, not the one before it" */
+
+block([L],
+  L : [%e^(1/(x^2+1))*(x+1)/(x^2+1)^2,
+       x^2*%e^(c*x)/(b+a*x),
+       x^4*%e^(c*x)/(b+a*x)^2,
+       x*%e^(d*x)/(c+b*x+a*x^2)^2],
+  map(lambda([f], ratsimp(diff(risch(f, x), x) - f)), L));
+[0, 0, 0, 0]$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/risch.lisp            | 17 ++++++++++--
 tests/rtest_integrate.mac | 70 +++++++++++++++++++++++++++++------------------
 3 files changed, 58 insertions(+), 30 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  129972e293540043facd367bee871e8c2da822d6 (commit)
      from  f1cf5a0f463b9bd1f4212996aec2b1c9847978af (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 129972e293540043facd367bee871e8c2da822d6
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 18:37:00 2026 +0200

    Make risch() try harder on nested exponentials
    
    risch() gave up on some integrals of an exponential of an exponential,
    although they have an elementary antiderivative. For example, it left
    diff(x*%e^(1/x)*%e^(x*%e^(-1/x)), x) unevaluated instead of returning
    x*%e^(1/x)*%e^(x*%e^(-1/x)).
    
    On the way, the Risch integrator looks for an answer in a form whose
    coefficients are constants. When there is none, it gave up, although
    an answer whose coefficients depend on x may exist. Now it looks for
    that too, with a method it already uses elsewhere.
    
    This fixes bug #5271.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 7c76e3018..6582f2f8f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -104,6 +104,7 @@ Bug fixes for numbered bugs:
 * \#5268 risch gives up on nested exponential
 * \#5269 risch loses the constant of integration in an exponential tower
 * \#5270 risch stops when the coefficient system is short of an equation
+* \#5271 risch stops when the inner monomial's coefficients are not constants
 * \#5274 solve(tan(x) = sqrt(2)+1) -> []
 * \#5277 newcontext(foo)$ newcontext(bar)$ is(y>0)$ killcontext(foo); -> Lisp error
 * \#5278 killcontext doesn't always refuse to kill an active context
diff --git a/src/risch.lisp b/src/risch.lisp
index 867c74de1..b1307f4ff 100644
--- a/src/risch.lisp
+++ b/src/risch.lisp
@@ -977,7 +977,7 @@
 	 temp risch-cary risch-nogood vector aarray rmu rrmu rarray
 	 risch-beta
 	 risch-lians
-	 constmu secondpass firsty firsttt startmu starttt)
+	 constmu secondpass firsty firsttt startmu starttt fromlsa)
      (desetq (expg n eta logeta logdiff) l)
      (cond ((or (pzerop a)
 		(pzerop (car a)))
@@ -1103,7 +1103,11 @@
      (setq startmu mu)
      (setq starttt tt)
    recurseloop
-     (setq f (r+ (ratqu (polcoef s risch-gamma risch-var) rbeta)
+     ;; The equation for the coefficient of the MU-th power only involves
+     ;; the term of S of the same degree as the leading term of R. That is
+     ;; the leading term of S when both degrees agree, and zero when S has
+     ;; the lower degree, as it can when LSACALL falls back on this walk.
+     (setq f (r+ (ratqu (polcoef s risch-beta risch-var) rbeta)
 		 (if expexpflag
 		     (r* mu (spderivative eta risch-mainvar))
 		     0)))
@@ -1168,6 +1172,12 @@
        (setq tt firsttt))
      (return (cond ((null flag)
 		    nil)
+		   ;; LSACALL returns either a complete answer or none.
+		   ;; A walk entered from there does the same and leaves the integral
+		   ;; unevaluated when it does not finish.
+		   (fromlsa
+		    (cxerfarg (rzero) expg n a risch-ratform risch-intvar
+			      risch-mainvar))
 		   (t
 		    (cxerfarg (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
 			      expg
@@ -1229,11 +1239,14 @@
      (setq rarray (reverse aarray))
      (setq temp
        (lsa rarray))
+     ;; LSACALL only looks for coefficients that are constants. When LSA finds
+     ;; no single such solution, let the RECURSE walk try. It also finds
+     ;; coefficients that depend on the variable of integration.
      (when (or (eq temp 'singular)
 	       (eq temp 'inconsistent))
-       (return (if (null flag)
-		   nil
-		   (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar))))
+       (setq mu rrmu)
+       (setq fromlsa t)
+       (go recurse))
      (setq temp (reverse  (cdr temp)))
      (setq rmu 0)
      (setq risch-y 0)
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 7a5555c19..5a95b9df8 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7179,6 +7179,15 @@ risch(diff(%e^(x*%e^(-1/x))/(x^2+1), x), x);
 %e^(%e^-(1/x)*x)/(x^2+1)$
 
 
+/* Bug #5271: "risch stops when the inner monomial's coefficients are not constants" */
+
+risch(diff(x*%e^(1/x)*%e^(x*%e^(-1/x)), x), x);
+%e^(%e^-(1/x)*x+1/x)*x$
+
+risch(diff((x+1)*%e^(1/x)*%e^(x*%e^(-1/x))/x, x), x);
+(%e^(%e^-(1/x)*x+1/x)*(x+1))/x$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/risch.lisp            | 23 ++++++++++++++++++-----
 tests/rtest_integrate.mac |  9 +++++++++
 3 files changed, 28 insertions(+), 5 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  f1cf5a0f463b9bd1f4212996aec2b1c9847978af (commit)
      from  c88fd44fa3f8b6acc204724f6d5b76238373436d (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit f1cf5a0f463b9bd1f4212996aec2b1c9847978af
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 18:30:20 2026 +0200

    Fix risch() giving up on some nested exponentials
    
    risch() gave up on some integrals of an exponential of an exponential
    times a rational function, although they have an elementary
    antiderivative. For example, it left diff(%e^(x*%e^(-1/x))/(x+1), x)
    unevaluated instead of returning %e^(x*%e^(-1/x))/(x+1).
    
    On the way, the Risch integrator solves a linear system that here has
    more than one solution. Any of them would do, because they only change
    the result by a constant of integration, but risch() gave up instead.
    Now it picks one.
    
    This fixes bug #5270.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 7e1de17e9..7c76e3018 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -103,6 +103,7 @@ Bug fixes for numbered bugs:
 * \#5267 solve(x^4-5*x^2+4) vs. solve(sin(x)^4-5*sin(x)^2+4)
 * \#5268 risch gives up on nested exponential
 * \#5269 risch loses the constant of integration in an exponential tower
+* \#5270 risch stops when the coefficient system is short of an equation
 * \#5274 solve(tan(x) = sqrt(2)+1) -> []
 * \#5277 newcontext(foo)$ newcontext(bar)$ is(y>0)$ killcontext(foo); -> Lisp error
 * \#5278 killcontext doesn't always refuse to kill an active context
diff --git a/src/risch.lisp b/src/risch.lisp
index 0cf513bf4..867c74de1 100644
--- a/src/risch.lisp
+++ b/src/risch.lisp
@@ -730,7 +730,7 @@
       (setq ans (rischadd w ans)))))
 
 (defun rischexpvar (flag l risch-ratform risch-intvar risch-y risch-var risch-mainvar)
-  (prog (lcm risch-m p risch-alphar risch-gamma delta r s
+  (prog (lcm risch-m p risch-alphar risch-gamma delta r s fromlsa
 	 tt denom k wl wv i ytemp ttemp yalpha f a expg n yn yd
 	 risch-beta)
      (desetq (f a expg n) l)
@@ -792,12 +792,15 @@
 	    (go l2)))
      (multiple-value-setq (risch-y risch-m)
        (lsa wl))
+     ;; LSA calls the system singular when it has more than one solution.
+     ;; Here the antiderivatives those solutions give differ only by a
+     ;; constant, so any one of them will do, and the DOWN2 walk finds one.
+     ;; The walk answers only when it solves the equation exactly, so it is
+     ;; also safe to let it try when LSA calls the system inconsistent.
      (if (or (eq risch-y 'singular)
 	     (eq risch-y 'inconsistent))
-	 (cond ((null flag)
-		(return nil))
-	       (t
-		(return (cxerfarg (rzero) expg n a risch-ratform risch-intvar risch-mainvar)))))
+	 (progn (setq fromlsa t)
+		(go down2)))
      (setq k 0)
      (setq lcm 0)
      (setq risch-y (cdr risch-y))
@@ -868,6 +871,12 @@
 				0)))))
 	      ((null flag)
 	       (return nil))
+	      ;; The LSA branch returns either a complete answer or none. A walk
+	      ;; entered from there does the same and leaves the integral
+	      ;; unevaluated when it does not finish.
+	      (fromlsa
+	       (return (cxerfarg (rzero) expg n a risch-ratform risch-intvar
+				 risch-mainvar)))
 	      ((and (risch-constp (setq ttemp (ratqu ttemp lcm)) risch-mainvar)
 		    $erfflag
 		    (equal (pdegree (car (get expg 'rischarg)) risch-mainvar) 2)
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index cc3aac76e..7a5555c19 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7167,7 +7167,16 @@ block([f, F, dF],
   F : risch(f, x),
   dF : diff(F, x),
   if freeof('integrate, F) then ratsimp(dF - f) else 'fail);
-0;
+0$
+
+
+/* Bug #5270: "risch stops when the coefficient system is short of an equation" */
+
+risch(diff(%e^(x*%e^(-1/x))/(x+1), x), x);
+%e^(%e^-(1/x)*x)/(x+1)$
+
+risch(diff(%e^(x*%e^(-1/x))/(x^2+1), x), x);
+%e^(%e^-(1/x)*x)/(x^2+1)$
 
 
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 |  1 +
 src/risch.lisp            | 19 ++++++++++++++-----
 tests/rtest_integrate.mac | 11 ++++++++++-
 3 files changed, 25 insertions(+), 6 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  c88fd44fa3f8b6acc204724f6d5b76238373436d (commit)
      from  2e7760550b8fd49601b574d520dad492d5603333 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit c88fd44fa3f8b6acc204724f6d5b76238373436d
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 18:22:09 2026 +0200

    Let risch() integrate more nested exponentials
    
    risch() left some integrals of an exponential of an exponential partly
    unevaluated, although they have an elementary antiderivative. For
    example, it did not integrate diff((3*x-2)*%e^(x*%e^(-1/x)), x) back
    to (3*x-2)*%e^(x*%e^(-1/x)).
    
    The Risch integrator works through such a problem one step at a time.
    One step has a solution that is fixed only up to a constant, and only
    a later step can tell which constant is right. risch() kept the first
    solution it found. Now it tries a second constant and works out the
    right one from the two attempts.
    
    An integral in rtest_abs_integrate that used to stay unevaluated now
    gets a closed form, so its expected result changes.
    
    This fixes bug #5269.
    
    AI-Assisted-By: Claude Opus 5.5

diff --git a/ChangeLog b/ChangeLog
index 5832b5a96..7e1de17e9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -102,6 +102,7 @@ Bug fixes for numbered bugs:
 * \#5264 sign(a+b) -> pnz or error, depending on whether a or b is declared imaginary
 * \#5267 solve(x^4-5*x^2+4) vs. solve(sin(x)^4-5*sin(x)^2+4)
 * \#5268 risch gives up on nested exponential
+* \#5269 risch loses the constant of integration in an exponential tower
 * \#5274 solve(tan(x) = sqrt(2)+1) -> []
 * \#5277 newcontext(foo)$ newcontext(bar)$ is(y>0)$ killcontext(foo); -> Lisp error
 * \#5278 killcontext doesn't always refuse to kill an active context
diff --git a/share/contrib/integration/rtest_abs_integrate.mac b/share/contrib/integration/rtest_abs_integrate.mac
index 50cff736f..da7cd58ae 100644
--- a/share/contrib/integration/rtest_abs_integrate.mac
+++ b/share/contrib/integration/rtest_abs_integrate.mac
@@ -841,10 +841,12 @@ integrate(abs(cos(x)),x,0,%pi);
 (ii : integrate(diff((tan(x)+x)*exp(tan(x)),x),x),0);
 0$
 
-/* of course a non-nounform is better, but a nounform is better
-than an infinite loop.*/
-is(inpart(ii,0)= nounify('integrate));
-true$
+block([f, F, dF],
+  f : diff((tan(x)+x)*exp(tan(x)), x),
+  F : ii,
+  dF : diff(F, x),
+  if freeof('integrate, F) then ratsimp(exponentialize(dF - f)) else 'fail);
+0$
 
 (remvalue(ii),0);
 0$ 
diff --git a/src/risch.lisp b/src/risch.lisp
index cd1cdfc50..0cf513bf4 100644
--- a/src/risch.lisp
+++ b/src/risch.lisp
@@ -967,7 +967,8 @@
 	 mu r s tt denom ymu rbeta expg n eta logeta logdiff
 	 temp risch-cary risch-nogood vector aarray rmu rrmu rarray
 	 risch-beta
-	 risch-lians)
+	 risch-lians
+	 constmu secondpass firsty firsttt startmu starttt)
      (desetq (expg n eta logeta logdiff) l)
      (cond ((or (pzerop a)
 		(pzerop (car a)))
@@ -1090,11 +1091,19 @@
    recurse
      (setq rbeta (polcoef r risch-beta risch-var))
      (setq risch-y '(0 . 1))
+     (setq startmu mu)
+     (setq starttt tt)
    recurseloop
      (setq f (r+ (ratqu (polcoef s risch-gamma risch-var) rbeta)
 		 (if expexpflag
 		     (r* mu (spderivative eta risch-mainvar))
 		     0)))
+     ;; F is affine in MU with slope ETA', so for an exponential over an
+     ;; exponential it vanishes at no more than one MU. There the coefficient
+     ;; equation is the quadrature Y' = A, and EXPPOLYCONTROL answers with one
+     ;; of its solutions. The others differ from that one by a constant.
+     (when (and expexpflag (null constmu) (rzerop f))
+       (setq constmu mu))
      (setq ymu (exppolycontrol nil
 			       f
 			       (ratqu (polcoef (ratnumerator tt)
@@ -1109,12 +1118,9 @@
 			       risch-degree
 			       risch-mainvar))
      (when (null ymu)
-       (return (cond ((null flag)
-		      nil)
-		     (t
-		      (return (cxerfarg (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
-					expg n (ratqu tt lcm)
-					risch-ratform risch-intvar risch-mainvar))))))
+       (go recursefail))
+     (when (and secondpass (eql mu constmu))
+       (setq ymu (r+ ymu (cons 1 1))))
      (setq risch-y (r+ risch-y (setq ymu (r* ymu (pexpt (list logeta 1 1) mu)))))
      (setq tt (r- tt
 		  (r* s ymu)
@@ -1123,21 +1129,49 @@
      (cond
        ((not (< mu 0))
 	(go recurseloop))
-       ((not flag)
-	(return (cond ((rzerop tt)
-		       (ratqu risch-y p))
-		      (t nil))))
        ((rzerop tt)
-	(return (cons (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
-		      '(0))))
-       (t
-	(return (cxerfarg (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
-			  expg
-			  n
-			  (ratqu tt lcm)
-			  risch-ratform
-			  risch-intvar
-			  risch-mainvar))))
+	(go recursedone))
+       ((and constmu (null secondpass))
+	;; The constant left open at CONSTMU is fixed by one of the coefficient
+	;; equations below it, which this walk has already passed. Walk again
+	;; with that constant one higher.
+	(setq firsty risch-y)
+	(setq firsttt tt)
+	(setq secondpass t)
+	(setq risch-y '(0 . 1))
+	(setq mu startmu)
+	(setq tt starttt)
+	(go recurseloop))
+       (secondpass
+	;; Both walks keep R*Y' + S*Y + TT equal to STARTTT, so the combination
+	;; with a constant CC leaves FIRSTTT - CC*SPREAD = 0 and solves the
+	;; equation.
+	(let* ((spread (r- firsttt tt))
+	       (cc (and (not (rzerop spread))
+			(ratqu firsttt spread))))
+	  (when (and cc (risch-constp cc risch-mainvar))
+	    (setq risch-y (r+ firsty (r* cc (r- risch-y firsty))))
+	    (go recursedone)))))
+   recursefail
+     ;; A second walk that came to nothing leaves the first walk's answer.
+     (when secondpass
+       (setq risch-y firsty)
+       (setq tt firsttt))
+     (return (cond ((null flag)
+		    nil)
+		   (t
+		    (cxerfarg (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
+			      expg
+			      n
+			      (ratqu tt lcm)
+			      risch-ratform
+			      risch-intvar
+			      risch-mainvar))))
+   recursedone
+     (return (if (null flag)
+		 (ratqu risch-y p)
+		 (cons (ratqu (r* risch-y (ratexpt (cons (list expg 1 1) 1) n)) p)
+		       '(0))))
    lsacall
      (setq rrmu mu)
    muloop
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 97864921f..cc3aac76e 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7154,6 +7154,22 @@ block([x],
 true$
 
 
+/* Bug #5269: "risch loses the constant of integration in an exponential tower" */
+
+risch(diff((3*x-2)*%e^(x*%e^(-1/x)), x), x);
+%e^(%e^-(1/x)*x)*(3*x-2);
+
+risch(diff((x^2+x+7)*%e^(x*%e^(-1/x)), x), x);
+%e^(%e^-(1/x)*x)*(x^2+x+7);
+
+block([f, F, dF],
+  f : diff(%e^(%e^-(1/x)*x-1/x)*(x-%e^(1/x)), x),
+  F : risch(f, x),
+  dF : diff(F, x),
+  if freeof('integrate, F) then ratsimp(dF - f) else 'fail);
+0;
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                                         |  1 +
 share/contrib/integration/rtest_abs_integrate.mac | 10 +--
 src/risch.lisp                                    | 76 ++++++++++++++++-------
 tests/rtest_integrate.mac                         | 16 +++++
 4 files changed, 78 insertions(+), 25 deletions(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  2e7760550b8fd49601b574d520dad492d5603333 (commit)
      from  f5ee506770a79d53a46f3e9c7f8df58e865165fd (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 2e7760550b8fd49601b574d520dad492d5603333
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 11:47:10 2026 +0200

    ratout.lisp: fix PGATHERCOEF for a constant primitive part
    
    When an operand of PGCDP is univariate in the main variable, PMODCONTENT
    returns it as the content, and its primitive part is the number 1, of
    which PGATHERCOEF then took CAR. So gcd:mod gave a Lisp error, for
    example for gcd(expand((y+1)*(y*x+2)), y^2-1). PGATHERCOEF now returns 1
    for a constant, as for a polynomial in another variable.
    
    This fixes bug #5502.
    
    Found and fixed by Claude Opus 5.5.

diff --git a/ChangeLog b/ChangeLog
index 1458953c2..5832b5a96 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -158,6 +158,7 @@ Bug fixes for numbered bugs:
 * \#5490 poisexpt(sin(x), -1) never returns
 * \#5491 poissubst(w, v, cos(u-v), e, 3) gives the cos(w-u) term the wrong sign
 * \#5492 intopois(x^2.5) gives a Lisp error instead of a Poisson coefficient
+* \#5502 gcd:mod signals a Lisp error when one polynomial is univariate in the main variable
 * \#5507 resultant(x^3+1, x, x) gives 1 instead of -1: wrong sign when the second polynomial has a factor c*x^k
 * \#5508 risch() and integrate() return wrong antiderivatives for some exp(rational) integrands
 * \#5517 sign(inf^2-inf) returns neg and is(inf^2-inf>0) returns false
diff --git a/src/ratout.lisp b/src/ratout.lisp
index cc263d591..fbce42611 100644
--- a/src/ratout.lisp
+++ b/src/ratout.lisp
@@ -46,7 +46,7 @@
 	   (pgcdm ratout-x ratout-y)))))
 
 (defun pgathercoef (p ratout-*chk ratout-*res)
-  (if (not (eq (car p) ratout-*chk))
+  (if (or (pcoefp p) (not (eq (car p) ratout-*chk)))
       1
       (labels
 	  ((pgath2 (p vmax)
diff --git a/tests/rtest_gcd.mac b/tests/rtest_gcd.mac
index d7da54dbf..082adf577 100644
--- a/tests/rtest_gcd.mac
+++ b/tests/rtest_gcd.mac
@@ -174,3 +174,9 @@ true;
 
 ratsimp(expand((y^2+1)*(x+2))/expand((y^2+x)*(x+2)));
 (y^2+1)/(y^2+x)$
+
+
+/* Bug #5502: "gcd:mod signals a Lisp error when one polynomial is univariate in the main variable" */
+
+block([gcd : 'mod, x, y], gcd(expand((y + 1) * (y*x + 2)), y^2 - 1));
+y+1$

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog           | 1 +
 src/ratout.lisp     | 2 +-
 tests/rtest_gcd.mac | 6 ++++++
 3 files changed, 8 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  f5ee506770a79d53a46f3e9c7f8df58e865165fd (commit)
      from  e7ebfaacecce524febeece7ce06283877e4aef6a (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit f5ee506770a79d53a46f3e9c7f8df58e865165fd
Author: David Scherfgen <d.s...@go...>
Date:   Wed Oct 14 11:34:32 2026 +0200

    defmfun-check.lisp: fix keyword options given as strings
    
    DEFMFUN-KEYWORDS accepted an option name given as a string, since
    $VERBIFY turns it into a valid keyword, but then called KEYWORDIFY on
    the string itself, so quad_qags(x^2, x, 0, 1, "epsrel" = 1d-8) signaled
    a Lisp error. Make the keyword from the verbified name.
    
    This fixes bug #5317.
    
    Found and fixed by Claude Opus 5.5.

diff --git a/ChangeLog b/ChangeLog
index ec40738e7..1458953c2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -120,6 +120,7 @@ Bug fixes for numbered bugs:
 * \#5310 sign() gives a wrong definite sign when SIGN-SHIFT collapses a sum to a product
 * \#5311 asksign asks about, and records an answer for, a different expression after SIGN-SHIFT
 * \#5312 lambert_w is not treated as complex: csign(lambert_w(x)) gives pnz
+* \#5317 quad_qags(x^2, x, 0, 1, "epsrel" = 1d-8) signals a Lisp error
 * \#5358 inverse_erf(1.0*%i) and other complex points with abs(z) = 1 report undefined 
 * \#5359 inverse_erfc of any complex number signals a Lisp error
 * \#5360 log_gamma of a complex bigfloat with negative real part leaves log(sin(...)) unevaluated
diff --git a/src/defmfun-check.lisp b/src/defmfun-check.lisp
index a8fd56f81..e8c61b879 100644
--- a/src/defmfun-check.lisp
+++ b/src/defmfun-check.lisp
@@ -237,7 +237,7 @@
 			  (member ($verbify opt) valid-keywords))
 		      (flet ((keywordify (x)
 			       (intern (subseq (symbol-name x) 1) :keyword)))
-			(list (keywordify opt) val))
+			(list (keywordify ($verbify opt)) val))
 		      (merror (intl:gettext "~M: Unrecognized keyword: ~M")
 			      fname opt))))
 	    options)))
diff --git a/tests/rtest_integrate.mac b/tests/rtest_integrate.mac
index 36692ebf4..97864921f 100644
--- a/tests/rtest_integrate.mac
+++ b/tests/rtest_integrate.mac
@@ -7147,6 +7147,13 @@ block([F, G, n, q, u, domain : complex],
 true$
 
 
+/* Bug #5317: "quad_qags(x^2, x, 0, 1, "epsrel" = 1d-8) signals a Lisp error" */
+
+block([x],
+  is(quad_qags(x, x, 0, 1, "epsrel" = 1d-10) = quad_qags(x, x, 0, 1, 'epsrel = 1d-10)));
+true$
+
+
 /**************************************/
 /* Leave this at the end of the file! */
 /**************************************/

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                 | 1 +
 src/defmfun-check.lisp    | 2 +-
 tests/rtest_integrate.mac | 7 +++++++
 3 files changed, 9 insertions(+), 1 deletion(-)


hooks/post-receive
-- 
Maxima CAS