|
| 1 | +#![allow(clippy::unreadable_literal)] |
| 2 | + |
| 3 | +//! An MT19937 Mersenne Twister rng implementation, with the goal of being |
| 4 | +//! compatible with CPython's `_random` module. |
| 5 | +//! |
| 6 | +//! This crate was translated from the original |
| 7 | +//! [implementation](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html) |
| 8 | +//! by a team at Hiroshima University. The original content of the header of |
| 9 | +//! their implementation, along with the license, is left intact below. |
| 10 | +//! |
| 11 | +//! # mt19937ar.c header |
| 12 | +
|
| 13 | +/*! |
| 14 | + A C-program for MT19937, with initialization improved 2002/1/26. |
| 15 | + Coded by Takuji Nishimura and Makoto Matsumoto. |
| 16 | +
|
| 17 | + Before using, initialize the state by using init_genrand(seed) |
| 18 | + or init_by_array(init_key, key_length). |
| 19 | +
|
| 20 | + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, |
| 21 | + All rights reserved. |
| 22 | +
|
| 23 | + Redistribution and use in source and binary forms, with or without |
| 24 | + modification, are permitted provided that the following conditions |
| 25 | + are met: |
| 26 | +
|
| 27 | + 1. Redistributions of source code must retain the above copyright |
| 28 | + notice, this list of conditions and the following disclaimer. |
| 29 | +
|
| 30 | + 2. Redistributions in binary form must reproduce the above copyright |
| 31 | + notice, this list of conditions and the following disclaimer in the |
| 32 | + documentation and/or other materials provided with the distribution. |
| 33 | +
|
| 34 | + 3. The names of its contributors may not be used to endorse or promote |
| 35 | + products derived from this software without specific prior written |
| 36 | + permission. |
| 37 | +
|
| 38 | + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| 39 | + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| 40 | + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
| 41 | + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR |
| 42 | + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| 43 | + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| 44 | + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
| 45 | + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
| 46 | + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
| 47 | + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| 48 | + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 49 | +
|
| 50 | +
|
| 51 | + Any feedback is very welcome. |
| 52 | + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html |
| 53 | + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) |
| 54 | +*/ |
| 55 | + |
| 56 | +// this was translated from c; all rights go to copyright holders listed above |
| 57 | +// https://gist.github.com/coolreader18/b56d510f1b0551d2954d74ad289f7d2e |
| 58 | + |
| 59 | +/* Period parameters */ |
| 60 | +const N: usize = 624; |
| 61 | +const M: usize = 397; |
| 62 | +const MATRIX_A: u32 = 0x9908b0dfu32; /* constant vector a */ |
| 63 | +const UPPER_MASK: u32 = 0x80000000u32; /* most significant w-r bits */ |
| 64 | +const LOWER_MASK: u32 = 0x7fffffffu32; /* least significant r bits */ |
| 65 | + |
| 66 | +pub struct MT19937 { |
| 67 | + mt: [u32; N], /* the array for the state vector */ |
| 68 | + mti: usize, /* mti==N+1 means mt[N] is not initialized */ |
| 69 | +} |
| 70 | +impl Default for MT19937 { |
| 71 | + fn default() -> Self { |
| 72 | + MT19937 { |
| 73 | + mt: [0; N], |
| 74 | + mti: N + 1, |
| 75 | + } |
| 76 | + } |
| 77 | +} |
| 78 | +impl std::fmt::Debug for MT19937 { |
| 79 | + fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { |
| 80 | + f.pad("MT19937") |
| 81 | + } |
| 82 | +} |
| 83 | + |
| 84 | +impl MT19937 { |
| 85 | + pub fn new_with_slice_seed(init_key: &[u32]) -> Self { |
| 86 | + let mut state = Self::default(); |
| 87 | + state.seed_slice(init_key); |
| 88 | + state |
| 89 | + } |
| 90 | + |
| 91 | + /** initializes self.mt[N] with a seed */ |
| 92 | + fn seed(&mut self, s: u32) { |
| 93 | + self.mt[0] = s; |
| 94 | + self.mti = 1; |
| 95 | + while self.mti < N { |
| 96 | + self.mt[self.mti] = 1812433253u32 |
| 97 | + .wrapping_mul(self.mt[self.mti - 1] ^ (self.mt[self.mti - 1] >> 30)) |
| 98 | + + self.mti as u32; |
| 99 | + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ |
| 100 | + /* In the previous versions, MSBs of the seed affect */ |
| 101 | + /* only MSBs of the array self.mt[]. */ |
| 102 | + /* 2002/01/09 modified by Makoto Matsumoto */ |
| 103 | + self.mti += 1; |
| 104 | + } |
| 105 | + } |
| 106 | + |
| 107 | + /** initialize by an array with array-length */ |
| 108 | + /** init_key is the array for initializing keys */ |
| 109 | + /** key_length is its length */ |
| 110 | + /** slight change for C++, 2004/2/26 */ |
| 111 | + pub fn seed_slice(&mut self, init_key: &[u32]) { |
| 112 | + let mut i; |
| 113 | + let mut j; |
| 114 | + let mut k; |
| 115 | + self.seed(19650218); |
| 116 | + i = 1; |
| 117 | + j = 0; |
| 118 | + k = if N > init_key.len() { |
| 119 | + N |
| 120 | + } else { |
| 121 | + init_key.len() |
| 122 | + }; |
| 123 | + while k != 0 { |
| 124 | + self.mt[i] = (self.mt[i] |
| 125 | + ^ ((self.mt[i - 1] ^ (self.mt[i - 1] >> 30)).wrapping_mul(1664525u32))) |
| 126 | + + init_key[j] |
| 127 | + + j as u32; /* non linear */ |
| 128 | + self.mt[i] &= 0xffffffffu32; /* for WORDSIZE > 32 machines */ |
| 129 | + i += 1; |
| 130 | + j += 1; |
| 131 | + if i >= N { |
| 132 | + self.mt[0] = self.mt[N - 1]; |
| 133 | + i = 1; |
| 134 | + } |
| 135 | + if j >= init_key.len() { |
| 136 | + j = 0; |
| 137 | + } |
| 138 | + k -= 1; |
| 139 | + } |
| 140 | + k = N - 1; |
| 141 | + while k != 0 { |
| 142 | + self.mt[i] = (self.mt[i] |
| 143 | + ^ ((self.mt[i - 1] ^ (self.mt[i - 1] >> 30)).wrapping_mul(1566083941u32))) |
| 144 | + - i as u32; /* non linear */ |
| 145 | + self.mt[i] &= 0xffffffffu32; /* for WORDSIZE > 32 machines */ |
| 146 | + i += 1; |
| 147 | + if i >= N { |
| 148 | + self.mt[0] = self.mt[N - 1]; |
| 149 | + i = 1; |
| 150 | + } |
| 151 | + k -= 1; |
| 152 | + } |
| 153 | + |
| 154 | + self.mt[0] = 0x80000000u32; /* MSB is 1; assuring non-zero initial array */ |
| 155 | + } |
| 156 | + |
| 157 | + /** generates a random number on [0,0xffffffff]-interval */ |
| 158 | + fn gen_u32(&mut self) -> u32 { |
| 159 | + let mut y: u32; |
| 160 | + let mag01 = |x| if (x & 0x1) == 1 { MATRIX_A } else { 0 }; |
| 161 | + /* mag01[x] = x * MATRIX_A for x=0,1 */ |
| 162 | + |
| 163 | + if self.mti >= N { |
| 164 | + /* generate N words at one time */ |
| 165 | + |
| 166 | + if self.mti == N + 1 |
| 167 | + /* if seed() has not been called, */ |
| 168 | + { |
| 169 | + self.seed(5489u32); |
| 170 | + } /* a default initial seed is used */ |
| 171 | + |
| 172 | + for kk in 0..N - M { |
| 173 | + y = (self.mt[kk] & UPPER_MASK) | (self.mt[kk + 1] & LOWER_MASK); |
| 174 | + self.mt[kk] = self.mt[kk + M] ^ (y >> 1) ^ mag01(y); |
| 175 | + } |
| 176 | + for kk in N - M..N - 1 { |
| 177 | + y = (self.mt[kk] & UPPER_MASK) | (self.mt[kk + 1] & LOWER_MASK); |
| 178 | + self.mt[kk] = self.mt[kk.wrapping_add(M.wrapping_sub(N))] ^ (y >> 1) ^ mag01(y); |
| 179 | + } |
| 180 | + y = (self.mt[N - 1] & UPPER_MASK) | (self.mt[0] & LOWER_MASK); |
| 181 | + self.mt[N - 1] = self.mt[M - 1] ^ (y >> 1) ^ mag01(y); |
| 182 | + |
| 183 | + self.mti = 0; |
| 184 | + } |
| 185 | + |
| 186 | + y = self.mt[self.mti]; |
| 187 | + self.mti += 1; |
| 188 | + |
| 189 | + /* Tempering */ |
| 190 | + y ^= y >> 11; |
| 191 | + y ^= (y << 7) & 0x9d2c5680u32; |
| 192 | + y ^= (y << 15) & 0xefc60000u32; |
| 193 | + y ^= y >> 18; |
| 194 | + |
| 195 | + y |
| 196 | + } |
| 197 | +} |
| 198 | + |
| 199 | +/** generates a random number on [0,1) with 53-bit resolution*/ |
| 200 | +/// This generates a float with the same algorithm that CPython uses; calling |
| 201 | +/// it with an `MT19937` with a given seed returns the same as it would in CPython. |
| 202 | +/// |
| 203 | +/// e.g.: |
| 204 | +/// ``` |
| 205 | +/// let mut m = mt19937::MT19937::new_with_slice_seed(&[12345]); |
| 206 | +/// let expected: f64 = 0.416619872545341163316834354191087186336517333984375; |
| 207 | +/// assert_eq!(mt19937::gen_res53(&mut m), expected); |
| 208 | +/// ``` |
| 209 | +/// and in Python: |
| 210 | +/// ```python |
| 211 | +/// import random |
| 212 | +/// random.seed(12345) |
| 213 | +/// expected = 0.416619872545341163316834354191087186336517333984375 |
| 214 | +/// assert random.random() == expected |
| 215 | +/// ``` |
| 216 | +/// (note that CPython converts ints to slices by taking the native endian ordering |
| 217 | +/// of the underlying "BigInt" implementation, but for seeds < u32::max_value(), |
| 218 | +/// just `&[seed]` should be fine.) |
| 219 | +/// Original mt19937ar.c attribution: |
| 220 | +/** These real versions are due to Isaku Wada, 2002/01/09 added */ |
| 221 | +pub fn gen_res53<R: rand::RngCore>(rng: &mut R) -> f64 { |
| 222 | + let a = rng.next_u32() >> 5; |
| 223 | + let b = rng.next_u32() >> 6; |
| 224 | + (a as f64 * 67108864.0 + b as f64) * (1.0 / 9007199254740992.0) |
| 225 | +} |
| 226 | + |
| 227 | +impl rand::RngCore for MT19937 { |
| 228 | + fn next_u32(&mut self) -> u32 { |
| 229 | + self.gen_u32() |
| 230 | + } |
| 231 | + fn next_u64(&mut self) -> u64 { |
| 232 | + rand_core::impls::next_u64_via_u32(self) |
| 233 | + } |
| 234 | + fn fill_bytes(&mut self, dest: &mut [u8]) { |
| 235 | + rand_core::impls::fill_bytes_via_next(self, dest) |
| 236 | + } |
| 237 | + fn try_fill_bytes(&mut self, dest: &mut [u8]) -> Result<(), rand::Error> { |
| 238 | + self.fill_bytes(dest); |
| 239 | + Ok(()) |
| 240 | + } |
| 241 | +} |
| 242 | + |
| 243 | +/// Very big seed, but this is the size that CPython uses as well |
| 244 | +pub struct Seed([u32; N]); |
| 245 | +impl Default for Seed { |
| 246 | + fn default() -> Self { |
| 247 | + Seed([0; N]) |
| 248 | + } |
| 249 | +} |
| 250 | +impl AsMut<[u8]> for Seed { |
| 251 | + fn as_mut(&mut self) -> &mut [u8] { |
| 252 | + // in the rare chance that we aren't aligned well, just ignore it; the |
| 253 | + // max entropy we could lose is 48 bits |
| 254 | + unsafe { self.0.align_to_mut().1 } |
| 255 | + } |
| 256 | +} |
| 257 | +impl rand::SeedableRng for MT19937 { |
| 258 | + type Seed = Seed; |
| 259 | + fn from_seed(seed: Self::Seed) -> Self { |
| 260 | + Self::new_with_slice_seed(&seed.0) |
| 261 | + } |
| 262 | +} |
0 commit comments