Raymond Toy pushed to branch issue-461-cdiv-use-4-doubles at cmucl / cmucl Commits: 336b9322 by Raymond Toy at 2026-01-15T09:13:30-08:00 Update cmucl.pot for new docstrings Forgot to update this when accurate complex division was committed. [skip-ci] - - - - - 1b688607 by Raymond Toy at 2026-01-18T15:38:33-08:00 Fix #463: Handle double-double-float and double-float comparisons - - - - - 67dc5c9d by Raymond Toy at 2026-01-18T15:38:34-08:00 Merge branch 'issue-463-double-double-float-comparison' into 'master' Fix #463: Handle double-double-float and double-float comparisons Closes #463 See merge request cmucl/cmucl!342 - - - - - 06df185b by Raymond Toy at 2026-01-19T13:33:37-08:00 Fix #459: More accurate complex division for double-double-float - - - - - 75324b6d by Raymond Toy at 2026-01-19T13:33:37-08:00 Merge branch 'issue-459-more-accurate-dd-complex-div' into 'master' Fix #459: More accurate complex division for double-double-float Closes #459 and #456 See merge request cmucl/cmucl!338 - - - - - edf61a9d by Raymond Toy at 2026-01-19T14:47:22-08:00 Merge branch 'master' into issue-461-cdiv-use-4-doubles Required fixing merge conflicts with master. - - - - - fae66d26 by Raymond Toy at 2026-01-19T14:59:34-08:00 Fix bad merge Forgot to make internal take 4 args. - - - - - 2 changed files: - src/code/numbers.lisp - src/compiler/float-tran-dd.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,46 +665,51 @@ ,@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 (a b c d) - (declare (,type a b c d) - ,@opt) + (declare (,type a b c d) + ,@opt) (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)) + (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+)) + (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+))))) + (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)) @@ -781,48 +732,82 @@ (multiple-value-bind (e f) (,subinternal b a d c) (values e (- f)))))))) - (defun ,name (a b c d) - ,docstring - (declare (,type a b c d) - ,@opt) - (let* ((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 a (* a 0.5d0) - b (* b 0.5d0) - s (* s 2))) - ;; If c or d is big, scale down c and d. - (when (>= cd ,+rbig+) - (setf c (* c 0.5d0) - d (* d 0.5d0) - s (* s 0.5d0))) - ;; If a or b is tiny, scale up a and b. - (when (<= ab (* ,+rmin+ ,+2/eps+)) - (setf a (* a ,+be+) - b (* b ,+be+) - s (/ s ,+be+))) - ;; If c or d is tiny, scale up c and d. - (when (<= cd (* ,+rmin+ ,+2/eps+)) - (setf c (* c ,+be+) - d (* d ,+be+) - s (* s ,+be+))) - (multiple-value-bind (e f) - (,internal a b c d) - (complex (* s e) - (* s f)))))))) + (defun ,name (xr xi yr yi) + ,docstring + (declare (,type xr xi yr yi) + ,@opt) + (let ((+rmin+ ,rmin) + (+rbig+ ,rbig) + (+2/eps+ (/ 2 ,eps)) + (+be+ (/ 2 (* ,eps ,eps)))) + (let* ((max-x (max (abs xr) (abs xi))) + (max-y (max (abs yr) (abs yi))) + ;; S is always multipled by xr power of 2. It's + ;; multiplied by either 2, 0.5, or +be+, which is + ;; also xr power of two. (Either 2^105 or 2^209). + ;; Hence, we can just use xr double-float instead + ;; of xr double-double-float because the + ;; double-double always has 0d0 for the lo part. + (s 1d0)) + (declare (double-float s)) + ;; If xr or xi is big, scale down xr and xi. + (when (>= max-x +rbig+) + (setf xr (* xr 0.5d0) + xi (* xi 0.5d0) + s (* s 2))) + ;; If yr or yi is big, scale down yr and yi. + (when (>= max-y +rbig+) + (setf yr (* yr 0.5d0) + yi (* yi 0.5d0) + s (* s 0.5d0))) + ;; If xr or xi is tiny, scale up xr and xi. + (when (<= max-x (* +rmin+ +2/eps+)) + (setf xr (* xr +be+) + xi (* xi +be+) + s (/ s +be+))) + ;; If yr or yi is tiny, scale up yr and yi. + (when (<= max-y (* +rmin+ +2/eps+)) + (setf yr (* yr +be+) + yi (* yi +be+) + s (* s +be+))) + (multiple-value-bind (e f) + (,internal xr xi yr yi) + (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. ===================================== src/compiler/float-tran-dd.lisp ===================================== @@ -681,6 +681,18 @@ (kernel:double-double-hi b) (kernel:double-double-lo b))) +(deftransform = ((a b) (vm::double-double-float double-float) *) + `(dd= (kernel:double-double-hi a) + (kernel:double-double-lo a) + b + 0d0)) + +(deftransform = ((a b) (double-float vm::double-double-float) *) + `(dd= a + 0d0 + (kernel:double-double-hi b) + (kernel:double-double-lo b))) + (deftransform < ((a b) (vm::double-double-float vm::double-double-float) *) `(dd< (kernel:double-double-hi a) @@ -688,10 +700,34 @@ (kernel:double-double-hi b) (kernel:double-double-lo b))) +(deftransform < ((a b) (vm::double-double-float double-float) *) + `(dd< (kernel:double-double-hi a) + (kernel:double-double-lo a) + b + 0d0)) + +(deftransform < ((a b) (double-float vm::double-double-float) *) + `(dd< a + 0d0 + (kernel:double-double-hi b) + (kernel:double-double-lo b))) + (deftransform > ((a b) (vm::double-double-float vm::double-double-float) *) `(dd> (kernel:double-double-hi a) (kernel:double-double-lo a) (kernel:double-double-hi b) (kernel:double-double-lo b))) + +(deftransform > ((a b) (vm::double-double-float double-float) *) + `(dd> (kernel:double-double-hi a) + (kernel:double-double-lo a) + b + 0d0)) + +(deftransform > ((a b) (double-float vm::double-double-float) *) + `(dd> a + 0d0 + (kernel:double-double-hi b) + (kernel:double-double-lo b))) ) ; end progn View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/eb37d9a9b059843f874542c... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/eb37d9a9b059843f874542c... You're receiving this email because of your account on gitlab.common-lisp.net.