packages feed

mersenne-random-pure64 (empty) → 0.1

raw patch · 13 files changed

+946/−0 lines, 13 filesdep +basedep +old-timedep +randomsetup-changed

Dependencies added: base, old-time, random

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) Don Stewart <dons@galois.com> 2008++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ LICENSE.mt19937-64 view
@@ -0,0 +1,54 @@+/* +   A C-program for MT19937-64 (2004/9/29 version).+   Coded by Takuji Nishimura and Makoto Matsumoto.++   This is a 64-bit version of Mersenne Twister pseudorandom number+   generator.++   Before using, initialize the state by using init_genrand64(seed)  +   or init_by_array64(init_key, key_length).++   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,+   All rights reserved.                          ++   Redistribution and use in source and binary forms, with or without+   modification, are permitted provided that the following conditions+   are met:++     1. Redistributions of source code must retain the above copyright+        notice, this list of conditions and the following disclaimer.++     2. Redistributions in binary form must reproduce the above copyright+        notice, this list of conditions and the following disclaimer in the+        documentation and/or other materials provided with the distribution.++     3. The names of its contributors may not be used to endorse or promote +        products derived from this software without specific prior written +        permission.++   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.++   References:+   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''+     ACM Transactions on Modeling and +     Computer Simulation 10. (2000) 348--357.+   M. Matsumoto and T. Nishimura,+     ``Mersenne Twister: a 623-dimensionally equidistributed+       uniform pseudorandom number generator''+     ACM Transactions on Modeling and +     Computer Simulation 8. (Jan. 1998) 3--30.++   Any feedback is very welcome.+   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)+*/
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ System/Random/Mersenne/Pure64.hsc view
@@ -0,0 +1,198 @@+{-# LANGUAGE CPP, ForeignFunctionInterface #-}+--------------------------------------------------------------------+-- |+-- Module     : System.Random.Mersenne.Pure64+-- 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+--+-- A purely functional binding 64 bit binding to the classic mersenne+-- twister random number generator. This is more flexible than the +-- impure 'mersenne-random' library, at the cost of being much slower.+-- This generator is however, many times faster than System.Random,+-- and yields high quality randoms with a long period.+--+-- This generator may be used with System.Random, however, that is+-- likely to be slower than using it directly.+--+module System.Random.Mersenne.Pure64 (++    -- * The random number generator+    PureMT          -- abstract: RandomGen++    -- * Introduction+    , pureMT        -- :: Word -> PureMT+    , newPureMT     -- :: IO PureMT++    -- $instance++    -- * Low level access to the generator++    -- $notes+    , randomInt     -- :: PureMT -> (Int   ,PureMT)+    , randomWord    -- :: PureMT -> (Word  ,PureMT)+    , randomInt64   -- :: PureMT -> (Int64 ,PureMT)+    , randomWord64  -- :: PureMT -> (Word64,PureMT)+    , randomDouble  -- :: PureMT -> (Double,PureMT)++    ) where++------------------------------------------------------------------------++#include "mt19937-64.h"++------------------------------------------------------------------------++import System.Random.Mersenne.Pure64.Base++#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 605+import GHC.ForeignPtr           (mallocPlainForeignPtrBytes)+#else+import Foreign.ForeignPtr       (mallocForeignPtrBytes)+#endif++import Foreign.C.Types+import Foreign++import System.CPUTime	( getCPUTime )+import System.Time+import System.IO.Unsafe+import Control.Monad++import System.Random++import Data.Word+import Data.Int++------------------------------------------------------------------------++-- | Create a PureMT generator from a 'Word' seed.+pureMT :: Word -> PureMT+pureMT = init_genrand64 . fromIntegral++-- | Create a new PureMT generator, using the clocktime as the base for the seed.+newPureMT :: IO PureMT+newPureMT = do+    ct             <- getCPUTime+    (TOD sec psec) <- getClockTime+    return $ pureMT (fromIntegral $ sec * 1013904242 + psec + ct)++------------------------------------------------------------------------+-- System.Random interface.++-- $instance +-- +-- Being purely functional, the PureMT generator is an instance of+-- RandomGen. However, it doesn't support 'split' yet.++instance RandomGen PureMT where+   next  = randomInt+   split = error "System.Random.Mersenne.Pure: unable to split the mersenne twister"++------------------------------------------------------------------------+-- Direct access to Int, Word and Double types++-- | Yield a new 'Int' value from the generator, returning a new+-- generator and that 'Int'. The full 64 bits will be used on a 64 bit machine.+randomInt :: PureMT -> (Int,PureMT)+randomInt g = (fromIntegral i, g')+        where (i,g') = genrand64_int64 g++-- | Yield a new 'Word' value from the generator, returning a new+-- generator and that 'Word'. +randomWord :: PureMT -> (Word,PureMT)+randomWord g = (fromIntegral i, g')+        where (i,g') = genrand64_int64 g++-- | Yield a new 'Int64' value from the generator, returning a new+-- generator and that 'Int64'. +randomInt64 :: PureMT -> (Int64,PureMT)+randomInt64 g = (fromIntegral i, g')+        where (i,g') = genrand64_int64 g++-- | Yield a new 'Word64' value from the generator, returning a new+-- generator and that 'Word64'. +randomWord64 :: PureMT -> (Word64,PureMT)+randomWord64 g = (fromIntegral i, g')+        where (i,g') = genrand64_int64 g++-- | Efficiently yield a new 53-bit precise 'Double' value, and a new generator.+randomDouble :: PureMT -> (Double,PureMT)+randomDouble g = (realToFrac i, g')+        where (i,g') = genrand64_real2 g++------------------------------------------------------------------------+-- 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.+--++-- | 'PureMT', a pure mersenne twister pseudo-random number generator+--+newtype PureMT  = PureMT (ForeignPtr MTState)++instance Show PureMT where+    show _ = show "<PureMT>"++sizeof_MTState :: Int+sizeof_MTState = (#const (sizeof (struct mt_state_t))) -- 2504 bytes++------------------------------------------------------------------------+-- Low level interface++init_genrand64 :: UInt64 -> PureMT+init_genrand64 seed = unsafePerformIO $ do+    fp <- mallocFastBytes sizeof_MTState+    withForeignPtr fp $ \p -> c_init_genrand64 p seed -- fill it+    return $ PureMT fp++genrand64_int64 :: PureMT -> (UInt64, PureMT)+genrand64_int64 o   = unsafePerformIO $ do+    PureMT n <- copyPureMT o+    v        <- withForeignPtr n $ c_genrand64_int64+    n `seq` v `seq`+        return (v,PureMT n)+{-# INLINE genrand64_int64 #-}++genrand64_real2 :: PureMT -> (CDouble, PureMT)+genrand64_real2 o   = unsafePerformIO $ do+    PureMT n <- copyPureMT o+    v        <- withForeignPtr n $ c_genrand64_real2+    n `seq` v `seq`+        return (v,PureMT n)+{-# INLINE genrand64_real2 #-}++-- ---------------------------------------------------------------------+-- Memory management++--+-- | Wrapper of mallocForeignPtrBytes with faster implementation+-- Allocate pinned heap memory for the state.+--+mallocFastBytes :: Int -> IO (ForeignPtr a)+mallocFastBytes s = do+#if __GLASGOW_HASKELL__ >= 605 && !defined(SLOW_FOREIGN_PTR)+    mallocPlainForeignPtrBytes s+#else+    mallocForeignPtrBytes s+#endif+{-# INLINE mallocFastBytes #-}++--+-- Duplicate a random state, prior to permuting it.+--+copyPureMT :: PureMT -> IO PureMT+copyPureMT (PureMT oldfp) = do+    newfp <- mallocFastBytes sizeof_MTState+    withForeignPtr newfp $ \p ->+        withForeignPtr oldfp $ \q ->+            c_memcpy (castPtr p) (castPtr q) (fromIntegral sizeof_MTState)+    return $! PureMT newfp+{-# INLINE copyPureMT #-}++
+ System/Random/Mersenne/Pure64/Base.hsc view
@@ -0,0 +1,55 @@+{-# LANGUAGE EmptyDataDecls, CPP, ForeignFunctionInterface #-}+--------------------------------------------------------------------+-- |+-- Module     : System.Random.Mersenne.Pure64+-- Copyright  : Copyright (c) 2008, Don Stewart <dons@galois.com>+-- License    : BSD3+-- Maintainer : Don Stewart <dons@galois.com>+-- Stability  : experimental+-- Portability: CPP, FFI, EmptyDataDecls+-- Tested with: GHC 6.8.2+--+-- A purely functional binding 64 bit binding to the classic mersenne+-- twister random number generator. This is more flexible than the +-- impure 'mersenne-random' library, at the cost of being much slower.+-- This generator is however, many times faster than System.Random,+-- and yields high quality randoms with a long period.+--+module System.Random.Mersenne.Pure64.Base where++#include "mt19937-64.h"+#include "mt19937-64-unsafe.h"++import Foreign.C.Types+import Foreign++data MTState++type UInt64 = CULLong++------------------------------------------------------------------------+-- pure version:++foreign import ccall unsafe "init_genrand64"+    c_init_genrand64  :: Ptr MTState -> UInt64 -> IO ()++foreign import ccall unsafe "genrand64_int64"+    c_genrand64_int64 :: Ptr MTState -> IO UInt64++foreign import ccall unsafe "genrand64_real2"+    c_genrand64_real2 :: Ptr MTState -> IO CDouble++------------------------------------------------------------------------+-- model: (for testing purposes)++foreign import ccall unsafe "init_genrand64_unsafe"+    c_init_genrand64_unsafe :: UInt64 -> IO ()++foreign import ccall unsafe "genrand64_int64_unsafe"+    c_genrand64_int64_unsafe :: IO UInt64++foreign import ccall unsafe "genrand64_real2_unsafe"+    c_genrand64_real2_unsafe :: IO CDouble++foreign import ccall unsafe "string.h memcpy"+    c_memcpy :: Ptr Word8 -> Ptr Word8 -> CSize -> IO (Ptr Word8)
+ TODO view
@@ -0,0 +1,5 @@+-fvia-C seems to break things!+Check sequences are repeatable+check behaviour on 32 bits+profile/benchmark+tests driver
+ cbits/mt19937-64-unsafe.c view
@@ -0,0 +1,122 @@+/* +   A C-program for MT19937-64 (2004/9/29 version).+   Coded by Takuji Nishimura and Makoto Matsumoto.++   This is a 64-bit version of Mersenne Twister pseudorandom number+   generator.++   Before using, initialize the state by using init_genrand64(seed)  +   or init_by_array64(init_key, key_length).++   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,+   All rights reserved.                          ++   Redistribution and use in source and binary forms, with or without+   modification, are permitted provided that the following conditions+   are met:++     1. Redistributions of source code must retain the above copyright+        notice, this list of conditions and the following disclaimer.++     2. Redistributions in binary form must reproduce the above copyright+        notice, this list of conditions and the following disclaimer in the+        documentation and/or other materials provided with the distribution.++     3. The names of its contributors may not be used to endorse or promote +        products derived from this software without specific prior written +        permission.++   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.++   References:+   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''+     ACM Transactions on Modeling and +     Computer Simulation 10. (2000) 348--357.+   M. Matsumoto and T. Nishimura,+     ``Mersenne Twister: a 623-dimensionally equidistributed+       uniform pseudorandom number generator''+     ACM Transactions on Modeling and +     Computer Simulation 8. (Jan. 1998) 3--30.++   Any feedback is very welcome.+   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)+*/+++#include <stdio.h>+#include "mt19937-64-unsafe.h"++#define NN 312+#define MM 156+#define MATRIX_A 0xB5026F5AA96619E9ULL+#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */+#define LM 0x7FFFFFFFULL /* Least significant 31 bits */+++/* The array for the state vector */+static unsigned long long mt[NN]; +/* mti==NN+1 means mt[NN] is not initialized */+static int mti=NN+1; ++/* initializes mt[NN] with a seed */+void init_genrand64_unsafe(unsigned long long seed)+{+    mt[0] = seed;+    for (mti=1; mti<NN; mti++) +        mt[mti] =  (6364136223846793005ULL * (mt[mti-1] ^ (mt[mti-1] >> 62)) + mti);+}++/* generates a random number on [0, 2^64-1]-interval */+unsigned long long genrand64_int64_unsafe(void)+{+    int i;+    unsigned long long x;+    static unsigned long long mag01[2]={0ULL, MATRIX_A};++    if (mti >= NN) { /* generate NN words at one time */++        /* if init_genrand64() has not been called, */+        /* a default initial seed is used     */+        if (mti == NN+1) +            init_genrand64_unsafe(5489ULL); ++        for (i=0;i<NN-MM;i++) {+            x = (mt[i]&UM)|(mt[i+1]&LM);+            mt[i] = mt[i+MM] ^ (x>>1) ^ mag01[(int)(x&1ULL)];+        }+        for (;i<NN-1;i++) {+            x = (mt[i]&UM)|(mt[i+1]&LM);+            mt[i] = mt[i+(MM-NN)] ^ (x>>1) ^ mag01[(int)(x&1ULL)];+        }+        x = (mt[NN-1]&UM)|(mt[0]&LM);+        mt[NN-1] = mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)];++        mti = 0;+    }+  +    x = mt[mti++];++    x ^= (x >> 29) & 0x5555555555555555ULL;+    x ^= (x << 17) & 0x71D67FFFEDA60000ULL;+    x ^= (x << 37) & 0xFFF7EEE000000000ULL;+    x ^= (x >> 43);++    return x;+}++/* generates a random number on [0,1)-real-interval */+double genrand64_real2_unsafe(void)+{+    return (genrand64_int64_unsafe() >> 11) * (1.0/9007199254740992.0);+}
+ cbits/mt19937-64.c view
@@ -0,0 +1,116 @@+/* +   A C-program for MT19937-64 (2004/9/29 version).+   Coded by Takuji Nishimura and Makoto Matsumoto.++   This is a 64-bit version of Mersenne Twister pseudorandom number+   generator.++   Before using, initialize the state by using init_genrand64(seed)  +   or init_by_array64(init_key, key_length).++   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,+   All rights reserved.                          ++   Redistribution and use in source and binary forms, with or without+   modification, are permitted provided that the following conditions+   are met:++     1. Redistributions of source code must retain the above copyright+        notice, this list of conditions and the following disclaimer.++     2. Redistributions in binary form must reproduce the above copyright+        notice, this list of conditions and the following disclaimer in the+        documentation and/or other materials provided with the distribution.++     3. The names of its contributors may not be used to endorse or promote +        products derived from this software without specific prior written +        permission.++   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.++   References:+   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''+     ACM Transactions on Modeling and +     Computer Simulation 10. (2000) 348--357.+   M. Matsumoto and T. Nishimura,+     ``Mersenne Twister: a 623-dimensionally equidistributed+       uniform pseudorandom number generator''+     ACM Transactions on Modeling and +     Computer Simulation 8. (Jan. 1998) 3--30.++   Any feedback is very welcome.+   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)+*/+++#include <stdio.h>+#include <stdlib.h>+#include "mt19937-64.h"++/*+ * When passed an empty state, initialize the states' mt[NN] with a seed+ */+void init_genrand64(mt_state st, unsigned long long seed)+{+    st->mti=NN+1;++    st->mt[0] = seed;+    for (st->mti=1; st->mti<NN; st->mti++) +        st->mt[st->mti] =  (6364136223846793005ULL * +                                (st->mt[st->mti-1] ^ (st->mt[st->mti-1] >> 62)) + st->mti);+}++/* generates a random number on [0, 2^64-1]-interval */+unsigned long long genrand64_int64(mt_state st)+{+    int i;+    unsigned long long x;+    static unsigned long long mag01[2]={0ULL, MATRIX_A};++    if (st->mti >= NN) { /* generate NN words at one time */++        /* if init_genrand64() has not been called, */+        /* a default initial seed is used     */+        if (st->mti == NN+1) +            init_genrand64(st, 5489ULL); ++        for (i=0;i<NN-MM;i++) {+            x = (st->mt[i]&UM)|(st->mt[i+1]&LM);+            st->mt[i] = st->mt[i+MM] ^ (x>>1) ^ mag01[(int)(x&1ULL)];+        }+        for (;i<NN-1;i++) {+            x = (st->mt[i]&UM)|(st->mt[i+1]&LM);+            st->mt[i] = st->mt[i+(MM-NN)] ^ (x>>1) ^ mag01[(int)(x&1ULL)];+        }+        x = (st->mt[NN-1]&UM)|(st->mt[0]&LM);+        st->mt[NN-1] = st->mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)];++        st->mti = 0;+    }+  +    x = st->mt[st->mti++];++    x ^= (x >> 29) & 0x5555555555555555ULL;+    x ^= (x << 17) & 0x71D67FFFEDA60000ULL;+    x ^= (x << 37) & 0xFFF7EEE000000000ULL;+    x ^= (x >> 43);++    return x;+}++/* generates a random number on [0,1)-real-interval */+double genrand64_real2(mt_state st)+{+    return (genrand64_int64(st) >> 11) * (1.0/9007199254740992.0);+}
+ include/mt19937-64-unsafe.h view
@@ -0,0 +1,64 @@+/* +   A C-program for MT19937-64 (2004/9/29 version).+   Coded by Takuji Nishimura and Makoto Matsumoto.++   This is a 64-bit version of Mersenne Twister pseudorandom number+   generator.++   Before using, initialize the state by using init_genrand64(seed)  +   or init_by_array64(init_key, key_length).++   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,+   All rights reserved.                          ++   Redistribution and use in source and binary forms, with or without+   modification, are permitted provided that the following conditions+   are met:++     1. Redistributions of source code must retain the above copyright+        notice, this list of conditions and the following disclaimer.++     2. Redistributions in binary form must reproduce the above copyright+        notice, this list of conditions and the following disclaimer in the+        documentation and/or other materials provided with the distribution.++     3. The names of its contributors may not be used to endorse or promote +        products derived from this software without specific prior written +        permission.++   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.++   References:+   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''+     ACM Transactions on Modeling and +     Computer Simulation 10. (2000) 348--357.+   M. Matsumoto and T. Nishimura,+     ``Mersenne Twister: a 623-dimensionally equidistributed+       uniform pseudorandom number generator''+     ACM Transactions on Modeling and +     Computer Simulation 8. (Jan. 1998) 3--30.++   Any feedback is very welcome.+   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)+*/+++/* initializes mt[NN] with a seed */+void init_genrand64_unsafe(unsigned long long seed);++/* generates a random number on [0, 2^64-1]-interval */+unsigned long long genrand64_int64_unsafe(void);++/* generates a random number on [0,1)-real-interval */+double genrand64_real2_unsafe(void);
+ include/mt19937-64.h view
@@ -0,0 +1,76 @@+/* +   A C-program for MT19937-64 (2004/9/29 version).+   Coded by Takuji Nishimura and Makoto Matsumoto.++   This is a 64-bit version of Mersenne Twister pseudorandom number+   generator.++   Before using, initialize the state by using init_genrand64(seed)  +   or init_by_array64(init_key, key_length).++   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,+   All rights reserved.                          ++   Redistribution and use in source and binary forms, with or without+   modification, are permitted provided that the following conditions+   are met:++     1. Redistributions of source code must retain the above copyright+        notice, this list of conditions and the following disclaimer.++     2. Redistributions in binary form must reproduce the above copyright+        notice, this list of conditions and the following disclaimer in the+        documentation and/or other materials provided with the distribution.++     3. The names of its contributors may not be used to endorse or promote +        products derived from this software without specific prior written +        permission.++   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.++   References:+   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''+     ACM Transactions on Modeling and +     Computer Simulation 10. (2000) 348--357.+   M. Matsumoto and T. Nishimura,+     ``Mersenne Twister: a 623-dimensionally equidistributed+       uniform pseudorandom number generator''+     ACM Transactions on Modeling and +     Computer Simulation 8. (Jan. 1998) 3--30.++   Any feedback is very welcome.+   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)+*/+++#define NN 312+#define MM 156+#define MATRIX_A 0xB5026F5AA96619E9ULL+#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */+#define LM 0x7FFFFFFFULL /* Least significant 31 bits */++/* The array for the state vector */+typedef struct mt_state_t {+    unsigned long long mt[NN]; +    int mti;+} mt_state_struct;++typedef struct mt_state_t *mt_state;++/* initializes mt[NN] with a seed */+void init_genrand64(mt_state st, unsigned long long seed);++/* generates a random number on [0, 2^64-1]-interval */+unsigned long long genrand64_int64(mt_state st);+double genrand64_real2(mt_state st);
+ mersenne-random-pure64.cabal view
@@ -0,0 +1,57 @@+name:            mersenne-random-pure64+version:         0.1+homepage:        http://code.haskell.org/~dons/code/mt19937-random+synopsis:        Generate high quality pseudorandom numbers purely using a Mersenne Twister+description:+    The Mersenne twister is a pseudorandom number generator developed by+    Makoto Matsumoto and Takuji Nishimura that is based on a matrix linear+    recurrence over a finite binary field. It provides for fast generation+    of very high quality pseudorandom numbers. The source for the C code+    can be found here:+    .+    <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html>+    .+    This library provides a purely functional binding to the 64 bit+    classic mersenne twister, along with instances of RandomGen, so the+    generator can be used with System.Random. The generator should+    typically be a few times faster than the default StdGen (but much+    slower than the impure 'mersenne-random' library based on SIMD+    instructions and destructive state updates.+    .+category:        Math, System+license:         BSD3+license-file:    LICENSE+copyright:       (c) 2008. Don Stewart <dons@galois.com>+author:          Don Stewart+maintainer:      Don Stewart <dons@galois.com>+cabal-version: >= 1.2.0+build-type:      Configure+tested-with:     GHC ==6.8.2+build-type:      Configure++flag small_base+  description: Build with new smaller base library+  default:     False++library+    exposed-modules: System.Random.Mersenne.Pure64+                     System.Random.Mersenne.Pure64.Base+    extensions:      CPP, ForeignFunctionInterface++    if flag(small_base)+        build-depends: base  < 3+    else+        build-depends: base >= 3, old-time, random++    cc-options:+        -O3 -finline-functions -fomit-frame-pointer+        -fno-strict-aliasing --param max-inline-insns-single=1800++    ghc-options:     -Wall -O2 -fvia-C -fexcess-precision++    c-sources:        cbits/mt19937-64.c+                      cbits/mt19937-64-unsafe.c+    include-dirs:     include+    includes:+    install-includes: mt19937-64.h+                      mt19937-64-unsafe.h
+ tests/Unit.hs view
@@ -0,0 +1,151 @@+{-# OPTIONS -fbang-patterns -fglasgow-exts #-}++import Control.Exception+import Control.Monad+import Data.Int+import Data.Typeable+import Data.Word+import System.CPUTime+import System.Environment+import System.IO+import Text.Printf+import qualified System.Random as Old+import qualified System.Random.Mersenne as Unsafe++import System.Random.Mersenne.Pure64+import System.Random.Mersenne.Pure64.Base+import Control.Concurrent+import Control.Concurrent.MVar++time :: IO t -> IO t+time a = do+    start <- getCPUTime+    v <- a+    end   <- getCPUTime+    let diff = (fromIntegral (end - start)) / (10^12)+    printf "Computation time: %0.3f sec\n" (diff :: Double)+    return v++seed = 7++main = do++    c_init_genrand64_unsafe seed+    let g = pureMT (fromIntegral seed)++    ------------------------------------------------------------------------+    -- calibrate+    s <- newMVar 0 :: IO (MVar Int)+    putStr "Callibrating ... " >> hFlush stdout++    tid <- forkIO $ do+        let go !i !g = do+                let (!_,!g') = randomWord64 g+                x <- swapMVar s i+                x `seq` go (i+1) g'+        go 0 g++    threadDelay (1000 * 1000)+    killThread tid+    lim <- readMVar s -- 1 sec worth of generation+    putStrLn $ "done. Using N="++ show lim++    time $ do+        let m = 2*lim+        putStr $ "Checking against released mt19937-64.c to depth " ++ show m ++ " "+        hFlush stdout+        equivalent g m++    speed lim ++    return ()++------------------------------------------------------------------------++equivalent !g !n | n > 0 = do++    i'      <- c_genrand64_int64_unsafe+    d'      <- c_genrand64_real2_unsafe++    let (i,g')  = randomWord64 g+        (d,g'') = randomDouble g'++    if i == fromIntegral i' && d == realToFrac d'+        then do when (n `rem` 500000 == 0) $ putChar '.' >> hFlush stdout+                equivalent g'' (n-1)++        else do print $ "Failed! " ++ show ((i,i') , (d,d'))+                return g''++equivalent g _ = do putStrLn "Matches model!"+                    return g++------------------------------------------------------------------------+-- compare with System.Random++-- overhead cause by random's badness+speed lim = do++ time $ do+    putStrLn $ "System.Random"+    let g = Old.mkStdGen 5+    let go :: Old.StdGen -> Int -> Int -> Int+        go !g !n !acc+            | n >= lim = acc+            | otherwise     =+                    let (a,g') = Old.random g+                    in go g' (n+1) (if a > acc then a else acc)+    print (go g 0 0)++ time $ do+    putStrLn $ "System.Random with our generator"+    let g = pureMT 5+    let go :: PureMT -> Int -> Int -> Int+        go !g !n !acc+            | n >= lim = acc+            | otherwise     =+                    let (a,g') = Old.random g+                    in go g' (n+1) (if a > acc then a else acc)+    print (go g 0 0)++ time $ do+    putStrLn $ "System.Random.Mersenne.Pure"+    let g = pureMT 5+    let go :: PureMT -> Int -> Int -> Int+        go !g !n !acc+            | n >= lim = acc+            | otherwise     =+                    let (a',g') = randomWord64 g+                        a = fromIntegral a'+                    in go g' (n+1) (if a > acc then a else acc)+    print (go g 0 0)++ time $ do+    putStrLn $ "System.Random.Mersenne.Pure (unique state)"+    c_init_genrand64_unsafe 5+    let go :: Int -> Int -> IO Int+        go !n !acc+            | n >= lim = return acc+            | otherwise     = do+                    a' <- c_genrand64_int64_unsafe+                    let a = fromIntegral a'+                    go (n+1::Int) (if a > acc then a else acc)+    print =<< go 0 0++ time $ do+    putStrLn $ "System.Random.Mersenne.Unsafe"+    g <- Unsafe.newMTGen (Just 5)++    let go :: Int -> Int -> IO Int+        go !n !acc+            | n >= lim = return acc+            | otherwise     = do+                    a <- Unsafe.random g+                    go (n+1::Int) (if a > acc then a else acc)++    print =<< go 0 0+++--    printf "MT is %s times faster generating %s\n" (show $x`div`y) (show (typeOf ty))+--    return ()+
+ tests/copy.hs view
@@ -0,0 +1,15 @@+{-# LANGUAGE BangPatterns #-}++import System.Random.Mersenne.Pure64++main = do++    let g = pureMT 7++    let go :: Int -> PureMT -> IO ()+        go 0 !p = return ()+        go n !p  = do+            let (_,q) = randomWord p+            go (n-1) q++    go 10000000 g -- 10s constant space