- cmucl-cvs - mailman3.common-lisp.net

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] Clean up code more
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 5b00b5d8 by Raymond Toy at 2025-12-12T07:04:17-08:00 Clean up code more Remove least-positive-double-float from *test-cases* because we don't use ulp anymore. Remove unused var `diff` for the same reason. - - - - - 1 changed file: - tests/float.lisp Changes: ===================================== tests/float.lisp ===================================== @@ -375,65 +375,65 @@ (complex 1d0 (scale-float 1d0 1023)) (complex (scale-float 1d0 -1023) (scale-float -1d0 -1023)) - 53 least-positive-double-float) + 53) ;; 2 (list (complex 1d0 1d0) (complex (scale-float 1d0 -1023) (scale-float 1d0 -1023)) (complex (scale-float 1d0 1023) 0) - 53 least-positive-double-float) + 53) ;; 3 (list (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 least-positive-double-float) + 53) ;; 4 (list (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) (complex 1d0 1d0) (complex (scale-float 1d0 1023) 0) - 53 least-positive-double-float) + 53) ;; 5 (list (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 least-positive-double-float) + 53) ;; 6 (list (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 least-positive-double-float) + 53) ;; 7 (list (complex (scale-float 1d0 -347) (scale-float 1d0 -54)) (complex (scale-float 1d0 -1037) (scale-float 1d0 -1058)) (complex 3.898125604559113300d289 8.174961907852353577d295) - 53 least-positive-double-float) + 53) ;; 8 (list (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) (complex 0.6d0 0.2d0) - 53 least-positive-double-float) + 53) ;; 9 (list (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) (complex 0.001953125d0 -0.001953125d0) - 53 least-positive-double-float) + 53) ;; 10 (list (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 least-positive-double-float) + 53) ;; 11 ;; From Maxima (list #c(5.43d-10 1.13d-100) #c(1.2d-311 5.7d-312) #c(3.691993880674614517999740937026568563794896024143749539711267954d301 -1.753697093319947872394996242210428954266103103602859195409591583d301) - 52 least-positive-double-float) + 52) ;; 12 ;; Found by ansi tests. z/z should be exactly 1. (list #c(1.565640716292489d19 0.0d0) #c(1.565640716292489d19 0.0d0) #c(1d0 0) - 53 least-positive-double-float) + 53) ;; 13 ;; Iteration 1. Without this, we would instead return ;; @@ -443,25 +443,25 @@ (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")) - 53 least-positive-double-float) + 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")) - 52 least-positive-double-float) + 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")) - 53 least-positive-double-float) + 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")) - 52 least-positive-double-float) + 52) )) ;; Relative error in terms of bits of accuracy. This is the @@ -482,12 +482,10 @@ (loop for k from 1 for test in *test-cases* do - (destructuring-bind (x y z-true expected-rel max-ulp) + (destructuring-bind (x y z-true expected-rel) test (let* ((z (/ x y)) - (diff (- z z-true)) (rel (rel-err z z-true))) - (assert-equal expected-rel rel k x y z z-true diff rel))))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/5b00b5d83514f7da9cb65aa… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/5b00b5d83514f7da9cb65aa… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] 2 commits: Add reference and comments.
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 8ca7b5cb by Raymond Toy at 2025-12-12T06:15:32-08:00 Add reference and comments. - - - - - 97e0e433 by Raymond Toy at 2025-12-12T06:44:13-08:00 Clean up impl. Remove the code for computing ulp; I don't want to use that anymore. We'll use the bits of accuracy used by Baudin and Smith. Add some comments. - - - - - 2 changed files: - src/code/numbers.lisp - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -596,6 +596,12 @@ ;; An implementation of Baudin and Smith's robust complex division for ;; double-floats. This is a pretty straightforward translation of the ;; original in https://arxiv.org/pdf/1210.4539. +;; +;; This also includes improvements mentioned in +;; https://lpc.events/event/11/contributions/1005/attachments/856/1625/Complex…. +;; In particular iteration 1 and 3 are added. Iteration 2 and 4 were +;; not added. The test examples from iteration 2 and 4 didn't change +;; with or without changes added. (let* (;; This is the value of Scilab's %eps variable. (eps (scale-float 1d0 -52)) (rmin least-positive-normalized-double-float) ===================================== tests/float.lisp ===================================== @@ -343,6 +343,9 @@ (assert-true (typep new-mode 'x86::float-modes)) (assert-equal new-mode (setf (x86::x87-floating-point-modes) new-mode)))) +;; 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) (let* ((sign (if (char= (aref string 0) #\-) -1 @@ -358,9 +361,13 @@ 1d0) (- exp 52))))) -;; Tests for complex division. From Baudin and Smith, and one test -;; from Maxima. Each test is a list of values: x, y, z-true, -;; bits-of-accuracy, and max-ulp. +;; 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. Tests 13-16 are for examples for improvement +;; iterations 1-4 from McGehearty. +;; +;; Each test is a list of values: x, y, z-true (the value of x/y), and +;; the bits of accuracy. (defparameter *test-cases* (list ;; 1 @@ -457,6 +464,10 @@ 52 least-positive-double-float) )) +;; 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)))) @@ -466,46 +477,6 @@ (min (rerr (realpart computed) (realpart expected)) (rerr (imagpart computed) (imagpart expected))))) -(defconstant +most-negative-normalized-float-exponent+ - (nth-value 1 (decode-float least-positive-normalized-double-float))) - -(defun power-of-two-p (n) - (cond ((integerp n) - (and (> n 0) - (= 0 (logand (abs n) (+ (abs n) -1))))) - ((floatp n) - (and (> n 0.0) - (= 0.5 (decode-float n)))))) - -(defun ulp (f) - (cond - ((= f 0.0) - least-positive-double-float) - (t - (multiple-value-bind (significand expon) - (decode-float f) - ;; If the exponent is smaller than the smallest - ;; normalized exponent, the ULP is the smallest float. - ;; Otherwise, the ULP has an exponent that is - ;; float-digits smaller, except when the fraction is a - ;; power of two, where we have to increase the exponent - ;; by 1. - (if (> expon +most-negative-normalized-float-exponent+) - (scale-float 0.5d0 - (- expon - (float-digits f) - (if (power-of-two-p significand) - 0 - -1))) - least-positive-double-float))))) - -(defun complex-ulp (computed expected) - (flet ((real-ulp (f) - (ulp (abs f)))) - (let ((diff (- computed expected))) - (min (real-ulp (realpart diff)) - (real-ulp (imagpart diff)))))) - (define-test complex-division.double (:tag :issues) (loop for k from 1 @@ -515,14 +486,11 @@ test (let* ((z (/ x y)) (diff (- z z-true)) - (rel (rel-err z z-true)) - (ulp (complex-ulp z z-true))) + (rel (rel-err z z-true))) (assert-equal expected-rel rel - k x y z z-true diff rel) - (assert-true (<= ulp max-ulp) - k x y z z-true diff ulp max-ulp))))) + k x y z z-true diff rel))))) (define-test complex-division.misc (:tag :issue) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3d75dda21fa51f25d6f5cf… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3d75dda21fa51f25d6f5cf… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] 7 commits: Add improvements from Patrick McGehearty
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 949dbaa2 by Raymond Toy at 2025-12-10T20:29:59-08:00 Add improvements from Patrick McGehearty Add suggestion for iteration 1 and add tests for sample cases for the other iterations. This change fixes the test case for iteration 1, which would have failed without this change. For whatever reason, the test from iteration 2 and 4 were already working to at least 52 bits of accuracy. I don't know why. Iteration 3 needs some work. - - - - - 43d732df by Raymond Toy at 2025-12-11T08:22:38-08:00 Add iteration 3 from Patrick McGehearty Very straightforward implementation. This could be cleaned up, but I'm making a checkpoint for now. The tests pass. (Well, the ULP test doesn't for test 16 because it's too tight.). - - - - - a75064f2 by Raymond Toy at 2025-12-11T17:01:59-08:00 Refactor iteration 4 code The two sets of code for |c|<=|d| or |d|<=|c| are basically identical except they're testing for c or d. Put that all into a new flet function taking c or d as the arg. - - - - - 11b159dd by Raymond Toy at 2025-12-11T17:10:24-08:00 Add rmin value to match McGehearty code - - - - - 5728d771 by Raymond Toy at 2025-12-11T17:13:25-08:00 Get rid of var ov and un in favor of rbig and rmin - - - - - 8292e387 by Raymond Toy at 2025-12-11T17:16:54-08:00 Some manual common subexpression elimination The arg for `maybe-scale` takes the absolute value so we don't have to keep taking the abs value. Also add local vars to hold the abs value of a and b so we don't have to call abs many times. - - - - - 3d75dda2 by Raymond Toy at 2025-12-11T17:29:33-08:00 Refactor some constants Move be to the list of constants. Get rid of bb which is always 2. This got absorbed into the constants be and 2/eps (new). - - - - - 2 changed files: - src/code/numbers.lisp - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -596,10 +596,15 @@ ;; An implementation of Baudin and Smith's robust complex division for ;; double-floats. This is a pretty straightforward translation of the ;; original in https://arxiv.org/pdf/1210.4539. -(let ((ov most-positive-double-float) - (un least-positive-normalized-double-float) - ;; This is the value of Scilab's %eps variable. - (eps (scale-float 1d0 -52))) +(let* (;; This is the value of Scilab's %eps variable. + (eps (scale-float 1d0 -52)) + (rmin least-positive-normalized-double-float) + (rbig (/ most-positive-double-float 2)) + (rmin2 (scale-float 1d0 -53)) + (rminscal (scale-float 1d0 51)) + (rmax2 (* rbig rmin2)) + (be (/ 2 (* eps eps))) + (2/eps (/ 2 eps))) (defun cdiv-double-float (x y) (declare (type (complex double-float) x y)) (labels @@ -623,7 +628,7 @@ (progn (format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d) (format t " r, tt = ~A ~A~%" r tt)) - (cond ((/= r 0) + (cond ((>= (abs r) rmin) (let ((br (* b r))) #+nil (format t "br = ~A~%" br) @@ -660,54 +665,84 @@ (f (internal-compreal b (- a) c d r tt))) (values e f)))) + (robust-internal (x y) (declare (type (complex double-float) x y)) (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) (d (imagpart y))) - (cond - ((<= (abs d) (abs c)) - ;; |d| <= |c|, so we can use robust-subinternal to - ;; perform the division. - (multiple-value-bind (e f) - (robust-subinternal a b c d) - (complex e f))) - (t - ;; |d| > |c|. So, instead compute - ;; - ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) - ;; - ;; Compare this to (a+i*b)/(c+i*d) and we see that - ;; realpart of the former is the same, but the - ;; imagpart of the former is the negative of the - ;; desired division. - (multiple-value-bind (e f) - (robust-subinternal b a d c) - (complex e (- f)))))))) + (flet ((maybe-scale (abs-tst) + ;; This implements McGehearty's iteration 3 to + ;; handle the case when some values are too big + ;; and should be scaled down. Also if some + ;; values are too tiny, scale them up. + (let ((abs-a (abs a)) + (abs-b (abs b))) + (if (or (> abs-tst rbig) + (> abs-a rbig) + (> abs-b rbig)) + (setf a (* a 0.5d0) + b (* b 0.5d0) + c (* c 0.5d0) + d (* d 0.5d0)) + (if (< abs-tst rmin2) + (setf a (* a rminscal) + b (* b rminscal) + c (* c rminscal) + d (* d rminscal)) + (if (or (and (< abs-a rmin) + (< abs-b rmax2) + (< abs-tst rmax2)) + (and (< abs-b rmin) + (< abs-a rmax2) + (< abs-tst rmax2))) + (setf a (* a rminscal) + b (* b rminscal) + c (* c rminscal) + d (* d rminscal)))))))) + (cond + ((<= (abs d) (abs c)) + ;; |d| <= |c|, so we can use robust-subinternal to + ;; perform the division. + (maybe-scale (abs c)) + (multiple-value-bind (e f) + (robust-subinternal a b c d) + (complex e f))) + (t + ;; |d| > |c|. So, instead compute + ;; + ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) + ;; + ;; Compare this to (a+i*b)/(c+i*d) and we see that + ;; realpart of the former is the same, but the + ;; imagpart of the former is the negative of the + ;; desired division. + (maybe-scale (abs d)) + (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))) - (bb 2d0) - (s 1d0) - (be (/ bb (* eps eps)))) + (s 1d0)) ;; If a or b is big, scale down a and b. - (when (>= ab (/ ov 2)) + (when (>= ab rbig) (setf x (/ x 2) s (* s 2))) ;; If c or d is big, scale down c and d. - (when (>= cd (/ ov 2)) + (when (>= cd rbig) (setf y (/ y 2) s (/ s 2))) ;; If a or b is tiny, scale up a and b. - (when (<= ab (* un (/ bb eps))) + (when (<= ab (* rmin 2/eps)) (setf x (* x be) s (/ s be))) ;; If c or d is tiny, scale up c and d. - (when (<= cd (* un (/ bb eps))) + (when (<= cd (* rmin 2/eps)) (setf y (* y be) s (* s be))) (* s ===================================== tests/float.lisp ===================================== @@ -343,6 +343,21 @@ (assert-true (typep new-mode 'x86::float-modes)) (assert-equal new-mode (setf (x86::x87-floating-point-modes) new-mode)))) +(defun parse-%a (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))))) + ;; Tests for complex division. From Baudin and Smith, and one test ;; from Maxima. Each test is a list of values: x, y, z-true, ;; bits-of-accuracy, and max-ulp. @@ -412,6 +427,34 @@ #c(1.565640716292489d19 0.0d0) #c(1d0 0) 53 least-positive-double-float) + ;; 13 + ;; Iteration 1. Without this, we would instead return + ;; + ;; (complex (parse-%a "0x1.ba8df8075bceep+155") (parse-%a "-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")) + 53 least-positive-double-float) + ;; 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")) + 52 least-positive-double-float) + ;; 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")) + 53 least-positive-double-float) + ;; 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")) + 52 least-positive-double-float) )) (defun rel-err (computed expected) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/1aa7c11f6cf9da0be58523… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/1aa7c11f6cf9da0be58523… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div-improve] 6 commits: Add iteration 3 from Patrick McGehearty
by Raymond Toy (＠rtoy) 12 Dec '25

12 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div-improve at cmucl / cmucl Commits: 43d732df by Raymond Toy at 2025-12-11T08:22:38-08:00 Add iteration 3 from Patrick McGehearty Very straightforward implementation. This could be cleaned up, but I'm making a checkpoint for now. The tests pass. (Well, the ULP test doesn't for test 16 because it's too tight.). - - - - - a75064f2 by Raymond Toy at 2025-12-11T17:01:59-08:00 Refactor iteration 4 code The two sets of code for |c|<=|d| or |d|<=|c| are basically identical except they're testing for c or d. Put that all into a new flet function taking c or d as the arg. - - - - - 11b159dd by Raymond Toy at 2025-12-11T17:10:24-08:00 Add rmin value to match McGehearty code - - - - - 5728d771 by Raymond Toy at 2025-12-11T17:13:25-08:00 Get rid of var ov and un in favor of rbig and rmin - - - - - 8292e387 by Raymond Toy at 2025-12-11T17:16:54-08:00 Some manual common subexpression elimination The arg for `maybe-scale` takes the absolute value so we don't have to keep taking the abs value. Also add local vars to hold the abs value of a and b so we don't have to call abs many times. - - - - - 3d75dda2 by Raymond Toy at 2025-12-11T17:29:33-08:00 Refactor some constants Move be to the list of constants. Get rid of bb which is always 2. This got absorbed into the constants be and 2/eps (new). - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -596,10 +596,15 @@ ;; An implementation of Baudin and Smith's robust complex division for ;; double-floats. This is a pretty straightforward translation of the ;; original in https://arxiv.org/pdf/1210.4539. -(let ((ov most-positive-double-float) - (un least-positive-normalized-double-float) - ;; This is the value of Scilab's %eps variable. - (eps (scale-float 1d0 -52))) +(let* (;; This is the value of Scilab's %eps variable. + (eps (scale-float 1d0 -52)) + (rmin least-positive-normalized-double-float) + (rbig (/ most-positive-double-float 2)) + (rmin2 (scale-float 1d0 -53)) + (rminscal (scale-float 1d0 51)) + (rmax2 (* rbig rmin2)) + (be (/ 2 (* eps eps))) + (2/eps (/ 2 eps))) (defun cdiv-double-float (x y) (declare (type (complex double-float) x y)) (labels @@ -623,7 +628,7 @@ (progn (format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d) (format t " r, tt = ~A ~A~%" r tt)) - (cond ((>= (abs r) un) + (cond ((>= (abs r) rmin) (let ((br (* b r))) #+nil (format t "br = ~A~%" br) @@ -660,54 +665,84 @@ (f (internal-compreal b (- a) c d r tt))) (values e f)))) + (robust-internal (x y) (declare (type (complex double-float) x y)) (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) (d (imagpart y))) - (cond - ((<= (abs d) (abs c)) - ;; |d| <= |c|, so we can use robust-subinternal to - ;; perform the division. - (multiple-value-bind (e f) - (robust-subinternal a b c d) - (complex e f))) - (t - ;; |d| > |c|. So, instead compute - ;; - ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) - ;; - ;; Compare this to (a+i*b)/(c+i*d) and we see that - ;; realpart of the former is the same, but the - ;; imagpart of the former is the negative of the - ;; desired division. - (multiple-value-bind (e f) - (robust-subinternal b a d c) - (complex e (- f)))))))) + (flet ((maybe-scale (abs-tst) + ;; This implements McGehearty's iteration 3 to + ;; handle the case when some values are too big + ;; and should be scaled down. Also if some + ;; values are too tiny, scale them up. + (let ((abs-a (abs a)) + (abs-b (abs b))) + (if (or (> abs-tst rbig) + (> abs-a rbig) + (> abs-b rbig)) + (setf a (* a 0.5d0) + b (* b 0.5d0) + c (* c 0.5d0) + d (* d 0.5d0)) + (if (< abs-tst rmin2) + (setf a (* a rminscal) + b (* b rminscal) + c (* c rminscal) + d (* d rminscal)) + (if (or (and (< abs-a rmin) + (< abs-b rmax2) + (< abs-tst rmax2)) + (and (< abs-b rmin) + (< abs-a rmax2) + (< abs-tst rmax2))) + (setf a (* a rminscal) + b (* b rminscal) + c (* c rminscal) + d (* d rminscal)))))))) + (cond + ((<= (abs d) (abs c)) + ;; |d| <= |c|, so we can use robust-subinternal to + ;; perform the division. + (maybe-scale (abs c)) + (multiple-value-bind (e f) + (robust-subinternal a b c d) + (complex e f))) + (t + ;; |d| > |c|. So, instead compute + ;; + ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) + ;; + ;; Compare this to (a+i*b)/(c+i*d) and we see that + ;; realpart of the former is the same, but the + ;; imagpart of the former is the negative of the + ;; desired division. + (maybe-scale (abs d)) + (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))) - (bb 2d0) - (s 1d0) - (be (/ bb (* eps eps)))) + (s 1d0)) ;; If a or b is big, scale down a and b. - (when (>= ab (/ ov 2)) + (when (>= ab rbig) (setf x (/ x 2) s (* s 2))) ;; If c or d is big, scale down c and d. - (when (>= cd (/ ov 2)) + (when (>= cd rbig) (setf y (/ y 2) s (/ s 2))) ;; If a or b is tiny, scale up a and b. - (when (<= ab (* un (/ bb eps))) + (when (<= ab (* rmin 2/eps)) (setf x (* x be) s (/ s be))) ;; If c or d is tiny, scale up c and d. - (when (<= cd (* un (/ bb eps))) + (when (<= cd (* rmin 2/eps)) (setf y (* y be) s (* s be))) (* s View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/949dbaa242c2a77687bcea… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/949dbaa242c2a77687bcea… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl] Pushed new branch issue-456-more-accurate-complex-div-improve
by Raymond Toy (＠rtoy) 11 Dec '25

11 Dec '25

Raymond Toy pushed new branch issue-456-more-accurate-complex-div-improve at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-456-more-accurate-c… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] 3 commits: Add comments and remove debugging.
by Raymond Toy (＠rtoy) 11 Dec '25

11 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 22b67f1e by Raymond Toy at 2025-12-10T17:48:16-08:00 Add comments and remove debugging. Add comments to cdiv-double-float based on the text of the paper and also add some derivations to make it a bit easier to follow what's going on. Remove the debugging prints too. - - - - - 7f08a0de by Raymond Toy at 2025-12-10T18:19:38-08:00 Fix round-off error for z/z not returning 1. In the ansi tests, there's a test case for (/ z z) for many random complex numbers z. This is expected to return #c(1 0), but we weren't. This is caused by the cdiv-double-float algorithm computing tt = 1/(c + d*r) and then multiplying by tt in various places. However, this causes an extra rounding operation. Instead, modify the code to use tt = (c + d*r) and divide by tt instead of multiply. This gets rid of the extra rounding. - - - - - 1aa7c11f by Raymond Toy at 2025-12-10T19:31:06-08:00 Update tests. Test 8 now has full accuracy, but test 11 from Maxima only has 52 bits. Add a test from the failed ansi-test that was testing (/ z z) being exactly 1. Fix typo in complex-division.single. - - - - - 2 changed files: - src/code/numbers.lisp - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -605,42 +605,85 @@ (labels ((internal-compreal (a b c d r tt) (declare (double-float a b c d r tt)) + ;; Compute the real part of the complex division + ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r). + ;; + ;; The realpart is (a*c+b*d)/(c^2+d^2). + ;; + ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r) + ;; + ;; Then + ;; + ;; (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r)) + ;; = (a + b*d/c)/(c+d*r) + ;; = (a + b*r)/(c + d*r). + ;; + ;; Thus tt = (c + d*r). + #+nil + (progn + (format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d) + (format t " r, tt = ~A ~A~%" r tt)) (cond ((/= r 0) (let ((br (* b r))) + #+nil + (format t "br = ~A~%" br) (if (/= br 0) - (* (+ a br) tt) - (+ (* a tt) - (* (* b tt) + (/ (+ a br) tt) + ;; b*r underflows. Instead, compute + ;; + ;; (a + b*r)*tt = a*tt + b*tt*r + ;; = a*tt + (b*tt)*r + ;; (a + b*r)/tt = a/tt + b/tt*r + ;; = a*tt + (b*tt)*r + (+ (/ a tt) + (* (/ b tt) r))))) (t - (* (+ a (* d (/ b c))) + ;; r = 0 so d is very tiny compared to c. + ;; + ;; (a + b*r)/tt = (a + b*(d/c))/tt + #+nil + (progn + (format t "r = 0~%") + (format t "a*tt = ~A~%" (* a tt))) + (/ (+ a (* d (/ b c))) tt)))) (robust-subinternal (a b c d) (declare (double-float a b c d)) (let* ((r (/ d c)) - (tt (/ (+ c (* d r))))) + (tt (+ c (* d r)))) + ;; e is the real part and f is the imaginary part. We + ;; can use internal-compreal for the imaginary part by + ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is + ;; the same as the real part of (b-i*a)/(c+i*d). (let ((e (internal-compreal a b c d r tt)) (f (internal-compreal b (- a) c d r tt))) - #+nil(format t "subint: e f = ~A ~A~%" e f) (values e f)))) (robust-internal (x y) (declare (type (complex double-float) x y)) - #+nil(format t "x = ~A, y = ~A~%" x y) (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) (d (imagpart y))) (cond ((<= (abs d) (abs c)) + ;; |d| <= |c|, so we can use robust-subinternal to + ;; perform the division. (multiple-value-bind (e f) (robust-subinternal a b c d) - #+nil(format t "1: e f = ~A ~A~%" e f) (complex e f))) (t + ;; |d| > |c|. So, instead compute + ;; + ;; (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2) + ;; + ;; Compare this to (a+i*b)/(c+i*d) and we see that + ;; realpart of the former is the same, but the + ;; imagpart of the former is the negative of the + ;; desired division. (multiple-value-bind (e f) (robust-subinternal b a d c) - #+nil(format t "2: e f = ~A ~A~%" e f) (complex e (- f)))))))) (let* ((a (realpart x)) (b (imagpart x)) @@ -648,19 +691,23 @@ (d (imagpart y)) (ab (max (abs a) (abs b))) (cd (max (abs c) (abs d))) - (b 2d0) + (bb 2d0) (s 1d0) - (be (/ b (* eps eps)))) + (be (/ bb (* eps eps)))) + ;; If a or b is big, scale down a and b. (when (>= ab (/ ov 2)) (setf x (/ x 2) s (* s 2))) + ;; If c or d is big, scale down c and d. (when (>= cd (/ ov 2)) (setf y (/ y 2) s (/ s 2))) - (when (<= ab (* un (/ b eps))) + ;; If a or b is tiny, scale up a and b. + (when (<= ab (* un (/ bb eps))) (setf x (* x be) s (/ s be))) - (when (<= cd (* un (/ b eps))) + ;; If c or d is tiny, scale up c and d. + (when (<= cd (* un (/ bb eps))) (setf y (* y be) s (* s be))) (* s ===================================== tests/float.lisp ===================================== @@ -388,7 +388,7 @@ (list (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) (complex 0.6d0 0.2d0) - 52 least-positive-double-float) + 53 least-positive-double-float) ;; 9 (list (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) @@ -399,11 +399,18 @@ (complex (scale-float 1d0 -343) (scale-float 1d0 -798)) (complex 1.02951151789360578d-84 6.97145987515076231d-220) 53 least-positive-double-float) + ;; 11 ;; From Maxima (list #c(5.43d-10 1.13d-100) #c(1.2d-311 5.7d-312) #c(3.691993880674614517999740937026568563794896024143749539711267954d301 -1.753697093319947872394996242210428954266103103602859195409591583d301) + 52 least-positive-double-float) + ;; 12 + ;; Found by ansi tests. z/z should be exactly 1. + (list #c(1.565640716292489d19 0.0d0) + #c(1.565640716292489d19 0.0d0) + #c(1d0 0) 53 least-positive-double-float) )) @@ -494,10 +501,10 @@ (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)))) + (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 View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/bb5c9e106d89594a1bc302… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/bb5c9e106d89594a1bc302… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] Clean up removing old code.
by Raymond Toy (＠rtoy) 10 Dec '25

10 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: bb5c9e10 by Raymond Toy at 2025-12-10T15:31:47-08:00 Clean up removing old code. - - - - - 2 changed files: - src/code/numbers.lisp - src/compiler/float-tran.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -593,180 +593,18 @@ (build-ratio x y) (build-ratio (truncate x gcd) (truncate y gcd)))))) - -#+nil -(defun two-arg-/ (x y) - (number-dispatch ((x number) (y number)) - (float-contagion / x y (ratio integer)) - - ((complex complex) - (let* ((rx (realpart x)) - (ix (imagpart x)) - (ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (+ ry (* r iy)))) - (canonical-complex (/ (+ rx (* ix r)) dn) - (/ (- ix (* rx r)) dn))) - (let* ((r (/ ry iy)) - (dn (+ iy (* r ry)))) - (canonical-complex (/ (+ (* rx r) ix) dn) - (/ (- (* ix r) rx) dn)))))) - (((foreach integer ratio single-float double-float) complex) - (let* ((ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (* ry (+ 1 (* r r))))) - (canonical-complex (/ x dn) - (/ (- (* x r)) dn))) - (let* ((r (/ ry iy)) - (dn (* iy (+ 1 (* r r))))) - (canonical-complex (/ (* x r) dn) - (/ (- x) dn)))))) - ((complex (or rational float)) - (canonical-complex (/ (realpart x) y) - (/ (imagpart x) y))) - - ((ratio ratio) - (let* ((nx (numerator x)) - (dx (denominator x)) - (ny (numerator y)) - (dy (denominator y)) - (g1 (gcd nx ny)) - (g2 (gcd dx dy))) - (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2)) - (* (maybe-truncate dx g2) (maybe-truncate ny g1))))) - - ((integer integer) - (integer-/-integer x y)) - - ((integer ratio) - (if (zerop x) - 0 - (let* ((ny (numerator y)) - (dy (denominator y)) - (gcd (gcd x ny))) - (build-ratio (* (maybe-truncate x gcd) dy) - (maybe-truncate ny gcd))))) - - ((ratio integer) - (let* ((nx (numerator x)) - (gcd (gcd nx y))) - (build-ratio (maybe-truncate nx gcd) - (* (maybe-truncate y gcd) (denominator x))))))) - -;; An implementation of Baudin and Smith's robust complex division. -;; This is a pretty straightforward translation of the original in -;; https://arxiv.org/pdf/1210.4539. -#+nil -(eval-when (:compile-toplevel :load-toplevel :execute) -(macrolet - ((frob (type max-float min-float eps-float) - (let ((name (symbolicate "CDIV-" type))) - `(progn - (let ((ov ,max-float) - (un ,min-float) - (eps ,eps-float)) - (defun ,name (x y) - (declare (type (complex ,type) x y)) - (labels - ((internal-compreal (a b c d r tt) - (declare (double-float a b c d r tt)) - ;; Iteration 1: compare r against DBL_MIN instead of 0 - (cond ((/= r 0) - (let ((br (* b r))) - (if (/= br 0) - (* (+ a br) tt) - (+ (* a tt) - (* (* b tt) - r))))) - (t - (* (+ a (* d (/ b c))) - tt)))) - - (robust-subinternal (a b c d) - (declare (double-float a b c d)) - (let* ((r (/ d c)) - (tt (/ (+ c (* d r))))) - (let ((e (internal-compreal a b c d r tt)) - (f (internal-compreal b (- a) c d r tt))) - #+nil(format t "subint: e f = ~A ~A~%" e f) - (values e - f)))) - - (robust-internal (x y) - (declare (type (complex ,type) x y)) - #+nil(format t "x = ~A, y = ~A~%" x y) - (let ((a (realpart x)) - (b (imagpart x)) - (c (realpart y)) - (d (imagpart y))) - (cond - ((<= (abs d) (abs c)) - (multiple-value-bind (e f) - (robust-subinternal a b c d) - #+nil(format t "1: e f = ~A ~A~%" e f) - (complex e f))) - (t - (multiple-value-bind (e f) - (robust-subinternal b a d c) - #+nil(format t "2: e f = ~A ~A~%" e f) - (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))) - (b 2d0) - (s 1d0) - (be (/ b (* eps eps)))) - (when (>= ab (/ ov 2)) - (setf x (/ x 2) - s (* s 2))) - (when (>= cd (/ ov 2)) - (setf y (/ y 2) - s (/ s 2))) - (when (<= ab (* un (/ b eps))) - (setf x (* x be) - s (/ s be))) - (when (<= cd (* un (/ b eps))) - (setf y (* y be) - s (* s be))) - ;; Iteration 2 - #+nil - (when (< cd eps) - (setf x (/ x eps)) - (setf y (/ y eps))) - - ;; Iteration 4 - #+nil - (when (> cd (/ ov 2)) - (setf x (/ x 2) - y (/ y 2))) - - (* s - (robust-internal x y)))))))))) - ;; The eps value for double-float is determined by looking at the - ;; value of %eps in Scilab. - (frob double-float most-positive-double-float least-positive-normalized-double-float - (scale-float 1d0 -52)) - ;; The eps value here is a guess. - (frob single-float most-positive-single-float least-positive-normalized-single-float - (scale-float 1f0 -22))) -) - +;; An implementation of Baudin and Smith's robust complex division for +;; double-floats. This is a pretty straightforward translation of the +;; original in https://arxiv.org/pdf/1210.4539. (let ((ov most-positive-double-float) (un least-positive-normalized-double-float) + ;; This is the value of Scilab's %eps variable. (eps (scale-float 1d0 -52))) (defun cdiv-double-float (x y) (declare (type (complex double-float) x y)) (labels ((internal-compreal (a b c d r tt) (declare (double-float a b c d r tt)) - ;; Iteration 1: compare r against DBL_MIN instead of 0 (cond ((/= r 0) (let ((br (* b r))) (if (/= br 0) @@ -777,7 +615,6 @@ (t (* (+ a (* d (/ b c))) tt)))) - (robust-subinternal (a b c d) (declare (double-float a b c d)) (let* ((r (/ d c)) @@ -787,7 +624,6 @@ #+nil(format t "subint: e f = ~A ~A~%" e f) (values e f)))) - (robust-internal (x y) (declare (type (complex double-float) x y)) #+nil(format t "x = ~A, y = ~A~%" x y) @@ -827,21 +663,11 @@ (when (<= cd (* un (/ b eps))) (setf y (* y be) s (* s be))) - ;; Iteration 2 - #+nil - (when (< cd eps) - (setf x (/ x eps)) - (setf y (/ y eps))) - - ;; Iteration 4 - #+nil - (when (> cd (/ ov 2)) - (setf x (/ x 2) - y (/ y 2))) - (* s (robust-internal x y)))))) +;; Smith's algorithm for complex division for (complex single-float). +;; We convert the parts to double-floats before computing the result. (defun cdiv-single-float (x y) (declare (type (complex single-float) x y)) (let ((a (float (realpart x) 1d0)) @@ -877,12 +703,15 @@ (cdiv-single-float (coerce x '(complex single-float)) y)) - (((foreach integer ratio single-float double-float (complex rational) (complex single-float)) + (((foreach integer ratio single-float double-float (complex rational) + (complex single-float)) (complex double-float)) (cdiv-double-float (coerce x '(complex double-float)) y)) (((complex double-double-float) - (foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float))) + (foreach (complex rational) (complex single-float) (complex double-float) + (complex double-double-float))) + ;; We should do something better for double-double floats. (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) @@ -935,6 +764,8 @@ (/ (- x) dn)))))) (((complex rational) (complex rational)) + ;; We probably don't need to do Smith's algorithm for rationals. + ;; A naive implementation of coplex division has no issues. (let ((a (realpart x)) (b (imagpart x)) (c (realpart y)) @@ -952,7 +783,8 @@ (f (/ (- b (* a r)) denom))) (canonical-complex e f)))))) - (((foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float)) + (((foreach (complex rational) (complex single-float) (complex double-float) + (complex double-double-float)) (or rational float)) (canonical-complex (/ (realpart x) y) (/ (imagpart x) y))) @@ -962,108 +794,6 @@ (cdiv-double-float (coerce x '(complex double-float)) (coerce y '(complex double-float)))) - #+nil - (((foreach integer ratio single-float double-float) complex) - (let* ((ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (* ry (+ 1 (* r r))))) - (canonical-complex (/ x dn) - (/ (- (* x r)) dn))) - (let* ((r (/ ry iy)) - (dn (* iy (+ 1 (* r r))))) - (canonical-complex (/ (* x r) dn) - (/ (- x) dn)))))) - #+nil - ((complex (or rational float)) - (canonical-complex (/ (realpart x) y) - (/ (imagpart x) y))) - - ((ratio ratio) - (let* ((nx (numerator x)) - (dx (denominator x)) - (ny (numerator y)) - (dy (denominator y)) - (g1 (gcd nx ny)) - (g2 (gcd dx dy))) - (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2)) - (* (maybe-truncate dx g2) (maybe-truncate ny g1))))) - - ((integer integer) - (integer-/-integer x y)) - - ((integer ratio) - (if (zerop x) - 0 - (let* ((ny (numerator y)) - (dy (denominator y)) - (gcd (gcd x ny))) - (build-ratio (* (maybe-truncate x gcd) dy) - (maybe-truncate ny gcd))))) - - ((ratio integer) - (let* ((nx (numerator x)) - (gcd (gcd nx y))) - (build-ratio (maybe-truncate nx gcd) - (* (maybe-truncate y gcd) (denominator x))))))) - - -#+nil -(defun two-arg-/ (x y) - (number-dispatch ((x number) (y number)) - (float-contagion / x y (ratio integer)) - - (((complex double-float) (complex double-float)) - (cdiv-double-float x y)) - (((complex double-float) (complex single-float)) - (cdiv-double-float x (coerce y '(complex double-float)))) - - (((complex single-float) (complex double-float)) - (cdiv-double-float (coerce y '(complex double-float)) - y)) - (((complex single-float) (complex single-float)) - (cdiv-single-float x y)) - - #+nil - (((foreach (complex single-float) (complex double-float)) - (complex double-float)) - (cdiv-double-float (complex (float (realpart x) 1d0) - (float (imagpart x) 1d0)) - y)) - (((foreach (complex rational) (complex double-double-float)) - (foreach (complex rational) (complex double-double-float))) - (let* ((rx (realpart x)) - (ix (imagpart x)) - (ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (+ ry (* r iy)))) - (canonical-complex (/ (+ rx (* ix r)) dn) - (/ (- ix (* rx r)) dn))) - (let* ((r (/ ry iy)) - (dn (+ iy (* r ry)))) - (canonical-complex (/ (+ (* rx r) ix) dn) - (/ (- (* ix r) rx) dn)))))) - (((foreach integer ratio single-float double-float) - (foreach (complex rational) (complex double-double-float))) - (let* ((ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (* ry (+ 1 (* r r))))) - (canonical-complex (/ x dn) - (/ (- (* x r)) dn))) - (let* ((r (/ ry iy)) - (dn (* iy (+ 1 (* r r))))) - (canonical-complex (/ (* x r) dn) - (/ (- x) dn)))))) - (((foreach (complex rational) (complex double-double-float)) - (or rational float)) - (canonical-complex (/ (realpart x) y) - (/ (imagpart x) y))) - ((ratio ratio) (let* ((nx (numerator x)) (dx (denominator x)) @@ -1092,7 +822,6 @@ (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/compiler/float-tran.lisp ===================================== @@ -1804,119 +1804,6 @@ #+double-double (frob double-double-float)) -#+(and nil sse2 complex-fp-vops) -(macrolet - ((frob (type one) - `(deftransform / ((x y) (,type ,type) * - :policy (> speed space)) - ;; Divide a complex by a complex - - ;; Here's how we do a complex division - ;; - ;; Compute (xr + i*xi)/(yr + i*yi) - ;; - ;; Assume |yi| < |yr|. Then - ;; - ;; (xr + i*xi) (xr + i*xi) - ;; ----------- = ----------------- - ;; (yr + i*yi) yr*(1 + i*(yi/yr)) - ;; - ;; (xr + i*xi)*(1 - i*(yi/yr)) - ;; = --------------------------- - ;; yr*(1 + (yi/yr)^2) - ;; - ;; (xr + i*xi)*(1 - i*(yi/yr)) - ;; = --------------------------- - ;; yr + (yi/yr)*yi - ;; - ;; This allows us to use a fast complex multiply followed by - ;; a real division. - '(let* ((ry (realpart y)) - (iy (imagpart y))) - (if (> (abs ry) (abs iy)) - (let* ((r (/ iy ry)) - (dn (+ ry (* r iy)))) - (/ (* x (complex ,one r)) - dn)) - (let* ((r (/ ry iy)) - (dn (+ iy (* r ry)))) - (/ (* x (complex r ,(- one))) - dn))))))) - (frob (complex single-float) 1f0) - (frob (complex double-float) 1d0)) - -#+(and nil sse2 complex-fp-vops) -(macrolet - ((frob (type one) - `(deftransform / ((x y) (,type ,type) * - :policy (> speed space)) - ;; Divide a complex by a complex - - ;; Here's how we do a complex division - ;; - ;; Compute (xr + i*xi)/(yr + i*yi) - ;; - ;; Assume |xi| < |xr| and |yi| < |yr|. Then - ;; - ;; (xr + i*xi) xr*(1 + i*(xi/xr)) - ;; ----------- = ----------------- - ;; (yr + i*yi) yr*(1 + i*(yi/yr)) - ;; - ;; xr*(1+i(xi/xr))*(1-i*(yi/yr)) - ;; = ----------------------------- - ;; yr*(1+(yi/yr)^2) - ;; - ;; xr*(1+i(xi/xr))*(1-i*(yi/yr)) - ;; = ----------------------------- - ;; yr + (yi/yr)*yi - ;; - ;; This allows us to use a fast complex multiply followed by - ;; a real division. - '(let* ((rx (realpart x)) - (ix (imagpart x)) - (ry (realpart y)) - (iy (imagpart y))) - (if (>= (abs rx) (abs ix)) - (if (>= (abs ry) (abs iy)) - ;; |xr| >= |xi| and |yr| >= |yi| - ;; - ;; xr*(1+i*(xi/xr))*(1-i*(yi/yr)) - ;; w = ------------------------------ - ;; yr + (yi/yr)*yi - ;; - ;; (1+i*p)*(1-i*q) = 1+p*q + i*(p-q) - (let* ((x-div (/ ix rx)) - (y-div (/ iy ry)) - (dn (+ ry (* y-div iy))) - (factor (/ xr dn))) - (* factor - (complex (1+ (* x-div y-div)) - (- x-div y-div)))) - ;; |xr| >= |xi| and |yr| < |yi| - ;; - ;; xr*(1+i*(xi/xr))*(yr/yi-i) - ;; w = ------------------------------ - ;; yi*((yr/yi)^2+1) - ;; - ;; xr*(1+i*(xi/xr))*(yr/yi-i) - ;; w = ------------------------------ - ;; (yr/yi)*yi + yi - ;; - ;; (1+i*p)*(q-i) = p+q + i*(p*q-1) - (let* ((x-div (/ ix rx)) - (y-div (/ ry iy)) - (dn (+ yi (* y-div iy))) - (factor (/ xr dn))) - (* factor - (complex (+ x-div y-div) - (1- (* x-div y-div)))))) - (let* ((r (/ ry iy)) - (dn (+ iy (* r ry)))) - (/ (* x (complex r ,(- one))) - dn))))))) - (frob (complex single-float) 1f0) - (frob (complex double-float) 1d0)) - (macrolet ((frob (type name) `(deftransform / ((x y) ((complex ,type) (complex ,type)) *) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/bb5c9e106d89594a1bc3020… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/bb5c9e106d89594a1bc3020… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-456-more-accurate-complex-div] 3 commits: Move cdiv functions from float-tran to numbers
by Raymond Toy (＠rtoy) 10 Dec '25

10 Dec '25

Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 270c3daf by Raymond Toy at 2025-12-10T10:14:57-08:00 Move cdiv functions from float-tran to numbers Move the cdiv functions from float-tran to numbers because the interpreter needs to use these when dispatching /. The deftransforms are fine with using the version from the KERNEL package. The tests now pass except for one where the precision is 52 bits instead of 53. More cleanup needed. - - - - - e1e1fb47 by Raymond Toy at 2025-12-10T14:58:33-08:00 Clean up implementation and fix contagion issues First, just define `cdiv-double-float` that implements Baudin and Smith for (complex double-float). For (complex single-float), convert to double-float and use Smith's algorithm instead of the more complex Baudin and Smith. Fix up the numeric contagion so we call these two routines appropriately. (This is quite messy!) Update the deftransforms to call kernel::cdiv-double-float and kernel::cdiv-single-float. - - - - - 825f7d96 by Raymond Toy at 2025-12-10T15:19:06-08:00 Fix bug in cdiv-single-float and add/update tests cdiv-single-float needs to convert single-float to double-float so we don't overflow. Update the accuracy of the one test to be 52 instead of 53. The equivalent C code has the same accuracy. Add a test that we don't overflow for division for complex single-float number. Finally, add tests to make sure we got the all the numeric contagion cases covered for division. Quite useful for getting the contagion cases covered. - - - - - 3 changed files: - src/code/numbers.lisp - src/compiler/float-tran.lisp - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -594,6 +594,7 @@ (build-ratio (truncate x gcd) (truncate y gcd)))))) +#+nil (defun two-arg-/ (x y) (number-dispatch ((x number) (y number)) (float-contagion / x y (ratio integer)) @@ -656,6 +657,441 @@ (build-ratio (maybe-truncate nx gcd) (* (maybe-truncate y gcd) (denominator x))))))) +;; An implementation of Baudin and Smith's robust complex division. +;; This is a pretty straightforward translation of the original in +;; https://arxiv.org/pdf/1210.4539. +#+nil +(eval-when (:compile-toplevel :load-toplevel :execute) +(macrolet + ((frob (type max-float min-float eps-float) + (let ((name (symbolicate "CDIV-" type))) + `(progn + (let ((ov ,max-float) + (un ,min-float) + (eps ,eps-float)) + (defun ,name (x y) + (declare (type (complex ,type) x y)) + (labels + ((internal-compreal (a b c d r tt) + (declare (double-float a b c d r tt)) + ;; Iteration 1: compare r against DBL_MIN instead of 0 + (cond ((/= r 0) + (let ((br (* b r))) + (if (/= br 0) + (* (+ a br) tt) + (+ (* a tt) + (* (* b tt) + r))))) + (t + (* (+ a (* d (/ b c))) + tt)))) + + (robust-subinternal (a b c d) + (declare (double-float a b c d)) + (let* ((r (/ d c)) + (tt (/ (+ c (* d r))))) + (let ((e (internal-compreal a b c d r tt)) + (f (internal-compreal b (- a) c d r tt))) + #+nil(format t "subint: e f = ~A ~A~%" e f) + (values e + f)))) + + (robust-internal (x y) + (declare (type (complex ,type) x y)) + #+nil(format t "x = ~A, y = ~A~%" x y) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (cond + ((<= (abs d) (abs c)) + (multiple-value-bind (e f) + (robust-subinternal a b c d) + #+nil(format t "1: e f = ~A ~A~%" e f) + (complex e f))) + (t + (multiple-value-bind (e f) + (robust-subinternal b a d c) + #+nil(format t "2: e f = ~A ~A~%" e f) + (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))) + (b 2d0) + (s 1d0) + (be (/ b (* eps eps)))) + (when (>= ab (/ ov 2)) + (setf x (/ x 2) + s (* s 2))) + (when (>= cd (/ ov 2)) + (setf y (/ y 2) + s (/ s 2))) + (when (<= ab (* un (/ b eps))) + (setf x (* x be) + s (/ s be))) + (when (<= cd (* un (/ b eps))) + (setf y (* y be) + s (* s be))) + ;; Iteration 2 + #+nil + (when (< cd eps) + (setf x (/ x eps)) + (setf y (/ y eps))) + + ;; Iteration 4 + #+nil + (when (> cd (/ ov 2)) + (setf x (/ x 2) + y (/ y 2))) + + (* s + (robust-internal x y)))))))))) + ;; The eps value for double-float is determined by looking at the + ;; value of %eps in Scilab. + (frob double-float most-positive-double-float least-positive-normalized-double-float + (scale-float 1d0 -52)) + ;; The eps value here is a guess. + (frob single-float most-positive-single-float least-positive-normalized-single-float + (scale-float 1f0 -22))) +) + +(let ((ov most-positive-double-float) + (un least-positive-normalized-double-float) + (eps (scale-float 1d0 -52))) + (defun cdiv-double-float (x y) + (declare (type (complex double-float) x y)) + (labels + ((internal-compreal (a b c d r tt) + (declare (double-float a b c d r tt)) + ;; Iteration 1: compare r against DBL_MIN instead of 0 + (cond ((/= r 0) + (let ((br (* b r))) + (if (/= br 0) + (* (+ a br) tt) + (+ (* a tt) + (* (* b tt) + r))))) + (t + (* (+ a (* d (/ b c))) + tt)))) + + (robust-subinternal (a b c d) + (declare (double-float a b c d)) + (let* ((r (/ d c)) + (tt (/ (+ c (* d r))))) + (let ((e (internal-compreal a b c d r tt)) + (f (internal-compreal b (- a) c d r tt))) + #+nil(format t "subint: e f = ~A ~A~%" e f) + (values e + f)))) + + (robust-internal (x y) + (declare (type (complex double-float) x y)) + #+nil(format t "x = ~A, y = ~A~%" x y) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (cond + ((<= (abs d) (abs c)) + (multiple-value-bind (e f) + (robust-subinternal a b c d) + #+nil(format t "1: e f = ~A ~A~%" e f) + (complex e f))) + (t + (multiple-value-bind (e f) + (robust-subinternal b a d c) + #+nil(format t "2: e f = ~A ~A~%" e f) + (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))) + (b 2d0) + (s 1d0) + (be (/ b (* eps eps)))) + (when (>= ab (/ ov 2)) + (setf x (/ x 2) + s (* s 2))) + (when (>= cd (/ ov 2)) + (setf y (/ y 2) + s (/ s 2))) + (when (<= ab (* un (/ b eps))) + (setf x (* x be) + s (/ s be))) + (when (<= cd (* un (/ b eps))) + (setf y (* y be) + s (* s be))) + ;; Iteration 2 + #+nil + (when (< cd eps) + (setf x (/ x eps)) + (setf y (/ y eps))) + + ;; Iteration 4 + #+nil + (when (> cd (/ ov 2)) + (setf x (/ x 2) + y (/ y 2))) + + (* s + (robust-internal x y)))))) + +(defun cdiv-single-float (x y) + (declare (type (complex single-float) x y)) + (let ((a (float (realpart x) 1d0)) + (b (float (imagpart x) 1d0)) + (c (float (realpart y) 1d0)) + (d (float (imagpart y) 1d0))) + (cond ((< (abs c) (abs d)) + (let* ((r (/ c d)) + (denom (+ (* c r) d)) + (e (float (/ (+ (* a r) b) denom) 1f0)) + (f (float (/ (- (* b r) a) denom) 1f0))) + (complex e f))) + (t + (let* ((r (/ d c)) + (denom (+ c (* d r))) + (e (float (/ (+ a (* b r)) denom) 1f0)) + (f (float (/ (- b (* a r)) denom) 1f0))) + (complex e f)))))) + +(defun two-arg-/ (x y) + (number-dispatch ((x number) (y number)) + (float-contagion / x y (ratio integer)) + + (((complex single-float) + (foreach (complex rational) (complex single-float))) + (cdiv-single-float x (coerce y '(complex single-float)))) + (((complex double-float) + (foreach (complex rational) (complex single-float) (complex double-float))) + (cdiv-double-float x (coerce y '(complex double-float)))) + + (((foreach integer ratio single-float (complex rational)) + (complex single-float)) + (cdiv-single-float (coerce x '(complex single-float)) + y)) + + (((foreach integer ratio single-float double-float (complex rational) (complex single-float)) + (complex double-float)) + (cdiv-double-float (coerce x '(complex double-float)) + y)) + (((complex double-double-float) + (foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float))) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (cond ((< (abs c) (abs d)) + (let* ((r (/ c d)) + (denom (+ (* c r) d)) + (e (/ (+ (* a r) b) denom)) + (f (/ (- (* b r) a) denom))) + (canonical-complex e f))) + (t + (let* ((r (/ d c)) + (denom (+ c (* d r))) + (e (/ (+ a (* b r)) denom)) + (f (/ (- b (* a r)) denom))) + (canonical-complex e f)))))) + + (((foreach integer ratio single-float double-float double-double-float + (complex rational) (complex single-float) (complex double-float)) + (complex double-double-float)) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (cond ((< (abs c) (abs d)) + (let* ((r (/ c d)) + (denom (+ (* c r) d)) + (e (/ (+ (* a r) b) denom)) + (f (/ (- (* b r) a) denom))) + (canonical-complex e f))) + (t + (let* ((r (/ d c)) + (denom (+ c (* d r))) + (e (/ (+ a (* b r)) denom)) + (f (/ (- b (* a r)) denom))) + (canonical-complex e f)))))) + + (((foreach integer ratio single-float double-float double-double-float) + (complex rational)) + (let* ((ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (* ry (+ 1 (* r r))))) + (canonical-complex (/ x dn) + (/ (- (* x r)) dn))) + (let* ((r (/ ry iy)) + (dn (* iy (+ 1 (* r r))))) + (canonical-complex (/ (* x r) dn) + (/ (- x) dn)))))) + (((complex rational) + (complex rational)) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (cond ((< (abs c) (abs d)) + (let* ((r (/ c d)) + (denom (+ (* c r) d)) + (e (/ (+ (* a r) b) denom)) + (f (/ (- (* b r) a) denom))) + (canonical-complex e f))) + (t + (let* ((r (/ d c)) + (denom (+ c (* d r))) + (e (/ (+ a (* b r)) denom)) + (f (/ (- b (* a r)) denom))) + (canonical-complex e f)))))) + + (((foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float)) + (or rational float)) + (canonical-complex (/ (realpart x) y) + (/ (imagpart x) y))) + + ((double-float + (complex single-float)) + (cdiv-double-float (coerce x '(complex double-float)) + (coerce y '(complex double-float)))) + + #+nil + (((foreach integer ratio single-float double-float) complex) + (let* ((ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (* ry (+ 1 (* r r))))) + (canonical-complex (/ x dn) + (/ (- (* x r)) dn))) + (let* ((r (/ ry iy)) + (dn (* iy (+ 1 (* r r))))) + (canonical-complex (/ (* x r) dn) + (/ (- x) dn)))))) + #+nil + ((complex (or rational float)) + (canonical-complex (/ (realpart x) y) + (/ (imagpart x) y))) + + ((ratio ratio) + (let* ((nx (numerator x)) + (dx (denominator x)) + (ny (numerator y)) + (dy (denominator y)) + (g1 (gcd nx ny)) + (g2 (gcd dx dy))) + (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2)) + (* (maybe-truncate dx g2) (maybe-truncate ny g1))))) + + ((integer integer) + (integer-/-integer x y)) + + ((integer ratio) + (if (zerop x) + 0 + (let* ((ny (numerator y)) + (dy (denominator y)) + (gcd (gcd x ny))) + (build-ratio (* (maybe-truncate x gcd) dy) + (maybe-truncate ny gcd))))) + + ((ratio integer) + (let* ((nx (numerator x)) + (gcd (gcd nx y))) + (build-ratio (maybe-truncate nx gcd) + (* (maybe-truncate y gcd) (denominator x))))))) + + +#+nil +(defun two-arg-/ (x y) + (number-dispatch ((x number) (y number)) + (float-contagion / x y (ratio integer)) + + (((complex double-float) (complex double-float)) + (cdiv-double-float x y)) + (((complex double-float) (complex single-float)) + (cdiv-double-float x (coerce y '(complex double-float)))) + + (((complex single-float) (complex double-float)) + (cdiv-double-float (coerce y '(complex double-float)) + y)) + (((complex single-float) (complex single-float)) + (cdiv-single-float x y)) + + #+nil + (((foreach (complex single-float) (complex double-float)) + (complex double-float)) + (cdiv-double-float (complex (float (realpart x) 1d0) + (float (imagpart x) 1d0)) + y)) + (((foreach (complex rational) (complex double-double-float)) + (foreach (complex rational) (complex double-double-float))) + (let* ((rx (realpart x)) + (ix (imagpart x)) + (ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (+ ry (* r iy)))) + (canonical-complex (/ (+ rx (* ix r)) dn) + (/ (- ix (* rx r)) dn))) + (let* ((r (/ ry iy)) + (dn (+ iy (* r ry)))) + (canonical-complex (/ (+ (* rx r) ix) dn) + (/ (- (* ix r) rx) dn)))))) + (((foreach integer ratio single-float double-float) + (foreach (complex rational) (complex double-double-float))) + (let* ((ry (realpart y)) + (iy (imagpart y))) + (if (> (abs ry) (abs iy)) + (let* ((r (/ iy ry)) + (dn (* ry (+ 1 (* r r))))) + (canonical-complex (/ x dn) + (/ (- (* x r)) dn))) + (let* ((r (/ ry iy)) + (dn (* iy (+ 1 (* r r))))) + (canonical-complex (/ (* x r) dn) + (/ (- x) dn)))))) + (((foreach (complex rational) (complex double-double-float)) + (or rational float)) + (canonical-complex (/ (realpart x) y) + (/ (imagpart x) y))) + + ((ratio ratio) + (let* ((nx (numerator x)) + (dx (denominator x)) + (ny (numerator y)) + (dy (denominator y)) + (g1 (gcd nx ny)) + (g2 (gcd dx dy))) + (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2)) + (* (maybe-truncate dx g2) (maybe-truncate ny g1))))) + + ((integer integer) + (integer-/-integer x y)) + + ((integer ratio) + (if (zerop x) + 0 + (let* ((ny (numerator y)) + (dy (denominator y)) + (gcd (gcd x ny))) + (build-ratio (* (maybe-truncate x gcd) dy) + (maybe-truncate ny gcd))))) + + ((ratio integer) + (let* ((nx (numerator x)) + (gcd (gcd nx y))) + (build-ratio (maybe-truncate nx gcd) + (* (maybe-truncate y gcd) (denominator x))))))) + (defun %negate (n) (number-dispatch ((n number)) ===================================== src/compiler/float-tran.lisp ===================================== @@ -1917,104 +1917,13 @@ (frob (complex single-float) 1f0) (frob (complex double-float) 1d0)) -;; An implementation of Baudin and Smith's robust complex division. -;; This is a pretty straightforward translation of the original in -;; https://arxiv.org/pdf/1210.4539. (macrolet - ((frob (type max-float min-float eps-float) - (let ((name (symbolicate "CDIV-" type))) - `(progn - (let ((ov ,max-float) - (un ,min-float) - (eps ,eps-float)) - (defun ,name (x y) - (declare (type (complex ,type) x y)) - (labels - ((internal-compreal (a b c d r tt) - (declare (,type a b c d r tt)) - (cond ((zerop r) - (let ((br (* b r))) - (if (/= br 0) - (* (+ a br) tt) - (+ (* a tt) - (* (* b tt) - r))))) - (t - (* (+ a (* d (/ b c))) - tt)))) - - (robust-subinternal (a b c d) - (declare (,type a b c d)) - (let* ((r (/ d c)) - (tt (/ (+ c (* d r))))) - (values (internal-compreal a b c d r tt) - (internal-compreal b (- a) c d r tt)))) - - (robust-internal (x y) - (declare (type (complex ,type) x y)) - (let ((a (realpart x)) - (b (imagpart x)) - (c (realpart y)) - (d (imagpart y))) - (cond - ((<= (abs d) (abs c)) - (multiple-value-bind (e f) - (robust-subinternal a b c d) - (complex e f))) - (t - (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))) - (b 2d0) - (s 1d0) - (be (/ b (* eps eps)))) - (when (>= ab (/ ov 2)) - (setf x (/ x 2) - s (* s 2))) - (when (>= cd (/ ov 2)) - (setf y (/ y 2) - s (/ s 2))) - (when (<= ab (* un (/ b eps))) - (setf x (* x be) - s (/ s be))) - (when (<= cd (* un (/ b eps))) - (setf y (* y be) - s (* s be))) - ;; Iteration 2 - #+nil - (when (< cd eps) - (setf x (/ x eps)) - (setf y (/ y eps))) - - ;; Iteration 4 - #+nil - (when (> cd (/ ov 2)) - (setf x (/ x 2) - y (/ y 2))) - - (* s - (robust-internal x y)))))))))) - ;; The eps value for double-float is determined by looking at the - ;; value of %eps in Scilab. - (frob double-float most-positive-double-float least-positive-normalized-double-float - (scale-float 1d0 -52)) - ;; The eps value here is a guess. - (frob single-float most-positive-single-float least-positive-normalized-single-float - (scale-float 1f0 -22))) - -(macrolet - ((frob (type) - (let ((name (symbolicate "CDIV-" type))) + ((frob (type name) `(deftransform / ((x y) ((complex ,type) (complex ,type)) *) - (,name x y))))) - (frob double-float) - (frob single-float)) + (,name x y)))) + (frob double-float kernel::cdiv-double-float) + (frob single-float kernel::cdiv-single-float)) + ;;;; Complex contagion: ===================================== tests/float.lisp ===================================== @@ -388,7 +388,7 @@ (list (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) (complex 0.6d0 0.2d0) - 53 least-positive-double-float) + 52 least-positive-double-float) ;; 9 (list (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) @@ -456,7 +456,7 @@ (min (real-ulp (realpart diff)) (real-ulp (imagpart diff)))))) -(define-test complex-division +(define-test complex-division.double (:tag :issues) (loop for k from 1 for test in *test-cases* @@ -473,3 +473,35 @@ k x y z z-true diff rel) (assert-true (<= ulp max-ulp) k x y z z-true diff ulp max-ulp))))) + +(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/e11ed3eb2b57116cf92220… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/e11ed3eb2b57116cf92220… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl] Pushed new branch issue-456-more-accurate-complex-div
by Raymond Toy (＠rtoy) 10 Dec '25

10 Dec '25

Raymond Toy pushed new branch issue-456-more-accurate-complex-div at cmucl / cmucl -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/tree/issue-456-more-accurate-c… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][master] 2 commits: Fix #455: Allow manual pipeline runs
by Raymond Toy (＠rtoy) 05 Dec '25

05 Dec '25

Raymond Toy pushed to branch master at cmucl / cmucl Commits: 25cd44a2 by Raymond Toy at 2025-12-04T17:57:49-08:00 Fix #455: Allow manual pipeline runs - - - - - ae1bdde7 by Raymond Toy at 2025-12-04T17:57:50-08:00 Merge branch 'issue-455-manual-pipeline-run' into 'master' Fix #455: Allow manual pipeline runs Closes #455 See merge request cmucl/cmucl!334 - - - - - 1 changed file: - .gitlab-ci.yml Changes: ===================================== .gitlab-ci.yml ===================================== @@ -133,6 +133,7 @@ linux:install: - if: $CI_PIPELINE_SOURCE == "push" - if: $CI_PIPELINE_SOURCE == "merge_request_event" - if: $CI_PIPELINE_SOURCE == "branch" + - if: $CI_PIPELINE_SOURCE == "web" linux:build: <<: *build_configuration @@ -148,6 +149,7 @@ linux:build: - if: $CI_PIPELINE_SOURCE == "schedule" - if: $CI_PIPELINE_SOURCE == "push" - if: $CI_PIPELINE_SOURCE == "merge_request_event" + - if: $CI_PIPELINE_SOURCE == "web" linux:cross-build: stage: build View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/16c334b0116235288a7ea7… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/16c334b0116235288a7ea7… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0