primes.c revision 91027
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 38static 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 45static char sccsid[] = "@(#)primes.c 8.5 (Berkeley) 5/10/95"; 46#endif 47static const char rcsid[] = 48 "$FreeBSD: head/games/primes/primes.c 91027 2002-02-21 18:13:31Z billf $"; 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 [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 <memory.h> 74#include <stdio.h> 75#include <stdlib.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 */ 89static char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 90 91/* 92 * prime[i] is the (i-1)th prime. 93 * 94 * We are able to sieve 2^32-1 because this byte table yields all primes 95 * up to 65537 and 65537^2 > 2^32-1. 96 */ 97extern ubig prime[]; 98extern ubig *pr_limit; /* largest prime in the prime array */ 99 100/* 101 * To avoid excessive sieves for small factors, we use the table below to 102 * setup our sieve blocks. Each element represents a odd number starting 103 * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13. 104 */ 105extern char pattern[]; 106extern size_t pattern_size; /* length of pattern array */ 107 108static int hflag; 109 110static void primes(ubig, ubig); 111static ubig read_num_buf(void); 112static void usage(void); 113 114int 115main(int argc, char *argv[]) 116{ 117 ubig start; /* where to start generating */ 118 ubig stop; /* don't generate at or above this value */ 119 int ch; 120 char *p; 121 122 while ((ch = getopt(argc, argv, "h")) != -1) 123 switch (ch) { 124 case 'h': 125 hflag++; 126 break; 127 case '?': 128 default: 129 usage(); 130 } 131 argc -= optind; 132 argv += optind; 133 134 start = 0; 135 stop = BIG; 136 137 /* 138 * Convert low and high args. Strtoul(3) sets errno to 139 * ERANGE if the number is too large, but, if there's 140 * a leading minus sign it returns the negation of the 141 * result of the conversion, which we'd rather disallow. 142 */ 143 switch (argc) { 144 case 2: 145 /* Start and stop supplied on the command line. */ 146 if (argv[0][0] == '-' || argv[1][0] == '-') 147 errx(1, "negative numbers aren't permitted."); 148 149 errno = 0; 150 start = strtoul(argv[0], &p, 0); 151 if (errno) 152 err(1, "%s", argv[0]); 153 if (*p != '\0') 154 errx(1, "%s: illegal numeric format.", argv[0]); 155 156 errno = 0; 157 stop = strtoul(argv[1], &p, 0); 158 if (errno) 159 err(1, "%s", argv[1]); 160 if (*p != '\0') 161 errx(1, "%s: illegal numeric format.", argv[1]); 162 break; 163 case 1: 164 /* Start on the command line. */ 165 if (argv[0][0] == '-') 166 errx(1, "negative numbers aren't permitted."); 167 168 errno = 0; 169 start = strtoul(argv[0], &p, 0); 170 if (errno) 171 err(1, "%s", argv[0]); 172 if (*p != '\0') 173 errx(1, "%s: illegal numeric format.", argv[0]); 174 break; 175 case 0: 176 start = read_num_buf(); 177 break; 178 default: 179 usage(); 180 } 181 182 if (start > stop) 183 errx(1, "start value must be less than stop value."); 184 primes(start, stop); 185 return (0); 186} 187 188/* 189 * read_num_buf -- 190 * This routine returns a number n, where 0 <= n && n <= BIG. 191 */ 192static ubig 193read_num_buf(void) 194{ 195 ubig val; 196 char *p, buf[100]; /* > max number of digits. */ 197 198 for (;;) { 199 if (fgets(buf, sizeof(buf), stdin) == NULL) { 200 if (ferror(stdin)) 201 err(1, "stdin"); 202 exit(0); 203 } 204 for (p = buf; isblank(*p); ++p); 205 if (*p == '\n' || *p == '\0') 206 continue; 207 if (*p == '-') 208 errx(1, "negative numbers aren't permitted."); 209 errno = 0; 210 val = strtoul(buf, &p, 0); 211 if (errno) 212 err(1, "%s", buf); 213 if (*p != '\n') 214 errx(1, "%s: illegal numeric format.", buf); 215 return (val); 216 } 217} 218 219/* 220 * primes - sieve and print primes from start up to and but not including stop 221 */ 222static void 223primes(ubig start, ubig stop) 224{ 225 char *q; /* sieve spot */ 226 ubig factor; /* index and factor */ 227 char *tab_lim; /* the limit to sieve on the table */ 228 ubig *p; /* prime table pointer */ 229 ubig fact_lim; /* highest prime for current block */ 230 231 /* 232 * A number of systems can not convert double values into unsigned 233 * longs when the values are larger than the largest signed value. 234 * We don't have this problem, so we can go all the way to BIG. 235 */ 236 if (start < 3) { 237 start = (ubig)2; 238 } 239 if (stop < 3) { 240 stop = (ubig)2; 241 } 242 if (stop <= start) { 243 return; 244 } 245 246 /* 247 * be sure that the values are odd, or 2 248 */ 249 if (start != 2 && (start&0x1) == 0) { 250 ++start; 251 } 252 if (stop != 2 && (stop&0x1) == 0) { 253 ++stop; 254 } 255 256 /* 257 * quick list of primes <= pr_limit 258 */ 259 if (start <= *pr_limit) { 260 /* skip primes up to the start value */ 261 for (p = &prime[0], factor = prime[0]; 262 factor < stop && p <= pr_limit; factor = *(++p)) { 263 if (factor >= start) { 264 printf(hflag ? "0x%lx\n" : "%lu\n", factor); 265 } 266 } 267 /* return early if we are done */ 268 if (p <= pr_limit) { 269 return; 270 } 271 start = *pr_limit+2; 272 } 273 274 /* 275 * we shall sieve a bytemap window, note primes and move the window 276 * upward until we pass the stop point 277 */ 278 while (start < stop) { 279 /* 280 * factor out 3, 5, 7, 11 and 13 281 */ 282 /* initial pattern copy */ 283 factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 284 memcpy(table, &pattern[factor], pattern_size-factor); 285 /* main block pattern copies */ 286 for (fact_lim=pattern_size-factor; 287 fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 288 memcpy(&table[fact_lim], pattern, pattern_size); 289 } 290 /* final block pattern copy */ 291 memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 292 293 /* 294 * sieve for primes 17 and higher 295 */ 296 /* note highest useful factor and sieve spot */ 297 if (stop-start > TABSIZE+TABSIZE) { 298 tab_lim = &table[TABSIZE]; /* sieve it all */ 299 fact_lim = (int)sqrt( 300 (double)(start)+TABSIZE+TABSIZE+1.0); 301 } else { 302 tab_lim = &table[(stop-start)/2]; /* partial sieve */ 303 fact_lim = (int)sqrt((double)(stop)+1.0); 304 } 305 /* sieve for factors >= 17 */ 306 factor = 17; /* 17 is first prime to use */ 307 p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 308 do { 309 /* determine the factor's initial sieve point */ 310 q = (char *)(start%factor); /* temp storage for mod */ 311 if ((long)q & 0x1) { 312 q = &table[(factor-(long)q)/2]; 313 } else { 314 q = &table[q ? factor-((long)q/2) : 0]; 315 } 316 /* sive for our current factor */ 317 for ( ; q < tab_lim; q += factor) { 318 *q = '\0'; /* sieve out a spot */ 319 } 320 } while ((factor=(ubig)(*(p++))) <= fact_lim); 321 322 /* 323 * print generated primes 324 */ 325 for (q = table; q < tab_lim; ++q, start+=2) { 326 if (*q) { 327 printf(hflag ? "0x%lx\n" : "%lu\n", start); 328 } 329 } 330 } 331} 332 333static void 334usage(void) 335{ 336 (void)fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); 337 exit(1); 338} 339