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