mini-gmp: mpz_congruent_p and mpz_probab_prime_p

[Niels Möller]

bodrato at mail.dm.unipi.it writes:

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

I hadn't heard about that before.

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

Hmm, why? y == 0 can happen if n divides a power of a (the base for the
test), in which case n is composite (except in corner cases, see below),
and in which case the test will eventually fail without explicitly
checking for y == 0.

There are some corner cases with the choice of bases, since we may get
n-1 <= a. We have n >= 31^2 = 961, so if we use a_k = k^2 + k + 41, so
this can happen if reps > 30.

Maybe it's good enough to stop the base loop when we get a >= n-1, even
if the caller specifies a large reps.

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

Makes sense to me.

Regards,
/Niels

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