Back to index

lightning-sunbird  0.9+nobinonly
sieve.c
Go to the documentation of this file.
00001 /*
00002  * sieve.c
00003  *
00004  * Finds prime numbers using the Sieve of Eratosthenes
00005  *
00006  * This implementation uses a bitmap to represent all odd integers in a
00007  * given range.  We iterate over this bitmap, crossing off the
00008  * multiples of each prime we find.  At the end, all the remaining set
00009  * bits correspond to prime integers.
00010  *
00011  * Here, we make two passes -- once we have generated a sieve-ful of
00012  * primes, we copy them out, reset the sieve using the highest
00013  * generated prime from the first pass as a base.  Then we cross out
00014  * all the multiples of all the primes we found the first time through,
00015  * and re-sieve.  In this way, we get double use of the memory we
00016  * allocated for the sieve the first time though.  Since we also
00017  * implicitly ignore multiples of 2, this amounts to 4 times the
00018  * values.
00019  *
00020  * This could (and probably will) be generalized to re-use the sieve a
00021  * few more times.
00022  *
00023  * ***** BEGIN LICENSE BLOCK *****
00024  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00025  *
00026  * The contents of this file are subject to the Mozilla Public License Version
00027  * 1.1 (the "License"); you may not use this file except in compliance with
00028  * the License. You may obtain a copy of the License at
00029  * http://www.mozilla.org/MPL/
00030  *
00031  * Software distributed under the License is distributed on an "AS IS" basis,
00032  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00033  * for the specific language governing rights and limitations under the
00034  * License.
00035  *
00036  * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library.
00037  *
00038  * The Initial Developer of the Original Code is
00039  * Michael J. Fromberger.
00040  * Portions created by the Initial Developer are Copyright (C) 1998
00041  * the Initial Developer. All Rights Reserved.
00042  *
00043  * Contributor(s):
00044  *
00045  * Alternatively, the contents of this file may be used under the terms of
00046  * either the GNU General Public License Version 2 or later (the "GPL"), or
00047  * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00048  * in which case the provisions of the GPL or the LGPL are applicable instead
00049  * of those above. If you wish to allow use of your version of this file only
00050  * under the terms of either the GPL or the LGPL, and not to allow others to
00051  * use your version of this file under the terms of the MPL, indicate your
00052  * decision by deleting the provisions above and replace them with the notice
00053  * and other provisions required by the GPL or the LGPL. If you do not delete
00054  * the provisions above, a recipient may use your version of this file under
00055  * the terms of any one of the MPL, the GPL or the LGPL.
00056  *
00057  * ***** END LICENSE BLOCK ***** */
00058 /* $Id: sieve.c,v 1.3 2004/04/27 23:04:37 gerv%gerv.net Exp $ */
00059 
00060 #include <stdio.h>
00061 #include <stdlib.h>
00062 #include <limits.h>
00063 
00064 typedef unsigned char  byte;
00065 
00066 typedef struct {
00067   int   size;
00068   byte *bits;
00069   long  base;
00070   int   next;
00071   int   nbits;
00072 } sieve;
00073 
00074 void sieve_init(sieve *sp, long base, int nbits);
00075 void sieve_grow(sieve *sp, int nbits);
00076 long sieve_next(sieve *sp);
00077 void sieve_reset(sieve *sp, long base);
00078 void sieve_cross(sieve *sp, long val);
00079 void sieve_clear(sieve *sp);
00080 
00081 #define S_ISSET(S, B)  (((S)->bits[(B)/CHAR_BIT]>>((B)%CHAR_BIT))&1)
00082 #define S_SET(S, B)    ((S)->bits[(B)/CHAR_BIT]|=(1<<((B)%CHAR_BIT)))
00083 #define S_CLR(S, B)    ((S)->bits[(B)/CHAR_BIT]&=~(1<<((B)%CHAR_BIT)))
00084 #define S_VAL(S, B)    ((S)->base+(2*(B)))
00085 #define S_BIT(S, V)    (((V)-((S)->base))/2)
00086 
00087 int main(int argc, char *argv[])
00088 {
00089   sieve   s;
00090   long    pr, *p;
00091   int     c, ix, cur = 0;
00092 
00093   if(argc < 2) {
00094     fprintf(stderr, "Usage: %s <width>\n", argv[0]);
00095     return 1;
00096   }
00097 
00098   c = atoi(argv[1]);
00099   if(c < 0) c = -c;
00100 
00101   fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c);
00102 
00103   sieve_init(&s, 3, c);
00104 
00105   c = 0;
00106   while((pr = sieve_next(&s)) > 0) {
00107     ++c;
00108   }
00109 
00110   p = calloc(c, sizeof(long));
00111   if(!p) {
00112     fprintf(stderr, "%s: out of memory after first half\n", argv[0]);
00113     sieve_clear(&s);
00114     exit(1);
00115   }
00116 
00117   fprintf(stderr, "%s: half done ... \n", argv[0]);
00118 
00119   for(ix = 0; ix < s.nbits; ix++) {
00120     if(S_ISSET(&s, ix)) {
00121       p[cur] = S_VAL(&s, ix);
00122       printf("%ld\n", p[cur]);
00123       ++cur;
00124     }
00125   }
00126 
00127   sieve_reset(&s, p[cur - 1]);
00128   fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur);
00129   for(ix = 0; ix < cur; ix++) {
00130     sieve_cross(&s, p[ix]);
00131     if(!(ix % 1000))
00132       fputc('.', stderr);
00133   }
00134   fputc('\n', stderr);
00135 
00136   free(p);
00137 
00138   fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]);
00139   c = 0;
00140   while((pr = sieve_next(&s)) > 0) {
00141     ++c;
00142   }
00143   
00144   fprintf(stderr, "%s: done!\n", argv[0]);
00145   for(ix = 0; ix < s.nbits; ix++) {
00146     if(S_ISSET(&s, ix)) {
00147       printf("%ld\n", S_VAL(&s, ix));
00148     }
00149   }
00150 
00151   sieve_clear(&s);
00152 
00153   return 0;
00154 }
00155 
00156 void sieve_init(sieve *sp, long base, int nbits)
00157 {
00158   sp->size = (nbits / CHAR_BIT);
00159 
00160   if(nbits % CHAR_BIT)
00161     ++sp->size;
00162 
00163   sp->bits = calloc(sp->size, sizeof(byte));
00164   memset(sp->bits, UCHAR_MAX, sp->size);
00165   if(!(base & 1))
00166     ++base;
00167   sp->base = base;
00168   
00169   sp->next = 0;
00170   sp->nbits = sp->size * CHAR_BIT;
00171 }
00172 
00173 void sieve_grow(sieve *sp, int nbits)
00174 {
00175   int  ns = (nbits / CHAR_BIT);
00176 
00177   if(nbits % CHAR_BIT)
00178     ++ns;
00179 
00180   if(ns > sp->size) {
00181     byte   *tmp;
00182     int     ix;
00183 
00184     tmp = calloc(ns, sizeof(byte));
00185     if(tmp == NULL) {
00186       fprintf(stderr, "Error: out of memory in sieve_grow\n");
00187       return;
00188     }
00189 
00190     memcpy(tmp, sp->bits, sp->size);
00191     for(ix = sp->size; ix < ns; ix++) {
00192       tmp[ix] = UCHAR_MAX;
00193     }
00194 
00195     free(sp->bits);
00196     sp->bits = tmp;
00197     sp->size = ns;
00198 
00199     sp->nbits = sp->size * CHAR_BIT;
00200   }
00201 }
00202 
00203 long sieve_next(sieve *sp)
00204 {
00205   long out;
00206   int  ix = 0;
00207   long val;
00208 
00209   if(sp->next > sp->nbits)
00210     return -1;
00211 
00212   out = S_VAL(sp, sp->next);
00213 #ifdef DEBUG
00214   fprintf(stderr, "Sieving %ld\n", out);
00215 #endif
00216 
00217   /* Sieve out all multiples of the current prime */
00218   val = out;
00219   while(ix < sp->nbits) {
00220     val += out;
00221     ix = S_BIT(sp, val);
00222     if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */
00223       S_CLR(sp, ix);
00224 #ifdef DEBUG
00225       fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix);
00226 #endif
00227     }
00228   }
00229 
00230   /* Scan ahead to the next prime */
00231   ++sp->next;
00232   while(sp->next < sp->nbits && !S_ISSET(sp, sp->next))
00233     ++sp->next;
00234 
00235   return out;
00236 }
00237 
00238 void sieve_cross(sieve *sp, long val)
00239 {
00240   int  ix = 0;
00241   long cur = val;
00242 
00243   while(cur < sp->base)
00244     cur += val;
00245 
00246   ix = S_BIT(sp, cur);
00247   while(ix < sp->nbits) {
00248     if(cur & 1) 
00249       S_CLR(sp, ix);
00250     cur += val;
00251     ix = S_BIT(sp, cur);
00252   }
00253 }
00254 
00255 void sieve_reset(sieve *sp, long base)
00256 {
00257   memset(sp->bits, UCHAR_MAX, sp->size);
00258   sp->base = base;
00259   sp->next = 0;
00260 }
00261 
00262 void sieve_clear(sieve *sp)
00263 {
00264   if(sp->bits) 
00265     free(sp->bits);
00266 
00267   sp->bits = NULL;
00268 }