Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: ef99f009 by Raymond Toy at 2026-01-02T09:40:54-08:00 Address (most) review comments - - - - - 3 changed files: - src/code/numbers.lisp - src/general-info/release-22a.md - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -635,7 +635,7 @@ ;; 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 ;; @@ -726,28 +726,26 @@ (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 @@ -909,6 +907,7 @@ (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)) ===================================== src/general-info/release-22a.md ===================================== @@ -35,6 +35,8 @@ public domain. 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: ===================================== tests/float.lisp ===================================== @@ -354,7 +354,7 @@ ;; 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)) @@ -445,30 +445,43 @@ ;; 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) )) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/ef99f00954a34ece65ebdb8b... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/ef99f00954a34ece65ebdb8b... You're receiving this email because of your account on gitlab.common-lisp.net.