github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/rand/randmts.c (about)

     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  }