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