The Tor Browser / file comparison
comparison: security/nss/lib/freebl/mpi/utils/sieve.c
security/nss/lib/freebl/mpi/utils/sieve.c
- changeset 0
- 6474c204b198
equal
deleted
inserted
replaced
| -1:000000000000 | 0:e4ab40b753b2 |
|---|---|
| 1 /* | |
| 2 * sieve.c | |
| 3 * | |
| 4 * Finds prime numbers using the Sieve of Eratosthenes | |
| 5 * | |
| 6 * This implementation uses a bitmap to represent all odd integers in a | |
| 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 | |
| 9 * bits correspond to prime integers. | |
| 10 * | |
| 11 * Here, we make two passes -- once we have generated a sieve-ful of | |
| 12 * primes, we copy them out, reset the sieve using the highest | |
| 13 * generated prime from the first pass as a base. Then we cross out | |
| 14 * all the multiples of all the primes we found the first time through, | |
| 15 * and re-sieve. In this way, we get double use of the memory we | |
| 16 * allocated for the sieve the first time though. Since we also | |
| 17 * implicitly ignore multiples of 2, this amounts to 4 times the | |
| 18 * values. | |
| 19 * | |
| 20 * This could (and probably will) be generalized to re-use the sieve a | |
| 21 * few more times. | |
| 22 * | |
| 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 | |
| 25 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ | |
| 26 | |
| 27 #include <stdio.h> | |
| 28 #include <stdlib.h> | |
| 29 #include <limits.h> | |
| 30 | |
| 31 typedef unsigned char byte; | |
| 32 | |
| 33 typedef struct { | |
| 34 int size; | |
| 35 byte *bits; | |
| 36 long base; | |
| 37 int next; | |
| 38 int nbits; | |
| 39 } sieve; | |
| 40 | |
| 41 void sieve_init(sieve *sp, long base, int nbits); | |
| 42 void sieve_grow(sieve *sp, int nbits); | |
| 43 long sieve_next(sieve *sp); | |
| 44 void sieve_reset(sieve *sp, long base); | |
| 45 void sieve_cross(sieve *sp, long val); | |
| 46 void sieve_clear(sieve *sp); | |
| 47 | |
| 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))) | |
| 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))) | |
| 52 #define S_BIT(S, V) (((V)-((S)->base))/2) | |
| 53 | |
| 54 int main(int argc, char *argv[]) | |
| 55 { | |
| 56 sieve s; | |
| 57 long pr, *p; | |
| 58 int c, ix, cur = 0; | |
| 59 | |
| 60 if(argc < 2) { | |
| 61 fprintf(stderr, "Usage: %s <width>\n", argv[0]); | |
| 62 return 1; | |
| 63 } | |
| 64 | |
| 65 c = atoi(argv[1]); | |
| 66 if(c < 0) c = -c; | |
| 67 | |
| 68 fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c); | |
| 69 | |
| 70 sieve_init(&s, 3, c); | |
| 71 | |
| 72 c = 0; | |
| 73 while((pr = sieve_next(&s)) > 0) { | |
| 74 ++c; | |
| 75 } | |
| 76 | |
| 77 p = calloc(c, sizeof(long)); | |
| 78 if(!p) { | |
| 79 fprintf(stderr, "%s: out of memory after first half\n", argv[0]); | |
| 80 sieve_clear(&s); | |
| 81 exit(1); | |
| 82 } | |
| 83 | |
| 84 fprintf(stderr, "%s: half done ... \n", argv[0]); | |
| 85 | |
| 86 for(ix = 0; ix < s.nbits; ix++) { | |
| 87 if(S_ISSET(&s, ix)) { | |
| 88 p[cur] = S_VAL(&s, ix); | |
| 89 printf("%ld\n", p[cur]); | |
| 90 ++cur; | |
| 91 } | |
| 92 } | |
| 93 | |
| 94 sieve_reset(&s, p[cur - 1]); | |
| 95 fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur); | |
| 96 for(ix = 0; ix < cur; ix++) { | |
| 97 sieve_cross(&s, p[ix]); | |
| 98 if(!(ix % 1000)) | |
| 99 fputc('.', stderr); | |
| 100 } | |
| 101 fputc('\n', stderr); | |
| 102 | |
| 103 free(p); | |
| 104 | |
| 105 fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]); | |
| 106 c = 0; | |
| 107 while((pr = sieve_next(&s)) > 0) { | |
| 108 ++c; | |
| 109 } | |
| 110 | |
| 111 fprintf(stderr, "%s: done!\n", argv[0]); | |
| 112 for(ix = 0; ix < s.nbits; ix++) { | |
| 113 if(S_ISSET(&s, ix)) { | |
| 114 printf("%ld\n", S_VAL(&s, ix)); | |
| 115 } | |
| 116 } | |
| 117 | |
| 118 sieve_clear(&s); | |
| 119 | |
| 120 return 0; | |
| 121 } | |
| 122 | |
| 123 void sieve_init(sieve *sp, long base, int nbits) | |
| 124 { | |
| 125 sp->size = (nbits / CHAR_BIT); | |
| 126 | |
| 127 if(nbits % CHAR_BIT) | |
| 128 ++sp->size; | |
| 129 | |
| 130 sp->bits = calloc(sp->size, sizeof(byte)); | |
| 131 memset(sp->bits, UCHAR_MAX, sp->size); | |
| 132 if(!(base & 1)) | |
| 133 ++base; | |
| 134 sp->base = base; | |
| 135 | |
| 136 sp->next = 0; | |
| 137 sp->nbits = sp->size * CHAR_BIT; | |
| 138 } | |
| 139 | |
| 140 void sieve_grow(sieve *sp, int nbits) | |
| 141 { | |
| 142 int ns = (nbits / CHAR_BIT); | |
| 143 | |
| 144 if(nbits % CHAR_BIT) | |
| 145 ++ns; | |
| 146 | |
| 147 if(ns > sp->size) { | |
| 148 byte *tmp; | |
| 149 int ix; | |
| 150 | |
| 151 tmp = calloc(ns, sizeof(byte)); | |
| 152 if(tmp == NULL) { | |
| 153 fprintf(stderr, "Error: out of memory in sieve_grow\n"); | |
| 154 return; | |
| 155 } | |
| 156 | |
| 157 memcpy(tmp, sp->bits, sp->size); | |
| 158 for(ix = sp->size; ix < ns; ix++) { | |
| 159 tmp[ix] = UCHAR_MAX; | |
| 160 } | |
| 161 | |
| 162 free(sp->bits); | |
| 163 sp->bits = tmp; | |
| 164 sp->size = ns; | |
| 165 | |
| 166 sp->nbits = sp->size * CHAR_BIT; | |
| 167 } | |
| 168 } | |
| 169 | |
| 170 long sieve_next(sieve *sp) | |
| 171 { | |
| 172 long out; | |
| 173 int ix = 0; | |
| 174 long val; | |
| 175 | |
| 176 if(sp->next > sp->nbits) | |
| 177 return -1; | |
| 178 | |
| 179 out = S_VAL(sp, sp->next); | |
| 180 #ifdef DEBUG | |
| 181 fprintf(stderr, "Sieving %ld\n", out); | |
| 182 #endif | |
| 183 | |
| 184 /* Sieve out all multiples of the current prime */ | |
| 185 val = out; | |
| 186 while(ix < sp->nbits) { | |
| 187 val += out; | |
| 188 ix = S_BIT(sp, val); | |
| 189 if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */ | |
| 190 S_CLR(sp, ix); | |
| 191 #ifdef DEBUG | |
| 192 fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix); | |
| 193 #endif | |
| 194 } | |
| 195 } | |
| 196 | |
| 197 /* Scan ahead to the next prime */ | |
| 198 ++sp->next; | |
| 199 while(sp->next < sp->nbits && !S_ISSET(sp, sp->next)) | |
| 200 ++sp->next; | |
| 201 | |
| 202 return out; | |
| 203 } | |
| 204 | |
| 205 void sieve_cross(sieve *sp, long val) | |
| 206 { | |
| 207 int ix = 0; | |
| 208 long cur = val; | |
| 209 | |
| 210 while(cur < sp->base) | |
| 211 cur += val; | |
| 212 | |
| 213 ix = S_BIT(sp, cur); | |
| 214 while(ix < sp->nbits) { | |
| 215 if(cur & 1) | |
| 216 S_CLR(sp, ix); | |
| 217 cur += val; | |
| 218 ix = S_BIT(sp, cur); | |
| 219 } | |
| 220 } | |
| 221 | |
| 222 void sieve_reset(sieve *sp, long base) | |
| 223 { | |
| 224 memset(sp->bits, UCHAR_MAX, sp->size); | |
| 225 sp->base = base; | |
| 226 sp->next = 0; | |
| 227 } | |
| 228 | |
| 229 void sieve_clear(sieve *sp) | |
| 230 { | |
| 231 if(sp->bits) | |
| 232 free(sp->bits); | |
| 233 | |
| 234 sp->bits = NULL; | |
| 235 } |
