prime sums

output

bclement:primes bclement$ gcc -O666 primes.c
bclement:primes bclement$ time ./a.out
There is no sum for 0!
There is no sum for 1!
There is no sum for 2!
There is no sum for 3!
There is no sum for 4!
There is no sum for 5!
There is no sum for 6!
There is no sum for 7!

real	3m56.149s
user	27m43.326s
sys	0m11.647s
bclement:primes bclement$ time ./a.out
There is no sum for 0!
There is no sum for 1!
There is no sum for 2!
There is no sum for 3!
There is no sum for 4!
There is no sum for 5!
There is no sum for 6!
There is no sum for 7!

real	0m0.593s
user	0m0.588s
sys	0m0.003s

primality.c

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

//binary expontentiation
unsigned long modPow(unsigned long base, unsigned long pow, unsigned long mod) {
    unsigned long res = 1;

    while (pow) {
        if (pow % 2 == 1) {
            res *= base;
            res %= mod;
        }
        base *= base;
        base %= mod;
        pow >>= 1;
    }

    return res;
}

//due to some statements made about finite fields we know that a^(c - 1)
//= 1 and there are no non-trivial roots of unity mod p (i.e. an element
//f such that f^2 = 1 mod p).
//for more in depth discussion read here;
//https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Concepts
unsigned long check(unsigned long base, unsigned long maxPower, unsigned long basePower, unsigned long mod) {
    unsigned long shouldBeTrivial = modPow(base, basePower, mod);
    // test to see the base case is 1
    if (shouldBeTrivial == 1) {
        return 1;
    }
    // the power of 2 in the expontent
    unsigned long exponentPower = 0;
    // test all squares for non-triviality
    while (exponentPower < maxPower) {
        if (shouldBeTrivial == mod - 1) {
            return 1;
        }
        shouldBeTrivial = (shouldBeTrivial * shouldBeTrivial) % mod;
        exponentPower += 1;
    }
    return 0;
}

unsigned long isPrime(unsigned long n) {
    if(n == 2){
        return 1;
    }

    // factor n - 1, which should be even, into the form (2^s)*d where d is
    // an odd factor
    unsigned long d = n - 1;
    unsigned long s = 0;

    while(d % 2 == 0){
        s += 1;
        d /= 2;
    }

    // the list of factors you have to bound by are listed out here:
    // https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants_of_the_test
    return check(2, s, d, n) && check(7, s, d, n) && check(61, s, d, n);

}

unsigned long * getPrimePairFor(unsigned long n, unsigned long * primes, unsigned long numPrimes) {
    if (n < 8) {
        return 0;
    }

    unsigned long * primeList = malloc(4*sizeof(unsigned long));
    if (n % 2 == 0) {
        primeList[0] = 2;
        primeList[1] = 2;
        n -= 4;
    } else {
        primeList[0] = 2;
        primeList[1] = 3;
        n -= 5;
    }

    unsigned long i;
    unsigned long prime;
    for (i = 0; i < numPrimes; i++) {
        prime = primes[i];
        if (isPrime(n - prime)) {
            break;
        }
    }

    primeList[2] = prime;
    primeList[3] = n - prime;

    return primeList;

}

unsigned long * getNPrimes(unsigned long n) {
    unsigned long *primes = malloc(sizeof(unsigned long) * n);
    memset(primes, 0, sizeof(unsigned long) * n);
    primes[0] = 2;
    int numPrimes = 1;
    int candidate = 3;
    int isPrime = 1;
    while(numPrimes < n) {
        for (int i = 0; i < numPrimes; i++) {
            if (candidate % primes[i] == 0) {
                isPrime = 0;
                break;
            }
        }
        if (isPrime) {
            primes[numPrimes++] = candidate;
        }
        isPrime = 1;
        candidate += 2;
    }
    return primes;
}

int main() {
    unsigned long n = 1000000000;

    unsigned long numPrimes = ceil(sqrt(sqrt(n)));
    unsigned long * primes = getNPrimes(numPrimes);

    for (unsigned long i = 0; i <= n; i++) {
        unsigned long * primeList = getPrimePairFor(i, primes, numPrimes);
        if (!primeList) {
            printf("There is no sum for %lu!\n", i);
        } else {
            //printf("%lu + %lu + %lu + %lu == %lu\n", primeList[0],primeList[1],primeList[2],primeList[3],i);
        }
    }

    free(primes);
}

