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