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