primes.c

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <pthread.h>
#include <math.h>

#define NUM_THREADS  8

typedef struct primePair primePair;

struct primePair {
   unsigned first;
   unsigned second;
};

typedef struct primeWorkerArgs primeWorkerArgs;

struct primeWorkerArgs {
    unsigned * primes;
    unsigned numPrimes;
    unsigned n;
    unsigned offset;
    primePair * primeSums;
};

unsigned * genPrimesUnder(unsigned n, unsigned numPrimes) {
    unsigned *candidates = malloc(sizeof(unsigned) * n + 1);
    memset(candidates, 1, sizeof(unsigned) * n + 1);
    unsigned limit = ceil(sqrt(n));
    for (unsigned i = 2; i < limit; i++) {
        if (candidates[i]) {
            for (unsigned j = 2*i; j < n; j += i) {
                candidates[j] = 0;
            }
        }
    }
    unsigned *primes = malloc(sizeof(unsigned) * numPrimes);
    memset(primes, 0, sizeof(unsigned) * numPrimes);
    unsigned currentPrime = 0;
    for (unsigned i = 2; i < n; i++) {
        if (candidates[i]) {
            primes[currentPrime++] = i;
        }
    }
    free(candidates);
    return primes;
}

void primeSumWorker(void *input) {
    primeWorkerArgs *args = (primeWorkerArgs *)(input);
    for (unsigned i = args->offset; i < args->numPrimes; i += NUM_THREADS) {
        if (args->primes[i] == 0) {
            break;
        }
        for (unsigned j = i; j < args->numPrimes; j++) {
            if (args->primes[j] == 0) {
                break;
            }
            unsigned sum = args->primes[i] + args->primes[j];
            if (sum > args->n) {
                break;
            }
            primePair primeSum = {.first = args->primes[i], .second = args->primes[j]};
            args->primeSums[sum] = primeSum;
        }
    }
}

primePair * genPrimeSums(unsigned * primes, unsigned numPrimes, unsigned n) {
    primePair * primeSums = malloc(sizeof(primePair) * (n + 1));
    memset(primeSums, 0, sizeof(primePair) * (n + 1));
    pthread_t threads[NUM_THREADS];
    int rc;
    primeWorkerArgs workerArgList[NUM_THREADS];
    for(unsigned t=0;t<NUM_THREADS;t++){
        primeWorkerArgs args = {primes, numPrimes, n, t, primeSums};
        workerArgList[t] = args;
        rc = pthread_create(&threads[t], NULL, (void *) &primeSumWorker, (void *) &workerArgList[t]);
        if (rc){
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
    }
    for(unsigned t=0;t<NUM_THREADS;t++){
        pthread_join(threads[t], NULL);
    }

    return primeSums;
}

unsigned * getPrimePairFor(unsigned n, primePair * primeSums) {
    if (n < 8) {
        return 0;
    }

    unsigned * primeList = malloc(4*sizeof(unsigned));
    if (n % 2 == 0) {
        primeList[0] = 2;
        primeList[1] = 2;
        n -= 4;
    } else {
        primeList[0] = 2;
        primeList[1] = 3;
        n -= 5;
    }

    primePair pair = primeSums[n];
    if (pair.first == 0) {
        return 0;
    }

    primeList[2] = pair.first;
    primeList[3] = pair.second;

    return primeList;

}

int main() {
    unsigned n = 10000000;
    unsigned numPrimes = ceil(n/log(n)*(1 + 1.2762f/log(n)));
    unsigned * primes = genPrimesUnder(n, numPrimes);
    primePair * primeSums = genPrimeSums(primes, numPrimes, n);
    for (unsigned i = 0; i <= n; i++) {
        unsigned * primeList = getPrimePairFor(i, primeSums);
        if (!primeList) {
            printf("There is no sum for %u!\n", i);
        }
    }
    free(primes);
    free(primeSums);
}