Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: b29c585b by Raymond Toy at 2026-01-19T11:07:59-08:00 Let bind the constants as needed in each routine This also means the macro needs args to define the values. Address review comments. - - - - - 184720ee by Raymond Toy at 2026-01-19T11:18:04-08:00 Remove unneeded defconstants and update/clarify some comments - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -602,53 +602,11 @@ ;; In particular iteration 1 and 3 are added. Iteration 2 and 4 were ;; not added. The test examples from iteration 2 and 4 didn't change ;; with or without changes added. -(defconstant +cdiv-rmin+ least-positive-normalized-double-float) -(defconstant +cdiv-rbig+ (/ most-positive-double-float 2)) -(defconstant +cdiv-rmin2+ (scale-float 1d0 -53)) -(defconstant +cdiv-rminscal+ (scale-float 1d0 51)) -(defconstant +cdiv-rmax2+ (* +cdiv-rbig+ +cdiv-rmin2+)) -;; This is the value of %eps from Scilab -(defconstant +cdiv-eps+ (scale-float 1d0 -52)) -(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 -;; 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 (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-float) -(defconstant +cdiv-dd-rbig+ - (/ most-positive-double-float 2)) -(defconstant +cdiv-dd-rmin2+ - (scale-float 1d0 -106)) -(defconstant +cdiv-dd-rminscal+ - (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 1d0 -104)) -(defconstant +cdiv-dd-be+ - (/ 2 (* +cdiv-dd-eps+ +cdiv-dd-eps+))) -(defconstant +cdiv-dd-2/eps+ - (/ 2 +cdiv-dd-eps+)) - - -;; Make these functions accessible. cdiv-double-float and -;; cdiv-single-float are used by deftransforms. Of course, two-arg-/ -;; is the interface to division. cdiv-generic isn't used anywhere -;; else. +;; Make these functions accessible outside the block compilation. +;; CDIV-DOUBLE-FLOAT and CDIV-SINGLE-FLOAT are used by deftransforms. +;; Of course, TWO-ARG-/ is the interface to division. CDIV-GENERIC +;; isn't used anywhere else. (declaim (ext:start-block cdiv-double-float cdiv-single-float cdiv-double-double-float two-arg-/)) @@ -657,27 +615,15 @@ ;; implementation is basically identical except for the constants. ;; use a macro to generate both versions. (eval-when (:compile-toplevel) -(defmacro define-cdiv (type) - (let* ((dd (if (eq type 'double-float) - "" - "DD-")) - (opt '((optimize (speed 3) (safety 0)))) +(defmacro define-cdiv (type &key rmin rbig rmin2 rminscal eps) + (let* ((opt '((optimize (speed 3) (safety 0)))) (name (symbolicate "CDIV-" type)) (compreal (symbolicate "CDIV-" type "-INTERNAL-COMPREAL")) (subinternal (symbolicate "CDIV-" type "-ROBUST-SUBINTERNAL")) (internal (symbolicate "CDIV-" type "-ROBUST-INTERNAL")) (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+")) - (+be+ (symbolicate "+CDIV-" dd "BE+")) - (+2/eps+ (symbolicate "+CDIV-" dd "2/EPS+"))) + type)))) `(progn (defun ,compreal (a b c d r tt) (declare (,type a b c d r tt) @@ -699,7 +645,7 @@ ;; = (a + b*r)/(c + d*r). ;; ;; Thus tt = (c + d*r). - (cond ((>= (abs r) ,+rmin+) + (cond ((>= (abs r) ,rmin) (let ((br (* b r))) (if (/= br 0) (/ (+ a br) tt) @@ -719,115 +665,154 @@ ,@opt) (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). + ;; e is the real part and f is the imaginary part. We can + ;; use 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 (,compreal a b c d r tt)) (f (,compreal b (- a) c d r tt))) (values e f)))) (defun ,internal (x y) - (declare (type (complex ,type) x y) - ,@opt) - (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) - (,subinternal a b c d))) - (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) - (,subinternal b a d c) - (values e (- f))))))))) - (defun ,name (x y) - ,docstring - (declare (type (complex ,type) x y) - ,@opt) - (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 is always multipled by a power of 2. It's either - ;; 2, 0.5, or +be+, which is also a power of two. - ;; (Either 2^105 or 2^209). Hence, we can just use a - ;; double-float instead of a double-double-float - ;; because the double-double always has 0d0 for the lo - ;; part. - (s 1d0)) - (declare (double-float 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+))) - (multiple-value-bind (e f) - (,internal x y) - (complex (* s e) - (* s f)))))))) + (declare (type (complex ,type) x y) + ,@opt) + (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 ((+rbig+ ,rbig) + (+rmin2+ ,rmin2) + (+rminscal+ ,rminscal) + (+rmax2+ (* ,rbig ,rmin2)) + (+rmin+ ,rmin) + (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) + (,subinternal a b c d))) + (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) + (,subinternal b a d c) + (values e (- f))))))))) + (defun ,name (x y) + ,docstring + (declare (type (complex ,type) x y) + ,@opt) + (let ((+rmin+ ,rmin) + (+rbig+ ,rbig) + (+2/eps+ (/ 2 ,eps)) + (+be+ (/ 2 (* ,eps ,eps)))) + (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 is always multipled by a power of 2. It's either + ;; 2, 0.5, or +be+, which is also a power of two. + ;; (Either 2^105 or 2^209). Hence, we can just use a + ;; double-float instead of a double-double-float + ;; because the double-double always has 0d0 for the lo + ;; part. + (s 1d0)) + (declare (double-float 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+))) + (multiple-value-bind (e f) + (,internal x y) + (complex (* s e) + (* s f))))))))) ) -(define-cdiv double-float) -(define-cdiv double-double-float) +(define-cdiv double-float + :rmin least-positive-normalized-double-float + :rbig (/ most-positive-double-float 2) + :rmin2 (scale-float 1d0 -53) + :rminscal (scale-float 1d0 51) + ;; This is the value of %eps from Scilab + :eps (scale-float 1d0 -52)) + +;; Same constants but for DOUBLE-DOUBLE-FLOATS. Some of these aren't +;; well-defined for DOUBLE-DOUBLE-FLOATS (like eps) so we make our +;; best guess at what they might be. Since double-doubles have +;; twice as many 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 since both +;; LEAST-POSITIVE-NORMALIZED-DOUBLE-DOUBLE-FLOAT and +;; MOST-POSITIVE-DOUBLE-DOUBLE-FLOAT have a low component of 0, +;; 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. This +;; also makes some of the operations a bit simpler since operations +;; between a DOUBLE-DOUBLE-FLOAT and a DOUBLE-FLOAT are simpler than +;; if both were DOUBLE-DOUBLE-FLOATs. +(define-cdiv double-double-float + :rmin least-positive-normalized-double-float + :rbig (/ most-positive-double-float 2) + :rmin2 (scale-float 1d0 -106) + :rminscal (scale-float 1d0 102) + :eps (scale-float 1d0 -104)) + ;; Smith's algorithm for complex division for (complex single-float). ;; We convert the parts to double-floats before computing the result. View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/7ef6184520ac5201dd878e9... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/7ef6184520ac5201dd878e9... You're receiving this email because of your account on gitlab.common-lisp.net.