-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Random.Generator.MT19937
-- Copyright   :  (c) Matt Harden 1999
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- A Haskell program for MT19937 pseudorandom number generator
--
-----------------------------------------------------------------------------

-- The original source was found at
--
-- http://members.primary.net/~matth/mt19937.hs
--
-- but I can't get to the site anymore.  As much as the orginal
-- formatting has been retained as possible. --mpd

-- TODO: Make an instance of RandomGen class

{-
   Function genrand generates an infinite list of pseudorandom
   unsigned integers (32bit) which are uniformly distributed
   among 0 to 2^32-1.  sgenrand(seed) uses an algorithm of Knuth
   to provide 624 initial values to genrand().

   Rewritten in Haskell by Matt Harden
      from original code in C by Takuji Nishimura.

   This program relies upon the GHC/Hugs extensions to Haskell.
   These are very likely to be available in any Haskell
   environment, and performance would suffer greatly without them.
-}

{-
   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Library General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later
   version.
   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   See the GNU Library General Public License for more details.
   You should have received a copy of the GNU Library General
   Public License along with this library; if not, write to the
   Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307  USA
-}

-- Copyright (C) 1999 Matt Harden
-- The original C code contained the following notice:
--   When you use this, send an email to: matumoto@math.keio.ac.jp
--   with an appropriate reference to your work.

{- REFERENCE -
   M. Matsumoto and T. Nishimura,
   "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
   Pseudo-Random Number Generator",
   ACM Transactions on Modeling and Computer Simulation,
   Vol. 8, No. 1, January 1998, pp 3--30.
-}

module Numeric.Random.Generator.MT19937 (W, genrand, test) where

import Data.Word (Word32)
import Data.Bits (Bits, shiftL, shiftR, complement, xor, bit, (.|.), (.&.))

infixl 8 .<<., .>>.

(.<<.), (.>>.) :: (Bits a) => (a -> Int -> a)
.<<. :: forall a. Bits a => a -> Int -> a
(.<<.) = forall a. Bits a => a -> Int -> a
shiftL
.>>. :: forall a. Bits a => a -> Int -> a
(.>>.) = forall a. Bits a => a -> Int -> a
shiftR

type W = Word32

-- Period parameters
parmN :: Int
parmN :: Int
parmN = Int
624
parmM :: Int
parmM :: Int
parmM = Int
397
parmA :: W
parmA :: W
parmA = W
0x9908b0df
upperMask :: W
upperMask :: W
upperMask = (forall a. Bits a => Int -> a
bit Int
31)
lowerMask :: W
lowerMask :: W
lowerMask = (forall a. Bits a => a -> a
complement W
upperMask)

-- Tempering parameters
temperingMaskB :: W -> W
temperingMaskB :: W -> W
temperingMaskB = (forall a. Bits a => a -> a -> a
.&. W
0x9d2c5680)
temperingMaskC :: W -> W
temperingMaskC :: W -> W
temperingMaskC = (forall a. Bits a => a -> a -> a
.&. W
0xefc60000)
temperingShiftU :: W -> W
temperingShiftU :: W -> W
temperingShiftU = (forall a. Bits a => a -> Int -> a
.>>. Int
11)
temperingShiftS :: W -> W
temperingShiftS :: W -> W
temperingShiftS = (forall a. Bits a => a -> Int -> a
.<<.  Int
7)
temperingShiftT :: W -> W
temperingShiftT :: W -> W
temperingShiftT = (forall a. Bits a => a -> Int -> a
.<<. Int
15)
temperingShiftL :: W -> W
temperingShiftL :: W -> W
temperingShiftL = (forall a. Bits a => a -> Int -> a
.>>. Int
18)

-- A Knuth algorithm just to seed the seed...
-- Line 25 of table 1
-- in [KNUTH 1981, The Art of Computer Programming Vol. 2 (2nd Ed.), pp102]
sgenrand :: W -> [W]
sgenrand :: W -> [W]
sgenrand W
0 = W -> [W]
sgenrand W
4357   -- 0 not acceptable.  Why 4357?  I dunno.
sgenrand W
seed = forall a. Int -> [a] -> [a]
take Int
parmN (forall a. (a -> a) -> a -> [a]
iterate (W
69069 forall a. Num a => a -> a -> a
*) W
seed)

mag01 :: Bool -> W
mag01 :: Bool -> W
mag01 Bool
False = W
0
mag01 Bool
True  = W
parmA

tempering :: W -> W
tempering :: W -> W
tempering = let ^= :: t -> (t -> t) -> t
(^=) t
x t -> t
f = forall a. Bits a => a -> a -> a
xor t
x (t -> t
f t
x) in
   (forall {t}. Bits t => t -> (t -> t) -> t
^= (W -> W
temperingShiftL)) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (forall {t}. Bits t => t -> (t -> t) -> t
^= (W -> W
temperingMaskC forall b c a. (b -> c) -> (a -> b) -> a -> c
. W -> W
temperingShiftT)) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (forall {t}. Bits t => t -> (t -> t) -> t
^= (W -> W
temperingMaskB forall b c a. (b -> c) -> (a -> b) -> a -> c
. W -> W
temperingShiftS)) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (forall {t}. Bits t => t -> (t -> t) -> t
^= (W -> W
temperingShiftU))

-- parameter to rand MUST be a list of (_N) words!
rand :: [W] -> [W]
rand :: [W] -> [W]
rand [W]
seed = forall a b. (a -> b) -> [a] -> [b]
map W -> W
tempering [W]
r2 where
   r :: [W]
r = [W]
seed forall a. [a] -> [a] -> [a]
++ [W]
r2
   r2 :: [W]
r2 = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Bits a => a -> a -> a
xor (forall a b. (a -> b) -> [a] -> [b]
map W -> W
f [W]
r3) (forall a. Int -> [a] -> [a]
drop Int
parmM [W]
r)
   r3 :: [W]
r3 = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\W
x W
y -> (W
x forall a. Bits a => a -> a -> a
.&. W
upperMask) forall a. Bits a => a -> a -> a
.|. (W
y forall a. Bits a => a -> a -> a
.&. W
lowerMask)) [W]
r (forall a. [a] -> [a]
tail [W]
r)
   f :: W -> W
f W
y = (W
y forall a. Bits a => a -> Int -> a
.>>. Int
1) forall a. Bits a => a -> a -> a
`xor` (Bool -> W
mag01 (forall a. Integral a => a -> Bool
odd W
y))

genrand :: W -> [W]
genrand :: W -> [W]
genrand = [W] -> [W]
rand forall b c a. (b -> c) -> (a -> b) -> a -> c
. W -> [W]
sgenrand

test :: IO ()
test :: IO ()
test = forall (t :: * -> *) (m :: * -> *) a.
(Foldable t, Monad m) =>
t (m a) -> m ()
sequence_ forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a. Show a => a -> IO ()
print forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
take Int
1000 forall a b. (a -> b) -> a -> b
$ W -> [W]
genrand W
4357