6
\$\begingroup\$

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.

\$\endgroup\$

1 Answer 1

2
\$\begingroup\$

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?

\$\endgroup\$
3
  • \$\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
    Commented 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
    Commented 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
    Commented Apr 14, 2011 at 20:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.