mini-gmp: mpz_congruent_p and mpz_probab_prime_p

[Niels Möller]

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

> I'm not sure what to do about mpz_probab_prime_p. One could implement it
> fairly simply as a Miller-Rabin test with no trial division (and some
> deterministically selected bases?!), or do some very limited trial
> division, say, checking for even numbers, and collecting the odd primes
> that fit in 32-bits.

Below is an attempt at implementing it (compile-tested only). Per
Kurlberg suggested using bases k^2 + k + 41 instead of the sequence I
came up with (3 + (2^k mod 509)). Maybe there are some other clever
choices (ideally, I'd like to have a sequence of bases with no
pseudoprimes smaller than a few thousand bits, but I have no idea if
that is possible. I haven't tried to research that). 

If I remember correctly, if the Miller-Rabin test passes for bases a and
b, it's useless to also try the base a * b. Hence, when using small
bases like here, it's desirable to have prime bases, or at least have
(almost) all bases include some new prime factor not present earlier in
the sequence of bases.

It's 63 non-comment lines, smaller than (plain) gcd, and smaller than
powm. So maybe not too unwieldy?

The very limited trial division should identify 84% of random inputs as
composite, I think, so average time to identify a composite should be
about 15% of a single powm, while the time to (probably) identify a
prime is reps times a single powm.

Regards,
/Niels

PS. And I already pushed a mpz_congruent_p. I forgot the test program in
my first push, but it should also be in now.

/* 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)
{
  mp_bitcnt_t i;

  /* Caller must initialize y to the base. */
  mpz_powm (y, y, q, n);

  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
    return 1;

  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;
    }
  return 0;
}

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

/* Bit p is set, for each of the above primes */
#define GMP_PRIME_MASK 0x208a28a8UL

int
mpz_probab_prime_p (const mpz_t n, int reps)
{
  int cmp;
  unsigned long gcd;
  unsigned x;
  mpz_t nm1;
  mpz_t q;
  mpz_t y;
  mp_bitcnt_t k;
  
  cmp = mpz_cmp_ui (n, 2);
  
  if (cmp <= 0)
    /* n < 2  ==>  non prime, n == 2  ==>  prime */
    return cmp + 1;

  if (mpz_tstbit (n, 0) == 0)
    /* Even number, and not 2 */
    return 0;

  gcd = mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT);
  if (gcd > 1)
    {
      /* If n is a prime, then gcd == n, and gcd is one of the above
	 primes */
      if (gcd > 29 || mpz_cmp_ui (n, gcd) != 0)
	return 0;
      return (GMP_PRIME_MASK >> gcd) & 1;
    }

  /* Ok, now n is odd, and all factors are >= 31. */
  if (mpz_cmp_ui (n, 31*31) < 0)    
    return 1;
  
  /* Use Miller-Rabin, with a deterministic sequence of bases, based
     on 2^k (mod 509), where conveniently 512 = 3 (mod 509). Adding 3
     to this sequence gives us a sequence of bases starting with
     4,5,7,11,19,35,67,131,259,6,... Since n > 31*31, all bases are in
     the range [4, n-2], so we always avoid -1, 0 and 1 (mod n). */ 

  mpz_init (nm1);
  mpz_init (q);
  mpz_init (y);

  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
  mpz_sub_ui (nm1, n, 1);
  k = mpz_scan1 (nm1, 0);
  mpz_tdiv_q_2exp (q, nm1, k);
  
  for (x = 1; reps > 0;
       reps--, x = ((x & 0xff) << 1) + 3*(x >> 8))
    {
      mpz_set_ui (y, 3 + x);
      if (!gmp_millerrabin (n, nm1, y, q, k))
	break;
    }
  mpz_clear (nm1);
  mpz_clear (q);
  mpz_clear (y);
    
  return reps == 0;
}

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