| 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) %G%"; |
| 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 | } |