mini-gmp: mpz_congruent_p and mpz_probab_prime_p

[bodrato at mail.dm.unipi.it]

Ciao,

Il Lun, 3 Marzo 2014 9:16 pm, Niels M�ller ha scritto:
> Below is an attempt at implementing it (compile-tested only). Per
> Kurlberg suggested using bases k^2 + k + 41 instead of the sequence I

The good old Euler's prime-generating polynomial... http://oeis.org/A007634

> /* Primality testing */
>
> static int
> gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
> 		 const mpz_t q, mp_bitcnt_t k)
...
>   for (i = 1; i < k; i++)
>     {
>       mpz_powm_ui (y, y, 2, n);
>       if (mpz_cmp (y, nm1) == 0)
> 	return 1;
>       if (mpz_cmp_ui (y, 1) == 0)
> 	return 0;

May I suggest to change the second comparison to mpz_cmp_ui (y, 1) <= 0?

> /* This product is 0xc0cfd797, and fits in 32 bits. */
> #define GMP_PRIME_PRODUCT \
>   (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)

Then I suggest to simplify the first few lines of prob_prime with:

/* Bit (p-1)/2 is set, for each odd prime below 67 */
#define GMP_PRIME_MASK 0x64b4cb6eUL

int
mpz_probab_prime_p (const mpz_t n, int reps)
{
  unsigned x;
  mpz_t nm1;
  mpz_t q;
  mpz_t y;
  mp_bitcnt_t k;

  if (mpz_even_p (n))
    return mpz_cmpabs_ui (n, 2) ? 0 : 2;

  if (mpz_cmpabs_ui (n, 66) < 0)
    return 2 * (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1) & 1);

  if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
    return 0;

  /* Ok, now n is odd, and all factors are >= 31. */
  if (mpz_cmpabs_ui (n, 31*31) < 0)
    return 2;

>   /* Use Miller-Rabin, with a deterministic sequence of bases, based
>      on 2^k (mod 509), where conveniently 512 = 3 (mod 509). Adding 3
etc...

I "kind of" tested it by substituting the isprime function in bootstrap.c
with the following:

int
isprime (unsigned long int t)
{
  mpz_t z;
  return mpz_probab_prime_p (mpz_roinit_n (z, &t, 1), 25);
}

and it worked...

Regards,
m

-- 
http://bodrato.it/