stream.generator.SFMT19937 Maven / Gradle / Ivy
The newest version!
/*
* streams library
*
* Copyright (C) 2011-2014 by Christian Bockermann, Hendrik Blom
*
* streams is a library, API and runtime environment for processing high
* volume data streams. It is composed of three submodules "stream-api",
* "stream-core" and "stream-runtime".
*
* The streams library (and its submodules) is free software: you can
* redistribute it and/or modify it under the terms of the
* GNU Affero General Public License as published by the Free Software
* Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* The stream.ai library (and its submodules) is distributed in the hope
* that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see http://www.gnu.org/licenses/.
*/
package stream.generator;
/**
* An adapation of SFMT
* (SIMD oriented Fast Mersenne Twister) version 1.3 by Mutsuo Saito
* (Hiroshima University) and Makoto Matsumoto (Hiroshima University). This
* adapation supports only the period 219937 − 1; the original
* supports some longer and shorter periods.
*
*
* The license (a modified BSD License) for the original C code from which this
* code is adapted:
*
*
* 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.
*
*
* @author Adrian King (ceroxylon at hotmail.com)
*/
public class SFMT19937 {
// Static variables
// ////////////////////////////////////////////////////////////
/**
* Mersenne Exponent. The period of the sequence is a multiple of 2
* MEXP − 1. If you adapt this code to support a
* different exponent, you must change many of the other constants here as
* well; consult the original C code.
*/
public final static int MEXP = 19937;
/**
* The SFMT generator has an internal state array of 128-bit integers, and
* N is its size.
*/
final static int N = MEXP / 128 + 1;
/**
* N32 is the size of internal state array when regarded as an
* array of 32-bit integers.
*/
public final static int N32 = N * 4;
/**
* The pick up position of the array.
*/
final static int POS1 = 122;
/**
* The parameter of shift left as four 32-bit registers.
*/
final static int SL1 = 18;
/**
* The parameter of shift left as one 128-bit register. The 128-bit integer
* is shifted by SL2 * 8 bits.
*/
final static int SL2 = 1;
/**
* The parameter of shift right as four 32-bit registers.
*/
final static int SR1 = 11;
/**
* The parameter of shift right as one 128-bit register. The 128-bit integer
* is shifted by SL2 * 8 bits.
*/
final static int SR2 = 1;
/**
* A bitmask parameter, used in the recursion to break symmetry of SIMD.
*/
final static int MSK1 = 0xdfffffef;
/**
* A bitmask parameter, used in the recursion to break symmetry of SIMD.
*/
final static int MSK2 = 0xddfecb7f;
/**
* A bitmask parameter, used in the recursion to break symmetry of SIMD.
*/
final static int MSK3 = 0xbffaffff;
/**
* A bitmask parameter, used in the recursion to break symmetry of SIMD.
*/
final static int MSK4 = 0xbffffff6;
/**
* Part of a 128-bit period certification vector.
*/
final static int PARITY1 = 0x00000001;
/**
* Part of a 128-bit period certification vector.
*/
final static int PARITY2 = 0x00000000;
/**
* Part of a 128-bit period certification vector.
*/
final static int PARITY3 = 0x00000000;
/**
* Part of a 128-bit period certification vector.
*/
final static int PARITY4 = 0x13c9e684;
/**
* A parity check vector which certifies the period of 2{@link #MEXP}
* .
*/
final static int parity[] = { PARITY1, PARITY2, PARITY3, PARITY4 };
/**
* A number mixed with the time of day to provide a unique seed to each
* generator of this type allocated.
*/
static long uniquifier = 314159265358979L;
// Instance variables
// //////////////////////////////////////////////////////////
/**
* The internal state array. Blocks of four consecutive ints
* are often treated as a single 128-bit integer that is
* little-endian—that is, its low-order bits are at lower indices in
* the array, and high-order bits at higher indices.
*/
final int[] sfmt = new int[N32];
/**
* Index counter of the next int to return from {@link #next}.
*/
int idx = N32;
// Constructors
// ////////////////////////////////////////////////////////////////
/**
* Constructor that initially uses a seed based on the time of day.
*/
public SFMT19937() {
this(System.nanoTime() + uniquifier++);
}
/**
* Constructor that initializes the internal state array by calling
* {@link #setSeed} with the specified argument.
*
* @param seed
* initial seed for this generator.
*/
public SFMT19937(long seed) {
setSeed(seed);
}
// Instance methods
// ////////////////////////////////////////////////////////////
/**
* Applies the recursion formula.
*
* @param r
* output array.
* @param rI
* index in r.
* @param a
* state array.
* @param aI
* index in a.
* @param b
* state array.
* @param bI
* index in b.
* @param c
* state array.
* @param cI
* index in c.
* @param d
* state array.
* @param dI
* index in d.
*/
void doRecursion(int[] r, int rI, int[] a, int aI, int[] b, int bI,
int[] c, int cI, int[] d, int dI) {
// 128-bit shift: x = a << SL2 * 8:
final int lShift = SL2 * 8;
int a0 = a[aI], a1 = a[aI + 1], a2 = a[aI + 2], a3 = a[aI + 3];
// for SL2 <= 3, this is more concise, but possibly not as fast (haven't
// timed it):
// int x0 = a0 << lShift,
// x1 = (a1 << lShift) | (a0 >>> (32 - lShift)),
// x2 = (a2 << lShift) | (a1 >>> (32 - lShift)),
// x3 = (a3 << lShift) | (a2 >>> (32 - lShift));
long hi = ((long) a3 << 32) | (a2 & (-1L >>> 32)), lo = ((long) a1 << 32)
| (a0 & (-1L >>> 32)), outLo = lo << lShift, outHi = (hi << lShift)
| (lo >>> (64 - lShift));
int x0 = (int) outLo, x1 = (int) (outLo >>> 32), x2 = (int) outHi, x3 = (int) (outHi >>> 32);
// 128-bit shift: y = c >>> SR2 * 8:
final int rShift = SR2 * 8;
hi = ((long) c[cI + 3] << 32) | (c[cI + 2] & (-1L >>> 32));
lo = ((long) c[cI + 1] << 32) | (c[cI] & (-1L >>> 32));
outHi = hi >>> rShift;
outLo = (lo >>> rShift) | (hi << (64 - rShift));
int y0 = (int) outLo, y1 = (int) (outLo >>> 32), y2 = (int) outHi, y3 = (int) (outHi >>> 32);
// rest of forumula:
r[rI] = a0 ^ x0 ^ ((b[bI] >>> SR1) & MSK1) ^ y0 ^ (d[dI] << SL1);
r[rI + 1] = a1 ^ x1 ^ ((b[bI + 1] >>> SR1) & MSK2) ^ y1
^ (d[dI + 1] << SL1);
r[rI + 2] = a2 ^ x2 ^ ((b[bI + 2] >>> SR1) & MSK3) ^ y2
^ (d[dI + 2] << SL1);
r[rI + 3] = a3 ^ x3 ^ ((b[bI + 3] >>> SR1) & MSK4) ^ y3
^ (d[dI + 3] << SL1);
}
/**
* Fills the internal state array with pseudorandom integers.
*/
void genRandAll() {
int i = 0, r1 = 4 * (N - 2), r2 = 4 * (N - 1);
for (; i < 4 * (N - POS1); i += 4) {
doRecursion(sfmt, i, sfmt, i, sfmt, i + 4 * POS1, sfmt, r1, sfmt,
r2);
r1 = r2;
r2 = i;
}
for (; i < 4 * N; i += 4) {
doRecursion(sfmt, i, sfmt, i, sfmt, i + 4 * (POS1 - N), sfmt, r1,
sfmt, r2);
r1 = r2;
r2 = i;
}
}
/**
* Fills a user-specified array with pseudorandom integers.
*
* @param array
* 128-bit array to be filled with pseudorandom numbers.
* @param size
* number of elements of array to fill.
* @throws IllegalArgumentException
* if size is greater than the length of
* array, or if size is less than
* {@link #N32}, or is not a multiple of 4.
*/
void genRandArray(int[] array, int size) {
if (size > array.length)
throw new IllegalArgumentException("Given size " + size
+ " exceeds array length " + array.length);
if (size < N32)
throw new IllegalArgumentException("Size must be at least " + N32
+ ", but is " + size);
if (size % 4 != 0)
throw new IllegalArgumentException("Size must be a multiple of 4: "
+ size);
int i = 0, j = 0, r1I = 4 * (N - 2), r2I = 4 * (N - 1);
int[] r1 = sfmt, r2 = sfmt;
for (; i < 4 * (N - POS1); i += 4) {
doRecursion(array, i, sfmt, i, sfmt, i + 4 * POS1, r1, r1I, r2, r2I);
r1 = r2;
r1I = r2I;
r2 = array;
r2I = i;
}
for (; i < 4 * N; i += 4) {
doRecursion(array, i, sfmt, i, array, i + 4 * (POS1 - N), r1, r1I,
r2, r2I);
assert r1 == r2;
r1I = r2I;
assert r2 == array;
r2I = i;
}
for (; i < size - 4 * N; i += 4) {
doRecursion(array, i, array, i - 4 * N, array, i + 4 * (POS1 - N),
r1, r1I, r2, r2I);
assert r1 == r2;
r1I = r2I;
assert r2 == array;
r2I = i;
}
for (; j < 4 * 2 * N - size; j++)
sfmt[j] = array[j + size - 4 * N];
for (; i < size; i += 4, j += 4) {
doRecursion(array, i, array, i - 4 * N, array, i + 4 * (POS1 - N),
r1, r1I, r2, r2I);
assert r1 == r2;
r1I = r2I;
assert r2 == array;
r2I = i;
sfmt[j] = array[i];
sfmt[j + 1] = array[i + 1];
sfmt[j + 2] = array[i + 2];
sfmt[j + 3] = array[i + 3];
}
}
/**
* Used by {@link #initByArray}.
*
* @param x
* 32-bit integer.
* @return 32-bit integer.
*/
int func1(int x) {
return (x ^ (x >>> 27)) * 1664525;
}
/**
* Used by {@link #initByArray}.
*
* @param x
* 32-bit integer.
* @return 32-bit integer.
*/
int func2(int x) {
return (x ^ (x >>> 27)) * 1566083941;
}
/**
* Certifies the period of 2{@link #MEXP}.
*/
void periodCertification() {
int inner = 0;
for (int i = 0; i < 4; i++)
inner ^= sfmt[i] & parity[i];
for (int i = 16; i > 0; i >>= 1)
inner ^= inner >> i;
if ((inner & 1) != 0) // check OK
return;
for (int i = 0; i < 4; i++) {
int work = 1;
for (int j = 0; j < 32; j++) {
if ((work & parity[i]) != 0) {
sfmt[i] ^= work;
return;
}
work <<= 1;
}
}
}
/**
* Returns the smallest array size allowed as an argument to
* {@link #fillArray}.
*/
public int minArraySize() {
return N32;
}
/**
* Generates and returns the next 32-bit pseudorandom number.
*
* @return next number.
*/
public int next() {
if (idx >= N32) {
genRandAll();
idx = 0;
}
return sfmt[idx++];
}
/**
* Fills the given array with pseudorandom 32-bit integers. Equivalent to
* {@link #fillArray(int[],int)} applied to
* (array,array.length).
*
* @param array
* array to fill.
*/
public void fillArray(int[] array) {
genRandArray(array, array.length);
}
/**
* Fills the given array with the specified number of pseudorandom 32-bit
* integers. The specified number of elements must be a multiple of four.
*
* @param array
* array to fill.
* @param elems
* the number of elements of array (starting at
* index zero) to fill; subsequent elements are not modified.
* @throws IllegalArgumentException
* if elems is greater than the length of
* array, or is less than {@link #N32}, or is not a
* multiple of 4.
*/
public void fillArray(int[] array, int elems) {
genRandArray(array, elems);
idx = N32;
}
/**
* Initializes the internal state array with a 32-bit seed.
*
* @param seed
* 32-bit seed.
*/
public void setIntSeed(int seed) {
sfmt[0] = seed;
for (int i = 1; i < N32; i++) {
int prev = sfmt[i - 1];
sfmt[i] = 1812433253 * (prev ^ (prev >>> 30)) + i;
}
periodCertification();
idx = N32;
}
/**
* Initializes the internal state array with a 64-bit seed.
*
* @param seed
* 64-bit seed.
*/
public void setSeed(long seed) {
initByArray((int) seed, (int) (seed >>> 32));
}
/**
* Initializes the internal state array with an array of 32-bit integers.
*
* @param key
* array of 32-bit integers, used as a seed.
*/
public void initByArray(int... key) {
@SuppressWarnings("unused")
int lag = N32 >= 623 ? 11 : N32 >= 68 ? 7 : N32 >= 39 ? 5 : 3, mid = (N32 - lag) / 2;
for (int i = sfmt.length - 1; i >= 0; i--)
sfmt[i] = 0x8b8b8b8b;
int count = key.length >= N32 ? key.length : N32 - 1, r = func1(0x8b8b8b8b);
sfmt[mid] += r;
r += key.length;
sfmt[mid + lag] += r;
sfmt[0] = r;
int i = 1, j = 0;
for (; j < count && j < key.length; j++) {
r = func1(sfmt[i] ^ sfmt[(i + mid) % N32]
^ sfmt[(i + N32 - 1) % N32]);
sfmt[(i + mid) % N32] += r;
r += key[j] + i;
sfmt[(i + mid + lag) % N32] += r;
sfmt[i] = r;
i = (i + 1) % N32;
}
for (; j < count; j++) {
r = func1(sfmt[i] ^ sfmt[(i + mid) % N32]
^ sfmt[(i + N32 - 1) % N32]);
sfmt[(i + mid) % N32] += r;
r += i;
sfmt[(i + mid + lag) % N32] += r;
sfmt[i] = r;
i = (i + 1) % N32;
}
for (j = 0; j < N32; j++) {
r = func2(sfmt[i] + sfmt[(i + mid) % N32]
+ sfmt[(i + N32 - 1) % N32]);
sfmt[(i + mid) % N32] ^= r;
r -= i;
sfmt[(i + mid + lag) % N32] ^= r;
sfmt[i] = r;
i = (i + 1) % N32;
}
periodCertification();
idx = N32;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy