[Maxima-commits] [git] Maxima CAS branch master updated. branch-5_48-base-310-g1c6a095f8

[willisbl <wil...@us...>]

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  1c6a095f833c3b3e15d13880c2476065a09bf9cf (commit)
      from  0be2b4e4d3801c78bfa69f43529bf06721e17790 (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 1c6a095f833c3b3e15d13880c2476065a09bf9cf
Author: Barton Willis <wi...@un...>
Date:   Mon Dec 8 09:57:16 2025 -0600

    Fix #4638 cabs/carg/polarform overflow/underflow and #4642 sign error
    
    - Add new function for the sign of a complex float at the top of `sign1` (compar.lisp). This avoids some floating point over and underflow problems that can happen later in the code.
    
    - In `bfloat-erf` and `complex-bfloat-erf` (gamma.lisp), do not redefine the special variable `1//2`, use `bfhalf` instead.
    
    - In `ellipt.lisp`, replace `pi` (which may be a long float in ECL) with `(float pi 1d0)`.
    
    - In `rpart.lisp`, introduce new `hypotenuse` function to handle double-float case without overflow.
    
    - Append tests from #4638 to `rtest16.mac`.
    
    - Append tests from #4642 to `rtest_sign.mac`.
    
    No unexpected failures in testsuite or share testsuite with SBCL 2.4.7 or Clozure CL 1.13.

diff --git a/src/compar.lisp b/src/compar.lisp
index 0f22d8ce5..51b2cc36d 100644
--- a/src/compar.lisp
+++ b/src/compar.lisp
@@ -1237,9 +1237,34 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
   `(multiple-value-bind (,lhs ,rhs) (compsplt ,expr)
      ,@forms))
 
+(defun sign-float (x)
+  "Maxima sign of a double-float x. Returns one of '$zero, '$pos, '$neg, or '$und (for NaN)"
+  (cond ((not (= x x)) '$und)   ; NaN check
+        ((zerop x) '$zero) ; ignores signed zero
+        ((minusp x) '$neg) ; x < 0
+        ((plusp x) '$pos)  ; x > 0
+        (t (merror "sign-float called on ~M ~%" x))))
+
+(defun sign-complex-float (x)
+  (multiple-value-bind (is-complex re im) (complex-number-p x #'floatp)
+    (cond
+      (is-complex  ;including the pure real case too
+       (cond
+         ((zerop im) (sign-float re)) ;imaginary part is zero, return sign of real part
+         ((and *complexsign* (zerop re)) '$imaginary) ; real part is zero, imaginary part nonzero, return imaginary
+         (*complexsign* '$complex) ;input is complex, return complex
+         (t (imag-err x)))) ;throw error for sign of complex
+      (t nil)))) ;x isn't a float or a complex float, return nil
+
 (defun sign1 (x)
   (setq x (specrepcheck x))
   (setq x (infsimp* x))
+
+  (let ((sgn (sign-complex-float x)))
+   (when sgn 
+    (setq sign sgn)
+     (return-from sign1 sgn))
+
   (when (and *complexsign* (atom x) (eq x '$infinity))
     ;; In Complex Mode the sign of infinity is complex.
     (when *debug-compar* (format t "~& in sign1 detect $infinity.~%"))
@@ -1278,7 +1303,7 @@ TDNEG TDZERO TDPN) to store it, and also sets SIGN."
                                    (setq sign s odds o evens e minus m)))
                            t)))))
 	 (sign exp))
-     (return sign)))
+     (return sign))))
 
 (defun numer (x)
   (let (($numer t) ; currently, no effect on $float, but proposed to
diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index 82ddc47fa..8b67ad950 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -1711,7 +1711,7 @@ first kind:
     ;;
     ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
     (let ((period (round (realpart phi) pi)))
-      (+ (base (- phi (* pi period)) m)
+      (+ (base (- phi (* (float pi 1d0) period)) m)
 	 (* 2 period (elliptic-ec m))))))
 
 ;; Complete version
diff --git a/src/gamma.lisp b/src/gamma.lisp
index f9fb3ab72..c7a0d0f2e 100644
--- a/src/gamma.lisp
+++ b/src/gamma.lisp
@@ -2241,30 +2241,29 @@
       (t
         result))))
 
+;; bfhalf is a special variable--it automatically updates when $fpprec is changed
 (defun bfloat-erf (z)
   ;; Warning!  This has round-off problems when abs(z) is very small.
-  (let ((1//2 ($bfloat '((rat simp) 1 2))))
   ;; The argument is real, the result is real too
     ($realpart
       (mul
         (simplify (list '(%signum) z))
         (sub 1
           (mul 
-            (div 1 (power ($bfloat '$%pi) 1//2))
-            (bfloat-gamma-incomplete 1//2 ($bfloat (power z 2)))))))))
+            (div 1 (power ($bfloat '$%pi) bfhalf))
+            (bfloat-gamma-incomplete bfhalf ($bfloat (power z 2))))))))
 
 (defun complex-bfloat-erf (z)
   ;; Warning!  This has round-off problems when abs(z) is very small.
   (let* (($ratprint nil)
-         (1//2 ($bfloat '((rat simp) 1 2)))
          (result
            (cmul
-             (cdiv (cpower (cpower z 2) 1//2) z)
+             (cdiv (cpower (cpower z 2) bfhalf) z)
              (sub 1
                (cmul 
-                 (div 1 (power ($bfloat '$%pi) 1//2))
+                 (div 1 (power ($bfloat '$%pi) bfhalf))
                  (complex-bfloat-gamma-incomplete 
-                   1//2
+                   bfhalf
                    ($bfloat (cpower z 2))))))))
     (cond
       ((zerop1 ($imagpart z))
diff --git a/src/rpart.lisp b/src/rpart.lisp
index 2f30fa270..09396e1f3 100644
--- a/src/rpart.lisp
+++ b/src/rpart.lisp
@@ -756,10 +756,35 @@
 	    (if absflag 0 (genatan (cdr ris) (car ris)))
 	    (cond ((equal (car ris) 0) (absarg-mabs (cdr ris)))
 		  ((equal (cdr ris) 0) (absarg-mabs (car ris)))
-		  (t (powers ($expand (add (powers (car ris) 2)
-					   (powers (cdr ris) 2))
-				      1 0)
-			     (half)))))))))
+		  (t (hypotenuse (car ris) (cdr ris)))))))))
+
+(defun hypotenuse-numerical (re im)
+ "Dispatch the CL abs function to return |re + %i im|. The inputs re and im should be floating point numbers.
+  We trust the compiler to work correctly for all double floats, including denormalized floats, and not needlessly
+  over or underflow."
+  (cond ((zerop im) (abs im))
+        ((zerop re) (abs re))
+        (t (abs (complex re im)))))
+
+(defun hypotenuse (re im)
+ (flet ((hypotenuse-default (re im) ;ok to use when no worries about floating point over/underflow
+          ;; For mixed binary64 and symbolic cases, computing re^2 or im^2 can cause a 
+          ;; floating point overflow. When an error happens, we'll punt to an abs nounform.
+
+          ;; I'd prefer to eliminate the following calls to expand, but doing so causes
+          ;; some testsuite failures.
+          (setq re ($expand re 1 0)
+                im ($expand im 1 0))
+          (or (ignore-errors (ftake 'mexpt (add (mul re re) (mul im im)) 1//2))
+              (ftake 'mabs (add re (mul '$%i im))))))
+
+    (cond ((or ($bfloatp re) ($bfloatp im)) ; at least one bigfloat
+            (hypotenuse-default ($bfloat re) ($bfloat im)))
+
+          ((or (and (floatp re) (mnump im)) (and (mnump re) (floatp im))) ;at least one float part
+             (hypotenuse-numerical ($float re) ($float im)))
+
+          (t (hypotenuse-default re im))))) ;fall back
 
 (defun genatan (num den)
   (let ((arg (take '(%atan2) num den)))
diff --git a/tests/rtest16.mac b/tests/rtest16.mac
index fc26b5c59..4d1edf279 100644
--- a/tests/rtest16.mac
+++ b/tests/rtest16.mac
@@ -3199,6 +3199,61 @@ kbateman[2](x);
 makelist(kbateman[n](0), n, 0, 5);
 [0,2/%pi,0,-2/(3*%pi), 0, 2/(5*%pi)];
 
+/* \#4642 sign(1.0e-310*%i) gives error because 1e-310*x/1e-310 fails*/
+errcatch(sign(1.0e-310*%i));
+[]$
+
+csign(1.0e-310*%i);
+imaginary$
+
+csign(1.0e-301);
+pos$
+
+csign(-1.0e-301);
+neg$
+
+csign(1.0e-301 + %i * 1.0e-301);
+complex$
+
+errcatch(sign(1.0b-310*%i));
+[]$
+
+csign(1.0b-310*%i);
+imaginary$
+
+csign(1.0b-301);
+pos$
+
+csign(-1.b-301);
+neg$
+
+csign(1.0b-301 + %i * 1.0e-301);
+complex$
+
+/* \#4638 cabs/carg/polarform overflow and underflow*/
+
+cabs(1e-170 + %i*1e-170);
+1.4142135623730951e-170$
+
+cabs(1e170 + %i*1e+170);
+1.4142135623730952e170$
+
+carg(1e170 + %i*1e+170);
+0.7853981633974483$
+
+carg(1e-310 + %i*1e-310);
+0.7853981633974483$
+
+/* Note: 1.4142135623730952e170*%e^(0.7853981633974483*%i) is not simplified! */
+(%emode : false,0);
+0$
+
+block([ans : polarform(1e170 + %i*1e+170)], [inpart(ans,0),inpart(ans,1),inpart(ans,2)]);
+["*",1.4142135623730952e170,%e^(0.7853981633974483*%i)]$
+
+(reset(%emode),0);
+0$
+
 kill(x);
 done$
 
diff --git a/tests/rtest_sign.mac b/tests/rtest_sign.mac
index acb6b01d2..4baae804f 100644
--- a/tests/rtest_sign.mac
+++ b/tests/rtest_sign.mac
@@ -1566,6 +1566,23 @@ pos$
 sign(x^(3/4));
 pz$
 
+/* \#4642 sign(1.0e-310*%i) gives error because 1e-310*x/1e-310 fails */
+errcatch(sign(1.0e-310*%i));
+[]$
+
+sign(1.0e-310);
+pos$
+
+sign(-1.0e-310);
+neg$
+
+csign(1.0e-310*%i);
+imaginary$
+
+csign(1.2e-310 + 1.0e-310*%i);
+complex$
+
+
 /* Clean up!*/
 
 (kill(all),0);

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

Summary of changes:
 src/compar.lisp      | 27 +++++++++++++++++++++++++-
 src/ellipt.lisp      |  2 +-
 src/gamma.lisp       | 13 ++++++-------
 src/rpart.lisp       | 33 +++++++++++++++++++++++++++----
 tests/rtest16.mac    | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 tests/rtest_sign.mac | 17 ++++++++++++++++
 6 files changed, 134 insertions(+), 13 deletions(-)


hooks/post-receive
-- 
Maxima CAS