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

diff --git a/core/math/integers/integers-tests.factor b/core/math/integers/integers-tests.factor
index e83aeffb9ec79a42ea1ba06d99031c8aaf2cf456..ddb6c908db2de6c9658db304db41ed2514e940fe 100644 (file)
--- a/core/math/integers/integers-tests.factor
+++ b/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

diff --git a/core/math/integers/integers.factor b/core/math/integers/integers.factor
index e675db2f6865dd543748389bcc7112267d428e2b..07a2e2bbd2e53a9944f4e4f0c11260fa973bdee4 100644 (file)
--- a/core/math/integers/integers.factor
+++ b/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,47 +139,67 @@ 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' )
-    [ 2over /i log2 53 > ]
-    [ [ (2/-with-epsilon) ] [ ] [ 1 + ] tri* ] while
+: mantissa-and-guard ( epsilon? num den scale -- epsilon?' mantissa-and-guard rem scale' )
+    2over /i log2 53 >
+    [ [ (2/-with-epsilon) ] [ ] [ 1 + ] tri* ] when
     [ /mod ] dip ; inline
 
-! Third step: post-scaling
+! Third step: rounding
+!
+! if the guard bit is 0, round down
+! else if the guard bit is 1 and (rem != 0 or epsilon is true), round up
+! else break the tie by alternating rounding down or up to avoid accumulating errors
+!
+! The epsilon trick works because epsilon is true if numerator bits were discarded.
+! Mathematically, (num+epsilon)/denom = (num/denum) + (epsilon/denom)
+! We have actually computed the "num/denum" part and use the "epsilon/denom"
+! to choose the correct rounding.
+!
+! Note that rounding down means doing nothing because we will
+! discard the guard bit after this
+: round-to-nearest ( epsilon? mantissa-and-guard rem -- mantissa-and-guard' )
+    over odd?
+    [
+        zero? [
+            dup 2 bitand zero? not rot or [ 1 + ] when
+        ] [ nip 1 + ] if
+    ] [ 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 ] }
+        { [ 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 over log2 52 > [ [ 2/ ] [ 1 + ] bi* ] when
+    [ 2/ ] dip ! drop guard bit
+    over 53 2^ = [ [ 2/ ] [ 1 + ] bi* ] when
     scale-float ; inline
 
-: round-to-nearest ( epsilon? fraction-and-guard rem -- fraction-and-guard' )
-    over odd?
-    [
-        zero? [
-            dup 2 bitand zero? not rot or [ 1 + ] when
-        ] [ nip 1 + ] if
-    ] [ drop nip ] if ;
-    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