Raymond Toy pushed to branch rtoy-print-using-ryu at cmucl / cmucl Commits: f65ff9eb by Raymond Toy at 2026-05-25T17:34:15-07:00 Fix #504: Do correct rounding in scale-float-maybe-underflow Previously, when dealing with denormals, we right shifted the bits so we were truncating instead of rounding to even as expected. Instead of trying to round the bits ourselves, build a normal-range float that is not a denormal by increasing the exponent outside the denormal range. `scale-float` is fine for this and doesn't introduce any round-off. Then multiply by a power of two to get the correct denormal. Since it's a power of two, there's no rounding except when we get a denormal. This is the rounding we want. The kludge in `float-ratio-float` isn't needed because it was trying to fix the denormal case. It incorrectly assumed that `scale-float` was rounding correctly. As we know now, it wasn't. Updated the test `float-ratio.double` which had the incorrect expected value. It should round to twice the least-positive-double-float. Because we've corrected the reader, cdiv.maxima-case was getting slightly different floats. Compute the expected value by converting the args to complex rationals and doing a rational division to get the fully accurate result. Added additional test for `scale-float-maybe-underflow` to make sure we get the rounding correct. And verify that we read 1.1e-44 and 1.1d-322 correctly. The former value showed up in an ansi-test failure before this fix. - - - - - 63c190db by Raymond Toy at 2026-05-25T17:54:42-07:00 Add some declarations and simplify comment Add some declarations to scale-float-maybe-underflow to help the compiler generate a little bit better code. Simplify the comment about how we do scaling for denormals. - - - - - 06ea24af by Raymond Toy at 2026-05-26T07:22:37-07:00 Change declaration of exp arg In `scale-float-maybe-underflow`, the exponent (`exp`) isn't limited to the double-float-exponent range. It can basically be any fixnum. This can be seen with `(scale-float 1d0 -9999)`. `scale-float-maybe-underflow` sees -9999, well outside the range of double-float-exponent type. - - - - - a5b2bbcc by Raymond Toy at 2026-05-26T07:37:03-07:00 Update release notes. No other changes. - - - - - af544ca1 by Raymond Toy at 2026-05-26T07:38:14-07:00 Merge branch 'issue-504-read-denormals-with-rounding' into rtoy-print-using-ryu - - - - - 3 changed files: - src/code/float.lisp - src/general-info/release-22a.md - tests/float.lisp Changes: ===================================== src/code/float.lisp ===================================== @@ -885,40 +885,57 @@ ;;; denormalized or underflows to 0. ;;; (defun scale-float-maybe-underflow (x exp) + (declare (type (or single-float double-float) x) + (fixnum exp)) (multiple-value-bind (sig old-exp) - (integer-decode-float x) + (integer-decode-float x) (let* ((digits (float-digits x)) + (1+digits (1+ digits)) (new-exp (+ exp old-exp digits (etypecase x (single-float vm:single-float-bias) (double-float vm:double-float-bias)))) (sign (if (minusp (float-sign x)) 1 0))) (cond - ((< new-exp - (etypecase x - (single-float vm:single-float-normal-exponent-min) - (double-float vm:double-float-normal-exponent-min))) - (when (vm:current-float-trap :inexact) - (error 'floating-point-inexact :operation 'scale-float - :operands (list x exp))) - (when (vm:current-float-trap :underflow) - (error 'floating-point-underflow :operation 'scale-float - :operands (list x exp))) - (let ((shift (1- new-exp))) - ;; Is it necessary to have this IF here? Is there any case - ;; where (ash sig shift) won't return 0 when - ;; shift < -(digits-1)? - (if (< shift (- (1- digits))) + ((< new-exp + (etypecase x + (single-float vm:single-float-normal-exponent-min) + (double-float vm:double-float-normal-exponent-min))) + (when (vm:current-float-trap :inexact) + (error 'floating-point-inexact :operation 'scale-float + :operands (list x exp))) + (when (vm:current-float-trap :underflow) + (error 'floating-point-underflow :operation 'scale-float + :operands (list x exp))) + ;; To round correctly, let the hardware multiplier do the + ;; rounding: build a normal float whose stored exponent is + ;; bumped up by 1+DIGITS (which puts it safely in the normal + ;; range), then multiply by 2^-(1+DIGITS). The multiplier is + ;; an exact power of two, so the multiplication is exact + ;; apart from the unavoidable rounding step that expresses + ;; the product as a denormal, which the FPU performs in the + ;; current rounding mode. If the bumped exponent is zero or + ;; negative the bumped float would itself be a denormal -- + ;; losing the implicit 1 bit of SIG -- so handle that case + ;; explicitly by returning signed zero. + (let ((bumped-exp (+ new-exp 1+digits))) + (cond + ((<= bumped-exp 0) (etypecase x (single-float (single-from-bits sign 0 0)) - (double-float (double-from-bits sign 0 0))) + (double-float (double-from-bits sign 0 0)))) + (t (etypecase x - (single-float (single-from-bits sign 0 (ash sig shift))) - (double-float (double-from-bits sign 0 (ash sig shift))))))) - (t - (etypecase x - (single-float (single-from-bits sign new-exp sig)) - (double-float (double-from-bits sign new-exp sig)))))))) + (single-float + (* (single-from-bits sign bumped-exp sig) + (scale-float 1f0 (- 1+digits)))) + (double-float + (* (double-from-bits sign bumped-exp sig) + (scale-float 1d0 (- 1+digits))))))))) + (t + (etypecase x + (single-float (single-from-bits sign new-exp sig)) + (double-float (double-from-bits sign new-exp sig)))))))) ;;; SCALE-FLOAT-MAYBE-OVERFLOW -- Internal @@ -1150,27 +1167,11 @@ (format t "2: f0, f1 = ~A~%" f0) (format t " scale = ~A~%" scale) (format t "scale-float f0 = ~A~%" (scale-float f0 scale))) - (let ((min-exponent - ;; Compute the min (unbiased) exponent - (ecase format - (single-float - (- vm:single-float-normal-exponent-min - vm:single-float-bias - vm:single-float-digits)) - (double-float - (- vm:double-float-normal-exponent-min - vm:double-float-bias - vm:double-float-digits))))) - ;; F0 is always between 0.5 and 1. If - ;; SCALE is the min exponent, we have a - ;; denormal number just less than the - ;; least-positive float. We want to - ;; return the least-positive-float so - ;; multiply F0 by 2 (without adjusting - ;; SCALE) to get the nearest float. - (if (= scale min-exponent) - (scale-float (* 2 f0) scale) - (scale-float f0 scale)))))))) + ;; SCALE-FLOAT-MAYBE-UNDERFLOW correctly rounds + ;; to nearest when constructing a denormal result, + ;; so just call SCALE-FLOAT here without any + ;; boundary-case fixup. + (scale-float f0 scale)))))) (floatit (bits) (let ((sign (if plusp 0 1))) (case format ===================================== src/general-info/release-22a.md ===================================== @@ -58,6 +58,8 @@ public domain. * #463: `double-double-float` is missing comparison operations between `double-double-float` and `double-float` * #474: Add functions to print and parse C-style hex floats. + * #504: Do correct rounding in `scale-float-maybe-underflow`. + This was causing some denormals to be read incorrectly. * Other changes: * Improvements to the PCL implementation of CLOS: * Changes to building procedure: ===================================== tests/float.lisp ===================================== @@ -179,7 +179,9 @@ (kernel::float-ratio-float (* 4 expo) 'double-float)) (assert-equal least-positive-double-float (kernel::float-ratio-float (* 494/100 expo) 'double-float)) - (assert-equal least-positive-double-float + ;; 988/100*10^-324 is very close to 2*least-positive (the exact ratio + ;; is 1.9997 * least-positive), so it rounds to 2*least-positive. + (assert-equal (* 2 least-positive-double-float) (kernel::float-ratio-float (* 988/100 expo) 'double-float))))) (define-test reader-error.small-single-floats @@ -678,8 +680,11 @@ (frob cdiv.maxima-case #c(5.43d-10 1.13d-100) #c(1.2d-311 5.7d-312) - #c(3.691993880674614517999740937026568563794896024143749539711267954d301 - -1.753697093319947872394996242210428954266103103602859195409591583d301) + ;; Compute the expected value using rational arithmetic after + ;; converting the complex numbers above to the equivalent + ;; complex rationals. + (/ (complex (rational 5.43d-10) (rational 1.13d-100)) + (complex (rational 1.2d-311) (rational 5.7d-312))) 52) ;; 12 ;; @@ -766,3 +771,85 @@ (coerce y '(complex single-float))) x y))) + +(define-test scale-float-underflow-rounding.single + (:tag :issues) + ;; SCALE-FLOAT into the denormal range must round to nearest, ties to + ;; even, instead of truncating the discarded bits. Each (X EXP BITS) + ;; triple gives a normal X, an exponent EXP to scale by, and the IEEE + ;; bits of the expected single-float result. + (ext:with-float-traps-masked (:underflow :inexact) + (dolist (case (list + ;; 1.7f0 * 2^-149: between denormals 1 and 2, closer to 2. + (list 1.7f0 -149 #x00000002) + ;; 1.5f0 * 2^-149: halfway between denormals 1 and 2; + ;; ties round to even -> 2. + (list 1.5f0 -149 #x00000002) + ;; 1.1f0 * 2^-149: closer to denormal 1. + (list 1.1f0 -149 #x00000001) + ;; 1.0001f0 * 2^-150: just above halfway between 0 and + ;; smallest denormal; rounds up to 1. + (list 1.0001f0 -150 #x00000001) + ;; 1.0f0 * 2^-150: exactly halfway between 0 and smallest + ;; denormal; ties round to even -> 0. + (list 1.0f0 -150 #x00000000) + ;; Largest single < 2 scaled by 2^-127: rounding carries + ;; into the implicit-1 position and produces the smallest + ;; normal number. + (list (kernel:make-single-float #x3fffffff) + -127 #x00800000))) + (destructuring-bind (x exp bits) case + (let ((result (scale-float x exp))) + (assert-equal bits (kernel:single-float-bits result) + x exp result)))))) + +(define-test scale-float-underflow-rounding.double + (:tag :issues) + ;; Like SCALE-FLOAT-UNDERFLOW-ROUNDING.SINGLE but for double-floats. + ;; Each (X EXP HI LO) gives a normal X, an exponent EXP, and the IEEE + ;; high and low bits of the expected double-float result. + (ext:with-float-traps-masked (:underflow :inexact) + (dolist (case (list + ;; 1.7d0 * 2^-1074: between denormals 1 and 2, closer to 2. + (list 1.7d0 -1074 0 2) + ;; 1.5d0 * 2^-1074: tie, rounds to even -> 2. + (list 1.5d0 -1074 0 2) + ;; 1.1d0 * 2^-1074: closer to denormal 1. + (list 1.1d0 -1074 0 1) + ;; 1.0001d0 * 2^-1075: just above halfway, rounds up. + (list 1.0001d0 -1075 0 1) + ;; 1.0d0 * 2^-1075: tie at the bottom, rounds to even -> 0. + (list 1.0d0 -1075 0 0) + ;; Largest double < 2 scaled by 2^-1023: rounding carries + ;; into the implicit-1 position and produces the smallest + ;; normal number. + (list (kernel:make-double-float #x3fffffff #xffffffff) + -1023 #x00100000 0))) + (destructuring-bind (x exp hi lo) case + (let ((result (scale-float x exp))) + (assert-equal hi (kernel:double-float-high-bits result) + x exp result) + (assert-equal lo (kernel:double-float-low-bits result) + x exp result)))))) + +(define-test scale-float-underflow-rounding.reader-single + (:tag :issues) + ;; The reader uses FLOAT-RATIO-FLOAT, which calls SCALE-FLOAT and + ;; hence SCALE-FLOAT-MAYBE-UNDERFLOW on denormal results. Reading + ;; small float literals must therefore also round to nearest. + (ext:with-float-traps-masked (:underflow :inexact) + ;; 1.1e-44 is closer to 8 * least-positive-single-float than to 7. + (assert-equal #x00000008 + (kernel:single-float-bits 1.1e-44)) + ;; 1.121e-44 (essentially the IEEE representation of 1.1e-44) reads + ;; as the same denormal. + (assert-equal #x00000008 + (kernel:single-float-bits 1.121e-44)))) + +(define-test scale-float-underflow-rounding.reader-double + (:tag :issues) + ;; Like the single-float reader test, in the double denormal range. + ;; 1.1d-322 is closest to denormal #x16 (= 22). + (ext:with-float-traps-masked (:underflow :inexact) + (assert-equal 0 (kernel:double-float-high-bits 1.1d-322)) + (assert-equal #x16 (kernel:double-float-low-bits 1.1d-322)))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/b71a67e281f3eab0f6e950b... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/b71a67e281f3eab0f6e950b... You're receiving this email because of your account on gitlab.common-lisp.net. Manage all notifications: https://gitlab.common-lisp.net/-/profile/notifications | Help: https://gitlab.common-lisp.net/help