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