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