Added Random.SFMT

.gitignore

@@ -0,0 +1,9 @@
+.bmx
+
+*.bak
+
+*.i
+
+*.i2
+
+*.a

sfmt.mod/LICENSE.txt

@@ -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.

sfmt.mod/common.bmx

@@ -0,0 +1,52 @@
+' Copyright (c) 2007-2020, Bruce A Henderson
+' 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 Bruce A Henderson 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 Bruce A Henderson ``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 Bruce A Henderson 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.
+'
+SuperStrict
+
+Import "src/*.h"
+Import "src/SFMT.c"
+
+
+Import "sfmtglue.c"
+
+Extern
+	Function bmx_sfmt_init_gen_rand:Byte Ptr(handle:Byte Ptr, seed:Int)
+	Function bmx_sfmt_free(handle:Byte Ptr)
+?bmxng
+	Function sfmt_genrand_uint32:UInt(handle:Byte Ptr)
+	Function bmx_gen_rand64(handle:Byte Ptr, v:ULong Ptr)
+?Not bmxng
+	Function sfmt_genrand_uint32:Int(handle:Byte Ptr)
+	Function bmx_gen_rand64(handle:Byte Ptr, v:Long Ptr)
+?
+	
+	Function bmx_genrand_real1:Double(handle:Byte Ptr)
+	Function bmx_genrand_real2:Double(handle:Byte Ptr)
+	Function bmx_genrand_real3:Double(handle:Byte Ptr)
+	Function bmx_genrand_res53:Double(handle:Byte Ptr)
+	
+End Extern
+
+

sfmt.mod/sfmt.bmx

@@ -0,0 +1,154 @@
+' Copyright (c) 2007-2020, Bruce A Henderson
+' 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 Bruce A Henderson 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 Bruce A Henderson ``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 Bruce A Henderson 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.
+'
+SuperStrict
+
+Rem
+bbdoc: Random Numbers - SFMT
+End Rem
+Module Random.SFMT
+
+ModuleInfo "Version: 1.06"
+ModuleInfo "License: BSD"
+ModuleInfo "Copyright: SFMT - 2006-2017 Mutsuo Saito, Makoto Matsumoto and Hiroshima"
+ModuleInfo "Copyright: Wrapper - 2007-2020 Bruce A Henderson"
+ModuleInfo "Modserver: BRL"
+
+ModuleInfo "History: 1.06"
+ModuleInfo "History: Integrated into new BRL.Random API."
+ModuleInfo "History: 1.05"
+ModuleInfo "History: Created TRandom to support multiple generators."
+ModuleInfo "History: Aligned pointer creation."
+ModuleInfo "History: 1.04"
+ModuleInfo "History: Update to SFMT 1.5.1"
+ModuleInfo "History: Refactored."
+ModuleInfo "History: 1.03"
+ModuleInfo "History: Update to SFMT 1.4.1"
+ModuleInfo "History: Updated for NG."
+ModuleInfo "History: 1.02"
+ModuleInfo "History: Fix for PPC Mac compile."
+ModuleInfo "History: 1.01"
+ModuleInfo "History: Automatically initializes via SeedRnd() if required."
+ModuleInfo "History: 1.00"
+ModuleInfo "History: Initial Version (SFMT 1.2)"
+
+ModuleInfo "CC_OPTS: -DMEXP=19937"
+ModuleInfo "CC_OPTS: -fno-strict-aliasing"
+ModuleInfo "CC_OPTS: -std=c11"
+
+Import BRL.Random
+
+?x86
+ModuleInfo "CC_OPTS: -msse2 -DHAVE_SSE2"
+?x64
+ModuleInfo "CC_OPTS: -msse2 -DHAVE_SSE2"
+?
+
+Import "common.bmx"
+
+Rem
+bbdoc: An instance of a random number generator.
+End Rem
+Type TSFMTRandom Extends TRandom
+
+	Field sfmtPtr:Byte Ptr
+	Field rnd_seed:Int
+	
+	Method New()
+		sfmtPtr = bmx_sfmt_init_gen_rand(Null, GenerateSeed())
+	End Method
+
+	Method New(seed:Int)
+		rnd_seed = seed
+		sfmtPtr = bmx_sfmt_init_gen_rand(Null, seed)
+	End Method
+	
+	Method SeedRnd(seed:Int)
+		rnd_seed = seed
+		bmx_sfmt_init_gen_rand(sfmtPtr, seed)
+	End Method
+
+	Method RndSeed:Int()
+		Return rnd_seed
+	End Method
+
+	Method Rand:Int( min_value:Int, max_value:Int = 1 )
+		Local Range:Double = max_value - min_value
+		If Range > 0 Return Int( bmx_genrand_res53(sfmtPtr)*(1:Double+Range) )+min_value
+		Return Int( bmx_genrand_res53(sfmtPtr)*(1:Double-Range) )+max_value
+	End Method
+
+	Method Rnd:Double( min_value!=1,max_value!=0 )
+		If max_value > min_value Return RndDouble() * (max_value - min_value) + min_value
+		Return RndDouble() * (min_value - max_value) + max_value
+	End Method
+
+	Method Rand64:Long( min_value:Long, max_value:Long = 1 )
+		Local Range:Long = max_value - min_value
+		If Range > 0 Return Long( bmx_genrand_res53(sfmtPtr) * (1:Long + Range) ) + min_value
+		Return Long( bmx_genrand_res53(sfmtPtr) * (1:Long - Range) ) + max_value
+	End Method
+	
+	Method RndFloat:Float()
+		Return Float(bmx_genrand_real3(sfmtPtr))
+	End Method
+	
+	Method RndDouble:Double()
+		Return bmx_genrand_res53(sfmtPtr)
+	End Method
+	
+	Method Delete()
+		If sfmtPtr Then
+			bmx_sfmt_free(sfmtPtr)
+			sfmtPtr = Null
+		End If
+	End Method
+	
+End Type
+
+
+Private
+Type TSFMTRandomFactory Extends TRandomFactory
+	
+	Method New()
+		Super.New()
+		Init()
+	End Method
+	
+	Method GetName:String()
+		Return "SFMTRandom"
+	End Method
+	
+	Method Create:TRandom(seed:Int)
+		Return New TSFMTRandom(seed)
+	End Method
+
+	Method Create:TRandom()
+		Return New TSFMTRandom()
+	End Method
+		
+End Type
+
+New TSFMTRandomFactory

sfmt.mod/sfmtglue.c

