mersenne-random (empty) → 0.1

raw patch · 28 files changed

+2709/−0 lines, 28 filesdep +basedep +old-timesetup-changed

Dependencies added: base, old-time

Files

LICENSE +30/−0
LICENSE.sfmt +29/−0
README +5/−0
Setup.lhs +3/−0
System/Random/Mersenne.hs +512/−0
TODO +27/−0
cbits/SFMT.c +628/−0
include/SFMT-alti.h +156/−0
include/SFMT-params.h +97/−0
include/SFMT-params11213.h +48/−0
include/SFMT-params1279.h +46/−0
include/SFMT-params132049.h +46/−0
include/SFMT-params19937.h +46/−0
include/SFMT-params216091.h +46/−0
include/SFMT-params2281.h +46/−0
include/SFMT-params4253.h +46/−0
include/SFMT-params44497.h +46/−0
include/SFMT-params607.h +46/−0
include/SFMT-params86243.h +46/−0
include/SFMT-sse2.h +121/−0
include/SFMT.h +168/−0
mersenne-random.cabal +109/−0
tests/Makefile +11/−0
tests/Unit.hs +264/−0
tests/number-fast.hs +13/−0
tests/pi.hs +52/−0
tests/runtests.sh +5/−0
tests/sum.c +17/−0

+ 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.sfmt view

@@ -0,0 +1,29 @@+Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima+University. All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are+met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.+    * 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.+    * Neither the name of the Hiroshima University nor the names of+      its contributors may 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.

+ README view

@@ -0,0 +1,5 @@+To build with sse2 support (core2 and other intels)++    $ runhaskell Setup.lhs -f use_sse2++Similar flags for altivec on the mac

+ Setup.lhs view

@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain

+ System/Random/Mersenne.hs view

