enumFromThenTo for Doubles
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/a90085bd45239fffd65c01c24752a9bbcef346f1/lib... 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.
Turns out the accumulated error is even worse:
Prelude> let old x y z = let eftt i j = i : eftt j (j+j-i) in let d =
y - x in maximum $ takeWhile (<= z + d) $ eftt x y
Prelude> old 0.0 0.1 86400.0
86400.0000005062
Prelude> let new x y z = let d = y - x in let go i = i : go (i + d) in
maximum $ takeWhile (<= z + d) $ go x
Prelude> new 0.0 0.1 86400.0
86400.00000054126
Sorry to spam the list. :-P Floating point is hard.
On Tue, Aug 9, 2016 at 8: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/a90085bd45239fffd65c01c24752a9bbcef346f1/lib...
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.
Sounds somewhat plausible. By all means give it a try. S -----Original Message----- From: ghc-devs [mailto:ghc-devs-bounces@haskell.org] On Behalf Of Andrew Farmer Sent: 10 August 2016 04:22 To: ghc-devs@haskell.org Subject: enumFromThenTo for Doubles 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/a90085bd45239fffd65c01c24752a9bbcef346f1/lib... 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 https://na01.safelinks.protection.outlook.com/?url=http%3a%2f%2fmail.haskell.org%2fcgi-bin%2fmailman%2flistinfo%2fghc-devs&data=02%7c01%7csimonpj%40microsoft.com%7c65f112900fb44b7186a408d3c0cd8631%7c72f988bf86f141af91ab2d7cd011db47%7c1%7c0%7c636063961357184976&sdata=Gz0DQ%2fGEUIyfHtmAjdbxpBt3YEnxbpoKKiygnCb%2fhYo%3d
Way back when I started with haskell I noticed this, and switched to using
this:
-- | Enumerate an inclusive range. Uses multiplication instead of
successive
-- addition to avoid loss of precision.
--
-- Also it doesn't require an Enum instance.
range :: (Num a, Ord a) => a -> a -> a -> [a]
range start end step = go 0
where
go i
| step >= 0 && val > end = []
| step < 0 && val < end = []
| otherwise = val : go (i+1)
where val = start + (i*step)
It's always seemed better in every way, except syntax convenience.
Wouldn't any approach with successive addition lose precision?
On Tue, Aug 9, 2016 at 8: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
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
participants (4)
-
Andrew Farmer -
Edward Kmett -
Evan Laforge -
Simon Peyton Jones