Annotation of OpenXM_contrib/gmp/demos/primes.c, Revision 1.1.1.1
1.1 maekawa 1: /* List and count primes.
2:
3: Copyright (C) 2000 Free Software Foundation, Inc.
4:
5: This program is free software; you can redistribute it and/or modify it under
6: the terms of the GNU General Public License as published by the Free Software
7: Foundation; either version 2 of the License, or (at your option) any later
8: version.
9:
10: This program is distributed in the hope that it will be useful, but WITHOUT ANY
11: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
12: PARTICULAR PURPOSE. See the GNU General Public License for more details.
13:
14: You should have received a copy of the GNU General Public License along with
15: this program; if not, write to the Free Software Foundation, Inc., 59 Temple
16: Place - Suite 330, Boston, MA 02111-1307, USA. */
17:
18: #include <stdlib.h>
19: #include <stdio.h>
20: #include "gmp.h"
21:
22: /* Sieve table size */
23: #define ST_SIZE 30000
24: /* Largest prime to sieve with */
25: #define MAX_S_PRIME 1000
26:
27: main (int argc, char **argv)
28: {
29: char *progname = argv[0];
30: mpz_t r0, r1; /* range */
31: mpz_t cur;
32: unsigned char *st;
33: unsigned long i, ii;
34: unsigned long nprimes = 0;
35: unsigned long last;
36: int flag_print = 1;
37: int flag_count = 0;
38:
39: st = malloc (ST_SIZE);
40:
41: while (argc != 1)
42: {
43: if (strcmp (argv[1], "-c") == 0)
44: {
45: flag_count = 1;
46: argv++;
47: argc--;
48: }
49: else if (strcmp (argv[1], "-p") == 0)
50: {
51: flag_print = 2;
52: argv++;
53: argc--;
54: }
55: else
56: break;
57: }
58:
59: if (flag_count)
60: flag_print--; /* clear unless an explicit -p */
61:
62: mpz_init (r0);
63: mpz_init (r1);
64: mpz_init (cur);
65:
66: if (argc == 2)
67: {
68: mpz_set_ui (r0, 2);
69: mpz_set_str (r1, argv[1], 0);
70: }
71: else if (argc == 3)
72: {
73: mpz_set_str (r0, argv[1], 0);
74: if (argv[2][0] == '+')
75: {
76: mpz_set_str (r1, argv[2] + 1, 0);
77: mpz_add (r1, r1, r0);
78: }
79: else
80: mpz_set_str (r1, argv[2], 0);
81: }
82: else
83: {
84: fprintf (stderr, "usage: %s [-c] [-g] [from [+]]to\n", progname);
85: exit (1);
86: }
87:
88: if (mpz_cmp_ui (r0, 2) < 0)
89: mpz_set_ui (r0, 2);
90:
91: if ((mpz_get_ui (r0) & 1) == 0)
92: {
93: if (mpz_cmp_ui (r0, 2) == 0)
94: {
95: if (flag_print)
96: puts ("2");
97: nprimes++;
98: }
99: mpz_add_ui (r0, r0, 1);
100: }
101:
102: mpz_set (cur, r0);
103:
104: while (mpz_cmp (cur, r1) <= 0)
105: {
106: memset (st, 1, ST_SIZE);
107: for (i = 3; i < MAX_S_PRIME; i += 2)
108: {
109: unsigned long start;
110: start = mpz_tdiv_ui (cur, i);
111: if (start != 0)
112: start = i - start;
113:
114: if (mpz_cmp_ui (cur, i - start) == 0)
115: start += i;
116: for (ii = start; ii < ST_SIZE; ii += i)
117: st[ii] = 0;
118: }
119: last = 0;
120: for (ii = 0; ii < ST_SIZE; ii += 2)
121: {
122: if (st[ii] != 0)
123: {
124: mpz_add_ui (cur, cur, ii - last);
125: last = ii;
126: if (mpz_cmp (cur, r1) > 0)
127: goto done;
128: if (mpz_probab_prime_p (cur, 3))
129: {
130: nprimes++;
131: if (flag_print)
132: {
133: mpz_out_str (stdout, 10, cur); puts ("");
134: }
135: }
136: }
137: }
138: mpz_add_ui (cur, cur, ST_SIZE - last);
139: }
140: done:
141: if (flag_count)
142: printf ("Pi(interval) = %lu\n", nprimes);
143:
144: exit (0);
145: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>