New code for primality testing

[Adrien Prost-Boucle]

Hi Marco,

Thank you for this improvement of the primality testing :-)

I was working on my side on implementing BPSW with general aim at proposing this for GMP.
So I'm very happy to see interest in this from other people too.
And to see a certainly-better-than-mine implementation of these Lucas stuff xD

And I like your approach that keeps backward compatibility with better perf.
(even though the rationale about why reps-24 MR test becomes awkward)

But ...,
I would like the API to better :
- provide flexibility for users to use and experiment primality testing,
- enable implementation of deterministic primality tests.

In that direction it would be useful to keep an API function
that explicitly performs one or several Miller-Rabin tests only.

Overall it would be great to extend the API with functions like these (names are just illustration) :
  mpz_prime_mr_p(const mpz_t n, const mpz_t a)  // Miller-Rabin, just one test base a
  mpz_prime_mr2_p(const mpz_t n, int reps)      // Miller-Rabin, perform reps tests
  mpz_prime_bpsw_p(const mpz_t n)               // BPSW test
And one or both like these for the fast Lucas-Lehmer test for Mersenne numbers :
  mpz_prime_mersenne_p(const mpz_t n)  // input is the Mersenne number
  mpz_prime_mersenne2_p(int p)         // the Mersenne number is 2^p-1
And this for deterministic (= proven) primality test :
  mpz_prime_p(const mpz_t n)

Rationale :
- mpz_prime_p() would do deterministic primality testing
  - The primary interest is to provide a placeholder for deterministic primality testing API.
  - A first-draft approach could be : BPSW + the deterministic variant of Miller-Rabin test.
    (basically all is already there in GMP)
  - With a fallback to Lucas-Lehmer test if the number is a Mersenne one (check is easy).
  - Future optimizations could consider AKR, APR, ECC tests
- 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)

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

Regards,
Adrien


On Fri, 2018-11-16 at 20:03 +0100, Marco Bodrato wrote:
> Dear developers,
> 
> in the last ten day or so, I pushed some new code for primality testing 
> to our repository. The code is used by mpz_probab_prime_p() and 
> mpz_nextprime().
> 
> The library, before, tested primality by "some trial divisions, then 
> <reps> Miller-Rabin probabilistic primality tests". With a suggested 
> range for <reps> "between 15 and 50" [see 
> https://gmplib.org/manual/Number-Theoretic-Functions.html ].
> mpz_nextprime() uses <reps> = 25.
> 
> The main idea of the new code is to replace the first few Miller-Rabin 
> tests with the BPSW test. I.e. a Miller-Rabin test with the fixed base 
> 2, and a strong Lucas' test with the parameters suggested by the BPSW 
> test. I decided to replace the first 24 MR tests... but it's easy to 
> change this!
> 
> I pushed two implementations of the test.
>   - One for mini-gmp, https://gmplib.org/repo/gmp/rev/1dd99c5d0bf6 .
>   - Another one for the main library:
>     + from https://gmplib.org/repo/gmp/rev/81976c309e55 ...
>     + to https://gmplib.org/repo/gmp/rev/d417bcef18e6
> 
> They both choose the parameter D in the sequence 5, -7, 9, -11, ... such 
> that (D/n)=-1, then compute the Lucas' sequence with P=1 and Q=(1-D)/4.
> But the formulas used to compute the sequence are different.
> 
> mini-gmp uses the simplest:
>   U_{2k} <- U_k * V_k
>   V_{2k} <- V_k^2 - 2Q^k
>   Q^{2k} = (Q^k)^2
> (1 multiplication, 2 squares, 3 modular reductions per step)
> 
> GMP uses another one (adapted from "Elementary Number Theory" by Peter 
> Hackman, as suggested by Niels)
>   U_{2k} = U_{k+1}^2 - |U_{k+1} - U_k|^2
>   U_{2k+1} = U_{k+1}^2 - Q*U_k^2
> (3 squares, 2 modular reductions per step)
> 
> ... when D=5, i.e. the U sequence is Fibonacci's sequence, a special 
> function is used already ported to the mpn layer [see 
> https://gmplib.org/repo/gmp/rev/81976c309e55 ]. I optimised this 
> because:
>   - it was not too difficult to do (ideas inherited from the existing 
> fib2 code);
>   - it's worth doing (I expect, to see this used half of the times, for 
> random inputs; line 80 in 
> https://gmplib.org/devel/lcov/shell/gmp/mpz/stronglucas.c.gcov.html 
> seems to confirm this).
> 
> Maybe the functions are not easy to follow. I focused on writing helper 
> functions for the BPSW test, the functions are not "general". Moreover, 
> they sometimes play with the signs, because we have to detect 
> zero/non-zero, we always compute squares, so we don't really need to 
> keep track of the correct "sign".
> By the way, I wrote some comments in the code. If anyone is interested 
> in reading the implementation, I'll be happy to read any comment.
> 
> The implementation can be improved in many ways.
> 
> Should also mpz_lucas_mod be ported to the mpn layer? Should it return 
> V_n^2 instead of V_n ? testing can be improved ... and so on...
> 
> But the main point that SHOULD be improved is the repeated use of 
> mpz_tdiv_r () or mpn_tdiv_qr () to reduce the many intermediate results, 
> always by the same modulus...
> 
> I would also be curious to explore the possible implementation of some 
> tests for numbers with some special form. E.g. Proth's theorem, once we 
> have found a number D such that (D/n)=-1, if n is a Proth number, using 
> a Lucas' test is a waste of time...
> 
> Ĝis baldaŭ,
> m
>