1 /* $NetBSD: primes.c,v 1.15 2008/02/02 18:15:14 matt Exp $ */
4 * Copyright (c) 1989, 1993
5 * The Regents of the University of California. All rights reserved.
7 * This code is derived from software contributed to Berkeley by
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
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.
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
35 #include <sys/cdefs.h>
37 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
38 The Regents of the University of California. All rights reserved.");
43 static char sccsid
[] = "@(#)primes.c 8.5 (Berkeley) 5/10/95";
45 __RCSID("$NetBSD: primes.c,v 1.15 2008/02/02 18:15:14 matt Exp $");
50 * primes - generate a table of primes between two values
52 * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo
54 * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
57 * primes [start [stop]]
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.
63 * validation check: there are 664579 primes between 0 and 10^7
79 * Eratosthenes sieve table
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.
85 * We make TABSIZE large to reduce the overhead of inner loop setup.
87 char table
[TABSIZE
]; /* Eratosthenes sieve of odd numbers */
90 * prime[i] is the (i-1)th prime.
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.
95 extern const ubig prime
[];
96 extern const ubig
*pr_limit
; /* largest prime in the prime array */
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.
103 extern const char pattern
[];
104 extern const int pattern_size
; /* length of pattern array */
108 int main(int, char *[]);
109 void primes(ubig
, ubig
);
110 ubig
read_num_buf(void);
111 void usage(void) __dead
;
114 main(int argc
, char *argv
[])
116 ubig start
; /* where to start generating */
117 ubig stop
; /* don't generate at or above this value */
121 while ((ch
= getopt(argc
, argv
, "d")) != -1)
137 * Convert low and high args. Strtoul(3) sets errno to
138 * ERANGE if the number is too large, but, if there's
139 * a leading minus sign it returns the negation of the
140 * result of the conversion, which we'd rather disallow.
144 /* Start and stop supplied on the command line. */
145 if (argv
[0][0] == '-' || argv
[1][0] == '-')
146 errx(1, "negative numbers aren't permitted.");
149 start
= strtoul(argv
[0], &p
, 10);
151 err(1, "%s", argv
[0]);
153 errx(1, "%s: illegal numeric format.", argv
[0]);
156 stop
= strtoul(argv
[1], &p
, 10);
158 err(1, "%s", argv
[1]);
160 errx(1, "%s: illegal numeric format.", argv
[1]);
163 /* Start on the command line. */
164 if (argv
[0][0] == '-')
165 errx(1, "negative numbers aren't permitted.");
168 start
= strtoul(argv
[0], &p
, 10);
170 err(1, "%s", argv
[0]);
172 errx(1, "%s: illegal numeric format.", argv
[0]);
175 start
= read_num_buf();
182 errx(1, "start value must be less than stop value.");
189 * This routine returns a number n, where 0 <= n && n <= BIG.
195 char *p
, buf
[100]; /* > max number of digits. */
198 if (fgets(buf
, sizeof(buf
), stdin
) == NULL
) {
203 for (p
= buf
; isblank(*p
); ++p
);
204 if (*p
== '\n' || *p
== '\0')
207 errx(1, "negative numbers aren't permitted.");
209 val
= strtoul(buf
, &p
, 10);
213 errx(1, "%s: illegal numeric format.", buf
);
219 * primes - sieve and print primes from start up to and but not including stop
221 * start where to start generating
222 * stop don't generate at or above this value
225 primes(ubig start
, ubig stop
)
227 char *q
; /* sieve spot */
228 ubig factor
; /* index and factor */
229 char *tab_lim
; /* the limit to sieve on the table */
230 const ubig
*p
; /* prime table pointer */
231 ubig fact_lim
; /* highest prime for current block */
232 ubig mod
; /* temp storage for mod */
236 * A number of systems can not convert double values into unsigned
237 * longs when the values are larger than the largest signed value.
238 * We don't have this problem, so we can go all the way to BIG.
251 * be sure that the values are odd, or 2
253 if (start
!= 2 && (start
&0x1) == 0) {
256 if (stop
!= 2 && (stop
&0x1) == 0) {
261 * quick list of primes <= pr_limit
263 if (start
<= *pr_limit
) {
264 /* skip primes up to the start value */
265 for (p
= &prime
[0], factor
= prime
[0];
266 factor
< stop
&& p
<= pr_limit
; factor
= *(++p
)) {
267 if (factor
>= start
) {
268 printf("%lu", (unsigned long) factor
);
271 (unsigned long) factor
- prev
);
277 /* return early if we are done */
285 * we shall sieve a bytemap window, note primes and move the window
286 * upward until we pass the stop point
288 while (start
< stop
) {
290 * factor out 3, 5, 7, 11 and 13
292 /* initial pattern copy */
293 factor
= (start
%(2*3*5*7*11*13))/2; /* starting copy spot */
294 memcpy(table
, &pattern
[factor
], pattern_size
-factor
);
295 /* main block pattern copies */
296 for (fact_lim
=pattern_size
-factor
;
297 fact_lim
+pattern_size
<=TABSIZE
; fact_lim
+=pattern_size
) {
298 memcpy(&table
[fact_lim
], pattern
, pattern_size
);
300 /* final block pattern copy */
301 memcpy(&table
[fact_lim
], pattern
, TABSIZE
-fact_lim
);
304 * sieve for primes 17 and higher
306 /* note highest useful factor and sieve spot */
307 if (stop
-start
> TABSIZE
+TABSIZE
) {
308 tab_lim
= &table
[TABSIZE
]; /* sieve it all */
309 fact_lim
= (int)sqrt(
310 (double)(start
)+TABSIZE
+TABSIZE
+1.0);
312 tab_lim
= &table
[(stop
-start
)/2]; /* partial sieve */
313 fact_lim
= (int)sqrt((double)(stop
)+1.0);
315 /* sieve for factors >= 17 */
316 factor
= 17; /* 17 is first prime to use */
317 p
= &prime
[7]; /* 19 is next prime, pi(19)=7 */
319 /* determine the factor's initial sieve point */
322 q
= &table
[(factor
-mod
)/2];
324 q
= &table
[mod
? factor
-(mod
/2) : 0];
326 /* sieve for our current factor */
327 for ( ; q
< tab_lim
; q
+= factor
) {
328 *q
= '\0'; /* sieve out a spot */
330 } while ((factor
=(ubig
)(*(p
++))) <= fact_lim
);
333 * print generated primes
335 for (q
= table
; q
< tab_lim
; ++q
, start
+=2) {
337 printf("%lu", (unsigned long) start
);
340 (unsigned long) start
- prev
);
352 (void)fprintf(stderr
, "usage: primes [-d] [start [stop]]\n");