Fix typo in floatToDigits

[ghc-base.git] / GHC / Float.lhs

diff --git a/GHC/Float.lhs b/GHC/Float.lhs
index bb9aa20..e0cc415 100644 (file)
--- a/GHC/Float.lhs
+++ b/GHC/Float.lhs
@@ -595,30 +595,46 @@ floatToDigits base x =
   -- will have an impossibly low exponent.  Adjust for this.
   (f, e) =
    let n = minExp - e0 in
-   if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
+   if n > 0 then (f0 `quot` (expt b n), e0+n) else (f0, e0)
   (r, s, mUp, mDn) =
    if e >= 0 then
-    let be = b^ e in
-    if f == b^(p-1) then
-      (f*be*b*2, 2*b, be*b, b)
+    let be = expt b e in
+    if f == expt b (p-1) then
+      (f*be*b*2, 2*b, be*b, be)     -- according to Burger and Dybvig
     else
       (f*be*2, 2, be, be)
    else
-    if e > minExp && f == b^(p-1) then
-      (f*b*2, b^(-e+1)*2, b, 1)
+    if e > minExp && f == expt b (p-1) then
+      (f*b*2, expt b (-e+1)*2, b, 1)
     else
-      (f*2, b^(-e)*2, 1, 1)
+      (f*2, expt b (-e)*2, 1, 1)
   k :: Int
   k =
    let
     k0 :: Int
     k0 =
      if b == 2 && base == 10 then
-        -- logBase 10 2 is slightly bigger than 3/10 so
-        -- the following will err on the low side.  Ignoring
-        -- the fraction will make it err even more.
-        -- Haskell promises that p-1 <= logBase b f < p.
-        (p - 1 + e0) * 3 `div` 10
+        -- logBase 10 2 is very slightly larger than 8651/28738
+        -- (about 5.3558e-10), so if log x >= 0, the approximation
+        -- k1 is too small, hence we add one and need one fixup step less.
+        -- If log x < 0, the approximation errs rather on the high side.
+        -- That is usually more than compensated for by ignoring the
+        -- fractional part of logBase 2 x, but when x is a power of 1/2
+        -- or slightly larger and the exponent is a multiple of the
+        -- denominator of the rational approximation to logBase 10 2,
+        -- k1 is larger than logBase 10 x. If k1 > 1 + logBase 10 x,
+        -- we get a leading zero-digit we don't want.
+        -- With the approximation 3/10, this happened for
+        -- 0.5^1030, 0.5^1040, ..., 0.5^1070 and values close above.
+        -- The approximation 8651/28738 guarantees k1 < 1 + logBase 10 x
+        -- for IEEE-ish floating point types with exponent fields
+        -- <= 17 bits and mantissae of several thousand bits, earlier
+        -- convergents to logBase 10 2 would fail for long double.
+        -- Using quot instead of div is a little faster and requires
+        -- fewer fixup steps for negative lx.
+        let lx = p - 1 + e0
+            k1 = (lx * 8651) `quot` 28738
+        in if lx >= 0 then k1 + 1 else k1
      else
        -- f :: Integer, log :: Float -> Float,
         --               ceiling :: Float -> Int
@@ -637,7 +653,7 @@ floatToDigits base x =
 
   gen ds rn sN mUpN mDnN =
    let
-    (dn, rn') = (rn * base) `divMod` sN
+    (dn, rn') = (rn * base) `quotRem` sN
     mUpN' = mUpN * base
     mDnN' = mDnN * base
    in