randmts.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
源码类别:

数学计算

开发平台:

Unix_Linux

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