=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/demos/Attic/primes.c,v retrieving revision 1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/demos/Attic/primes.c 2000/09/09 14:13:19 1.1 +++ OpenXM_contrib/gmp/demos/Attic/primes.c 2003/08/25 16:06:03 1.1.1.2 @@ -1,6 +1,8 @@ /* List and count primes. + Written by tege while on holiday in Rodupp, August 2001. + Between 10 and 500 times faster than previous program. -Copyright (C) 2000 Free Software Foundation, Inc. +Copyright 2001 Free Software Foundation, Inc. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software @@ -17,27 +19,94 @@ Place - Suite 330, Boston, MA 02111-1307, USA. */ #include #include +#include +#include +#include + +/* IDEAS: + * Do not fill primes[] with real primes when the range [fr,to] is small, + when fr,to are relatively large. Fill primes[] with odd numbers instead. + [Probably a bad idea, since the primes[] array would become very large.] + * Separate small primes and large primes when sieving. Either the Montgomery + way (i.e., having a large array a multiple of L1 cache size), or just + separate loops for primes <= S and primes > S. The latter primes do not + require an inner loop, since they will touch the sieving array at most once. + * Pre-fill sieving array with an appropriately aligned ...00100100... pattern, + then omit 3 from primes array. (May require similar special handling of 3 + as we now have for 2.) + * A large SIEVE_LIMIT currently implies very large memory usage, mainly due + to the sieving array in make_primelist, but also because of the primes[] + array. We might want to stage the program, using sieve_region/find_primes + to build primes[]. Make report() a function pointer, as part of achieving + this. + * Store primes[] as two arrays, one array with primes represented as delta + values using just 8 bits (if gaps are too big, store bogus primes!) + and one array with "rem" values. The latter needs 32-bit values. + * A new entry point, mpz_probab_prime_likely_p, would be useful. + * Improve command line syntax and versatility. "primes -f FROM -t TO", + allow either to be omitted for open interval. (But disallow + "primes -c -f FROM" since that would be infinity.) Allow printing a + limited *number* of primes using syntax like "primes -f FROM -n NUMBER". + * When looking for maxgaps, we should not perform any primality testing until + we find possible record gaps. Should speed up the searches tremendously. + */ + #include "gmp.h" -/* Sieve table size */ -#define ST_SIZE 30000 -/* Largest prime to sieve with */ -#define MAX_S_PRIME 1000 +struct primes +{ + unsigned int prime; + signed int rem; +}; -main (int argc, char **argv) +struct primes *primes; +unsigned long n_primes; + +void find_primes (unsigned char *, mpz_t, unsigned long, mpz_t); +void sieve_region (unsigned char *, mpz_t, unsigned long); +void make_primelist (unsigned long); + +int flag_print = 1; +int flag_count = 0; +int flag_maxgap = 0; +unsigned long maxgap = 0; +unsigned long total_primes = 0; + +void +report (mpz_t prime) { + total_primes += 1; + if (flag_print) + { + mpz_out_str (stdout, 10, prime); + printf ("\n"); + } + if (flag_maxgap) + { + static unsigned long prev_prime_low = 0; + unsigned long gap; + if (prev_prime_low != 0) + { + gap = mpz_get_ui (prime) - prev_prime_low; + if (maxgap < gap) + maxgap = gap; + } + prev_prime_low = mpz_get_ui (prime); + } +} + +int +main (int argc, char *argv[]) +{ char *progname = argv[0]; - mpz_t r0, r1; /* range */ - mpz_t cur; - unsigned char *st; - unsigned long i, ii; - unsigned long nprimes = 0; - unsigned long last; - int flag_print = 1; - int flag_count = 0; + mpz_t fr, to; + mpz_t fr2, to2; + unsigned long sieve_lim; + unsigned long est_n_primes; + unsigned char *s; + mpz_t tmp; + mpz_t siev_sqr_lim; - st = malloc (ST_SIZE); - while (argc != 1) { if (strcmp (argv[1], "-c") == 0) @@ -52,94 +121,268 @@ main (int argc, char **argv) argv++; argc--; } + else if (strcmp (argv[1], "-g") == 0) + { + flag_maxgap = 1; + argv++; + argc--; + } else break; } - if (flag_count) + if (flag_count || flag_maxgap) flag_print--; /* clear unless an explicit -p */ - mpz_init (r0); - mpz_init (r1); - mpz_init (cur); + mpz_init (fr); + mpz_init (to); + mpz_init (fr2); + mpz_init (to2); - if (argc == 2) + if (argc == 3) { - mpz_set_ui (r0, 2); - mpz_set_str (r1, argv[1], 0); - } - else if (argc == 3) - { - mpz_set_str (r0, argv[1], 0); + mpz_set_str (fr, argv[1], 0); if (argv[2][0] == '+') { - mpz_set_str (r1, argv[2] + 1, 0); - mpz_add (r1, r1, r0); + mpz_set_str (to, argv[2] + 1, 0); + mpz_add (to, to, fr); } else - mpz_set_str (r1, argv[2], 0); + mpz_set_str (to, argv[2], 0); } + else if (argc == 2) + { + mpz_set_ui (fr, 0); + mpz_set_str (to, argv[1], 0); + } else { - fprintf (stderr, "usage: %s [-c] [-g] [from [+]]to\n", progname); + fprintf (stderr, "usage: %s [-c] [-p] [-g] [from [+]]to\n", progname); exit (1); } - if (mpz_cmp_ui (r0, 2) < 0) - mpz_set_ui (r0, 2); + mpz_set (fr2, fr); + if (mpz_cmp_ui (fr2, 3) < 0) + { + mpz_set_ui (fr2, 2); + report (fr2); + mpz_set_ui (fr2, 3); + } + mpz_setbit (fr2, 0); /* make odd */ + mpz_sub_ui (to2, to, 1); + mpz_setbit (to2, 0); /* make odd */ - if ((mpz_get_ui (r0) & 1) == 0) + mpz_init (tmp); + mpz_init (siev_sqr_lim); + + mpz_sqrt (tmp, to2); +#define SIEVE_LIMIT 10000000 + if (mpz_cmp_ui (tmp, SIEVE_LIMIT) < 0) { - if (mpz_cmp_ui (r0, 2) == 0) + sieve_lim = mpz_get_ui (tmp); + } + else + { + sieve_lim = SIEVE_LIMIT; + mpz_sub (tmp, to2, fr2); + if (mpz_cmp_ui (tmp, sieve_lim) < 0) + sieve_lim = mpz_get_ui (tmp); /* limit sieving for small ranges */ + } + mpz_set_ui (siev_sqr_lim, sieve_lim + 1); + mpz_mul_ui (siev_sqr_lim, siev_sqr_lim, sieve_lim + 1); + + est_n_primes = (size_t) (sieve_lim / log((double) sieve_lim) * 1.13) + 10; + primes = malloc (est_n_primes * sizeof primes[0]); + make_primelist (sieve_lim); + assert (est_n_primes >= n_primes); + +#if DEBUG + printf ("sieve_lim = %lu\n", sieve_lim); + printf ("n_primes = %lu (3..%u)\n", + n_primes, primes[n_primes - 1].prime); +#endif + +#define S (1 << 15) /* FIXME: Figure out L1 cache size */ + s = malloc (S/2); + while (mpz_cmp (fr2, to2) <= 0) + { + unsigned long rsize; + rsize = S; + mpz_add_ui (tmp, fr2, rsize); + if (mpz_cmp (tmp, to2) > 0) { - if (flag_print) - puts ("2"); - nprimes++; + mpz_sub (tmp, to2, fr2); + rsize = mpz_get_ui (tmp) + 2; } - mpz_add_ui (r0, r0, 1); +#if DEBUG + printf ("Sieving region ["); mpz_out_str (stdout, 10, fr2); + printf (","); mpz_add_ui (tmp, fr2, rsize - 2); + mpz_out_str (stdout, 10, tmp); printf ("]\n"); +#endif + sieve_region (s, fr2, rsize); + find_primes (s, fr2, rsize / 2, siev_sqr_lim); + + mpz_add_ui (fr2, fr2, S); } + free (s); - mpz_set (cur, r0); + if (flag_count) + printf ("Pi(interval) = %lu\n", total_primes); - while (mpz_cmp (cur, r1) <= 0) + if (flag_maxgap) + printf ("max gap: %lu\n", maxgap); + + return 0; +} + +/* Find primes in region [fr,fr+rsize). Requires that fr is odd and that + rsize is even. The sieving array s should be aligned for "long int" and + have rsize/2 entries, rounded up to the nearest multiple of "long int". */ +void +sieve_region (unsigned char *s, mpz_t fr, unsigned long rsize) +{ + unsigned long ssize = rsize / 2; + unsigned long start, start2, prime; + unsigned long i; + mpz_t tmp; + + mpz_init (tmp); + +#if 0 + /* initialize sieving array */ + for (ii = 0; ii < (ssize + sizeof (long) - 1) / sizeof (long); ii++) + ((long *) s) [ii] = ~0L; +#else + { + signed long k; + long *se = (long *) (s + ((ssize + sizeof (long) - 1) & -sizeof (long))); + for (k = -((ssize + sizeof (long) - 1) / sizeof (long)); k < 0; k++) + se[k] = ~0L; + } +#endif + + for (i = 0; i < n_primes; i++) { - memset (st, 1, ST_SIZE); - for (i = 3; i < MAX_S_PRIME; i += 2) + prime = primes[i].prime; + + if (primes[i].rem >= 0) { - unsigned long start; - start = mpz_tdiv_ui (cur, i); - if (start != 0) - start = i - start; - - if (mpz_cmp_ui (cur, i - start) == 0) - start += i; - for (ii = start; ii < ST_SIZE; ii += i) - st[ii] = 0; + start2 = primes[i].rem; } - last = 0; - for (ii = 0; ii < ST_SIZE; ii += 2) + else { - if (st[ii] != 0) + mpz_set_ui (tmp, prime); + mpz_mul_ui (tmp, tmp, prime); + if (mpz_cmp (fr, tmp) <= 0) { - mpz_add_ui (cur, cur, ii - last); - last = ii; - if (mpz_cmp (cur, r1) > 0) - goto done; - if (mpz_probab_prime_p (cur, 3)) + mpz_sub (tmp, tmp, fr); + if (mpz_cmp_ui (tmp, 2 * ssize) > 0) + break; /* avoid overflow at next line, also speedup */ + start = mpz_get_ui (tmp); + } + else + { + start = (prime - mpz_tdiv_ui (fr, prime)) % prime; + if (start % 2 != 0) + start += prime; /* adjust if even divisable */ + } + start2 = start / 2; + } + +#if 0 + for (ii = start2; ii < ssize; ii += prime) + s[ii] = 0; + primes[i].rem = ii - ssize; +#else + { + signed long k; + unsigned char *se = s + ssize; /* point just beyond sieving range */ + for (k = start2 - ssize; k < 0; k += prime) + se[k] = 0; + primes[i].rem = k; + } +#endif + } + mpz_clear (tmp); +} + +/* Find primes in region [fr,fr+rsize), using the previously sieved s[]. */ +void +find_primes (unsigned char *s, mpz_t fr, unsigned long ssize, + mpz_t siev_sqr_lim) +{ + unsigned long j, ij; + mpz_t tmp; + + mpz_init (tmp); + for (j = 0; j < (ssize + sizeof (long) - 1) / sizeof (long); j++) + { + if (((long *) s) [j] != 0) + { + for (ij = 0; ij < sizeof (long); ij++) + { + if (s[j * sizeof (long) + ij] != 0) { - nprimes++; - if (flag_print) - { - mpz_out_str (stdout, 10, cur); puts (""); - } + if (j * sizeof (long) + ij >= ssize) + goto out; + mpz_add_ui (tmp, fr, (j * sizeof (long) + ij) * 2); + if (mpz_cmp (tmp, siev_sqr_lim) < 0 || + mpz_probab_prime_p (tmp, 3)) + report (tmp); } } } - mpz_add_ui (cur, cur, ST_SIZE - last); } - done: - if (flag_count) - printf ("Pi(interval) = %lu\n", nprimes); + out: + mpz_clear (tmp); +} - exit (0); +/* Generate a lits of primes and store in the global array primes[]. */ +void +make_primelist (unsigned long maxprime) +{ +#if 1 + unsigned char *s; + unsigned long ssize = maxprime / 2; + unsigned long i, ii, j; + + s = malloc (ssize); + memset (s, ~0, ssize); + for (i = 3; ; i += 2) + { + unsigned long isqr = i * i; + if (isqr >= maxprime) + break; + if (s[i * i / 2 - 1] == 0) + continue; /* only sieve with primes */ + for (ii = i * i / 2 - 1; ii < ssize; ii += i) + s[ii] = 0; + } + n_primes = 0; + for (j = 0; j < ssize; j++) + { + if (s[j] != 0) + { + primes[n_primes].prime = j * 2 + 3; + primes[n_primes].rem = -1; + n_primes++; + } + } + /* FIXME: This should not be needed if fencepost errors were fixed... */ + if (primes[n_primes - 1].prime > maxprime) + n_primes--; + free (s); +#else + unsigned long i; + n_primes = 0; + for (i = 3; i <= maxprime; i += 2) + { + if (i < 7 || (i % 3 != 0 && i % 5 != 0 && i % 7 != 0)) + { + primes[n_primes].prime = i; + primes[n_primes].rem = -1; + n_primes++; + } + } +#endif }