|
Packit |
5c3484 |
/* gmp_nextprime -- generate small primes reasonably efficiently for internal
|
|
Packit |
5c3484 |
GMP needs.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Contributed to the GNU project by Torbjorn Granlund. Miscellaneous
|
|
Packit |
5c3484 |
improvements by Martin Boij.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
|
|
Packit |
5c3484 |
SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
|
|
Packit |
5c3484 |
GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Copyright 2009 Free Software Foundation, Inc.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
This file is part of the GNU MP Library.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
Packit |
5c3484 |
it under the terms of either:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU Lesser General Public License as published by the Free
|
|
Packit |
5c3484 |
Software Foundation; either version 3 of the License, or (at your
|
|
Packit |
5c3484 |
option) any later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU General Public License as published by the Free Software
|
|
Packit |
5c3484 |
Foundation; either version 2 of the License, or (at your option) any
|
|
Packit |
5c3484 |
later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or both in parallel, as here.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is distributed in the hope that it will be useful, but
|
|
Packit |
5c3484 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
Packit |
5c3484 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
5c3484 |
for more details.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
You should have received copies of the GNU General Public License and the
|
|
Packit |
5c3484 |
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
Packit |
5c3484 |
see https://www.gnu.org/licenses/. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/*
|
|
Packit |
5c3484 |
Optimisation ideas:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
1. Unroll the sieving loops. Should reach 1 write/cycle. That would be a 2x
|
|
Packit |
5c3484 |
improvement.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
2. Separate sieving with primes p < SIEVESIZE and p >= SIEVESIZE. The latter
|
|
Packit |
5c3484 |
will need at most one write, and thus not need any inner loop.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
3. For primes p >= SIEVESIZE, i.e., typically the majority of primes, we
|
|
Packit |
5c3484 |
perform more than one division per sieving write. That might dominate the
|
|
Packit |
5c3484 |
entire run time for the nextprime function. A incrementally initialised
|
|
Packit |
5c3484 |
remainder table of Pi(65536) = 6542 16-bit entries could replace that
|
|
Packit |
5c3484 |
division.
|
|
Packit |
5c3484 |
*/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#include "gmp.h"
|
|
Packit |
5c3484 |
#include "gmp-impl.h"
|
|
Packit |
5c3484 |
#include <string.h> /* for memset */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
unsigned long int
|
|
Packit |
5c3484 |
gmp_nextprime (gmp_primesieve_t *ps)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
unsigned long p, d, pi;
|
|
Packit |
5c3484 |
unsigned char *sp;
|
|
Packit |
5c3484 |
static unsigned char addtab[] =
|
|
Packit |
5c3484 |
{ 2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,
|
|
Packit |
5c3484 |
2,4,6,2,6,4,2,4,2,10,2,10 };
|
|
Packit |
5c3484 |
unsigned char *addp = addtab;
|
|
Packit |
5c3484 |
unsigned long ai;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Look for already sieved primes. A sentinel at the end of the sieving
|
|
Packit |
5c3484 |
area allows us to use a very simple loop here. */
|
|
Packit |
5c3484 |
d = ps->d;
|
|
Packit |
5c3484 |
sp = ps->s + d;
|
|
Packit |
5c3484 |
while (*sp != 0)
|
|
Packit |
5c3484 |
sp++;
|
|
Packit |
5c3484 |
if (sp != ps->s + SIEVESIZE)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
d = sp - ps->s;
|
|
Packit |
5c3484 |
ps->d = d + 1;
|
|
Packit |
5c3484 |
return ps->s0 + 2 * d;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Handle the number 2 separately. */
|
|
Packit |
5c3484 |
if (ps->s0 < 3)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
ps->s0 = 3 - 2 * SIEVESIZE; /* Tricky */
|
|
Packit |
5c3484 |
return 2;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Exhausted computed primes. Resieve, then call ourselves recursively. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if 0
|
|
Packit |
5c3484 |
for (sp = ps->s; sp < ps->s + SIEVESIZE; sp++)
|
|
Packit |
5c3484 |
*sp = 0;
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
memset (ps->s, 0, SIEVESIZE);
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ps->s0 += 2 * SIEVESIZE;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Update sqrt_s0 as needed. */
|
|
Packit |
5c3484 |
while ((ps->sqrt_s0 + 1) * (ps->sqrt_s0 + 1) <= ps->s0 + 2 * SIEVESIZE - 1)
|
|
Packit |
5c3484 |
ps->sqrt_s0++;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
pi = ((ps->s0 + 3) / 2) % 3;
|
|
Packit |
5c3484 |
if (pi > 0)
|
|
Packit |
5c3484 |
pi = 3 - pi;
|
|
Packit |
5c3484 |
if (ps->s0 + 2 * pi <= 3)
|
|
Packit |
5c3484 |
pi += 3;
|
|
Packit |
5c3484 |
sp = ps->s + pi;
|
|
Packit |
5c3484 |
while (sp < ps->s + SIEVESIZE)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
*sp = 1, sp += 3;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
pi = ((ps->s0 + 5) / 2) % 5;
|
|
Packit |
5c3484 |
if (pi > 0)
|
|
Packit |
5c3484 |
pi = 5 - pi;
|
|
Packit |
5c3484 |
if (ps->s0 + 2 * pi <= 5)
|
|
Packit |
5c3484 |
pi += 5;
|
|
Packit |
5c3484 |
sp = ps->s + pi;
|
|
Packit |
5c3484 |
while (sp < ps->s + SIEVESIZE)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
*sp = 1, sp += 5;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
pi = ((ps->s0 + 7) / 2) % 7;
|
|
Packit |
5c3484 |
if (pi > 0)
|
|
Packit |
5c3484 |
pi = 7 - pi;
|
|
Packit |
5c3484 |
if (ps->s0 + 2 * pi <= 7)
|
|
Packit |
5c3484 |
pi += 7;
|
|
Packit |
5c3484 |
sp = ps->s + pi;
|
|
Packit |
5c3484 |
while (sp < ps->s + SIEVESIZE)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
*sp = 1, sp += 7;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
p = 11;
|
|
Packit |
5c3484 |
ai = 0;
|
|
Packit |
5c3484 |
while (p <= ps->sqrt_s0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
pi = ((ps->s0 + p) / 2) % p;
|
|
Packit |
5c3484 |
if (pi > 0)
|
|
Packit |
5c3484 |
pi = p - pi;
|
|
Packit |
5c3484 |
if (ps->s0 + 2 * pi <= p)
|
|
Packit |
5c3484 |
pi += p;
|
|
Packit |
5c3484 |
sp = ps->s + pi;
|
|
Packit |
5c3484 |
while (sp < ps->s + SIEVESIZE)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
*sp = 1, sp += p;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
p += addp[ai];
|
|
Packit |
5c3484 |
ai = (ai + 1) % 48;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
ps->d = 0;
|
|
Packit |
5c3484 |
return gmp_nextprime (ps);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
void
|
|
Packit |
5c3484 |
gmp_init_primesieve (gmp_primesieve_t *ps)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
ps->s0 = 0;
|
|
Packit |
5c3484 |
ps->sqrt_s0 = 0;
|
|
Packit |
5c3484 |
ps->d = SIEVESIZE;
|
|
Packit |
5c3484 |
ps->s[SIEVESIZE] = 0; /* sentinel */
|
|
Packit |
5c3484 |
}
|