Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: 23c3c67a by Raymond Toy at 2026-01-11T17:11:18-08:00 Remove old macrolet and don't make internal routines available Don't need the old macrolet anymore so remove it. The internal routines aren't used anywhere else so don't expose them outside the block compilation. - - - - - 3b75ab81 by Raymond Toy at 2026-01-14T08:51:03-08:00 Use double-floats for constants for dd cdiv; update comments All of the constants used for the double-double float version are powers of two and can be exactly represented as doubles instead of double-doubles. This allows slightly faster code we can use simpler routines between a double-double and a double. Add comment about this and also fix up some of the comments that had typos or were incorrect. - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -612,25 +612,33 @@ (defconstant +cdiv-be+ (/ 2 (* +cdiv-eps+ +cdiv-eps+))) (defconstant +cdiv-2/eps+ (/ 2 +cdiv-eps+)) -;; Same constants but for double-double-floats. Some of these aren't -;; well-defined for double-double-floats so we make our best guess at +;; Same constants but for DOUBLE-DOUBLE-FLOATS. Some of these aren't +;; well-defined for DOUBLE-DOUBLE-FLOATS so we make our best guess at ;; what they might be. Since double-doubles have about twice as many -;; bits of precision as a double-float, we generally just double the -;; exponent of the corresponding double-float values above. +;; bits of precision as a DOUBLE-FLOAT, we generally just double the +;; exponent (for SCALE-FLOAT) of the corresponding DOUBLE-FLOAT values +;; above. +;; +;; Note also that both LEAST-POSITIVE-NORMALIZED-DOUBLE-DOUBLE-FLOAT +;; and MOST-POSITIVE-DOUBLE-DOUBLE-FLOAT have a low component of 0, so +;; there's no loss of precision if we use the corresponding +;; DOUBLE-FLOAT values. Likewise, all the other constants here are +;; powers of 2, so the DOUBLE-DOUBLE-FLOAT values can be represented +;; exactly as a DOUBLE-FLOAT. We can use DOUBLE-FLOAT values. (defconstant +cdiv-dd-rmin+ - least-positive-normalized-double-double-float) + least-positive-normalized-double-float) (defconstant +cdiv-dd-rbig+ - (/ most-positive-double-double-float 2)) + (/ most-positive-double-float 2)) (defconstant +cdiv-dd-rmin2+ - (scale-float 1w0 -106)) + (scale-float 1d0 -106)) (defconstant +cdiv-dd-rminscal+ - (scale-float 1w0 102)) + (scale-float 1d0 102)) (defconstant +cdiv-dd-rmax2+ (* +cdiv-dd-rbig+ +cdiv-dd-rmin2+)) ;; Epsilon for double-doubles isn't really well-defined because things ;; like (+ 1w0 1w-200) is a valid double-double float. (defconstant +cdiv-dd-eps+ - (scale-float 1w0 -104)) + (scale-float 1d0 -104)) (defconstant +cdiv-dd-be+ (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+))) (defconstant +cdiv-dd-2/eps+ @@ -643,179 +651,11 @@ ;; else. (declaim (ext:start-block cdiv-double-float cdiv-single-float cdiv-double-double-float - two-arg-/ - cdiv-double-float-internal-compreal - cdiv-double-float-robust-subinternal - cdiv-double-float-robust-internal - cdiv-double-double-float-internal-compreal - cdiv-double-double-float-robust-subinternal - cdiv-double-double-float-robust-internal)) + two-arg-/)) ;; 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. -#+nil -(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 0.5d0) - s (* s 2))) - ;; If c or d is big, scale down c and d. - (when (>= cd ,+rbig+) - (setf y (* y 0.5d0) - s (* s 0.5d0))) - ;; 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)) - (eval-when (:compile-toplevel) (defmacro define-cdiv (type) (let* ((dd (if (eq type 'double-float) @@ -842,11 +682,14 @@ (defun ,compreal (a b c d r tt) (declare (,type a b c d r tt) ,@opt) - ;; Compute the real part of the complex division - ;; (a+ib)/(c+id), assuming |c| <= |d|. r = d/c and tt = 1/(c+d*r). + ;; Computes the real part of the complex division + ;; (a+ib)/(c+id), assuming |d| <= |c|. Then let r = d/c and + ;; tt = (c+d*r). ;; ;; The realpart is (a*c+b*d)/(c^2+d^2). ;; + ;; The denominator is + ;; ;; c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r) ;; ;; Then @@ -862,15 +705,12 @@ (/ (+ 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. + ;; r is very small so d is very tiny compared to c. ;; (/ (+ a (* d (/ b c))) tt)))) @@ -931,9 +771,7 @@ ;; perform the division. (multiple-value-bind (a b c d) (maybe-scale (abs c) a b c d) - (multiple-value-bind (e f) - (,subinternal a b c d) - (complex e f)))) + (,subinternal a b c d))) (t ;; |d| > |c|. So, instead compute ;; @@ -947,7 +785,7 @@ (maybe-scale (abs d) a b c d) (multiple-value-bind (e f) (,subinternal b a d c) - (complex e (- f))))))))) + (values e (- f))))))))) (defun ,name (x y) ,docstring (declare (type (complex ,type) x y) @@ -982,14 +820,15 @@ (when (<= cd (* ,+rmin+ ,+2/eps+)) (setf y (* y ,+be+) s (* s ,+be+))) - (* s - (,internal x y))))))) + (multiple-value-bind (e f) + (,internal x y) + (complex (* s e) + (* s f)))))))) ) (define-cdiv double-float) (define-cdiv double-double-float) - ;; 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) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/e1a0c3802ba8a4079d384e6... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/e1a0c3802ba8a4079d384e6... You're receiving this email because of your account on gitlab.common-lisp.net.