@@ -0,0 +1,512 @@+{-# 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(..)++    -- * 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++------------------------------------------------------------------------+-- $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)+-}++------------------------------------------------------------------------+-- Efficient basic instances++instance MTRandom Word    where+    random !_          = randomWord+--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Word64  where+    random !_ = fmap fromIntegral randomWord64+--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Word32  where+    random !_ = fmap fromIntegral randomWord+--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Word16  where+    random !_ = fmap fromIntegral randomWord+--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Word8   where+    random !_          = fmap fromIntegral randomWord+--    randomR (lo,hi) !g = randomIvalWord g (fromIntegral lo) (fromIntegral hi)++------------------------------------------------------------------------++instance MTRandom Double where+    random !_ = randomDouble+--    randomR (lo,hi) g = randomIvalDouble g lo hi id++------------------------------------------------------------------------++instance MTRandom Int     where+    random !_          = randomInt+--    randomR (lo,hi) !g = randomIvalInt g lo hi++instance MTRandom Int64   where+    random !_ = fmap fromIntegral randomInt64+--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Int32   where+    random !_ = fmap fromIntegral randomInt+--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Int16   where+    random !_ = fmap fromIntegral randomInt+--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Int8    where+    random !_ = fmap fromIntegral randomInt+--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)++instance MTRandom Integer where+    random !_ = fmap fromIntegral randomInt+--    randomR (lo,hi) !g = randomIvalInt g (fromIntegral lo) (fromIntegral hi)++------------------------------------------------------------------------++instance MTRandom Bool    where+    random !_ = do x <- randomWord; return $! x .&. 1 /= 0++{-+    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)+-}++------------------------------------------------------------------------++-- | 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.h gen_rand64_mix"+    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.h genrand_real2"+    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.h genrand_res53"+    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!+--

+ TODO view

@@ -0,0 +1,27 @@++-- support mulitple generators++    Want to save the state, and later continue the generation.+    In the case of mt19937ar.c, if one saves mt[N] (624 words) and the+    counter mti (1 word), then one can continue the computation by loading+    these values.++    * too much copying...+    * supporting multiple, impure generators is possible though.+++-- Also provide a RandomGen instance++-- eager evaluation+    -- fill arrays in chunks+    -- possible api for this?++-- support fast range bounds+    -- too slow on haskell side currently.++-- statistical correctness, DieHard suite.+    -- too complex.++-- coverage tests+    -- not working+

+ cbits/SFMT.c view

@@ -0,0 +1,628 @@+/** + * @file  SFMT.c+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)+ *+ * @author Mutsuo Saito (Hiroshima University)+ * @author Makoto Matsumoto (Hiroshima University)+ *+ * Copyright (C) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima+ * University. All rights reserved.+ *+ * The new BSD License is applied to this software, see LICENSE.txt+ */+#include <string.h>+#include <assert.h>+#include "SFMT.h"+#include "SFMT-params.h"++#if defined(__BIG_ENDIAN__) && !defined(__amd64) && !defined(BIG_ENDIAN64)+#define BIG_ENDIAN64 1+#endif+#if defined(HAVE_ALTIVEC) && !defined(BIG_ENDIAN64)+#define BIG_ENDIAN64 1+#endif+#if defined(ONLY64) && !defined(BIG_ENDIAN64)+  #if defined(__GNUC__)+    #error "-DONLY64 must be specified with -DBIG_ENDIAN64"+  #endif+#undef ONLY64+#endif+/*------------------------------------------------------+  128-bit SIMD data type for Altivec, SSE2 or standard C+  ------------------------------------------------------*/+#if defined(HAVE_ALTIVEC)+  #if !defined(__APPLE__)+    #include <altivec.h>+  #endif+/** 128-bit data structure */+union W128_T {+    vector unsigned int s;+    uint32_t u[4];+};+/** 128-bit data type */+typedef union W128_T w128_t;++#elif defined(HAVE_SSE2)+  #include <emmintrin.h>++/** 128-bit data structure */+union W128_T {+    __m128i si;+    uint32_t u[4];+};+/** 128-bit data type */+typedef union W128_T w128_t;++#else++/** 128-bit data structure */+struct W128_T {+    uint32_t u[4];+};+/** 128-bit data type */+typedef struct W128_T w128_t;++#endif++/*--------------------------------------+  FILE GLOBAL VARIABLES+  internal state, index counter and flag +  --------------------------------------*/++/** the 128-bit internal state array */+static w128_t sfmt[N];++/** the 32bit integer pointer to the 128-bit internal state array */+static uint32_t *psfmt32 = &sfmt[0].u[0];++#if !defined(BIG_ENDIAN64) || defined(ONLY64)+/** the 64bit integer pointer to the 128-bit internal state array */+static uint64_t *psfmt64 = (uint64_t *)&sfmt[0].u[0];+#endif++/** index counter to the 32-bit internal state array */+static int idx;++/** a flag: it is 0 if and only if the internal state is not yet+ * initialized. */+static int initialized = 0;++/** UNUSED: a parity check vector which certificate the period of 2^{MEXP} */+static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4};++/*----------------+  STATIC FUNCTIONS+  ----------------*/+inline static int idxof(int i);+inline static void rshift128(w128_t *out,  w128_t const *in, int shift);+inline static void lshift128(w128_t *out,  w128_t const *in, int shift);+inline static void gen_rand_all(void);+inline static void gen_rand_array(w128_t *array, int size);+inline static uint32_t func1(uint32_t x);+inline static uint32_t func2(uint32_t x);+static void period_certification(void);+#if defined(BIG_ENDIAN64) && !defined(ONLY64)+inline static void swap(w128_t *array, int size);+#endif++#if defined(HAVE_ALTIVEC)+  #include "SFMT-alti.h"+#elif defined(HAVE_SSE2)+  #include "SFMT-sse2.h"+#endif++/**+ * This function simulate a 64-bit index of LITTLE ENDIAN + * in BIG ENDIAN machine.+ */+#ifdef ONLY64+inline static int idxof(int i) {+    return i ^ 1;+}+#else+inline static int idxof(int i) {+    return i;+}+#endif+/**+ * This function simulates SIMD 128-bit right shift by the standard C.+ * The 128-bit integer given in in is shifted by (shift * 8) bits.+ * This function simulates the LITTLE ENDIAN SIMD.+ * @param out the output of this function+ * @param in the 128-bit data to be shifted+ * @param shift the shift value+ */+#ifdef ONLY64+inline static void rshift128(w128_t *out, w128_t const *in, int shift) {+    uint64_t th, tl, oh, ol;++    th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);+    tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);++    oh = th >> (shift * 8);+    ol = tl >> (shift * 8);+    ol |= th << (64 - shift * 8);+    out->u[0] = (uint32_t)(ol >> 32);+    out->u[1] = (uint32_t)ol;+    out->u[2] = (uint32_t)(oh >> 32);+    out->u[3] = (uint32_t)oh;+}+#else+inline static void rshift128(w128_t *out, w128_t const *in, int shift) {+    uint64_t th, tl, oh, ol;++    th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);+    tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);++    oh = th >> (shift * 8);+    ol = tl >> (shift * 8);+    ol |= th << (64 - shift * 8);+    out->u[1] = (uint32_t)(ol >> 32);+    out->u[0] = (uint32_t)ol;+    out->u[3] = (uint32_t)(oh >> 32);+    out->u[2] = (uint32_t)oh;+}+#endif+/**+ * This function simulates SIMD 128-bit left shift by the standard C.+ * The 128-bit integer given in in is shifted by (shift * 8) bits.+ * This function simulates the LITTLE ENDIAN SIMD.+ * @param out the output of this function+ * @param in the 128-bit data to be shifted+ * @param shift the shift value+ */+#ifdef ONLY64+inline static void lshift128(w128_t *out, w128_t const *in, int shift) {+    uint64_t th, tl, oh, ol;++    th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);+    tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);++    oh = th << (shift * 8);+    ol = tl << (shift * 8);+    oh |= tl >> (64 - shift * 8);+    out->u[0] = (uint32_t)(ol >> 32);+    out->u[1] = (uint32_t)ol;+    out->u[2] = (uint32_t)(oh >> 32);+    out->u[3] = (uint32_t)oh;+}+#else+inline static void lshift128(w128_t *out, w128_t const *in, int shift) {+    uint64_t th, tl, oh, ol;++    th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);+    tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);++    oh = th << (shift * 8);+    ol = tl << (shift * 8);+    oh |= tl >> (64 - shift * 8);+    out->u[1] = (uint32_t)(ol >> 32);+    out->u[0] = (uint32_t)ol;+    out->u[3] = (uint32_t)(oh >> 32);+    out->u[2] = (uint32_t)oh;+}+#endif++/**+ * This function represents the recursion formula.+ * @param r output+ * @param a a 128-bit part of the internal state array+ * @param b a 128-bit part of the internal state array+ * @param c a 128-bit part of the internal state array+ * @param d a 128-bit part of the internal state array+ */+#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))+#ifdef ONLY64+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,+				w128_t *d) {+    w128_t x;+    w128_t y;++    lshift128(&x, a, SL2);+    rshift128(&y, c, SR2);+    r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK2) ^ y.u[0] +	^ (d->u[0] << SL1);+    r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK1) ^ y.u[1] +	^ (d->u[1] << SL1);+    r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK4) ^ y.u[2] +	^ (d->u[2] << SL1);+    r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK3) ^ y.u[3] +	^ (d->u[3] << SL1);+}+#else+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,+				w128_t *d) {+    w128_t x;+    w128_t y;++    lshift128(&x, a, SL2);+    rshift128(&y, c, SR2);+    r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0] +	^ (d->u[0] << SL1);+    r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1] +	^ (d->u[1] << SL1);+    r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2] +	^ (d->u[2] << SL1);+    r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3] +	^ (d->u[3] << SL1);+}+#endif+#endif++#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))+/**+ * This function fills the internal state array with pseudorandom+ * integers.+ */+inline static void gen_rand_all(void) {+    int i;+    w128_t *r1, *r2;++    r1 = &sfmt[N - 2];+    r2 = &sfmt[N - 1];+    for (i = 0; i < N - POS1; i++) {+	do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1], r1, r2);+	r1 = r2;+	r2 = &sfmt[i];+    }+    for (; i < N; i++) {+	do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1 - N], r1, r2);+	r1 = r2;+	r2 = &sfmt[i];+    }+}++/**+ * This function fills the user-specified array with pseudorandom+ * integers.+ *+ * @param array an 128-bit array to be filled by pseudorandom numbers.  + * @param size number of 128-bit pseudorandom numbers to be generated.+ */+inline static void gen_rand_array(w128_t *array, int size) {+    int i, j;+    w128_t *r1, *r2;++    r1 = &sfmt[N - 2];+    r2 = &sfmt[N - 1];+    for (i = 0; i < N - POS1; i++) {+	do_recursion(&array[i], &sfmt[i], &sfmt[i + POS1], r1, r2);+	r1 = r2;+	r2 = &array[i];+    }+    for (; i < N; i++) {+	do_recursion(&array[i], &sfmt[i], &array[i + POS1 - N], r1, r2);+	r1 = r2;+	r2 = &array[i];+    }+    for (; i < size - N; i++) {+	do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2);+	r1 = r2;+	r2 = &array[i];+    }+    for (j = 0; j < 2 * N - size; j++) {+	sfmt[j] = array[j + size - N];+    }+    for (; i < size; i++, j++) {+	do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2);+	r1 = r2;+	r2 = &array[i];+	sfmt[j] = array[i];+    }+}+#endif++#if defined(BIG_ENDIAN64) && !defined(ONLY64) && !defined(HAVE_ALTIVEC)+inline static void swap(w128_t *array, int size) {+    int i;+    uint32_t x, y;++    for (i = 0; i < size; i++) {+	x = array[i].u[0];+	y = array[i].u[2];+	array[i].u[0] = array[i].u[1];+	array[i].u[2] = array[i].u[3];+	array[i].u[1] = x;+	array[i].u[3] = y;+    }+}+#endif+/**+ * This function represents a function used in the initialization+ * by init_by_array+ * @param x 32-bit integer+ * @return 32-bit integer+ */+static uint32_t func1(uint32_t x) {+    return (x ^ (x >> 27)) * (uint32_t)1664525UL;+}++/**+ * This function represents a function used in the initialization+ * by init_by_array+ * @param x 32-bit integer+ * @return 32-bit integer+ */+static uint32_t func2(uint32_t x) {+    return (x ^ (x >> 27)) * (uint32_t)1566083941UL;+}++/**+ * This function certificate the period of 2^{MEXP}+ */+static void period_certification(void) {+    int inner = 0;+    int i, j;+    uint32_t work;++    for (i = 0; i < 4; i++)+	inner ^= psfmt32[idxof(i)] & parity[i];+    for (i = 16; i > 0; i >>= 1)+	inner ^= inner >> i;+    inner &= 1;+    /* check OK */+    if (inner == 1) {+	return;+    }+    /* check NG, and modification */+    for (i = 0; i < 4; i++) {+	work = 1;+	for (j = 0; j < 32; j++) {+	    if ((work & parity[i]) != 0) {+		psfmt32[idxof(i)] ^= work;+		return;+	    }+	    work = work << 1;+	}+    }+}++/*----------------+  PUBLIC FUNCTIONS+  ----------------*/+/**+ * This function returns the identification string.+ * The string shows the word size, the Mersenne exponent,+ * and all parameters of this generator.+ */+const char *get_idstring(void) {+    return IDSTR;+}++int get_initialized(void) { return initialized; }++/**+ * This function returns the minimum size of array used for \b+ * fill_array32() function.+ * @return minimum size of array used for fill_array32() function.+ */+int get_min_array_size32(void) {+    return N32;+}++/**+ * This function returns the minimum size of array used for \b+ * fill_array64() function.+ * @return minimum size of array used for fill_array64() function.+ */+int get_min_array_size64(void) {+    return N64;+}++#ifndef ONLY64+/**+ * This function generates and returns 32-bit pseudorandom number.+ * init_gen_rand or init_by_array must be called before this function.+ * @return 32-bit pseudorandom number+ */+uint32_t gen_rand32(void) {+    uint32_t r;++    assert(initialized);+    if (idx >= N32) {+	gen_rand_all();+	idx = 0;+    }+    r = psfmt32[idx++];+    return r;+}+#endif+/**+ * 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. + * @return 64-bit pseudorandom number+ */+uint64_t gen_rand64(void) {+#if defined(BIG_ENDIAN64) && !defined(ONLY64)+    uint32_t r1, r2;+#else+    uint64_t r;+#endif++    assert(initialized);+    assert(idx % 2 == 0);++    if (idx >= N32) {+	gen_rand_all();+	idx = 0;+    }+#if defined(BIG_ENDIAN64) && !defined(ONLY64)+    r1 = psfmt32[idx];+    r2 = psfmt32[idx + 1];+    idx += 2;+    return ((uint64_t)r2 << 32) | r1;+#else+    r = psfmt64[idx / 2];+    idx += 2;+    return r;+#endif+}++#ifndef ONLY64+/**+ * This function generates pseudorandom 32-bit integers in the+ * specified array[] by one call. The number of pseudorandom integers+ * is specified by the argument size, which must be at least 624 and a+ * multiple of four.  The generation by this function is much faster+ * than the following gen_rand function.+ *+ * For initialization, init_gen_rand or init_by_array must be called+ * before the first call of this function. This function can not be+ * used after calling gen_rand function, without initialization.+ *+ * @param array an array where pseudorandom 32-bit integers are filled+ * by this function.  The pointer to the array must be \b "aligned"+ * (namely, must be a multiple of 16) in the SIMD version, since it+ * refers to the address of a 128-bit integer.  In the standard C+ * version, the pointer is arbitrary.+ *+ * @param size the number of 32-bit pseudorandom integers to be+ * generated.  size must be a multiple of 4, and greater than or equal+ * to (MEXP / 128 + 1) * 4.+ *+ * @note \b memalign or \b posix_memalign is available to get aligned+ * memory. Mac OSX doesn't have these functions, but \b malloc of OSX+ * returns the pointer to the aligned memory block.+ */+void fill_array32(uint32_t *array, int size) {+    assert(initialized);+    assert(idx == N32);+    assert(size % 4 == 0);+    assert(size >= N32);++    gen_rand_array((w128_t *)array, size / 4);+    idx = N32;+}+#endif++/**+ * This function generates pseudorandom 64-bit integers in the+ * specified array[] by one call. The number of pseudorandom integers+ * is specified by the argument size, which must be at least 312 and a+ * multiple of two.  The generation by this function is much faster+ * than the following gen_rand function.+ *+ * For initialization, init_gen_rand or init_by_array must be called+ * before the first call of this function. This function can not be+ * used after calling gen_rand function, without initialization.+ *+ * @param array an array where pseudorandom 64-bit integers are filled+ * by this function.  The pointer to the array must be "aligned"+ * (namely, must be a multiple of 16) in the SIMD version, since it+ * refers to the address of a 128-bit integer.  In the standard C+ * version, the pointer is arbitrary.+ *+ * @param size the number of 64-bit pseudorandom integers to be+ * generated.  size must be a multiple of 2, and greater than or equal+ * to (MEXP / 128 + 1) * 2+ *+ * @note \b memalign or \b posix_memalign is available to get aligned+ * memory. Mac OSX doesn't have these functions, but \b malloc of OSX+ * returns the pointer to the aligned memory block.+ */+void fill_array64(uint64_t *array, int size) {+    assert(initialized);+    assert(idx == N32);+    assert(size % 2 == 0);+    assert(size >= N64);++    gen_rand_array((w128_t *)array, size / 2);+    idx = N32;++#if defined(BIG_ENDIAN64) && !defined(ONLY64)+    swap((w128_t *)array, size /2);+#endif+}++/**+ * This function initializes the internal state array with a 32-bit+ * integer seed.+ *+ * @param seed a 32-bit integer used as the seed.+ */+void init_gen_rand(uint32_t seed) {+    int i;++    psfmt32[idxof(0)] = seed;+    for (i = 1; i < N32; i++) {+	psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)] +					    ^ (psfmt32[idxof(i - 1)] >> 30))+	    + i;+    }+    idx = N32;+    period_certification();+    initialized = 1;+}++/**+ * This function initializes the internal state array,+ * with an array of 32-bit integers used as the seeds+ * @param init_key the array of 32-bit integers, used as a seed.+ * @param key_length the length of init_key.+ */+void init_by_array(uint32_t *init_key, int key_length) {+    int i, j, count;+    uint32_t r;+    int lag;+    int mid;+    int size = N * 4;++    if (size >= 623) {+	lag = 11;+    } else if (size >= 68) {+	lag = 7;+    } else if (size >= 39) {+	lag = 5;+    } else {+	lag = 3;+    }+    mid = (size - lag) / 2;++    memset(sfmt, 0x8b, sizeof(sfmt));+    if (key_length + 1 > N32) {+	count = key_length + 1;+    } else {+	count = N32;+    }+    r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)] +	      ^ psfmt32[idxof(N32 - 1)]);+    psfmt32[idxof(mid)] += r;+    r += key_length;+    psfmt32[idxof(mid + lag)] += r;+    psfmt32[idxof(0)] = r;++    count--;+    for (i = 1, j = 0; (j < count) && (j < key_length); j++) {+	r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] +		  ^ psfmt32[idxof((i + N32 - 1) % N32)]);+	psfmt32[idxof((i + mid) % N32)] += r;+	r += init_key[j] + i;+	psfmt32[idxof((i + mid + lag) % N32)] += r;+	psfmt32[idxof(i)] = r;+	i = (i + 1) % N32;+    }+    for (; j < count; j++) {+	r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] +		  ^ psfmt32[idxof((i + N32 - 1) % N32)]);+	psfmt32[idxof((i + mid) % N32)] += r;+	r += i;+	psfmt32[idxof((i + mid + lag) % N32)] += r;+	psfmt32[idxof(i)] = r;+	i = (i + 1) % N32;+    }+    for (j = 0; j < N32; j++) {+	r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % N32)] +		  + psfmt32[idxof((i + N32 - 1) % N32)]);+	psfmt32[idxof((i + mid) % N32)] ^= r;+	r -= i;+	psfmt32[idxof((i + mid + lag) % N32)] ^= r;+	psfmt32[idxof(i)] = r;+	i = (i + 1) % N32;+    }++    idx = N32;+    period_certification();+    initialized = 1;+}

+ include/SFMT-alti.h view

@@ -0,0 +1,156 @@+/** + * @file SFMT-alti.h + *+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)+ * pseudorandom number generator+ *+ * @author Mutsuo Saito (Hiroshima University)+ * @author Makoto Matsumoto (Hiroshima University)+ *+ * Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima+ * University. All rights reserved.+ *+ * The new BSD License is applied to this software.+ * see LICENSE.txt+ */++#ifndef SFMT_ALTI_H+#define SFMT_ALTI_H++inline static vector unsigned int vec_recursion(vector unsigned int a,+						vector unsigned int b,+						vector unsigned int c,+						vector unsigned int d)+    ALWAYSINLINE;++/**+ * This function represents the recursion formula in AltiVec and BIG ENDIAN.+ * @param a a 128-bit part of the interal state array+ * @param b a 128-bit part of the interal state array+ * @param c a 128-bit part of the interal state array+ * @param d a 128-bit part of the interal state array+ * @return output+ */+inline static vector unsigned int vec_recursion(vector unsigned int a,+						vector unsigned int b,+						vector unsigned int c,+						vector unsigned int d) {++    const vector unsigned int sl1 = ALTI_SL1;+    const vector unsigned int sr1 = ALTI_SR1;+#ifdef ONLY64+    const vector unsigned int mask = ALTI_MSK64;+    const vector unsigned char perm_sl = ALTI_SL2_PERM64;+    const vector unsigned char perm_sr = ALTI_SR2_PERM64;+#else+    const vector unsigned int mask = ALTI_MSK;+    const vector unsigned char perm_sl = ALTI_SL2_PERM;+    const vector unsigned char perm_sr = ALTI_SR2_PERM;+#endif+    vector unsigned int v, w, x, y, z;+    x = vec_perm(a, (vector unsigned int)perm_sl, perm_sl);+    v = a;+    y = vec_sr(b, sr1);+    z = vec_perm(c, (vector unsigned int)perm_sr, perm_sr);+    w = vec_sl(d, sl1);+    z = vec_xor(z, w);+    y = vec_and(y, mask);+    v = vec_xor(v, x);+    z = vec_xor(z, y);+    z = vec_xor(z, v);+    return z;+}++/**+ * This function fills the internal state array with pseudorandom+ * integers.+ */+inline static void gen_rand_all(void) {+    int i;+    vector unsigned int r, r1, r2;++    r1 = sfmt[N - 2].s;+    r2 = sfmt[N - 1].s;+    for (i = 0; i < N - POS1; i++) {+	r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);+	sfmt[i].s = r;+	r1 = r2;+	r2 = r;+    }+    for (; i < N; i++) {+	r = vec_recursion(sfmt[i].s, sfmt[i + POS1 - N].s, r1, r2);+	sfmt[i].s = r;+	r1 = r2;+	r2 = r;+    }+}++/**+ * This function fills the user-specified array with pseudorandom+ * integers.+ *+ * @param array an 128-bit array to be filled by pseudorandom numbers.  + * @param size number of 128-bit pesudorandom numbers to be generated.+ */+inline static void gen_rand_array(w128_t *array, int size) {+    int i, j;+    vector unsigned int r, r1, r2;++    r1 = sfmt[N - 2].s;+    r2 = sfmt[N - 1].s;+    for (i = 0; i < N - POS1; i++) {+	r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);+	array[i].s = r;+	r1 = r2;+	r2 = r;+    }+    for (; i < N; i++) {+	r = vec_recursion(sfmt[i].s, array[i + POS1 - N].s, r1, r2);+	array[i].s = r;+	r1 = r2;+	r2 = r;+    }+    /* main loop */+    for (; i < size - N; i++) {+	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);+	array[i].s = r;+	r1 = r2;+	r2 = r;+    }+    for (j = 0; j < 2 * N - size; j++) {+	sfmt[j].s = array[j + size - N].s;+    }+    for (; i < size; i++) {+	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);+	array[i].s = r;+	sfmt[j++].s = r;+	r1 = r2;+	r2 = r;+    }+}++#ifndef ONLY64+#if defined(__APPLE__)+#define ALTI_SWAP (vector unsigned char) \+	(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11)+#else+#define ALTI_SWAP {4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11}+#endif+/**+ * This function swaps high and low 32-bit of 64-bit integers in user+ * specified array.+ *+ * @param array an 128-bit array to be swaped.+ * @param size size of 128-bit array.+ */+inline static void swap(w128_t *array, int size) {+    int i;+    const vector unsigned char perm = ALTI_SWAP;++    for (i = 0; i < size; i++) {+	array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);+    }+}+#endif++#endif

+ include/SFMT-params.h view

@@ -0,0 +1,97 @@+#ifndef SFMT_PARAMS_H+#define SFMT_PARAMS_H++#if !defined(MEXP)+#ifdef __GNUC__+  #warning "MEXP is not defined. I assume MEXP is 19937."+#endif+  #define MEXP 19937+#endif+/*-----------------+  BASIC DEFINITIONS+  -----------------*/+/** Mersenne Exponent. The period of the sequence + *  is a multiple of 2^MEXP-1.+ * #define MEXP 19937 */+/** SFMT generator has an internal state array of 128-bit integers,+ * and N is its size. */+#define N (MEXP / 128 + 1)+/** N32 is the size of internal state array when regarded as an array+ * of 32-bit integers.*/+#define N32 (N * 4)+/** N64 is the size of internal state array when regarded as an array+ * of 64-bit integers.*/+#define N64 (N * 2)++/*----------------------+  the parameters of SFMT+  following definitions are in paramsXXXX.h file.+  ----------------------*/+/** the pick up position of the array.+#define POS1 122 +*/++/** the parameter of shift left as four 32-bit registers.+#define SL1 18+ */++/** the parameter of shift left as one 128-bit register. + * The 128-bit integer is shifted by (SL2 * 8) bits. +#define SL2 1 +*/++/** the parameter of shift right as four 32-bit registers.+#define SR1 11+*/++/** the parameter of shift right as one 128-bit register. + * The 128-bit integer is shifted by (SL2 * 8) bits. +#define SR2 1 +*/++/** A bitmask, used in the recursion.  These parameters are introduced+ * to break symmetry of SIMD.+#define MSK1 0xdfffffefU+#define MSK2 0xddfecb7fU+#define MSK3 0xbffaffffU+#define MSK4 0xbffffff6U +*/++/** These definitions are part of a 128-bit period certification vector.+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0x00000000U+#define PARITY4	0xc98e126aU+*/++#if MEXP == 607+  #include "SFMT-params607.h"+#elif MEXP == 1279+  #include "SFMT-params1279.h"+#elif MEXP == 2281+  #include "SFMT-params2281.h"+#elif MEXP == 4253+  #include "SFMT-params4253.h"+#elif MEXP == 11213+  #include "SFMT-params11213.h"+#elif MEXP == 19937+  #include "SFMT-params19937.h"+#elif MEXP == 44497+  #include "SFMT-params44497.h"+#elif MEXP == 86243+  #include "SFMT-params86243.h"+#elif MEXP == 132049+  #include "SFMT-params132049.h"+#elif MEXP == 216091+  #include "SFMT-params216091.h"+#else+#ifdef __GNUC__+  #error "MEXP is not valid."+  #undef MEXP+#else+  #undef MEXP+#endif++#endif++#endif /* SFMT_PARAMS_H */

+ include/SFMT-params11213.h view

@@ -0,0 +1,48 @@+#ifndef SFMT_PARAMS11213_H+#define SFMT_PARAMS11213_H++#include "SFMT.h"++#define POS1	68+#define SL1	14+#define SL2	3+#define SR1	7+#define SR2	3+#define MSK1	0xeffff7fbU+#define MSK2	0xffffffefU+#define MSK3	0xdfdfbfffU+#define MSK4	0x7fffdbfdU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0xe8148000U+#define PARITY4	0xd0c7afa3U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}+    #define ALTI_SL2_PERM64	{3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}+    #define ALTI_SR2_PERM	{5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}+    #define ALTI_SR2_PERM64	{13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}+#endif	/* For OSX */+#define IDSTR	"SFMT-11213:68-14-3-7-3:effff7fb-ffffffef-dfdfbfff-7fffdbfd"++#endif /* SFMT_PARAMS11213_H */

+ include/SFMT-params1279.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS1279_H+#define SFMT_PARAMS1279_H++#define POS1	7+#define SL1	14+#define SL2	3+#define SR1	5+#define SR2	1+#define MSK1	0xf7fefffdU+#define MSK2	0x7fefcfffU+#define MSK3	0xaff3ef3fU+#define MSK4	0xb5ffff7fU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0x00000000U+#define PARITY4	0x20000000U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}+    #define ALTI_SL2_PERM64	{3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-1279:7-14-3-5-1:f7fefffd-7fefcfff-aff3ef3f-b5ffff7f"++#endif /* SFMT_PARAMS1279_H */

+ include/SFMT-params132049.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS132049_H+#define SFMT_PARAMS132049_H++#define POS1	110+#define SL1	19+#define SL2	1+#define SR1	21+#define SR2	1+#define MSK1	0xffffbb5fU+#define MSK2	0xfb6ebf95U+#define MSK3	0xfffefffaU+#define MSK4	0xcff77fffU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0xcb520000U+#define PARITY4	0xc7e91c7dU+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}+    #define ALTI_SL2_PERM64	{1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-132049:110-19-1-21-1:ffffbb5f-fb6ebf95-fffefffa-cff77fff"++#endif /* SFMT_PARAMS132049_H */

+ include/SFMT-params19937.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS19937_H+#define SFMT_PARAMS19937_H++#define POS1	122+#define SL1	18+#define SL2	1+#define SR1	11+#define SR2	1+#define MSK1	0xdfffffefU+#define MSK2	0xddfecb7fU+#define MSK3	0xbffaffffU+#define MSK4	0xbffffff6U+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0x00000000U+#define PARITY4	0x13c9e684U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}+    #define ALTI_SL2_PERM64	{1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6"++#endif /* SFMT_PARAMS19937_H */

+ include/SFMT-params216091.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS216091_H+#define SFMT_PARAMS216091_H++#define POS1	627+#define SL1	11+#define SL2	3+#define SR1	10+#define SR2	1+#define MSK1	0xbff7bff7U+#define MSK2	0xbfffffffU+#define MSK3	0xbffffa7fU+#define MSK4	0xffddfbfbU+#define PARITY1	0xf8000001U+#define PARITY2	0x89e80709U+#define PARITY3	0x3bd2b64bU+#define PARITY4	0x0c64b1e4U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}+    #define ALTI_SL2_PERM64	{3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-216091:627-11-3-10-1:bff7bff7-bfffffff-bffffa7f-ffddfbfb"++#endif /* SFMT_PARAMS216091_H */

+ include/SFMT-params2281.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS2281_H+#define SFMT_PARAMS2281_H++#define POS1	12+#define SL1	19+#define SL2	1+#define SR1	5+#define SR2	1+#define MSK1	0xbff7ffbfU+#define MSK2	0xfdfffffeU+#define MSK3	0xf7ffef7fU+#define MSK4	0xf2f7cbbfU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0x00000000U+#define PARITY4	0x41dfa600U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}+    #define ALTI_SL2_PERM64	{1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-2281:12-19-1-5-1:bff7ffbf-fdfffffe-f7ffef7f-f2f7cbbf"++#endif /* SFMT_PARAMS2281_H */

+ include/SFMT-params4253.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS4253_H+#define SFMT_PARAMS4253_H++#define POS1	17+#define SL1	20+#define SL2	1+#define SR1	7+#define SR2	1+#define MSK1	0x9f7bffffU+#define MSK2	0x9fffff5fU+#define MSK3	0x3efffffbU+#define MSK4	0xfffff7bbU+#define PARITY1	0xa8000001U+#define PARITY2	0xaf5390a3U+#define PARITY3	0xb740b3f8U+#define PARITY4	0x6c11486dU+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}+    #define ALTI_SL2_PERM64	{1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-4253:17-20-1-7-1:9f7bffff-9fffff5f-3efffffb-fffff7bb"++#endif /* SFMT_PARAMS4253_H */

+ include/SFMT-params44497.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS44497_H+#define SFMT_PARAMS44497_H++#define POS1	330+#define SL1	5+#define SL2	3+#define SR1	9+#define SR2	3+#define MSK1	0xeffffffbU+#define MSK2	0xdfbebfffU+#define MSK3	0xbfbf7befU+#define MSK4	0x9ffd7bffU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0xa3ac4000U+#define PARITY4	0xecc1327aU+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}+    #define ALTI_SL2_PERM64	{3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}+    #define ALTI_SR2_PERM	{5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}+    #define ALTI_SR2_PERM64	{13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}+#endif	/* For OSX */+#define IDSTR	"SFMT-44497:330-5-3-9-3:effffffb-dfbebfff-bfbf7bef-9ffd7bff"++#endif /* SFMT_PARAMS44497_H */

+ include/SFMT-params607.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS607_H+#define SFMT_PARAMS607_H++#define POS1	2+#define SL1	15+#define SL2	3+#define SR1	13+#define SR2	3+#define MSK1	0xfdff37ffU+#define MSK2	0xef7f3f7dU+#define MSK3	0xff777b7dU+#define MSK4	0x7ff7fb2fU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0x00000000U+#define PARITY4	0x5986f054U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}+    #define ALTI_SL2_PERM64	{3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}+    #define ALTI_SR2_PERM	{5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}+    #define ALTI_SR2_PERM64	{13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}+#endif	/* For OSX */+#define IDSTR	"SFMT-607:2-15-3-13-3:fdff37ff-ef7f3f7d-ff777b7d-7ff7fb2f"++#endif /* SFMT_PARAMS607_H */

+ include/SFMT-params86243.h view

@@ -0,0 +1,46 @@+#ifndef SFMT_PARAMS86243_H+#define SFMT_PARAMS86243_H++#define POS1	366+#define SL1	6+#define SL2	7+#define SR1	19+#define SR2	1+#define MSK1	0xfdbffbffU+#define MSK2	0xbff7ff3fU+#define MSK3	0xfd77efffU+#define MSK4	0xbf9ff3ffU+#define PARITY1	0x00000001U+#define PARITY2	0x00000000U+#define PARITY3	0x00000000U+#define PARITY4	0xe9528d85U+++/* PARAMETERS FOR ALTIVEC */+#if defined(__APPLE__)	/* For OSX */+    #define ALTI_SL1	(vector unsigned int)(SL1, SL1, SL1, SL1)+    #define ALTI_SR1	(vector unsigned int)(SR1, SR1, SR1, SR1)+    #define ALTI_MSK	(vector unsigned int)(MSK1, MSK2, MSK3, MSK4)+    #define ALTI_MSK64 \+	(vector unsigned int)(MSK2, MSK1, MSK4, MSK3)+    #define ALTI_SL2_PERM \+	(vector unsigned char)(25,25,25,25,3,25,25,25,7,0,1,2,11,4,5,6)+    #define ALTI_SL2_PERM64 \+	(vector unsigned char)(7,25,25,25,25,25,25,25,15,0,1,2,3,4,5,6)+    #define ALTI_SR2_PERM \+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)+    #define ALTI_SR2_PERM64 \+	(vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)+#else	/* For OTHER OSs(Linux?) */+    #define ALTI_SL1	{SL1, SL1, SL1, SL1}+    #define ALTI_SR1	{SR1, SR1, SR1, SR1}+    #define ALTI_MSK	{MSK1, MSK2, MSK3, MSK4}+    #define ALTI_MSK64	{MSK2, MSK1, MSK4, MSK3}+    #define ALTI_SL2_PERM	{25,25,25,25,3,25,25,25,7,0,1,2,11,4,5,6}+    #define ALTI_SL2_PERM64	{7,25,25,25,25,25,25,25,15,0,1,2,3,4,5,6}+    #define ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}+    #define ALTI_SR2_PERM64	{15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}+#endif	/* For OSX */+#define IDSTR	"SFMT-86243:366-6-7-19-1:fdbffbff-bff7ff3f-fd77efff-bf9ff3ff"++#endif /* SFMT_PARAMS86243_H */

+ include/SFMT-sse2.h view

@@ -0,0 +1,121 @@+/** + * @file  SFMT-sse2.h+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) for Intel SSE2+ *+ * @author Mutsuo Saito (Hiroshima University)+ * @author Makoto Matsumoto (Hiroshima University)+ *+ * @note We assume LITTLE ENDIAN in this file+ *+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima+ * University. All rights reserved.+ *+ * The new BSD License is applied to this software, see LICENSE.txt+ */++#ifndef SFMT_SSE2_H+#define SFMT_SSE2_H++PRE_ALWAYS static __m128i mm_recursion(__m128i *a, __m128i *b, __m128i c,+				   __m128i d, __m128i mask) ALWAYSINLINE;++/**+ * This function represents the recursion formula.+ * @param a a 128-bit part of the interal state array+ * @param b a 128-bit part of the interal state array+ * @param c a 128-bit part of the interal state array+ * @param d a 128-bit part of the interal state array+ * @param mask 128-bit mask+ * @return output+ */+PRE_ALWAYS static __m128i mm_recursion(__m128i *a, __m128i *b, +				   __m128i c, __m128i d, __m128i mask) {+    __m128i v, x, y, z;+    +    x = _mm_load_si128(a);+    y = _mm_srli_epi32(*b, SR1);+    z = _mm_srli_si128(c, SR2);+    v = _mm_slli_epi32(d, SL1);+    z = _mm_xor_si128(z, x);+    z = _mm_xor_si128(z, v);+    x = _mm_slli_si128(x, SL2);+    y = _mm_and_si128(y, mask);+    z = _mm_xor_si128(z, x);+    z = _mm_xor_si128(z, y);+    return z;+}++/**+ * This function fills the internal state array with pseudorandom+ * integers.+ */+inline static void gen_rand_all(void) {+    int i;+    __m128i r, r1, r2, mask;+    mask = _mm_set_epi32(MSK4, MSK3, MSK2, MSK1);++    r1 = _mm_load_si128(&sfmt[N - 2].si);+    r2 = _mm_load_si128(&sfmt[N - 1].si);+    for (i = 0; i < N - POS1; i++) {+	r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1].si, r1, r2, mask);+	_mm_store_si128(&sfmt[i].si, r);+	r1 = r2;+	r2 = r;+    }+    for (; i < N; i++) {+	r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1 - N].si, r1, r2, mask);+	_mm_store_si128(&sfmt[i].si, r);+	r1 = r2;+	r2 = r;+    }+}++/**+ * This function fills the user-specified array with pseudorandom+ * integers.+ *+ * @param array an 128-bit array to be filled by pseudorandom numbers.  + * @param size number of 128-bit pesudorandom numbers to be generated.+ */+inline static void gen_rand_array(w128_t *array, int size) {+    int i, j;+    __m128i r, r1, r2, mask;+    mask = _mm_set_epi32(MSK4, MSK3, MSK2, MSK1);++    r1 = _mm_load_si128(&sfmt[N - 2].si);+    r2 = _mm_load_si128(&sfmt[N - 1].si);+    for (i = 0; i < N - POS1; i++) {+	r = mm_recursion(&sfmt[i].si, &sfmt[i + POS1].si, r1, r2, mask);+	_mm_store_si128(&array[i].si, r);+	r1 = r2;+	r2 = r;+    }+    for (; i < N; i++) {+	r = mm_recursion(&sfmt[i].si, &array[i + POS1 - N].si, r1, r2, mask);+	_mm_store_si128(&array[i].si, r);+	r1 = r2;+	r2 = r;+    }+    /* main loop */+    for (; i < size - N; i++) {+	r = mm_recursion(&array[i - N].si, &array[i + POS1 - N].si, r1, r2,+			 mask);+	_mm_store_si128(&array[i].si, r);+	r1 = r2;+	r2 = r;+    }+    for (j = 0; j < 2 * N - size; j++) {+	r = _mm_load_si128(&array[j + size - N].si);+	_mm_store_si128(&sfmt[j].si, r);+    }+    for (; i < size; i++) {+	r = mm_recursion(&array[i - N].si, &array[i + POS1 - N].si, r1, r2,+			 mask);+	_mm_store_si128(&array[i].si, r);+	_mm_store_si128(&sfmt[j++].si, r);+	r1 = r2;+	r2 = r;+    }+}++#endif

+ include/SFMT.h view

@@ -0,0 +1,168 @@+/** + * @file SFMT.h + *+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom+ * number generator+ *+ * @author Mutsuo Saito (Hiroshima University)+ * @author Makoto Matsumoto (Hiroshima University)+ *+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima+ * University. All rights reserved.+ *+ * The new BSD License is applied to this software.+ * see LICENSE.txt+ *+ * @note We assume that your system has inttypes.h.  If your system+ * doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,+ * and you have to define PRIu64 and PRIx64 in this file as follows:+ * @verbatim+ typedef unsigned int uint32_t+ typedef unsigned long long uint64_t  + #define PRIu64 "llu"+ #define PRIx64 "llx"+@endverbatim+ * uint32_t must be exactly 32-bit unsigned integer type (no more, no+ * less), and uint64_t must be exactly 64-bit unsigned integer type.+ * PRIu64 and PRIx64 are used for printf function to print 64-bit+ * unsigned int and 64-bit unsigned int in hexadecimal format.+ */++#ifndef SFMT_H+#define SFMT_H++#include <stdio.h>++#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)+  #include <inttypes.h>+#elif defined(_MSC_VER) || defined(__BORLANDC__)+  typedef unsigned int uint32_t;+  typedef unsigned __int64 uint64_t;+  #define inline __inline+#else+  #include <inttypes.h>+  #if defined(__GNUC__)+    #define inline __inline__+  #endif+#endif++#ifndef PRIu64+  #if defined(_MSC_VER) || defined(__BORLANDC__)+    #define PRIu64 "I64u"+    #define PRIx64 "I64x"+  #else+    #define PRIu64 "llu"+    #define PRIx64 "llx"+  #endif+#endif++#if defined(__GNUC__)+#define ALWAYSINLINE __attribute__((always_inline))+#else+#define ALWAYSINLINE+#endif++#if defined(_MSC_VER)+  #if _MSC_VER >= 1200+    #define PRE_ALWAYS __forceinline+  #else+    #define PRE_ALWAYS inline+  #endif+#else+  #define PRE_ALWAYS inline+#endif++uint32_t gen_rand32(void);+uint64_t gen_rand64(void);+void fill_array32(uint32_t *array, int size);+void fill_array64(uint64_t *array, int size);+void init_gen_rand(uint32_t seed);+void init_by_array(uint32_t *init_key, int key_length);+const char *get_idstring(void);+int get_min_array_size32(void);+int get_min_array_size64(void);+int get_initialized(void);++/* These real versions are due to Isaku Wada */+/** generates a random number on [0,1]-real-interval */+inline static double to_real1(uint32_t v)+{+    return v * (1.0/4294967295.0); +    /* divided by 2^32-1 */ +}++/** generates a random number on [0,1]-real-interval */+inline static double genrand_real1(void)+{+    return to_real1(gen_rand32());+}++/** generates a random number on [0,1)-real-interval */+inline static double to_real2(uint32_t v)+{+    return v * (1.0/4294967296.0); +    /* divided by 2^32 */+}++/** generates a random number on [0,1)-real-interval */+inline static double genrand_real2(void)+{+    return to_real2(gen_rand32());+}++/** generates a random number on (0,1)-real-interval */+inline static double to_real3(uint32_t v)+{+    return (((double)v) + 0.5)*(1.0/4294967296.0); +    /* divided by 2^32 */+}++/** generates a random number on (0,1)-real-interval */+inline static double genrand_real3(void)+{+    return to_real3(gen_rand32());+}+/** These real versions are due to Isaku Wada */++/** generates a random number on [0,1) with 53-bit resolution*/+inline static double to_res53(uint64_t v) +{ +    return v * (1.0/18446744073709551616.0L);+}++/** generates a random number on [0,1) with 53-bit resolution from two+ * 32 bit integers */+inline static double to_res53_mix(uint32_t x, uint32_t y) +{ +    return to_res53(x | ((uint64_t)y << 32));+}++/** generates a random number on [0,1) with 53-bit resolution+ */+inline static double genrand_res53(void) +{ +    return to_res53(gen_rand64());+} ++/** generates a random number on [0,1) with 53-bit resolution+    using 32bit integer.+ */+inline static double genrand_res53_mix(void) +{ +    uint32_t x, y;++    x = gen_rand32();+    y = gen_rand32();+    return to_res53_mix(x, y);+} ++/* build a 64 bit integer from two 32 bit ones */+inline static uint64_t gen_rand64_mix(void)+{ +    uint32_t x, y;++    x = gen_rand32();+    y = gen_rand32();+    return (x | ((uint64_t)y << 32));+}+#endif

+ mersenne-random.cabal view

@@ -0,0 +1,109 @@+name:            mersenne-random+version:         0.1+homepage:        http://code.haskell.org/~dons/code/mersenne-random+synopsis:        Generate high quality pseudorandom numbers using a SIMD Fast 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+    .+    This library uses SFMT, the SIMD-oriented Fast Mersenne Twister, a+    variant of Mersenne Twister that is much faster than the original. +    It is designed to be fast when it runs on 128-bit SIMD. It can be+    compiled with either SSE2 and PowerPC AltiVec support, to take+    advantage of these instructions.+    .+    By default the period of the function is 2^19937-1, however, you can+    compile in other defaults. Note that this algorithm on its own+    is not cryptographically secure.+    .+    For more information about the algorithm and implementation, see+    the SFMT homepage,+    .+    <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>+    .+    and, Mutsuo Saito and Makoto Matsumoto,+    /SIMD-oriented Fast Mersenne Twister: a 128-bit Pseudorandom Number+    Generator/, in the Proceedings of MCQMC2006, here:+    .+    <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/sfmt.pdf>+    .+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, Hugs ==2005+build-type:      Configure++flag small_base+  description: Build with new smaller base library+  default:     False++flag use_sse2+  description: Build with SSE2 support.+  default:     False++flag use_altivec+  description: Build with Altivec support.+  default:     False++flag big_endian64+  description: Build for a big endian 64 bit machine.+  default:     False++library+    exposed-modules: System.Random.Mersenne+    extensions:      CPP, ForeignFunctionInterface, BangPatterns++    if flag(small_base)+        build-depends: base  < 3+    else+        build-depends: base >= 3, old-time++    -- For information on how to set different periods, or tune +    -- for your arch, see,+    --+    -- <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/howto-compile.html>+    --+    -- SSE2 supported on: Pentium M, Pentium 4, Core, Core 2 etc.+    -- See: http://en.wikipedia.org/wiki/SSE2#CPUs_supporting_SSE2+    --+    -- Enable use_sse2 flag if you have one of those archs.+    --+    -- Works well on core 2 duo.+    --+    -- Enable use_altivec flag to use smid on powerpc.+    -- Enable big_endian64 flag on a big endian machine 64 bit machine+    --  (e.g. UltraSparc)+    --+    cc-options:+        -DMEXP=19937+        -DNDEBUG +        -O3 -finline-functions -fomit-frame-pointer+        -fno-strict-aliasing --param max-inline-insns-single=1800++    if flag(use_sse2)+        cc-options:+            -msse2 +            -DHAVE_SSE2++    if flag(big_endian64)+        cc-options:+            -DBIG_ENDIAN64++    if flag(use_altivec)+        cc-options:+            -DHAVE_ALTIVEC++    ghc-options:     -Wall -O2 -fvia-C -fexcess-precision++    c-sources:        cbits/SFMT.c+    include-dirs:     include+    includes:         SFMT.h+    install-includes: SFMT.h+    extra-libraries:  m

+ tests/Makefile view

@@ -0,0 +1,11 @@+all:+	ghc -O2 -ddump-simpl-stats -no-recomp Unit.hs --make+	./Unit+	ghc -O2 pi.hs --make+	./pi 100000000++hpc:+	ghc --make -i. -i../ -optc-DMEXP=19937 -optc-DNDEBUG -optc-O3 -optc-finline-functions -optc-fomit-frame-pointer -optc-fno-strict-aliasing -optc--param -optcmax-inline-insns-single=1800 -#include "SFMT.h" -package base-3.0.1.0 -package old-time-1.0.0.0 -O -Wall -O2 -fvia-C -XCPP -XForeignFunctionInterface -XBangPatterns Unit.hs -o Unit -fhpc -no-recomp  ../dist/build/cbits/SFMT.o -I../include++clean:+	rm -f *~ Unit *.hi *.o

+ tests/Unit.hs view

@@ -0,0 +1,264 @@+{- LANGUAGE BangPatterns, PatternSignatures, ScopedTypeVariables -}+{-# OPTIONS -fbang-patterns -fglasgow-exts -#include "SFMT.h" #-}++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 System.Random.Mersenne++main = do+    print  version++    g <- newMTGen (Just 5)++    time $ gen g+    ranges_strict g+--    ranges_strict_range g+    ranges_ty g+--    ranges g++    speed (undefined :: Int) g+    speed (undefined :: Integer) g+    speed (undefined :: Double) g++    time $ sums g+    time $ sum_lazy g++gen g = do+    forM_ [0 .. 10000] $ \i -> do+        x <- random g :: IO Word+        when (i < 1000) $ do+            printf "%12u " (fromIntegral x :: Int)+            when (i `rem` 4 == 3) $+                putChar '\n'+++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+++-- overhead cause by random's badness+sums g = do++    let lim = 100000000+    printf "Generating %d randoms ...\n" lim+++    let go :: Int -> Int -> IO Int+        go !n !acc+            | n >= lim = return acc+            | otherwise     = do+                    a <- random g+                    go (n+1) (if a < acc then a else acc)+    print =<< go 0 0+++sum_lazy g = do+    let lim = 100000000+    printf "Generating %d randoms lazily ...\n" lim++    xs <- randoms g :: IO [Int]+    print (minimum (take lim xs)) -- faster when it fuses.++++{-+-- overhead cause by random's badness+ranges g = do+    let n = 10000000+    test g n (undefined :: Bool)+    test g n (undefined :: Word8)+    test g n (undefined :: Word16)+    test g n (undefined :: Word32)+    test g n (undefined :: Word64)+    test g n (undefined :: Word)+    test g n (undefined :: Int)+    test g n (undefined :: Int64)+    test g n (undefined :: Int32)+    test g n (undefined :: Int16)+    test g n (undefined :: Int8)+    test g n (undefined :: Integer)+-}++-- overhead cause by random's badness+ranges_ty g = do+    let n = 10000000+    test_type g n (undefined :: Bool)+    test_type g n (undefined :: Word8)+    test_type g n (undefined :: Word16)+    test_type g n (undefined :: Word32)+    test_type g n (undefined :: Word64)+    test_type g n (undefined :: Word)+    test_type g n (undefined :: Int)+    test_type g n (undefined :: Int64)+    test_type g n (undefined :: Int32)+    test_type g n (undefined :: Int16)+    test_type g n (undefined :: Int8)++-- overhead cause by random's badness+ranges_strict g = do+    let n = 100000000+    test_strict g n (undefined :: Word)+    test_strict g n (undefined :: Word8)+    test_strict g n (undefined :: Word16)+    test_strict g n (undefined :: Word32)+    test_strict g n (undefined :: Word64)+    test_strict g n (undefined :: Int)+    test_strict g n (undefined :: Int64)+    test_strict g n (undefined :: Int32)+    test_strict g n (undefined :: Int16)+    test_strict g n (undefined :: Int8)+    test_strict g n (undefined :: Double)+    test_strict g n (undefined :: Integer)+    test_strict g n (undefined :: Bool)++{-+-- overhead cause by random's badness+ranges_strict_range g = do+    let n = 100000000+    test_strict_range g n (undefined :: Word)+    test_strict_range g n (undefined :: Word8)+    test_strict_range g n (undefined :: Word16)+    test_strict_range g n (undefined :: Word32)+    test_strict_range g n (undefined :: Word64)+    test_strict_range g n (undefined :: Int)+    test_strict_range g n (undefined :: Int64)+    test_strict_range g n (undefined :: Int32)+    test_strict_range g n (undefined :: Int16)+    test_strict_range g n (undefined :: Int8)+    test_strict_range g n (undefined :: Double)+    test_strict_range g n (undefined :: Integer)+    test_strict_range g n (undefined :: Bool)+    -}+++------------------------------------------------------------------------+-- check values are in range for randomRs++{-+test :: forall a . (Show a, Ord a, Typeable a, MTRandom a) => MTGen -> Int -> a -> IO ()+test g n ty = do+    a' <- random g :: IO a+    b' <- random g :: IO a+    let (a,b) = (min a' b', max a' b')++    printf "%d bounded :: %s ...\t" n (show $  typeOf ty)+    hFlush stdout++    time $ do+        xs <- randomRs (a,b) g :: IO [a]++        sequence_+            [ if x >= a && x <= b+                    then return ()+                    else error $ "Fail " ++ show (x,a,b)+            | x <- take n xs ]++--      printf " all good."++-}++test_type :: forall a . (Bounded a, Show a, Ord a, Typeable a, MTRandom a) => MTGen -> Int -> a -> IO ()+test_type g n ty = do++    printf "lazy   generation of %d  :: %s ...\t" n (show $  typeOf ty)+    hFlush stdout++    time $ do+        xs <- randoms g :: IO [a]++        sequence_+            [ if x >= minBound && x <= maxBound+                    then return ()+                    else error $ "Fail " ++ show x+            | x <- take n xs ]++--      printf "all good. "++test_strict :: forall a . (Show a, Ord a, Typeable a, MTRandom a) => MTGen -> Int -> a -> IO ()+test_strict g n ty = do++    printf "strict generation of %d :: type %s ...\t" n (show $  typeOf ty)+    hFlush stdout++    time $ do+        let go i+                | i > n = return ()+                | otherwise = do+                    x <- random g :: IO a+                    x `seq` go (i+1)++        go 0++--      printf "all good. "++{-+test_strict_range :: forall a . (Show a, Ord a, Typeable a, MTRandom a) => MTGen -> Int -> a -> IO ()+test_strict_range g n ty = do+    a' <- random g :: IO a+    b' <- random g :: IO a+    let (a,b) = (min a' b', max a' b')++    printf "strict, ranged generation of %d :: type %s ...\t" n (show $  typeOf ty)+    hFlush stdout++    time $ do+        let go i+                | i > n = return ()+                | otherwise = do+                    x <- randomR (a,b) g :: IO a+                    if x >= a && x <= b+                        then go (i+1)+                        else error $ "test_strict_range failed " ++ show (a,b,x)++        go 0++        printf "all good. "+-}++------------------------------------------------------------------------+-- compare with System.Random++-- overhead cause by random's badness+speed :: forall a . (Show a, Ord a, Typeable a, Old.Random a, Num a, MTRandom a)+                            => a -> MTGen -> IO ()+speed ty g  = do++    let lim = 5000000++    time $ do+        putStrLn $ "System.Random: " ++ show lim ++ " " ++ show (show $ typeOf ty)+        let g = Old.mkStdGen 5++        let go :: Old.StdGen -> Int -> a -> a+            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: " ++ show lim ++ " " ++ show (show $ typeOf ty)+        let go !n !acc+                | n >= lim = return acc+                | otherwise     = do+                        a <- random g :: IO a+                        go (n+1::Int) (if a > acc then a else acc)+        print =<< go 0 0++

+ tests/number-fast.hs view

@@ -0,0 +1,13 @@++import System.Environment+import System.Random.Mersenne+import Control.Monad+import Data.Char++main = do+    [i, n] <- map read `fmap` getArgs+    g <- newMTGen Nothing+    replicateM_ i $ do+          j <- randomR (1,n) g+          putStrLn . map chr . take j =<< randomRs (ord 'a', ord 'z') g+

+ tests/pi.hs view

@@ -0,0 +1,52 @@+{-# OPTIONS -cpp -fglasgow-exts -fvia-C -optc-O2 -optc-msse2 -optc-march=pentium4 -optc-ffast-math -fexcess-precision #-}++#if defined(FAST)+import System.Random.Mersenne+#else+import System.Random+#endif+import System.Environment++#if defined(FAST)++main = do+    [lim] <- mapM readIO =<< getArgs+    g <- newMTGen Nothing++    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++    print =<< go 0 0++{-+    $ time ./pi 100000000+    3.14166916+    ./pi 100000000  4.31s user 0.00s system 99% cpu 4.311 total+-}++#else++main = do+    [lim] <- mapM readIO =<< getArgs+    g <- newStdGen++    let go :: StdGen -> Int -> Int -> IO Double+        go g throws ins+            | throws >= lim  = return ((4 * fromIntegral ins) / (fromIntegral throws))+            | otherwise = do+                let (x::Double,g1) = random g+                    (y::Double,g2) = random g1+                if x * x + y * y < 1+                    then go g2 (throws+1) $! ins + 1+                    else go g2 (throws+1) ins++    print =<< go g 0 0++#endif

+ tests/runtests.sh view

@@ -0,0 +1,5 @@+#!/bin/sh++echo "Testing C ... "+gcc sum.c -o sum-c -I../include ../dist/build/cbits/SFMT.o+/usr/bin/time ./sum-c

+ tests/sum.c view

@@ -0,0 +1,17 @@+#include "SFMT.h"++int main() {+    long i;+    long n, acc = 0;+    long lim =100000000; +    printf("Generating %d randoms ...\n", lim);++    init_gen_rand(5);+    for (i=0;  i < lim; i++) {+        n = gen_rand64();+        if (n < acc) { acc = n; }+    }+    printf("%ld\n", acc);++    return 0;+}