1/*	$NetBSD: primes.c,v 1.18 2010/05/13 17:52:12 tnozaki 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. Neither the name of the University nor the names of its contributors
19 *    may be used to endorse or promote products derived from this software
20 *    without specific prior written permission.
21 *
22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32 * SUCH DAMAGE.
33 */
34
35#include <sys/cdefs.h>
36#ifndef lint
37__COPYRIGHT("@(#) Copyright (c) 1989, 1993\
38 The Regents of the University of California.  All rights reserved.");
39#endif /* not lint */
40
41#ifndef lint
42#if 0
43static char sccsid[] = "@(#)primes.c	8.5 (Berkeley) 5/10/95";
44#else
45__RCSID("$NetBSD: primes.c,v 1.18 2010/05/13 17:52:12 tnozaki Exp $");
46#endif
47#endif /* not lint */
48
49/*
50 * primes - generate a table of primes between two values
51 *
52 * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo
53 *
54 * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
55 *
56 * usage:
57 *	primes [start [stop]]
58 *
59 *	Print primes >= start and < stop.  If stop is omitted,
60 *	the value 4294967295 (2^32-1) is assumed.  If start is
61 *	omitted, start is read from standard input.
62 *
63 * validation check: there are 664579 primes between 0 and 10^7
64 */
65
66#include <ctype.h>
67#include <err.h>
68#include <errno.h>
69#include <limits.h>
70#include <math.h>
71#include <memory.h>
72#include <stdio.h>
73#include <stdlib.h>
74#include <unistd.h>
75
76#include "primes.h"
77
78/*
79 * Eratosthenes sieve table
80 *
81 * We only sieve the odd numbers.  The base of our sieve windows are always
82 * odd.  If the base of table is 1, table[i] represents 2*i-1.  After the
83 * sieve, table[i] == 1 if and only iff 2*i-1 is prime.
84 *
85 * We make TABSIZE large to reduce the overhead of inner loop setup.
86 */
87static char table[TABSIZE];	 /* Eratosthenes sieve of odd numbers */
88
89/*
90 * prime[i] is the (i-1)th prime.
91 *
92 * We are able to sieve 2^32-1 because this byte table yields all primes
93 * up to 65537 and 65537^2 > 2^32-1.
94 */
95extern const ubig prime[];
96extern const ubig *pr_limit;		/* largest prime in the prime array */
97
98/*
99 * To avoid excessive sieves for small factors, we use the table below to
100 * setup our sieve blocks.  Each element represents a odd number starting
101 * with 1.  All non-zero elements are factors of 3, 5, 7, 11 and 13.
102 */
103extern const char pattern[];
104extern const int pattern_size;	/* length of pattern array */
105
106static int dflag;
107
108static void primes(ubig, ubig);
109static ubig read_num_buf(void);
110static void usage(void) __dead;
111
112int
113main(int argc, char *argv[])
114{
115	ubig start;		/* where to start generating */
116	ubig stop;		/* don't generate at or above this value */
117	int ch;
118	char *p;
119
120	while ((ch = getopt(argc, argv, "d")) != -1)
121		switch (ch) {
122		case 'd':
123			dflag++;
124			break;
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 */
190ubig
191read_num_buf(void)
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((unsigned char)*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 *	start	where to start generating
221 *	stop	don't generate at or above this value
222 */
223void
224primes(ubig start, ubig stop)
225{
226	char *q;		/* sieve spot */
227	ubig factor;		/* index and factor */
228	char *tab_lim;		/* the limit to sieve on the table */
229	const ubig *p;		/* prime table pointer */
230	ubig fact_lim;		/* highest prime for current block */
231	ubig mod;		/* temp storage for mod */
232	ubig prev = 0;
233
234	/*
235	 * A number of systems can not convert double values into unsigned
236	 * longs when the values are larger than the largest signed value.
237	 * We don't have this problem, so we can go all the way to BIG.
238	 */
239	if (start < 3) {
240		start = (ubig)2;
241	}
242	if (stop < 3) {
243		stop = (ubig)2;
244	}
245	if (stop <= start) {
246		return;
247	}
248
249	/*
250	 * be sure that the values are odd, or 2
251	 */
252	if (start != 2 && (start&0x1) == 0) {
253		++start;
254	}
255	if (stop != 2 && (stop&0x1) == 0) {
256		++stop;
257	}
258
259	/*
260	 * quick list of primes <= pr_limit
261	 */
262	if (start <= *pr_limit) {
263		/* skip primes up to the start value */
264		for (p = &prime[0], factor = prime[0];
265		    factor < stop && p <= pr_limit; factor = *(++p)) {
266			if (factor >= start) {
267				printf("%lu", (unsigned long) factor);
268				if (dflag) {
269					printf(" (%lu)",
270					    (unsigned long) factor - prev);
271				}
272				putchar('\n');
273			}
274			prev = factor;
275		}
276		/* return early if we are done */
277		if (p <= pr_limit) {
278			return;
279		}
280		start = *pr_limit+2;
281	}
282
283	/*
284	 * we shall sieve a bytemap window, note primes and move the window
285	 * upward until we pass the stop point
286	 */
287	while (start < stop) {
288		/*
289		 * factor out 3, 5, 7, 11 and 13
290		 */
291		/* initial pattern copy */
292		factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
293		memcpy(table, &pattern[factor], pattern_size-factor);
294		/* main block pattern copies */
295		for (fact_lim=pattern_size-factor;
296		    fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) {
297			memcpy(&table[fact_lim], pattern, pattern_size);
298		}
299		/* final block pattern copy */
300		memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
301
302		/*
303		 * sieve for primes 17 and higher
304		 */
305		/* note highest useful factor and sieve spot */
306		if (stop-start > TABSIZE+TABSIZE) {
307			tab_lim = &table[TABSIZE]; /* sieve it all */
308			fact_lim = sqrt((double)(start)+TABSIZE+TABSIZE+1.0);
309		} else {
310			tab_lim = &table[(stop-start)/2]; /* partial sieve */
311			fact_lim = sqrt((double)(stop)+1.0);
312		}
313		/* sieve for factors >= 17 */
314		factor = 17;	/* 17 is first prime to use */
315		p = &prime[7];	/* 19 is next prime, pi(19)=7 */
316		do {
317			/* determine the factor's initial sieve point */
318			mod = start%factor;
319			if (mod & 0x1) {
320				q = &table[(factor-mod)/2];
321			} else {
322				q = &table[mod ? factor-(mod/2) : 0];
323			}
324			/* sieve for our current factor */
325			for ( ; q < tab_lim; q += factor) {
326				*q = '\0'; /* sieve out a spot */
327			}
328		} while ((factor=(ubig)(*(p++))) <= fact_lim);
329
330		/*
331		 * print generated primes
332		 */
333		for (q = table; q < tab_lim; ++q, start+=2) {
334			if (*q) {
335				printf("%lu", (unsigned long) start);
336				if (dflag) {
337					printf(" (%lu)",
338					    (unsigned long) start - prev);
339					prev = start;
340				}
341				putchar('\n');
342			}
343		}
344	}
345}
346
347void
348usage(void)
349{
350	(void)fprintf(stderr, "usage: primes [-d] [start [stop]]\n");
351	exit(1);
352}
353