The Tor Browser / file revision

 /*

  * 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.

  *

  * This Source Code Form is subject to the terms of the Mozilla Public

  * License, v. 2.0. If a copy of the MPL was not distributed with this

  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

 #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;

 }