]> git.cameronkatri.com Git - bsdgames-darwin.git/blob - primes/primes.c
merge with Lite, new RCS id conventions, etc.
[bsdgames-darwin.git] / primes / primes.c
1 /* $NetBSD: primes.c,v 1.4 1995/03/23 08:35:55 cgd Exp $ */
2
3 /*
4 * Copyright (c) 1989, 1993
5 * The Regents of the University of California. All rights reserved.
6 *
7 * This code is derived from software contributed to Berkeley by
8 * Landon Curt Noll.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in the
17 * documentation and/or other materials provided with the distribution.
18 * 3. All advertising materials mentioning features or use of this software
19 * must display the following acknowledgement:
20 * This product includes software developed by the University of
21 * California, Berkeley and its contributors.
22 * 4. Neither the name of the University nor the names of its contributors
23 * may be used to endorse or promote products derived from this software
24 * without specific prior written permission.
25 *
26 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
27 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
29 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
30 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
31 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
32 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
33 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
34 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
35 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
36 * SUCH DAMAGE.
37 */
38
39 #ifndef lint
40 static char copyright[] =
41 "@(#) Copyright (c) 1989, 1993\n\
42 The Regents of the University of California. All rights reserved.\n";
43 #endif /* not lint */
44
45 #ifndef lint
46 #if 0
47 static char sccsid[] = "@(#)primes.c 8.4 (Berkeley) 3/21/94";
48 #else
49 static char rcsid[] = "$NetBSD: primes.c,v 1.4 1995/03/23 08:35:55 cgd Exp $";
50 #endif
51 #endif /* not lint */
52
53 /*
54 * primes - generate a table of primes between two values
55 *
56 * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo
57 *
58 * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
59 *
60 * usage:
61 * primes [start [stop]]
62 *
63 * Print primes >= start and < stop. If stop is omitted,
64 * the value 4294967295 (2^32-1) is assumed. If start is
65 * omitted, start is read from standard input.
66 *
67 * validation check: there are 664579 primes between 0 and 10^7
68 */
69
70 #include <ctype.h>
71 #include <err.h>
72 #include <errno.h>
73 #include <limits.h>
74 #include <math.h>
75 #include <memory.h>
76 #include <stdio.h>
77 #include <stdlib.h>
78
79 #include "primes.h"
80
81 /*
82 * Eratosthenes sieve table
83 *
84 * We only sieve the odd numbers. The base of our sieve windows are always
85 * odd. If the base of table is 1, table[i] represents 2*i-1. After the
86 * sieve, table[i] == 1 if and only iff 2*i-1 is prime.
87 *
88 * We make TABSIZE large to reduce the overhead of inner loop setup.
89 */
90 char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */
91
92 /*
93 * prime[i] is the (i-1)th prime.
94 *
95 * We are able to sieve 2^32-1 because this byte table yields all primes
96 * up to 65537 and 65537^2 > 2^32-1.
97 */
98 extern ubig prime[];
99 extern ubig *pr_limit; /* largest prime in the prime array */
100
101 /*
102 * To avoid excessive sieves for small factors, we use the table below to
103 * setup our sieve blocks. Each element represents a odd number starting
104 * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13.
105 */
106 extern char pattern[];
107 extern int pattern_size; /* length of pattern array */
108
109 void primes __P((ubig, ubig));
110 ubig read_num_buf __P((void));
111 void usage __P((void));
112
113 int
114 main(argc, argv)
115 int argc;
116 char *argv[];
117 {
118 ubig start; /* where to start generating */
119 ubig stop; /* don't generate at or above this value */
120 int ch;
121 char *p;
122
123 while ((ch = getopt(argc, argv, "")) != EOF)
124 switch (ch) {
125 case '?':
126 default:
127 usage();
128 }
129 argc -= optind;
130 argv += optind;
131
132 start = 0;
133 stop = BIG;
134
135 /*
136 * Convert low and high args. Strtoul(3) sets errno to
137 * ERANGE if the number is too large, but, if there's
138 * a leading minus sign it returns the negation of the
139 * result of the conversion, which we'd rather disallow.
140 */
141 switch (argc) {
142 case 2:
143 /* Start and stop supplied on the command line. */
144 if (argv[0][0] == '-' || argv[1][0] == '-')
145 errx(1, "negative numbers aren't permitted.");
146
147 errno = 0;
148 start = strtoul(argv[0], &p, 10);
149 if (errno)
150 err(1, "%s", argv[0]);
151 if (*p != '\0')
152 errx(1, "%s: illegal numeric format.", argv[0]);
153
154 errno = 0;
155 stop = strtoul(argv[1], &p, 10);
156 if (errno)
157 err(1, "%s", argv[1]);
158 if (*p != '\0')
159 errx(1, "%s: illegal numeric format.", argv[1]);
160 break;
161 case 1:
162 /* Start on the command line. */
163 if (argv[0][0] == '-')
164 errx(1, "negative numbers aren't permitted.");
165
166 errno = 0;
167 start = strtoul(argv[0], &p, 10);
168 if (errno)
169 err(1, "%s", argv[0]);
170 if (*p != '\0')
171 errx(1, "%s: illegal numeric format.", argv[0]);
172 break;
173 case 0:
174 start = read_num_buf();
175 break;
176 default:
177 usage();
178 }
179
180 if (start > stop)
181 errx(1, "start value must be less than stop value.");
182 primes(start, stop);
183 exit(0);
184 }
185
186 /*
187 * read_num_buf --
188 * This routine returns a number n, where 0 <= n && n <= BIG.
189 */
190 ubig
191 read_num_buf()
192 {
193 ubig val;
194 char *p, buf[100]; /* > max number of digits. */
195
196 for (;;) {
197 if (fgets(buf, sizeof(buf), stdin) == NULL) {
198 if (ferror(stdin))
199 err(1, "stdin");
200 exit(0);
201 }
202 for (p = buf; isblank(*p); ++p);
203 if (*p == '\n' || *p == '\0')
204 continue;
205 if (*p == '-')
206 errx(1, "negative numbers aren't permitted.");
207 errno = 0;
208 val = strtoul(buf, &p, 10);
209 if (errno)
210 err(1, "%s", buf);
211 if (*p != '\n')
212 errx(1, "%s: illegal numeric format.", buf);
213 return (val);
214 }
215 }
216
217 /*
218 * primes - sieve and print primes from start up to and but not including stop
219 */
220 void
221 primes(start, stop)
222 ubig start; /* where to start generating */
223 ubig stop; /* don't generate at or above this value */
224 {
225 register char *q; /* sieve spot */
226 register ubig factor; /* index and factor */
227 register char *tab_lim; /* the limit to sieve on the table */
228 register ubig *p; /* prime table pointer */
229 register ubig fact_lim; /* highest prime for current block */
230
231 /*
232 * A number of systems can not convert double values into unsigned
233 * longs when the values are larger than the largest signed value.
234 * We don't have this problem, so we can go all the way to BIG.
235 */
236 if (start < 3) {
237 start = (ubig)2;
238 }
239 if (stop < 3) {
240 stop = (ubig)2;
241 }
242 if (stop <= start) {
243 return;
244 }
245
246 /*
247 * be sure that the values are odd, or 2
248 */
249 if (start != 2 && (start&0x1) == 0) {
250 ++start;
251 }
252 if (stop != 2 && (stop&0x1) == 0) {
253 ++stop;
254 }
255
256 /*
257 * quick list of primes <= pr_limit
258 */
259 if (start <= *pr_limit) {
260 /* skip primes up to the start value */
261 for (p = &prime[0], factor = prime[0];
262 factor < stop && p <= pr_limit; factor = *(++p)) {
263 if (factor >= start) {
264 printf("%u\n", factor);
265 }
266 }
267 /* return early if we are done */
268 if (p <= pr_limit) {
269 return;
270 }
271 start = *pr_limit+2;
272 }
273
274 /*
275 * we shall sieve a bytemap window, note primes and move the window
276 * upward until we pass the stop point
277 */
278 while (start < stop) {
279 /*
280 * factor out 3, 5, 7, 11 and 13
281 */
282 /* initial pattern copy */
283 factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
284 memcpy(table, &pattern[factor], pattern_size-factor);
285 /* main block pattern copies */
286 for (fact_lim=pattern_size-factor;
287 fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) {
288 memcpy(&table[fact_lim], pattern, pattern_size);
289 }
290 /* final block pattern copy */
291 memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
292
293 /*
294 * sieve for primes 17 and higher
295 */
296 /* note highest useful factor and sieve spot */
297 if (stop-start > TABSIZE+TABSIZE) {
298 tab_lim = &table[TABSIZE]; /* sieve it all */
299 fact_lim = (int)sqrt(
300 (double)(start)+TABSIZE+TABSIZE+1.0);
301 } else {
302 tab_lim = &table[(stop-start)/2]; /* partial sieve */
303 fact_lim = (int)sqrt((double)(stop)+1.0);
304 }
305 /* sieve for factors >= 17 */
306 factor = 17; /* 17 is first prime to use */
307 p = &prime[7]; /* 19 is next prime, pi(19)=7 */
308 do {
309 /* determine the factor's initial sieve point */
310 q = (char *)(start%factor); /* temp storage for mod */
311 if ((int)q & 0x1) {
312 q = &table[(factor-(int)q)/2];
313 } else {
314 q = &table[q ? factor-((int)q/2) : 0];
315 }
316 /* sive for our current factor */
317 for ( ; q < tab_lim; q += factor) {
318 *q = '\0'; /* sieve out a spot */
319 }
320 } while ((factor=(ubig)(*(p++))) <= fact_lim);
321
322 /*
323 * print generated primes
324 */
325 for (q = table; q < tab_lim; ++q, start+=2) {
326 if (*q) {
327 printf("%u\n", start);
328 }
329 }
330 }
331 }
332
333 void
334 usage()
335 {
336 (void)fprintf(stderr, "usage: primes [start [stop]]\n");
337 exit(1);
338 }