Improvement to mpz_nextprime

[Décio Luiz Gazzoni Filho]

On Jun 19, 2006, at 12:04 PM, Torbjorn Granlund wrote:
> Décio Luiz Gazzoni Filho <decio at decpp.net> writes:
>
>   Has anyone considered implementing a sieve for mpz_nextprime?
>
> The commented-out code in mpz/nextprime.c is a somewhat backwards
> sieve.
>
> I have made an unfinished effort at implementing a real sieve.  I
> have attached my stuff at the end.  Perhaps it is not useful.
>
>   Estimate the distance to the next prime (crude estimate is log n,
>   would have to look up better estimates), add some experimentally-
>   determined fudge factor, sieve to a reasonable amount and apply
>   mpz_probab_prime() to the survivors. If no primes are found,
>   either  keep sieving (probably with a smaller range), or revert
>   to the  original sequential algorithm, whichever works best in
>   practice.
>
> Makes sense.  Make sure that, if you resieving with a smaller size,
> you don't keep decreasing it so that no forward progress is made.
>
> I'd suggest something along these lines:
>
> sievesize = 1.5 * log (n)
> small_sievesize = sievesize / 2
> if (sievesize > LIMIT_TO_CONSERVE_MEMORY)
>   {
>     sievesize = LIMIT_TO_CONSERVE_MEMORY;
>     small_sievesize = LIMIT_TO_CONSERVE_MEMORY;
>   }
> success = try_sieve (p, n, sievesize);
> if (success)
>   return p;
> n += sievesize;
> for (;;)
>   {
>     success = try_sieve (p, n, small_sievesize);
>     if (success)
>       return p;
>     n += small_sievesize;
>   }
>
> The next question is how large primes to include in the sieve, and how
> to compute these primes.  There have been plans to include a prime
> list in GMP; such a list will have to be quite limited in size.

There are some storage tricks, of course. For instance, using a  
bitmap and a wheel of 2*3*5*7 = 210, one could store primes up to  
100k in less than 3 KB of memory. I have a gut feeling that's enough  
for the sieve, as long as people stick to sensible prime sizes -- for  
really huge stuff they're better off using e.g. PFGW.

> For small enough n, one would include an n dependent number of sieve
> primes, over a limit one need to cap the primes and then run a pseudo
> prime test for the remaining numbers.

Here's an idea for determining the cutoff point of the sieve: about 1/ 
p of surviving composites are likely to be sieved out by a given  
prime p; multiply by the time it would take to perform the sieve, and  
whenever this is less than the time a pseudoprime test in that range  
would take, stop sieving and begin PrP testing. Since sieve lengths  
are unlkely to be larger than the L2 cache, this time could be  
estimated with decent precision.

The only remaining question is how to obtain the starting residues  
for the sieve (e.g. in a call to mpz_nextprime(10^1000), one would  
have to compute 10^1000 mod 2, 3, 5, 7, 11, etc.) If subquadratic  
division is implemented, then the best strategy is, I believe:

-Compute m = 2*3*5*7*11*... up until the number is about half the  
size of 10^1000;
-Compute r = 10^1000 mod m;
-Break up m in two similarly sized halves m_1 and m_2;
-Compute r_1 = r mod m_1 and r_2 = r mod m_2;
-Proceed recursively until the modulos are prime.

The only problem I see with this scheme is memory usage when  
accumulating products. Would that be a concern?

> Also mpz_probab_prime_p needs an overhaul, perhaps at the same time.
> As Paul suggests, we want a variant of mpz_nextprime with an
> repetition argument.  Or, perhaps it is better to add new variants of
> mpz_probab_prime_p and mpz_nextprime that accepts some sort of
> probability argument, allowing us to more smoothly change from Miller
> Rabin to a stronger underlying test.

Also, perhaps implement a stronger test like the Frobenius test?

Before you ask, I'm not volunteering for this part, just the sieve. (:

Décio


-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 186 bytes
Desc: This is a digitally signed message part
Url : http://gmplib.org/list-archives/gmp-discuss/attachments/20060619/e9f17afa/PGP.bin