findprimes.c

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

#define PRIMESIZE sizeof(uint64_t)

static uint64_t * primes_found = NULL;
static uint64_t   prime_count  = 0;

/*
 *   findPrimes(maxPrime, maxCount)
 *
 *  Either finds all primes up to maxPrime, or finds as many primes as
 *  maxCount, whichever is reached first. 0 in one means "not a factor".
 *  0 in both is disallowed, and returns 0.
 *
 *  Upon normal completion, it will return the amount of primes that
 *  have been found in all executions.
 */

uint64_t findPrimes(uint64_t maxPrime, uint64_t maxCount)
{
    if (maxPrime == 0 && maxCount == 0) { return 0; }

    char primeCond = maxPrime != 0;
    char countCond = maxCount != 0;

    if (primes_found == NULL)
    {
        primes_found = malloc(PRIMESIZE);
        primes_found[0] = 2;
        prime_count     = 1;
    }

    uint64_t prime = primes_found[prime_count-1] + 1;

    long double maxTestNum_f;
    uint32_t    maxTestNum_i;

    uint64_t    curIndex;
    uint64_t    curPrime;
    uint8_t     isPrime;

    long double divResult_f;
    uint64_t    divResult_i;

    while (1)
    {
        // Exit conditions
        if (primeCond && maxPrime < prime)        { break; }
        if (countCond && maxCount <= prime_count) { break; }

        // No number has a prime factor larger than the square root of it
        maxTestNum_f = sqrtl((long double)prime);
        maxTestNum_i = (uint32_t)maxTestNum_f;

        isPrime = 1;

        // Test against all previous primes

        for (curIndex = 0; curIndex < prime_count; curIndex++)
        {
            curPrime = primes_found[curIndex];
            if (curPrime > maxTestNum_i) { break; }

            divResult_f = (long double)prime / (long double)curPrime;
            divResult_i = prime / curPrime;

            if (divResult_f == (long double)divResult_i)
            {
                isPrime = 0;
                break;
            }
        }

        if (isPrime)
        {
            // Extend by one prime
            primes_found = realloc(primes_found, PRIMESIZE * (++prime_count));
            primes_found[prime_count - 1] = prime;
        }

        prime++;
    }

    return prime_count;
}

/*   nthPrime(index)
 *
 *  Finds the nth prime, assuming it's been calculated. If it hasn't
 *  been calculated, return 0.
 *
 *  For reference: nthPrime(1) = 2, nthPrime(2) = 3, nthPrime(3) = 5
 */

uint64_t nthPrime(uint64_t index)
{
    if (--index >= prime_count) { return 0; }
    return primes_found[index];
}

/*   primeIndex(checknum)
 *
 *  Returns the index of checknum in primes_found plus one.
 *  A value of 0 means checknum is 0, 1, or composite.
 */

uint64_t primeIndex(uint64_t checknum)
{
    // Already known cases - don't waste time on them
    if (checknum <  2) { return 0; }
    if (checknum == 2) { return 1; }
    if (checknum == primes_found[prime_count-1]) { return prime_count; }

    if (checknum > primes_found[prime_count-1])
    {
        findPrimes(checknum, 0);
        if (checknum != primes_found[prime_count-1]) { return 0; }
    }

    uint64_t bottomSlice = 0;
    uint64_t topSlice    = prime_count - 1;
    uint64_t sliceRange;
    uint64_t middle;
    uint64_t middleItem;
    uint64_t ret = 0;

    // Binary search for the index

    while (1)
    {
        // Once our range goes to 1, the only value left is either our check
        // number or our check num is composite

        if (topSlice == bottomSlice)
        {
            if (primes_found[topSlice] == checknum) { ret = topSlice; break; }
            else { return 0; }
        }

        if (primes_found[topSlice] == checknum)    { ret = topSlice; break; }
        if (primes_found[bottomSlice] == checknum) { ret = bottomSlice; break; }

        sliceRange = topSlice - bottomSlice;
        middle = (topSlice + bottomSlice) / 2;
        middleItem = primes_found[middle];

        if (middleItem == checknum) { ret = middle; break; }

        if (middleItem > checknum) { topSlice -= (sliceRange / 2); }
        else { bottomSlice += (sliceRange / 2); }
    }

    return ret + 1;
}


/*   primeFactors(tofactor, **retfactors)
 *
 *  Returns a value equal to the amount of prime factors in the number
 *  tofactor. The prime factors are stored in the uint64_t* address
 *  given to us. The values in retfactors, when multiplied together,
 *  will get you tofactor.
 *
 *  Example:
 *
 *  uint64_t * factors;
 *  uint32_t factorcount = primeFactors(29374, &factors);
 *
 *  gets you:
 *
 *  factors     = {2, 19, 773}
 *  factorcount = 3
 */

uint32_t primeFactors(uint64_t tofactor, uint64_t ** retfactors)
{
    uint64_t * factors = NULL;
    uint32_t factorcount = 0;
    uint32_t i, curPrime;

    long double maxTestNum_f;
    uint32_t    maxTestNum_i;

    long double divResult_f;
    uint64_t    divResult_i;

    uint8_t     noMoreFactors = 0;

    while (1)
    {
        // No number has a prime factor larger than the square root of it
        maxTestNum_f = sqrtl((long double)tofactor);
        maxTestNum_i = (uint32_t)maxTestNum_f;

        i = 1;

        while (1)
        {
            curPrime = nthPrime(i);
            if (curPrime == 0) { findPrimes(0, i); }

            // Once we're looking above the square root, we've exhausted all factors
            if (curPrime > maxTestNum_i) { noMoreFactors = 1; break; }

            divResult_f = (long double)tofactor / (long double)curPrime;
            divResult_i = tofactor / curPrime;

            if (divResult_f == (long double)divResult_i)
            {
                tofactor = divResult_i;
                break;
            }

            i++;
        }

        if (noMoreFactors)
        {
            factors = realloc(factors, sizeof(uint64_t) * (++factorcount));
            factors[factorcount-1] = tofactor;
            break;
        }
        else
        {
            factors = realloc(factors, sizeof(uint64_t) * (++factorcount));
            factors[factorcount-1] = curPrime;
        }
    }

    *retfactors = factors;
    return factorcount;
};

int main(void)
{
    uint64_t gotoPrime = 10001;
    uint64_t testPrime = 1298909;  // known prime
    uint64_t testNum   = 29374;  // known composite
    uint64_t index;
    uint64_t * factors;
    uint32_t factorcount;

    printf("Found %lu primes\n", findPrimes(0, gotoPrime));
    printf("Prime %lu is %lu\n", gotoPrime, nthPrime(gotoPrime));

    index = primeIndex(testPrime);
    printf("Prime index of %lu is %lu\n", testPrime, index);
    printf("Reverse lookup - %lu = nthPrime(%lu) = %lu\n", testPrime, index, nthPrime(index));

    factorcount = primeFactors(testNum, &factors);

    printf("%lu has %u prime factors:\n", testNum, factorcount);

    uint32_t i;

    for (i = 0; i < factorcount; i++)
    {
        printf(" * %lu\n", factors[i]);
    }

    return 0;
}