diff --git a/src/internal/mersenne.ts b/src/internal/mersenne.ts index 82a0915baf9..392336e6f38 100644 --- a/src/internal/mersenne.ts +++ b/src/internal/mersenne.ts @@ -1,325 +1,190 @@ +// Based on: https://github.com/dubzzz/pure-rand/blob/5da623a4eebf47f01c5b87c43c9fb4d43ad949bd/src/generator/MersenneTwister.ts + +const N = 624; +const M = 397; +const R = 31; +const A = 0x9908b0df; +const F = 1812433253; +const U = 11; +const S = 7; +const B = 0x9d2c5680; +const T = 15; +const C = 0xefc60000; +const L = 18; +const MASK_LOWER = 2 ** R - 1; +const MASK_UPPER = 2 ** R; +const HIGH_MULTIPLIER_53 = 67108864.0; +const FLOATIFY_32 = 1.0 / 4294967296.0; +const FLOATIFY_53 = 1.0 / 9007199254740992.0; + +const { imul, trunc } = Math; + /** - * Copyright (c) 2022-2023 Faker - * - * This is a version of the original source code migrated to TypeScript and - * modified by the Faker team. - * - * Check LICENSE for more details on copyright. - * - * ----------------------------------------------------------------------------- - * - * Copyright (c) 2006 Y. Okada - * - * This program is a JavaScript version of Mersenne Twister, with concealment - * and encapsulation in class, an almost straight conversion from the original - * program, mt19937ar.c, translated by Y. Okada on July 17, 2006, and modified - * a little at July 20, 2006, but there are not any substantial differences. - * - * In this program, procedure descriptions and comments of original source code - * were not removed. - * - * Lines commented with //c// were originally descriptions of c procedure. - * And a few following lines are appropriate JavaScript descriptions. - * - * Lines commented with /* and *\/ are original comments. - * Lines commented with // are additional comments in this JavaScript version. - * - * Before using this version, create at least one instance of - * MersenneTwister19937 class, and initialize the each state, given below - * in C comments, of all the instances. - * - * ----------------------------------------------------------------------------- + * Generates the random state from the given number or array. * - * A C-program for MT19937, with initialization improved 2002/1/26. - * Coded by Takuji Nishimura and Makoto Matsumoto. - * - * Before using, initialize the state by using init_genrand(seed) - * or init_by_array(init_key, key_length). - * - * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions are met: - * - * 1. Redistributions of source code must retain the above copyright notice, - * this list of conditions and the following disclaimer. - * - * 2. Redistributions in binary form must reproduce the above copyright notice, - * this list of conditions and the following disclaimer in the documentation - * and/or other materials provided with the distribution. - * - * 3. The names of its contributors may not be used to endorse or promote - * products derived from this software without specific prior written - * permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" - * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, - * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS - * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR - * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT - * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS - * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN - * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) - * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF - * THE POSSIBILITY OF SUCH DAMAGE. - * - * Any feedback is very welcome. - * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html - * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) - * - * @internal + * @param seed The seed to generate the random state from. */ -export class MersenneTwister19937 { - private readonly N = 624; - private readonly M = 397; - private readonly MATRIX_A = 0x9908b0df; // constant vector a - private readonly UPPER_MASK = 0x80000000; // most significant w-r bits - private readonly LOWER_MASK = 0x7fffffff; // least significant r bits - private mt: number[] = Array.from({ length: this.N }); // the array for the state vector - private mti = this.N + 1; // mti==N+1 means mt[N] is not initialized +function seedFrom(seed: number | number[]): number[] { + return typeof seed === 'number' ? numberSeeded(seed) : arraySeeded(seed); +} - /** - * Returns a 32-bits unsigned integer from an operand to which applied a bit - * operator. - * - * @param n1 number to unsign - */ - private unsigned32(n1: number): number { - return n1 < 0 ? (n1 ^ this.UPPER_MASK) + this.UPPER_MASK : n1; - } +/** + * Generates the random state from the given number. + * + * @param seed The seed to generate the random state from. + */ +function numberSeeded(seed: number): number[] { + const out = Array.from({ length: N }); + out[0] = seed; - /** - * Emulates lowerflow of a c 32-bits unsigned integer variable, instead of - * the operator -. These both arguments must be non-negative integers - * expressible using unsigned 32 bits. - * - * @param n1 dividend - * @param n2 divisor - */ - private subtraction32(n1: number, n2: number): number { - return n1 < n2 - ? this.unsigned32((0x100000000 - (n2 - n1)) & 0xffffffff) - : n1 - n2; + for (let idx = 1; idx !== N; ++idx) { + const xored = out[idx - 1] ^ (out[idx - 1] >>> 30); + out[idx] = trunc(imul(F, xored) + idx); } - /** - * Emulates overflow of a c 32-bits unsigned integer variable, instead of the operator +. - * these both arguments must be non-negative integers expressible using unsigned 32 bits. - * - * @param n1 number one for addition - * @param n2 number two for addition - */ - private addition32(n1: number, n2: number): number { - return this.unsigned32((n1 + n2) & 0xffffffff); - } + return out; +} - /** - * Emulates overflow of a c 32-bits unsigned integer variable, instead of the operator *. - * These both arguments must be non-negative integers expressible using unsigned 32 bits. - * - * @param n1 number one for multiplication - * @param n2 number two for multiplication - */ - private multiplication32(n1: number, n2: number): number { - let sum = 0; - for (let i = 0; i < 32; ++i) { - if ((n1 >>> i) & 0x1) { - sum = this.addition32(sum, this.unsigned32(n2 << i)); - } +/** + * Generates the random state from the given array. + * + * @param seed The seed to generate the random state from. + */ +function arraySeeded(seed: number[]): number[] { + const out = numberSeeded(19650218); + + let idxOut = 1; + let idxSeed = 0; + + for (let iteration = Math.max(N, seed.length); iteration !== 0; --iteration) { + const xored = out[idxOut - 1] ^ (out[idxOut - 1] >>> 30); + out[idxOut] = trunc( + (out[idxOut] ^ imul(xored, 1664525)) + seed[idxSeed] + idxSeed + ); + idxOut++; + idxSeed++; + + if (idxOut >= N) { + out[0] = out[N - 1]; + idxOut = 1; } - return sum; + if (idxSeed >= seed.length) { + idxSeed = 0; + } } - /** - * Initializes mt[N] with a seed. - * - * @param seed the seed to use - */ - initGenrand(seed: number): void { - this.mt[0] = this.unsigned32(seed & 0xffffffff); - for (this.mti = 1; this.mti < this.N; this.mti++) { - this.mt[this.mti] = - //(1812433253 * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); - this.addition32( - this.multiplication32( - 1812433253, - this.unsigned32( - this.mt[this.mti - 1] ^ (this.mt[this.mti - 1] >>> 30) - ) - ), - this.mti - ); - - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array mt[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ - this.mt[this.mti] = this.unsigned32(this.mt[this.mti] & 0xffffffff); + for (let iteration = N - 1; iteration !== 0; iteration--) { + out[idxOut] = trunc( + (out[idxOut] ^ + imul(out[idxOut - 1] ^ (out[idxOut - 1] >>> 30), 1566083941)) - + idxOut + ); + idxOut++; + + if (idxOut >= N) { + out[0] = out[N - 1]; + idxOut = 1; } } - /** - * Initialize by an array with array-length. - * - * @param initKey is the array for initializing keys - * @param keyLength is its length - */ - initByArray(initKey: number[], keyLength: number): void { - this.initGenrand(19650218); - let i = 1; - let j = 0; - let k = Math.max(this.N, keyLength); - for (; k; k--) { - // mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525)) + init_key[j] + j; - this.mt[i] = this.addition32( - this.addition32( - this.unsigned32( - this.mt[i] ^ - this.multiplication32( - this.unsigned32(this.mt[i - 1] ^ (this.mt[i - 1] >>> 30)), - 1664525 - ) - ), - initKey[j] - ), - j - ); - // mt[i] &= 0xffffffff; for WORDSIZE > 32 machines - this.mt[i] = this.unsigned32(this.mt[i] & 0xffffffff); - i++; - j++; - if (i >= this.N) { - this.mt[0] = this.mt[this.N - 1]; - i = 1; - } - - if (j >= keyLength) { - j = 0; - } - } + out[0] = 0x80000000; + return out; +} - for (k = this.N - 1; k; k--) { - // mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941)) - i - this.mt[i] = this.subtraction32( - this.unsigned32( - this.mt[i] ^ - this.multiplication32( - this.unsigned32(this.mt[i - 1] ^ (this.mt[i - 1] >>> 30)), - 1566083941 - ) - ), - i - ); - // mt[i] &= 0xffffffff; for WORDSIZE > 32 machines - this.mt[i] = this.unsigned32(this.mt[i] & 0xffffffff); - i++; - if (i >= this.N) { - this.mt[0] = this.mt[this.N - 1]; - i = 1; - } - } +/** + * Twists the given randomness states. + * + * @param states The randomness states to twist. + * + * @returns The given states. + */ +function twist(states: number[]): number[] { + for (let idx = 0; idx !== N - M; ++idx) { + const y = (states[idx] & MASK_UPPER) + (states[idx + 1] & MASK_LOWER); + states[idx] = states[idx + M] ^ (y >>> 1) ^ (-(y & 1) & A); + } - this.mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array + for (let idx = N - M; idx !== N - 1; ++idx) { + const y = (states[idx] & MASK_UPPER) + (states[idx + 1] & MASK_LOWER); + states[idx] = states[idx + M - N] ^ (y >>> 1) ^ (-(y & 1) & A); } - // moved outside of genrandInt32() by jwatte 2010-11-17; generate less garbage - private mag01 = [0x0, this.MATRIX_A]; + const y = (states[N - 1] & MASK_UPPER) + (states[0] & MASK_LOWER); + states[N - 1] = states[M - 1] ^ (y >>> 1) ^ (-(y & 1) & A); + return states; +} +/** + * Mersenne Twister based random number generator. + */ +export class MersenneTwister19937 { /** - * Generates a random number on [0,2^32]-interval + * Creates a new Mersenne Twister random number generator based on the given seed. + * + * @param seed The seed to generate the random state from. Defaults to a random seed. */ - genrandInt32(): number { - let y: number; - - if (this.mti >= this.N) { - // generate N words at one time - let kk: number; - - // if initGenrand() has not been called a default initial seed is used - if (this.mti === this.N + 1) { - this.initGenrand(5489); - } - - for (kk = 0; kk < this.N - this.M; kk++) { - y = this.unsigned32( - (this.mt[kk] & this.UPPER_MASK) | (this.mt[kk + 1] & this.LOWER_MASK) - ); - - this.mt[kk] = this.unsigned32( - this.mt[kk + this.M] ^ (y >>> 1) ^ this.mag01[y & 0x1] - ); - } - - for (; kk < this.N - 1; kk++) { - y = this.unsigned32( - (this.mt[kk] & this.UPPER_MASK) | (this.mt[kk + 1] & this.LOWER_MASK) - ); - - this.mt[kk] = this.unsigned32( - this.mt[kk + (this.M - this.N)] ^ (y >>> 1) ^ this.mag01[y & 0x1] - ); - } - - y = this.unsigned32( - (this.mt[this.N - 1] & this.UPPER_MASK) | (this.mt[0] & this.LOWER_MASK) - ); - - this.mt[this.N - 1] = this.unsigned32( - this.mt[this.M - 1] ^ (y >>> 1) ^ this.mag01[y & 0x1] - ); - - this.mti = 0; - } - - y = this.mt[this.mti++]; - - // Tempering - y = this.unsigned32(y ^ (y >>> 11)); - y = this.unsigned32(y ^ ((y << 7) & 0x9d2c5680)); - y = this.unsigned32(y ^ ((y << 15) & 0xefc60000)); - y = this.unsigned32(y ^ (y >>> 18)); - - return y; - } + constructor(seed?: number | number[]); + /** + * Creates a new Mersenne Twister random number generator based on the given seed. + * + * @param seed The seed to generate the random state from. Defaults to a random seed. + * @param states The states to use, must be an array of 624 32-bit integers. + * @param index The index to use, must be a number between 0 and 623. + */ + constructor( + seed: number | number[] = Math.random() * Number.MAX_SAFE_INTEGER, + private states: number[] = twist(seedFrom(seed)), + private index: number = 0 + ) {} /** - * Generates a random number on [0,2^32]-interval + * Generates the next 32-bit unsigned integer. */ - genrandInt31(): number { - return this.genrandInt32() >>> 1; + nextU32(): number { + let y = this.states[this.index]; + y ^= this.states[this.index] >>> U; + y ^= (y << S) & B; + y ^= (y << T) & C; + y ^= y >>> L; + if (++this.index >= N) { + this.states = twist(this.states); + this.index = 0; + } + + return y >>> 0; } /** - * Generates a random number on [0,1]-real-interval + * Generates the next 32-bit float. */ - genrandReal1(): number { - return this.genrandInt32() * (1.0 / 4294967295.0); // divided by 2^32-1 + nextF32(): number { + return this.nextU32() * FLOATIFY_32; } /** - * Generates a random number on [0,1)-real-interval + * Generates the next 53-bit unsigned integer. */ - genrandReal2(): number { - return this.genrandInt32() * (1.0 / 4294967296.0); // divided by 2^32 + nextU53(): number { + const high = this.nextU32() >>> 5; + const low = this.nextU32() >>> 6; + return high * HIGH_MULTIPLIER_53 + low; } /** - * Generates a random number on (0,1)-real-interval + * Generates the next 53-bit float. */ - genrandReal3(): number { - return (this.genrandInt32() + 0.5) * (1.0 / 4294967296.0); // divided by 2^32 + nextF53(): number { + return this.nextU53() * FLOATIFY_53; } /** - * Generates a random number on [0,1) with 53-bit resolution + * Sets the seed of the random number generator. + * + * @param seed The seed to use. */ - genrandRes53(): number { - const a = this.genrandInt32() >>> 5, - b = this.genrandInt32() >>> 6; - return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); + seed(seed: number | number[]): void { + this.states = twist(seedFrom(seed)); + this.index = 0; } - // These real versions are due to Isaku Wada, 2002/01/09 } diff --git a/src/utils/mersenne.ts b/src/utils/mersenne.ts index 585f84b5454..be8ce2b9f7b 100644 --- a/src/utils/mersenne.ts +++ b/src/utils/mersenne.ts @@ -22,20 +22,14 @@ import type { Randomizer } from '../randomizer'; export function generateMersenne32Randomizer( seed: number = randomSeed() ): Randomizer { - const twister = new MersenneTwister19937(); - - twister.initGenrand(seed); + const twister = new MersenneTwister19937(seed); return { next(): number { - return twister.genrandReal2(); + return twister.nextF32(); }, seed(seed: number | number[]): void { - if (typeof seed === 'number') { - twister.initGenrand(seed); - } else if (Array.isArray(seed)) { - twister.initByArray(seed, seed.length); - } + twister.seed(seed); }, }; } @@ -60,20 +54,14 @@ export function generateMersenne32Randomizer( export function generateMersenne53Randomizer( seed: number = randomSeed() ): Randomizer { - const twister = new MersenneTwister19937(); - - twister.initGenrand(seed); + const twister = new MersenneTwister19937(seed); return { next(): number { - return twister.genrandRes53(); + return twister.nextF53(); }, seed(seed: number | number[]): void { - if (typeof seed === 'number') { - twister.initGenrand(seed); - } else if (Array.isArray(seed)) { - twister.initByArray(seed, seed.length); - } + twister.seed(seed); }, }; } diff --git a/test/utils/mersenne.spec.ts b/test/utils/mersenne.spec.ts index fe1b4bbc2cb..a5484b5c395 100644 --- a/test/utils/mersenne.spec.ts +++ b/test/utils/mersenne.spec.ts @@ -16,69 +16,63 @@ import { const NON_SEEDED_BASED_RUN = 25; -function newTwister(seed: number = randomSeed()): MersenneTwister19937 { - const twister = new MersenneTwister19937(); - twister.initGenrand(seed); - return twister; -} - describe('MersenneTwister19937', () => { - describe('genrandInt32()', () => { + describe('nextU32()', () => { it('should be able to return 0', () => { - const twister = newTwister(257678572); + const twister = new MersenneTwister19937(257678572); // There is no single value seed that can produce 0 in the first call for (let i = 0; i < 5; i++) { - twister.genrandInt32(); + twister.nextU32(); } - const actual = twister.genrandInt32(); + const actual = twister.nextU32(); expect(actual).toBe(0); }); it('should be able to return 2^32-1', () => { - const twister = newTwister(2855577693); - const actual = twister.genrandInt32(); + const twister = new MersenneTwister19937(2855577693); + const actual = twister.nextU32(); expect(actual).toBe(2 ** 32 - 1); }); }); - describe('genrandReal2()', () => { + describe('nextF32()', () => { it('should be able to return 0', () => { - const twister = newTwister(); + const twister = new MersenneTwister19937(); // shortcut to return minimal value // the test above shows that it is possible to return 0 - twister.genrandInt32 = () => 0; - const actual = twister.genrandReal2(); + twister.nextU32 = () => 0; + const actual = twister.nextF32(); expect(actual).toBe(0); }); it('should be able to return almost 1', () => { - const twister = newTwister(); + const twister = new MersenneTwister19937(); // shortcut to return maximal value // the test above shows that it is possible to return 2^32-1 - twister.genrandInt32 = () => 2 ** 32 - 1; - const actual = twister.genrandReal2(); + twister.nextU32 = () => 2 ** 32 - 1; + const actual = twister.nextF32(); expect(actual).toBe(TWISTER_32CO_MAX_VALUE); }); }); - describe('genrandRes53()', () => { + describe('nextF53()', () => { it('should be able to return 0', () => { - const twister = newTwister(); + const twister = new MersenneTwister19937(); // shortcut to return minimal value // the test above shows that it is possible to return 0 - twister.genrandInt32 = () => 0; - const actual = twister.genrandRes53(); + twister.nextU32 = () => 0; + const actual = twister.nextF53(); expect(actual).toBe(0); }); it('should be able to return almost 1', () => { - const twister = newTwister(); + const twister = new MersenneTwister19937(); // shortcut to return maximal value // the test above shows that it is possible to return 2^32-1 - twister.genrandInt32 = () => 2 ** 32 - 1; - const actual = twister.genrandRes53(); + twister.nextU32 = () => 2 ** 32 - 1; + const actual = twister.nextF53(); expect(actual).toBe(TWISTER_53CO_MAX_VALUE); }); });