Ciao,
Il 2022-05-18 19:32 Marco Bodrato ha scritto:
Il Mer, 18 Maggio 2022 7:48 am, Niels Möller ha scritto:
Seth Troisi <brain...@gmail.com> writes:
I was reading the stronglucas code
Great!
- if (((*PTR (n) & 6) == 0) && UNLIKELY (mpz_perfect_square_p
(n)))
Our test-suite did not trigger that branch.
Now it does: https://gmplib.org/repo/gmp/rev/7ecb57e1bb82
I took the list of base-2 pseudo-primes by Jan Feitsma, published at
http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html .
Trowing the numbers in the list through the current implementation, it
was easy to find examples triggering the different branches!
Very useful!
the whole "((*PTR (n) & 6) == 0) &&" code is an optimization,
[...]
Should we avoid to repeat the check here and call the
function?
Maybe we can directly call the function...
On the other side, maybe we should avoid some calls to jacobi...
The current search for an odd D such that the Jacobi symbol (n/|D|) is
negative is performed by the following code:
D = GMP_NUMB_BITS % 16 == 0 ? (GMP_NUMB_BITS % 32 == 0 ? 17 : 15) : 5;
do
{
if (UNLIKELY (D >= maxD))
return 1;
D += 2;
jac = mpz_oddjacobi_ui (n, D);
}
while (jac == 1);
If we already checked 15, we may skip all the composite D. Should we at
least use a code that skips the multiple of 3? Something like:
unsigned Ddiff = 2;
#if GMP_NUMB_BITS % 16 == 0
const unsigned D2 = 6;
#if GMP_NUMB_BITS % 32 == 0
D = 19;
Ddiff = 4;
#else
D = 17;
#endif
#else
const unsigned D2 = 4;
D = 7;
#endif
for (;;)
{
jac = mpz_oddjacobi_ui (n, D);
if (jac != 1)
break;
if (UNLIKELY (D >= maxD))
return 1;
D += Ddiff;
Ddiff = D2 - Ddiff;
}
In the common case (GMP_NUMB_BITS % 32 == 0), I'd expect that about 50%
of the numbers entering stronglucas will use D=5 (Fibonacci), 25% D=7,
1/8 D=11, 1/16 D=13, 1/32 D=15, 1/64 D=17, so that only 1/64 should use
this code. I'd expect that D=19 will be used for one half of them,
but... should we check D=21 for the other half?
Ĝis,
m
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
