/* Mersenne Twister pseudo-random number generator functions.

Copyright 2002, 2003, 2013, 2014 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 either:

  * 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.

or

  * the GNU General Public License as published by the Free Software

    Foundation; either version 2 of the License, or (at your option) any

    later version.

or both in parallel, as here.

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 General Public License

for more details.

You should have received copies of the GNU General Public License and the

GNU Lesser General Public License along with the GNU MP Library.  If not,

see https://www.gnu.org/licenses/.  */

#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_t          t, b;

  unsigned long  e = 0x40118124;

  unsigned long  bit = 0x20000000;

  mpz_init2 (t, 19937L);

  mpz_init_set (b, r);

  do

    {

      mpz_mul (r, r, r);

    reduce:

      for (;;)

        {

          mpz_tdiv_q_2exp (t, r, 19937L);

          if (SIZ (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.

   Moreover, c) is not guaranted, there are many seeds yielding to the

   same sequence, because gcd (1074888996, 2^19937 - 20023 - 1) = 12.

   E.g. x and x'=x*19^((2^19937-20023-1) / 12) mod (2^19937-20023), if

   chosen as seed1, generate the same seed2, for every x.

   Similarly x" can be obtained from x', obtaining 12 different

   values.

 */

static void

randseed_mt (gmp_randstate_ptr 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_init2 (mod, 19938L);

  mpz_init2 (seed1, 19937L);

  mpz_setbit (mod, 19937L);

  mpz_sub_ui (mod, mod, 20027L);

  mpz_mod (seed1, seed, mod);  /* Reduce `seed' modulo `mod'.  */

  mpz_clear (mod);

  mpz_add_ui (seed1, seed1, 2L);  /* seed1 is now ready.  */

  mangle_seed (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);

  mpz_clear (seed1);

  cnt++;

  ASSERT (cnt <= N);

  while (cnt < N)

    p->mt[cnt++] = 0;

  /* 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_ptr rstate)

{

  __gmp_randinit_mt_noseed (rstate);

  RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;

}