1 /* gmp_nextprime -- generate small primes reasonably efficiently for internal
4 Contributed to the GNU project by Torbjorn Granlund. Miscellaneous
5 improvements by Martin Boij.
7 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
8 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
9 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
11 Copyright 2009 Free Software Foundation, Inc.
13 This file is part of the GNU MP Library.
15 The GNU MP Library is free software; you can redistribute it and/or modify
16 it under the terms of the GNU Lesser General Public License as published by
17 the Free Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
20 The GNU MP Library is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 License for more details.
25 You should have received a copy of the GNU Lesser General Public License
26 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
31 1. Unroll the sieving loops. Should reach 1 write/cycle. That would be a 2x
34 2. Separate sieving with primes p < SIEVESIZE and p >= SIEVESIZE. The latter
35 will need at most one write, and thus not need any inner loop.
37 3. For primes p >= SIEVESIZE, i.e., typically the majority of primes, we
38 perform more than one division per sieving write. That might dominate the
39 entire run time for the nextprime function. A incrementally initialised
40 remainder table of Pi(65536) = 6542 16-bit entries could replace that
46 #include <string.h> /* for memset */
50 gmp_nextprime (gmp_primesieve_t *ps)
52 unsigned long p, d, pi;
54 static unsigned char addtab[] =
55 { 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,
56 2,4,6,2,6,4,2,4,2,10,2,10 };
57 unsigned char *addp = addtab;
60 /* Look for already sieved primes. A sentinel at the end of the sieving
61 area allows us to use a very simple loop here. */
66 if (sp != ps->s + SIEVESIZE)
70 return ps->s0 + 2 * d;
73 /* Handle the number 2 separately. */
76 ps->s0 = 3 - 2 * SIEVESIZE; /* Tricky */
80 /* Exhausted computed primes. Resieve, then call ourselves recursively. */
83 for (sp = ps->s; sp < ps->s + SIEVESIZE; sp++)
86 memset (ps->s, 0, SIEVESIZE);
89 ps->s0 += 2 * SIEVESIZE;
91 /* Update sqrt_s0 as needed. */
92 while ((ps->sqrt_s0 + 1) * (ps->sqrt_s0 + 1) <= ps->s0 + 2 * SIEVESIZE - 1)
95 pi = ((ps->s0 + 3) / 2) % 3;
98 if (ps->s0 + 2 * pi <= 3)
101 while (sp < ps->s + SIEVESIZE)
106 pi = ((ps->s0 + 5) / 2) % 5;
109 if (ps->s0 + 2 * pi <= 5)
112 while (sp < ps->s + SIEVESIZE)
117 pi = ((ps->s0 + 7) / 2) % 7;
120 if (ps->s0 + 2 * pi <= 7)
123 while (sp < ps->s + SIEVESIZE)
130 while (p <= ps->sqrt_s0)
132 pi = ((ps->s0 + p) / 2) % p;
135 if (ps->s0 + 2 * pi <= p)
138 while (sp < ps->s + SIEVESIZE)
146 return gmp_nextprime (ps);
150 gmp_init_primesieve (gmp_primesieve_t *ps)
155 ps->s[SIEVESIZE] = 0; /* sentinel */