1
/* List and count primes.
2
Written by tege while on holiday in Rodupp, August 2001.
3
Between 10 and 500 times faster than previous program.
5
Copyright 2001 Free Software Foundation, Inc.
7
This program is free software; you can redistribute it and/or modify it under
8
the terms of the GNU General Public License as published by the Free Software
9
Foundation; either version 2 of the License, or (at your option) any later
12
This program is distributed in the hope that it will be useful, but WITHOUT ANY
13
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
14
PARTICULAR PURPOSE. See the GNU General Public License for more details.
16
You should have received a copy of the GNU General Public License along with
17
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
18
Place - Suite 330, Boston, MA 02111-1307, USA. */
27
* Do not fill primes[] with real primes when the range [fr,to] is small,
28
when fr,to are relatively large. Fill primes[] with odd numbers instead.
29
[Probably a bad idea, since the primes[] array would become very large.]
30
* Separate small primes and large primes when sieving. Either the Montgomery
31
way (i.e., having a large array a multiple of L1 cache size), or just
32
separate loops for primes <= S and primes > S. The latter primes do not
33
require an inner loop, since they will touch the sieving array at most once.
34
* Pre-fill sieving array with an appropriately aligned ...00100100... pattern,
35
then omit 3 from primes array. (May require similar special handling of 3
36
as we now have for 2.)
37
* A large SIEVE_LIMIT currently implies very large memory usage, mainly due
38
to the sieving array in make_primelist, but also because of the primes[]
39
array. We might want to stage the program, using sieve_region/find_primes
40
to build primes[]. Make report() a function pointer, as part of achieving
42
* Store primes[] as two arrays, one array with primes represented as delta
43
values using just 8 bits (if gaps are too big, store bogus primes!)
44
and one array with "rem" values. The latter needs 32-bit values.
45
* A new entry point, mpz_probab_prime_likely_p, would be useful.
46
* Improve command line syntax and versatility. "primes -f FROM -t TO",
47
allow either to be omitted for open interval. (But disallow
48
"primes -c -f FROM" since that would be infinity.) Allow printing a
49
limited *number* of primes using syntax like "primes -f FROM -n NUMBER".
50
* When looking for maxgaps, we should not perform any primality testing until
51
we find possible record gaps. Should speed up the searches tremendously.
62
struct primes *primes;
63
unsigned long n_primes;
65
void find_primes (unsigned char *, mpz_t, unsigned long, mpz_t);
66
void sieve_region (unsigned char *, mpz_t, unsigned long);
67
void make_primelist (unsigned long);
72
unsigned long maxgap = 0;
73
unsigned long total_primes = 0;
81
mpz_out_str (stdout, 10, prime);
86
static unsigned long prev_prime_low = 0;
88
if (prev_prime_low != 0)
90
gap = mpz_get_ui (prime) - prev_prime_low;
94
prev_prime_low = mpz_get_ui (prime);
99
main (int argc, char *argv[])
101
char *progname = argv[0];
104
unsigned long sieve_lim;
105
unsigned long est_n_primes;
112
if (strcmp (argv[1], "-c") == 0)
118
else if (strcmp (argv[1], "-p") == 0)
124
else if (strcmp (argv[1], "-g") == 0)
134
if (flag_count || flag_maxgap)
135
flag_print--; /* clear unless an explicit -p */
144
mpz_set_str (fr, argv[1], 0);
145
if (argv[2][0] == '+')
147
mpz_set_str (to, argv[2] + 1, 0);
148
mpz_add (to, to, fr);
151
mpz_set_str (to, argv[2], 0);
156
mpz_set_str (to, argv[1], 0);
160
fprintf (stderr, "usage: %s [-c] [-p] [-g] [from [+]]to\n", progname);
165
if (mpz_cmp_ui (fr2, 3) < 0)
171
mpz_setbit (fr2, 0); /* make odd */
172
mpz_sub_ui (to2, to, 1);
173
mpz_setbit (to2, 0); /* make odd */
176
mpz_init (siev_sqr_lim);
179
#define SIEVE_LIMIT 10000000
180
if (mpz_cmp_ui (tmp, SIEVE_LIMIT) < 0)
182
sieve_lim = mpz_get_ui (tmp);
186
sieve_lim = SIEVE_LIMIT;
187
mpz_sub (tmp, to2, fr2);
188
if (mpz_cmp_ui (tmp, sieve_lim) < 0)
189
sieve_lim = mpz_get_ui (tmp); /* limit sieving for small ranges */
191
mpz_set_ui (siev_sqr_lim, sieve_lim + 1);
192
mpz_mul_ui (siev_sqr_lim, siev_sqr_lim, sieve_lim + 1);
194
est_n_primes = (size_t) (sieve_lim / log((double) sieve_lim) * 1.13) + 10;
195
primes = malloc (est_n_primes * sizeof primes[0]);
196
make_primelist (sieve_lim);
197
assert (est_n_primes >= n_primes);
200
printf ("sieve_lim = %lu\n", sieve_lim);
201
printf ("n_primes = %lu (3..%u)\n",
202
n_primes, primes[n_primes - 1].prime);
205
#define S (1 << 15) /* FIXME: Figure out L1 cache size */
207
while (mpz_cmp (fr2, to2) <= 0)
211
mpz_add_ui (tmp, fr2, rsize);
212
if (mpz_cmp (tmp, to2) > 0)
214
mpz_sub (tmp, to2, fr2);
215
rsize = mpz_get_ui (tmp) + 2;
218
printf ("Sieving region ["); mpz_out_str (stdout, 10, fr2);
219
printf (","); mpz_add_ui (tmp, fr2, rsize - 2);
220
mpz_out_str (stdout, 10, tmp); printf ("]\n");
222
sieve_region (s, fr2, rsize);
223
find_primes (s, fr2, rsize / 2, siev_sqr_lim);
225
mpz_add_ui (fr2, fr2, S);
230
printf ("Pi(interval) = %lu\n", total_primes);
233
printf ("max gap: %lu\n", maxgap);
238
/* Find primes in region [fr,fr+rsize). Requires that fr is odd and that
239
rsize is even. The sieving array s should be aligned for "long int" and
240
have rsize/2 entries, rounded up to the nearest multiple of "long int". */
242
sieve_region (unsigned char *s, mpz_t fr, unsigned long rsize)
244
unsigned long ssize = rsize / 2;
245
unsigned long start, start2, prime;
252
/* initialize sieving array */
253
for (ii = 0; ii < (ssize + sizeof (long) - 1) / sizeof (long); ii++)
254
((long *) s) [ii] = ~0L;
258
long *se = (long *) (s + ((ssize + sizeof (long) - 1) & -sizeof (long)));
259
for (k = -((ssize + sizeof (long) - 1) / sizeof (long)); k < 0; k++)
264
for (i = 0; i < n_primes; i++)
266
prime = primes[i].prime;
268
if (primes[i].rem >= 0)
270
start2 = primes[i].rem;
274
mpz_set_ui (tmp, prime);
275
mpz_mul_ui (tmp, tmp, prime);
276
if (mpz_cmp (fr, tmp) <= 0)
278
mpz_sub (tmp, tmp, fr);
279
if (mpz_cmp_ui (tmp, 2 * ssize) > 0)
280
break; /* avoid overflow at next line, also speedup */
281
start = mpz_get_ui (tmp);
285
start = (prime - mpz_tdiv_ui (fr, prime)) % prime;
287
start += prime; /* adjust if even divisable */
293
for (ii = start2; ii < ssize; ii += prime)
295
primes[i].rem = ii - ssize;
299
unsigned char *se = s + ssize; /* point just beyond sieving range */
300
for (k = start2 - ssize; k < 0; k += prime)
309
/* Find primes in region [fr,fr+rsize), using the previously sieved s[]. */
311
find_primes (unsigned char *s, mpz_t fr, unsigned long ssize,
318
for (j = 0; j < (ssize + sizeof (long) - 1) / sizeof (long); j++)
320
if (((long *) s) [j] != 0)
322
for (ij = 0; ij < sizeof (long); ij++)
324
if (s[j * sizeof (long) + ij] != 0)
326
if (j * sizeof (long) + ij >= ssize)
328
mpz_add_ui (tmp, fr, (j * sizeof (long) + ij) * 2);
329
if (mpz_cmp (tmp, siev_sqr_lim) < 0 ||
330
mpz_probab_prime_p (tmp, 3))
340
/* Generate a lits of primes and store in the global array primes[]. */
342
make_primelist (unsigned long maxprime)
346
unsigned long ssize = maxprime / 2;
347
unsigned long i, ii, j;
350
memset (s, ~0, ssize);
351
for (i = 3; ; i += 2)
353
unsigned long isqr = i * i;
354
if (isqr >= maxprime)
356
if (s[i * i / 2 - 1] == 0)
357
continue; /* only sieve with primes */
358
for (ii = i * i / 2 - 1; ii < ssize; ii += i)
362
for (j = 0; j < ssize; j++)
366
primes[n_primes].prime = j * 2 + 3;
367
primes[n_primes].rem = -1;
371
/* FIXME: This should not be needed if fencepost errors were fixed... */
372
if (primes[n_primes - 1].prime > maxprime)
378
for (i = 3; i <= maxprime; i += 2)
380
if (i < 7 || (i % 3 != 0 && i % 5 != 0 && i % 7 != 0))
382
primes[n_primes].prime = i;
383
primes[n_primes].rem = -1;
added
removed
gmp3/demos/primes.c
