January 2026 - cmucl-cvs - mailman3.common-lisp.net

[Git][cmucl/cmucl][master] 2 commits: Fix #428: Add some constants for common values in irrat.lisp
by Raymond Toy (＠rtoy) 20 Jan '26

20 Jan '26

Raymond Toy pushed to branch master at cmucl / cmucl Commits: 6d264493 by Raymond Toy at 2026-01-20T07:02:39-08:00 Fix #428: Add some constants for common values in irrat.lisp - - - - - ea927341 by Raymond Toy at 2026-01-20T07:02:40-08:00 Merge branch 'issue-428-add-constants-for-irrat' into 'master' Fix #428: Add some constants for common values in irrat.lisp Closes #428 See merge request cmucl/cmucl!344 - - - - - 3 changed files: - src/code/irrat.lisp - src/general-info/release-22a.md - src/i18n/locale/cmucl.pot Changes: ===================================== src/code/irrat.lisp ===================================== @@ -24,6 +24,18 @@ (defconstant pi 3.14159265358979323846264338327950288419716939937511L0) ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0) +(defconstant +log-2d0+ + (log 2d0) + "Natural log of 2d0") + +(defconstant +log-2w0+ + (log 2w0) + "Natural log of 2w0") + +(defconstant +log2-e+ + (log (exp 1d0) 2d0) + "Base 2 logarithm of e") + ;;; Make these INLINE, since the call to C is at least as compact as a Lisp ;;; call, and saves number consing to boot. ;;; @@ -488,14 +500,14 @@ complex) (if (and (zerop base) (plusp (realpart power))) (* base power) - (exp (* power (* (log2 base 1w0) (log 2w0)))))) + (exp (* power (* (log2 base 1w0) +log-2w0+))))) (((foreach fixnum (or bignum ratio) single-float double-float) (foreach (complex double-float))) ;; Result should have double-float accuracy. Use log2 in ;; case the base won't fit in a double-float. (if (and (zerop base) (plusp (realpart power))) (* base power) - (exp (* power (* (log2 base) (log 2d0)))))) + (exp (* power (* (log2 base) +log-2d0+))))) ((double-float (foreach (complex rational) (complex single-float))) (if (and (zerop base) (plusp (realpart power))) @@ -508,7 +520,7 @@ ;; in case the base won't fit in a double-float. (if (and (zerop base) (plusp (realpart power))) (* base power) - (exp (* power (* (log2 base 1w0) (log 2w0)))))) + (exp (* power (* (log2 base 1w0) +log-2w0+))))) (((foreach fixnum (or bignum ratio) single-float) (foreach (complex rational) (complex single-float))) (if (and (zerop base) (plusp (realpart power))) @@ -566,18 +578,18 @@ ;; log(2), with the precision specified by the type of F (number-dispatch ((f real)) ((double-float) - #.(log 2d0)) + +log-2d0+) #+double-double ((double-double-float) - #.(log 2w0)))) + +log-2w0+))) (log-2-pi (f) ;; log(pi), with the precision specified by the type of F (number-dispatch ((f real)) ((double-float) - #.(/ pi (log 2d0))) + #.(/ pi +log-2d0+)) #+double-double ((double-double-float) - #.(/ dd-pi (log 2w0))))) + #.(/ dd-pi +log-2w0+)))) (log1p (x) ;; log(1+x), with the precision specified by the type of ;; X @@ -827,7 +839,7 @@ (if (minusp number) (complex (coerce (log (- number)) 'single-float) (coerce pi 'single-float)) - (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float))) + (coerce (/ (log2 number) +log2-e+) 'single-float))) ((ratio) (if (minusp number) (complex (coerce (log (- number)) 'single-float) @@ -852,7 +864,7 @@ (coerce (%log1p (coerce (- number 1) 'double-float)) 'single-float) (coerce (/ (- (log2 top) (log2 bot)) - #.(log (exp 1d0) 2d0)) + +log2-e+) 'single-float))))) (((foreach single-float double-float)) ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)? @@ -1393,7 +1405,7 @@ This is for use with J /= 0 only when |z| is huge." (let ((t0 #.(/ 1 (sqrt 2.0d0))) (t1 1.2d0) (t2 3d0) - (ln2 #.(log 2d0)) + (ln2 +log-2d0+) (x (float (realpart z) 1.0d0)) (y (float (imagpart z) 1.0d0))) (multiple-value-bind (rho k) ===================================== src/general-info/release-22a.md ===================================== @@ -28,6 +28,7 @@ public domain. * #373: Add macros to handle temp files and directories * #387: Update CLX from upstream to version from 2024-09-11 * #408: Add C functions to get user name and home dir + * #428: Add some constants for common values in irrat.lisp * #437: Use `realpath` to find the lisp executable * #444: `unix-stat` and friends return wrong timestamp * #449: Touch 64-bit timestamp.txt with desired time for testing ===================================== src/i18n/locale/cmucl.pot ===================================== @@ -5008,6 +5008,18 @@ msgid "" " (= x (float (rationalize x) x))" msgstr "" +#: src/code/irrat.lisp +msgid "Natural log of 2d0" +msgstr "" + +#: src/code/irrat.lisp +msgid "Natural log of 2w0" +msgstr "" + +#: src/code/irrat.lisp +msgid "Base 2 logarithm of e" +msgstr "" + #: src/code/irrat.lisp msgid "Return e raised to the power NUMBER." msgstr "" View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3f733bae60b20dadee5f4f… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3f733bae60b20dadee5f4f… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-461-cdiv-use-4-doubles] 4 commits: Fix bug in deftransform / of complex floats
by Raymond Toy (＠rtoy) 20 Jan '26

