diff options
Diffstat (limited to 'mozilla/security/nss/lib/freebl/mpi/utils/sieve.c')
| -rw-r--r-- | mozilla/security/nss/lib/freebl/mpi/utils/sieve.c | 268 |
1 files changed, 268 insertions, 0 deletions
diff --git a/mozilla/security/nss/lib/freebl/mpi/utils/sieve.c b/mozilla/security/nss/lib/freebl/mpi/utils/sieve.c new file mode 100644 index 0000000..3c8d18a --- /dev/null +++ b/mozilla/security/nss/lib/freebl/mpi/utils/sieve.c @@ -0,0 +1,268 @@ +/* + * sieve.c + * + * Finds prime numbers using the Sieve of Eratosthenes + * + * This implementation uses a bitmap to represent all odd integers in a + * given range. We iterate over this bitmap, crossing off the + * multiples of each prime we find. At the end, all the remaining set + * bits correspond to prime integers. + * + * Here, we make two passes -- once we have generated a sieve-ful of + * primes, we copy them out, reset the sieve using the highest + * generated prime from the first pass as a base. Then we cross out + * all the multiples of all the primes we found the first time through, + * and re-sieve. In this way, we get double use of the memory we + * allocated for the sieve the first time though. Since we also + * implicitly ignore multiples of 2, this amounts to 4 times the + * values. + * + * This could (and probably will) be generalized to re-use the sieve a + * few more times. + * + * ***** BEGIN LICENSE BLOCK ***** + * Version: MPL 1.1/GPL 2.0/LGPL 2.1 + * + * The contents of this file are subject to the Mozilla Public License Version + * 1.1 (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * Software distributed under the License is distributed on an "AS IS" basis, + * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License + * for the specific language governing rights and limitations under the + * License. + * + * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library. + * + * The Initial Developer of the Original Code is + * Michael J. Fromberger. + * Portions created by the Initial Developer are Copyright (C) 1998 + * the Initial Developer. All Rights Reserved. + * + * Contributor(s): + * + * Alternatively, the contents of this file may be used under the terms of + * either the GNU General Public License Version 2 or later (the "GPL"), or + * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), + * in which case the provisions of the GPL or the LGPL are applicable instead + * of those above. If you wish to allow use of your version of this file only + * under the terms of either the GPL or the LGPL, and not to allow others to + * use your version of this file under the terms of the MPL, indicate your + * decision by deleting the provisions above and replace them with the notice + * and other provisions required by the GPL or the LGPL. If you do not delete + * the provisions above, a recipient may use your version of this file under + * the terms of any one of the MPL, the GPL or the LGPL. + * + * ***** END LICENSE BLOCK ***** */ +/* $Id: sieve.c,v 1.3 2004/04/27 23:04:37 gerv%gerv.net Exp $ */ + +#include <stdio.h> +#include <stdlib.h> +#include <limits.h> + +typedef unsigned char byte; + +typedef struct { + int size; + byte *bits; + long base; + int next; + int nbits; +} sieve; + +void sieve_init(sieve *sp, long base, int nbits); +void sieve_grow(sieve *sp, int nbits); +long sieve_next(sieve *sp); +void sieve_reset(sieve *sp, long base); +void sieve_cross(sieve *sp, long val); +void sieve_clear(sieve *sp); + +#define S_ISSET(S, B) (((S)->bits[(B)/CHAR_BIT]>>((B)%CHAR_BIT))&1) +#define S_SET(S, B) ((S)->bits[(B)/CHAR_BIT]|=(1<<((B)%CHAR_BIT))) +#define S_CLR(S, B) ((S)->bits[(B)/CHAR_BIT]&=~(1<<((B)%CHAR_BIT))) +#define S_VAL(S, B) ((S)->base+(2*(B))) +#define S_BIT(S, V) (((V)-((S)->base))/2) + +int main(int argc, char *argv[]) +{ + sieve s; + long pr, *p; + int c, ix, cur = 0; + + if(argc < 2) { + fprintf(stderr, "Usage: %s <width>\n", argv[0]); + return 1; + } + + c = atoi(argv[1]); + if(c < 0) c = -c; + + fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c); + + sieve_init(&s, 3, c); + + c = 0; + while((pr = sieve_next(&s)) > 0) { + ++c; + } + + p = calloc(c, sizeof(long)); + if(!p) { + fprintf(stderr, "%s: out of memory after first half\n", argv[0]); + sieve_clear(&s); + exit(1); + } + + fprintf(stderr, "%s: half done ... \n", argv[0]); + + for(ix = 0; ix < s.nbits; ix++) { + if(S_ISSET(&s, ix)) { + p[cur] = S_VAL(&s, ix); + printf("%ld\n", p[cur]); + ++cur; + } + } + + sieve_reset(&s, p[cur - 1]); + fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur); + for(ix = 0; ix < cur; ix++) { + sieve_cross(&s, p[ix]); + if(!(ix % 1000)) + fputc('.', stderr); + } + fputc('\n', stderr); + + free(p); + + fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]); + c = 0; + while((pr = sieve_next(&s)) > 0) { + ++c; + } + + fprintf(stderr, "%s: done!\n", argv[0]); + for(ix = 0; ix < s.nbits; ix++) { + if(S_ISSET(&s, ix)) { + printf("%ld\n", S_VAL(&s, ix)); + } + } + + sieve_clear(&s); + + return 0; +} + +void sieve_init(sieve *sp, long base, int nbits) +{ + sp->size = (nbits / CHAR_BIT); + + if(nbits % CHAR_BIT) + ++sp->size; + + sp->bits = calloc(sp->size, sizeof(byte)); + memset(sp->bits, UCHAR_MAX, sp->size); + if(!(base & 1)) + ++base; + sp->base = base; + + sp->next = 0; + sp->nbits = sp->size * CHAR_BIT; +} + +void sieve_grow(sieve *sp, int nbits) +{ + int ns = (nbits / CHAR_BIT); + + if(nbits % CHAR_BIT) + ++ns; + + if(ns > sp->size) { + byte *tmp; + int ix; + + tmp = calloc(ns, sizeof(byte)); + if(tmp == NULL) { + fprintf(stderr, "Error: out of memory in sieve_grow\n"); + return; + } + + memcpy(tmp, sp->bits, sp->size); + for(ix = sp->size; ix < ns; ix++) { + tmp[ix] = UCHAR_MAX; + } + + free(sp->bits); + sp->bits = tmp; + sp->size = ns; + + sp->nbits = sp->size * CHAR_BIT; + } +} + +long sieve_next(sieve *sp) +{ + long out; + int ix = 0; + long val; + + if(sp->next > sp->nbits) + return -1; + + out = S_VAL(sp, sp->next); +#ifdef DEBUG + fprintf(stderr, "Sieving %ld\n", out); +#endif + + /* Sieve out all multiples of the current prime */ + val = out; + while(ix < sp->nbits) { + val += out; + ix = S_BIT(sp, val); + if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */ + S_CLR(sp, ix); +#ifdef DEBUG + fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix); +#endif + } + } + + /* Scan ahead to the next prime */ + ++sp->next; + while(sp->next < sp->nbits && !S_ISSET(sp, sp->next)) + ++sp->next; + + return out; +} + +void sieve_cross(sieve *sp, long val) +{ + int ix = 0; + long cur = val; + + while(cur < sp->base) + cur += val; + + ix = S_BIT(sp, cur); + while(ix < sp->nbits) { + if(cur & 1) + S_CLR(sp, ix); + cur += val; + ix = S_BIT(sp, cur); + } +} + +void sieve_reset(sieve *sp, long base) +{ + memset(sp->bits, UCHAR_MAX, sp->size); + sp->base = base; + sp->next = 0; +} + +void sieve_clear(sieve *sp) +{ + if(sp->bits) + free(sp->bits); + + sp->bits = NULL; +} |