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