20 Jan '26

Raymond Toy pushed to branch issue-461-cdiv-use-4-doubles at cmucl / cmucl Commits: 3353b960 by Raymond Toy at 2026-01-19T16:36:43-08:00 Fix bug in deftransform / of complex floats We forgot to quote the body of the code that is the expansion of the deftransform. This doesn't show up until you try to compile some code that calls the deftransform. - - - - - 96037b55 by Raymond Toy at 2026-01-19T16:37:59-08:00 Need to make cdiv-single-float take 4 floats instead of 2 complex - - - - - 3f733bae by Raymond Toy at 2026-01-19T16:54:46-08:00 Fix critical typo in deftransform for / for complex floats We inadvertently left off the quote for the code that the deftransform is supposed to use. Without that the compiler produces an error when the deftransform would be applicable. This bug was introduced when the MR for #459 was landed. - - - - - f7c31886 by Raymond Toy at 2026-01-19T16:56:53-08:00 Merge branch 'master' into issue-461-cdiv-use-4-doubles - - - - - 2 changed files: - src/code/numbers.lisp - src/compiler/float-tran.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -811,13 +811,13 @@ ;; 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) +(defun cdiv-single-float (xr xi yr yi) "Accurate division of complex single-float numbers x and y." - (declare (type (complex single-float) x y)) - (let ((a (float (realpart x) 1d0)) - (b (float (imagpart x) 1d0)) - (c (float (realpart y) 1d0)) - (d (float (imagpart y) 1d0))) + (declare (single-float xr xi yr yi)) + (let ((a (float xr 1d0)) + (b (float xi 1d0)) + (c (float yr 1d0)) + (d (float yi 1d0))) (cond ((< (abs c) (abs d)) (let* ((r (/ c d)) (denom (+ (* c r) d)) @@ -858,7 +858,10 @@ (((complex single-float) (foreach (complex rational) (complex single-float))) - (cdiv-single-float x (coerce y '(complex single-float)))) + (cdiv-single-float (realpart x) + (imagpart x) + (coerce (realpart y) 'single-float) + (coerce (imagpart y) 'single-float))) (((complex double-float) (foreach (complex rational) (complex single-float) (complex double-float))) (cdiv-double-float (realpart x) @@ -868,8 +871,10 @@ (((foreach integer ratio single-float (complex rational)) (complex single-float)) - (cdiv-single-float (coerce x '(complex single-float)) - y)) + (cdiv-single-float (coerce (realpart x) 'single-float) + (coerce (imagpart x) 'single-float) + (realpart y) + (imagpart y))) (((foreach integer ratio single-float double-float (complex rational) (complex single-float)) ===================================== src/compiler/float-tran.lisp ===================================== @@ -1807,10 +1807,10 @@ (macrolet ((frob (type name) `(deftransform / ((x y) ((complex ,type) (complex ,type)) *) - (,name (realpart x) - (imagpart x) - (realpart y) - (imagpart y))))) + '(,name (realpart x) + (imagpart x) + (realpart y) + (imagpart y))))) (frob double-float kernel::cdiv-double-float) (frob single-float kernel::cdiv-single-float)) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/fae66d26c66d2632e670b8… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/fae66d26c66d2632e670b8… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][master] Fix critical typo in deftransform for / for complex floats
by Raymond Toy (＠rtoy) 20 Jan '26

