randmts.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
- /* Mersenne Twister pseudo-random number generator functions.
- Copyright 2002, 2003 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 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.
- 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 Lesser General Public
- License for more details.
- You should have received a copy of the GNU Lesser General Public License
- along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
- #include "gmp.h"
- #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_srcptr b_orig)
- {
- mpz_t t, b;
- unsigned long e = 0x40118124;
- unsigned long bit = 0x20000000;
- mpz_init (t);
- mpz_init_set (b, b_orig); /* in case r==b_orig */
- mpz_set (r, b);
- do
- {
- mpz_mul (r, r, r);
- reduce:
- for (;;)
- {
- mpz_tdiv_q_2exp (t, r, 19937L);
- if (mpz_sgn (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.
- */
- static void
- randseed_mt (gmp_randstate_t 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_init (mod);
- mpz_init (seed1);
- mpz_set_ui (mod, 0L);
- mpz_setbit (mod, 19937L);
- mpz_sub_ui (mod, mod, 20027L);
- mpz_mod (seed1, seed, mod); /* Reduce `seed' modulo `mod'. */
- mpz_add_ui (seed1, seed1, 2L); /* seed1 is now ready. */
- mangle_seed (seed1, 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);
- cnt++;
- ASSERT (cnt <= N);
- while (cnt < N)
- p->mt[cnt++] = 0;
- mpz_clear (mod);
- mpz_clear (seed1);
- /* 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_t rstate)
- {
- __gmp_randinit_mt_noseed (rstate);
- RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
- }