GitLab

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.

@@ -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)))))))

Raymond Toy pushed to branch issue-504-read-denormals-with-rounding at cmucl / cmucl

Commits:

1 changed file:

Changes: