1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/security/nss/lib/freebl/mpi/utils/sieve.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,235 @@
1.4 +/*
1.5 + * sieve.c
1.6 + *
1.7 + * Finds prime numbers using the Sieve of Eratosthenes
1.8 + *
1.9 + * This implementation uses a bitmap to represent all odd integers in a
1.10 + * given range. We iterate over this bitmap, crossing off the
1.11 + * multiples of each prime we find. At the end, all the remaining set
1.12 + * bits correspond to prime integers.
1.13 + *
1.14 + * Here, we make two passes -- once we have generated a sieve-ful of
1.15 + * primes, we copy them out, reset the sieve using the highest
1.16 + * generated prime from the first pass as a base. Then we cross out
1.17 + * all the multiples of all the primes we found the first time through,
1.18 + * and re-sieve. In this way, we get double use of the memory we
1.19 + * allocated for the sieve the first time though. Since we also
1.20 + * implicitly ignore multiples of 2, this amounts to 4 times the
1.21 + * values.
1.22 + *
1.23 + * This could (and probably will) be generalized to re-use the sieve a
1.24 + * few more times.
1.25 + *
1.26 + * This Source Code Form is subject to the terms of the Mozilla Public
1.27 + * License, v. 2.0. If a copy of the MPL was not distributed with this
1.28 + * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
1.29 +
1.30 +#include <stdio.h>
1.31 +#include <stdlib.h>
1.32 +#include <limits.h>
1.33 +
1.34 +typedef unsigned char byte;
1.35 +
1.36 +typedef struct {
1.37 + int size;
1.38 + byte *bits;
1.39 + long base;
1.40 + int next;
1.41 + int nbits;
1.42 +} sieve;
1.43 +
1.44 +void sieve_init(sieve *sp, long base, int nbits);
1.45 +void sieve_grow(sieve *sp, int nbits);
1.46 +long sieve_next(sieve *sp);
1.47 +void sieve_reset(sieve *sp, long base);
1.48 +void sieve_cross(sieve *sp, long val);
1.49 +void sieve_clear(sieve *sp);
1.50 +
1.51 +#define S_ISSET(S, B) (((S)->bits[(B)/CHAR_BIT]>>((B)%CHAR_BIT))&1)
1.52 +#define S_SET(S, B) ((S)->bits[(B)/CHAR_BIT]|=(1<<((B)%CHAR_BIT)))
1.53 +#define S_CLR(S, B) ((S)->bits[(B)/CHAR_BIT]&=~(1<<((B)%CHAR_BIT)))
1.54 +#define S_VAL(S, B) ((S)->base+(2*(B)))
1.55 +#define S_BIT(S, V) (((V)-((S)->base))/2)
1.56 +
1.57 +int main(int argc, char *argv[])
1.58 +{
1.59 + sieve s;
1.60 + long pr, *p;
1.61 + int c, ix, cur = 0;
1.62 +
1.63 + if(argc < 2) {
1.64 + fprintf(stderr, "Usage: %s <width>\n", argv[0]);
1.65 + return 1;
1.66 + }
1.67 +
1.68 + c = atoi(argv[1]);
1.69 + if(c < 0) c = -c;
1.70 +
1.71 + fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c);
1.72 +
1.73 + sieve_init(&s, 3, c);
1.74 +
1.75 + c = 0;
1.76 + while((pr = sieve_next(&s)) > 0) {
1.77 + ++c;
1.78 + }
1.79 +
1.80 + p = calloc(c, sizeof(long));
1.81 + if(!p) {
1.82 + fprintf(stderr, "%s: out of memory after first half\n", argv[0]);
1.83 + sieve_clear(&s);
1.84 + exit(1);
1.85 + }
1.86 +
1.87 + fprintf(stderr, "%s: half done ... \n", argv[0]);
1.88 +
1.89 + for(ix = 0; ix < s.nbits; ix++) {
1.90 + if(S_ISSET(&s, ix)) {
1.91 + p[cur] = S_VAL(&s, ix);
1.92 + printf("%ld\n", p[cur]);
1.93 + ++cur;
1.94 + }
1.95 + }
1.96 +
1.97 + sieve_reset(&s, p[cur - 1]);
1.98 + fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur);
1.99 + for(ix = 0; ix < cur; ix++) {
1.100 + sieve_cross(&s, p[ix]);
1.101 + if(!(ix % 1000))
1.102 + fputc('.', stderr);
1.103 + }
1.104 + fputc('\n', stderr);
1.105 +
1.106 + free(p);
1.107 +
1.108 + fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]);
1.109 + c = 0;
1.110 + while((pr = sieve_next(&s)) > 0) {
1.111 + ++c;
1.112 + }
1.113 +
1.114 + fprintf(stderr, "%s: done!\n", argv[0]);
1.115 + for(ix = 0; ix < s.nbits; ix++) {
1.116 + if(S_ISSET(&s, ix)) {
1.117 + printf("%ld\n", S_VAL(&s, ix));
1.118 + }
1.119 + }
1.120 +
1.121 + sieve_clear(&s);
1.122 +
1.123 + return 0;
1.124 +}
1.125 +
1.126 +void sieve_init(sieve *sp, long base, int nbits)
1.127 +{
1.128 + sp->size = (nbits / CHAR_BIT);
1.129 +
1.130 + if(nbits % CHAR_BIT)
1.131 + ++sp->size;
1.132 +
1.133 + sp->bits = calloc(sp->size, sizeof(byte));
1.134 + memset(sp->bits, UCHAR_MAX, sp->size);
1.135 + if(!(base & 1))
1.136 + ++base;
1.137 + sp->base = base;
1.138 +
1.139 + sp->next = 0;
1.140 + sp->nbits = sp->size * CHAR_BIT;
1.141 +}
1.142 +
1.143 +void sieve_grow(sieve *sp, int nbits)
1.144 +{
1.145 + int ns = (nbits / CHAR_BIT);
1.146 +
1.147 + if(nbits % CHAR_BIT)
1.148 + ++ns;
1.149 +
1.150 + if(ns > sp->size) {
1.151 + byte *tmp;
1.152 + int ix;
1.153 +
1.154 + tmp = calloc(ns, sizeof(byte));
1.155 + if(tmp == NULL) {
1.156 + fprintf(stderr, "Error: out of memory in sieve_grow\n");
1.157 + return;
1.158 + }
1.159 +
1.160 + memcpy(tmp, sp->bits, sp->size);
1.161 + for(ix = sp->size; ix < ns; ix++) {
1.162 + tmp[ix] = UCHAR_MAX;
1.163 + }
1.164 +
1.165 + free(sp->bits);
1.166 + sp->bits = tmp;
1.167 + sp->size = ns;
1.168 +
1.169 + sp->nbits = sp->size * CHAR_BIT;
1.170 + }
1.171 +}
1.172 +
1.173 +long sieve_next(sieve *sp)
1.174 +{
1.175 + long out;
1.176 + int ix = 0;
1.177 + long val;
1.178 +
1.179 + if(sp->next > sp->nbits)
1.180 + return -1;
1.181 +
1.182 + out = S_VAL(sp, sp->next);
1.183 +#ifdef DEBUG
1.184 + fprintf(stderr, "Sieving %ld\n", out);
1.185 +#endif
1.186 +
1.187 + /* Sieve out all multiples of the current prime */
1.188 + val = out;
1.189 + while(ix < sp->nbits) {
1.190 + val += out;
1.191 + ix = S_BIT(sp, val);
1.192 + if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */
1.193 + S_CLR(sp, ix);
1.194 +#ifdef DEBUG
1.195 + fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix);
1.196 +#endif
1.197 + }
1.198 + }
1.199 +
1.200 + /* Scan ahead to the next prime */
1.201 + ++sp->next;
1.202 + while(sp->next < sp->nbits && !S_ISSET(sp, sp->next))
1.203 + ++sp->next;
1.204 +
1.205 + return out;
1.206 +}
1.207 +
1.208 +void sieve_cross(sieve *sp, long val)
1.209 +{
1.210 + int ix = 0;
1.211 + long cur = val;
1.212 +
1.213 + while(cur < sp->base)
1.214 + cur += val;
1.215 +
1.216 + ix = S_BIT(sp, cur);
1.217 + while(ix < sp->nbits) {
1.218 + if(cur & 1)
1.219 + S_CLR(sp, ix);
1.220 + cur += val;
1.221 + ix = S_BIT(sp, cur);
1.222 + }
1.223 +}
1.224 +
1.225 +void sieve_reset(sieve *sp, long base)
1.226 +{
1.227 + memset(sp->bits, UCHAR_MAX, sp->size);
1.228 + sp->base = base;
1.229 + sp->next = 0;
1.230 +}
1.231 +
1.232 +void sieve_clear(sieve *sp)
1.233 +{
1.234 + if(sp->bits)
1.235 + free(sp->bits);
1.236 +
1.237 + sp->bits = NULL;
1.238 +}
