Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: 3416f6d5 by Raymond Toy at 2026-01-08T08:13:59-08:00 Fix #456: More accurate complex division - - - - - f8d90cd0 by Raymond Toy at 2026-01-08T08:13:59-08:00 Merge branch 'issue-456-more-accurate-complex-div' into 'master' Fix #456: More accurate complex division Closes #456 See merge request cmucl/cmucl!335 - - - - - 118bae2c by Raymond Toy at 2026-01-08T08:22:47-08:00 Merge branch 'master' into issue-459-more-accurate-dd-complex-div - - - - - 1 changed file: - tests/float.lisp Changes: ===================================== tests/float.lisp ===================================== @@ -591,3 +591,211 @@ (:tag :issues) (assert-equal -2w300 (* -2w300 1w0))) + + + +;; 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-hex-float (string) + (let* ((sign (if (char= (aref string 0) #\-) + -1 + 1)) + (dot-posn (position #\. string)) + (p-posn (position #\p string)) + (lead (parse-integer string :start (1- dot-posn) :end dot-posn)) + (frac (parse-integer string :start (1+ dot-posn) :end p-posn :radix 16)) + (exp (parse-integer string :start (1+ p-posn)))) + (* sign + (scale-float (float (+ (ash lead 52) + frac) + 1d0) + (- exp 52))))) + +;; Relative error in terms of bits of accuracy. This is the +;; definition used by Baudin and Smith. A result of 53 means the two +;; numbers have identical bits. For complex numbers, we use the min +;; of the bits of accuracy of the real and imaginary parts. +(defun rel-err (computed expected) + (flet ((rerr (c e) + (let ((diff (abs (- c e)))) + (if (zerop diff) + (float-digits diff) + (floor (- (log (/ diff (abs e)) 2d0))))))) + (min (rerr (realpart computed) (realpart expected)) + (rerr (imagpart computed) (imagpart expected))))) + +(defun do-cdiv-test (x y z-true expected-rel) + (let* ((z (/ x y)) + (rel (rel-err z z-true))) + (assert-equal expected-rel + rel + x y z z-true rel))) +;; Issue #456: improve accuracy of division of complex double-floats. +;; +;; Tests for complex division. Tests 1-10 are from Baudin and Smith. +;; Test 11 is an example from Maxima. Test 12 is an example from the +;; ansi-tests where (/ z z) didn't produce exactly 1. Tests 13-16 are +;; for examples for improvement iterations 1-4 from McGehearty. +(macrolet + ((frob (name x y z-true rel) + `(define-test ,name + (:tag :issues) + (do-cdiv-test ,x ,y ,z-true ,rel)))) + ;; First cases are from Baudin and Smith + ;; 1 + (frob cdiv.baudin-case.1 + (complex 1d0 1d0) + (complex 1d0 (scale-float 1d0 1023)) + (complex (scale-float 1d0 -1023) + (scale-float -1d0 -1023)) + 53) + ;; 2 + (frob cdiv.baudin-case.2 + (complex 1d0 1d0) + (complex (scale-float 1d0 -1023) (scale-float 1d0 -1023)) + (complex (scale-float 1d0 1023) 0) + 53) + ;; 3 + (frob cdiv.baudin-case.3 + (complex (scale-float 1d0 1023) (scale-float 1d0 -1023)) + (complex (scale-float 1d0 677) (scale-float 1d0 -677)) + (complex (scale-float 1d0 346) (scale-float -1d0 -1008)) + 53) + ;; 4 + (frob cdiv.baudin-case.4.overflow + (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) + (complex 1d0 1d0) + (complex (scale-float 1d0 1023) 0) + 53) + ;; 5 + (frob cdiv.baudin-case.5.underflow-ratio + (complex (scale-float 1d0 1020) (scale-float 1d0 -844)) + (complex (scale-float 1d0 656) (scale-float 1d0 -780)) + (complex (scale-float 1d0 364) (scale-float -1d0 -1072)) + 53) + ;; 6 + (frob cdiv.baudin-case.6.underflow-realpart + (complex (scale-float 1d0 -71) (scale-float 1d0 1021)) + (complex (scale-float 1d0 1001) (scale-float 1d0 -323)) + (complex (scale-float 1d0 -1072) (scale-float 1d0 20)) + 53) + ;; 7 + (frob cdiv.baudin-case.7.overflow-both-parts + (complex (scale-float 1d0 -347) (scale-float 1d0 -54)) + (complex (scale-float 1d0 -1037) (scale-float 1d0 -1058)) + (complex 3.898125604559113300d289 8.174961907852353577d295) + 53) + ;; 8 + (frob cdiv.baudin-case.8 + (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) + (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) + (complex 0.6d0 0.2d0) + 53) + ;; 9 + (frob cdiv.baudin-case.9 + (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) + (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) + (complex 0.001953125d0 -0.001953125d0) + 53) + ;; 10 + (frob cdiv.baudin-case.10.improve-imagpart-accuracy + (complex (scale-float 1d0 -622) (scale-float 1d0 -1071)) + (complex (scale-float 1d0 -343) (scale-float 1d0 -798)) + (complex 1.02951151789360578d-84 6.97145987515076231d-220) + 53) + ;; 11 + ;; + ;; From Maxima. This was from a (private) email where Maxima used + ;; CL:/ to compute the ratio but was not very accurate. + (frob cdiv.maxima-case + #c(5.43d-10 1.13d-100) + #c(1.2d-311 5.7d-312) + #c(3.691993880674614517999740937026568563794896024143749539711267954d301 + -1.753697093319947872394996242210428954266103103602859195409591583d301) + 52) + ;; 12 + ;; + ;; Found by ansi tests. z/z should be exactly 1. + (frob cdiv.ansi-test-z/z + #c(1.565640716292489d19 0.0d0) + #c(1.565640716292489d19 0.0d0) + #c(1d0 0) + 53) + ;; 13 + ;; Iteration 1. Without this, we would instead return + ;; + ;; (complex (parse-hex-float "0x1.ba8df8075bceep+155") + ;; (parse-hex-float "-0x1.a4ad6329485f0p-895")) + ;; + ;; whose imaginary part is quite a bit off. + (frob cdiv.mcgehearty-iteration.1 + (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. + (frob cdiv.mcgehearty-iteration.2 + (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 + (frob cdiv.mcgehearty-iteration.3 + (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 + (frob cdiv.mcgehearty-iteration.4 + (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)) + +(define-test complex-division.misc + (:tag :issue) + (let ((num '(1 + 1/2 + 1.0 + 1d0 + #c(1 2) + #c(1.0 2.0) + #c(1d0 2d0) + #c(1w0 2w0)))) + ;; Try all combinations of divisions of different types. This is + ;; primarily to test that we got all the numeric contagion cases + ;; for division in CL:/. + (dolist (x num) + (dolist (y num) + (assert-true (/ x y) + x y))))) + +(define-test complex-division.single + (:tag :issues) + (let* ((x #c(1 2)) + (y (complex (expt 2 127) (expt 2 127))) + (expected (coerce (/ x y) + '(complex single-float)))) + ;; A naive implementation of complex division would cause an + ;; overflow in computing the denominator. + (assert-equal expected + (/ (coerce x '(complex single-float)) + (coerce y '(complex single-float))) + x + y))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/78e7201838cc9b7f0de878c... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/78e7201838cc9b7f0de878c... You're receiving this email because of your account on gitlab.common-lisp.net.