Problem with gmp_randinit_set

[Marco Bodrato]

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/