Raymond Toy pushed to branch rtoy-print-using-ryu at cmucl / cmucl Commits: 50d9df5c by Raymond Toy at 2026-05-26T18:47:45-07:00 Fix double rounding when reading denormal float literals `float-ratio-float` was always rounding the rational full precision before handing the result to `scale-float`. When the result was a denormal, `scale-float-maybe-underflow` could then round the number again to remove some bits to fit in a denormal. This causes double rounding as seen with "7.290983e-39". To fix this `float-ratio-float` re-rounds directly to denormal range. The number of extra bits in the denormal is detected and the extra bits are dropped using round-to-nearest-even. The denormal is emitted directly so `scale-float` isn't involved. The normal path is unchanged and `scale-float-maybe-underflow` still needs correctly rounding in case `scale-float` would produce a denormal. - - - - - bc6c560e by Raymond Toy at 2026-05-26T18:56:29-07:00 Add tests to exercise float-ratio-float Each test calls `float-ratio-float` with an exact rational to test a specific branch of the new helpers. - - - - - 94858a23 by Raymond Toy at 2026-05-26T20:41:16-07:00 Merge branch 'issue-504-read-denormals-with-rounding' into rtoy-print-using-ryu - - - - - 2 changed files: - src/code/float.lisp - tests/float.lisp Changes: ===================================== src/code/float.lisp ===================================== @@ -1152,26 +1152,74 @@ (assert (= len (the fixnum (1+ digits)))) (multiple-value-bind (f0) (floatit (ash bits -1)) - #+nil - (progn - (format t "x = ~A~%" x) - (format t "1: f0, f1 = ~A~%" f0) - (format t " scale = ~A~%" (1+ scale))) - (scale-float f0 (1+ scale)))) (t (multiple-value-bind (f0) (floatit bits) - #+nil - (progn - (format t "2: f0, f1 = ~A~%" f0) - (format t " scale = ~A~%" scale) - (format t "scale-float f0 = ~A~%" (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)))))) + (denormal-excess () + ;; How many bits of precision the result loses by being + ;; denormal instead of normal. A normal-precision return + ;; would be BITS*2^(SCALE-DIGITS) with BITS having DIGITS + ;; bits. Once it's known that this representation will + ;; produce a denormal -- equivalently, that SCALE-FLOAT- + ;; MAYBE-UNDERFLOW would take the underflow branch -- + ;; (1 - BIAS) - SCALE bits of the mantissa fall below the + ;; denormal's narrower storage and must be rounded off. + ;; Zero in the normal range. + (let ((bias + (ecase format + (single-float vm:single-float-bias) + (double-float vm:double-float-bias)))) + (declare (fixnum bias)) + (max 0 (the fixnum + (- (the fixnum (- 1 bias)) + scale))))) + (round-denormal (fraction-and-guard rem excess) + ;; FRACTION-AND-GUARD has (1+ DIGITS) bits with one guard + ;; bit; round it to (- DIGITS EXCESS) bits using round-to- + ;; nearest, ties to even, with REM as the sticky tail. + ;; Drops EXCESS+1 low bits in a single step. This is the + ;; one rounding the denormal result undergoes; no + ;; subsequent SCALE-FLOAT call is needed, so there is no + ;; double rounding. + (declare (type unsigned-byte fraction-and-guard rem) + (fixnum excess)) + (let* ((shift (1+ excess)) + (low (ldb (byte shift 0) fraction-and-guard)) + (quot (ash fraction-and-guard (- shift))) + (halfway (ash 1 excess))) + (declare (fixnum shift)) + (cond ((< low halfway) quot) + ((> low halfway) (1+ quot)) + ((not (zerop rem)) (1+ quot)) + ((oddp quot) (1+ quot)) + (t quot)))) + (denormal-from-bits (mantissa excess) + ;; MANTISSA has at most (- DIGITS EXCESS) bits and is the + ;; stored significand of a denormal result. Denormal + ;; storage holds (1- DIGITS) bits, so rounding can carry + ;; into the smallest normal only when EXCESS = 1, in + ;; which case MANTISSA can be exactly (ASH 1 (1- DIGITS)). + (declare (fixnum excess)) + (let ((sign (if plusp 0 1))) + (case format + (single-float + (cond ((and (= excess 1) + (= mantissa + (ash 1 (1- vm:single-float-digits)))) + (single-from-bits + sign vm:single-float-normal-exponent-min 0)) + (t + (single-from-bits sign 0 mantissa)))) + (double-float + (cond ((and (= excess 1) + (= mantissa + (ash 1 (1- vm:double-float-digits)))) + (double-from-bits + sign vm:double-float-normal-exponent-min 0)) + (t + (double-from-bits sign 0 mantissa))))))) (floatit (bits) (let ((sign (if plusp 0 1))) (case format @@ -1189,16 +1237,36 @@ (declare (fixnum extra)) (cond ((/= extra 1) (assert (> extra 1))) - ((oddp fraction-and-guard) - (return - (if (zerop rem) - (float-and-scale - (if (zerop (logand fraction-and-guard 2)) - fraction-and-guard - (1+ fraction-and-guard))) - (float-and-scale (1+ fraction-and-guard))))) (t - (return (float-and-scale fraction-and-guard))))) + (return + (let ((excess (denormal-excess))) + (cond + ((zerop excess) + ;; Normal result: original odd/even tie-break. + (cond ((oddp fraction-and-guard) + (if (zerop rem) + (float-and-scale + (if (zerop + (logand fraction-and-guard 2)) + fraction-and-guard + (1+ fraction-and-guard))) + (float-and-scale + (1+ fraction-and-guard)))) + (t + (float-and-scale fraction-and-guard)))) + (t + ;; Denormal result: re-round directly to the + ;; denormal's narrower precision so the only + ;; rounding step happens here. Rounding to + ;; DIGITS first and re-rounding via + ;; SCALE-FLOAT-MAYBE-UNDERFLOW would double- + ;; round (e.g. 7.290983e-39 would land on an + ;; artifical tie at the 24-bit boundary). + (let ((mantissa + (round-denormal fraction-and-guard rem + excess))) + (declare (type unsigned-byte mantissa)) + (denormal-from-bits mantissa excess))))))))) (setq shifted-num (ash shifted-num -1)) (incf scale))))))) ===================================== tests/float.lisp ===================================== @@ -853,3 +853,176 @@ (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)))) + +(define-test float-ratio-float.denormal-double-rounding.single + (:tag :issues) + ;; Regression for the double-rounding bug exposed by reading + ;; "7.290983e-39". The exact rational lies at 5203019.332 * 2^-149, + ;; below halfway between denormals #x4f644b and #x4f644c. An earlier + ;; fix rounded to 24 bits first (lifting it to the artificial tie + ;; 5203019.5 * 2^-149) and then re-rounded ties-to-even to #x4f644c. + ;; FLOAT-RATIO-FLOAT now re-rounds directly to denormal precision in a + ;; single step, yielding the correctly-rounded #x4f644b. + (assert-equal #x4f644b + (kernel:single-float-bits 7.290983e-39))) + +(define-test float-ratio-float.denormal-low-below-halfway.single + (:tag :issues) + ;; Exercise ROUND-DENORMAL's `low < halfway' branch via FLOAT-RATIO- + ;; FLOAT. Value 11/10 * least-positive is 1.1 denormal-units; the + ;; loop sees several bits of low and the comparison must take the + ;; round-down path so the mantissa stays at 1. + (let ((x (* 11/10 (expt 2 -149)))) + (assert-equal #x00000001 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x))) + +(define-test float-ratio-float.denormal-low-above-halfway.single + (:tag :issues) + ;; ROUND-DENORMAL's `low > halfway' branch: 17/10 * least-positive + ;; (= 1.7 denormal-units), past halfway between 1 and 2, rounds up. + (let ((x (* 17/10 (expt 2 -149)))) + (assert-equal #x00000002 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x))) + +(define-test float-ratio-float.denormal-tie-to-even.single + (:tag :issues) + ;; ROUND-DENORMAL's tie path, REM = 0: exact halfway between two + ;; denormals rounds to even. 3/2 * least-positive ties between + ;; denormals 1 and 2; the even neighbour is 2. + (let ((x (* 3/2 (expt 2 -149)))) + (assert-equal #x00000002 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; 5/2 * least-positive ties between 2 and 3; even is 2. + (let ((x (* 5/2 (expt 2 -149)))) + (assert-equal #x00000002 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; 7/2 * least-positive ties between 3 and 4; even is 4. + (let ((x (* 7/2 (expt 2 -149)))) + (assert-equal #x00000004 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; 1/2 * least-positive ties between 0 and 1; even is 0. + (let ((x (* 1/2 (expt 2 -149)))) + (assert-equal #x00000000 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x))) + +(define-test float-ratio-float.denormal-tie-with-sticky.single + (:tag :issues) + ;; ROUND-DENORMAL's tie path with REM != 0: when the loop's + ;; FRACTION-AND-GUARD reaches an exact halfway pattern but the + ;; division has a nonzero remainder, the rounding must go up because + ;; the original rational is strictly above halfway. 3/2 * least- + ;; positive + a tiny fraction of least-positive lands just past tie 1-2, + ;; rounds up to 2. + (let ((x (+ (* 3/2 (expt 2 -149)) + (* (expt 2 -149) 1/1000000000)))) + (assert-equal #x00000002 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; 1/2 + tiny: just past tie 0-1, rounds up to 1. + (let ((x (+ (* 1/2 (expt 2 -149)) + (* (expt 2 -149) 1/1000000000)))) + (assert-equal #x00000001 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x))) + +(define-test float-ratio-float.denormal-excess.single + (:tag :issues) + ;; Exercise DENORMAL-EXCESS over a range of magnitudes. EXCESS + ;; varies inversely with the result's magnitude; each case has the + ;; result land on a chosen denormal so an off-by-one in the EXCESS + ;; computation would shift the result by a factor of two. + ;; EXCESS = 1: value near smallest normal; pick (2^23 - 1) * 2^-149, + ;; the largest denormal. + (let ((x (* (1- (expt 2 23)) (expt 2 -149)))) + (assert-equal #x007fffff + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; EXCESS = 4: small denormal, mantissa 2^19. + (let ((x (expt 2 -130))) + (assert-equal (ash 1 19) + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; EXCESS = 21: very small denormal, mantissa 8. + (let ((x (* 8 (expt 2 -149)))) + (assert-equal #x00000008 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; EXCESS = 23 (the maximum): only the bottom three denormals are + ;; reachable. 3 * least-positive gives mantissa 3. + (let ((x (* 3 (expt 2 -149)))) + (assert-equal #x00000003 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x))) + +(define-test float-ratio-float.denormal-carry-to-normal.single + (:tag :issues) + ;; DENORMAL-FROM-BITS's carry-into-smallest-normal branch. Rounding + ;; promotes the denormal mantissa to 2^(DIGITS-1), which doesn't fit + ;; in the denormal's stored mantissa width; the result must be the + ;; smallest normal (stored exponent = 1, mantissa = 0). + ;; + ;; (2^24 - 1)/2 * 2^-149 is an exact tie between the largest denormal + ;; (mantissa 2^23 - 1) and the smallest normal (2^-126); the latter + ;; has stored mantissa 0 which is even, so the tie rounds to it. + (let ((x (* (1- (expt 2 24)) 1/2 (expt 2 -149)))) + (assert-equal #x00800000 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; Just past that tie, also rounds up to smallest normal. + (let ((x (+ (* (1- (expt 2 24)) 1/2 (expt 2 -149)) + (* (expt 2 -149) 1/1000)))) + (assert-equal #x00800000 + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x)) + ;; Just below that tie, rounds down to the largest denormal. + (let ((x (- (* (1- (expt 2 24)) 1/2 (expt 2 -149)) + (* (expt 2 -149) 1/1000)))) + (assert-equal #x007fffff + (kernel:single-float-bits + (kernel::float-ratio-float x 'single-float)) + x))) + +(define-test float-ratio-float.denormal-double-rounding.double + (:tag :issues) + ;; Double-float equivalents. Pick a rational that lands strictly + ;; below halfway between two denormals after 53-bit rounding would + ;; otherwise lift to the halfway position. Verify a few denormal + ;; cases also work end-to-end through FLOAT-RATIO-FLOAT. + ;; + ;; 1.1d-322 -> mantissa #x16 (= 22). + (assert-equal 0 (kernel:double-float-high-bits 1.1d-322)) + (assert-equal #x16 (kernel:double-float-low-bits 1.1d-322)) + ;; Tie-to-even, double: 3/2 * least-positive ties between denormals + ;; 1 and 2; even neighbour is 2. + (let ((x (* 3/2 (expt 2 -1074)))) + (assert-equal 0 (kernel:double-float-high-bits + (kernel::float-ratio-float x 'double-float))) + (assert-equal 2 (kernel:double-float-low-bits + (kernel::float-ratio-float x 'double-float)))) + ;; Carry into smallest normal: (2^53 - 1)/2 * 2^-1074 exactly halfway + ;; between the largest double denormal and the smallest normal; + ;; rounds up to the smallest normal #x00100000_00000000. + (let* ((x (* (1- (expt 2 53)) 1/2 (expt 2 -1074))) + (r (kernel::float-ratio-float x 'double-float))) + (assert-equal #x00100000 (kernel:double-float-high-bits r) x) + (assert-equal 0 (kernel:double-float-low-bits r) x))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/af544ca1ef157e4bc78eddb... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/af544ca1ef157e4bc78eddb... 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