diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/LICENSE.sfmt b/LICENSE.sfmt
new file mode 100644
--- /dev/null
+++ b/LICENSE.sfmt
@@ -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.
diff --git a/README b/README
new file mode 100644
--- /dev/null
+++ b/README
@@ -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
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/System/Random/Mersenne.hs b/System/Random/Mersenne.hs
new file mode 100644
--- /dev/null
+++ b/System/Random/Mersenne.hs
@@ -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!
+--
diff --git a/TODO b/TODO
new file mode 100644
--- /dev/null
+++ b/TODO
@@ -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
+
diff --git a/cbits/SFMT.c b/cbits/SFMT.c
new file mode 100644
--- /dev/null
+++ b/cbits/SFMT.c
@@ -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;
+}
diff --git a/include/SFMT-alti.h b/include/SFMT-alti.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-alti.h
@@ -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
diff --git a/include/SFMT-params.h b/include/SFMT-params.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params.h
@@ -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 */
diff --git a/include/SFMT-params11213.h b/include/SFMT-params11213.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params11213.h
@@ -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 */
diff --git a/include/SFMT-params1279.h b/include/SFMT-params1279.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params1279.h
@@ -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 */
diff --git a/include/SFMT-params132049.h b/include/SFMT-params132049.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params132049.h
@@ -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 */
diff --git a/include/SFMT-params19937.h b/include/SFMT-params19937.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params19937.h
@@ -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 */
diff --git a/include/SFMT-params216091.h b/include/SFMT-params216091.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params216091.h
@@ -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 */
diff --git a/include/SFMT-params2281.h b/include/SFMT-params2281.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params2281.h
@@ -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 */
diff --git a/include/SFMT-params4253.h b/include/SFMT-params4253.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params4253.h
@@ -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 */
diff --git a/include/SFMT-params44497.h b/include/SFMT-params44497.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params44497.h
@@ -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 */
diff --git a/include/SFMT-params607.h b/include/SFMT-params607.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params607.h
@@ -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 */
diff --git a/include/SFMT-params86243.h b/include/SFMT-params86243.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-params86243.h
@@ -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 */
diff --git a/include/SFMT-sse2.h b/include/SFMT-sse2.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT-sse2.h
@@ -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
diff --git a/include/SFMT.h b/include/SFMT.h
new file mode 100644
--- /dev/null
+++ b/include/SFMT.h
@@ -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
diff --git a/mersenne-random.cabal b/mersenne-random.cabal
new file mode 100644
--- /dev/null
+++ b/mersenne-random.cabal
@@ -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
diff --git a/tests/Makefile b/tests/Makefile
new file mode 100644
--- /dev/null
+++ b/tests/Makefile
@@ -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
diff --git a/tests/Unit.hs b/tests/Unit.hs
new file mode 100644
--- /dev/null
+++ b/tests/Unit.hs
@@ -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
+
+
diff --git a/tests/number-fast.hs b/tests/number-fast.hs
new file mode 100644
--- /dev/null
+++ b/tests/number-fast.hs
@@ -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
+
diff --git a/tests/pi.hs b/tests/pi.hs
new file mode 100644
--- /dev/null
+++ b/tests/pi.hs
@@ -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
diff --git a/tests/runtests.sh b/tests/runtests.sh
new file mode 100644
--- /dev/null
+++ b/tests/runtests.sh
@@ -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
diff --git a/tests/sum.c b/tests/sum.c
new file mode 100644
--- /dev/null
+++ b/tests/sum.c
@@ -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;
+}