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