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