module System.Random.MWC
(
Gen
, GenIO
, GenST
, Seed
, fromSeed
, toSeed
, Variate(..)
, normal
, create
, initialize
, withSystemRandom
, save
, restore
, uniformVector
) where
#if defined(__GLASGOW_HASKELL__) && !defined(__HADDOCK__)
#include "MachDeps.h"
#endif
import Control.Exception (IOException, catch)
import Control.Monad (ap, liftM, unless)
import Control.Monad.Primitive (PrimMonad, PrimState, unsafePrimToIO)
import Control.Monad.ST (ST)
import Data.Bits (Bits, (.&.), (.|.), shiftL, shiftR, xor)
import Data.Int (Int8, Int16, Int32, Int64)
import Data.IORef (atomicModifyIORef, newIORef)
import Data.Ratio ((%), numerator)
import Data.Time.Clock.POSIX (getPOSIXTime)
import Data.Typeable (Typeable)
import Data.Vector.Generic (Vector, unsafeFreeze)
import Data.Word (Word, Word8, Word16, Word32, Word64)
import Foreign.Marshal.Alloc (allocaBytes)
import Foreign.Marshal.Array (peekArray)
import Prelude hiding (catch)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM
import qualified Data.Vector.Unboxed as I
import qualified Data.Vector.Unboxed.Mutable as M
import System.CPUTime (cpuTimePrecision, getCPUTime)
import System.IO (IOMode(..), hGetBuf, hPutStrLn, stderr, withBinaryFile)
import System.IO.Unsafe (unsafePerformIO)
class M.Unbox a => Variate a where
uniform :: (PrimMonad m) => Gen (PrimState m) -> m a
uniformR :: (PrimMonad m) => (a,a) -> Gen (PrimState m) -> m a
instance Variate Int8 where
uniform = uniform1 fromIntegral
uniformR = uniformRange
instance Variate Int16 where
uniform = uniform1 fromIntegral
uniformR = uniformRange
instance Variate Int32 where
uniform = uniform1 fromIntegral
uniformR = uniformRange
instance Variate Int64 where
uniform = uniform2 wordsTo64Bit
uniformR = uniformRange
instance Variate Word8 where
uniform = uniform1 fromIntegral
uniformR = uniformRange
instance Variate Word16 where
uniform = uniform1 fromIntegral
uniformR = uniformRange
instance Variate Word32 where
uniform = uniform1 fromIntegral
uniformR = uniformRange
instance Variate Word64 where
uniform = uniform2 wordsTo64Bit
uniformR = uniformRange
instance Variate Bool where
uniform = uniform1 wordToBool
uniformR (False,True) g = uniform g
uniformR (False,False) _ = return False
uniformR (True,True) _ = return True
uniformR (True,False) g = uniform g
instance Variate Float where
uniform = uniform1 wordToFloat
uniformR (x1,x2) = uniform1 (\w -> x1 + (x2x1) * wordToFloat w)
instance Variate Double where
uniform = uniform2 wordsToDouble
uniformR (x1,x2) = uniform2 (\w1 w2 -> x1 + (x2x1) * wordsToDouble w1 w2)
instance Variate Int where
#if WORD_SIZE_IN_BITS < 64
uniform = uniform1 fromIntegral
#else
uniform = uniform2 wordsTo64Bit
#endif
uniformR = uniformRange
instance Variate Word where
#if WORD_SIZE_IN_BITS < 64
uniform = uniform1 fromIntegral
#else
uniform = uniform2 wordsTo64Bit
#endif
uniformR = uniformRange
instance (Variate a, Variate b) => Variate (a,b) where
uniform g = (,) `liftM` uniform g `ap` uniform g
uniformR ((x1,y1),(x2,y2)) g = (,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g
instance (Variate a, Variate b, Variate c) => Variate (a,b,c) where
uniform g = (,,) `liftM` uniform g `ap` uniform g `ap` uniform g
uniformR ((x1,y1,z1),(x2,y2,z2)) g =
(,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` uniformR (z1,z2) g
instance (Variate a, Variate b, Variate c, Variate d) => Variate (a,b,c,d) where
uniform g = (,,,) `liftM` uniform g `ap` uniform g `ap` uniform g
`ap` uniform g
uniformR ((x1,y1,z1,t1),(x2,y2,z2,t2)) g =
(,,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap`
uniformR (z1,z2) g `ap` uniformR (t1,t2) g
wordsTo64Bit :: (Integral a) => Word32 -> Word32 -> a
wordsTo64Bit x y =
fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64)
wordToBool :: Word32 -> Bool
wordToBool i = (i .&. 1) /= 0
wordToFloat :: Word32 -> Float
wordToFloat x = (fromIntegral i * m_inv_32) + 0.5 + m_inv_33
where m_inv_33 = 1.16415321826934814453125e-10
m_inv_32 = 2.3283064365386962890625e-10
i = fromIntegral x :: Int32
wordsToDouble :: Word32 -> Word32 -> Double
wordsToDouble x y = (fromIntegral u * m_inv_32 + (0.5 + m_inv_53) +
fromIntegral (v .&. 0xFFFFF) * m_inv_52)
where m_inv_52 = 2.220446049250313080847263336181640625e-16
m_inv_53 = 1.1102230246251565404236316680908203125e-16
m_inv_32 = 2.3283064365386962890625e-10
u = fromIntegral x :: Int32
v = fromIntegral y :: Int32
newtype Gen s = Gen (M.MVector s Word32)
type GenIO = Gen (PrimState IO)
type GenST s = Gen (PrimState (ST s))
ioff, coff :: Int
ioff = 256
coff = 257
create :: PrimMonad m => m (Gen (PrimState m))
create = initialize defaultSeed
initialize :: (PrimMonad m, Vector v Word32) =>
v Word32 -> m (Gen (PrimState m))
initialize seed = do
q <- M.unsafeNew 258
fill q
if fini == 258
then do
M.unsafeWrite q ioff $ G.unsafeIndex seed ioff .&. 255
M.unsafeWrite q coff $ G.unsafeIndex seed coff
else do
M.unsafeWrite q ioff 255
M.unsafeWrite q coff 362436
return (Gen q)
where fill q = go 0 where
go i | i == 256 = return ()
| otherwise = M.unsafeWrite q i s >> go (i+1)
where s | i >= fini = if fini == 0
then G.unsafeIndex defaultSeed i
else G.unsafeIndex defaultSeed i `xor`
G.unsafeIndex seed (i `mod` fini)
| otherwise = G.unsafeIndex seed i
fini = G.length seed
newtype Seed = Seed {
fromSeed :: I.Vector Word32
}
deriving (Eq, Show, Typeable)
toSeed :: (Vector v Word32) => v Word32 -> Seed
toSeed v = Seed $ I.create $ do { Gen q <- initialize v; return q }
safeFreeze :: (PrimMonad m, Vector v a) => G.Mutable v (PrimState m) a -> m (v a)
safeFreeze v = do
v' <- GM.unsafeNew (GM.length v)
GM.unsafeCopy v' v
unsafeFreeze v'
save :: PrimMonad m => Gen (PrimState m) -> m Seed
save (Gen q) = Seed `liftM` safeFreeze q
restore :: PrimMonad m => Seed -> m (Gen (PrimState m))
restore (Seed s) = M.unsafeNew n >>= fill
where fill q = go 0 where
go !i | i >= n = return $! Gen q
| otherwise = M.unsafeWrite q i (I.unsafeIndex s i) >> go (i+1)
n = I.length s
withTime :: (PrimMonad m) => (Gen (PrimState m) -> m a) -> IO a
withTime act = do
c <- (numerator . (%cpuTimePrecision)) `liftM` getCPUTime
t <- toRational `liftM` getPOSIXTime
let n = fromIntegral (numerator t) :: Word64
seed = [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)]
unsafePrimToIO $ initialize (I.fromList seed) >>= act
withSystemRandom :: PrimMonad m => (Gen (PrimState m) -> m a) -> IO a
withSystemRandom act = tryRandom `catch` \(_::IOException) -> do
seen <- atomicModifyIORef warned ((,) True)
unless seen $ do
hPutStrLn stderr ("Warning: Couldn't open " ++ show random)
hPutStrLn stderr ("Warning: using system clock for seed instead " ++
"(quality will be lower)")
withTime act
where tryRandom = do
let nbytes = 1024
ws <- allocaBytes nbytes $ \buf -> do
nread <- withBinaryFile random ReadMode $
\h -> hGetBuf h buf nbytes
peekArray (nread `div` 4) buf
unsafePrimToIO $ initialize (I.fromList ws) >>= act
random = "/dev/urandom"
warned = unsafePerformIO $ newIORef False
nextIndex :: Integral a => a -> Int
nextIndex i = fromIntegral j
where j = fromIntegral (i+1) :: Word8
a :: Word64
a = 1540315826
uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32
uniformWord32 (Gen q) = do
i <- nextIndex `liftM` M.unsafeRead q ioff
c <- fromIntegral `liftM` M.unsafeRead q coff
qi <- fromIntegral `liftM` M.unsafeRead q i
let t = a * qi + c
c' = fromIntegral (t `shiftR` 32)
x = fromIntegral t + c'
(# x', c'' #) | x < c' = (# x + 1, c' + 1 #)
| otherwise = (# x, c' #)
M.unsafeWrite q i x'
M.unsafeWrite q ioff (fromIntegral i)
M.unsafeWrite q coff (fromIntegral c'')
return x'
uniform1 :: PrimMonad m => (Word32 -> a) -> Gen (PrimState m) -> m a
uniform1 f gen = do
i <- uniformWord32 gen
return $! f i
uniform2 :: PrimMonad m => (Word32 -> Word32 -> a) -> Gen (PrimState m) -> m a
uniform2 f (Gen q) = do
i <- nextIndex `liftM` M.unsafeRead q ioff
let j = nextIndex i
c <- fromIntegral `liftM` M.unsafeRead q coff
qi <- fromIntegral `liftM` M.unsafeRead q i
qj <- fromIntegral `liftM` M.unsafeRead q j
let t = a * qi + c
c' = fromIntegral (t `shiftR` 32)
x = fromIntegral t + c'
(# x', c'' #) | x < c' = (# x + 1, c' + 1 #)
| otherwise = (# x, c' #)
u = a * qj + fromIntegral c''
d' = fromIntegral (u `shiftR` 32)
y = fromIntegral u + d'
(# y', d'' #) | y < d' = (# y + 1, d' + 1 #)
| otherwise = (# y, d' #)
M.unsafeWrite q i x'
M.unsafeWrite q j y'
M.unsafeWrite q ioff (fromIntegral j)
M.unsafeWrite q coff (fromIntegral d'')
return $! f x' y'
type family Unsigned a :: *
type instance Unsigned Int8 = Word8
type instance Unsigned Int16 = Word16
type instance Unsigned Int32 = Word32
type instance Unsigned Int64 = Word64
type instance Unsigned Int = Word
type instance Unsigned Word8 = Word8
type instance Unsigned Word16 = Word16
type instance Unsigned Word32 = Word32
type instance Unsigned Word64 = Word64
type instance Unsigned Word = Word
sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a
sub x y = fromIntegral x fromIntegral y
add :: (Integral a, Integral (Unsigned a)) => a -> Unsigned a -> a
add m x = m + fromIntegral x
unsignedRange :: (PrimMonad m, Integral a, Bounded a) => a -> m a -> m a
unsignedRange n rnd = go
where
buckets = maxBound `div` n
maxN = buckets * n
go = do x <- rnd
if x < maxN then return (x `div` buckets)
else go
uniformRange :: ( PrimMonad m
, Integral a, Bounded a, Variate a
, Integral (Unsigned a), Bounded (Unsigned a), Variate (Unsigned a))
=> (a,a) -> Gen (PrimState m) -> m a
uniformRange (x1,x2) g
| x1 == minBound && x2 == maxBound = uniform g
| otherwise = do x <- unsignedRange (sub x2 x1 + 1) (uniform g)
return $! add x1 x
uniformVector :: (PrimMonad m, Variate a, Vector v a)
=> Gen (PrimState m) -> Int -> m (v a)
uniformVector gen n = G.replicateM n (uniform gen)
data T = T !Double !Double
normal :: PrimMonad m => Gen (PrimState m) -> m Double
normal gen = loop
where
loop = do
u <- (subtract 1 . (*2)) `liftM` uniform gen
ri <- uniform gen
let i = fromIntegral ((ri :: Word32) .&. 127)
bi = I.unsafeIndex blocks i
bj = I.unsafeIndex blocks (i+1)
if abs u < I.unsafeIndex ratios i
then return $! u * bi
else if i == 0
then normalTail (u < 0)
else do
let x = u * bi
xx = x * x
d = exp (0.5 * (bi * bi xx))
e = exp (0.5 * (bj * bj xx))
c <- uniform gen
if e + c * (d e) < 1
then return x
else loop
blocks = let f = exp (0.5 * r * r)
in (`I.snoc` 0) . I.cons (v/f) . I.cons r .
I.unfoldrN 126 go $! T r f
where
go (T b g) = let !u = T h (exp (0.5 * h * h))
h = sqrt (2 * log (v / b + g))
in Just (h, u)
v = 9.91256303526217e-3
r = 3.442619855899
ratios = I.zipWith (/) (I.tail blocks) blocks
normalTail neg = tailing
where tailing = do
x <- ((/r) . log) `liftM` uniform gen
y <- log `liftM` uniform gen
if y * (2) < x * x
then tailing
else return $! if neg then x r else r x
defaultSeed :: I.Vector Word32
defaultSeed = I.fromList [
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 ]