-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathgetprime.c
332 lines (282 loc) · 8.95 KB
/
getprime.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
/* Dynamic Eratosthenes sieve.
Copyright 2001, 2002, 2003, 2005, 2006, 2007, 2008, 2009, 2012
Paul Zimmermann, Alexander Kruppa, Dave Newman.
This file is part of the ECM Library.
The ECM Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The ECM Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the ECM Library; see the file COPYING.LIB. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#include <math.h>
#include <stdlib.h>
#ifdef OUTSIDE_LIBECM
# include "ecm-ecm.h"
#else
# include "ecm-impl.h"
#endif
/* This function returns successive odd primes, starting with 3.
To perform a loop over all primes <= B1, do the following
(compile this file with -DMAIN to count primes):
for (p = 2.0; p <= B1; p = getprime ())
{
...
}
It is slightly less efficient (1.5 to 2 times) than Dan Bernstein's
primegen library (http://cr.yp.to/primegen.html), however it is
fast enough for our usage here.
*/
/* We allow primes up to 2^53. This means len, the primes in primes[] etc.
will stay well below 2^32. */
static double offset = 0.0; /* offset for current primes, must be 0 or odd */
static int current = -1; /* index of previous prime */
static unsigned int *primes = NULL; /* table of small primes up to sqrt(p) */
static unsigned int nprimes = 0; /* length of primes[] */
static unsigned char *sieve = NULL; /* sieving table */
static int len = 0; /* length of sieving table, WITHOUT sentinel */
static unsigned int *moduli = NULL; /* offset for small primes,
moduli[i] = offset mod primes[i] */
/* sieve[i] == 1 if offset+2*i is a prime, otherwise sieve[i] == 0.
sieve has len + 1 bytes allocated, the last byte is always 1 (a sentinel).
This allows us avoid testing for the array end in the loop that looks
for the next prime in sieve[]. */
/* The last prime returned by getprime is offset + 2*current */
/* primes[] contains small primes to needed to sieve out composites in
sieve, i.e. all primes <= sqrt(offset + 2 * (len - 1)).
moduli[i] contains the smallest k so that offset+2*(len+k) is divisible
by primes[i], i.e. after advancing the sieve array by len,
sieve[moduli[i]] is divisible by primes[i].
*/
void
getprime_clear ()
{
offset = 0.0;
current = -1;
free (primes);
primes = NULL;
nprimes = 0;
free (sieve);
sieve = NULL;
len = 0;
free (moduli);
moduli = NULL;
}
/* For p > 1, return 1 if p is prime and 0 if p is not prime.
Requires that all primes <= sqrt(p) are in *primes */
static int
isprime_ui (unsigned int p, unsigned int *primes)
{
int i;
for (i = 0; primes[i] * primes[i] <= p; i++)
if (p % primes[i] == 0)
return 0;
return 1;
}
double
getprime ()
{
/* the following complex block is equivalent to:
while ((++current < len) && (sieve[current] == 0));
but is faster.
*/
if (len > 0L)
{
unsigned char *ptr = sieve + current;
while (*(++ptr) == 0);
current = ptr - sieve;
}
else
current = len;
if (current < len) /* most calls will end here */
return offset + 2.0 * (double) current;
/* otherwise we have to advance the sieve */
offset += 2.0 * (double) len;
/* first enlarge sieving table if too small */
if ((double) len * (double) len < offset && len > 0)
{
free (sieve);
len *= 2;
sieve = (unsigned char *) malloc ((len + 1) * sizeof (unsigned char));
/* assume this "small" malloc will not fail in normal usage */
if (sieve == NULL)
{
fprintf (stderr, "Cannot allocate memory in getprime\n");
exit (1);
}
}
/* now enlarge small prime table if too small */
if ((nprimes == 0) || (primes[nprimes-1] < sqrt(offset + 2*len)))
{
if (nprimes == 0) /* initialization */
{
nprimes = 1;
primes = (unsigned int *) malloc (nprimes * sizeof(unsigned int));
/* assume this "small" malloc will not fail in normal usage */
ASSERT(primes != NULL);
moduli = (unsigned int *) malloc (nprimes * sizeof(unsigned int));
/* assume this "small" malloc will not fail in normal usage */
ASSERT(moduli != NULL);
len = 1;
sieve = (unsigned char *) malloc((len + 1) *
sizeof(unsigned char)); /* len=1 here */
/* assume this "small" malloc will not fail in normal usage */
ASSERT(sieve != NULL);
offset = 5.0;
sieve[0] = 1; /* corresponding to 5 */
sieve[1] = 1; /* place the sentinel */
primes[0] = 3;
moduli[0] = 1; /* After we advance sieve[], sieve[0] will
correspond to 7 and sieve[1] to 9, which is
the smallest odd multiple of 3 */
current = -1;
return 3.0;
}
else
{
/* extend the existing table of small primes */
unsigned int i, j;
i = nprimes;
nprimes *= 2;
primes = (unsigned int *) realloc (primes, nprimes *
sizeof(unsigned int));
moduli = (unsigned int *) realloc (moduli, nprimes *
sizeof(unsigned int));
/* assume those "small" realloc's will not fail in normal usage */
ASSERT_ALWAYS(primes != NULL && moduli != NULL);
for (; i < nprimes; i++)
{
unsigned int p;
/* find next (odd) prime */
for (p = primes[i - 1] + 2; !isprime_ui (p, primes); p += 2);
primes[i] = p;
/* moduli[i] is the smallest m such that offset + 2*m = k*p */
j = (unsigned long) fmod (offset, (double) p);
j = (j == 0) ? j : p - j; /* -offset mod p */
if ((j % 2) != 0)
j += p; /* ensure j is even */
moduli[i] = j / 2;
}
}
}
/* now sieve for new primes */
{
int i, p;
unsigned int j;
/* Set sieve (including sentinel at the end) to 1 */
for (i = 0; i < len + 1; i++)
sieve[i] = 1;
for (j = 0; j < nprimes; j++)
{
p = primes[j];
for (i = moduli[j]; i < len; i += p)
sieve[i] = 0;
moduli[j] = i - len; /* for next sieving array */
}
}
current = -1;
while (sieve[++current] == 0);
ASSERT(current < len); /* otherwise we found a prime gap >= sqrt(x) around x */
return offset + 2.0 * (double) current;
}
/* Skips forward or backward in the sieve so that the next call to getprime
returns the smallest prime >= pp */
void
getprime_seek (double pp)
{
int i, p;
unsigned int j;
if (pp <= 3.)
{
getprime_clear ();
return;
}
offset = floor (pp / 2.) * 2. + 1.; /* make sure offset is odd */
/* Choose a large enough sieve array length */
for (i = 2; (double) i * (double) i < offset; i *= 2);
/* Now allocate sieving table */
if (len > 0)
free (sieve);
len = i;
sieve = (unsigned char *) malloc ((len + 1) * sizeof (unsigned char));
/* assume this "small" malloc will not fail in normal usage */
ASSERT_ALWAYS(sieve != NULL);
j = 1;
/* Find out how many small odd primes we'll need */
for (p = 5; (double)p*(double)p <= offset + (double)(2*len); p += 2)
{
for (i = 3; i*i <= p && p % i != 0; i += 2);
if (i*i <= p)
continue;
if ((double)p*(double)p < offset + (double)len)
j++;
}
/* Allocate memory for small primes */
if (nprimes != 0)
{
free (primes);
free (moduli);
}
nprimes = j;
primes = (unsigned int *) malloc (nprimes * sizeof(unsigned int));
moduli = (unsigned int *) malloc (nprimes * sizeof(unsigned int));
ASSERT_ALWAYS(primes != NULL && moduli != NULL);
/* Fill small primes and moduli arrays */
for (p = 3, j = 0; j < nprimes; p += 2)
{
for (i = 3; i*i <= p && p % i != 0; i += 2);
if (i*i <= p)
continue;
primes[j] = p;
i = (unsigned int) fmod (offset, (double)p);
i = (i == 0) ? i : p - i; /* -offset mod p */
if (i % 2 != 0)
i += p; /* ensure i is even */
moduli[j] = i / 2;
j++;
}
/* now sieve for new primes */
for (i = 0; i < len + 1; i++)
sieve[i] = 1;
for (j = 0; j < nprimes; j++)
{
p = primes[j];
for (i = moduli[j]; i < len; i += p)
sieve[i] = 0;
moduli[j] = i - len; /* for next sieving array */
}
current = -1;
}
#ifdef MAIN
int
main (int argc, char *argv[])
{
double p, B1, B2;
unsigned long pi = 0;
if (argc != 3)
{
fprintf (stderr, "Usage: getprime <bound> <bound>\n");
exit (EXIT_FAILURE);
}
B1 = atof (argv[1]);
B2 = atof (argv[2]);
if (B1 > 0.)
getprime_seek (B1);
p = 0;
if (B1 <= 2)
{
printf("2\n");
pi++;
}
for (p = getprime (); p <= B2; p = getprime (), pi++)
printf("%1.0f\n", p);
/* printf ("pi(%1.0f) - pi(%1.0f - 1) = %lu\n", B2, B1, pi); */
getprime_clear ();
return 0;
}
#endif