# New code for primality testing

Marco Bodrato bodrato at mail.dm.unipi.it
Fri Nov 16 19:03:55 UTC 2018

```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

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/
```