# mini-gmp: mpz_congruent_p and mpz_probab_prime_p

bodrato at mail.dm.unipi.it bodrato at mail.dm.unipi.it
Tue Mar 4 19:35:02 UTC 2014

```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 */

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/

```