OLD | NEW |
1 /* | 1 /* |
2 * sieve.c | 2 * sieve.c |
3 * | 3 * |
4 * Finds prime numbers using the Sieve of Eratosthenes | 4 * Finds prime numbers using the Sieve of Eratosthenes |
5 * | 5 * |
6 * This implementation uses a bitmap to represent all odd integers in a | 6 * This implementation uses a bitmap to represent all odd integers in a |
7 * given range. We iterate over this bitmap, crossing off the | 7 * given range. We iterate over this bitmap, crossing off the |
8 * multiples of each prime we find. At the end, all the remaining set | 8 * multiples of each prime we find. At the end, all the remaining set |
9 * bits correspond to prime integers. | 9 * bits correspond to prime integers. |
10 * | 10 * |
(...skipping 10 matching lines...) Expand all Loading... |
21 * few more times. | 21 * few more times. |
22 * | 22 * |
23 * This Source Code Form is subject to the terms of the Mozilla Public | 23 * This Source Code Form is subject to the terms of the Mozilla Public |
24 * License, v. 2.0. If a copy of the MPL was not distributed with this | 24 * License, v. 2.0. If a copy of the MPL was not distributed with this |
25 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ | 25 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ |
26 | 26 |
27 #include <stdio.h> | 27 #include <stdio.h> |
28 #include <stdlib.h> | 28 #include <stdlib.h> |
29 #include <limits.h> | 29 #include <limits.h> |
30 | 30 |
31 typedef unsigned char byte; | 31 typedef unsigned char byte; |
32 | 32 |
33 typedef struct { | 33 typedef struct { |
34 int size; | 34 int size; |
35 byte *bits; | 35 byte *bits; |
36 long base; | 36 long base; |
37 int next; | 37 int next; |
38 int nbits; | 38 int nbits; |
39 } sieve; | 39 } sieve; |
40 | 40 |
41 void sieve_init(sieve *sp, long base, int nbits); | 41 void sieve_init(sieve *sp, long base, int nbits); |
42 void sieve_grow(sieve *sp, int nbits); | 42 void sieve_grow(sieve *sp, int nbits); |
43 long sieve_next(sieve *sp); | 43 long sieve_next(sieve *sp); |
44 void sieve_reset(sieve *sp, long base); | 44 void sieve_reset(sieve *sp, long base); |
45 void sieve_cross(sieve *sp, long val); | 45 void sieve_cross(sieve *sp, long val); |
46 void sieve_clear(sieve *sp); | 46 void sieve_clear(sieve *sp); |
47 | 47 |
48 #define S_ISSET(S, B) (((S)->bits[(B)/CHAR_BIT]>>((B)%CHAR_BIT))&1) | 48 #define S_ISSET(S, B) (((S)->bits[(B) / CHAR_BIT] >> ((B) % CHAR_BIT)) & 1) |
49 #define S_SET(S, B) ((S)->bits[(B)/CHAR_BIT]|=(1<<((B)%CHAR_BIT))) | 49 #define S_SET(S, B) ((S)->bits[(B) / CHAR_BIT] |= (1 << ((B) % CHAR_BIT))) |
50 #define S_CLR(S, B) ((S)->bits[(B)/CHAR_BIT]&=~(1<<((B)%CHAR_BIT))) | 50 #define S_CLR(S, B) ((S)->bits[(B) / CHAR_BIT] &= ~(1 << ((B) % CHAR_BIT))) |
51 #define S_VAL(S, B) ((S)->base+(2*(B))) | 51 #define S_VAL(S, B) ((S)->base + (2 * (B))) |
52 #define S_BIT(S, V) (((V)-((S)->base))/2) | 52 #define S_BIT(S, V) (((V) - ((S)->base)) / 2) |
53 | 53 |
54 int main(int argc, char *argv[]) | 54 int main(int argc, char *argv[]) { |
55 { | 55 sieve s; |
56 sieve s; | 56 long pr, *p; |
57 long pr, *p; | 57 int c, ix, cur = 0; |
58 int c, ix, cur = 0; | |
59 | 58 |
60 if(argc < 2) { | 59 if (argc < 2) { |
61 fprintf(stderr, "Usage: %s <width>\n", argv[0]); | 60 fprintf(stderr, "Usage: %s <width>\n", argv[0]); |
62 return 1; | 61 return 1; |
63 } | 62 } |
64 | 63 |
65 c = atoi(argv[1]); | 64 c = atoi(argv[1]); |
66 if(c < 0) c = -c; | 65 if (c < 0) c = -c; |
67 | 66 |
68 fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c); | 67 fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c); |
69 | 68 |
70 sieve_init(&s, 3, c); | 69 sieve_init(&s, 3, c); |
71 | 70 |
72 c = 0; | 71 c = 0; |
73 while((pr = sieve_next(&s)) > 0) { | 72 while ((pr = sieve_next(&s)) > 0) { |
74 ++c; | 73 ++c; |
75 } | 74 } |
76 | 75 |
77 p = calloc(c, sizeof(long)); | 76 p = calloc(c, sizeof(long)); |
78 if(!p) { | 77 if (!p) { |
79 fprintf(stderr, "%s: out of memory after first half\n", argv[0]); | 78 fprintf(stderr, "%s: out of memory after first half\n", argv[0]); |
80 sieve_clear(&s); | 79 sieve_clear(&s); |
81 exit(1); | 80 exit(1); |
82 } | 81 } |
83 | 82 |
84 fprintf(stderr, "%s: half done ... \n", argv[0]); | 83 fprintf(stderr, "%s: half done ... \n", argv[0]); |
85 | 84 |
86 for(ix = 0; ix < s.nbits; ix++) { | 85 for (ix = 0; ix < s.nbits; ix++) { |
87 if(S_ISSET(&s, ix)) { | 86 if (S_ISSET(&s, ix)) { |
88 p[cur] = S_VAL(&s, ix); | 87 p[cur] = S_VAL(&s, ix); |
89 printf("%ld\n", p[cur]); | 88 printf("%ld\n", p[cur]); |
90 ++cur; | 89 ++cur; |
91 } | 90 } |
92 } | 91 } |
93 | 92 |
94 sieve_reset(&s, p[cur - 1]); | 93 sieve_reset(&s, p[cur - 1]); |
95 fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur); | 94 fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur); |
96 for(ix = 0; ix < cur; ix++) { | 95 for (ix = 0; ix < cur; ix++) { |
97 sieve_cross(&s, p[ix]); | 96 sieve_cross(&s, p[ix]); |
98 if(!(ix % 1000)) | 97 if (!(ix % 1000)) fputc('.', stderr); |
99 fputc('.', stderr); | |
100 } | 98 } |
101 fputc('\n', stderr); | 99 fputc('\n', stderr); |
102 | 100 |
103 free(p); | 101 free(p); |
104 | 102 |
105 fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]); | 103 fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]); |
106 c = 0; | 104 c = 0; |
107 while((pr = sieve_next(&s)) > 0) { | 105 while ((pr = sieve_next(&s)) > 0) { |
108 ++c; | 106 ++c; |
109 } | 107 } |
110 | 108 |
111 fprintf(stderr, "%s: done!\n", argv[0]); | 109 fprintf(stderr, "%s: done!\n", argv[0]); |
112 for(ix = 0; ix < s.nbits; ix++) { | 110 for (ix = 0; ix < s.nbits; ix++) { |
113 if(S_ISSET(&s, ix)) { | 111 if (S_ISSET(&s, ix)) { |
114 printf("%ld\n", S_VAL(&s, ix)); | 112 printf("%ld\n", S_VAL(&s, ix)); |
115 } | 113 } |
116 } | 114 } |
117 | 115 |
118 sieve_clear(&s); | 116 sieve_clear(&s); |
119 | 117 |
120 return 0; | 118 return 0; |
121 } | 119 } |
122 | 120 |
123 void sieve_init(sieve *sp, long base, int nbits) | 121 void sieve_init(sieve *sp, long base, int nbits) { |
124 { | |
125 sp->size = (nbits / CHAR_BIT); | 122 sp->size = (nbits / CHAR_BIT); |
126 | 123 |
127 if(nbits % CHAR_BIT) | 124 if (nbits % CHAR_BIT) ++sp->size; |
128 ++sp->size; | |
129 | 125 |
130 sp->bits = calloc(sp->size, sizeof(byte)); | 126 sp->bits = calloc(sp->size, sizeof(byte)); |
131 memset(sp->bits, UCHAR_MAX, sp->size); | 127 memset(sp->bits, UCHAR_MAX, sp->size); |
132 if(!(base & 1)) | 128 if (!(base & 1)) ++base; |
133 ++base; | |
134 sp->base = base; | 129 sp->base = base; |
135 | 130 |
136 sp->next = 0; | 131 sp->next = 0; |
137 sp->nbits = sp->size * CHAR_BIT; | 132 sp->nbits = sp->size * CHAR_BIT; |
138 } | 133 } |
139 | 134 |
140 void sieve_grow(sieve *sp, int nbits) | 135 void sieve_grow(sieve *sp, int nbits) { |
141 { | 136 int ns = (nbits / CHAR_BIT); |
142 int ns = (nbits / CHAR_BIT); | |
143 | 137 |
144 if(nbits % CHAR_BIT) | 138 if (nbits % CHAR_BIT) ++ns; |
145 ++ns; | |
146 | 139 |
147 if(ns > sp->size) { | 140 if (ns > sp->size) { |
148 byte *tmp; | 141 byte *tmp; |
149 int ix; | 142 int ix; |
150 | 143 |
151 tmp = calloc(ns, sizeof(byte)); | 144 tmp = calloc(ns, sizeof(byte)); |
152 if(tmp == NULL) { | 145 if (tmp == NULL) { |
153 fprintf(stderr, "Error: out of memory in sieve_grow\n"); | 146 fprintf(stderr, "Error: out of memory in sieve_grow\n"); |
154 return; | 147 return; |
155 } | 148 } |
156 | 149 |
157 memcpy(tmp, sp->bits, sp->size); | 150 memcpy(tmp, sp->bits, sp->size); |
158 for(ix = sp->size; ix < ns; ix++) { | 151 for (ix = sp->size; ix < ns; ix++) { |
159 tmp[ix] = UCHAR_MAX; | 152 tmp[ix] = UCHAR_MAX; |
160 } | 153 } |
161 | 154 |
162 free(sp->bits); | 155 free(sp->bits); |
163 sp->bits = tmp; | 156 sp->bits = tmp; |
164 sp->size = ns; | 157 sp->size = ns; |
165 | 158 |
166 sp->nbits = sp->size * CHAR_BIT; | 159 sp->nbits = sp->size * CHAR_BIT; |
167 } | 160 } |
168 } | 161 } |
169 | 162 |
170 long sieve_next(sieve *sp) | 163 long sieve_next(sieve *sp) { |
171 { | |
172 long out; | 164 long out; |
173 int ix = 0; | 165 int ix = 0; |
174 long val; | 166 long val; |
175 | 167 |
176 if(sp->next > sp->nbits) | 168 if (sp->next > sp->nbits) return -1; |
177 return -1; | |
178 | 169 |
179 out = S_VAL(sp, sp->next); | 170 out = S_VAL(sp, sp->next); |
180 #ifdef DEBUG | 171 #ifdef DEBUG |
181 fprintf(stderr, "Sieving %ld\n", out); | 172 fprintf(stderr, "Sieving %ld\n", out); |
182 #endif | 173 #endif |
183 | 174 |
184 /* Sieve out all multiples of the current prime */ | 175 /* Sieve out all multiples of the current prime */ |
185 val = out; | 176 val = out; |
186 while(ix < sp->nbits) { | 177 while (ix < sp->nbits) { |
187 val += out; | 178 val += out; |
188 ix = S_BIT(sp, val); | 179 ix = S_BIT(sp, val); |
189 if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */ | 180 if ((val & 1) && ix < sp->nbits) {/* && S_ISSET(sp, ix)) { */ |
190 S_CLR(sp, ix); | 181 S_CLR(sp, ix); |
191 #ifdef DEBUG | 182 #ifdef DEBUG |
192 fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix); | 183 fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix); |
193 #endif | 184 #endif |
194 } | 185 } |
195 } | 186 } |
196 | 187 |
197 /* Scan ahead to the next prime */ | 188 /* Scan ahead to the next prime */ |
198 ++sp->next; | 189 ++sp->next; |
199 while(sp->next < sp->nbits && !S_ISSET(sp, sp->next)) | 190 while (sp->next < sp->nbits && !S_ISSET(sp, sp->next)) ++sp->next; |
200 ++sp->next; | |
201 | 191 |
202 return out; | 192 return out; |
203 } | 193 } |
204 | 194 |
205 void sieve_cross(sieve *sp, long val) | 195 void sieve_cross(sieve *sp, long val) { |
206 { | 196 int ix = 0; |
207 int ix = 0; | |
208 long cur = val; | 197 long cur = val; |
209 | 198 |
210 while(cur < sp->base) | 199 while (cur < sp->base) cur += val; |
211 cur += val; | |
212 | 200 |
213 ix = S_BIT(sp, cur); | 201 ix = S_BIT(sp, cur); |
214 while(ix < sp->nbits) { | 202 while (ix < sp->nbits) { |
215 if(cur & 1)· | 203 if (cur & 1) S_CLR(sp, ix); |
216 S_CLR(sp, ix); | |
217 cur += val; | 204 cur += val; |
218 ix = S_BIT(sp, cur); | 205 ix = S_BIT(sp, cur); |
219 } | 206 } |
220 } | 207 } |
221 | 208 |
222 void sieve_reset(sieve *sp, long base) | 209 void sieve_reset(sieve *sp, long base) { |
223 { | |
224 memset(sp->bits, UCHAR_MAX, sp->size); | 210 memset(sp->bits, UCHAR_MAX, sp->size); |
225 sp->base = base; | 211 sp->base = base; |
226 sp->next = 0; | 212 sp->next = 0; |
227 } | 213 } |
228 | 214 |
229 void sieve_clear(sieve *sp) | 215 void sieve_clear(sieve *sp) { |
230 { | 216 if (sp->bits) free(sp->bits); |
231 if(sp->bits)· | |
232 free(sp->bits); | |
233 | 217 |
234 sp->bits = NULL; | 218 sp->bits = NULL; |
235 } | 219 } |
OLD | NEW |