New code for primality testing

[Marco Bodrato]

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

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