Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: a26a83cd by Raymond Toy at 2026-01-08T20:17:09-08:00 Macroize cdiv for double-float and double-double-float Since the code for double-float and double-double-float is basically identical except for the constants used, create a macro to generate both versions with the appropriate constants. - - - - - 32ca6718 by Raymond Toy at 2026-01-08T20:20:54-08:00 Remove old impl that has been replaced by the macro. - - - - - 74c910e5 by Raymond Toy at 2026-01-08T20:21:11-08:00 Update due to changed docstrings. - - - - - 2 changed files: - src/code/numbers.lisp - src/i18n/locale/cmucl.pot Changes: ===================================== src/code/numbers.lisp ===================================== @@ -645,144 +645,169 @@ cdiv-double-double-float two-arg-/)) -(defun cdiv-double-float (x y) - "Accurate division of complex double-float numbers x and y." - (declare (type (complex double-float) x y) - (optimize (speed 3) (safety 0))) - (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). - (cond ((>= (abs r) +cdiv-rmin+) - (let ((br (* b r))) - (if (/= br 0) - (/ (+ 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 - ;; r = 0 so d is very tiny compared to c. - ;; - (/ (+ 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)))) - ;; 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))) - (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))) - (declare (double-float a b c d)) - (flet ((maybe-scale (abs-tst a b c d) - (declare (double-float a b c d)) - ;; 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 +cdiv-rbig+) - (> abs-a +cdiv-rbig+) - (> abs-b +cdiv-rbig+)) - (setf a (* a 0.5d0) - b (* b 0.5d0) - c (* c 0.5d0) - d (* d 0.5d0)) - (if (< abs-tst +cdiv-rmin2+) - (setf a (* a +cdiv-rminscal+) - b (* b +cdiv-rminscal+) - c (* c +cdiv-rminscal+) - d (* d +cdiv-rminscal+)) - (if (or (and (< abs-a +cdiv-rmin+) - (< abs-b +cdiv-rmax2+) - (< abs-tst +cdiv-rmax2+)) - (and (< abs-b +cdiv-rmin+) - (< abs-a +cdiv-rmax2+) - (< abs-tst +cdiv-rmax2+))) - (setf a (* a +cdiv-rminscal+) - b (* b +cdiv-rminscal+) - c (* c +cdiv-rminscal+) - d (* d +cdiv-rminscal+))))) - (values a b c d)))) - (cond - ((<= (abs d) (abs c)) - ;; |d| <= |c|, so we can use robust-subinternal to - ;; perform the division. - (multiple-value-bind (a b c d) - (maybe-scale (abs c) a b c d) - (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 (a b c d) - (maybe-scale (abs d) a b c 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))) - (s 1d0)) - (declare (double-float s)) - ;; If a or b is big, scale down a and b. - (when (>= ab +cdiv-rbig+) - (setf x (/ x 2) - s (* s 2))) - ;; If c or d is big, scale down c and d. - (when (>= cd +cdiv-rbig+) - (setf y (/ y 2) - s (/ s 2))) - ;; If a or b is tiny, scale up a and b. - (when (<= ab (* +cdiv-rmin+ +cdiv-2/eps+)) - (setf x (* x +cdiv-be+) - s (/ s +cdiv-be+))) - ;; If c or d is tiny, scale up c and d. - (when (<= cd (* +cdiv-rmin+ +cdiv-2/eps+)) - (setf y (* y +cdiv-be+) - s (* s +cdiv-be+))) - (* s - (robust-internal x y))))) +;; The code for both the double-float and double-double-float +;; implementation is basically identical except for the constants. +;; use a macro to generate both versions. +(macrolet + ((frob (type) + (let* ((dd (if (eq type 'double-float) + "" + "DD-")) + (name (symbolicate "CDIV-" type)) + (docstring (let ((*print-case* :downcase)) + (format nil "Accurate division of complex ~A numbers x and y." + type))) + ;; Create the correct symbols for the constants we need, + ;; as defined above. + (+rmin+ (symbolicate "+CDIV-" dd "RMIN+")) + (+rbig+ (symbolicate "+CDIV-" dd "RBIG+")) + (+rmin2+ (symbolicate "+CDIV-" dd "RMIN2+")) + (+rminscal+ (symbolicate "+CDIV-" dd "RMINSCAL+")) + (+rmax2+ (symbolicate "+CDIV-" dd "RMAX2+")) + (+eps+ (symbolicate "+CDIV-" dd "EPS+")) + (+be+ (symbolicate "+CDIV-" dd "BE+")) + (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+"))) + `(defun ,name (x y) + ,docstring + (declare (type (complex ,type) x y) + (optimize (speed 3) (safety 0))) + (labels + ((internal-compreal (a b c d r tt) + (declare (,type 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). + (cond ((>= (abs r) ,+rmin+) + (let ((br (* b r))) + (if (/= br 0) + (/ (+ 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 + ;; r = 0 so d is very tiny compared to c. + ;; + (/ (+ 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)))) + ;; 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))) + (values e + f)))) + (robust-internal (x y) + (declare (type (complex ,type) x y)) + (let ((a (realpart x)) + (b (imagpart x)) + (c (realpart y)) + (d (imagpart y))) + (declare (,type a b c d)) + (flet ((maybe-scale (abs-tst a b c d) + (declare (,type a b c d)) + ;; 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+))))) + (values a b c d)))) + (cond + ((<= (abs d) (abs c)) + ;; |d| <= |c|, so we can use robust-subinternal to + ;; perform the division. + (multiple-value-bind (a b c d) + (maybe-scale (abs c) a b c d) + (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 (a b c d) + (maybe-scale (abs d) a b c 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))) + (s (coerce 1 ',type))) + (declare (,type s)) + ;; If a or b is big, scale down a and b. + (when (>= ab ,+rbig+) + (setf x (/ x 2) + s (* s 2))) + ;; If c or d is big, scale down c and d. + (when (>= cd ,+rbig+) + (setf y (/ y 2) + s (/ s 2))) + ;; If a or b is tiny, scale up a and b. + (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 (* ,+rmin+ ,+2/eps+)) + (setf y (* y ,+be+) + s (* s ,+be+))) + (* s + (robust-internal x y)))))))) + (frob double-float) + #+double-double + (frob double-double-float)) ;; Smith's algorithm for complex division for (complex single-float). ;; We convert the parts to double-floats before computing the result. @@ -806,146 +831,6 @@ (f (float (/ (- b (* a r)) denom) 1f0))) (complex e f)))))) -(defun cdiv-double-double-float (x y) - "Accurate division of complex double-double-float numbers x and y." - (declare (type (complex double-double-float) x y) - (optimize (speed 2) (safety 0))) - (labels - ((internal-compreal (a b c d r tt) - (declare (double-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). - (cond ((>= (abs r) +cdiv-dd-rmin+) - (let ((br (* b r))) - (if (/= br 0) - (/ (+ 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 - ;; r = 0 so d is very tiny compared to c. - ;; - (/ (+ a (* d (/ b c))) - tt)))) - (robust-subinternal (a b c d) - (declare (double-double-float a b c d)) - (let* ((r (/ d c)) - (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))) - (values e - f)))) - - (robust-internal (x y) - (declare (type (complex double-double-float) x y)) - (let ((a (realpart x)) - (b (imagpart x)) - (c (realpart y)) - (d (imagpart y))) - (declare (double-double-float a b c d)) - (flet ((maybe-scale (abs-tst a b c d) - (declare (double-double-float a b c d)) - ;; 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 +cdiv-dd-rbig+) - (> abs-a +cdiv-dd-rbig+) - (> abs-b +cdiv-dd-rbig+)) - (setf a (* a 0.5d0) - b (* b 0.5d0) - c (* c 0.5d0) - d (* d 0.5d0)) - (if (< abs-tst +cdiv-dd-rmin2+) - (setf a (* a +cdiv-dd-rminscal+) - b (* b +cdiv-dd-rminscal+) - c (* c +cdiv-dd-rminscal+) - d (* d +cdiv-dd-rminscal+)) - (if (or (and (< abs-a +cdiv-dd-rmin+) - (< abs-b +cdiv-dd-rmax2+) - (< abs-tst +cdiv-dd-rmax2+)) - (and (< abs-b +cdiv-dd-rmin+) - (< abs-a +cdiv-dd-rmax2+) - (< abs-tst +cdiv-dd-rmax2+))) - (setf a (* a +cdiv-dd-rminscal+) - b (* b +cdiv-dd-rminscal+) - c (* c +cdiv-dd-rminscal+) - d (* d +cdiv-dd-rminscal+))))) - (values a b c d)))) - (cond - ((<= (abs d) (abs c)) - ;; |d| <= |c|, so we can use robust-subinternal to - ;; perform the division. - (multiple-value-bind (a b c d) - (maybe-scale (abs c) a b c d) - (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 (a b c d) - (maybe-scale (abs d) a b c 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))) - (s 1w0)) - (declare (double-double-float s)) - ;; If a or b is big, scale down a and b. - (when (>= ab +cdiv-dd-rbig+) - (setf x (/ x 2) - s (* s 2))) - ;; If c or d is big, scale down c and d. - (when (>= cd +cdiv-dd-rbig+) - (setf y (/ y 2) - s (/ s 2))) - ;; If a or b is tiny, scale up a and b. - (when (<= ab (* +cdiv-dd-rmin+ +cdiv-dd-2/eps+)) - (setf x (* x +cdiv-dd-be+) - s (/ s +cdiv-dd-be+))) - ;; If c or d is tiny, scale up c and d. - (when (<= cd (* +cdiv-dd-rmin+ +cdiv-dd-2/eps+)) - (setf y (* y +cdiv-dd-be+) - s (* s +cdiv-dd-be+))) - (* s - (robust-internal x y))))) - ;; Generic implementation of Smith's algorithm. (defun cdiv-generic (x y) "Complex division of generic numbers x and y. One of x or y should be ===================================== src/i18n/locale/cmucl.pot ===================================== @@ -4390,11 +4390,11 @@ msgid "Accurate division of complex double-float numbers x and y." msgstr "" #: src/code/numbers.lisp -msgid "Accurate division of complex single-float numbers x and y." +msgid "Accurate division of complex double-double-float numbers x and y." msgstr "" #: src/code/numbers.lisp -msgid "Accurate division of complex double-double-float numbers x and y." +msgid "Accurate division of complex single-float numbers x and y." msgstr "" #: src/code/numbers.lisp View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/f1d99fbc00e0f789b610cc6... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/f1d99fbc00e0f789b610cc6... You're receiving this email because of your account on gitlab.common-lisp.net.