State of PRNG code in GMP

[Marco Bodrato]

Ciao,

Il 2020-06-02 11:22 tg at gmplib.org ha scritto:
> I'd like to have a coherent interface which also provide reentrant
> random functions to the mpn layer.  And in no case should mpn pull in
> mpz.

Makes sense.

> Unfortunately, Mersenne Twister uses mpz, and I am not saying that that
> was a bad choice when you implemented it.  But in order to provide
> reentrant mpn PRNG functions, we either rewrite the relevant parts of 
> MT
> to use mpn, or exclude it from a new mpn PRNG interface.

Mersenne Twister only uses mpz for initialization. Moreover there is a 
little "bug" in the initialization procedure, so that the sequence can 
be the same even if the seed is different (in the range qhere it is 
supposed to generate different sequencese).

That's wht some years ago we started rewriting that init function.
Of course this will yield to different sequences too.

https://gmplib.org/list-archives/gmp-bugs/2017-March/004106.html

Here is the almost mpz-free init function using the xxtea scrambler.



#define MX (((z >> 5 ^ y << 2) + (y >> 3 ^ z << 4)) \
	    ^ ((sum ^ y) + (k[(p & 3) ^ e] ^ z)))
#define DELTA 0x9E3779B9

static void
mangle_buf (gmp_uint_least32_t *buf, int high_bit)
{
   /* Apply XXTEA decryption with a constant key.
      XXTEA offers good randomness and speed properties.  */
   static const gmp_uint_least32_t key[8] = {0x4BEDAF6D, 0x674DD5FB, 
0xB79D42BC, 0x94C371EA,
					    0x5443092C, 0xA67C9FE2, 0x31CC686A, 0xC41175D6};
   const gmp_uint_least32_t *k;
   gmp_uint_least32_t z;
   gmp_uint_least32_t y = buf[0] & 0xFFFFFFFF;
   gmp_uint_least32_t sum = ((DELTA << 1) + DELTA) << 1; /* DELTA*6 */
   int p, e;
   k = high_bit ? key + 4 : key;
   do {
     e = sum >> 2 & 3;
     p = N - 2;
     do {
       z = buf[p - 1] & 0xFFFFFFFF;
       y = (buf[p] -= MX) & 0xFFFFFFFF;
     } while (--p);
     z = buf[N - 2] & 0xFFFFFFFF;
     y = (buf[0] -= MX) & 0xFFFFFFFF;
   } while (sum -= DELTA);
}

/* Seeding function.  */
static void
xxtea_randseed_mt (gmp_randstate_t rstate, mpz_srcptr seed)
{
   int i, high_bit;
   size_t cnt;

   gmp_rand_mt_struct *p;
   mpz_t seed1;  /* Intermediate result.  */
   mp_limb_t mp_d[BITS_TO_LIMBS (19936)];

   high_bit = mpz_tstbit (seed, 19936);
   p = (gmp_rand_mt_struct *) RNG_STATE (rstate);

   PTR (seed1) = mp_d;
   ALLOC (seed1) = BITS_TO_LIMBS (19936);
   SIZ (seed1) = 0;

   mpz_fdiv_r_2exp (seed1, seed, 19936);
   mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0,
	      CHAR_BIT * sizeof (p->mt[1]) - 32, seed1);

   cnt++;
   ASSERT (cnt <= N);
   while (cnt < N)
     p->mt[cnt++] = 0;

   mangle_buf (&p->mt[1], high_bit);	/* Perform the mangling of the 
buffer. */

   p->mt[0] = 0x80000000;	/* Set the first bit. */
   if (!high_bit) {
     cnt = N - 1;
     do {
       if (p->mt[cnt] != 0) {
	p->mt[0] = 0;		/* Unset the first bit. */
	break;
       }
     } while (--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;
}


> Here is some code I have tinkered with.

> typedef enum {
>   GMP_PRNG_ALG_LC,  GMP_PRNG_ALG_MT,  GMP_PRNG_ALG_AES,
>   GMP_PRNG_ALG_LIM,  GMP_PRNG_ALG_DEFAULT = GMP_PRNG_ALG_AES
> } gmp_prng_alg;

What's LIM?

Ĝis,
m