Commit | Line | Data |
---|---|---|
e0bbfbf9 | 1 | /* |
31575405 KB |
2 | * Copyright (c) 1989, 1993 |
3 | * The Regents of the University of California. All rights reserved. | |
bc2e2efb KB |
4 | * |
5 | * This code is derived from software contributed to Berkeley by | |
6 | * Landon Curt Noll. | |
7 | * | |
ad787160 C |
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. | |
e0bbfbf9 DF |
35 | */ |
36 | ||
37 | #ifndef lint | |
31575405 KB |
38 | static char copyright[] = |
39 | "@(#) Copyright (c) 1989, 1993\n\ | |
40 | The Regents of the University of California. All rights reserved.\n"; | |
bc2e2efb | 41 | #endif /* not lint */ |
e0bbfbf9 | 42 | |
779464d0 | 43 | #ifndef lint |
ed554bc5 | 44 | static char sccsid[] = "@(#)primes.c 8.4 (Berkeley) 3/21/94"; |
bc2e2efb | 45 | #endif /* not lint */ |
779464d0 SL |
46 | |
47 | /* | |
bc2e2efb KB |
48 | * primes - generate a table of primes between two values |
49 | * | |
39bd0b9c | 50 | * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo |
779464d0 | 51 | * |
39bd0b9c | 52 | * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ |
779464d0 | 53 | * |
bc2e2efb KB |
54 | * usage: |
55 | * primes [start [stop]] | |
779464d0 | 56 | * |
18a3e3ca | 57 | * Print primes >= start and < stop. If stop is omitted, |
bc2e2efb KB |
58 | * the value 4294967295 (2^32-1) is assumed. If start is |
59 | * omitted, start is read from standard input. | |
60 | * | |
18a3e3ca | 61 | * validation check: there are 664579 primes between 0 and 10^7 |
779464d0 SL |
62 | */ |
63 | ||
62d14cd7 | 64 | #include <ctype.h> |
39bd0b9c KB |
65 | #include <err.h> |
66 | #include <errno.h> | |
39bd0b9c | 67 | #include <limits.h> |
779464d0 | 68 | #include <math.h> |
bc2e2efb | 69 | #include <memory.h> |
39bd0b9c KB |
70 | #include <stdio.h> |
71 | #include <stdlib.h> | |
62d14cd7 | 72 | |
bc2e2efb | 73 | #include "primes.h" |
779464d0 | 74 | |
bc2e2efb KB |
75 | /* |
76 | * Eratosthenes sieve table | |
77 | * | |
78 | * We only sieve the odd numbers. The base of our sieve windows are always | |
79 | * odd. If the base of table is 1, table[i] represents 2*i-1. After the | |
80 | * sieve, table[i] == 1 if and only iff 2*i-1 is prime. | |
81 | * | |
82 | * We make TABSIZE large to reduce the overhead of inner loop setup. | |
83 | */ | |
84 | char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ | |
779464d0 | 85 | |
bc2e2efb KB |
86 | /* |
87 | * prime[i] is the (i-1)th prime. | |
88 | * | |
89 | * We are able to sieve 2^32-1 because this byte table yields all primes | |
90 | * up to 65537 and 65537^2 > 2^32-1. | |
91 | */ | |
92 | extern ubig prime[]; | |
93 | extern ubig *pr_limit; /* largest prime in the prime array */ | |
779464d0 | 94 | |
bc2e2efb KB |
95 | /* |
96 | * To avoid excessive sieves for small factors, we use the table below to | |
97 | * setup our sieve blocks. Each element represents a odd number starting | |
98 | * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13. | |
99 | */ | |
100 | extern char pattern[]; | |
101 | extern int pattern_size; /* length of pattern array */ | |
779464d0 | 102 | |
62d14cd7 KB |
103 | void primes __P((ubig, ubig)); |
104 | ubig read_num_buf __P((void)); | |
105 | void usage __P((void)); | |
779464d0 | 106 | |
39bd0b9c | 107 | int |
779464d0 | 108 | main(argc, argv) |
39bd0b9c KB |
109 | int argc; |
110 | char *argv[]; | |
779464d0 | 111 | { |
39bd0b9c KB |
112 | ubig start; /* where to start generating */ |
113 | ubig stop; /* don't generate at or above this value */ | |
114 | int ch; | |
62d14cd7 | 115 | char *p; |
39bd0b9c KB |
116 | |
117 | while ((ch = getopt(argc, argv, "")) != EOF) | |
118 | switch (ch) { | |
119 | case '?': | |
120 | default: | |
121 | usage(); | |
122 | } | |
123 | argc -= optind; | |
124 | argv += optind; | |
bc2e2efb | 125 | |
bc2e2efb KB |
126 | start = 0; |
127 | stop = BIG; | |
62d14cd7 KB |
128 | |
129 | /* | |
130 | * Convert low and high args. Strtoul(3) sets errno to | |
131 | * ERANGE if the number is too large, but, if there's | |
132 | * a leading minus sign it returns the negation of the | |
133 | * result of the conversion, which we'd rather disallow. | |
134 | */ | |
135 | switch (argc) { | |
136 | case 2: | |
137 | /* Start and stop supplied on the command line. */ | |
138 | if (argv[0][0] == '-' || argv[1][0] == '-') | |
139 | errx(1, "negative numbers aren't permitted."); | |
140 | ||
141 | errno = 0; | |
142 | start = strtoul(argv[0], &p, 10); | |
143 | if (errno) | |
39bd0b9c | 144 | err(1, "%s", argv[0]); |
62d14cd7 KB |
145 | if (*p != '\0') |
146 | errx(1, "%s: illegal numeric format.", argv[0]); | |
147 | ||
148 | errno = 0; | |
149 | stop = strtoul(argv[1], &p, 10); | |
150 | if (errno) | |
151 | err(1, "%s", argv[1]); | |
152 | if (*p != '\0') | |
153 | errx(1, "%s: illegal numeric format.", argv[1]); | |
154 | break; | |
155 | case 1: | |
156 | /* Start on the command line. */ | |
157 | if (argv[0][0] == '-') | |
158 | errx(1, "negative numbers aren't permitted."); | |
159 | ||
160 | errno = 0; | |
161 | start = strtoul(argv[0], &p, 10); | |
162 | if (errno) | |
39bd0b9c | 163 | err(1, "%s", argv[0]); |
62d14cd7 KB |
164 | if (*p != '\0') |
165 | errx(1, "%s: illegal numeric format.", argv[0]); | |
166 | break; | |
167 | case 0: | |
168 | start = read_num_buf(); | |
169 | break; | |
170 | default: | |
171 | usage(); | |
bc2e2efb | 172 | } |
62d14cd7 | 173 | |
39bd0b9c | 174 | if (start > stop) |
62d14cd7 | 175 | errx(1, "start value must be less than stop value."); |
bc2e2efb KB |
176 | primes(start, stop); |
177 | exit(0); | |
178 | } | |
779464d0 | 179 | |
bc2e2efb | 180 | /* |
62d14cd7 KB |
181 | * read_num_buf -- |
182 | * This routine returns a number n, where 0 <= n && n <= BIG. | |
bc2e2efb | 183 | */ |
62d14cd7 KB |
184 | ubig |
185 | read_num_buf() | |
bc2e2efb | 186 | { |
62d14cd7 KB |
187 | ubig val; |
188 | char *p, buf[100]; /* > max number of digits. */ | |
bc2e2efb | 189 | |
62d14cd7 KB |
190 | for (;;) { |
191 | if (fgets(buf, sizeof(buf), stdin) == NULL) { | |
192 | if (ferror(stdin)) | |
193 | err(1, "stdin"); | |
7e439f5a | 194 | exit(0); |
bc2e2efb | 195 | } |
62d14cd7 KB |
196 | for (p = buf; isblank(*p); ++p); |
197 | if (*p == '\n' || *p == '\0') | |
bc2e2efb | 198 | continue; |
62d14cd7 KB |
199 | if (*p == '-') |
200 | errx(1, "negative numbers aren't permitted."); | |
201 | errno = 0; | |
202 | val = strtoul(buf, &p, 10); | |
203 | if (errno) | |
204 | err(1, "%s", buf); | |
205 | if (*p != '\n') | |
206 | errx(1, "%s: illegal numeric format.", buf); | |
207 | return (val); | |
208 | } | |
bc2e2efb | 209 | } |
779464d0 SL |
210 | |
211 | /* | |
bc2e2efb | 212 | * primes - sieve and print primes from start up to and but not including stop |
779464d0 | 213 | */ |
bc2e2efb KB |
214 | void |
215 | primes(start, stop) | |
216 | ubig start; /* where to start generating */ | |
217 | ubig stop; /* don't generate at or above this value */ | |
218 | { | |
219 | register char *q; /* sieve spot */ | |
220 | register ubig factor; /* index and factor */ | |
221 | register char *tab_lim; /* the limit to sieve on the table */ | |
222 | register ubig *p; /* prime table pointer */ | |
223 | register ubig fact_lim; /* highest prime for current block */ | |
224 | ||
779464d0 | 225 | /* |
39bd0b9c KB |
226 | * A number of systems can not convert double values into unsigned |
227 | * longs when the values are larger than the largest signed value. | |
228 | * We don't have this problem, so we can go all the way to BIG. | |
779464d0 | 229 | */ |
bc2e2efb KB |
230 | if (start < 3) { |
231 | start = (ubig)2; | |
232 | } | |
233 | if (stop < 3) { | |
234 | stop = (ubig)2; | |
235 | } | |
236 | if (stop <= start) { | |
237 | return; | |
238 | } | |
239 | ||
779464d0 | 240 | /* |
bc2e2efb | 241 | * be sure that the values are odd, or 2 |
779464d0 | 242 | */ |
bc2e2efb KB |
243 | if (start != 2 && (start&0x1) == 0) { |
244 | ++start; | |
245 | } | |
246 | if (stop != 2 && (stop&0x1) == 0) { | |
247 | ++stop; | |
779464d0 | 248 | } |
779464d0 | 249 | |
bc2e2efb KB |
250 | /* |
251 | * quick list of primes <= pr_limit | |
252 | */ | |
253 | if (start <= *pr_limit) { | |
254 | /* skip primes up to the start value */ | |
255 | for (p = &prime[0], factor = prime[0]; | |
39bd0b9c | 256 | factor < stop && p <= pr_limit; factor = *(++p)) { |
bc2e2efb KB |
257 | if (factor >= start) { |
258 | printf("%u\n", factor); | |
259 | } | |
260 | } | |
261 | /* return early if we are done */ | |
262 | if (p <= pr_limit) { | |
263 | return; | |
264 | } | |
265 | start = *pr_limit+2; | |
779464d0 | 266 | } |
779464d0 | 267 | |
bc2e2efb KB |
268 | /* |
269 | * we shall sieve a bytemap window, note primes and move the window | |
270 | * upward until we pass the stop point | |
271 | */ | |
272 | while (start < stop) { | |
273 | /* | |
274 | * factor out 3, 5, 7, 11 and 13 | |
275 | */ | |
276 | /* initial pattern copy */ | |
277 | factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ | |
278 | memcpy(table, &pattern[factor], pattern_size-factor); | |
279 | /* main block pattern copies */ | |
280 | for (fact_lim=pattern_size-factor; | |
39bd0b9c | 281 | fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { |
bc2e2efb KB |
282 | memcpy(&table[fact_lim], pattern, pattern_size); |
283 | } | |
284 | /* final block pattern copy */ | |
285 | memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); | |
286 | ||
287 | /* | |
288 | * sieve for primes 17 and higher | |
289 | */ | |
290 | /* note highest useful factor and sieve spot */ | |
291 | if (stop-start > TABSIZE+TABSIZE) { | |
292 | tab_lim = &table[TABSIZE]; /* sieve it all */ | |
293 | fact_lim = (int)sqrt( | |
294 | (double)(start)+TABSIZE+TABSIZE+1.0); | |
295 | } else { | |
296 | tab_lim = &table[(stop-start)/2]; /* partial sieve */ | |
297 | fact_lim = (int)sqrt((double)(stop)+1.0); | |
298 | } | |
299 | /* sieve for factors >= 17 */ | |
300 | factor = 17; /* 17 is first prime to use */ | |
301 | p = &prime[7]; /* 19 is next prime, pi(19)=7 */ | |
302 | do { | |
303 | /* determine the factor's initial sieve point */ | |
304 | q = (char *)(start%factor); /* temp storage for mod */ | |
305 | if ((int)q & 0x1) { | |
306 | q = &table[(factor-(int)q)/2]; | |
307 | } else { | |
308 | q = &table[q ? factor-((int)q/2) : 0]; | |
309 | } | |
310 | /* sive for our current factor */ | |
311 | for ( ; q < tab_lim; q += factor) { | |
312 | *q = '\0'; /* sieve out a spot */ | |
313 | } | |
314 | } while ((factor=(ubig)(*(p++))) <= fact_lim); | |
315 | ||
316 | /* | |
317 | * print generated primes | |
318 | */ | |
319 | for (q = table; q < tab_lim; ++q, start+=2) { | |
320 | if (*q) { | |
321 | printf("%u\n", start); | |
322 | } | |
323 | } | |
324 | } | |
779464d0 | 325 | } |
39bd0b9c KB |
326 | |
327 | void | |
328 | usage() | |
329 | { | |
62d14cd7 | 330 | (void)fprintf(stderr, "usage: primes [start [stop]]\n"); |
39bd0b9c KB |
331 | exit(1); |
332 | } |