math.integers, comment and simplify bignum/f

change the "while" that could only execute once to "when"
change the f/loop word name to "mantissa-and-guard" since it's what it
computes
change the check against 2^53 to be explicit

core/math/integers/integers-tests.factor

@@ -279,3 +279,6 @@ IN: math.integers.tests
 { 0x0.8p-1022 } [ 8 1026 2^ /f ] unit-test
 { 0x0.6p-1022 } [ 6 1026 2^ /f ] unit-test
 { 0x0.4p-1022 } [ 4 1026 2^ /f ] unit-test
+
+! rounding triggering special case in post-scale
+{ 1.0 } [ 300 2^ 1 - 300 2^ /f ] unit-test

core/math/integers/integers.factor

@@ -116,6 +116,16 @@ M: bignum (log2) bignum-log2 ; inline
 ! provided with absolutely no warranty."
 
 ! First step: pre-scaling
+! As an optimization to minimize the size of the operands of the bignum
+! divisions we will do, we start by stripping any trailing zeros from
+! the denominator and move it into the scale factor.
+! We want a result in ]2^54;2^53] to find the mantissa
+! in ]2^53,2^52] with 1 extra "guard" bit for rounding.
+! So we shift the numerator to get the result of the integer division
+! "num/den" in the range ]2^54; 2^53]; Our shift is only a guess
+! based on the magnitude of the inputs, so it
+! will actually give results in the range ]2^55; 2^53].
+! Note: epsilon is used for rounding in step 3.
 : twos ( x -- y ) dup 1 - bitxor log2 ; inline
 
 : scale-denonimator ( den -- scaled-den scale' )
@@ -129,33 +139,37 @@ M: bignum (log2) bignum-log2 ; inline
         54 + [ (epsilon?) ] [ shift ] 2bi
     ] keep ; inline
 
-: pre-scale ( num den -- epsilon? mantissa den' scale )
+: pre-scale ( num den -- epsilon? num' den' scale )
     scale-denonimator [
         [ scale-numerator ] keep swap
     ] dip swap - ; inline
 
-! Second step: loop
+! Second step: compute mantissa
+! "num/den" would be in the range ]2^55; 2^53]. After this step
+! it will be in the range ]2^54; 2^53]. Compute "num/den" and the
+! reminder used for rounding
 : (2/-with-epsilon) ( epsilon? num -- epsilon?' num' )
     [ 1 bitand zero? not or ] [ 2/ ] bi ; inline
 
-: /f-loop ( epsilon? mantissa den scale -- epsilon?' fraction-and-guard rem scale' )
+: mantissa-and-guard ( epsilon? num den scale -- epsilon?' mantissa-and-guard rem scale' )
-    [ 2over /i log2 53 > ]
+    2over /i log2 53 >
-    [ [ (2/-with-epsilon) ] [ ] [ 1 + ] tri* ] while
+    [ [ (2/-with-epsilon) ] [ ] [ 1 + ] tri* ] when
     [ /mod ] dip ; inline
 
-! Third step: post-scaling
+! Third step: rounding
-: scale-float ( mantissa scale -- float' )
+!
-    {
+! if the guard bit is 0, round down
-        { [ dup 1024 > ] [ 2drop 1/0. ] }
+! else if the guard bit is 1 and (rem != 0 or epsilon is true), round up
-        { [ dup -1021 < ] [ 1021 + shift bits>double ] }
+! else break the tie by alternating rounding down or up to avoid accumulating errors
-        [ [ 52 2^ 1 - bitand ] dip 1022 + 52 shift bitor bits>double ]
+!
-    } cond ; inline
+! The epsilon trick works because epsilon is true if numerator bits were discarded.
-
+! Mathematically, (num+epsilon)/denom = (num/denum) + (epsilon/denom)
-: post-scale ( mantissa scale -- n )
+! We have actually computed the "num/denum" part and use the "epsilon/denom"
-    [ 2/ ] dip over log2 52 > [ [ 2/ ] [ 1 + ] bi* ] when
+! to choose the correct rounding.
-    scale-float ; inline
+!
-
+! Note that rounding down means doing nothing because we will
-: round-to-nearest ( epsilon? fraction-and-guard rem -- fraction-and-guard' )
+! discard the guard bit after this
+: round-to-nearest ( epsilon? mantissa-and-guard rem -- mantissa-and-guard' )
     over odd?
     [
         zero? [
@@ -164,12 +178,28 @@ M: bignum (log2) bignum-log2 ; inline
     ] [ drop nip ] if ;
     inline
 
+! Fourth step: post-scaling
+! Because of rounding, our mantissa with guard bit is now in the
+! range [2^54;2^53], so we have to handle 2^54 specially.
+: scale-float ( mantissa scale -- float' )
+    ! At this point, the scale value is the exponent minus 1.
+    {
+        { [ dup 1024 > ] [ 2drop 1/0. ] }
+        { [ dup -1021 < ] [ 1021 + shift bits>double ] } ! subnormals and underflow
+        [ [ 52 2^ 1 - bitand ] dip 1022 + 52 shift bitor bits>double ]
+    } cond ; inline
+
+: post-scale ( mantissa scale -- n )
+    [ 2/ ] dip ! drop guard bit
+    over 53 2^ = [ [ 2/ ] [ 1 + ] bi* ] when
+    scale-float ; inline
+
 ! Main word
 : /f-abs ( m n -- f )
     over zero? [ nip zero? 0/0. 0.0 ? ] [
         [ drop 1/0. ] [
             pre-scale
-            /f-loop
+            mantissa-and-guard
             [ round-to-nearest ] dip
             post-scale
         ] if-zero