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/bb5c9e106d89594a1bc3020... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/bb5c9e106d89594a1bc3020... You're receiving this email because of your account on gitlab.common-lisp.net.