Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl Commits: 262c2323 by Raymond Toy at 2026-01-07T11:02:30-08:00 Factor out common code to cdiv-generic cdiv-generic implements Smith's algorithm for other cases with mixing different types of reals and complexes together. But we leave the case of a real divided by a complex as is, because it has a simpler implementation. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -774,6 +774,25 @@ (f (float (/ (- b (* a r)) denom) 1f0))) (complex e f)))))) +;; Generic implementation of Smith's algorithm. +(defun cdiv-generic (x y) + (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)))))) + (defun two-arg-/ (x y) (number-dispatch ((x number) (y number)) (float-contagion / x y (ratio integer)) @@ -799,45 +818,17 @@ (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)) - (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)))))) + (cdiv-generic x y)) (((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)))))) + (cdiv-generic x y)) (((foreach integer ratio single-float double-float double-double-float) (complex rational)) + ;; Smith's algorithm, but takes advantage of the fact that the + ;; numerator is a real number and not complex. (let* ((ry (realpart y)) (iy (imagpart y))) (if (> (abs ry) (abs iy)) @@ -853,22 +844,7 @@ (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)) - (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)))))) + (cdiv-generic x y)) (((foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float)) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/262c2323e4d18089405d76f5... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/262c2323e4d18089405d76f5... You're receiving this email because of your account on gitlab.common-lisp.net.