@@ -0,0 +1,78 @@
+/*
+ Copyright (c) 2007-2020, Bruce A Henderson
+ 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 Bruce A Henderson 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 Bruce A Henderson ``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 Bruce A Henderson 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.
+*/
+
+
+#include "brl.mod/blitz.mod/blitz.h"
+#include <stdlib.h>
+#include <malloc.h>
+#include "SFMT.h"
+
+sfmt_t * bmx_sfmt_init_gen_rand(sfmt_t * sfmt, int seed) {
+	if (sfmt == NULL) {
+#ifdef WIN32
+		sfmt = _aligned_malloc(sizeof(*sfmt), 16);
+#else
+		sfmt = aligned_malloc(16, sizeof(*sfmt));
+#endif
+	}
+	sfmt_init_gen_rand(sfmt, seed);
+	return sfmt;
+}
+
+void bmx_sfmt_free(sfmt_t * sfmt) {
+#ifdef WIN32
+	_aligned_free(sfmt);
+#else
+	free(sfmt);
+#endif
+}
+
+double bmx_genrand_real1(sfmt_t * sfmt) {
+	return sfmt_genrand_real1(sfmt);
+}
+
+double bmx_genrand_real2(sfmt_t * sfmt) {
+	return sfmt_genrand_real2(sfmt);
+}
+
+double bmx_genrand_real3(sfmt_t * sfmt) {
+	return sfmt_genrand_real3(sfmt);
+}
+
+#ifdef BMX_NG
+void bmx_gen_rand64(sfmt_t * sfmt, BBUInt64 * r) {
+	*r = sfmt_genrand_uint64(sfmt);
+}
+#else
+void bmx_gen_rand64(sfmt_t * sfmt, BBInt64 * r) {
+	*r = sfmt_genrand_uint64(sfmt);
+}
+#endif
+
+double bmx_genrand_res53(sfmt_t * sfmt) {
+    return sfmt_genrand_res53(sfmt);
+}

sfmt.mod/src/SFMT-alti.h

@@ -0,0 +1,156 @@
+#pragma once
+/**
+ * @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);
+
+/**
+ * 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 = SFMT_ALTI_SL1;
+    const vector unsigned int sr1 = SFMT_ALTI_SR1;
+#ifdef ONLY64
+    const vector unsigned int mask = SFMT_ALTI_MSK64;
+    const vector unsigned char perm_sl = SFMT_ALTI_SL2_PERM64;
+    const vector unsigned char perm_sr = SFMT_ALTI_SR2_PERM64;
+#else
+    const vector unsigned int mask = SFMT_ALTI_MSK;
+    const vector unsigned char perm_sl = SFMT_ALTI_SL2_PERM;
+    const vector unsigned char perm_sr = SFMT_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.
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    vector unsigned int r, r1, r2;
+
+    r1 = sfmt->state[N - 2].s;
+    r2 = sfmt->state[N - 1].s;
+    for (i = 0; i < N - POS1; i++) {
+	r = vec_recursion(sfmt->state[i].s, sfmt->state[i + POS1].s, r1, r2);
+	sfmt->state[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+    for (; i < N; i++) {
+	r = vec_recursion(sfmt->state[i].s, sfmt->state[i + POS1 - N].s, r1, r2);
+	sfmt->state[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(sfmt_t * sfmt, w128_t *array, int size) {
+    int i, j;
+    vector unsigned int r, r1, r2;
+
+    r1 = sfmt->state[N - 2].s;
+    r2 = sfmt->state[N - 1].s;
+    for (i = 0; i < N - POS1; i++) {
+	r = vec_recursion(sfmt->state[i].s, sfmt->state[i + POS1].s, r1, r2);
+	array[i].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+    for (; i < N; i++) {
+	r = vec_recursion(sfmt->state[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->state[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->state[j++].s = r;
+	r1 = r2;
+	r2 = r;
+    }
+}
+
+#ifndef ONLY64
+#if defined(__APPLE__)
+#define SFMT_ALTI_SWAP (vector unsigned char) \
+	(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11)
+#else
+#define SFMT_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 = SFMT_ALTI_SWAP;
+
+    for (i = 0; i < size; i++) {
+	array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);
+    }
+}
+#endif
+
+#endif

sfmt.mod/src/SFMT-common.h

@@ -0,0 +1,164 @@
+#pragma once
+/**
+ * @file SFMT-common.h
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
+ * number generator with jump function. This file includes common functions
+ * used in random number generation and jump.
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (The University of Tokyo)
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+ * University and The University of Tokyo.
+ * All rights reserved.
+ *
+ * The 3-clause BSD License is applied to this software, see
+ * LICENSE.txt
+ */
+#ifndef SFMT_COMMON_H
+#define SFMT_COMMON_H
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+#include "SFMT.h"
+
+inline static void do_recursion(w128_t * r, w128_t * a, w128_t * b,
+                                w128_t * c, w128_t * d);
+
+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);
+
+/**
+ * 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
+ */
+#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, SFMT_SL2);
+    rshift128(&y, c, SFMT_SR2);
+    r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SFMT_SR1) & SFMT_MSK2) ^ y.u[0]
+        ^ (d->u[0] << SFMT_SL1);
+    r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SFMT_SR1) & SFMT_MSK1) ^ y.u[1]
+        ^ (d->u[1] << SFMT_SL1);
+    r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SFMT_SR1) & SFMT_MSK4) ^ y.u[2]
+        ^ (d->u[2] << SFMT_SL1);
+    r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SFMT_SR1) & SFMT_MSK3) ^ y.u[3]
+        ^ (d->u[3] << SFMT_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, SFMT_SL2);
+    rshift128(&y, c, SFMT_SR2);
+    r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SFMT_SR1) & SFMT_MSK1)
+        ^ y.u[0] ^ (d->u[0] << SFMT_SL1);
+    r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SFMT_SR1) & SFMT_MSK2)
+        ^ y.u[1] ^ (d->u[1] << SFMT_SL1);
+    r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SFMT_SR1) & SFMT_MSK3)
+        ^ y.u[2] ^ (d->u[2] << SFMT_SL1);
+    r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SFMT_SR1) & SFMT_MSK4)
+        ^ y.u[3] ^ (d->u[3] << SFMT_SL1);
+}
+#endif
+#if defined(__cplusplus)
+}
+#endif
+
+#endif // SFMT_COMMON_H

sfmt.mod/src/SFMT-neon.h

