Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: 78e72018 by Raymond Toy at 2026-01-07T20:30:11-08:00 Add cdiv-dd test code Remove the old test code and add new test code that fits in with the current scheme used for cdiv-double-float. - - - - - 1 changed file: - tests/float.lisp Changes: ===================================== tests/float.lisp ===================================== @@ -374,37 +374,32 @@ (min (rerr (realpart computed) (realpart expected)) (rerr (imagpart computed) (imagpart expected))))) -(define-test complex-division.double-double - (:tag :issues) - (loop for k from 1 - for test in *test-cases* - do - (destructuring-bind (x y z-true expected-rel expected-rel-w) - test - (declare (ignore expected-rel z-true)) - (flet ((compute-true (a b) - ;; Convert a and b to complex rationals, do the - ;; division and convert back to get the true - ;; expected result. - (coerce - (/ (complex (rational (realpart a)) - (rational (imagpart a))) - (complex (rational (realpart b)) - (rational (imagpart b)))) - '(complex ext:double-double-float)))) - (let* ((z (/ (coerce x '(complex ext:double-double-float)) - (coerce y '(complex ext:double-double-float)))) - (z-true (compute-true x y)) - (rel (rel-err z z-true))) - (assert-equal expected-rel-w - rel - k x y z z-true rel)))))) (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))) + +(defun do-cdiv-dd-test (x y z-true expected-rel) + (flet ((compute-true (a b) + ;; Convert a and b to complex rationals, do the + ;; division and convert back to get the true + ;; expected result. + (coerce + (/ (complex (rational (realpart a)) + (rational (imagpart a))) + (complex (rational (realpart b)) + (rational (imagpart b)))) + '(complex ext:double-double-float)))) + (let* ((z (/ (coerce x '(complex ext:double-double-float)) + (coerce y '(complex ext:double-double-float)))) + (z-true (compute-true x y)) + (rel (rel-err z z-true))) + (assert-equal expected-rel + rel + k 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. @@ -412,10 +407,11 @@ ;; 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) + ((frob (name x y z-true rel dd-rel) `(define-test ,name (:tag :issues) - (do-cdiv-test ,x ,y ,z-true ,rel)))) + (do-cdiv-test ,x ,y ,z-true ,rel) + (do-cdiv-dd-test ,x ,y ,z-true ,dd-rel)))) ;; First cases are from Baudin and Smith ;; 1 (frob cdiv.baudin-case.1 @@ -423,61 +419,71 @@ (complex 1d0 (scale-float 1d0 1023)) (complex (scale-float 1d0 -1023) (scale-float -1d0 -1023)) - 53) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 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) + 53 + 106) ;; 11 ;; ;; From Maxima. This was from a (private) email where Maxima used @@ -487,7 +493,8 @@ #c(1.2d-311 5.7d-312) #c(3.691993880674614517999740937026568563794896024143749539711267954d301 -1.753697093319947872394996242210428954266103103602859195409591583d301) - 52) + 52 + 107) ;; 12 ;; ;; Found by ansi tests. z/z should be exactly 1. @@ -495,7 +502,8 @@ #c(1.565640716292489d19 0.0d0) #c(1.565640716292489d19 0.0d0) #c(1d0 0) - 53) + 53 + 106) ;; 13 ;; Iteration 1. Without this, we would instead return ;; @@ -510,7 +518,8 @@ (parse-hex-float "0x1.98b3fbc1677bbp-697")) (complex (parse-hex-float "0x1.BA8DF8075BCEEp+155") (parse-hex-float "-0x1.A4AD628DA5B74p-895")) - 53) + 53 + 106) ;; 14 ;; Iteration 2. (frob cdiv.mcgehearty-iteration.2 @@ -520,7 +529,8 @@ (parse-hex-float "0x1.417bd33105808p-256")) (complex (parse-hex-float "0x1.cde593daa4ffep-810") (parse-hex-float "-0x1.179b9a63df6d3p-799")) - 52) + 52 + 106) ;; 15 ;; Iteration 3 (frob cdiv.mcgehearty-iteration.3 @@ -530,7 +540,8 @@ (parse-hex-float "0x1.3d80439e9a119p-731")) (complex (parse-hex-float "-0x1.3b35ed806ae5ap-333") (parse-hex-float "-0x0.05e01bcbfd9f6p-1022")) - 53) + 53 + 106) ;; 16 ;; Iteration 4 (frob cdiv.mcgehearty-iteration.4 @@ -540,7 +551,8 @@ (parse-hex-float "0x1.5bd78c9335899p+1021")) (complex (parse-hex-float "-0x1.423c6ce00c73bp-710") (parse-hex-float "0x1.d9edcf45bcb0ep-708")) - 52)) + 52 + 106)) (define-test complex-division.misc (:tag :issue) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/78e7201838cc9b7f0de878ca... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/78e7201838cc9b7f0de878ca... You're receiving this email because of your account on gitlab.common-lisp.net.