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) 5/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: }