Factorial using "divide, swing and conquer" algorithm - https://oeis.org/A000142/a000142.pdf

factorial.c

/**
 * Factorial algorithm - 
 *  Divide, Swing & Conquer from https://oeis.org/A000142/a000142.pdf
 **/


/**
 * check if input number is prime. 
 * Return 1 if prime, otherwise 0
 **/
char is_prime(uint64_t n)
{
    if ((n & 0x01) == 0)
        return 0;
    for (int i = 3; i*i <= n+1; i++)
        if (n % i == 0)
            return 0;
    return 1;
}

/**
 * Return the next occuring prime number 
 **/
uint64_t get_next_prime(uint64_t start_num)
{
    #define MAX_PRIME 50000
    for (int i = start_num+1; i < MAX_PRIME; i++)
        if (is_prime(i))
            return i;
    return 0;
}

uint64_t product(uint64_t *list, uint64_t len, uint64_t *index)
{
    // static uint64_t index = 0;
    if (len == 0)
        return 1;
    if (len == 1)
        return list[index[0]++];
    
    uint64_t hlen = len >> 1;
    return product(list, len - hlen, index) * product(list, hlen, index);
}

uint64_t prime_swing(uint64_t n)
{
    uint64_t count = 0;
    static uint64_t factor_list[500] = {0};

    for (uint64_t prime = 2; prime <= n && prime > 0; prime = get_next_prime(prime))
    {
        uint64_t q = n;
        uint64_t p = 1;
        do {
            q = q / prime;
            if (q & 0x01)
                p *= prime; 
        } while (q != 0);
        if (p > 1) 
        {
            factor_list[count++] = p;
            // printf("Count: %ld\t Prime: %ld\t %ld\n", count, prime, p);
        }
    }

    uint64_t index  = 0;
    return product(factor_list, count, &index);
}

uint64_t factorial(uint64_t n)
{
    if (n < 2)
        return 1;
    uint64_t f1 = factorial(n >> 1);
    return f1 * f1 * prime_swing(n);
}