stronglucas.c

[Marco Bodrato]

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