# Problem with gmp_randinit_set

Marco Bodrato bodrato at mail.dm.unipi.it
Sun Feb 19 09:32:22 UTC 2017

```Ciao,

Il Dom, 19 Febbraio 2017 9:21 am, Niels Möller ha scritto:
> But shifting, as you suggest, might be simpler, using
>
>   2 B^623 = 20023 (mod p)

A possible generalization follows:

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

#define GSIZE (19937 / GMP_NUMB_BITS)
#define GSHIFT (19937 % GMP_NUMB_BITS)
#define SIZE (GSIZE + 1)
#define K 20023
/* B^312 mod p */
#define B312MODP ((mp_limb_t) K << (GMP_NUMB_BITS-GSHIFT))
#define USE_NIELS_SHORTCUT (((mp_limb_t) K >> GSHIFT) == 0)

/* Reduces r modulo p in place, result is SIZE limbs (non-canonical
representation). */
static void
reduce (mp_ptr rp, mp_size_t n)
{
mp_limb_t cy, hi;
assert (n > GSIZE);
assert (GSHIFT != 0); /* 19937 is prime */

if (USE_NIELS_SHORTCUT)
{
if (n == SIZE)
return;

while (n >= 2*SIZE)
{
rp[n - SIZE] = mpn_addmul_1 (rp + n - 2*SIZE,
rp + n - SIZE, SIZE, B312MODP);
n -= GSIZE;
}

cy = mpn_addmul_1 (rp, rp + SIZE, n - SIZE, B312MODP);
if (n < 2*SIZE - 1)
cy = mpn_add_1 (rp + n - SIZE, rp + n - SIZE, 2*SIZE - n - 1, cy);
hi = rp[GSIZE];
rp[GSIZE] = cy + (hi & (((mp_limb_t)1 << GSHIFT) - 1))
+ mpn_add_1 (rp, rp, SIZE - 1, (hi >> GSHIFT) * K);
}
else
{
while (n >= 2 * GSIZE)
{
hi = mpn_rshift (rp + n - GSIZE, rp + n - GSIZE, GSIZE, GSHIFT);
cy = mpn_addmul_1 (rp + n - 2*GSIZE, rp + n - GSIZE, GSIZE, K);
rp[n - GSIZE] = cy + (hi >> (GMP_NUMB_BITS - GSHIFT));
n -= GSIZE - 1;
}

hi = mpn_rshift (rp + GSIZE, rp + GSIZE, n - GSIZE, GSHIFT);
cy = mpn_addmul_1 (rp, rp + GSIZE, n - GSIZE, K);
cy = mpn_add_1 (rp + n - GSIZE, rp + n - GSIZE, 2*GSIZE - n, cy);
rp[GSIZE] = cy + (hi >> (GMP_NUMB_BITS - GSHIFT));
}
}

int
main (int argc, char **argv)
{
mpz_t p, x, r, ref;
int i;
gmp_randstate_t rands;

mpz_inits (p, x, r, ref, NULL);
gmp_randinit_default (rands);

mpz_setbit (p, 19937);
mpz_sub_ui (p, p, K);

for (i = 0;  i < 64000; i++)
{
mpz_urandomb (x, rands, 16);
mpz_rrandomb (x, rands, mpz_get_ui (x) + 20000);
mpz_tdiv_r (ref, x, p);
mpz_set (r, x);
reduce (mpz_limbs_modify (r, mpz_size (r)), mpz_size(r));
mpz_limbs_finish (r, SIZE);

if (!mpz_congruent_p (r, ref, p))
{
gmp_printf ("Failed (%i) for size %d: %Zx\n"
" got: %Zx\n"
" exp: %Zx\n",
i, mpz_size (x), x, r, ref);
exit (1);
}
}
mpz_clears (p, x, r, ref, NULL);
gmp_randclear(rands);
}

Best regards,
m

--
http://bodrato.it/

```