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.2 (Berkeley) 03/01/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 <err.h> 39 #include <errno.h> 40 #include <ctype.h> 41 #include <limits.h> 42 #include <math.h> 43 #include <memory.h> 44 #include <stdio.h> 45 #include <stdlib.h> 46 #include "primes.h" 47 48 /* 49 * Eratosthenes sieve table 50 * 51 * We only sieve the odd numbers. The base of our sieve windows are always 52 * odd. If the base of table is 1, table[i] represents 2*i-1. After the 53 * sieve, table[i] == 1 if and only iff 2*i-1 is prime. 54 * 55 * We make TABSIZE large to reduce the overhead of inner loop setup. 56 */ 57 char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 58 59 /* 60 * prime[i] is the (i-1)th prime. 61 * 62 * We are able to sieve 2^32-1 because this byte table yields all primes 63 * up to 65537 and 65537^2 > 2^32-1. 64 */ 65 extern ubig prime[]; 66 extern ubig *pr_limit; /* largest prime in the prime array */ 67 68 /* 69 * To avoid excessive sieves for small factors, we use the table below to 70 * setup our sieve blocks. Each element represents a odd number starting 71 * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13. 72 */ 73 extern char pattern[]; 74 extern int pattern_size; /* length of pattern array */ 75 76 #define MAX_LINE 255 /* max line allowed on stdin */ 77 78 void primes __P((ubig, ubig)); 79 char *read_num_buf __P((FILE *, char *)); 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 buf[MAX_LINE+1]; /* input buffer */ 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 /* 102 * parse args 103 */ 104 start = 0; 105 stop = BIG; 106 if (argc == 2) { 107 /* convert low and high args */ 108 if (read_num_buf(NULL, argv[0]) == NULL) 109 exit (1); 110 if (sscanf(argv[0], "%lu", &start) != 1) 111 err(1, "%s", argv[0]); 112 if (read_num_buf(NULL, argv[1]) == NULL) 113 exit (1); 114 if (sscanf(argv[1], "%lu", &stop) != 1) 115 err(1, "%s", argv[0]); 116 } else if (argc == 1) { 117 /* convert low arg */ 118 if (read_num_buf(NULL, argv[0]) == NULL) 119 exit (1); 120 if (sscanf(argv[0], "%lu", &start) != 1) 121 err(1, "%s", argv[0]); 122 } else { 123 /* read input until we get a good line */ 124 if (read_num_buf(stdin, buf) == NULL) 125 exit (1); 126 if (sscanf(buf, "%lu", &start) != 1) 127 err(1, "illegal value: %s", buf); 128 } 129 if (start > stop) 130 errx(1, "start value > stop value"); 131 primes(start, stop); 132 exit(0); 133 } 134 135 /* 136 * read_num_buf - read a number buffer from a stream 137 * 138 * Read a number on a line of the form: 139 * 140 * ^[ \t]*\(+?[0-9][0-9]\)*.*$ 141 * 142 * where ? is a 1-or-0 operator and the number is within \( \). 143 * 144 * If does not match the above pattern, it is ignored and a new 145 * line is read. If the number is too large or small, we will 146 * object and read a new line. 147 * 148 * We have to be very careful on how we check the magnitude of the 149 * input. We can not use numeric checks because of the need to 150 * check values against maximum numeric values. 151 * 152 * This routine will return a line containing a ascii number between 153 * 0 and BIG, or it will return NULL. 154 * 155 * If the stream is NULL then buf will be processed as if were 156 * a single line stream. 157 * 158 * returns: 159 * char * pointer to leading digit or + 160 * NULL EOF or error 161 */ 162 char * 163 read_num_buf(input, buf) 164 FILE *input; /* input stream or NULL */ 165 char *buf; /* input buffer */ 166 { 167 static char limit[MAX_LINE+1]; /* ascii value of BIG */ 168 static int limit_len; /* digit count of limit */ 169 int len; /* digits in input (excluding +/-) */ 170 char *s; /* line start marker */ 171 char *d; /* first digit, skip +/- */ 172 char *p; /* scan pointer */ 173 char *z; /* zero scan pointer */ 174 175 /* form the ascii value of BIG if needed */ 176 if (!isascii(limit[0]) || !isdigit(limit[0])) { 177 (void)sprintf(limit, "%lu", BIG); 178 limit_len = strlen(limit); 179 } 180 181 /* 182 * the search for a good line 183 */ 184 if (input != NULL && fgets(buf, MAX_LINE, input) == NULL) { 185 /* error or EOF */ 186 return (NULL); 187 } 188 do { 189 190 /* ignore leading whitespace */ 191 for (s=buf; *s && s < buf+MAX_LINE; ++s) { 192 if (!isascii(*s) || !isspace(*s)) { 193 break; 194 } 195 } 196 197 /* object if - */ 198 if (*s == '-') { 199 (void)fprintf(stderr, 200 "primes: numbers can't be negative.\n"); 201 continue; 202 } 203 204 /* skip over any leading + */ 205 if (*s == '+') { 206 d = s+1; 207 } else { 208 d = s; 209 } 210 211 /* note leading zeros */ 212 for (z=d; *z && z < buf+MAX_LINE; ++z) { 213 if (*z != '0') { 214 break; 215 } 216 } 217 218 /* scan for the first non-digit/non-plus/non-minus */ 219 for (p=d; *p && p < buf+MAX_LINE; ++p) { 220 if (!isascii(*p) || !isdigit(*p)) { 221 break; 222 } 223 } 224 225 /* ignore empty lines */ 226 if (p == d) { 227 continue; 228 } 229 *p = '\0'; 230 231 /* object if too many digits */ 232 len = strlen(z); 233 len = (len<=0) ? 1 : len; 234 /* accept if digit count is below limit */ 235 if (len < limit_len) { 236 /* we have good input */ 237 return (s); 238 239 /* 240 * reject very large numbers, carefully check against 241 * near limit numbers 242 */ 243 } else if (len > limit_len || strcmp(z, limit) > 0) { 244 errno = ERANGE; 245 warn("%s", z); 246 continue; 247 } 248 /* number is near limit, but is under it */ 249 return (s); 250 } while (input != NULL && fgets(buf, MAX_LINE, input) != NULL); 251 252 /* error or EOF */ 253 return (NULL); 254 } 255 256 /* 257 * primes - sieve and print primes from start up to and but not including stop 258 */ 259 void 260 primes(start, stop) 261 ubig start; /* where to start generating */ 262 ubig stop; /* don't generate at or above this value */ 263 { 264 register char *q; /* sieve spot */ 265 register ubig factor; /* index and factor */ 266 register char *tab_lim; /* the limit to sieve on the table */ 267 register ubig *p; /* prime table pointer */ 268 register ubig fact_lim; /* highest prime for current block */ 269 270 /* 271 * A number of systems can not convert double values into unsigned 272 * longs when the values are larger than the largest signed value. 273 * We don't have this problem, so we can go all the way to BIG. 274 */ 275 if (start < 3) { 276 start = (ubig)2; 277 } 278 if (stop < 3) { 279 stop = (ubig)2; 280 } 281 if (stop <= start) { 282 return; 283 } 284 285 /* 286 * be sure that the values are odd, or 2 287 */ 288 if (start != 2 && (start&0x1) == 0) { 289 ++start; 290 } 291 if (stop != 2 && (stop&0x1) == 0) { 292 ++stop; 293 } 294 295 /* 296 * quick list of primes <= pr_limit 297 */ 298 if (start <= *pr_limit) { 299 /* skip primes up to the start value */ 300 for (p = &prime[0], factor = prime[0]; 301 factor < stop && p <= pr_limit; factor = *(++p)) { 302 if (factor >= start) { 303 printf("%u\n", factor); 304 } 305 } 306 /* return early if we are done */ 307 if (p <= pr_limit) { 308 return; 309 } 310 start = *pr_limit+2; 311 } 312 313 /* 314 * we shall sieve a bytemap window, note primes and move the window 315 * upward until we pass the stop point 316 */ 317 while (start < stop) { 318 /* 319 * factor out 3, 5, 7, 11 and 13 320 */ 321 /* initial pattern copy */ 322 factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 323 memcpy(table, &pattern[factor], pattern_size-factor); 324 /* main block pattern copies */ 325 for (fact_lim=pattern_size-factor; 326 fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 327 memcpy(&table[fact_lim], pattern, pattern_size); 328 } 329 /* final block pattern copy */ 330 memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 331 332 /* 333 * sieve for primes 17 and higher 334 */ 335 /* note highest useful factor and sieve spot */ 336 if (stop-start > TABSIZE+TABSIZE) { 337 tab_lim = &table[TABSIZE]; /* sieve it all */ 338 fact_lim = (int)sqrt( 339 (double)(start)+TABSIZE+TABSIZE+1.0); 340 } else { 341 tab_lim = &table[(stop-start)/2]; /* partial sieve */ 342 fact_lim = (int)sqrt((double)(stop)+1.0); 343 } 344 /* sieve for factors >= 17 */ 345 factor = 17; /* 17 is first prime to use */ 346 p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 347 do { 348 /* determine the factor's initial sieve point */ 349 q = (char *)(start%factor); /* temp storage for mod */ 350 if ((int)q & 0x1) { 351 q = &table[(factor-(int)q)/2]; 352 } else { 353 q = &table[q ? factor-((int)q/2) : 0]; 354 } 355 /* sive for our current factor */ 356 for ( ; q < tab_lim; q += factor) { 357 *q = '\0'; /* sieve out a spot */ 358 } 359 } while ((factor=(ubig)(*(p++))) <= fact_lim); 360 361 /* 362 * print generated primes 363 */ 364 for (q = table; q < tab_lim; ++q, start+=2) { 365 if (*q) { 366 printf("%u\n", start); 367 } 368 } 369 } 370 } 371 372 void 373 usage() 374 { 375 (void)fprintf(stderr, "usage: primes [start stop]\n"); 376 exit(1); 377 } 378