Good catch. Adding and subtracting over and over relying on massive
cancellation over and over is a recipe for floating point disaster!
That said, so is adding small things to a large thing over and over, you'll
accumulate the rounding error from the addition of small numbers to the
large base.
An even better fix would be then to track the base and the current
accumulated total delta from the base and do a final addition at the end.
Then you only ever add like sized step sizes together before adding them to
the base. This would stage the additions such that you get another matissa
sized window of possible accumulation.
Something like:
enumFromThen n m = go (n - m) n 0 where
go !d !b !a = db `seq` db : go d b (a + d) where
db = d + b
which probably beats Kahan in practice, Kahan-Babuška-Neumaier should be
more stable still, and there are other techniques that go further into
accuracy at the cost of significant performance.
But we don't need a general number summation algorithm. All the numbers
except the base are the same, we have a final hammer available to us: just
do multiplication of the delta by the number of steps and add it to the
base.
That should be the most numerically stable thing possible, given that we
are forced to do at least the first massive cancellation between the from
and then steps by the API we have to meet, but I don't have benchmarks to
say which technique wins in practice.
-Edward
On Tue, Aug 9, 2016 at 11:22 PM, Andrew Farmer
Noticed this today:
ghci> let xs = [0.0,0.1 .. 86400.0] in maximum xs 86400.0000005062
enumFromThenTo is implemented by numericEnumFromThenTo:
https://github.com/ghc/ghc/blob/a90085bd45239fffd65c01c24752a9 bbcef346f1/libraries/base/GHC/Real.hs#L227
Which probably accumulates error in numericEnumFromThen with the (m+m-n):
numericEnumFromThen n m = n `seq` m `seq` (n : numericEnumFromThen m (m+m-n))
Why not define numericEnumFromThen as:
numericEnumFromThen n m = let d = m - n in d `seq` go d n where go delta x = x `seq` (x : go delta (x + delta))
(or with BangPatterns)
numericEnumFromThen n m = go (m - n) n where go !delta !x = x : go delta (x + delta)
Seems like we'd save a lot of subtractions by using the worker function. _______________________________________________ ghc-devs mailing list ghc-devs@haskell.org http://mail.haskell.org/cgi-bin/mailman/listinfo/ghc-devs