security/nss/lib/freebl/mpi/utils/sieve.c

Wed, 31 Dec 2014 06:09:35 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Wed, 31 Dec 2014 06:09:35 +0100
changeset 0
6474c204b198
permissions
-rw-r--r--

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 }

mercurial