@@ -0,0 +1,111 @@
+/**
+ * @file  SFMT-neon.h
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) for ARM with 128b NEON 
+ *
+ * @author Masaki Ota 
+ *
+ * @note We assume LITTLE ENDIAN in this file
+ */
+
+#ifndef SFMT_NEON_H
+#define SFMT_NEON_H
+
+inline static void neon_recursion(uint32x4_t * r, uint32x4_t a, uint32x4_t b,
+                                uint32x4_t c, uint32x4_t d);
+
+
+
+/**
+ * This function represents the recursion formula.
+ * @param r an output
+ * @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
+ */
+inline static void neon_recursion(uint32x4_t * r, uint32x4_t a, uint32x4_t b,
+                                uint32x4_t c, uint32x4_t d)
+{
+    uint32x4_t v, x, y, z;
+    static const uint32x4_t vzero = {0,0,0,0};
+    static const uint32x4_t vmask = {SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4};
+
+#define rotate_bytes(A, B, C) vreinterpretq_u32_u8(vextq_u8(vreinterpretq_u8_u32(A),vreinterpretq_u8_u32(B),(C)))
+
+    y = vshrq_n_u32(b, SFMT_SR1);
+    z = rotate_bytes(c, vzero, SFMT_SR2);
+    v = vshlq_n_u32(d, SFMT_SL1);
+    z = veorq_u32(z, a);
+    z = veorq_u32(z, v);
+    x = rotate_bytes(vzero, a, 16-SFMT_SL2);
+    y = vandq_u32(y, vmask);
+    z = veorq_u32(z, x);
+    z = veorq_u32(z, y);
+    *r = z;
+}
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    uint32x4_t r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        neon_recursion(&pstate[i].si, pstate[i].si, pstate[i + SFMT_POS1].si, r1, r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        neon_recursion(&pstate[i].si, pstate[i].si, pstate[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
+{
+    int i, j;
+    uint32x4_t r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        neon_recursion(&array[i].si, pstate[i].si, pstate[i + SFMT_POS1].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        neon_recursion(&array[i].si, pstate[i].si, array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < size - SFMT_N; i++) {
+        neon_recursion(&array[i].si, array[i - SFMT_N].si, array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (j = 0; j < 2 * SFMT_N - size; j++) {
+        pstate[j] = array[j + size - SFMT_N];
+    }
+    for (; i < size; i++, j++) {
+        neon_recursion(&array[i].si, array[i - SFMT_N].si, array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = pstate[j].si = array[i].si;
+    }
+}
+
+#endif

sfmt.mod/src/SFMT-params.h

@@ -0,0 +1,98 @@
+#pragma once
+#ifndef SFMT_PARAMS_H
+#define SFMT_PARAMS_H
+
+#if !defined(SFMT_MEXP)
+#if defined(__GNUC__) && !defined(__ICC)
+  #warning "SFMT_MEXP is not defined. I assume MEXP is 19937."
+#endif
+  #define SFMT_MEXP 19937
+#endif
+/*-----------------
+  BASIC DEFINITIONS
+  -----------------*/
+/** Mersenne Exponent. The period of the sequence
+ *  is a multiple of 2^MEXP-1.
+ * #define SFMT_MEXP 19937 */
+/** SFMT generator has an internal state array of 128-bit integers,
+ * and N is its size. */
+#define SFMT_N (SFMT_MEXP / 128 + 1)
+/** N32 is the size of internal state array when regarded as an array
+ * of 32-bit integers.*/
+#define SFMT_N32 (SFMT_N * 4)
+/** N64 is the size of internal state array when regarded as an array
+ * of 64-bit integers.*/
+#define SFMT_N64 (SFMT_N * 2)
+
+/*----------------------
+  the parameters of SFMT
+  following definitions are in paramsXXXX.h file.
+  ----------------------*/
+/** the pick up position of the array.
+#define SFMT_POS1 122
+*/
+
+/** the parameter of shift left as four 32-bit registers.
+#define SFMT_SL1 18
+ */
+
+/** the parameter of shift left as one 128-bit register.
+ * The 128-bit integer is shifted by (SFMT_SL2 * 8) bits.
+#define SFMT_SL2 1
+*/
+
+/** the parameter of shift right as four 32-bit registers.
+#define SFMT_SR1 11
+*/
+
+/** the parameter of shift right as one 128-bit register.
+ * The 128-bit integer is shifted by (SFMT_SR2 * 8) bits.
+#define SFMT_SR2 1
+*/
+
+/** A bitmask, used in the recursion.  These parameters are introduced
+ * to break symmetry of SIMD.
+#define SFMT_MSK1 0xdfffffefU
+#define SFMT_MSK2 0xddfecb7fU
+#define SFMT_MSK3 0xbffaffffU
+#define SFMT_MSK4 0xbffffff6U
+*/
+
+/** These definitions are part of a 128-bit period certification vector.
+#define SFMT_PARITY1    0x00000001U
+#define SFMT_PARITY2    0x00000000U
+#define SFMT_PARITY3    0x00000000U
+#define SFMT_PARITY4    0xc98e126aU
+*/
+
+#if SFMT_MEXP == 607
+  #include "SFMT-params607.h"
+#elif SFMT_MEXP == 1279
+  #include "SFMT-params1279.h"
+#elif SFMT_MEXP == 2281
+  #include "SFMT-params2281.h"
+#elif SFMT_MEXP == 4253
+  #include "SFMT-params4253.h"
+#elif SFMT_MEXP == 11213
+  #include "SFMT-params11213.h"
+#elif SFMT_MEXP == 19937
+  #include "SFMT-params19937.h"
+#elif SFMT_MEXP == 44497
+  #include "SFMT-params44497.h"
+#elif SFMT_MEXP == 86243
+  #include "SFMT-params86243.h"
+#elif SFMT_MEXP == 132049
+  #include "SFMT-params132049.h"
+#elif SFMT_MEXP == 216091
+  #include "SFMT-params216091.h"
+#else
+#if defined(__GNUC__) && !defined(__ICC)
+  #error "SFMT_MEXP is not valid."
+  #undef SFMT_MEXP
+#else
+  #undef SFMT_MEXP
+#endif
+
+#endif
+
+#endif /* SFMT_PARAMS_H */

sfmt.mod/src/SFMT-params11213.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS11213_H
+#define SFMT_PARAMS11213_H
+
+#define SFMT_POS1	68
+#define SFMT_SL1	14
+#define SFMT_SL2	3
+#define SFMT_SR1	7
+#define SFMT_SR2	3
+#define SFMT_MSK1	0xeffff7fbU
+#define SFMT_MSK2	0xffffffefU
+#define SFMT_MSK3	0xdfdfbfffU
+#define SFMT_MSK4	0x7fffdbfdU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0xe8148000U
+#define SFMT_PARITY4	0xd0c7afa3U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
+    #define SFMT_ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
+    #define SFMT_ALTI_SR2_PERM	{5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}
+    #define SFMT_ALTI_SR2_PERM64 {13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-11213:68-14-3-7-3:effff7fb-ffffffef-dfdfbfff-7fffdbfd"
+
+#endif /* SFMT_PARAMS11213_H */

sfmt.mod/src/SFMT-params1279.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS1279_H
+#define SFMT_PARAMS1279_H
+
+#define SFMT_POS1	7
+#define SFMT_SL1	14
+#define SFMT_SL2	3
+#define SFMT_SR1	5
+#define SFMT_SR2	1
+#define SFMT_MSK1	0xf7fefffdU
+#define SFMT_MSK2	0x7fefcfffU
+#define SFMT_MSK3	0xaff3ef3fU
+#define SFMT_MSK4	0xb5ffff7fU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0x00000000U
+#define SFMT_PARITY4	0x20000000U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
+    #define SFMT_ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-1279:7-14-3-5-1:f7fefffd-7fefcfff-aff3ef3f-b5ffff7f"
+
+#endif /* SFMT_PARAMS1279_H */

sfmt.mod/src/SFMT-params132049.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS132049_H
+#define SFMT_PARAMS132049_H
+
+#define SFMT_POS1	110
+#define SFMT_SL1	19
+#define SFMT_SL2	1
+#define SFMT_SR1	21
+#define SFMT_SR2	1
+#define SFMT_MSK1	0xffffbb5fU
+#define SFMT_MSK2	0xfb6ebf95U
+#define SFMT_MSK3	0xfffefffaU
+#define SFMT_MSK4	0xcff77fffU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0xcb520000U
+#define SFMT_PARITY4	0xc7e91c7dU
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
+    #define SFMT_ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-132049:110-19-1-21-1:ffffbb5f-fb6ebf95-fffefffa-cff77fff"
+
+#endif /* SFMT_PARAMS132049_H */

sfmt.mod/src/SFMT-params19937.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS19937_H
+#define SFMT_PARAMS19937_H
+
+#define SFMT_POS1	122
+#define SFMT_SL1	18
+#define SFMT_SL2	1
+#define SFMT_SR1	11
+#define SFMT_SR2	1
+#define SFMT_MSK1	0xdfffffefU
+#define SFMT_MSK2	0xddfecb7fU
+#define SFMT_MSK3	0xbffaffffU
+#define SFMT_MSK4	0xbffffff6U
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0x00000000U
+#define SFMT_PARITY4	0x13c9e684U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
+    #define SFMT_ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6"
+
+#endif /* SFMT_PARAMS19937_H */

sfmt.mod/src/SFMT-params216091.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS216091_H
+#define SFMT_PARAMS216091_H
+
+#define SFMT_POS1	627
+#define SFMT_SL1	11
+#define SFMT_SL2	3
+#define SFMT_SR1	10
+#define SFMT_SR2	1
+#define SFMT_MSK1	0xbff7bff7U
+#define SFMT_MSK2	0xbfffffffU
+#define SFMT_MSK3	0xbffffa7fU
+#define SFMT_MSK4	0xffddfbfbU
+#define SFMT_PARITY1	0xf8000001U
+#define SFMT_PARITY2	0x89e80709U
+#define SFMT_PARITY3	0x3bd2b64bU
+#define SFMT_PARITY4	0x0c64b1e4U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
+    #define SFMT_ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-216091:627-11-3-10-1:bff7bff7-bfffffff-bffffa7f-ffddfbfb"
+
+#endif /* SFMT_PARAMS216091_H */

sfmt.mod/src/SFMT-params2281.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS2281_H
+#define SFMT_PARAMS2281_H
+
+#define SFMT_POS1	12
+#define SFMT_SL1	19
+#define SFMT_SL2	1
+#define SFMT_SR1	5
+#define SFMT_SR2	1
+#define SFMT_MSK1	0xbff7ffbfU
+#define SFMT_MSK2	0xfdfffffeU
+#define SFMT_MSK3	0xf7ffef7fU
+#define SFMT_MSK4	0xf2f7cbbfU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0x00000000U
+#define SFMT_PARITY4	0x41dfa600U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
+    #define SFMT_ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-2281:12-19-1-5-1:bff7ffbf-fdfffffe-f7ffef7f-f2f7cbbf"
+
+#endif /* SFMT_PARAMS2281_H */

sfmt.mod/src/SFMT-params4253.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS4253_H
+#define SFMT_PARAMS4253_H
+
+#define SFMT_POS1	17
+#define SFMT_SL1	20
+#define SFMT_SL2	1
+#define SFMT_SR1	7
+#define SFMT_SR2	1
+#define SFMT_MSK1	0x9f7bffffU
+#define SFMT_MSK2	0x9fffff5fU
+#define SFMT_MSK3	0x3efffffbU
+#define SFMT_MSK4	0xfffff7bbU
+#define SFMT_PARITY1	0xa8000001U
+#define SFMT_PARITY2	0xaf5390a3U
+#define SFMT_PARITY3	0xb740b3f8U
+#define SFMT_PARITY4	0x6c11486dU
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
+    #define SFMT_ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-4253:17-20-1-7-1:9f7bffff-9fffff5f-3efffffb-fffff7bb"
+
+#endif /* SFMT_PARAMS4253_H */

sfmt.mod/src/SFMT-params44497.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS44497_H
+#define SFMT_PARAMS44497_H
+
+#define SFMT_POS1	330
+#define SFMT_SL1	5
+#define SFMT_SL2	3
+#define SFMT_SR1	9
+#define SFMT_SR2	3
+#define SFMT_MSK1	0xeffffffbU
+#define SFMT_MSK2	0xdfbebfffU
+#define SFMT_MSK3	0xbfbf7befU
+#define SFMT_MSK4	0x9ffd7bffU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0xa3ac4000U
+#define SFMT_PARITY4	0xecc1327aU
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
+    #define SFMT_ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
+    #define SFMT_ALTI_SR2_PERM	{5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}
+    #define SFMT_ALTI_SR2_PERM64 {13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-44497:330-5-3-9-3:effffffb-dfbebfff-bfbf7bef-9ffd7bff"
+
+#endif /* SFMT_PARAMS44497_H */

sfmt.mod/src/SFMT-params607.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS607_H
+#define SFMT_PARAMS607_H
+
+#define SFMT_POS1	2
+#define SFMT_SL1	15
+#define SFMT_SL2	3
+#define SFMT_SR1	13
+#define SFMT_SR2	3
+#define SFMT_MSK1	0xfdff37ffU
+#define SFMT_MSK2	0xef7f3f7dU
+#define SFMT_MSK3	0xff777b7dU
+#define SFMT_MSK4	0x7ff7fb2fU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0x00000000U
+#define SFMT_PARITY4	0x5986f054U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{3,21,21,21,7,0,1,2,11,4,5,6,15,8,9,10}
+    #define SFMT_ALTI_SL2_PERM64 {3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
+    #define SFMT_ALTI_SR2_PERM	{5,6,7,0,9,10,11,4,13,14,15,8,19,19,19,12}
+    #define SFMT_ALTI_SR2_PERM64 {13,14,15,0,1,2,3,4,19,19,19,8,9,10,11,12}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-607:2-15-3-13-3:fdff37ff-ef7f3f7d-ff777b7d-7ff7fb2f"
+
+#endif /* SFMT_PARAMS607_H */

sfmt.mod/src/SFMT-params86243.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS86243_H
+#define SFMT_PARAMS86243_H
+
+#define SFMT_POS1	366
+#define SFMT_SL1	6
+#define SFMT_SL2	7
+#define SFMT_SR1	19
+#define SFMT_SR2	1
+#define SFMT_MSK1	0xfdbffbffU
+#define SFMT_MSK2	0xbff7ff3fU
+#define SFMT_MSK3	0xfd77efffU
+#define SFMT_MSK4	0xbf9ff3ffU
+#define SFMT_PARITY1	0x00000001U
+#define SFMT_PARITY2	0x00000000U
+#define SFMT_PARITY3	0x00000000U
+#define SFMT_PARITY4	0xe9528d85U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define SFMT_ALTI_SL1 \
+	(vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+    #define SFMT_ALTI_SR1 \
+	(vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+    #define SFMT_ALTI_MSK \
+	(vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+    #define SFMT_ALTI_MSK64 \
+	(vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+    #define SFMT_ALTI_SL2_PERM \
+	(vector unsigned char)(25,25,25,25,3,25,25,25,7,0,1,2,11,4,5,6)
+    #define SFMT_ALTI_SL2_PERM64 \
+	(vector unsigned char)(7,25,25,25,25,25,25,25,15,0,1,2,3,4,5,6)
+    #define SFMT_ALTI_SR2_PERM \
+	(vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+    #define SFMT_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 SFMT_ALTI_SL1	{SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+    #define SFMT_ALTI_SR1	{SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+    #define SFMT_ALTI_MSK	{SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+    #define SFMT_ALTI_MSK64	{SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+    #define SFMT_ALTI_SL2_PERM	{25,25,25,25,3,25,25,25,7,0,1,2,11,4,5,6}
+    #define SFMT_ALTI_SL2_PERM64 {7,25,25,25,25,25,25,25,15,0,1,2,3,4,5,6}
+    #define SFMT_ALTI_SR2_PERM	{7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+    #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif	/* For OSX */
+#define SFMT_IDSTR	"SFMT-86243:366-6-7-19-1:fdbffbff-bff7ff3f-fd77efff-bf9ff3ff"
+
+#endif /* SFMT_PARAMS86243_H */

sfmt.mod/src/SFMT-sse2-msc.h

@@ -0,0 +1,321 @@
+#pragma once
+/**
+ * @file  SFMT-sse2-msc.h
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) for Intel SSE2 for MSC
+ *
+ * @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.
+ * Copyright (C) 2013 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ *
+ * The new BSD License is applied to this software, see LICENSE.txt
+ */
+
+#ifndef SFMT_SSE2_MSC_H
+#define SFMT_SSE2_MSC_H
+
+/* This header file is included only when _MSC_VER is defined. */
+#if _MSC_VER > 1700
+
+inline static __m128i __vectorcall mm_recursion(__m128i a, __m128i b,
+                                                __m128i c, __m128i d);
+
+/**
+ * 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
+ * @return new value
+ */
+inline static __m128i __vectorcall mm_recursion(__m128i a, __m128i b,
+                                                __m128i c, __m128i d)
+{
+    __m128i v, x, y, z;
+
+    y = _mm_srli_epi32(b, SFMT_SR1);
+    z = _mm_srli_si128(c, SFMT_SR2);
+    v = _mm_slli_epi32(d, SFMT_SL1);
+    z = _mm_xor_si128(z, a);
+    z = _mm_xor_si128(z, v);
+    x = _mm_slli_si128(a, SFMT_SL2);
+    y = _mm_and_si128(y, sse2_param_mask.si);
+    z = _mm_xor_si128(z, x);
+    return _mm_xor_si128(z, y);
+}
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        pstate[i].si = mm_recursion(pstate[i].si,
+                                    pstate[i + SFMT_POS1].si, r1, r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        pstate[i].si = mm_recursion(pstate[i].si,
+                                    pstate[i + SFMT_POS1 - SFMT_N].si,
+                                    r1, r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
+{
+    int i, j;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        array[i].si = mm_recursion(pstate[i].si,
+                                   pstate[i + SFMT_POS1].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        array[i].si = mm_recursion(pstate[i].si,
+                                   array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < size - SFMT_N; i++) {
+        array[i].si = mm_recursion(array[i - SFMT_N].si,
+                                   array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (j = 0; j < 2 * SFMT_N - size; j++) {
+        pstate[j] = array[j + size - SFMT_N];
+    }
+    for (; i < size; i++, j++) {
+        array[i].si = mm_recursion(array[i - SFMT_N].si,
+                                   array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+        pstate[j] = array[i];
+    }
+}
+
+#elif defined(SFMT_USE_MACRO_FUNCTION_FOR_MSC)
+/**
+ * This function represents the recursion formula.
+ * @param r an output
+ * @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
+ */
+#define mm_recursion(r, a, b, c, d)           \
+do {                                          \
+    __m128i v, x, y, z;                       \
+                                              \
+    y = _mm_srli_epi32((b), SFMT_SR1);        \
+    z = _mm_srli_si128((c), SFMT_SR2);        \
+    v = _mm_slli_epi32((d), SFMT_SL1);        \
+    z = _mm_xor_si128(z, (a));                \
+    z = _mm_xor_si128(z, v);                  \
+    x = _mm_slli_si128((a), SFMT_SL2);        \
+    y = _mm_and_si128(y, sse2_param_mask.si); \
+    z = _mm_xor_si128(z, x);                  \
+    r = _mm_xor_si128(z, y);                  \
+} while (0)
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        mm_recursion(pstate[i].si, pstate[i].si,
+                     pstate[i + SFMT_POS1].si, r1, r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        mm_recursion(pstate[i].si, pstate[i].si,
+                     pstate[i + SFMT_POS1 - SFMT_N].si,
+                     r1, r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
+{
+    int i, j;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        mm_recursion(array[i].si, pstate[i].si,
+                     pstate[i + SFMT_POS1].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        mm_recursion(array[i].si, pstate[i].si,
+                     array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < size - SFMT_N; i++) {
+        mm_recursion(array[i].si, array[i - SFMT_N].si,
+                     array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (j = 0; j < 2 * SFMT_N - size; j++) {
+        pstate[j] = array[j + size - SFMT_N];
+    }
+    for (; i < size; i++, j++) {
+        mm_recursion(array[i].si, array[i - SFMT_N].si,
+                     array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+        r1 = r2;
+        r2 = array[i].si;
+        pstate[j] = array[i];
+    }
+}
+#else
+inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
+                                __m128i c, __m128i * d);
+
+/**
+ * This function represents the recursion formula.
+ * @param r an output
+ * @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
+ */
+inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
+                                __m128i c, __m128i * d)
+{
+    __m128i v, x, y, z;
+
+    y = _mm_srli_epi32(b, SFMT_SR1);
+    z = _mm_srli_si128(c, SFMT_SR2);
+    v = _mm_slli_epi32(*d, SFMT_SL1);
+    z = _mm_xor_si128(z, a);
+    z = _mm_xor_si128(z, v);
+    x = _mm_slli_si128(a, SFMT_SL2);
+    y = _mm_and_si128(y, sse2_param_mask.si);
+    z = _mm_xor_si128(z, x);
+    z = _mm_xor_si128(z, y);
+    *r = z;
+}
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        mm_recursion(&pstate[i].si, pstate[i].si,
+                     pstate[i + SFMT_POS1].si, r1, &r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        mm_recursion(&pstate[i].si, pstate[i].si,
+                     pstate[i + SFMT_POS1 - SFMT_N].si,
+                     r1, &r2);
+        r1 = r2;
+        r2 = pstate[i].si;
+    }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
+{
+    int i, j;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        mm_recursion(&array[i].si, pstate[i].si,
+                     pstate[i + SFMT_POS1].si, r1, &r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+        mm_recursion(&array[i].si, pstate[i].si,
+                     array[i + SFMT_POS1 - SFMT_N].si, r1, &r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (; i < size - SFMT_N; i++) {
+        mm_recursion(&array[i].si, array[i - SFMT_N].si,
+                     array[i + SFMT_POS1 - SFMT_N].si, r1, &r2);
+        r1 = r2;
+        r2 = array[i].si;
+    }
+    for (j = 0; j < 2 * SFMT_N - size; j++) {
+        pstate[j] = array[j + size - SFMT_N];
+    }
+    for (; i < size; i++, j++) {
+        mm_recursion(&array[i].si, array[i - SFMT_N].si,
+                     array[i + SFMT_POS1 - SFMT_N].si, r1, &r2);
+        r1 = r2;
+        r2 = array[i].si;
+        pstate[j] = array[i];
+    }
+}
+#endif
+
+#endif

sfmt.mod/src/SFMT-sse2.h

@@ -0,0 +1,121 @@
+#pragma once
+/**
+ * @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
+
+inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
+				__m128i c, __m128i d);
+
+/**
+ * This function represents the recursion formula.
+ * @param r an output
+ * @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
+ */
+inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
+				__m128i c, __m128i d)
+{
+    __m128i v, x, y, z;
+
+    y = _mm_srli_epi32(b, SFMT_SR1);
+    z = _mm_srli_si128(c, SFMT_SR2);
+    v = _mm_slli_epi32(d, SFMT_SL1);
+    z = _mm_xor_si128(z, a);
+    z = _mm_xor_si128(z, v);
+    x = _mm_slli_si128(a, SFMT_SL2);
+    y = _mm_and_si128(y, sse2_param_mask.si);
+    z = _mm_xor_si128(z, x);
+    z = _mm_xor_si128(z, y);
+    *r = z;
+}
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+	mm_recursion(&pstate[i].si, pstate[i].si,
+		     pstate[i + SFMT_POS1].si, r1, r2);
+	r1 = r2;
+	r2 = pstate[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+	mm_recursion(&pstate[i].si, pstate[i].si,
+		     pstate[i + SFMT_POS1 - SFMT_N].si,
+		     r1, r2);
+	r1 = r2;
+	r2 = pstate[i].si;
+    }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
+{
+    int i, j;
+    __m128i r1, r2;
+    w128_t * pstate = sfmt->state;
+
+    r1 = pstate[SFMT_N - 2].si;
+    r2 = pstate[SFMT_N - 1].si;
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+	mm_recursion(&array[i].si, pstate[i].si,
+		     pstate[i + SFMT_POS1].si, r1, r2);
+	r1 = r2;
+	r2 = array[i].si;
+    }
+    for (; i < SFMT_N; i++) {
+	mm_recursion(&array[i].si, pstate[i].si,
+		     array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+	r1 = r2;
+	r2 = array[i].si;
+    }
+    for (; i < size - SFMT_N; i++) {
+	mm_recursion(&array[i].si, array[i - SFMT_N].si,
+		     array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+	r1 = r2;
+	r2 = array[i].si;
+    }
+    for (j = 0; j < 2 * SFMT_N - size; j++) {
+	pstate[j] = array[j + size - SFMT_N];
+    }
+    for (; i < size; i++, j++) {
+	mm_recursion(&array[i].si, array[i - SFMT_N].si,
+		     array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+	r1 = r2;
+	r2 = array[i].si;
+	pstate[j] = array[i];
+    }
+}
+
+
+#endif

sfmt.mod/src/SFMT.c

@@ -0,0 +1,437 @@
+/**
+ * @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.
+ * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+ * University and The University of Tokyo.
+ * Copyright (C) 2013 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ * All rights reserved.
+ *
+ * The 3-clause BSD License is applied to this software, see
+ * LICENSE.txt
+ */
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+#include <string.h>
+#include <assert.h>
+#include "SFMT.h"
+#include "SFMT-params.h"
+#include "SFMT-common.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
+
+/*----------------
+  STATIC FUNCTIONS
+  ----------------*/
+inline static int idxof(int i);
+inline static void gen_rand_array(sfmt_t * sfmt, 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(sfmt_t * sfmt);
+#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)
+/**
+ * parameters used by sse2.
+ */
+  static const w128_t sse2_param_mask = {{SFMT_MSK1, SFMT_MSK2,
+                                          SFMT_MSK3, SFMT_MSK4}};
+  #if defined(_MSC_VER)
+    #include "SFMT-sse2-msc.h"
+  #else
+    #include "SFMT-sse2.h"
+  #endif
+#elif defined(HAVE_NEON)
+  #include "SFMT-neon.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
+
+#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2)) && (!defined(HAVE_NEON))
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ *
+ * @param sfmt SFMT internal state
+ * @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(sfmt_t * sfmt, w128_t *array, int size) {
+    int i, j;
+    w128_t *r1, *r2;
+
+    r1 = &sfmt->state[SFMT_N - 2];
+    r2 = &sfmt->state[SFMT_N - 1];
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        do_recursion(&array[i], &sfmt->state[i], &sfmt->state[i + SFMT_POS1], r1, r2);
+        r1 = r2;
+        r2 = &array[i];
+    }
+    for (; i < SFMT_N; i++) {
+        do_recursion(&array[i], &sfmt->state[i],
+                     &array[i + SFMT_POS1 - SFMT_N], r1, r2);
+        r1 = r2;
+        r2 = &array[i];
+    }
+    for (; i < size - SFMT_N; i++) {
+        do_recursion(&array[i], &array[i - SFMT_N],
+                     &array[i + SFMT_POS1 - SFMT_N], r1, r2);
+        r1 = r2;
+        r2 = &array[i];
+    }
+    for (j = 0; j < 2 * SFMT_N - size; j++) {
+        sfmt->state[j] = array[j + size - SFMT_N];
+    }
+    for (; i < size; i++, j++) {
+        do_recursion(&array[i], &array[i - SFMT_N],
+                     &array[i + SFMT_POS1 - SFMT_N], r1, r2);
+        r1 = r2;
+        r2 = &array[i];
+        sfmt->state[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}
+ * @param sfmt SFMT internal state
+ */
+static void period_certification(sfmt_t * sfmt) {
+    uint32_t inner = 0;
+    int i, j;
+    uint32_t work;
+    uint32_t *psfmt32 = &sfmt->state[0].u[0];
+    const uint32_t parity[4] = {SFMT_PARITY1, SFMT_PARITY2,
+                                SFMT_PARITY3, SFMT_PARITY4};
+
+    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
+  ----------------*/
+#define UNUSED_VARIABLE(x) (void)(x)
+/**
+ * This function returns the identification string.
+ * The string shows the word size, the Mersenne exponent,
+ * and all parameters of this generator.
+ * @param sfmt SFMT internal state
+ */
+const char *sfmt_get_idstring(sfmt_t * sfmt) {
+    UNUSED_VARIABLE(sfmt);
+    return SFMT_IDSTR;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array32() function.
+ * @param sfmt SFMT internal state
+ * @return minimum size of array used for fill_array32() function.
+ */
+int sfmt_get_min_array_size32(sfmt_t * sfmt) {
+    UNUSED_VARIABLE(sfmt);
+    return SFMT_N32;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array64() function.
+ * @param sfmt SFMT internal state
+ * @return minimum size of array used for fill_array64() function.
+ */
+int sfmt_get_min_array_size64(sfmt_t * sfmt) {
+    UNUSED_VARIABLE(sfmt);
+    return SFMT_N64;
+}
+
+#if !defined(HAVE_SSE2) && !defined(HAVE_ALTIVEC) && !defined(HAVE_NEON)
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+    int i;
+    w128_t *r1, *r2;
+
+    r1 = &sfmt->state[SFMT_N - 2];
+    r2 = &sfmt->state[SFMT_N - 1];
+    for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+        do_recursion(&sfmt->state[i], &sfmt->state[i],
+                     &sfmt->state[i + SFMT_POS1], r1, r2);
+        r1 = r2;
+        r2 = &sfmt->state[i];
+    }
+    for (; i < SFMT_N; i++) {
+        do_recursion(&sfmt->state[i], &sfmt->state[i],
+                     &sfmt->state[i + SFMT_POS1 - SFMT_N], r1, r2);
+        r1 = r2;
+        r2 = &sfmt->state[i];
+    }
+}
+#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 sfmt SFMT internal state
+ * @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 sfmt_fill_array32(sfmt_t * sfmt, uint32_t *array, int size) {
+    assert(sfmt->idx == SFMT_N32);
+    assert(size % 4 == 0);
+    assert(size >= SFMT_N32);
+
+    gen_rand_array(sfmt, (w128_t *)array, size / 4);
+    sfmt->idx = SFMT_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.
+ *
+ * @param sfmt SFMT internal state
+ * 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 sfmt_fill_array64(sfmt_t * sfmt, uint64_t *array, int size) {
+    assert(sfmt->idx == SFMT_N32);
+    assert(size % 2 == 0);
+    assert(size >= SFMT_N64);
+
+    gen_rand_array(sfmt, (w128_t *)array, size / 2);
+    sfmt->idx = SFMT_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 sfmt SFMT internal state
+ * @param seed a 32-bit integer used as the seed.
+ */
+void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed) {
+    int i;
+
+    uint32_t *psfmt32 = &sfmt->state[0].u[0];
+
+    psfmt32[idxof(0)] = seed;
+    for (i = 1; i < SFMT_N32; i++) {
+        psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)]
+                                            ^ (psfmt32[idxof(i - 1)] >> 30))
+            + i;
+    }
+    sfmt->idx = SFMT_N32;
+    period_certification(sfmt);
+}
+
+/**
+ * This function initializes the internal state array,
+ * with an array of 32-bit integers used as the seeds
+ * @param sfmt SFMT internal state
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ */
+void sfmt_init_by_array(sfmt_t * sfmt, uint32_t *init_key, int key_length) {
+    int i, j, count;
+    uint32_t r;
+    int lag;
+    int mid;
+    int size = SFMT_N * 4;
+    uint32_t *psfmt32 = &sfmt->state[0].u[0];
+
+    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_t));
+    if (key_length + 1 > SFMT_N32) {
+        count = key_length + 1;
+    } else {
+        count = SFMT_N32;
+    }
+    r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)]
+              ^ psfmt32[idxof(SFMT_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) % SFMT_N32)]
+                  ^ psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
+        psfmt32[idxof((i + mid) % SFMT_N32)] += r;
+        r += init_key[j] + i;
+        psfmt32[idxof((i + mid + lag) % SFMT_N32)] += r;
+        psfmt32[idxof(i)] = r;
+        i = (i + 1) % SFMT_N32;
+    }
+    for (; j < count; j++) {
+        r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % SFMT_N32)]
+                  ^ psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
+        psfmt32[idxof((i + mid) % SFMT_N32)] += r;
+        r += i;
+        psfmt32[idxof((i + mid + lag) % SFMT_N32)] += r;
+        psfmt32[idxof(i)] = r;
+        i = (i + 1) % SFMT_N32;
+    }
+    for (j = 0; j < SFMT_N32; j++) {
+        r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % SFMT_N32)]
+                  + psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
+        psfmt32[idxof((i + mid) % SFMT_N32)] ^= r;
+        r -= i;
+        psfmt32[idxof((i + mid + lag) % SFMT_N32)] ^= r;
+        psfmt32[idxof(i)] = r;
+        i = (i + 1) % SFMT_N32;
+    }
+
+    sfmt->idx = SFMT_N32;
+    period_certification(sfmt);
+}
+#if defined(__cplusplus)
+}
+#endif

sfmt.mod/src/SFMT.h

@@ -0,0 +1,304 @@
+#pragma once
+/**
+ * @file SFMT.h
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
+ * number generator using C structure.
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (The University of Tokyo)
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+ * University and The University of Tokyo.
+ * All rights reserved.
+ *
+ * The 3-clause 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 SFMTST_H
+#define SFMTST_H
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+#include <stdio.h>
+#include <assert.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
+
+#include "SFMT-params.h"
+
+/*------------------------------------------
+  128-bit SIMD like data type for 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];
+    uint64_t u64[2];
+};
+#elif defined(HAVE_NEON)
+  #include <arm_neon.h>
+
+/** 128-bit data structure */
+union W128_T {
+    uint32_t u[4];
+    uint64_t u64[2];
+    uint32x4_t si;
+};
+#elif defined(HAVE_SSE2)
+  #include <emmintrin.h>
+
+/** 128-bit data structure */
+union W128_T {
+    uint32_t u[4];
+    uint64_t u64[2];
+    __m128i si;
+};
+#else
+/** 128-bit data structure */
+union W128_T {
+    uint32_t u[4];
+    uint64_t u64[2];
+};
+#endif
+
+/** 128-bit data type */
+typedef union W128_T w128_t;
+
+/**
+ * SFMT internal state
+ */
+struct SFMT_T {
+    /** the 128-bit internal state array */
+    w128_t state[SFMT_N];
+    /** index counter to the 32-bit internal state array */
+    int idx;
+};
+
+typedef struct SFMT_T sfmt_t;
+
+void sfmt_fill_array32(sfmt_t * sfmt, uint32_t * array, int size);
+void sfmt_fill_array64(sfmt_t * sfmt, uint64_t * array, int size);
+void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed);
+void sfmt_init_by_array(sfmt_t * sfmt, uint32_t * init_key, int key_length);
+const char * sfmt_get_idstring(sfmt_t * sfmt);
+int sfmt_get_min_array_size32(sfmt_t * sfmt);
+int sfmt_get_min_array_size64(sfmt_t * sfmt);
+void sfmt_gen_rand_all(sfmt_t * sfmt);
+
+#ifndef ONLY64
+/**
+ * This function generates and returns 32-bit pseudorandom number.
+ * init_gen_rand or init_by_array must be called before this function.
+ * @param sfmt SFMT internal state
+ * @return 32-bit pseudorandom number
+ */
+inline static uint32_t sfmt_genrand_uint32(sfmt_t * sfmt) {
+    uint32_t r;
+    uint32_t * psfmt32 = &sfmt->state[0].u[0];
+
+    if (sfmt->idx >= SFMT_N32) {
+        sfmt_gen_rand_all(sfmt);
+        sfmt->idx = 0;
+    }
+    r = psfmt32[sfmt->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.
+ * @param sfmt SFMT internal state
+ * @return 64-bit pseudorandom number
+ */
+inline static uint64_t sfmt_genrand_uint64(sfmt_t * sfmt) {
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+    uint32_t * psfmt32 = &sfmt->state[0].u[0];
+    uint32_t r1, r2;
+#else
+    uint64_t r;
+#endif
+    uint64_t * psfmt64 = &sfmt->state[0].u64[0];
+    assert(sfmt->idx % 2 == 0);
+
+    if (sfmt->idx >= SFMT_N32) {
+        sfmt_gen_rand_all(sfmt);
+        sfmt->idx = 0;
+    }
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+    r1 = psfmt32[sfmt->idx];
+    r2 = psfmt32[sfmt->idx + 1];
+    sfmt->idx += 2;
+    return ((uint64_t)r2 << 32) | r1;
+#else
+    r = psfmt64[sfmt->idx / 2];
+    sfmt->idx += 2;
+    return r;
+#endif
+}
+
+/* =================================================
+   The following real versions are due to Isaku Wada
+   ================================================= */
+/**
+ * converts an unsigned 32-bit number to a double on [0,1]-real-interval.
+ * @param v 32-bit unsigned integer
+ * @return double on [0,1]-real-interval
+ */
+inline static double sfmt_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
+ * @param sfmt SFMT internal state
+ * @return double on [0,1]-real-interval
+ */
+inline static double sfmt_genrand_real1(sfmt_t * sfmt)
+{
+    return sfmt_to_real1(sfmt_genrand_uint32(sfmt));
+}
+
+/**
+ * converts an unsigned 32-bit integer to a double on [0,1)-real-interval.
+ * @param v 32-bit unsigned integer
+ * @return double on [0,1)-real-interval
+ */
+inline static double sfmt_to_real2(uint32_t v)
+{
+    return v * (1.0/4294967296.0);
+    /* divided by 2^32 */
+}
+
+/**
+ * generates a random number on [0,1)-real-interval
+ * @param sfmt SFMT internal state
+ * @return double on [0,1)-real-interval
+ */
+inline static double sfmt_genrand_real2(sfmt_t * sfmt)
+{
+    return sfmt_to_real2(sfmt_genrand_uint32(sfmt));
+}
+
+/**
+ * converts an unsigned 32-bit integer to a double on (0,1)-real-interval.
+ * @param v 32-bit unsigned integer
+ * @return double on (0,1)-real-interval
+ */
+inline static double sfmt_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
+ * @param sfmt SFMT internal state
+ * @return double on (0,1)-real-interval
+ */
+inline static double sfmt_genrand_real3(sfmt_t * sfmt)
+{
+    return sfmt_to_real3(sfmt_genrand_uint32(sfmt));
+}
+
+/**
+ * converts an unsigned 32-bit integer to double on [0,1)
+ * with 53-bit resolution.
+ * @param v 32-bit unsigned integer
+ * @return double on [0,1)-real-interval with 53-bit resolution.
+ */
+inline static double sfmt_to_res53(uint64_t v)
+{
+    return (v >> 11) * (1.0/9007199254740992.0);
+}
+
+/**
+ * generates a random number on [0,1) with 53-bit resolution
+ * @param sfmt SFMT internal state
+ * @return double on [0,1) with 53-bit resolution
+ */
+inline static double sfmt_genrand_res53(sfmt_t * sfmt)
+{
+    return sfmt_to_res53(sfmt_genrand_uint64(sfmt));
+}
+
+
+/* =================================================
+   The following function are added by Saito.
+   ================================================= */
+/**
+ * generates a random number on [0,1) with 53-bit resolution from two
+ * 32 bit integers
+ */
+inline static double sfmt_to_res53_mix(uint32_t x, uint32_t y)
+{
+    return sfmt_to_res53(x | ((uint64_t)y << 32));
+}
+
+/**
+ * generates a random number on [0,1) with 53-bit resolution
+ * using two 32bit integers.
+ * @param sfmt SFMT internal state
+ * @return double on [0,1) with 53-bit resolution
+ */
+inline static double sfmt_genrand_res53_mix(sfmt_t * sfmt)
+{
+    uint32_t x, y;
+
+    x = sfmt_genrand_uint32(sfmt);
+    y = sfmt_genrand_uint32(sfmt);
+    return sfmt_to_res53_mix(x, y);
+}
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif