I'm trying to speed up the following code for summing doubles with full precision, based on a Python recipe from Raymond Hettinger. Any thoughts on how to make this go faster?
-- Full precision summation based on http://code.activestate.com/recipes/393090/
msum :: [Double] -> Double
msum [] = 0
msum (x:xs) = sum $ foldl' (\ps x' -> reverse (partials x' ps)) [x] xs
partials :: Double -> [Double] -> [Double]
partials x = foldl' (\acc p -> case hilo x p of (hi, 0) -> hi:acc
(hi, lo) -> hi:lo:acc) []
hilo :: Double -> Double -> (Double, Double)
hilo x y | abs x < abs y = hilo y x
| otherwise = (hi, lo) where hi = x + y
lo = y - (hi - x)
This is already a second draft; you can see a previous version of the code in this gist.
1 Answer 1
You could write hilo
as hilo :: Double -> Double -> [Double]
, which returns [hi,lo]
for lo /= 0
and [hi]
for lo = 0
. Then you can write partials
as
partials x = foldl' (\acc p -> (hilo x p) ++ acc) -> hi:acc []
BTW: How could lo
ever become non-zero considering
lo = y - (hi - x) = y - (x + y - x) = y - y = 0
Rounding?
-
\$\begingroup\$ Did you benchmark that approach? For me it's ending up slightly slower. And, yes, lo only ends up nonzero as an artifact of IEEE754; the lack of precision in floating point is what this function is intended to address. The original Python recipe I linked to explains more and links to the paper that proves it works. \$\endgroup\$jfager– jfager2011年04月14日 16:30:49 +00:00Commented Apr 14, 2011 at 16:30
-
\$\begingroup\$ Also, I've got a version that passes the acc through to hilo and lets hilo do the appending itself. I'm really surprised to see that that version is also slightly slower than this in wall time on my benchmark, and ends up not doing the math right when run under the profiler (and before ghc 7.0.3, it didn't do the math right at all - I need to file a bug). \$\endgroup\$jfager– jfager2011年04月14日 16:38:32 +00:00Commented Apr 14, 2011 at 16:38
-
\$\begingroup\$ It was just an educated guess, no benchmarking. I don't see much other things you could try. Maybe you could get rid of the reverse in msum by building the result list in partials somehow backwards? \$\endgroup\$Landei– Landei2011年04月14日 20:29:50 +00:00Commented Apr 14, 2011 at 20:29