Primes by sieve with wheel factorization in C
/******************************************
/
/ wheelk2.c
/
/ Find the primes up to n, using a
/ Sieve of Eratosthenes, optimized by
/ wheel factorization
/
/ For "Exploring Small Primes with C"
/
/ Paul Soper
/
/ May, 2014
/
/******************************************/
/******************************************
/
/ This program uses a k=2 wheel, with M = 6
/
/******************************************/
#include "stdio.h"
#include "math.h" /* because we use sqrt() */
#include "stdlib.h" /* because we use malloc */
int main() {
/* we will use a k = 2 wheel, which will have
a stepper wheel of six elements */
int n = 200000; /* the size of the sieve */
/* allocate memory for the sieve */
int* si = (int*)malloc(sizeof(int)*(n+1));
/* now we need M for our wheel, using the first -
2 primes because this is a k=2 wheel */
int M = 2 * 3;
int i = 0; /* a general counter */
/* allocate memory for the stepper wheel
we need M bytes */
int* W = (int*)malloc(sizeof(int)*(M+1));
/* Define the stepper wheel using this algorithm:
If x is divisible by 2 or 3, W[x] = 0.
Otherwise, W[x] = the distance to the next integer
after x which is not divisible by 2 or 3.
We are hard-coding it here, to make it
easier to see, but the algorithm could
easily be written */
W[1] = 4; W[2] = 0; W[3] = 0;
W[4] = 0; W[5] = 2; W[6] = 0;
/* Now hard-code the first M elements of the
sieve, and si[0], with this algorithm:
If x is prime, si[x] = 1,
and otherwise, si[x] = 0 */
*si = 0;
*(si+1) = 0;
*(si+2) = 1;
*(si+3) = 1;
*(si+4) = 0;
*(si+5) = 1;
*(si+6) = 0;
/* Now define the rest of the sieve, with this
algorithm: If W[x mod M] = 0, si[x] = 0.
Otherwise, si[x] = 1. This is the same as
saying if x is not co-prime with 2 or 3,
we don't have to consider it as a candidate
prime, so start the sieve with only co-primes
with 2 or 3 marked as candidates */
for (i = 7; i <= n; ++i) {
if (*(W+(i%M)) == 0) {
*(si+i) = 0;
}
else {
*(si+i) = 1;
}
}
/* Now step through the sieve using the stepper wheel */
int p = 5;
int f;
for (p = 5; p <= (int)sqrt(n); ++p) {
if (*(si+p) == 1) {
for (f = p; f <= (n/p); f = f + *(W + (f%M))) {
*(si + p*f) = 0;
}
}
}
/* allocate memory for an array into which to
collect all of the primes - there are
certainly not more than 2n/log(n) of them */
int* pr = (int*)malloc(sizeof(int)*(2*n/log(n)));
int c = 0; /* a count of the primes */
/* step through the sieve, and if a prime has been
identified, store it in pr[], and count it */
for (i = 0; i < n; ++i) {
if (*(si + i) != 0) {
*(pr + c) = i;
++c;
}
}
/* print out the list of primes */
for (i = 0; i < c; ++i) {
printf("Prime %d: %d\n", i + 1, *(pr + i));
}
/* free the allocated memory */
free (pr);
free (W);
free (si);
return (0);
}
Like this:
Like Loading...