On the improvement of the the primality test function

[Adrien Prost-Boucle]

Hello,

I have recently discovered that there exists a simple and fast way of
deciding for sure whether some "small" numbers are primes. Small means
under 341,550,071,728,321, which fits in 64 bits.
This is according to the following web page and its references:
http://primes.utm.edu/prove/prove2_3.html

After observing the results of the GMP primality test
function, mpz_probab_prime_p(), I noticed that the results are
probabilistic above 1000000 and that the algorithm SPRP descibed in the
page above are not present.
So I decided to implement it and do some benchmarking :-)

First, I implemented 32-bit and 64-bit versions of the algorithm so
tests are done the fastest possible for numbers that fit in 32-bit and
64-bit machine registers.

Then, I investigated what was the appropriate new limit under which
testing divisibility by odd numbers can be done. I lowered the GMP
limit of 1000000 by 100000. The 32-bit SPRP implementation is fast
enough to take care of the rest with speedup, while returning proven
primality result.

And finally I use the 64-bit SPRP function for all numbers that fit
into an unsigned long on my machine.

The resulting code is attached, in file pprime_p.c.
It can be used a drop-in replacement for the corresponding file in
src/mpz/ of GMP (for test only).
Note that it is _far_ from being ready to put in GMP, for reasons
discussed below.

I conducted tests to measure average speedup and accuracy of results
compared to vanilla GMP.

Attached is file primes.c and Makefile for compilation.
It's for amd64-capable systems.
Run the following command to compare against GMP:
	./primes isprimeGMP-check 20000
Run the following command to have execution time:
	./primes isprimeGMP-bench 20000
Modify the macro SETVAL in both commands to run tests on different
ranges of numbers and with some randomness.

My observations are the following:

In the 32-bit range and up to the last number whose primality is
proven, 341,550,071,728,321, the speedup varies between 2x and 10x
(while returning exact primality, not probabilistic).

Above that limit however, there is some slowdown, in the order of 30%-
40%. I am unsure where this comes from... I suspect my implementation
of the powmod64 function is not optimal, the loop can run up to 64
times due to shift, maybe there is a more optimized solution using some
macros or ASM instructions?
Or simply do a few more trial divisions to rapidly find small factors,
like already done when numbers don't fit an unsigned long?
Surely the slowdown on this range can be solved. Feedback about that
would be much appreciated!

My code in pprime_p.c is _far_ from being ready to put in GMP, in fact
it is more of a preliminary implementation for test and as support for
discussion about algorithms.

My understanding of GMP internals is quite low.
In particular, I'm often unsure about how to use utility macros such as umul_ppmm and friends, and about all macros that define on which kind of target we're compiling what should be checked about their capabilities.
Advice about that is also most welcomed!

Last point of this message, about the tests.
During implementation of course I occasionally produced bugs which I had a hard time to corner. I have a few suggestions to improve tests related to primality.

tests/mpz/t-pprime_p.c
=> I suggest adding primality test of all values between one known prime and the next prime, and ensure returned primality is correct (because we do know actual primality). Maybe in or or 2 such select ranges. This is linked to other suggestion below.

tests/mpz/t-nextprime.c, function refmpz_nextprime()
=> Make it check if an abnormally high number of increments is done.
Because this function, as well as the mpz_nextprime() under test, use mpz_probab_prime_p(). So if there is a problem in mpz_probab_prime_p() the current code may run indefinitely, if e.g. all numbers are considered non prime by mpz_probab_prime_p() above a certain number:

void
refmpz_nextprime (mpz_ptr p, mpz_srcptr t)
{
  mpz_add_ui (p, t, 1L);
  unsigned c = 0;
  while (! mpz_probab_prime_p (p, 10))
  {
    mpz_add_ui (p, p, 1L);
    c ++;
    if(c > 10000)
    {
      gmp_printf ("Error: too many tested numbers at %Zu\n", p);
      abort();
    }
  }
}

tests/mpz/t-nextprime.c, function main():
Make it print the values before abort(), when values are different.
This eases investigation:

  if (mpz_cmp (nxtp, ref_nxtp) != 0)
  {
    gmp_printf ("From %Zu\n", x);
    gmp_printf ("got  %Zu\n", nxtp);
    gmp_printf ("want %Zu\n", ref_nxtp);
    abort ();
  }

Best regards,
Adrien