GitLab

at 2026-01-02T09:40:54-08:00

Address (most) review comments

 	   ;; Thus tt = (c + d*r).

 	   (cond ((>= (abs r) +rmin+)

 		  (let ((br (* b r)))

-		    (if (/= br 0)

+		    (if (/= br 0d0)

 			(/ (+ a br) tt)

 			;; b*r underflows.  Instead, compute

 			;;

 		    (multiple-value-bind (e f)

 			(robust-subinternal b a d c)

 		      (complex e (- f))))))))))

-      (let* ((a (realpart x))

-	     (b (imagpart x))

-	     (c (realpart y))

-	     (d (imagpart y))

-	     (ab (max (abs a) (abs b)))

-	     (cd (max (abs c) (abs d)))

+      (let* ((max-ab (max (abs (realpart x))

+			  (abs (imagpart x))))

+	     (max-cd (max (abs (realpart y))

+			  (abs (imagpart y))))

 	     (s 1d0))

 	(declare (double-float s))

 	;; If a or b is big, scale down a and b.

-	(when (>= ab +rbig+)

-	  (setf x (/ x 2)

-		s (* s 2)))

+	(when (>= max-ab +rbig+)

+	  (setf x (/ x 2d0)

+		s (* s 2d0)))

 	;; If c or d is big, scale down c and d.

-	(when (>= cd +rbig+)

-	  (setf y (/ y 2)

-		s (/ s 2)))

+	(when (>= max-cd +rbig+)

+	  (setf y (/ y 2d0)

+		s (/ s 2d0)))

 	;; If a or b is tiny, scale up a and b.

-	(when (<= ab (* +rmin+ +2/eps+))

+	(when (<= max-ab (* +rmin+ +2/eps+))

 	  (setf x (* x +be+)

 		s (/ s +be+)))

 	;; If c or d is tiny, scale up c and d.

-	(when (<= cd (* +rmin+ +2/eps+))

+	(when (<= max-cd (* +rmin+ +2/eps+))

 	  (setf y (* y +be+)

 		s (* s +be+)))

 	(* s

        (build-ratio (maybe-truncate nx gcd)

 		    (* (maybe-truncate y gcd) (denominator x)))))))

+

 (defun %negate (n)

   (number-dispatch ((n number))

     (((foreach fixnum single-float double-float #+long-float long-float))

             defpackage with errno symbols

     * #453: Use correct flags for analyzer and always save logs.

     * #456: Improve accuracy for division of complex double-floats

+            using Baudin and Smith's robust complex division algorithm

+            with improvements by Patrick McGehearty.

     * #458: Spurious overflow in double-double-float multiply

   * Other changes:

   * Improvements to the PCL implementation of CLOS:

 ;; Rudimentary code to read C %a formatted numbers that look like

 ;; "-0x1.c4dba4ba1ee79p-620".  We assume STRING is exactly in this

 ;; format.  No error-checking is done.

-(defun parse-%a (string)

+(defun parse-hex-float (string)

   (let* ((sign (if (char= (aref string 0) #\-)

 		   -1

 		   1))

    ;; 13

    ;; Iteration 1.  Without this, we would instead return

    ;;

-   ;;   (complex (parse-%a "0x1.ba8df8075bceep+155") (parse-%a "-0x1.a4ad6329485f0p-895"))

+   ;;   (complex (parse-hex-float "0x1.ba8df8075bceep+155")

+   ;;            (parse-hex-float "-0x1.a4ad6329485f0p-895"))

    ;;

    ;; whose imaginary part is quite a bit off.

-   (list (complex (parse-%a "0x1.73a3dac1d2f1fp+509") (parse-%a "-0x1.c4dba4ba1ee79p-620"))

-	 (complex (parse-%a "0x1.adf526c249cf0p+353") (parse-%a "0x1.98b3fbc1677bbp-697"))

-	 (complex (parse-%a "0x1.BA8DF8075BCEEp+155") (parse-%a "-0x1.A4AD628DA5B74p-895"))

+   (list (complex (parse-hex-float "0x1.73a3dac1d2f1fp+509")

+		  (parse-hex-float "-0x1.c4dba4ba1ee79p-620"))

+	 (complex (parse-hex-float "0x1.adf526c249cf0p+353")

+		  (parse-hex-float "0x1.98b3fbc1677bbp-697"))

+	 (complex (parse-hex-float "0x1.BA8DF8075BCEEp+155")

+		  (parse-hex-float "-0x1.A4AD628DA5B74p-895"))

 	 53)

    ;; 14

    ;; Iteration 2.

-   (list (complex (parse-%a "-0x0.000000008e4f8p-1022") (parse-%a "0x0.0000060366ba7p-1022"))

-	 (complex (parse-%a "-0x1.605b467369526p-245") (parse-%a "0x1.417bd33105808p-256"))

-	 (complex (parse-%a "0x1.cde593daa4ffep-810") (parse-%a "-0x1.179b9a63df6d3p-799"))

+   (list (complex (parse-hex-float "-0x0.000000008e4f8p-1022")

+		  (parse-hex-float "0x0.0000060366ba7p-1022"))

+	 (complex (parse-hex-float "-0x1.605b467369526p-245")

+		  (parse-hex-float "0x1.417bd33105808p-256"))

+	 (complex (parse-hex-float "0x1.cde593daa4ffep-810")

+		  (parse-hex-float "-0x1.179b9a63df6d3p-799"))

 	 52)

    ;; 15

    ;; Iteration 3

-   (list (complex (parse-%a "0x1.cb27eece7c585p-355 ") (parse-%a "0x0.000000223b8a8p-1022"))

-	 (complex (parse-%a "-0x1.74e7ed2b9189fp-22") (parse-%a "0x1.3d80439e9a119p-731"))

-	 (complex (parse-%a "-0x1.3b35ed806ae5ap-333") (parse-%a "-0x0.05e01bcbfd9f6p-1022"))

+   (list (complex (parse-hex-float "0x1.cb27eece7c585p-355 ")

+		  (parse-hex-float "0x0.000000223b8a8p-1022"))

+	 (complex (parse-hex-float "-0x1.74e7ed2b9189fp-22")

+		  (parse-hex-float "0x1.3d80439e9a119p-731"))

+	 (complex (parse-hex-float "-0x1.3b35ed806ae5ap-333")

+		  (parse-hex-float "-0x0.05e01bcbfd9f6p-1022"))

 	 53)

    ;; 16

    ;; Iteration 4

-   (list (complex (parse-%a "-0x1.f5c75c69829f0p-530") (parse-%a "-0x1.e73b1fde6b909p+316"))

-	 (complex (parse-%a "-0x1.ff96c3957742bp+1023") (parse-%a "0x1.5bd78c9335899p+1021"))

-	 (complex (parse-%a "-0x1.423c6ce00c73bp-710") (parse-%a "0x1.d9edcf45bcb0ep-708"))

+   (list (complex (parse-hex-float "-0x1.f5c75c69829f0p-530")

+		  (parse-hex-float "-0x1.e73b1fde6b909p+316"))

+	 (complex (parse-hex-float "-0x1.ff96c3957742bp+1023")

+		  (parse-hex-float "0x1.5bd78c9335899p+1021"))

+	 (complex (parse-hex-float "-0x1.423c6ce00c73bp-710")

+		  (parse-hex-float "0x1.d9edcf45bcb0ep-708"))

 	 52)

    ))

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl

Commits:

3 changed files:

Changes: