Problem with gmp_randinit_set

[Niels Möller]

nisse at lysator.liu.se (Niels Möller) writes:

> nisse at lysator.liu.se (Niels Möller) writes:
>
> Or if we want to take advantage of the structure, we need an mpn
> function to reduce numbers modulo 2^19937 - 20023.

Below is a sketch for the 64-bit case, not yet working. These things are
a bit tricky to get right, but it's not very complex code either.

Regards,
/Niels

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

#if GMP_NUMB_BITS != 64
#error unsupported limb size
#endif

#define SIZE 312
#define K 20023
/* B^312 mod p */
#define B312MODP ((mp_limb_t) K << 31)

/* Reduces r modulo p in place, result is 312 limbs (non-canonical
   representation). */
static void 
reduce (mp_ptr rp, mp_size_t n)
{
  mp_limb_t cy, hi;
  assert (n >= SIZE);
  if (n == SIZE)
    return;

  while (n >= 2*SIZE)
    {
      rp[n - SIZE] = mpn_addmul_1 (rp + n - 2*SIZE, 
				   rp + n - SIZE, SIZE, B312MODP);
      n -= (SIZE-1);
    }
  
  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[SIZE - 1];
  rp[SIZE - 1] = cy + (hi & (((mp_limb_t)1<<31) - 1))
    + mpn_add_1 (rp, rp, SIZE - 1, (hi >> 31) * K);
}

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_set_ui (p, 1);
  mpz_mul_2exp (p, p, 19937);
  mpz_sub_ui (p, p, K);

  for (i = 0;  i < 1000; 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 for size %d: %Zx\n"
		      " got: %Zx\n"
		      " exp: %Zx\n",
		      mpz_size (x), x, r, ref);
	  exit (1);
	}
    }
  mpz_clears (p, x, r, ref, NULL);
  gmp_randclear(rands);
}

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.