1 #ifndef lint 2 static char sccsid[] = "@(#)primes.c 4.1 (Wollongong) 06/13/83"; 3 #endif 4 5 /* 6 * primes [ number ] 7 * 8 * Print all primes greater than argument (or number read from stdin). 9 * 10 * A free translation of 'primes.s' 11 * 12 */ 13 14 #include <stdio.h> 15 #include <math.h> 16 17 #define TABSIZE 1000 /* size of sieve table */ 18 #define BIG 4294967296. /* largest unsigned int */ 19 20 char table[TABSIZE]; /* table for sieve of Eratosthenes */ 21 int tabbits = 8*TABSIZE; /* number of bits in table */ 22 23 float fstart; 24 unsigned start; /* lowest number to test for prime */ 25 char bittab[] = { /* bit positions (to save shifting) */ 26 01, 02, 04, 010, 020, 040, 0100, 0200 27 }; 28 29 unsigned pt[] = { /* primes < 100 */ 30 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 31 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 32 }; 33 34 unsigned factab[] = { /* difference between succesive trial factors */ 35 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 36 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 37 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2 38 }; 39 40 main(argc, argv) 41 int argc; 42 char **argv; 43 { 44 register unsigned *fp; 45 register char *p; 46 register int i; 47 unsigned quot; 48 unsigned factor, v; 49 50 if (argc >= 2) { /* get starting no. from arg */ 51 if (sscanf(argv[1], "%f", &fstart) != 1 52 || fstart < 0.0 || fstart >= BIG) { 53 ouch(); 54 exit(1); 55 } 56 } else { /* get starting no. from stdin */ 57 while ((i = scanf("%f", &fstart)) != 1 58 || fstart < 0.0 || fstart >= BIG) { 59 if (i == EOF) 60 exit(1); 61 ouch(); 62 } 63 } 64 start = (unsigned)fstart; 65 66 /* 67 * Quick list of primes < 100 68 */ 69 if (start <= 97) { 70 for (fp = pt; *fp < start; fp++) 71 ; 72 do 73 printf("%u\n", *fp); 74 while (++fp < &pt[sizeof(pt) / sizeof(*pt)]); 75 start = 100; 76 } 77 quot = start/2; 78 start = quot * 2 + 1; 79 80 /* 81 * Loop forever: 82 */ 83 for (;;) { 84 /* 85 * Generate primes via sieve of Eratosthenes 86 */ 87 for (p = table; p < &table[TABSIZE]; p++) /* clear sieve */ 88 *p = '\0'; 89 v = (unsigned)sqrt((float)(start + tabbits)); /* highest useful factor */ 90 sieve(3); 91 sieve(5); 92 sieve(7); 93 factor = 11; 94 fp = &factab[1]; 95 do { 96 sieve(factor); 97 factor += *fp; 98 if (++fp >= &factab[sizeof(factab) / sizeof(*factab)]) 99 fp = factab; 100 } while (factor <= v); 101 /* 102 * Print generated primes 103 */ 104 for (i = 0; i < 8*TABSIZE; i += 2) { 105 if ((table[i>>3] & bittab[i&07]) == 0) 106 printf("%u\n", start); 107 start += 2; 108 } 109 } 110 } 111 112 /* 113 * Insert all multiples of given factor into the sieve 114 */ 115 sieve(factor) 116 unsigned factor; 117 { 118 register int i; 119 unsigned off; 120 unsigned quot; 121 122 quot = start / factor; 123 off = (quot * factor) - start; 124 if ((int)off < 0) 125 off += factor; 126 while (off < tabbits ) { 127 i = (int)off; 128 table[i>>3] |= bittab[i&07]; 129 off += factor; 130 } 131 } 132 133 /* 134 * Error message 135 */ 136 ouch() 137 { 138 fprintf(stderr, "Ouch.\n"); 139 } 140