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, 2002, 2006, 2012 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 3 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, see https://www.gnu.org/licenses/. */
26 * Do not fill primes[] with real primes when the range [fr,to] is small,
27 when fr,to are relatively large. Fill primes[] with odd numbers instead.
28 [Probably a bad idea, since the primes[] array would become very large.]
29 * Separate small primes and large primes when sieving. Either the Montgomery
30 way (i.e., having a large array a multiple of L1 cache size), or just
31 separate loops for primes <= S and primes > S. The latter primes do not
32 require an inner loop, since they will touch the sieving array at most once.
33 * Pre-fill sieving array with an appropriately aligned ...00100100... pattern,
34 then omit 3 from primes array. (May require similar special handling of 3
35 as we now have for 2.)
36 * A large SIEVE_LIMIT currently implies very large memory usage, mainly due
37 to the sieving array in make_primelist, but also because of the primes[]
38 array. We might want to stage the program, using sieve_region/find_primes
39 to build primes[]. Make report() a function pointer, as part of achieving
41 * Store primes[] as two arrays, one array with primes represented as delta
42 values using just 8 bits (if gaps are too big, store bogus primes!)
43 and one array with "rem" values. The latter needs 32-bit values.
44 * A new entry point, mpz_probab_prime_likely_p, would be useful.
45 * Improve command line syntax and versatility. "primes -f FROM -t TO",
46 allow either to be omitted for open interval. (But disallow
47 "primes -c -f FROM" since that would be infinity.) Allow printing a
48 limited *number* of primes using syntax like "primes -f FROM -n NUMBER".
49 * When looking for maxgaps, we should not perform any primality testing until
50 we find possible record gaps. Should speed up the searches tremendously.
61 struct primes
*primes
;
62 unsigned long n_primes
;
64 void find_primes (unsigned char *, mpz_t
, unsigned long, mpz_t
);
65 void sieve_region (unsigned char *, mpz_t
, unsigned long);
66 void make_primelist (unsigned long);
71 unsigned long maxgap
= 0;
72 unsigned long total_primes
= 0;
80 mpz_out_str (stdout
, 10, prime
);
85 static unsigned long prev_prime_low
= 0;
87 if (prev_prime_low
!= 0)
89 gap
= mpz_get_ui (prime
) - prev_prime_low
;
93 prev_prime_low
= mpz_get_ui (prime
);
98 main (int argc
, char *argv
[])
100 char *progname
= argv
[0];
103 unsigned long sieve_lim
;
104 unsigned long est_n_primes
;
111 if (strcmp (argv
[1], "-c") == 0)
117 else if (strcmp (argv
[1], "-p") == 0)
123 else if (strcmp (argv
[1], "-g") == 0)
133 if (flag_count
|| flag_maxgap
)
134 flag_print
--; /* clear unless an explicit -p */
143 mpz_set_str (fr
, argv
[1], 0);
144 if (argv
[2][0] == '+')
146 mpz_set_str (to
, argv
[2] + 1, 0);
147 mpz_add (to
, to
, fr
);
150 mpz_set_str (to
, argv
[2], 0);
155 mpz_set_str (to
, argv
[1], 0);
159 fprintf (stderr
, "usage: %s [-c] [-p] [-g] [from [+]]to\n", progname
);
164 if (mpz_cmp_ui (fr2
, 3) < 0)
170 mpz_setbit (fr2
, 0); /* make odd */
171 mpz_sub_ui (to2
, to
, 1);
172 mpz_setbit (to2
, 0); /* make odd */
175 mpz_init (siev_sqr_lim
);
178 #define SIEVE_LIMIT 10000000
179 if (mpz_cmp_ui (tmp
, SIEVE_LIMIT
) < 0)
181 sieve_lim
= mpz_get_ui (tmp
);
185 sieve_lim
= SIEVE_LIMIT
;
186 mpz_sub (tmp
, to2
, fr2
);
187 if (mpz_cmp_ui (tmp
, sieve_lim
) < 0)
188 sieve_lim
= mpz_get_ui (tmp
); /* limit sieving for small ranges */
190 mpz_set_ui (siev_sqr_lim
, sieve_lim
+ 1);
191 mpz_mul_ui (siev_sqr_lim
, siev_sqr_lim
, sieve_lim
+ 1);
193 est_n_primes
= (size_t) (sieve_lim
/ log((double) sieve_lim
) * 1.13) + 10;
194 primes
= malloc (est_n_primes
* sizeof primes
[0]);
195 make_primelist (sieve_lim
);
196 assert (est_n_primes
>= n_primes
);
199 printf ("sieve_lim = %lu\n", sieve_lim
);
200 printf ("n_primes = %lu (3..%u)\n",
201 n_primes
, primes
[n_primes
- 1].prime
);
204 #define S (1 << 15) /* FIXME: Figure out L1 cache size */
206 while (mpz_cmp (fr2
, to2
) <= 0)
210 mpz_add_ui (tmp
, fr2
, rsize
);
211 if (mpz_cmp (tmp
, to2
) > 0)
213 mpz_sub (tmp
, to2
, fr2
);
214 rsize
= mpz_get_ui (tmp
) + 2;
217 printf ("Sieving region ["); mpz_out_str (stdout
, 10, fr2
);
218 printf (","); mpz_add_ui (tmp
, fr2
, rsize
- 2);
219 mpz_out_str (stdout
, 10, tmp
); printf ("]\n");
221 sieve_region (s
, fr2
, rsize
);
222 find_primes (s
, fr2
, rsize
/ 2, siev_sqr_lim
);
224 mpz_add_ui (fr2
, fr2
, S
);
229 printf ("Pi(interval) = %lu\n", total_primes
);
232 printf ("max gap: %lu\n", maxgap
);
237 /* Find primes in region [fr,fr+rsize). Requires that fr is odd and that
238 rsize is even. The sieving array s should be aligned for "long int" and
239 have rsize/2 entries, rounded up to the nearest multiple of "long int". */
241 sieve_region (unsigned char *s
, mpz_t fr
, unsigned long rsize
)
243 unsigned long ssize
= rsize
/ 2;
244 unsigned long start
, start2
, prime
;
251 /* initialize sieving array */
252 for (ii
= 0; ii
< (ssize
+ sizeof (long) - 1) / sizeof (long); ii
++)
253 ((long *) s
) [ii
] = ~0L;
257 long *se
= (long *) (s
+ ((ssize
+ sizeof (long) - 1) & -sizeof (long)));
258 for (k
= -((ssize
+ sizeof (long) - 1) / sizeof (long)); k
< 0; k
++)
263 for (i
= 0; i
< n_primes
; i
++)
265 prime
= primes
[i
].prime
;
267 if (primes
[i
].rem
>= 0)
269 start2
= primes
[i
].rem
;
273 mpz_set_ui (tmp
, prime
);
274 mpz_mul_ui (tmp
, tmp
, prime
);
275 if (mpz_cmp (fr
, tmp
) <= 0)
277 mpz_sub (tmp
, tmp
, fr
);
278 if (mpz_cmp_ui (tmp
, 2 * ssize
) > 0)
279 break; /* avoid overflow at next line, also speedup */
280 start
= mpz_get_ui (tmp
);
284 start
= (prime
- mpz_tdiv_ui (fr
, prime
)) % prime
;
286 start
+= prime
; /* adjust if even divisible */
292 for (ii
= start2
; ii
< ssize
; ii
+= prime
)
294 primes
[i
].rem
= ii
- ssize
;
298 unsigned char *se
= s
+ ssize
; /* point just beyond sieving range */
299 for (k
= start2
- ssize
; k
< 0; k
+= prime
)
308 /* Find primes in region [fr,fr+rsize), using the previously sieved s[]. */
310 find_primes (unsigned char *s
, mpz_t fr
, unsigned long ssize
,
317 for (j
= 0; j
< (ssize
+ sizeof (long) - 1) / sizeof (long); j
++)
319 if (((long *) s
) [j
] != 0)
321 for (ij
= 0; ij
< sizeof (long); ij
++)
323 if (s
[j
* sizeof (long) + ij
] != 0)
325 if (j
* sizeof (long) + ij
>= ssize
)
327 mpz_add_ui (tmp
, fr
, (j
* sizeof (long) + ij
) * 2);
328 if (mpz_cmp (tmp
, siev_sqr_lim
) < 0 ||
329 mpz_probab_prime_p (tmp
, 10))
339 /* Generate a list of primes and store in the global array primes[]. */
341 make_primelist (unsigned long maxprime
)
345 unsigned long ssize
= maxprime
/ 2;
346 unsigned long i
, ii
, j
;
349 memset (s
, ~0, ssize
);
350 for (i
= 3; ; i
+= 2)
352 unsigned long isqr
= i
* i
;
353 if (isqr
>= maxprime
)
355 if (s
[i
* i
/ 2 - 1] == 0)
356 continue; /* only sieve with primes */
357 for (ii
= i
* i
/ 2 - 1; ii
< ssize
; ii
+= i
)
361 for (j
= 0; j
< ssize
; j
++)
365 primes
[n_primes
].prime
= j
* 2 + 3;
366 primes
[n_primes
].rem
= -1;
370 /* FIXME: This should not be needed if fencepost errors were fixed... */
371 if (primes
[n_primes
- 1].prime
> maxprime
)
377 for (i
= 3; i
<= maxprime
; i
+= 2)
379 if (i
< 7 || (i
% 3 != 0 && i
% 5 != 0 && i
% 7 != 0))
381 primes
[n_primes
].prime
= i
;
382 primes
[n_primes
].rem
= -1;