1/* Mersenne Twister pseudo-random number generator functions. 2 3Copyright 2002, 2003, 2013, 2014 Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library. 6 7The GNU MP Library is free software; you can redistribute it and/or modify 8it under the terms of either: 9 10 * the GNU Lesser General Public License as published by the Free 11 Software Foundation; either version 3 of the License, or (at your 12 option) any later version. 13 14or 15 16 * the GNU General Public License as published by the Free Software 17 Foundation; either version 2 of the License, or (at your option) any 18 later version. 19 20or both in parallel, as here. 21 22The GNU MP Library is distributed in the hope that it will be useful, but 23WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 24or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 25for more details. 26 27You should have received copies of the GNU General Public License and the 28GNU Lesser General Public License along with the GNU MP Library. If not, 29see https://www.gnu.org/licenses/. */ 30 31#include "gmp-impl.h" 32#include "randmt.h" 33 34 35/* Calculate (b^e) mod (2^n-k) for e=1074888996, n=19937 and k=20023, 36 needed by the seeding function below. */ 37static void 38mangle_seed (mpz_ptr r) 39{ 40 mpz_t t, b; 41 unsigned long e = 0x40118124; 42 unsigned long bit = 0x20000000; 43 44 mpz_init2 (t, 19937L); 45 mpz_init_set (b, r); 46 47 do 48 { 49 mpz_mul (r, r, r); 50 51 reduce: 52 for (;;) 53 { 54 mpz_tdiv_q_2exp (t, r, 19937L); 55 if (SIZ (t) == 0) 56 break; 57 mpz_tdiv_r_2exp (r, r, 19937L); 58 mpz_addmul_ui (r, t, 20023L); 59 } 60 61 if ((e & bit) != 0) 62 { 63 e ^= bit; 64 mpz_mul (r, r, b); 65 goto reduce; 66 } 67 68 bit >>= 1; 69 } 70 while (bit != 0); 71 72 mpz_clear (t); 73 mpz_clear (b); 74} 75 76 77/* Seeding function. Uses powering modulo a non-Mersenne prime to obtain 78 a permutation of the input seed space. The modulus is 2^19937-20023, 79 which is probably prime. The power is 1074888996. In order to avoid 80 seeds 0 and 1 generating invalid or strange output, the input seed is 81 first manipulated as follows: 82 83 seed1 = seed mod (2^19937-20027) + 2 84 85 so that seed1 lies between 2 and 2^19937-20026 inclusive. Then the 86 powering is performed as follows: 87 88 seed2 = (seed1^1074888996) mod (2^19937-20023) 89 90 and then seed2 is used to bootstrap the buffer. 91 92 This method aims to give guarantees that: 93 a) seed2 will never be zero, 94 b) seed2 will very seldom have a very low population of ones in its 95 binary representation, and 96 c) every seed between 0 and 2^19937-20028 (inclusive) will yield a 97 different sequence. 98 99 CAVEATS: 100 101 The period of the seeding function is 2^19937-20027. This means that 102 with seeds 2^19937-20027, 2^19937-20026, ... the exact same sequences 103 are obtained as with seeds 0, 1, etc.; it also means that seed -1 104 produces the same sequence as seed 2^19937-20028, etc. 105 */ 106 107static void 108randseed_mt (gmp_randstate_t rstate, mpz_srcptr seed) 109{ 110 int i; 111 size_t cnt; 112 113 gmp_rand_mt_struct *p; 114 mpz_t mod; /* Modulus. */ 115 mpz_t seed1; /* Intermediate result. */ 116 117 p = (gmp_rand_mt_struct *) RNG_STATE (rstate); 118 119 mpz_init2 (mod, 19938L); 120 mpz_init2 (seed1, 19937L); 121 122 mpz_setbit (mod, 19937L); 123 mpz_sub_ui (mod, mod, 20027L); 124 mpz_mod (seed1, seed, mod); /* Reduce `seed' modulo `mod'. */ 125 mpz_clear (mod); 126 mpz_add_ui (seed1, seed1, 2L); /* seed1 is now ready. */ 127 mangle_seed (seed1); /* Perform the mangling by powering. */ 128 129 /* Copy the last bit into bit 31 of mt[0] and clear it. */ 130 p->mt[0] = (mpz_tstbit (seed1, 19936L) != 0) ? 0x80000000 : 0; 131 mpz_clrbit (seed1, 19936L); 132 133 /* Split seed1 into N-1 32-bit chunks. */ 134 mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0, 135 8 * sizeof (p->mt[1]) - 32, seed1); 136 mpz_clear (seed1); 137 cnt++; 138 ASSERT (cnt <= N); 139 while (cnt < N) 140 p->mt[cnt++] = 0; 141 142 /* Warm the generator up if necessary. */ 143 if (WARM_UP != 0) 144 for (i = 0; i < WARM_UP / N; i++) 145 __gmp_mt_recalc_buffer (p->mt); 146 147 p->mti = WARM_UP % N; 148} 149 150 151static const gmp_randfnptr_t Mersenne_Twister_Generator = { 152 randseed_mt, 153 __gmp_randget_mt, 154 __gmp_randclear_mt, 155 __gmp_randiset_mt 156}; 157 158/* Initialize MT-specific data. */ 159void 160gmp_randinit_mt (gmp_randstate_t rstate) 161{ 162 __gmp_randinit_mt_noseed (rstate); 163 RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator; 164} 165