roman10
Previous post covers two naive methods for primality testing. This post introduces an algorithm to generate prime numbers called Sieve of Eratosthenes. With this algorithm, an optimization can be applied to speed up the primality testing.
Sieve of Eratosthenes
Suppose we want to generate prime numbers up to N. The algorithm consists of two loops, and it can be described as below,
- Create a mask for all integers up to N.
- In the first loop, let x be 2, mark 2x, 3x, 4x for all n such that nx < N.
- finding an integer bigger than x and not marked, let it be the new x, and repeat the step above
- the first loop terminates when we cannot find an integer bigger than x and not marked
- In the second loop, starting from 2, collects the numbers that are not marked, which are the prime numbers up to N.
Below is an implementation of the Sieve of Eratosthenes algorithm,
/**
generate prime numbers
**/
#include <stdio.h>
#include <stdlib.h>
#define MAXV 10000000
#define ININ 1000
unsigned char bits[MAXV/8];
void setbit(int pos) {
bits[pos/8] |= (0x01 << (pos%8));
}
int testbit(int pos) {
return (bits[pos/8] & (0x01 << (pos%8)));
}
void genprime(int n, int **pr, int *npr) {
int i, j, k;
for (i = 2; i < n; ++i) {
if (!testbit(i)) {
for (j = 2; i*j < n; ++j) {
setbit(i*j);
}
}
}
k = ININ;
*pr = malloc(sizeof(int) * k);
*npr = 0;
for (i = 2, j = 0; i < n; ++i) {
if (!testbit(i)) {
(*pr)[j++] = i;
(*npr)++;
if (j == k) {
k += ININ;
*pr = realloc(*pr, sizeof(int)*k);
}
}
}
}
int main(int argc, char **argv) {
int n = atoi(argv[1]);
int *pr;
int npr;
int i;
struct timeval stTime, edTime;
if (n > MAXV) {
printf("the max n is %d\n", MAXV);
return 0;
}
gettimeofday(&stTime, NULL);
genprime(n, &pr, &npr);
printf("No. of prime numbers: %d\n", npr);
for (i = 0; i < npr; ++i) {
printf("%d ", pr[i]);
}
printf("\n");
gettimeofday(&edTime, NULL);
printf("time: %u:%u\n", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));
if (pr != NULL) {
free(pr);
}
return 0;
}
Note that the code above use a single bit to indicate if an integer is marked or not.
Save the code to prim2.c, and then compile the code with the command below,
gcc -o prim2 prim2.c
Below is a screenshot of some simple tests,
Figure 1. Simple Tests of Implementation for Sieve of Eratosthenes
Primality Testing with Sieve of Eratosthenes for Precomputation
One way to speed up primality testing is to precompute a list of prime numbers up to N. Then we first check if an input integer n is divisible by the prime numbers generated. If the number is not divisible by any of those prime numbers, we then proceed to the test mentioned in the naive approach post.
Below is an implementation of this precomputation approach,
/**
a primality test with precomputation
**/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAXV 10000000
#define ININ 1000
unsigned char bits[MAXV/8];
void setbit(int pos) {
bits[pos/8] |= (0x01 << (pos%8));
}
int testbit(int pos) {
return (bits[pos/8] & (0x01 << (pos%8)));
}
void genprime(int n, int **pr, int *npr) {
int i, j, k;
for (i = 2; i < n; ++i) {
if (!testbit(i)) {
for (j = 2; i*j < n; ++j) {
setbit(i*j);
}
}
}
k = ININ;
*pr = malloc(sizeof(int) * k);
*npr = 0;
for (i = 2, j = 0; i < n; ++i) {
if (!testbit(i)) {
(*pr)[j++] = i;
(*npr)++;
if (j == k) {
k += ININ;
*pr = realloc(*pr, sizeof(int)*k);
}
}
}
}
int prime(int n, int *pr, int npr, int x) {
int i;
int j;
for (i = 0; i < npr; ++i) {
if (n == pr[i]) {
return 1;
}
if (n%pr[i] == 0) {
printf("%d is divisble by %d\n", n, pr[i]);
return 0;
}
}
if (n == 2) {
return 1;
}
if (!(n&1)) {
printf("%d is divisible by 2\n", n);
return 0;
}
j = sqrt(n);
//to be safe, run to j+1
i = (x+1) > 3 ? (x+1):3;
if (!(i&1)) {
++i; //make sure i starts as odd number
}
for (; i <= j + 1; i+=2) {
if (n%i == 0) {
printf("%d is divisible by %d\n", n, i);
return 0;
}
}
return 1;
}
int main(int argc, char **argv) {
int m, n;
int *pr;
int npr;
int i;
struct timeval stTime, stTime2, edTime;
if (argc < 3) {
printf("./prim3 <primality test number> <prime number generation limit>\n");
return 0;
}
m = atoi(argv[1]);
n = atoi(argv[2]);
if (n > MAXV) {
printf("the max n is %d\n", MAXV);
return 0;
}
gettimeofday(&stTime, NULL);
genprime(n, &pr, &npr);
gettimeofday(&stTime2, NULL);
//printf("No. of prime numbers: %d\n", npr);
/* for (i = 0; i < npr; ++i) {
printf("%d\n", pr[i]);
}*/
if (prime(m, pr, npr, n)) {
printf("%d is prime number\n", m);
} else {
printf("%d is not prime number\n", m);
}
gettimeofday(&edTime, NULL);
printf("total time: %u:%u\n", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));
printf("test prime time: %u:%u\n", (unsigned int)(edTime.tv_sec - stTime2.tv_sec), (unsigned int)(edTime.tv_usec - stTime2.tv_usec));
if (pr != NULL) {
free(pr);
}
return 0;
}
Save the source code to prim3.c, and then compile the code with the command below,
gcc -o prim3 prim3.c –lm
Below are some simple tests for the compiled program,
roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2 10
2 is prime number
total time: 0:144
test prime time: 0:48
roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 3 10
3 is prime number
total time: 0:60
test prime time: 0:22
roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483647 1000
2147483647 is prime number
total time: 0:469
test prime time: 0:302
roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483643 1000
2147483643 is divisble by 3
2147483643 is not prime number
total time: 0:223
test prime time: 0:57
roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483641 1000
2147483641 is divisible by 2699
2147483641 is not prime number
total time: 0:250
test prime time: 0:80
References:
1. Wikipedia: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
If you enjoyed this article, please consider sharing it!
