Annotation of src/games/primes/primes.c, Revision 1.3
1.1 cgd 1: /*
2: * Copyright (c) 1989 The Regents of the University of California.
3: * All rights reserved.
4: *
5: * This code is derived from software contributed to Berkeley by
6: * Landon Curt Noll.
7: *
8: * Redistribution and use in source and binary forms, with or without
9: * modification, are permitted provided that the following conditions
10: * are met:
11: * 1. Redistributions of source code must retain the above copyright
12: * notice, this list of conditions and the following disclaimer.
13: * 2. Redistributions in binary form must reproduce the above copyright
14: * notice, this list of conditions and the following disclaimer in the
15: * documentation and/or other materials provided with the distribution.
16: * 3. All advertising materials mentioning features or use of this software
17: * must display the following acknowledgement:
18: * This product includes software developed by the University of
19: * California, Berkeley and its contributors.
20: * 4. Neither the name of the University nor the names of its contributors
21: * may be used to endorse or promote products derived from this software
22: * without specific prior written permission.
23: *
24: * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
25: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
26: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
27: * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
28: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
29: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
30: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
31: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
33: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
34: * SUCH DAMAGE.
35: */
36:
37: #ifndef lint
38: char copyright[] =
39: "@(#) Copyright (c) 1989 The Regents of the University of California.\n\
40: All rights reserved.\n";
41: #endif /* not lint */
42:
43: #ifndef lint
1.2 mycroft 44: /*static char sccsid[] = "from: @(#)primes.c 5.4 (Berkeley) 6/1/90";*/
1.3 ! cgd 45: static char rcsid[] = "$Id: primes.c,v 1.2 1993/08/01 18:53:04 mycroft Exp $";
1.1 cgd 46: #endif /* not lint */
47:
48: /*
49: * primes - generate a table of primes between two values
50: *
51: * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo
52: *
53: * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
54: *
55: * usage:
56: * primes [start [stop]]
57: *
58: * Print primes >= start and < stop. If stop is omitted,
59: * the value 4294967295 (2^32-1) is assumed. If start is
60: * omitted, start is read from standard input.
61: *
62: * Prints "ouch" if start or stop is bogus.
63: *
64: * validation check: there are 664579 primes between 0 and 10^7
65: */
66:
67: #include <stdio.h>
68: #include <math.h>
69: #include <memory.h>
70: #include <ctype.h>
1.3 ! cgd 71: #include <limits.h>
1.1 cgd 72: #include "primes.h"
73:
74: /*
75: * Eratosthenes sieve table
76: *
77: * We only sieve the odd numbers. The base of our sieve windows are always
78: * odd. If the base of table is 1, table[i] represents 2*i-1. After the
79: * sieve, table[i] == 1 if and only iff 2*i-1 is prime.
80: *
81: * We make TABSIZE large to reduce the overhead of inner loop setup.
82: */
83: char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */
84:
85: /*
86: * prime[i] is the (i-1)th prime.
87: *
88: * We are able to sieve 2^32-1 because this byte table yields all primes
89: * up to 65537 and 65537^2 > 2^32-1.
90: */
91: extern ubig prime[];
92: extern ubig *pr_limit; /* largest prime in the prime array */
93:
94: /*
95: * To avoid excessive sieves for small factors, we use the table below to
96: * setup our sieve blocks. Each element represents a odd number starting
97: * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13.
98: */
99: extern char pattern[];
100: extern int pattern_size; /* length of pattern array */
101:
102: #define MAX_LINE 255 /* max line allowed on stdin */
103:
104: char *read_num_buf(); /* read a number buffer */
105: void primes(); /* print the primes in range */
106: char *program; /* our name */
107:
108: main(argc, argv)
109: int argc; /* arg count */
110: char *argv[]; /* args */
111: {
112: char buf[MAX_LINE+1]; /* input buffer */
113: char *ret; /* return result */
114: ubig start; /* where to start generating */
115: ubig stop; /* don't generate at or above this value */
116:
117: /*
118: * parse args
119: */
120: program = argv[0];
121: start = 0;
122: stop = BIG;
123: if (argc == 3) {
124: /* convert low and high args */
125: if (read_num_buf(NULL, argv[1]) == NULL) {
126: fprintf(stderr, "%s: ouch\n", program);
127: exit(1);
128: }
129: if (read_num_buf(NULL, argv[2]) == NULL) {
130: fprintf(stderr, "%s: ouch\n", program);
131: exit(1);
132: }
1.3 ! cgd 133: if (sscanf(argv[1], "%lu", &start) != 1) {
1.1 cgd 134: fprintf(stderr, "%s: ouch\n", program);
135: exit(1);
136: }
1.3 ! cgd 137: if (sscanf(argv[2], "%lu", &stop) != 1) {
1.1 cgd 138: fprintf(stderr, "%s: ouch\n", program);
139: exit(1);
140: }
141:
142: } else if (argc == 2) {
143: /* convert low arg */
144: if (read_num_buf(NULL, argv[1]) == NULL) {
145: fprintf(stderr, "%s: ouch\n", program);
146: exit(1);
147: }
1.3 ! cgd 148: if (sscanf(argv[1], "%lu", &start) != 1) {
1.1 cgd 149: fprintf(stderr, "%s: ouch\n", program);
150: exit(1);
151: }
152:
153: } else {
154: /* read input until we get a good line */
155: if (read_num_buf(stdin, buf) != NULL) {
156:
157: /* convert the buffer */
1.3 ! cgd 158: if (sscanf(buf, "%lu", &start) != 1) {
1.1 cgd 159: fprintf(stderr, "%s: ouch\n", program);
160: exit(1);
161: }
162: } else {
163: exit(0);
164: }
165: }
166: if (start > stop) {
167: fprintf(stderr, "%s: ouch\n", program);
168: exit(1);
169: }
170: primes(start, stop);
171: exit(0);
172: }
173:
174: /*
175: * read_num_buf - read a number buffer from a stream
176: *
177: * Read a number on a line of the form:
178: *
179: * ^[ \t]*\(+?[0-9][0-9]\)*.*$
180: *
181: * where ? is a 1-or-0 operator and the number is within \( \).
182: *
183: * If does not match the above pattern, it is ignored and a new
184: * line is read. If the number is too large or small, we will
185: * print ouch and read a new line.
186: *
187: * We have to be very careful on how we check the magnitude of the
188: * input. We can not use numeric checks because of the need to
189: * check values against maximum numeric values.
190: *
191: * This routine will return a line containing a ascii number between
192: * 0 and BIG, or it will return NULL.
193: *
194: * If the stream is NULL then buf will be processed as if were
195: * a single line stream.
196: *
197: * returns:
198: * char * pointer to leading digit or +
199: * NULL EOF or error
200: */
201: char *
202: read_num_buf(input, buf)
203: FILE *input; /* input stream or NULL */
204: char *buf; /* input buffer */
205: {
206: static char limit[MAX_LINE+1]; /* ascii value of BIG */
207: static int limit_len; /* digit count of limit */
208: int len; /* digits in input (excluding +/-) */
209: char *s; /* line start marker */
210: char *d; /* first digit, skip +/- */
211: char *p; /* scan pointer */
212: char *z; /* zero scan pointer */
213:
1.3 ! cgd 214: /* form the ascii value of BIG if needed */
1.1 cgd 215: if (!isascii(limit[0]) || !isdigit(limit[0])) {
1.3 ! cgd 216: sprintf(limit, "%lu", BIG);
1.1 cgd 217: limit_len = strlen(limit);
218: }
219:
220: /*
221: * the search for a good line
222: */
223: if (input != NULL && fgets(buf, MAX_LINE, input) == NULL) {
224: /* error or EOF */
225: return NULL;
226: }
227: do {
228:
229: /* ignore leading whitespace */
230: for (s=buf; *s && s < buf+MAX_LINE; ++s) {
231: if (!isascii(*s) || !isspace(*s)) {
232: break;
233: }
234: }
235:
236: /* object if - */
237: if (*s == '-') {
1.3 ! cgd 238: fprintf(stderr, "%s: ouch for minuses\n", program);
1.1 cgd 239: continue;
240: }
241:
242: /* skip over any leading + */
243: if (*s == '+') {
244: d = s+1;
245: } else {
246: d = s;
247: }
248:
249: /* note leading zeros */
250: for (z=d; *z && z < buf+MAX_LINE; ++z) {
251: if (*z != '0') {
252: break;
253: }
254: }
255:
256: /* scan for the first non-digit/non-plus/non-minus */
257: for (p=d; *p && p < buf+MAX_LINE; ++p) {
258: if (!isascii(*p) || !isdigit(*p)) {
259: break;
260: }
261: }
262:
263: /* ignore empty lines */
264: if (p == d) {
265: continue;
266: }
267: *p = '\0';
268:
269: /* object if too many digits */
270: len = strlen(z);
271: len = (len<=0) ? 1 : len;
272: /* accept if digit count is below limit */
273: if (len < limit_len) {
274: /* we have good input */
275: return s;
276:
277: /* reject very large numbers */
278: } else if (len > limit_len) {
1.3 ! cgd 279: fprintf(stderr, "%s: %s too big\n", program, z);
1.1 cgd 280: continue;
281:
282: /* carefully check against near limit numbers */
283: } else if (strcmp(z, limit) > 0) {
1.3 ! cgd 284: fprintf(stderr, "%s: %s a bit too big\n", program, z);
1.1 cgd 285: continue;
286: }
287: /* number is near limit, but is under it */
288: return s;
289: } while (input != NULL && fgets(buf, MAX_LINE, input) != NULL);
290:
291: /* error or EOF */
292: return NULL;
293: }
294:
295: /*
296: * primes - sieve and print primes from start up to and but not including stop
297: */
298: void
299: primes(start, stop)
300: ubig start; /* where to start generating */
301: ubig stop; /* don't generate at or above this value */
302: {
303: register char *q; /* sieve spot */
304: register ubig factor; /* index and factor */
305: register char *tab_lim; /* the limit to sieve on the table */
306: register ubig *p; /* prime table pointer */
307: register ubig fact_lim; /* highest prime for current block */
308:
309: /*
1.3 ! cgd 310: * NetBSD has no problems with handling conversion
! 311: * between doubles and unsigned long, so we can go
! 312: * all the way to BIG.
1.1 cgd 313: */
314: if (start < 3) {
315: start = (ubig)2;
316: }
317: if (stop < 3) {
318: stop = (ubig)2;
319: }
320: if (stop <= start) {
321: return;
322: }
323:
324: /*
325: * be sure that the values are odd, or 2
326: */
327: if (start != 2 && (start&0x1) == 0) {
328: ++start;
329: }
330: if (stop != 2 && (stop&0x1) == 0) {
331: ++stop;
332: }
333:
334: /*
335: * quick list of primes <= pr_limit
336: */
337: if (start <= *pr_limit) {
338: /* skip primes up to the start value */
339: for (p = &prime[0], factor = prime[0];
340: factor < stop && p <= pr_limit;
341: factor = *(++p)) {
342: if (factor >= start) {
343: printf("%u\n", factor);
344: }
345: }
346: /* return early if we are done */
347: if (p <= pr_limit) {
348: return;
349: }
350: start = *pr_limit+2;
351: }
352:
353: /*
354: * we shall sieve a bytemap window, note primes and move the window
355: * upward until we pass the stop point
356: */
357: while (start < stop) {
358: /*
359: * factor out 3, 5, 7, 11 and 13
360: */
361: /* initial pattern copy */
362: factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
363: memcpy(table, &pattern[factor], pattern_size-factor);
364: /* main block pattern copies */
365: for (fact_lim=pattern_size-factor;
366: fact_lim+pattern_size<=TABSIZE;
367: fact_lim+=pattern_size) {
368: memcpy(&table[fact_lim], pattern, pattern_size);
369: }
370: /* final block pattern copy */
371: memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
372:
373: /*
374: * sieve for primes 17 and higher
375: */
376: /* note highest useful factor and sieve spot */
377: if (stop-start > TABSIZE+TABSIZE) {
378: tab_lim = &table[TABSIZE]; /* sieve it all */
379: fact_lim = (int)sqrt(
380: (double)(start)+TABSIZE+TABSIZE+1.0);
381: } else {
382: tab_lim = &table[(stop-start)/2]; /* partial sieve */
383: fact_lim = (int)sqrt((double)(stop)+1.0);
384: }
385: /* sieve for factors >= 17 */
386: factor = 17; /* 17 is first prime to use */
387: p = &prime[7]; /* 19 is next prime, pi(19)=7 */
388: do {
389: /* determine the factor's initial sieve point */
390: q = (char *)(start%factor); /* temp storage for mod */
391: if ((int)q & 0x1) {
392: q = &table[(factor-(int)q)/2];
393: } else {
394: q = &table[q ? factor-((int)q/2) : 0];
395: }
396: /* sive for our current factor */
397: for ( ; q < tab_lim; q += factor) {
398: *q = '\0'; /* sieve out a spot */
399: }
400: } while ((factor=(ubig)(*(p++))) <= fact_lim);
401:
402: /*
403: * print generated primes
404: */
405: for (q = table; q < tab_lim; ++q, start+=2) {
406: if (*q) {
407: printf("%u\n", start);
408: }
409: }
410: }
411: }
CVSweb <webmaster@jp.NetBSD.org>
![[BACK]](/icons/back.gif){kind=icon}
Return to primes.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [cvs.NetBSD.org] / src / games / primes