[Haskell-cafe] How to improve speed? (MersenneTwister is
several times slower than C version)
isto
isto.aho at dnainternet.net
Thu Nov 2 14:09:21 EST 2006
Hi & no problems (I didn't tell it clearly right away),
I modified the code along the comments given by Lemmih and things
improved a lot. mod-operator can be removed by two loops as in C
version, which still further improved the speed. I tried this with
the old version and the speed-up was next to nothing. The new
version can be found below.
Now, the results are comparable to the C version. On the second run
it was even better, but usually it is about 0.78 on my machine.
$ time ./testMT
Testing Mersenne Twister.
3063349438
real 0m0.791s
user 0m0.776s
sys 0m0.008s
$ time ./testMT
Testing Mersenne Twister.
3063349438
real 0m0.414s
user 0m0.400s
Can somebody tell, why this happens? There were similar timings
with ST-version and with the old version (so that it sometimes runs
almost twice as fast or twice as slow than normally). Is this something
system dependent (amd64/ubuntu)? C version seems to be within 10%
that is between 0.60 and 0.66.
Anyhow, this random generator now seems to randomly outperform this
particular C version :) C version is located at
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
About Lemmih's proposals. Did I follow them correctly? I'm using ghc
and there the inlining demanding seemed to work. With =100 there are
more often 0.4 timings than with =16, but also =100 gives 0.78 seconds
often (majority of runs).
And can something still be done? At this point I would leave FFI and
other languages out as one goal here is to learn Haskelling and somehow
solutions that require writing code only in one language sounds nicer.
br, Isto
to, 2006-11-02 kello 08:18 -0500, Lennart Augustsson kirjoitti:
> Oh, sorry, I thought your version was a rewritten version of mine. :)
> The names are so similar, after all.
-------------------------- Mersenne.hs (64bit version removed)
module Mersenne where
import Data.Bits
import Data.Word
import Data.Array.Base
import Data.Array.MArray
import Data.Array.IO
data MT32 = MT32 (IOUArray Int Word32) Int
last32bitsof :: Word32 -> Word32
last32bitsof a = a .&. 0xffffffff -- == (2^32-1)
lm32 = 0x7fffffff :: Word32
um32 = 0x80000000 :: Word32
mA32 = 0x9908b0df :: Word32 -- == 2567483615
-- Array of length 624.
initialiseGenerator32 :: Int -> IO MT32
initialiseGenerator32 seed = do
let s = last32bitsof (fromIntegral seed)::Word32
mt <- newArray (0,623) (0::Word32)
unsafeWrite mt 0 s
mtLoop mt s 1
generateNumbers32 mt
return (MT32 mt 0)
mtLoop :: (IOUArray Int Word32) -> Word32 -> Int -> IO (IOUArray Int
Word32)
mtLoop mt lastNro n = loop lastNro n
where
loop :: Word32 -> Int -> IO (IOUArray Int Word32)
loop lastNro n
| n == 624 = return mt
| otherwise = do
let n1 = lastNro `xor` (shiftR lastNro 30)
new = (1812433253 * n1 + (fromIntegral n)::Word32)
unsafeWrite mt n new
loop new $! (n+1)
generateNumbers32 :: (IOUArray Int Word32) -> IO ()
generateNumbers32 mt = do
gLoop1 0
gLoop2 227
wL <- unsafeRead mt 623
w0 <- unsafeRead mt 0
w396 <- unsafeRead mt 396
let y = (wL .&. um32) .|. (w0 .&. lm32) :: Word32
if even y
then unsafeWrite mt 623 (w396 `xor` (shiftR y 1))
else unsafeWrite mt 623 (w396 `xor` (shiftR y 1) `xor` mA32)
return ()
where
gLoop1 :: Int -> IO ()
gLoop1 227 = return ()
gLoop1 i = do
wi <- unsafeRead mt i
wi1 <- unsafeRead mt (i+1)
-- w3 <- unsafeRead mt ((i+397) `mod` 624)
w3 <- unsafeRead mt (i+397)
let y = (wi .&. um32) .|. (wi1 .&. lm32)
if even y
then unsafeWrite mt i (w3 `xor` (shiftR y 1))
else unsafeWrite mt i (w3 `xor` (shiftR y 1) `xor` mA32)
gLoop1 $! (i+1)
gLoop2 :: Int -> IO ()
gLoop2 623 = return ()
gLoop2 i = do
wi <- unsafeRead mt i
wi1 <- unsafeRead mt (i+1)
-- w3 <- unsafeRead mt ((i+397) `mod` 624)
w3 <- unsafeRead mt (i-227)
let y = (wi .&. um32) .|. (wi1 .&. lm32)
if even y
then unsafeWrite mt i (w3 `xor` (shiftR y 1))
else unsafeWrite mt i (w3 `xor` (shiftR y 1) `xor` mA32)
gLoop2 $! (i+1)
{-
gLoop :: Int -> IO ()
gLoop 623 = do
wL <- unsafeRead mt 623
w0 <- unsafeRead mt 0
w396 <- unsafeRead mt 396
let y = (wL .&. um32) .|. (w0 .&. lm32) :: Word32
if even y
then unsafeWrite mt 623 (w396 `xor` (shiftR y 1))
else unsafeWrite mt 623 (w396 `xor` (shiftR y 1) `xor` mA32)
return ()
gLoop i = do
wi <- unsafeRead mt i
wi1 <- unsafeRead mt (i+1)
w3 <- unsafeRead mt ((i+397) `mod` 624)
let y = (wi .&. um32) .|. (wi1 .&. lm32)
if even y
then unsafeWrite mt i (w3 `xor` (shiftR y 1))
else unsafeWrite mt i (w3 `xor` (shiftR y 1) `xor` mA32)
gLoop $! (i+1)
-}
-- next32 :: MT32 -> IO (Word32, MT32)
-- next32 (MT32 mt 624) = do
-- next32 (MT32 mt i) = do
next32 :: IOUArray Int Word32 -> Int -> IO (Word32, Int)
next32 mt 624 = do
generateNumbers32 mt
-- let m = MT32 mt 0
-- (w,m') <- next32 m
-- return (w,m')
(w,iN) <- next32 mt 0
return (w,iN)
next32 mt i = do
y <- unsafeRead mt i
let y1 = y `xor` (shiftR y 11)
y2 = y1 `xor` ((shiftL y1 7 ) .&. 0x9d2c5680) -- == 2636928640
y3 = y2 `xor` ((shiftL y2 15) .&. 0xefc60000) -- == 4022730752
y4 = y3 `xor` (shiftR y3 18)
return $ (y4, (i+1))
-- return $ (y4, MT32 mt (i+1))
--------------------------
-------------------------- testMT.hs
module Main where
-- Compile eg with
-- ghc -O3 -optc-O3 -optc-ffast-math -fexcess-precision
-- -funfolding-use-threshold=16 --make testMT
import Mersenne
import GHC.Prim
-- genRNums32 :: MT32 -> Int -> IO Int
-- genRNums32 mt nCnt = gRN mt nCnt
genRNums32 :: MT32 -> Int -> IO MT32
genRNums32 (MT32 mt pos) nCnt = gRN nCnt pos
where
gRN :: Int -> Int -> IO MT32
gRN 1 iCurr = do
(r,iNew) <- next32 mt iCurr
putStrLn $ (show r)
return (MT32 mt iNew)
gRN nCnt iCurr = do
(_,iNew) <- next32 mt iCurr
gRN (nCnt-1) $! iNew
main = do
putStrLn "Testing Mersenne Twister."
mt32 <- initialiseGenerator32 100
genRNums32 mt32 10000000
--------------------------
More information about the Haskell-Cafe
mailing list