Haskell repa : inner product result of two array is wrong

Question

import Data.Array.Repa
import Data.Vector.Unboxed hiding(zipWith)
import Prelude hiding(zipWith,replicate)

dotp :: Array U DIM1 Float -> Array U DIM1 Float -> IO Float
dotp x y =  sumAllP $ zipWith (*) x y

x = fromUnboxed (Z:.20000000) $ replicate 20000000 1 
y = fromUnboxed (Z:.20000000) $ replicate 20000000 1

main = (dotp x y) >>= print

Compiled with ghc -O2 -threaded test.hs

However, when I execute the executable with ./test +RTS -N1

the result is 1.6777216e7

while ./test +RTS -N2 or with more cores,

the result is correct : 2.0e7

What is wrong with my code?

Answer 1

You use Float instead of Double. Float doesn't have enough significant digits (~7) for your use case:

λ> (1.6777216e7 :: Float) + 1
1.6777216e7

Switch to Double. If all your values are integral, consider Int (or Intxx from Data.Int) or Word (from Data.Word).

See also:

• WP: Single-precision floating-point format
• What Every Computer Scientist Should Know About Floating-Point Arithmetic