[Haskell-cafe] How to improve speed? (MersenneTwister is several times slower than C version)

[isto]

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

--------------------------