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