Permalink
Browse files

Try an immutable MWC generator. Refactor the timing code a little.

  • Loading branch information...
1 parent 74ea0e8 commit 216ddcb7353e05e363ec7a0263b70a9e3103432a @bos committed Sep 20, 2009
Showing with 122 additions and 29 deletions.
  1. +78 −0 tests/MWC.hs
  2. +44 −29 tests/T.hs
View
@@ -0,0 +1,78 @@
+{-# LANGUAGE MagicHash, TypeOperators #-}
+
+module MWC (MWC, init, next) where
+
+import Data.Array.Vector
+--import Data.Bits
+import Data.IntMap
+import Data.Word
+import Prelude hiding (init)
+import GHC.Word
+import GHC.Base (Int(..))
+
+data MWC = MWC {-# UNPACK #-} !Word32 {-# UNPACK #-} !Word8 !(IntMap Word32)
+
+init :: Word32 -> MWC
+init k = MWC 362436 255 (insert 0 k seed)
+{-# INLINE init #-}
+
+-- | Unchecked 64-bit right shift.
+shiftR :: Word64 -> Int -> Word64
+shiftR (W64# x#) (I# i#) = W64# (x# `uncheckedShiftRL64#` i#)
+
+next :: MWC -> Word32 :*: MWC
+next (MWC c i m) = ti :*: MWC (fromIntegral (t `shiftR` 32))
+ j
+ (insert ji ti m)
+ where t = a * fromIntegral (m!ji) + fromIntegral c
+ a = 809430660 :: Word64
+ j = i + 1
+ ji = fromIntegral j
+ ti = fromIntegral t
+{-# INLINE next #-}
+
+seed = fromDistinctAscList $ zip [0..] [
+ 0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184,
+ 0xb3f8f692, 0x4114ea35, 0xb6af0230, 0xebb457d2, 0x47693630, 0x15bc0433,
+ 0x2e1e5b18, 0xbe91129c, 0xcc0815a0, 0xb1260436, 0xd6f605b1, 0xeaadd777,
+ 0x8f59f791, 0xe7149ed9, 0x72d49dd5, 0xd68d9ded, 0xe2a13153, 0x67648eab,
+ 0x48d6a1a1, 0xa69ab6d7, 0x236f34ec, 0x4e717a21, 0x9d07553d, 0x6683a701,
+ 0x19004315, 0x7b6429c5, 0x84964f99, 0x982eb292, 0x3a8be83e, 0xc1df1845,
+ 0x3cf7b527, 0xb66a7d3f, 0xf93f6838, 0x736b1c85, 0x5f0825c1, 0x37e9904b,
+ 0x724cd7b3, 0xfdcb7a46, 0xfdd39f52, 0x715506d5, 0xbd1b6637, 0xadabc0c0,
+ 0x219037fc, 0x9d71b317, 0x3bec717b, 0xd4501d20, 0xd95ea1c9, 0xbe717202,
+ 0xa254bd61, 0xd78a6c5b, 0x043a5b16, 0x0f447a25, 0xf4862a00, 0x48a48b75,
+ 0x1e580143, 0xd5b6a11b, 0x6fb5b0a4, 0x5aaf27f9, 0x668bcd0e, 0x3fdf18fd,
+ 0x8fdcec4a, 0x5255ce87, 0xa1b24dbf, 0x3ee4c2e1, 0x9087eea2, 0xa4131b26,
+ 0x694531a5, 0xa143d867, 0xd9f77c03, 0xf0085918, 0x1e85071c, 0x164d1aba,
+ 0xe61abab5, 0xb8b0c124, 0x84899697, 0xea022359, 0x0cc7fa0c, 0xd6499adf,
+ 0x746da638, 0xd9e5d200, 0xefb3360b, 0x9426716a, 0xabddf8c2, 0xdd1ed9e4,
+ 0x17e1d567, 0xa9a65000, 0x2f37dbc5, 0x9a4b8fd5, 0xaeb22492, 0x0ebe8845,
+ 0xd89dd090, 0xcfbb88c6, 0xb1325561, 0x6d811d90, 0x03aa86f4, 0xbddba397,
+ 0x0986b9ed, 0x6f4cfc69, 0xc02b43bc, 0xee916274, 0xde7d9659, 0x7d3afd93,
+ 0xf52a7095, 0xf21a009c, 0xfd3f795e, 0x98cef25b, 0x6cb3af61, 0x6fa0e310,
+ 0x0196d036, 0xbc198bca, 0x15b0412d, 0xde454349, 0x5719472b, 0x8244ebce,
+ 0xee61afc6, 0xa60c9cb5, 0x1f4d1fd0, 0xe4fb3059, 0xab9ec0f9, 0x8d8b0255,
+ 0x4e7430bf, 0x3a22aa6b, 0x27de22d3, 0x60c4b6e6, 0x0cf61eb3, 0x469a87df,
+ 0xa4da1388, 0xf650f6aa, 0x3db87d68, 0xcdb6964c, 0xb2649b6c, 0x6a880fa9,
+ 0x1b0c845b, 0xe0af2f28, 0xfc1d5da9, 0xf64878a6, 0x667ca525, 0x2114b1ce,
+ 0x2d119ae3, 0x8d29d3bf, 0x1a1b4922, 0x3132980e, 0xd59e4385, 0x4dbd49b8,
+ 0x2de0bb05, 0xd6c96598, 0xb4c527c3, 0xb5562afc, 0x61eeb602, 0x05aa192a,
+ 0x7d127e77, 0xc719222d, 0xde7cf8db, 0x2de439b8, 0x250b5f1a, 0xd7b21053,
+ 0xef6c14a1, 0x2041f80f, 0xc287332e, 0xbb1dbfd3, 0x783bb979, 0x9a2e6327,
+ 0x6eb03027, 0x0225fa2f, 0xa319bc89, 0x864112d4, 0xfe990445, 0xe5e2e07c,
+ 0xf7c6acb8, 0x1bc92142, 0x12e9b40e, 0x2979282d, 0x05278e70, 0xe160ba4c,
+ 0xc1de0909, 0x458b9bf4, 0xbfce9c94, 0xa276f72a, 0x8441597d, 0x67adc2da,
+ 0x6162b854, 0x7f9b2f4a, 0x0d995b6b, 0x193b643d, 0x399362b3, 0x8b653a4b,
+ 0x1028d2db, 0x2b3df842, 0x6eecafaf, 0x261667e9, 0x9c7e8cda, 0x46063eab,
+ 0x7ce7a3a1, 0xadc899c9, 0x017291c4, 0x528d1a93, 0x9a1ee498, 0xbb7d4d43,
+ 0x7837f0ed, 0x34a230cc, 0x614a628d, 0xb03f93b8, 0xd72e3b08, 0x604c98db,
+ 0x3cfacb79, 0x8b81646a, 0xc0f082fa, 0xd1f92388, 0xe5a91e39, 0xf95c756d,
+ 0x1177742f, 0xf8819323, 0x5c060b80, 0x96c1cd8f, 0x47d7b440, 0xbbb84197,
+ 0x35f749cc, 0x95b0e132, 0x8d90ad54, 0x5c3f9423, 0x4994005b, 0xb58f53b9,
+ 0x32df7348, 0x60f61c29, 0x9eae2f32, 0x85a3d398, 0x3b995dd4, 0x94c5e460,
+ 0x8e54b9f3, 0x87bc6e2a, 0x90bbf1ea, 0x55d44719, 0x2cbbfe6e, 0x439d82f0,
+ 0x4eb3782d, 0xc3f1e669, 0x61ff8d9e, 0x0909238d, 0xef406165, 0x09c1d762,
+ 0x705d184f, 0x188f2cc4, 0x9c5aa12a, 0xc7a5d70e, 0xbc78cb1b, 0x1d26ae62,
+ 0x23f96ae3, 0xd456bf32, 0xe4654f55, 0x31462bd8
+ ]
View
@@ -11,21 +11,26 @@ import Data.Array.Vector
import Statistics.RandomVariate
import qualified System.Random as R
import qualified System.Random.Mersenne as M
+import Text.Printf
+import System.IO
+import MWC as MWC
+
+type T = Word32
getTime :: IO Double
getTime = (fromRational . toRational) `fmap` getPOSIXTime
-time :: IO a -> IO (Double, a)
-time act = do
+time :: String -> IO (a,Int) -> IO ()
+time desc act = do
+ putStr (desc ++ ": ")
s <- getTime
- ret <- act
- e <- ret `seq` getTime
- return (e-s, ret)
+ (ret,n) <- act
+ e <- getTime
+ let t = e-s
+ hPrintf stdout "%.3g secs (%.3g M/sec)\n" t (fromIntegral n/t/1e6)
count = 100 * 1000000
-type T = Double
-
--summ :: [T] -> T
--summ = foldl' (+) 0
@@ -37,42 +42,52 @@ loop n act = go 0 0
d <- act
go (i+1) (s+d)
-system :: IO Double
+system :: IO (Double, Int)
system = do
gen <- R.getStdGen
let a :*: g = loop 0 0 gen
R.setStdGen g
- return a
+ a `seq` return (a, k)
where loop !i !a !g | i == k = a :*: g
| otherwise = let (b,h) = R.random g
in loop (i+1) (a+b) h
k = count `div` 100
+immut :: Word32 -> IO (Word32, Int)
+immut k = let r = loop (MWC.init k) 0 0
+ in r `seq` return (r,n)
+ where loop !g !i !a | i == n = a
+ | otherwise = let (b:*:h) = MWC.next g
+ in loop h (i+1) (a+b)
+ n = count `div` 10
+
-mwc :: Word32 -> IO T
-mwc k = return $! runST go where
- go = do
- gen <- initialize (singletonU k)
- loop count (uniform gen)
+mwc :: Word32 -> IO (T, Int)
+mwc k = let r = runST go
+ in r `seq` return (r,count)
+ where
+ go = do
+ gen <- initialize (singletonU k)
+ loop count (uniform gen)
-mwca :: Word32 -> IO T
-mwca k = return $! runST go where
- go = do
- gen <- initialize (singletonU k)
- sumU `fmap` uniformArray gen count
+mwca :: Word32 -> IO (T, Int)
+mwca k = let r = runST go
+ in r `seq` return (r,count)
+ where
+ go = do
+ gen <- initialize (singletonU k)
+ sumU `fmap` uniformArray gen count
-mersenne :: IO T
+mersenne :: IO (T, Int)
mersenne = do
gen <- M.getStdGen
- loop count (M.random gen)
+ r <- loop count (M.random gen)
+ return (r,count)
main = do
forM_ [1..3] $ \n -> do
- putStr "system: "
- print =<< time system
- putStr "mwc: "
- print =<< time (mwc n)
- putStr "mwca: "
- print =<< time (mwca n)
- putStr "mersenne: "
- print =<< time mersenne
+ time "system" system
+ time "immut" (immut n)
+ time "mwc" (mwc n)
+ time "mwca" (mwca n)
+ time "mersenne" mersenne

0 comments on commit 216ddcb

Please sign in to comment.