~aritra1911/prime

5d336d1936e8e9741091d9216f5d74575fcbe180 — Aritra Sarkar 2 years ago cb552ec
Use entire words for sieve of eratosthenes
3 files changed, 39 insertions(+), 11 deletions(-)

M .gitignore
M eratosthenes/Makefile
M eratosthenes/eratosthenes.c
M .gitignore => .gitignore +1 -0
@@ 2,3 2,4 @@
*.o
lucas-lehmer/lucas-lehmer
sha_consts/sha_consts
eratosthenes/pi

M eratosthenes/Makefile => eratosthenes/Makefile +1 -1
@@ 1,4 1,4 @@
CFLAGS = -Wall -Wextra -pedantic -std=c99 -g -O0
CFLAGS = -Wall -Wextra -pedantic -std=c99 -g -O3
LDFLAGS = -lm

SRCS = eratosthenes.c main.c

M eratosthenes/eratosthenes.c => eratosthenes/eratosthenes.c +37 -10
@@ 1,9 1,11 @@
#include <stdlib.h>
#include <math.h>

unsigned long int sieve_of_eratosthenes(const unsigned long int n,
unsigned long int sieve_of_eratosthenes(unsigned long int n,
                                        unsigned long int **primes)
{
    if (n <= 1) return 0;

    /* Sieve of Eratosthenes */

    /* Pseudocode :


@@ 12,8 14,16 @@ unsigned long int sieve_of_eratosthenes(const unsigned long int n,
     * Time Complexity : O(n log log n)  [according to Wikipedia]
     */

    /* TODO: We're wasting 7 bits here */
    char *a = calloc(n - 1, sizeof *a);
    /* Let A be an array of Boolean values, indexed by integers 0 to
     * (n - 2) for integers 2 through n. Hence the length of a is
     * (n - 1). */

    size_t *a;
    size_t num_bits = 8 * sizeof *a;
    size_t words = ceil((n - 2) / (double) num_bits);

    /* initially set all bits to 0 */
    a = calloc(words, sizeof *a);

    if ( primes ) {
        *primes = malloc((n - 1) * sizeof *primes);


@@ 21,20 31,37 @@ unsigned long int sieve_of_eratosthenes(const unsigned long int n,

    unsigned long int limit = sqrt(n);
    for (unsigned long int i = 2; i <= limit; i++) {
        if ( !a[i - 2] ) {
        // if ( !a[i - 2] ) {
        if ( !( a[(i - 2) / num_bits] & ((size_t) 1 << (i - 2) % num_bits) ) ) {
            for (unsigned long int j = i * i; j <= n; j += i) {
                a[j - 2] = 1;
                // a[j - 2] = 1;
                a[(j - 2) / num_bits] |= (size_t) 1 << (j - 2) % num_bits;
            }
        }
    }

    int p = 0;
    for (unsigned long int i = 0; i < (n - 1); i++) {
        if ( !a[i] ) {
            if ( primes ) {
                (*primes)[p] = i + 2;

    if ( !primes ) {
        /* Calculate the number of set bits */
        for (size_t i = 0; i < words; i++) {
            while ( a[i] ) {
                p++;

                /* Consume least significant power of 2 in a[i] */
                a[i] &= a[i] - 1;
            }
            p++;
        }

        free(a);

        /* # of unset bits = # of bits in total - # of set bits */
        return (n - 1) - p;
    }

    for (unsigned long int i = 0; i < (n - 1); i++) {
        if ( !( a[i / num_bits] & ((size_t) 1 << i % num_bits) ) ) {
	    (*primes)[p++] = i + 2;
        }
    }