Project Euler #27 (Quadratic Primes)

main.c

#include <assert.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

ptrdiff_t indexOfPrime(uint64_t num);

const uint64_t *findPrimes(size_t amount);

inline bool isPrime(uint64_t num) {
    return indexOfPrime(num) >= 0;
}

inline unsigned int primeAt(size_t index) {
    return findPrimes(index + 1)[index];
}

uint64_t quadraticPrimeMaxBound(int64_t a, int64_t b) {
    assert(isPrime(abs(a)) || a == 1);
    assert(b > 0 && isPrime(b));
    int64_t n = 0;
    int64_t result;
    do {
        result = (n * n) + (a * n) + b;
        n++;
    } while (result > 0 && isPrime((uint64_t) result));
    return (uint64_t) n - 1;
}

int main(void) {
    // Disable stdout buffering
    setvbuf(stdout, NULL, _IONBF, 0);

    int64_t highest_a = 0, highest_b = 0;
    uint64_t highest_bound = 0;
    int64_t a, b;
    for (a = -999; a < 1000; a++) {
        if (a != 1 && !isPrime(abs(a))) continue;
        for (b = 0; b <= 1000; b++) {
            if (!isPrime(abs(b))) continue;
            uint64_t bound = quadraticPrimeMaxBound(a, b);
            if (bound > highest_bound) {
                highest_bound = bound;
                highest_a = a;
                highest_b = b;
            }
        }
    }
    printf("n^2 + %lldn + %lld produces up to %llu primes\n", highest_a, highest_b, highest_bound);

    return 0;
}

//
// Prime numbers
//

static size_t numKnownPrimes = 0;
static uint64_t *knownPrimes = NULL;
static uint64_t highestKnownPrime = 0;

void computeMorePrimes(size_t amount) {
    assert(amount > 0);
    // Increase buffer size
    size_t desiredSize = numKnownPrimes + amount;
    if (desiredSize < 2) desiredSize = 2;
    knownPrimes = realloc(knownPrimes, sizeof(uint64_t) * desiredSize);
    if (knownPrimes == NULL) {
        fputs("Out of memory\n", stderr);
        exit(1);
    }
    if (numKnownPrimes < 2) {
        /*
         * Initialize our state with the primes '2' and '3'.
         * We can't just put it two, since the main-loop assumes that
         * 'highestKnownPrime' isn't an odd number.
         */
        knownPrimes[0] = 2;
        knownPrimes[1] = 3;
        highestKnownPrime = 3;
        numKnownPrimes = 2;
    }
    size_t index;
    for (index = numKnownPrimes; index < desiredSize; index++) {
        inline bool isNewPrime(uint64_t num) {
            // Start at the second prime number, since we assume 'num' is odd
            size_t previousIndex;
            for (previousIndex = 1; previousIndex < index; previousIndex++) {
                int previousPrime = knownPrimes[previousIndex];
                if (num % previousPrime == 0) {
                    return false;
                }
            }
            return true; // It's not divisible by any previous primes
        }
        uint64_t potentialPrime = highestKnownPrime;
        do {
            potentialPrime += 2;
        } while (!isNewPrime(potentialPrime));
        knownPrimes[index] = potentialPrime;
        highestKnownPrime = potentialPrime;
    }
    numKnownPrimes = desiredSize;
}

ptrdiff_t indexOfPrime(uint64_t num) {
    if (num < 2) return -1; // Safety check to prevent a NPE
    while (num > highestKnownPrime) {
        computeMorePrimes(25);
    }
    // Perform a binary search over the previous prime numbers
    // Since the prime numbers are sorted, it'll work fine to find the index
    size_t low = 0;
    size_t high = numKnownPrimes;
    while (low <= high) {
        size_t mid = (low + high) >> 1;
        uint64_t midVal = knownPrimes[mid];
        if (midVal < num) {
            low = mid + 1;
        } else if (midVal > num) {
            high = mid - 1;
        } else {
            return mid; // Found it
        }
    }
    return -1; // Not found
}

inline const uint64_t *findPrimes(size_t amount) {
    ptrdiff_t numUnknown = amount - numKnownPrimes;
    if (numUnknown > 0) {
        computeMorePrimes((size_t) numUnknown);
    }
    assert(numKnownPrimes >= amount);
    return knownPrimes;
}