Raymond Toy pushed to branch issue-504-read-denormals-with-rounding 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. - - - - - 1 changed file: - src/code/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))))))) View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/50d9df5cd5c68fd9760d5f21... -- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/-/commit/50d9df5cd5c68fd9760d5f21... 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