Mersenne.hs

{-# LANGUAGE CPP, ForeignFunctionInterface, BangPatterns #-}
--------------------------------------------------------------------
-- |
-- Module     : System.Random.Mersenne
-- Copyright  : Copyright (c) 2008, Don Stewart <dons@galois.com>
-- License    : BSD3
-- Maintainer : Don Stewart <dons@galois.com>
-- Stability  : experimental
-- Portability: CPP, FFI
-- Tested with: GHC 6.8.2
--
-- Generate pseudo-random numbers using the SIMD-oriented Fast Mersenne Twister(SFMT)
-- pseudorandom number generator. This is a /much/ faster generator than
-- the default 'System.Random' generator for Haskell (~50x faster
-- generation for Doubles, on a core 2 duo), however, it is not 
-- nearly as flexible.
--
-- This library may be compiled with the '-f use_sse2' or '-f
-- use_altivec' flags to configure, on intel and powerpc machines
-- respectively, to enable high performance vector instructions to be used.
-- This typically results in a 2-3x speedup in generation time.
--
-- This will work for newer intel chips such as Pentium 4s, and Core, Core2* chips.
--
module System.Random.Mersenne (

    -- * The random number generator
    MTGen

    -- ** Initialising the generator
    , newMTGen

    -- * Random values of various types

    -- $notes
    , MTRandom(..)

    -- $globalrng
    , getStdRandom
    , getStdGen
    , setStdGen

    -- * Miscellaneous
    , version

    -- $example

    ) where

#if defined(__GLASGOW_HASKELL__) && !defined(__HADDOCK__)
#include "MachDeps.h"
#endif

import Foreign.C.Types
import Foreign.C.String

import System.CPUTime	( getCPUTime )
import System.Time
import System.IO.Unsafe
-- import Control.Monad

import Data.Word
import Data.Int
import Data.Bits
-- import Data.Char
import Data.IORef

------------------------------------------------------------------------
-- $example
--
-- An example, calculation of pi via a monte carlo method:
--
-- > import System.Random.Mersenne
-- > import System.Environment
--
-- We'll roll the dice 'lim' times, 
--
-- > main = do
-- >    [lim] <- mapM readIO =<< getArgs
--
-- Now, define a loop that runs this many times, plotting a 'x' and 'y'
-- position, then working out if its in and outside the circle. 
-- The ratio of inside\/total points at then gives us an approximation
-- of pi.
--
-- > let go :: Int -> Int -> IO Double
-- >     go throws ins
-- >         | throws >= lim  = return ((4 * fromIntegral ins) / (fromIntegral throws))
-- >         | otherwise = do
-- >             x <- random g :: IO Double
-- >             y <- random g :: IO Double
-- >             if x * x + y * y < 1
-- >                 then go (throws+1) $! ins + 1
-- >                 else go (throws+1) ins
--
-- Compiling this, '-fexcess-precision', for accurate Doubles,
--
-- > $ ghc -fexcess-precision -fvia-C pi.hs -o pi
-- > $ ./pi 10000000                                                 
-- > 3.1417304
--

------------------------------------------------------------------------

-- | A single, global SIMD fast mersenne twister random number generator
-- This generator is evidence that you have initialised the generator,
--
data MTGen = MTGen

-- | Return an initialised SIMD Fast Mersenne Twister.
-- The generator is initialised based on the clock time, if Nothing
-- is passed as a seed. For deterministic behaviour, pass an explicit seed.
--
-- Due to the current SFMT library being vastly impure, currently only a single
-- generator is allowed per-program. Attempts to reinitialise it will fail.
--
newMTGen :: Maybe Word32 -> IO MTGen
newMTGen (Just n) = do
    dup <- c_get_initialized
    if dup == 0
       then do c_init_gen_rand (fromIntegral n)
               return MTGen

       else error $ "System.Random.Mersenne: " ++
                    "Only one mersenne twister generator can be created per process"

newMTGen Nothing = do
    ct             <- getCPUTime
    (TOD sec psec) <- getClockTime
    newMTGen (Just (fromIntegral $ sec * 1013904242 + psec + ct) )

------------------------------------------------------------------------

-- $notes
--
-- Instances MTRandom for Word, Word64, Word32, Word16, Word8
-- all return, quickly, a random inhabintant of that type, in its full
-- range. Similarly for Int types.
--
-- Int and Word will be 32 bits on a 32 bit machine, and 64 on a 64 bit
-- machine. The double precision will be 32 bits on a 32 bit machine,
-- and 53 on a 64 bit machine.
--
-- The MTRandom instance for Double returns a Double in the interval [0,1).
-- The Bool instance takes the lower bit off a random word.

-- | Given an initialised SFMT generator, the MTRandom
-- allows the programmer to extract values of a variety of 
-- types.
--
-- Minimal complete definition: 'random'.
--
class MTRandom a where

  -- | The same as 'randomR', but using a default range determined by the type:
  --
  -- * For bounded types (instances of 'Bounded', such as 'Char'),
  --   the range is normally the whole type.
  --
  -- * For fractional types, the range is normally the semi-closed interval
  -- @[0,1)@.
  --
  -- * For 'Integer', the range is (arbitrarily) the range of 'Int'.
  random :: MTGen -> IO a

  -- | Plural variant of 'random', producing an infinite list of
  -- random values instead of returning a new generator.
  randoms  :: MTGen -> IO [a]
  randoms !g = unsafeInterleaveIO $ do
                x  <- random g
                xs <- randoms g
                return (x : xs)
  -- There are real overheads here. Consider eagerly filling chunks
  -- and extracting elements piecewise.

{-
  -- | Takes a range /(lo,hi)/ and a random number generator
  -- /g/, and returns a random value uniformly distributed in the closed
  -- interval /[lo,hi]/, together with a new generator. It is unspecified
  -- what happens if /lo>hi/. For continuous types there is no requirement
  -- that the values /lo/ and /hi/ are ever produced, but they may be,
  -- depending on the implementation and the interval.
  randomR :: (a,a) -> MTGen -> IO a
-}

{-
  -- | Plural variant of 'random', producing an infinite list of
  -- random values instead of returning a new generator.
  randomRs  :: (a,a) -> MTGen -> IO [a]
  randomRs p !g = unsafeInterleaveIO $ do
                x  <- randomR p g
                xs <- randomRs p g
                return (x : xs)
-}

  -- | A variant of 'random' that uses the global random number generator
  -- (see "System.Random#globalrng").
  -- Essentially a convenience function if you're already in IO.
  --
  -- Note that there are performance penalties calling randomIO in an
  -- inner loop, rather than 'random' applied to a global generator. The
  -- cost comes in retrieving the random gen from an IORef, which is
  -- non-trivial. Expect a 3x slow down in speed of random generation.
  randomIO :: IO a
  randomIO = getStdRandom random
  {-# INLINE randomIO #-}

------------------------------------------------------------------------
-- Efficient basic instances

instance MTRandom Word    where
    random !_          = randomWord
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Word64  where
    random !_ = fmap fromIntegral randomWord64
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Word32  where
    random !_ = fmap fromIntegral randomWord
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Word16  where
    random !_ = fmap fromIntegral randomWord
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Word8   where
    random !_          = fmap fromIntegral randomWord
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)

------------------------------------------------------------------------

instance MTRandom Double where
    random !_ = randomDouble
    {-# INLINE random #-}
--    randomR (lo,hi) g = randomIvalDouble g lo hi id

------------------------------------------------------------------------

instance MTRandom Int     where
    random !_          = randomInt
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalInt g lo hi

instance MTRandom Int64   where
    random !_ = fmap fromIntegral randomInt64
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Int32   where
    random !_ = fmap fromIntegral randomInt
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Int16   where
    random !_ = fmap fromIntegral randomInt
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Int8    where
    random !_ = fmap fromIntegral randomInt
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)

instance MTRandom Integer where
    random !_ = fmap fromIntegral randomInt
    {-# INLINE random #-}
--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)

------------------------------------------------------------------------

instance MTRandom Bool    where
    random !_ = do x <- randomWord; return $! x .&. 1 /= 0
    {-# INLINE random #-}

{-
    randomR (a,b) !g = int2Bool `fmap` randomIvalInt g (bool2Int a) (bool2Int b)

        where
            bool2Int :: Bool -> Int
            bool2Int False = 0
            bool2Int True  = 1

            int2Bool :: Int -> Bool
            int2Bool 0 = False
            int2Bool _ = True
-}

------------------------------------------------------------------------

{-
randomIvalInt :: (MTRandom a, Num a) => MTGen -> Int -> Int -> IO a
randomIvalInt g l h
  | l > h     = randomIvalInt g h l
  | otherwise = do
        v <- f n 1
        return $ (fromIntegral (l + v `mod` k))
     where
       k = h - l + 1
       b = maxBound :: Int
       n = iLogBase b k

       f 0 acc = return acc
       f i acc = do
           x <- random g :: IO Int
           f (i-1) (fromIntegral x + acc * b)
{-# INLINE randomIvalInt #-}

iLogBase :: Int -> Int -> Int
iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
-}

------------------------------------------------------------------------

{-
randomIvalWord :: (MTRandom a, Num a) => MTGen -> Word -> Word -> IO a
randomIvalWord g l h
  | l > h     = randomIvalWord g h l
  | otherwise = do
        v <- f n 1
        return $ (fromIntegral (l + v `mod` k))
     where
       k = h - l + 1
       b = maxBound :: Word
       n = iLogBaseWord b k

       f 0 acc = return acc
       f i acc = do
           x <- random g :: IO Word
           f (i-1) (fromIntegral x + acc * b)
{-# INLINE randomIvalWord #-}

iLogBaseWord :: Word -> Word -> Word
iLogBaseWord b i = if i < b then 1 else 1 + iLogBaseWord b (i `div` b)
-}

------------------------------------------------------------------------

{-
--
-- Too slow:
--
randomIvalDouble :: (MTRandom a, Fractional a) => MTGen -> Double -> Double -> (Double -> a) -> IO a
randomIvalDouble g l h fromDouble
  | l > h     = randomIvalDouble g h l fromDouble
  | otherwise = do
       x <- random g :: IO Int
       return $ fromDouble ((l+h)/2) +
                fromDouble ((h-l) / realToFrac intRange) *
                fromIntegral x
{-# INLINE randomIvalDouble #-}

intRange :: Integer
intRange = toInteger (maxBound::Int) - toInteger (minBound::Int)
-}
------------------------------------------------------------------------
--
-- Using a single global random number generator
--

{- $globalrng #globalrng#

There is a single, implicit, global random number generator of type
'StdGen', held in some global variable maintained by the 'IO' monad. It is
initialised automatically in some system-dependent fashion. To get
deterministic behaviour, use 'setStdGen'.
-}

theStdGen :: IORef MTGen
theStdGen  = unsafePerformIO $ do
   rng <- newMTGen Nothing
   newIORef rng
{-# NOINLINE theStdGen #-}

-- |Sets the global random number generator.
setStdGen :: MTGen -> IO ()
setStdGen = writeIORef theStdGen

-- |Gets the global random number generator.
getStdGen :: IO MTGen
getStdGen  = readIORef theStdGen

-- | Uses the supplied function to get a value from the current global
-- random generator, and updates the global generator with the new
-- generator returned by the function. For example, @rollDice@ gets a
-- random integer between 1 and 6:
--
-- >  rollDice :: IO Int
-- >  rollDice = getMTRandom (randomR (1,6))
--
getStdRandom :: (MTGen -> IO a) -> IO a
getStdRandom f = do
    st <- readIORef theStdGen
    f st
{-# INLINE getStdRandom #-}

------------------------------------------------------------------------

-- | Returns the identification string for the SMFT version. 
-- The string shows the word size, the Mersenne exponent, and all parameters of this generator.
version :: String
version = unsafePerformIO (peekCString =<< c_get_idstring)

------------------------------------------------------------------------
-- Safe primitives: depend on the word size. It's generally not a 
-- good idea to mix generation of different types, unless you commit
-- to either 32 or 64 bits only.
--
-- So you should only use these functions for getting at randoms.

randomInt :: IO Int
randomInt = fmap fromIntegral
#if WORD_SIZE_IN_BITS < 64
    c_gen_rand32
#else
    c_gen_rand64
#endif

-- TODO randomWord64, for 32 bit machines

randomWord :: IO Word
randomWord = fmap fromIntegral
#if WORD_SIZE_IN_BITS < 64
    c_gen_rand32
#else
    c_gen_rand64
#endif

randomWord64 :: IO Word64
randomWord64 = fmap fromIntegral
#if WORD_SIZE_IN_BITS < 64
    c_gen_rand64_mix
#else
    c_gen_rand64
#endif

randomInt64 :: IO Int64
randomInt64 = fmap fromIntegral
#if WORD_SIZE_IN_BITS < 64
    c_gen_rand64_mix
#else
    c_gen_rand64
#endif

randomDouble :: IO Double
randomDouble = fmap realToFrac
#if WORD_SIZE_IN_BITS < 64
    c_genrand_real2
#else
    c_genrand_res53
#endif

------------------------------------------------------------------------
-- Generating chunks at a time.
--

{-
min_array_size :: Int
min_array_size = fromIntegral . unsafePerformIO $ -- constant
#if WORD_SIZE_IN_BITS < 64
    c_get_min_array_size32
#else
    c_get_min_array_size64
#endif

-- | Fill an array with 'n' random Ints
fill_array :: Ptr Int -> Int -> IO ()
fill_array p n =
#if WORD_SIZE_IN_BITS < 64
    c_fill_array32 (castPtr p) (fromIntegral n)
#else
    c_fill_array64 (castPtr p) (fromIntegral n)
#endif
-}

------------------------------------------------------------------------
-- We can have only one mersenne supply in a program.

-- You have to commit at initialisation time to call either
-- rand_gen32 or rand_gen64, and correspondingly, real2 or res53
-- for doubles.
--

type UInt32 = CUInt
type UInt64 = CULLong


-- | This function initializes the internal state array with a 32-bit integer seed.
foreign import ccall unsafe "SFMT.h init_gen_rand"
    c_init_gen_rand         :: UInt32 -> IO ()

-- Getting a random int

-- This function generates and returns 64-bit pseudorandom number.
-- init_gen_rand or init_by_array must be called before this function.
-- The function gen_rand64 should not be called after gen_rand32,
-- unless an initialization is again executed. 

#if WORD_SIZE_IN_BITS < 64

foreign import ccall unsafe "SFMT.h gen_rand32"
    c_gen_rand32            :: IO UInt32

foreign import ccall unsafe "SFMT_wrap.h gen_rand64_mix_wrap"
    c_gen_rand64_mix        :: IO UInt64

#else

foreign import ccall unsafe "SFMT.h gen_rand64"
    c_gen_rand64            :: IO UInt64

#endif

-- Getting a random double

-- | Generates a random number on [0,1)-real-interval
-- calls gen_rand32
-- | Generates a random number on [0,1) with 53-bit resolution. Fast on 64 bit machines.
-- calls gen_rand64

-- | generates a random number on [0,1) with 53-bit resolution using
-- 32bit integer

#if WORD_SIZE_IN_BITS < 64
foreign import ccall unsafe "SFMT_wrap.h genrand_real2_wrap"
    c_genrand_real2         :: IO CDouble

-- foreign import ccall unsafe "SFMT.h genrand_res53_mix"
--     c_genrand_res53_mix     :: IO CDouble
#else
foreign import ccall unsafe "SFMT_wrap.h genrand_res53_wrap"
    c_genrand_res53         :: IO CDouble
#endif

------------------------------------------------------------------------

{-
-- Generates a random number on [0,1]-real-interval
-- calls gen_rand32
foreign import ccall unsafe "SFMT.h genrand_real1"
    c_genrand_real1         :: IO CDouble

-- | Generates a random number on (0,1)-real-interval
-- calls gen_rand32
foreign import ccall unsafe "SFMT.h genrand_real3"
    c_genrand_real3         :: IO CDouble
-}

------------------------------------------------------------------------

{-
foreign import ccall unsafe "SFMT.h get_min_array_size32"
    c_get_min_array_size32  :: IO CInt

foreign import ccall unsafe "SFMT.h get_min_array_size64"
    c_get_min_array_size64  :: IO CInt

foreign import ccall unsafe "SFMT.h fill_array32"
    c_fill_array32          :: Ptr UInt32 -> CInt -> IO ()

foreign import ccall unsafe "SFMT.h fill_array64"
    c_fill_array64          :: Ptr UInt64 -> CInt -> IO ()
-}

------------------------------------------------------------------------

foreign import ccall unsafe "SFMT.h get_idstring"
    c_get_idstring          :: IO CString

foreign import ccall unsafe "SFMT.h get_initialized"
    c_get_initialized        :: IO CInt

--
-- Invariant: we can never call rand32 if we're in 64 bit land,
-- and never call rand64 if in 32 bit land.
-- 
-- audit this!
--