michael@0: /*
michael@0:  * sieve.c
michael@0:  *
michael@0:  * Finds prime numbers using the Sieve of Eratosthenes
michael@0:  *
michael@0:  * This implementation uses a bitmap to represent all odd integers in a
michael@0:  * given range.  We iterate over this bitmap, crossing off the
michael@0:  * multiples of each prime we find.  At the end, all the remaining set
michael@0:  * bits correspond to prime integers.
michael@0:  *
michael@0:  * Here, we make two passes -- once we have generated a sieve-ful of
michael@0:  * primes, we copy them out, reset the sieve using the highest
michael@0:  * generated prime from the first pass as a base.  Then we cross out
michael@0:  * all the multiples of all the primes we found the first time through,
michael@0:  * and re-sieve.  In this way, we get double use of the memory we
michael@0:  * allocated for the sieve the first time though.  Since we also
michael@0:  * implicitly ignore multiples of 2, this amounts to 4 times the
michael@0:  * values.
michael@0:  *
michael@0:  * This could (and probably will) be generalized to re-use the sieve a
michael@0:  * few more times.
michael@0:  *
michael@0:  * This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  * License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
michael@0: 
michael@0: #include <stdio.h>
michael@0: #include <stdlib.h>
michael@0: #include <limits.h>
michael@0: 
michael@0: typedef unsigned char  byte;
michael@0: 
michael@0: typedef struct {
michael@0:   int   size;
michael@0:   byte *bits;
michael@0:   long  base;
michael@0:   int   next;
michael@0:   int   nbits;
michael@0: } sieve;
michael@0: 
michael@0: void sieve_init(sieve *sp, long base, int nbits);
michael@0: void sieve_grow(sieve *sp, int nbits);
michael@0: long sieve_next(sieve *sp);
michael@0: void sieve_reset(sieve *sp, long base);
michael@0: void sieve_cross(sieve *sp, long val);
michael@0: void sieve_clear(sieve *sp);
michael@0: 
michael@0: #define S_ISSET(S, B)  (((S)->bits[(B)/CHAR_BIT]>>((B)%CHAR_BIT))&1)
michael@0: #define S_SET(S, B)    ((S)->bits[(B)/CHAR_BIT]|=(1<<((B)%CHAR_BIT)))
michael@0: #define S_CLR(S, B)    ((S)->bits[(B)/CHAR_BIT]&=~(1<<((B)%CHAR_BIT)))
michael@0: #define S_VAL(S, B)    ((S)->base+(2*(B)))
michael@0: #define S_BIT(S, V)    (((V)-((S)->base))/2)
michael@0: 
michael@0: int main(int argc, char *argv[])
michael@0: {
michael@0:   sieve   s;
michael@0:   long    pr, *p;
michael@0:   int     c, ix, cur = 0;
michael@0: 
michael@0:   if(argc < 2) {
michael@0:     fprintf(stderr, "Usage: %s <width>\n", argv[0]);
michael@0:     return 1;
michael@0:   }
michael@0: 
michael@0:   c = atoi(argv[1]);
michael@0:   if(c < 0) c = -c;
michael@0: 
michael@0:   fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c);
michael@0: 
michael@0:   sieve_init(&s, 3, c);
michael@0: 
michael@0:   c = 0;
michael@0:   while((pr = sieve_next(&s)) > 0) {
michael@0:     ++c;
michael@0:   }
michael@0: 
michael@0:   p = calloc(c, sizeof(long));
michael@0:   if(!p) {
michael@0:     fprintf(stderr, "%s: out of memory after first half\n", argv[0]);
michael@0:     sieve_clear(&s);
michael@0:     exit(1);
michael@0:   }
michael@0: 
michael@0:   fprintf(stderr, "%s: half done ... \n", argv[0]);
michael@0: 
michael@0:   for(ix = 0; ix < s.nbits; ix++) {
michael@0:     if(S_ISSET(&s, ix)) {
michael@0:       p[cur] = S_VAL(&s, ix);
michael@0:       printf("%ld\n", p[cur]);
michael@0:       ++cur;
michael@0:     }
michael@0:   }
michael@0: 
michael@0:   sieve_reset(&s, p[cur - 1]);
michael@0:   fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur);
michael@0:   for(ix = 0; ix < cur; ix++) {
michael@0:     sieve_cross(&s, p[ix]);
michael@0:     if(!(ix % 1000))
michael@0:       fputc('.', stderr);
michael@0:   }
michael@0:   fputc('\n', stderr);
michael@0: 
michael@0:   free(p);
michael@0: 
michael@0:   fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]);
michael@0:   c = 0;
michael@0:   while((pr = sieve_next(&s)) > 0) {
michael@0:     ++c;
michael@0:   }
michael@0:   
michael@0:   fprintf(stderr, "%s: done!\n", argv[0]);
michael@0:   for(ix = 0; ix < s.nbits; ix++) {
michael@0:     if(S_ISSET(&s, ix)) {
michael@0:       printf("%ld\n", S_VAL(&s, ix));
michael@0:     }
michael@0:   }
michael@0: 
michael@0:   sieve_clear(&s);
michael@0: 
michael@0:   return 0;
michael@0: }
michael@0: 
michael@0: void sieve_init(sieve *sp, long base, int nbits)
michael@0: {
michael@0:   sp->size = (nbits / CHAR_BIT);
michael@0: 
michael@0:   if(nbits % CHAR_BIT)
michael@0:     ++sp->size;
michael@0: 
michael@0:   sp->bits = calloc(sp->size, sizeof(byte));
michael@0:   memset(sp->bits, UCHAR_MAX, sp->size);
michael@0:   if(!(base & 1))
michael@0:     ++base;
michael@0:   sp->base = base;
michael@0:   
michael@0:   sp->next = 0;
michael@0:   sp->nbits = sp->size * CHAR_BIT;
michael@0: }
michael@0: 
michael@0: void sieve_grow(sieve *sp, int nbits)
michael@0: {
michael@0:   int  ns = (nbits / CHAR_BIT);
michael@0: 
michael@0:   if(nbits % CHAR_BIT)
michael@0:     ++ns;
michael@0: 
michael@0:   if(ns > sp->size) {
michael@0:     byte   *tmp;
michael@0:     int     ix;
michael@0: 
michael@0:     tmp = calloc(ns, sizeof(byte));
michael@0:     if(tmp == NULL) {
michael@0:       fprintf(stderr, "Error: out of memory in sieve_grow\n");
michael@0:       return;
michael@0:     }
michael@0: 
michael@0:     memcpy(tmp, sp->bits, sp->size);
michael@0:     for(ix = sp->size; ix < ns; ix++) {
michael@0:       tmp[ix] = UCHAR_MAX;
michael@0:     }
michael@0: 
michael@0:     free(sp->bits);
michael@0:     sp->bits = tmp;
michael@0:     sp->size = ns;
michael@0: 
michael@0:     sp->nbits = sp->size * CHAR_BIT;
michael@0:   }
michael@0: }
michael@0: 
michael@0: long sieve_next(sieve *sp)
michael@0: {
michael@0:   long out;
michael@0:   int  ix = 0;
michael@0:   long val;
michael@0: 
michael@0:   if(sp->next > sp->nbits)
michael@0:     return -1;
michael@0: 
michael@0:   out = S_VAL(sp, sp->next);
michael@0: #ifdef DEBUG
michael@0:   fprintf(stderr, "Sieving %ld\n", out);
michael@0: #endif
michael@0: 
michael@0:   /* Sieve out all multiples of the current prime */
michael@0:   val = out;
michael@0:   while(ix < sp->nbits) {
michael@0:     val += out;
michael@0:     ix = S_BIT(sp, val);
michael@0:     if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */
michael@0:       S_CLR(sp, ix);
michael@0: #ifdef DEBUG
michael@0:       fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix);
michael@0: #endif
michael@0:     }
michael@0:   }
michael@0: 
michael@0:   /* Scan ahead to the next prime */
michael@0:   ++sp->next;
michael@0:   while(sp->next < sp->nbits && !S_ISSET(sp, sp->next))
michael@0:     ++sp->next;
michael@0: 
michael@0:   return out;
michael@0: }
michael@0: 
michael@0: void sieve_cross(sieve *sp, long val)
michael@0: {
michael@0:   int  ix = 0;
michael@0:   long cur = val;
michael@0: 
michael@0:   while(cur < sp->base)
michael@0:     cur += val;
michael@0: 
michael@0:   ix = S_BIT(sp, cur);
michael@0:   while(ix < sp->nbits) {
michael@0:     if(cur & 1) 
michael@0:       S_CLR(sp, ix);
michael@0:     cur += val;
michael@0:     ix = S_BIT(sp, cur);
michael@0:   }
michael@0: }
michael@0: 
michael@0: void sieve_reset(sieve *sp, long base)
michael@0: {
michael@0:   memset(sp->bits, UCHAR_MAX, sp->size);
michael@0:   sp->base = base;
michael@0:   sp->next = 0;
michael@0: }
michael@0: 
michael@0: void sieve_clear(sieve *sp)
michael@0: {
michael@0:   if(sp->bits) 
michael@0:     free(sp->bits);
michael@0: 
michael@0:   sp->bits = NULL;
michael@0: }