source/libs/gmp/gmp-src/rand/randmts.c

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