content: 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 {-# UNPACK #-} !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 = (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
  --bConst = 4077358422479273989
  mConst = 18446383549859758079
  skipConst = 1073741824
  