New code for primality testing

[Marco Bodrato]

Ciao Adrien,

Nice to hear your comments!

Il Sab, 17 Novembre 2018 11:58 am, Adrien Prost-Boucle ha scritto:
> I was working on my side on implementing BPSW with general aim at

So, you have some possible implementation details in mind... If you want
to read the code I pushed, and ask any question like "why did you
implement this way that step?", I'll be happy to answer.

> And I like your approach that keeps backward compatibility with better
> perf.

It is intended for a GMP-6.2 release. For GMP-7 we should give different
functions, I think.

> (even though the rationale about why reps-24 MR test becomes awkward)

There isn't any good reason for the number 24, I'm ready to change it, if
other proposals arrive.

> Overall it would be great to extend the API with functions like these

>   mpz_prime_mr_p(const mpz_t n, const mpz_t a)  // Miller-Rabin, just one

Easy to do, with a wrapper around the function millerrabin already
available as an internal function in
https://gmplib.org/repo/gmp/file/17707/mpz/millerrabin.c#l130

We can assume we will always have a Miller-Rabin implementation, so that
we may decide to export it in the API....

>   mpz_prime_mr2_p(const mpz_t n, int reps)      // Miller-Rabin, perform
> reps tests

It should also take a gmp_randstate_t parameter, and chose bases in the
range [3..n/2]. But ... is this in the spirit of the API of the library?

The aim of this function is to test primality up to a given confidence, we
already have it, using MR or another algorithm is an implementation
detail.

>   mpz_prime_bpsw_p(const mpz_t n)               // BPSW test

This is equivalent to mpz_prime_mr_p (n,2) && mpz_stronglucas_p (n), where
the second function should be a wrapper around the internal function in
https://gmplib.org/repo/gmp/file/17707/mpz/stronglucas.c#l62

But again, will we keep on having BPSW also in the future releases?
Moreover... what is BPSW? Initial rial division or not? Up to which
number? Strong or weak PSP tests? I'm not sure we should export it.

> And one or both like these for the fast Lucas-Lehmer test for Mersenne

I don't agree. The set of Mersenne's numbers is too small to write special
code for it in a general library.
I wrote a comment about Lucas-Lehmer, because Mersenne's number require a
few special lines to be handled correctly...
https://gmplib.org/repo/gmp/file/17707/mpz/lucmod.c#l57
but I do not really thin it should be implemented.

We should maybe implement the more general Lucas–Lehmer–Riesel test!
A (small) brick to build it was just inserted in the code:
https://gmplib.org/repo/gmp/file/17707/mpn/generic/strongfibo.c#l69

>   mpz_prime_mersenne2_p(int p)         // the Mersenne number is 2^p-1

I'd suggest to hard-code a list of the already known primes, return 2 if
the value is in the list, 0 otherwise, and call the generic functions (at
first on the exponent, then on the number) if p is out of range :-) It
would be very fast up to some millions of digits! And too slow (as
expected) for larger numbers :-D

>   - Future optimizations could consider AKR, APR, ECC tests

...too complex for this library.
Maybe the GMP-ECM [https://gforge.inria.fr/projects/ecm/] project already
has the building blocks for Elliptic curve. GMP doesn't.

Some other free software projects already contain code for primality
proving. E.g. Pari/GP [http://pari.math.u-bordeaux.fr/] implement both
APR-CL and ECPP... I don't think GMP should duplicate that (good) work.

> - There is a lot of interest and activity around Mersenne numbers
>   And the baseline Lucas-Lehmer test is easy to implement (competing with
> GIMPS perf is not the goal)

A program to be distributed as demos/mersenne.c, maybe? There are some
demos using the mpz layer, this one might use mpn.

> What's your general point of view / intention about primality testing
> API ?

The probabilistic test should take a gmp_randstate_t parameter, so that
repeated calls to the test can give increasing confidence.
For this use-case we should also think about a way to skip the initial
checks (trial division? BPSW?) for the second calls. A flag or different
functions?

In my opinion, we should remove the prototype of the function
mpz_millerrabin from gmp.h . It was added years ago
https://gmplib.org/repo/gmp/rev/6041
it has never been documented, and was marked as internal a little later
https://gmplib.org/repo/gmp/rev/6290

I don't think the API should change much more than this...


For future implementations, it would be nice to look for some other tests
giving more confidence than Miller-Rabin...

I'd personally would like to have an implementation also for some tests
for special numbers. But not so special as Mersenne's or Fermat's...
E.g. d*2^n +/- 1 with d < 2^n

But currently I'll simply try to re-read, verify and refine the exiting
code. I'm not yet sure it's correct.

Ĝis,
m

-- 
http://bodrato.it/papers/