Improve numeric stability of numericEnumFrom for floating numbers

Fixes #15081.

Test Plan: cd libraries/base && make test TEST="enumNumeric"

Reviewers: hvr, bgamari

Reviewed By: bgamari

Subscribers: dfeuer, simonpj, thomie, carter

GHC Trac Issues: #15081

Differential Revision: https://phabricator.haskell.org/D4650

libraries/base/GHC/Real.hs

@@ -216,10 +216,19 @@ class  (Real a, Fractional a) => RealFrac a  where
 -- These 'numeric' enumerations come straight from the Report
 
 numericEnumFrom         :: (Fractional a) => a -> [a]
-numericEnumFrom n       =  n `seq` (n : numericEnumFrom (n + 1))
+numericEnumFrom n       = go 0
+  where
+    -- See Note [Numeric Stability of Enumerating Floating Numbers]
+    go !k = let !n' = n + k
+             in n' : go (k + 1)
 
 numericEnumFromThen     :: (Fractional a) => a -> a -> [a]
-numericEnumFromThen n m = n `seq` m `seq` (n : numericEnumFromThen m (m+m-n))
+numericEnumFromThen n m = go 0
+  where
+    step = m - n
+    -- See Note [Numeric Stability of Enumerating Floating Numbers]
+    go !k = let !n' = n + k * step
+             in n' : go (k + 1)
 
 numericEnumFromTo       :: (Ord a, Fractional a) => a -> a -> [a]
 numericEnumFromTo n m   = takeWhile (<= m + 1/2) (numericEnumFrom n)
@@ -232,6 +241,49 @@ numericEnumFromThenTo e1 e2 e3
                                  predicate | e2 >= e1  = (<= e3 + mid)
                                            | otherwise = (>= e3 + mid)
 
+{- Note [Numeric Stability of Enumerating Floating Numbers]
+-----------------------------------------------------------
+When enumerate floating numbers, we could add the increment to the last number
+at every run (as what we did previously):
+
+    numericEnumFrom n =  n `seq` (n : numericEnumFrom (n + 1))
+
+This approach is concise and really fast, only needs an addition operation.
+However when a floating number is large enough, for `n`, `n` and `n+1` will
+have the same binary representation. For example (all number has type
+`Double`):
+
+    9007199254740990                 is: 0x433ffffffffffffe
+    9007199254740990 + 1             is: 0x433fffffffffffff
+    (9007199254740990 + 1) + 1       is: 0x4340000000000000
+    ((9007199254740990 + 1) + 1) + 1 is: 0x4340000000000000
+
+When we evaluate ([9007199254740990..9007199254740991] :: Double), we would
+never reach the condition in `numericEnumFromTo`
+
+    9007199254740990 + 1 + 1 + ... > 9007199254740991 + 1/2
+
+We would fall into infinite loop (as reported in Trac #15081).
+
+To remedy the situation, we record the number of `1` that needed to be added
+to the start number, rather than increasing `1` at every time. This approach
+can improvement the numeric stability greatly at the cost of a multiplication.
+
+Furthermore, we use the type of the enumerated number, `Fractional a => a`,
+as the type of multiplier. In rare situations, the multiplier could be very
+large and will lead to the enumeration to infinite loop, too, which should
+be very rare. Consider the following example:
+
+    [1..9007199254740994]
+
+We could fix that by using an Integer as multiplier but we don't do that.
+The benchmark on T7954.hs shows that this approach leads to significant
+degeneration on performance (33% increase allocation and 300% increase on
+elapsed time).
+
+See Trac #15081 and Phab:D4650 for the related discussion about this problem.
+-}
+
 --------------------------------------------------------------
 -- Instances for Int
 --------------------------------------------------------------

libraries/base/tests/all.T

@@ -2,6 +2,7 @@
 test('readFloat', exit_code(1), compile_and_run, [''])
 test('enumDouble', normal, compile_and_run, [''])
 test('enumRatio', normal, compile_and_run, [''])
+test('enumNumeric', normal, compile_and_run, [''])
 test('tempfiles', normal, compile_and_run, [''])
 test('fixed', normal, compile_and_run, [''])
 test('quotOverflow', normal, compile_and_run, [''])

libraries/base/tests/enumNumeric.hs

+main :: IO ()
+main = do
+    print $ map (/2) ([5..6] :: [Double])
+    print $ ([9007199254740990..9007199254740991] :: [Double])
+    print $ map (/2) ([9007199254740990..9007199254740991] :: [Double])
+    print $ ([9007199254740989..9007199254740990] :: [Double])
+    print $ map (/2) ([9007199254740989..9007199254740990] :: [Double])

libraries/base/tests/enumNumeric.stdout

+[2.5,3.0]
+[9.00719925474099e15,9.007199254740991e15,9.007199254740992e15,9.007199254740992e15]
+[4.503599627370495e15,4.5035996273704955e15,4.503599627370496e15,4.503599627370496e15]
+[9.007199254740989e15,9.00719925474099e15]
+[4.5035996273704945e15,4.503599627370495e15]

testsuite/tests/perf/should_run/all.T

@@ -384,8 +384,9 @@ test('T7954',
                       [(wordsize(32), 920045264, 10),
               # some date:  1380051408    (64-bit Windows machine)
               # 2014-04-04:  920045264    (64-bit Windows machine)
-                       (wordsize(64), 1680051336, 10)]),
+                       (wordsize(64), 1280051632, 10)]),
               # 2014-02-10: 1680051336 (x86_64/Linux), call arity analysis
+              # 2018-05-03: 1280051632 (x86_64/Linux), refactor numericEnumFrom
       only_ways(['normal'])
       ],
      compile_and_run,