mini-gmp: mpz_congruent_p and mpz_probab_prime_p
Niels Möller
nisse at lysator.liu.se
Wed Mar 5 17:07:01 UTC 2014
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.
More information about the gmp-devel
mailing list