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