GitLab

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.

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

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

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

Raymond Toy pushed to branch rtoy-print-using-ryu at cmucl / cmucl

Commits:

3 changed files:

Changes: