|
// # License 1 |
|
// |
|
// - http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html |
|
// - https://gist.github.com/banksean/300494 |
|
// - https://github.com/boo1ean/mersenne-twister |
|
// - https://www.npmjs.com/package/mersenne-twister |
|
// |
|
// BSD-3 |
|
// |
|
// Coded by Takuji Nishimura and Makoto Matsumoto. |
|
// |
|
// 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) |
|
|
|
// # License 2 |
|
// |
|
// Zlib |
|
// |
|
// - https://create.stephan-brumme.com/mersenne-twister |
|
// - https://github.com/stbrumme/mersenne-twister |
|
// - https://gitlab.com/rockerest/fast-mersenne-twister |
|
// - https://www.npmjs.com/package/fast-mersenne-twister |
|
// |
|
// Additionally, the implementation in this file borrows some ideas from the |
|
// `fast-mersenne-twister` implementation, which is based on the Stephan |
|
// Brumme’s version. |
|
// |
|
// This software is provided 'as-is', without any express or implied warranty. |
|
// In no event will the authors be held liable for any damages arising from the |
|
// use of this software. |
|
// |
|
// Permission is granted to anyone to use this software for any purpose, |
|
// including commercial applications, and to alter it and redistribute it |
|
// freely, subject to the following restrictions: |
|
// |
|
// 1. The origin of this software must not be misrepresented; you must not claim |
|
// that you wrote the original software. If you use this software in a |
|
// product, an acknowledgment in the product documentation would be |
|
// appreciated but is not required. |
|
// 2. Altered source versions must be plainly marked as such, and must not be |
|
// misrepresented as being the original software. |
|
// 3. This notice may not be removed or altered from any source distribution. |
|
|
|
// # License 3 |
|
// |
|
// BSD-3 2021 Marina Miyaoka <miyaokamarina@gmail.com> (https://twitter.com/miyaokamarina) |
|
|
|
// it seems to work faster when extracted |
|
const next = (mt: Uint32Array, i: number, j: number, k: number) => { |
|
// ╭— most significant w-r bits |
|
// | ╭— least significant r bits |
|
// ↓ ↓ |
|
j = (mt[i]! & 0x80000000) | (mt[j]! & 0x7fffffff); |
|
// ↑ |
|
// ╰— reuse j to store the y value |
|
// ╭—————————————┬— this is y |
|
// | | ╭— constant vector a |
|
// ↓ ↓ ↓ |
|
mt[i] = mt[k]! ^ (j >>> 1) ^ (-(j & 0x1) & 0x9908b0df); |
|
// ↑ ↑ |
|
// ╰———————————————————————┴— instead of multiplication or array access, we’ll use a mask |
|
}; |
|
|
|
const twist = (mt: Uint32Array) => { |
|
// generate N words at one time |
|
let i = 0; |
|
|
|
while (i < 227) next(mt, i++, i, i + 396); |
|
while (i < 623) next(mt, i++, i, i - 228); |
|
|
|
next(mt, 623, 0, 396); |
|
}; |
|
|
|
export interface MersenneTwister { |
|
/** |
|
* Random 32-bit unsigned integer. |
|
*/ |
|
u32(): number; |
|
|
|
/** |
|
* Random $[0, 1]$ float with single precision. |
|
*/ |
|
f32_ii(): number; |
|
|
|
/** |
|
* Random $[0, 1)$ float with single precision. |
|
*/ |
|
f32_ix(): number; |
|
|
|
/** |
|
* Random $(0, 1)$ float with single precision. |
|
*/ |
|
f32_xx(): number; |
|
|
|
/** |
|
* Random 53-bit unsigned integer. |
|
*/ |
|
u53(): number; |
|
|
|
/** |
|
* Random $[0, 1)$ float with double precision. |
|
*/ |
|
f64_ix(): number; |
|
|
|
/** |
|
* Export the state. |
|
*/ |
|
save(): Uint32Array; |
|
} |
|
|
|
/** |
|
* Creates or restores a Mersenne Twister. |
|
* |
|
* @param seed Either a seed value, or a previously exported state. |
|
*/ |
|
export const MersenneTwister = (seed = Date.now() as number | Uint32Array): MersenneTwister => { |
|
let i = 1; |
|
let mt = new Uint32Array(624); |
|
|
|
const u32 = () => { |
|
if (i >= 624) { |
|
/* #__NOINLINE__ */ twist(mt); |
|
i = 0; |
|
} |
|
|
|
let y = mt[i++]!; |
|
|
|
/* Tempering */ |
|
y ^= y >>> 11; |
|
y ^= (y << 7) & 0x9d2c5680; |
|
y ^= (y << 15) & 0xefc60000; |
|
y ^= y >>> 18; |
|
|
|
return y >>> 0; |
|
}; |
|
|
|
// looks like constant division hacks make significant performance |
|
// difference, so don’t fuck with it, let v8 optimize everything for you |
|
const f32_ii = () => u32() / 0x0_ffff_ffff; |
|
const f32_ix = () => u32() / 0x1_0000_0000; |
|
const f32_xx = () => (u32() + 0.5) / 0x1_0000_0000; |
|
|
|
const u53 = () => (u32() >>> 5) * 67108864 + (u32() >>> 6); |
|
|
|
const f64_ix = () => u53() / 0x20_0000_0000_0000; // = max safe integer + 1 = 2^53 |
|
|
|
const save = () => { |
|
let dump = new Uint32Array(625); |
|
|
|
dump[0] = i; |
|
dump.set(mt, 1); |
|
|
|
return dump; |
|
}; |
|
|
|
if (typeof seed === 'number') { |
|
mt[0] = seed; |
|
|
|
while (i < 624) { |
|
seed = mt[i - 1]! ^ (mt[i - 1]! >>> 30); |
|
// ↑ |
|
// ╰— reuse seed to store the s value |
|
// ╭————————————————————————————————————┬— this is s |
|
// ↓ ↓ |
|
mt[i] = (((seed >>> 16) * 1812433253) << 16) + (seed & 0xffff) * 1812433253 + i++; |
|
// 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 |
|
} |
|
} else { |
|
i = seed[0]!; |
|
mt.set(seed.slice(1)); |
|
} |
|
|
|
return { u32, f32_ii, f32_ix, f32_xx, u53, f64_ix, save }; |
|
}; |