module Quant.RNG.MWC64X (
MWC64X(..)
, randomWord32
, randomWord64
, randomInt
, randomDouble
, randomInt64
, skip
) where
import Data.Int
import Data.Bits
import Data.Word
import Data.Random.Internal.Words
import System.Random
data MWC64X = MWC64X !Word64 deriving (Eq,Show)
randomWord32 :: MWC64X -> (Word32, MWC64X)
randomWord32 (MWC64X state) = (x `xor` c, MWC64X state')
where c = fromIntegral $ state `shiftR` 32 :: Word32
x = fromIntegral $ state .&. 0xFFFFFFFF :: Word32
state' = fromIntegral x * aConst + fromIntegral c
randomInt :: MWC64X -> (Int,MWC64X)
randomInt g = (fromIntegral i, g')
where (i, g') = randomWord64 g
randomWord64 :: MWC64X -> (Word64, MWC64X)
randomWord64 x = (buildWord64'' y1 y2, x'')
where (y1, x' ) = randomWord32 x
(y2, x'') = randomWord32 x'
randomDouble :: MWC64X -> (Double, MWC64X)
randomDouble x = (fromIntegral (val `div` 2048) / 9007199254740992, x')
where (val, x') = randomWord64 x
randomInt64 :: MWC64X -> (Int64,MWC64X)
randomInt64 g = (fromIntegral i, g')
where (i, g') = randomWord64 g
instance RandomGen MWC64X where
next g = (fromIntegral w, g')
where (w, g') = randomWord64 g
split g = (skip g skipConst, g)
addMod64 :: Word64 -> Word64 -> Word64 -> Word64
addMod64 a b m = (a+b) `mod` m
mulMod64 :: Word64 -> Word64 -> Word64 -> Word64
mulMod64 a b m = f 0 a b
where f r a1 b1
| a1 == 0 = r
| otherwise = f r' a' b'
where r' = if a1 .&. 1 == 1 then addMod64 r b1 m else r
b' = addMod64 b1 b1 m
a' = a `shiftR` 1
powMod64 :: Word64 -> Word64 -> Word64 -> Word64
powMod64 a e m = f a 1 e
where f sqr acc e1
| e1 == 0 = acc
| otherwise = f sqr' acc' e'
where acc' = if e1 .&. 1 == 1 then mulMod64 acc sqr m else acc
sqr' = mulMod64 sqr sqr m
e' = e `shiftR` 1
skip :: MWC64X -> Word64 -> MWC64X
skip (MWC64X st) d = MWC64X st'
where
m = powMod64 aConst d mConst
c = st `shiftR` 32
x = st .&. 0xFFFFFFFF
x' = mulMod64 (x * aConst + c) m mConst
x'' = fromIntegral $ x' `div` aConst :: Word32
c' = fromIntegral $ x' `mod` aConst :: Word32
st' = buildWord64'' c' x''
mkWord64 :: Word32 -> Word32 -> Word64
mkWord64 a b = (fromIntegral $ a `shiftL` 32) .&. fromIntegral b
aConst, mConst, skipConst :: Word64
aConst = 4294883355
mConst = 18446383549859758079
skipConst = 1073741824