stronglucas.c
Marco Bodrato
bodrato at mail.dm.unipi.it
Sun May 29 17:12:15 CEST 2022
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 <braintwo at 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
More information about the gmp-devel
mailing list