20 Jan '26

Raymond Toy pushed to branch master at cmucl / cmucl Commits: 3f733bae by Raymond Toy at 2026-01-19T16:54:46-08:00 Fix critical typo in deftransform for / for complex floats We inadvertently left off the quote for the code that the deftransform is supposed to use. Without that the compiler produces an error when the deftransform would be applicable. This bug was introduced when the MR for #459 was landed. - - - - - 1 changed file: - src/compiler/float-tran.lisp Changes: ===================================== src/compiler/float-tran.lisp ===================================== @@ -1807,7 +1807,7 @@ (macrolet ((frob (type name) `(deftransform / ((x y) ((complex ,type) (complex ,type)) *) - (,name x y)))) + '(,name x y)))) (frob double-float kernel::cdiv-double-float) (frob single-float kernel::cdiv-single-float)) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/3f733bae60b20dadee5f4f0… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/3f733bae60b20dadee5f4f0… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-461-cdiv-use-4-doubles] 7 commits: Update cmucl.pot for new docstrings
by Raymond Toy (＠rtoy) 19 Jan '26

19 Jan '26

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/eb37d9a9b059843f874542… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/eb37d9a9b059843f874542… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][master] 2 commits: Fix #459: More accurate complex division for double-double-float
by Raymond Toy (＠rtoy) 19 Jan '26

19 Jan '26

