Annotation of OpenXM_contrib/gmp/demos/primes.c, Revision 1.1.1.2
1.1 maekawa 1: /* List and count primes.
1.1.1.2 ! ohara 2: Written by tege while on holiday in Rodupp, August 2001.
! 3: Between 10 and 500 times faster than previous program.
1.1 maekawa 4:
1.1.1.2 ! ohara 5: Copyright 2001 Free Software Foundation, Inc.
1.1 maekawa 6:
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
10: version.
11:
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.
15:
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. */
19:
20: #include <stdlib.h>
21: #include <stdio.h>
1.1.1.2 ! ohara 22: #include <string.h>
! 23: #include <math.h>
! 24: #include <assert.h>
! 25:
! 26: /* IDEAS:
! 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
! 41: this.
! 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.
! 52: */
! 53:
1.1 maekawa 54: #include "gmp.h"
55:
1.1.1.2 ! ohara 56: struct primes
! 57: {
! 58: unsigned int prime;
! 59: signed int rem;
! 60: };
! 61:
! 62: struct primes *primes;
! 63: unsigned long n_primes;
! 64:
! 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);
! 68:
! 69: int flag_print = 1;
! 70: int flag_count = 0;
! 71: int flag_maxgap = 0;
! 72: unsigned long maxgap = 0;
! 73: unsigned long total_primes = 0;
1.1 maekawa 74:
1.1.1.2 ! ohara 75: void
! 76: report (mpz_t prime)
1.1 maekawa 77: {
1.1.1.2 ! ohara 78: total_primes += 1;
! 79: if (flag_print)
! 80: {
! 81: mpz_out_str (stdout, 10, prime);
! 82: printf ("\n");
! 83: }
! 84: if (flag_maxgap)
! 85: {
! 86: static unsigned long prev_prime_low = 0;
! 87: unsigned long gap;
! 88: if (prev_prime_low != 0)
! 89: {
! 90: gap = mpz_get_ui (prime) - prev_prime_low;
! 91: if (maxgap < gap)
! 92: maxgap = gap;
! 93: }
! 94: prev_prime_low = mpz_get_ui (prime);
! 95: }
! 96: }
1.1 maekawa 97:
1.1.1.2 ! ohara 98: int
! 99: main (int argc, char *argv[])
! 100: {
! 101: char *progname = argv[0];
! 102: mpz_t fr, to;
! 103: mpz_t fr2, to2;
! 104: unsigned long sieve_lim;
! 105: unsigned long est_n_primes;
! 106: unsigned char *s;
! 107: mpz_t tmp;
! 108: mpz_t siev_sqr_lim;
1.1 maekawa 109:
110: while (argc != 1)
111: {
112: if (strcmp (argv[1], "-c") == 0)
113: {
114: flag_count = 1;
115: argv++;
116: argc--;
117: }
118: else if (strcmp (argv[1], "-p") == 0)
119: {
120: flag_print = 2;
121: argv++;
122: argc--;
123: }
1.1.1.2 ! ohara 124: else if (strcmp (argv[1], "-g") == 0)
! 125: {
! 126: flag_maxgap = 1;
! 127: argv++;
! 128: argc--;
! 129: }
1.1 maekawa 130: else
131: break;
132: }
133:
1.1.1.2 ! ohara 134: if (flag_count || flag_maxgap)
1.1 maekawa 135: flag_print--; /* clear unless an explicit -p */
136:
1.1.1.2 ! ohara 137: mpz_init (fr);
! 138: mpz_init (to);
! 139: mpz_init (fr2);
! 140: mpz_init (to2);
1.1 maekawa 141:
1.1.1.2 ! ohara 142: if (argc == 3)
1.1 maekawa 143: {
1.1.1.2 ! ohara 144: mpz_set_str (fr, argv[1], 0);
1.1 maekawa 145: if (argv[2][0] == '+')
146: {
1.1.1.2 ! ohara 147: mpz_set_str (to, argv[2] + 1, 0);
! 148: mpz_add (to, to, fr);
1.1 maekawa 149: }
150: else
1.1.1.2 ! ohara 151: mpz_set_str (to, argv[2], 0);
! 152: }
! 153: else if (argc == 2)
! 154: {
! 155: mpz_set_ui (fr, 0);
! 156: mpz_set_str (to, argv[1], 0);
1.1 maekawa 157: }
158: else
159: {
1.1.1.2 ! ohara 160: fprintf (stderr, "usage: %s [-c] [-p] [-g] [from [+]]to\n", progname);
1.1 maekawa 161: exit (1);
162: }
163:
1.1.1.2 ! ohara 164: mpz_set (fr2, fr);
! 165: if (mpz_cmp_ui (fr2, 3) < 0)
! 166: {
! 167: mpz_set_ui (fr2, 2);
! 168: report (fr2);
! 169: mpz_set_ui (fr2, 3);
! 170: }
! 171: mpz_setbit (fr2, 0); /* make odd */
! 172: mpz_sub_ui (to2, to, 1);
! 173: mpz_setbit (to2, 0); /* make odd */
! 174:
! 175: mpz_init (tmp);
! 176: mpz_init (siev_sqr_lim);
! 177:
! 178: mpz_sqrt (tmp, to2);
! 179: #define SIEVE_LIMIT 10000000
! 180: if (mpz_cmp_ui (tmp, SIEVE_LIMIT) < 0)
! 181: {
! 182: sieve_lim = mpz_get_ui (tmp);
! 183: }
! 184: else
! 185: {
! 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 */
! 190: }
! 191: mpz_set_ui (siev_sqr_lim, sieve_lim + 1);
! 192: mpz_mul_ui (siev_sqr_lim, siev_sqr_lim, sieve_lim + 1);
1.1 maekawa 193:
1.1.1.2 ! ohara 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);
! 198:
! 199: #if DEBUG
! 200: printf ("sieve_lim = %lu\n", sieve_lim);
! 201: printf ("n_primes = %lu (3..%u)\n",
! 202: n_primes, primes[n_primes - 1].prime);
! 203: #endif
! 204:
! 205: #define S (1 << 15) /* FIXME: Figure out L1 cache size */
! 206: s = malloc (S/2);
! 207: while (mpz_cmp (fr2, to2) <= 0)
1.1 maekawa 208: {
1.1.1.2 ! ohara 209: unsigned long rsize;
! 210: rsize = S;
! 211: mpz_add_ui (tmp, fr2, rsize);
! 212: if (mpz_cmp (tmp, to2) > 0)
1.1 maekawa 213: {
1.1.1.2 ! ohara 214: mpz_sub (tmp, to2, fr2);
! 215: rsize = mpz_get_ui (tmp) + 2;
1.1 maekawa 216: }
1.1.1.2 ! ohara 217: #if DEBUG
! 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");
! 221: #endif
! 222: sieve_region (s, fr2, rsize);
! 223: find_primes (s, fr2, rsize / 2, siev_sqr_lim);
! 224:
! 225: mpz_add_ui (fr2, fr2, S);
1.1 maekawa 226: }
1.1.1.2 ! ohara 227: free (s);
! 228:
! 229: if (flag_count)
! 230: printf ("Pi(interval) = %lu\n", total_primes);
! 231:
! 232: if (flag_maxgap)
! 233: printf ("max gap: %lu\n", maxgap);
1.1 maekawa 234:
1.1.1.2 ! ohara 235: return 0;
! 236: }
1.1 maekawa 237:
1.1.1.2 ! ohara 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". */
! 241: void
! 242: sieve_region (unsigned char *s, mpz_t fr, unsigned long rsize)
! 243: {
! 244: unsigned long ssize = rsize / 2;
! 245: unsigned long start, start2, prime;
! 246: unsigned long i;
! 247: mpz_t tmp;
! 248:
! 249: mpz_init (tmp);
! 250:
! 251: #if 0
! 252: /* initialize sieving array */
! 253: for (ii = 0; ii < (ssize + sizeof (long) - 1) / sizeof (long); ii++)
! 254: ((long *) s) [ii] = ~0L;
! 255: #else
! 256: {
! 257: signed long k;
! 258: long *se = (long *) (s + ((ssize + sizeof (long) - 1) & -sizeof (long)));
! 259: for (k = -((ssize + sizeof (long) - 1) / sizeof (long)); k < 0; k++)
! 260: se[k] = ~0L;
! 261: }
! 262: #endif
! 263:
! 264: for (i = 0; i < n_primes; i++)
1.1 maekawa 265: {
1.1.1.2 ! ohara 266: prime = primes[i].prime;
! 267:
! 268: if (primes[i].rem >= 0)
1.1 maekawa 269: {
1.1.1.2 ! ohara 270: start2 = primes[i].rem;
1.1 maekawa 271: }
1.1.1.2 ! ohara 272: else
! 273: {
! 274: mpz_set_ui (tmp, prime);
! 275: mpz_mul_ui (tmp, tmp, prime);
! 276: if (mpz_cmp (fr, tmp) <= 0)
! 277: {
! 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);
! 282: }
! 283: else
! 284: {
! 285: start = (prime - mpz_tdiv_ui (fr, prime)) % prime;
! 286: if (start % 2 != 0)
! 287: start += prime; /* adjust if even divisable */
! 288: }
! 289: start2 = start / 2;
! 290: }
! 291:
! 292: #if 0
! 293: for (ii = start2; ii < ssize; ii += prime)
! 294: s[ii] = 0;
! 295: primes[i].rem = ii - ssize;
! 296: #else
! 297: {
! 298: signed long k;
! 299: unsigned char *se = s + ssize; /* point just beyond sieving range */
! 300: for (k = start2 - ssize; k < 0; k += prime)
! 301: se[k] = 0;
! 302: primes[i].rem = k;
! 303: }
! 304: #endif
! 305: }
! 306: mpz_clear (tmp);
! 307: }
! 308:
! 309: /* Find primes in region [fr,fr+rsize), using the previously sieved s[]. */
! 310: void
! 311: find_primes (unsigned char *s, mpz_t fr, unsigned long ssize,
! 312: mpz_t siev_sqr_lim)
! 313: {
! 314: unsigned long j, ij;
! 315: mpz_t tmp;
! 316:
! 317: mpz_init (tmp);
! 318: for (j = 0; j < (ssize + sizeof (long) - 1) / sizeof (long); j++)
! 319: {
! 320: if (((long *) s) [j] != 0)
1.1 maekawa 321: {
1.1.1.2 ! ohara 322: for (ij = 0; ij < sizeof (long); ij++)
1.1 maekawa 323: {
1.1.1.2 ! ohara 324: if (s[j * sizeof (long) + ij] != 0)
1.1 maekawa 325: {
1.1.1.2 ! ohara 326: if (j * sizeof (long) + ij >= ssize)
! 327: goto out;
! 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))
! 331: report (tmp);
1.1 maekawa 332: }
333: }
334: }
335: }
1.1.1.2 ! ohara 336: out:
! 337: mpz_clear (tmp);
! 338: }
! 339:
! 340: /* Generate a lits of primes and store in the global array primes[]. */
! 341: void
! 342: make_primelist (unsigned long maxprime)
! 343: {
! 344: #if 1
! 345: unsigned char *s;
! 346: unsigned long ssize = maxprime / 2;
! 347: unsigned long i, ii, j;
! 348:
! 349: s = malloc (ssize);
! 350: memset (s, ~0, ssize);
! 351: for (i = 3; ; i += 2)
! 352: {
! 353: unsigned long isqr = i * i;
! 354: if (isqr >= maxprime)
! 355: break;
! 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)
! 359: s[ii] = 0;
! 360: }
! 361: n_primes = 0;
! 362: for (j = 0; j < ssize; j++)
! 363: {
! 364: if (s[j] != 0)
! 365: {
! 366: primes[n_primes].prime = j * 2 + 3;
! 367: primes[n_primes].rem = -1;
! 368: n_primes++;
! 369: }
! 370: }
! 371: /* FIXME: This should not be needed if fencepost errors were fixed... */
! 372: if (primes[n_primes - 1].prime > maxprime)
! 373: n_primes--;
! 374: free (s);
! 375: #else
! 376: unsigned long i;
! 377: n_primes = 0;
! 378: for (i = 3; i <= maxprime; i += 2)
! 379: {
! 380: if (i < 7 || (i % 3 != 0 && i % 5 != 0 && i % 7 != 0))
! 381: {
! 382: primes[n_primes].prime = i;
! 383: primes[n_primes].rem = -1;
! 384: n_primes++;
! 385: }
! 386: }
! 387: #endif
1.1 maekawa 388: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to primes.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / demos