From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/rand/randmts.c | 171 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 gmp-6.3.0/rand/randmts.c (limited to 'gmp-6.3.0/rand/randmts.c') diff --git a/gmp-6.3.0/rand/randmts.c b/gmp-6.3.0/rand/randmts.c new file mode 100644 index 0000000..2b95a2a --- /dev/null +++ b/gmp-6.3.0/rand/randmts.c @@ -0,0 +1,171 @@ +/* Mersenne Twister pseudo-random number generator functions. + +Copyright 2002, 2003, 2013, 2014 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp-impl.h" +#include "randmt.h" + + +/* Calculate (b^e) mod (2^n-k) for e=1074888996, n=19937 and k=20023, + needed by the seeding function below. */ +static void +mangle_seed (mpz_ptr r) +{ + mpz_t t, b; + unsigned long e = 0x40118124; + unsigned long bit = 0x20000000; + + mpz_init2 (t, 19937L); + mpz_init_set (b, r); + + do + { + mpz_mul (r, r, r); + + reduce: + for (;;) + { + mpz_tdiv_q_2exp (t, r, 19937L); + if (SIZ (t) == 0) + break; + mpz_tdiv_r_2exp (r, r, 19937L); + mpz_addmul_ui (r, t, 20023L); + } + + if ((e & bit) != 0) + { + e ^= bit; + mpz_mul (r, r, b); + goto reduce; + } + + bit >>= 1; + } + while (bit != 0); + + mpz_clear (t); + mpz_clear (b); +} + + +/* Seeding function. Uses powering modulo a non-Mersenne prime to obtain + a permutation of the input seed space. The modulus is 2^19937-20023, + which is probably prime. The power is 1074888996. In order to avoid + seeds 0 and 1 generating invalid or strange output, the input seed is + first manipulated as follows: + + seed1 = seed mod (2^19937-20027) + 2 + + so that seed1 lies between 2 and 2^19937-20026 inclusive. Then the + powering is performed as follows: + + seed2 = (seed1^1074888996) mod (2^19937-20023) + + and then seed2 is used to bootstrap the buffer. + + This method aims to give guarantees that: + a) seed2 will never be zero, + b) seed2 will very seldom have a very low population of ones in its + binary representation, and + c) every seed between 0 and 2^19937-20028 (inclusive) will yield a + different sequence. + + CAVEATS: + + The period of the seeding function is 2^19937-20027. This means that + with seeds 2^19937-20027, 2^19937-20026, ... the exact same sequences + are obtained as with seeds 0, 1, etc.; it also means that seed -1 + produces the same sequence as seed 2^19937-20028, etc. + + Moreover, c) is not guaranted, there are many seeds yielding to the + same sequence, because gcd (1074888996, 2^19937 - 20023 - 1) = 12. + E.g. x and x'=x*19^((2^19937-20023-1) / 12) mod (2^19937-20023), if + chosen as seed1, generate the same seed2, for every x. + Similarly x" can be obtained from x', obtaining 12 different + values. + */ + +static void +randseed_mt (gmp_randstate_ptr rstate, mpz_srcptr seed) +{ + int i; + size_t cnt; + + gmp_rand_mt_struct *p; + mpz_t mod; /* Modulus. */ + mpz_t seed1; /* Intermediate result. */ + + p = (gmp_rand_mt_struct *) RNG_STATE (rstate); + + mpz_init2 (mod, 19938L); + mpz_init2 (seed1, 19937L); + + mpz_setbit (mod, 19937L); + mpz_sub_ui (mod, mod, 20027L); + mpz_mod (seed1, seed, mod); /* Reduce `seed' modulo `mod'. */ + mpz_clear (mod); + mpz_add_ui (seed1, seed1, 2L); /* seed1 is now ready. */ + mangle_seed (seed1); /* Perform the mangling by powering. */ + + /* Copy the last bit into bit 31 of mt[0] and clear it. */ + p->mt[0] = (mpz_tstbit (seed1, 19936L) != 0) ? 0x80000000 : 0; + mpz_clrbit (seed1, 19936L); + + /* Split seed1 into N-1 32-bit chunks. */ + mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0, + 8 * sizeof (p->mt[1]) - 32, seed1); + mpz_clear (seed1); + cnt++; + ASSERT (cnt <= N); + while (cnt < N) + p->mt[cnt++] = 0; + + /* Warm the generator up if necessary. */ + if (WARM_UP != 0) + for (i = 0; i < WARM_UP / N; i++) + __gmp_mt_recalc_buffer (p->mt); + + p->mti = WARM_UP % N; +} + + +static const gmp_randfnptr_t Mersenne_Twister_Generator = { + randseed_mt, + __gmp_randget_mt, + __gmp_randclear_mt, + __gmp_randiset_mt +}; + +/* Initialize MT-specific data. */ +void +gmp_randinit_mt (gmp_randstate_ptr rstate) +{ + __gmp_randinit_mt_noseed (rstate); + RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator; +} -- cgit v1.2.3