Raymond Toy pushed to branch master at cmucl / cmucl Commits: 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 - - - - - 5 changed files: - src/code/numbers.lisp - src/compiler/float-tran-dd.lisp - src/general-info/release-22a.md - src/i18n/locale/cmucl.pot - tests/float.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -602,34 +602,40 @@ ;; 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+)) - -;; 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. -(declaim (ext:start-block cdiv-double-float cdiv-single-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). +;; +;; 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-/)) + +;; 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. +(eval-when (:compile-toplevel) +(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)))) + `(progn + (defun ,compreal (a b c d r tt) + (declare (,type a b c d r tt) + ,@opt) + ;; 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 @@ -639,124 +645,174 @@ ;; = (a + b*r)/(c + d*r). ;; ;; Thus tt = (c + d*r). - (cond ((>= (abs r) +cdiv-rmin+) + (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. + ;; r is very small 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)) + (defun ,subinternal (a b c d) + (declare (,type a b c d) + ,@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). - (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))))) + ;; 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 ((+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 + ;; multiplied by 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 + :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. @@ -825,13 +881,14 @@ (((complex double-double-float) (foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float))) - ;; We should do something better for double-double floats. - (cdiv-generic x y)) + (cdiv-double-double-float x + (coerce y '(complex double-double-float)))) (((foreach integer ratio single-float double-float double-double-float (complex rational) (complex single-float) (complex double-float)) (complex double-double-float)) - (cdiv-generic x y)) + (cdiv-double-double-float (coerce x '(complex double-double-float)) + y)) (((foreach integer ratio single-float double-float double-double-float) (complex rational)) ===================================== src/compiler/float-tran-dd.lisp ===================================== @@ -576,6 +576,10 @@ (truly-the ,(type-specifier (node-derived-type node)) (kernel:%make-double-double-float hi lo)))) +(deftransform / ((a b) ((complex vm::double-double-float) (complex vm::double-double-float)) + *) + `(kernel::cdiv-double-double-float a b)) + (declaim (inline sqr-d)) (defun sqr-d (a) "Square" ===================================== src/general-info/release-22a.md ===================================== @@ -38,6 +38,9 @@ public domain. using Baudin and Smith's robust complex division algorithm with improvements by Patrick McGehearty. * #458: Spurious overflow in double-double-float multiply + * #459: Improve accuracy for division of complex + double-double-floats. The same algorithm is used here as + for #456. * Other changes: * Improvements to the PCL implementation of CLOS: * Changes to building procedure: ===================================== src/i18n/locale/cmucl.pot ===================================== @@ -4389,6 +4389,10 @@ msgstr "" msgid "Accurate division of complex double-float numbers x and y." msgstr "" +#: src/code/numbers.lisp +msgid "Accurate division of complex double-double-float numbers x and y." +msgstr "" + #: src/code/numbers.lisp msgid "Accurate division of complex single-float numbers x and y." msgstr "" ===================================== tests/float.lisp ===================================== @@ -343,6 +343,249 @@ (assert-true (typep new-mode 'x86::float-modes)) (assert-equal new-mode (setf (x86::x87-floating-point-modes) new-mode)))) +;; Rudimentary code to read C %a formatted numbers that look like +;; "-0x1.c4dba4ba1ee79p-620". We assume STRING is exactly in this +;; format. No error-checking is done. +(defun parse-hex-float (string) + (let* ((sign (if (char= (aref string 0) #\-) + -1 + 1)) + (dot-posn (position #\. string)) + (p-posn (position #\p string)) + (lead (parse-integer string :start (1- dot-posn) :end dot-posn)) + (frac (parse-integer string :start (1+ dot-posn) :end p-posn :radix 16)) + (exp (parse-integer string :start (1+ p-posn)))) + (* sign + (scale-float (float (+ (ash lead 52) + frac) + 1d0) + (- exp 52))))) + +;; Relative error in terms of bits of accuracy. This is the +;; definition used by Baudin and Smith. A result of 53 means the two +;; numbers have identical bits. For complex numbers, we use the min +;; of the bits of accuracy of the real and imaginary parts. +(defun rel-err (computed expected) + (flet ((rerr (c e) + (let ((diff (abs (- c e)))) + (if (zerop diff) + (float-digits diff) + (floor (- (log (/ diff (abs e)) 2d0))))))) + (min (rerr (realpart computed) (realpart expected)) + (rerr (imagpart computed) (imagpart expected))))) + +(defun do-cdiv-test (x y z-true expected-rel) + (let* ((z (/ x y)) + (rel (rel-err z z-true))) + (assert-equal expected-rel + rel + x y z z-true rel))) + +(defun do-cdiv-dd-test (x y z-true expected-rel) + (flet ((compute-true (a b) + ;; Convert a and b to complex rationals, do the + ;; division and convert back to get the true + ;; expected result. + (coerce + (/ (complex (rational (realpart a)) + (rational (imagpart a))) + (complex (rational (realpart b)) + (rational (imagpart b)))) + '(complex ext:double-double-float)))) + (let* ((z (/ (coerce x '(complex ext:double-double-float)) + (coerce y '(complex ext:double-double-float)))) + (z-true (compute-true x y)) + (rel (rel-err z z-true))) + (assert-equal expected-rel + rel + k x y z z-true rel)))) + +;; Issue #456: improve accuracy of division of complex double-floats. +;; +;; Tests for complex division. Tests 1-10 are from Baudin and Smith. +;; Test 11 is an example from Maxima. Test 12 is an example from the +;; ansi-tests where (/ z z) didn't produce exactly 1. Tests 13-16 are +;; for examples for improvement iterations 1-4 from McGehearty. +(macrolet + ((frob (name x y z-true rel dd-rel) + `(define-test ,name + (:tag :issues) + (do-cdiv-test ,x ,y ,z-true ,rel) + (do-cdiv-dd-test ,x ,y ,z-true ,dd-rel)))) + ;; First cases are from Baudin and Smith + ;; 1 + (frob cdiv.baudin-case.1 + (complex 1d0 1d0) + (complex 1d0 (scale-float 1d0 1023)) + (complex (scale-float 1d0 -1023) + (scale-float -1d0 -1023)) + 53 + 106) + ;; 2 + (frob cdiv.baudin-case.2 + (complex 1d0 1d0) + (complex (scale-float 1d0 -1023) (scale-float 1d0 -1023)) + (complex (scale-float 1d0 1023) 0) + 53 + 106) + ;; 3 + (frob cdiv.baudin-case.3 + (complex (scale-float 1d0 1023) (scale-float 1d0 -1023)) + (complex (scale-float 1d0 677) (scale-float 1d0 -677)) + (complex (scale-float 1d0 346) (scale-float -1d0 -1008)) + 53 + 106) + ;; 4 + (frob cdiv.baudin-case.4.overflow + (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) + (complex 1d0 1d0) + (complex (scale-float 1d0 1023) 0) + 53 + 106) + ;; 5 + (frob cdiv.baudin-case.5.underflow-ratio + (complex (scale-float 1d0 1020) (scale-float 1d0 -844)) + (complex (scale-float 1d0 656) (scale-float 1d0 -780)) + (complex (scale-float 1d0 364) (scale-float -1d0 -1072)) + 53 + 106) + ;; 6 + (frob cdiv.baudin-case.6.underflow-realpart + (complex (scale-float 1d0 -71) (scale-float 1d0 1021)) + (complex (scale-float 1d0 1001) (scale-float 1d0 -323)) + (complex (scale-float 1d0 -1072) (scale-float 1d0 20)) + 53 + 106) + ;; 7 + (frob cdiv.baudin-case.7.overflow-both-parts + (complex (scale-float 1d0 -347) (scale-float 1d0 -54)) + (complex (scale-float 1d0 -1037) (scale-float 1d0 -1058)) + (complex 3.898125604559113300d289 8.174961907852353577d295) + 53 + 106) + ;; 8 + (frob cdiv.baudin-case.8 + (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074)) + (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074)) + (complex 0.6d0 0.2d0) + 53 + 106) + ;; 9 + (frob cdiv.baudin-case.9 + (complex (scale-float 1d0 1015) (scale-float 1d0 -989)) + (complex (scale-float 1d0 1023) (scale-float 1d0 1023)) + (complex 0.001953125d0 -0.001953125d0) + 53 + 106) + ;; 10 + (frob cdiv.baudin-case.10.improve-imagpart-accuracy + (complex (scale-float 1d0 -622) (scale-float 1d0 -1071)) + (complex (scale-float 1d0 -343) (scale-float 1d0 -798)) + (complex 1.02951151789360578d-84 6.97145987515076231d-220) + 53 + 106) + ;; 11 + ;; + ;; From Maxima. This was from a (private) email where Maxima used + ;; CL:/ to compute the ratio but was not very accurate. + (frob cdiv.maxima-case + #c(5.43d-10 1.13d-100) + #c(1.2d-311 5.7d-312) + #c(3.691993880674614517999740937026568563794896024143749539711267954d301 + -1.753697093319947872394996242210428954266103103602859195409591583d301) + 52 + 107) + ;; 12 + ;; + ;; Found by ansi tests. z/z should be exactly 1. + (frob cdiv.ansi-test-z/z + #c(1.565640716292489d19 0.0d0) + #c(1.565640716292489d19 0.0d0) + #c(1d0 0) + 53 + 106) + ;; 13 + ;; Iteration 1. Without this, we would instead return + ;; + ;; (complex (parse-hex-float "0x1.ba8df8075bceep+155") + ;; (parse-hex-float "-0x1.a4ad6329485f0p-895")) + ;; + ;; whose imaginary part is quite a bit off. + (frob cdiv.mcgehearty-iteration.1 + (complex (parse-hex-float "0x1.73a3dac1d2f1fp+509") + (parse-hex-float "-0x1.c4dba4ba1ee79p-620")) + (complex (parse-hex-float "0x1.adf526c249cf0p+353") + (parse-hex-float "0x1.98b3fbc1677bbp-697")) + (complex (parse-hex-float "0x1.BA8DF8075BCEEp+155") + (parse-hex-float "-0x1.A4AD628DA5B74p-895")) + 53 + 106) + ;; 14 + ;; Iteration 2. + (frob cdiv.mcgehearty-iteration.2 + (complex (parse-hex-float "-0x0.000000008e4f8p-1022") + (parse-hex-float "0x0.0000060366ba7p-1022")) + (complex (parse-hex-float "-0x1.605b467369526p-245") + (parse-hex-float "0x1.417bd33105808p-256")) + (complex (parse-hex-float "0x1.cde593daa4ffep-810") + (parse-hex-float "-0x1.179b9a63df6d3p-799")) + 52 + 106) + ;; 15 + ;; Iteration 3 + (frob cdiv.mcgehearty-iteration.3 + (complex (parse-hex-float "0x1.cb27eece7c585p-355 ") + (parse-hex-float "0x0.000000223b8a8p-1022")) + (complex (parse-hex-float "-0x1.74e7ed2b9189fp-22") + (parse-hex-float "0x1.3d80439e9a119p-731")) + (complex (parse-hex-float "-0x1.3b35ed806ae5ap-333") + (parse-hex-float "-0x0.05e01bcbfd9f6p-1022")) + 53 + 106) + ;; 16 + ;; Iteration 4 + (frob cdiv.mcgehearty-iteration.4 + (complex (parse-hex-float "-0x1.f5c75c69829f0p-530") + (parse-hex-float "-0x1.e73b1fde6b909p+316")) + (complex (parse-hex-float "-0x1.ff96c3957742bp+1023") + (parse-hex-float "0x1.5bd78c9335899p+1021")) + (complex (parse-hex-float "-0x1.423c6ce00c73bp-710") + (parse-hex-float "0x1.d9edcf45bcb0ep-708")) + 52 + 106)) + +(define-test complex-division.misc + (:tag :issue) + (let ((num '(1 + 1/2 + 1.0 + 1d0 + #c(1 2) + #c(1.0 2.0) + #c(1d0 2d0) + #c(1w0 2w0)))) + ;; Try all combinations of divisions of different types. This is + ;; primarily to test that we got all the numeric contagion cases + ;; for division in CL:/. + (dolist (x num) + (dolist (y num) + (assert-true (/ x y) + x y))))) + +(define-test complex-division.single + (:tag :issues) + (let* ((x #c(1 2)) + (y (complex (expt 2 127) (expt 2 127))) + (expected (coerce (/ x y) + '(complex single-float)))) + ;; A naive implementation of complex division would cause an + ;; overflow in computing the denominator. + (assert-equal expected + (/ (coerce x '(complex single-float)) + (coerce y '(complex single-float))) + x + y))) + ;; Issue #458 (define-test dd-mult-overflow (:tag :issues) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/67dc5c9deb9a47cae250d0… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/67dc5c9deb9a47cae250d0… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] Clean up a comment
by Raymond Toy (＠rtoy) 19 Jan '26

19 Jan '26

Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div at cmucl / cmucl Commits: 466dee14 by Raymond Toy at 2026-01-19T11:26:01-08:00 Clean up a comment - - - - - 1 changed file: - src/code/numbers.lisp Changes: ===================================== src/code/numbers.lisp ===================================== @@ -751,12 +751,12 @@ (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 is always multipled by a power of 2. It's + ;; multiplied by 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. View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/466dee1428a1b6ed3a439e7… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/466dee1428a1b6ed3a439e7… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] 2 commits: Let bind the constants as needed in each routine
by Raymond Toy (＠rtoy) 19 Jan '26

19 Jan '26

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/7ef6184520ac5201dd878e… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/7ef6184520ac5201dd878e… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-459-more-accurate-dd-complex-div] 5 commits: Update cmucl.pot for new docstrings
by Raymond Toy (＠rtoy) 19 Jan '26

19 Jan '26

Raymond Toy pushed to branch issue-459-more-accurate-dd-complex-div 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 - - - - - bde45741 by Raymond Toy at 2026-01-18T15:53:04-08:00 Merge branch 'master' into issue-459-more-accurate-dd-complex-div - - - - - 7ef61845 by Raymond Toy at 2026-01-19T06:41:38-08:00 Update cmucl.pot - - - - - 1 changed file: - src/compiler/float-tran-dd.lisp Changes: ===================================== src/compiler/float-tran-dd.lisp ===================================== @@ -678,6 +678,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) @@ -685,10 +697,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/3b75ab818edd40c13b9e1f… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3b75ab818edd40c13b9e1f… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][master] 2 commits: Fix #463: Handle double-double-float and double-float comparisons
by Raymond Toy (＠rtoy) 18 Jan '26

18 Jan '26

Raymond Toy pushed to branch master at cmucl / cmucl Commits: 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 - - - - - 1 changed file: - src/compiler/float-tran-dd.lisp Changes: ===================================== src/compiler/float-tran-dd.lisp ===================================== @@ -674,6 +674,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) @@ -681,10 +693,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/336b9322330eee9a4b44d9… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/336b9322330eee9a4b44d9… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0

[Git][cmucl/cmucl][issue-428-add-constants-for-irrat] 2 commits: Add constant for log base 2 of e
by Raymond Toy (＠rtoy) 18 Jan '26

18 Jan '26

Raymond Toy pushed to branch issue-428-add-constants-for-irrat at cmucl / cmucl Commits: 46e88b6a by Raymond Toy at 2026-01-18T07:51:21-08:00 Add constant for log base 2 of e - - - - - 688490f6 by Raymond Toy at 2026-01-18T07:57:48-08:00 Update release notes and cmucl.pot for docstrings - - - - - 3 changed files: - src/code/irrat.lisp - src/general-info/release-22a.md - src/i18n/locale/cmucl.pot Changes: ===================================== src/code/irrat.lisp ===================================== @@ -32,7 +32,10 @@ (log 2w0) "Natural log of 2w0") - +(defconstant +log2-e+ + (log (exp 1d0) 2d0) + "Base 2 logarithm of e") + ;;; Make these INLINE, since the call to C is at least as compact as a Lisp ;;; call, and saves number consing to boot. ;;; @@ -836,7 +839,7 @@ (if (minusp number) (complex (coerce (log (- number)) 'single-float) (coerce pi 'single-float)) - (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float))) + (coerce (/ (log2 number) +log2-e+) 'single-float))) ((ratio) (if (minusp number) (complex (coerce (log (- number)) 'single-float) @@ -861,7 +864,7 @@ (coerce (%log1p (coerce (- number 1) 'double-float)) 'single-float) (coerce (/ (- (log2 top) (log2 bot)) - #.(log (exp 1d0) 2d0)) + +log2-e+) 'single-float))))) (((foreach single-float double-float)) ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)? ===================================== src/general-info/release-22a.md ===================================== @@ -28,6 +28,7 @@ public domain. * #373: Add macros to handle temp files and directories * #387: Update CLX from upstream to version from 2024-09-11 * #408: Add C functions to get user name and home dir + * #428: Add some constants for common values in irrat.lisp * #437: Use `realpath` to find the lisp executable * #444: `unix-stat` and friends return wrong timestamp * #449: Touch 64-bit timestamp.txt with desired time for testing ===================================== src/i18n/locale/cmucl.pot ===================================== @@ -5012,6 +5012,10 @@ msgstr "" msgid "Natural log of 2w0" msgstr "" +#: src/code/irrat.lisp +msgid "Base 2 logarithm of e" +msgstr "" + #: src/code/irrat.lisp msgid "Return e raised to the power NUMBER." msgstr "" View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3451a113880c01328881c9… -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/3451a113880c01328881c9… You're receiving this email because of your account on gitlab.common-lisp.net.

1 0