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/
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
