Find coprime binomial coefficents

coprime.c

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <gmp.h>
#include <assert.h>
#include "ff_asm_primes.h"
#define STB_DS_IMPLEMENTATION
#include "stb_ds.h"
typedef struct binomial_struct *Binomial;
struct binomial_struct
{
	uint32_t n;
	uint32_t k;
	uint32_t *prime;
	uint32_t *index;
	double bits;
}; 

double GammaApproximation(double n, double k, double base)
{
	double result = 0;
	result = (lgamma(n+1) -lgamma(n-k+1) -lgamma(k+1)) / log(base);
	return result;
}

Binomial CreateBinomial()
{
	Binomial b = malloc(sizeof(struct binomial_struct));
	b->n = 0;
	b->k = 0;
	b->prime = NULL;
	b->index = NULL;
	b->bits  = 0.0f;
	return b; 
}


void DestroyBinomial(Binomial b)
{
	if(b)
	{
		arrfree(b->prime);
		arrfree(b->index);
		free(b);
	}
}

int KummerTheorem(int n, int k, int primeNumber)
{
	int carry = 0;
	int carryCount = 0;
	int digit0 = k;
	int digit1 = n - k;
	
	while(digit0 > 0 || digit1 > 0)
	{
		int digitA = digit0 % primeNumber;
        	int digitB = digit1 % primeNumber;
        	
        	//Count carries during addition
        	if(digitA + digitB + carry >= primeNumber)
        	{
        		carryCount += 1;
        		carry = 1;
        	}
        	else
        	{
        		carry = 0;
        	}
        	digit0 /= primeNumber;
        	digit1 /= primeNumber;
	}
	return carryCount;
}

void FillBinomial(Binomial b, int n, int k)
{
	assert(b != NULL);
	b->n = n;
	b->k = k;
	b->bits = GammaApproximation(n,k,2);
	int currentPrimeIndex = 0;
	while(currentPrimeIndex < 5259 && first5239[currentPrimeIndex] <= n)
	{
		int carries = KummerTheorem(n, k, first5239[currentPrimeIndex]);
		if(carries > 0)
		{
			arrput(b->prime, first5239[currentPrimeIndex]);
			arrput(b->index, carries);
		}
		currentPrimeIndex += 1;
	}
	assert(arrlen(b->prime) == arrlen(b->index));	
}

void PrintBinomial(Binomial b)
{
	if(b)
	{
		assert(arrlen(b->prime) == arrlen(b->index));	
		printf("%3d, %3d :", b->n, b->k);
		for(size_t i = 0; i < arrlen(b->prime); i++)
		{
			printf("(%2d ^%d)",b->prime[i],b->index[i]);
		}
		printf("\n");
	}
}

int AreBinomialsCoprime(Binomial b1, Binomial b2) 
{
	if(arrlen(b1->prime) == 0 || arrlen(b2->prime) == 0){return 1;}
	size_t i = 0, j = 0;
	while(i < arrlen(b1->prime) && j < arrlen(b2->prime))
	{
		if(b1->prime[i] == b2->prime[j]){return 0;}
		else if(b1->prime[i] < b2->prime[j]){i++;}
		else{j++;}
	}
	return 1;
}
int AreAllCoprime(Binomial *list, size_t size)
{
	for(size_t i = 0; i < size; i++)
	{
		for(size_t j = i + 1; j < size; j++) {if(!AreBinomialsCoprime(list[i], list[j])) {return 0;}}
	}
	return 1;
}

void FindCoprimeCombinationsRecursive(Binomial *coefficients, size_t n, Binomial *current, size_t m, size_t start, size_t depth,int *count)
{
	if (depth == m)
	{
		if(AreAllCoprime(current, m))
		{
			(*count)++;
			double bits = 0;
			printf("Combination %d: ", *count);
			for(size_t i = 0; i < m; i++) 
			{
				printf("(%d,%d)", current[i]->n, current[i]->k);
				bits += current[i]->bits;
				if(i < m - 1){printf(", ");}
			}
			printf("(%.3f)\n", bits);
		}
		return;
	}

	for(size_t i = start; i < n; i++)
	{
		current[depth] = coefficients[i];
		FindCoprimeCombinationsRecursive(coefficients, n, current, m, i + 1, depth + 1, count);
	}
}

void FindCoprimeCombinations(Binomial *coefficients, size_t n, size_t m)
{
	if (m == 0 || m > n)
	{
		printf("Invalid combination size.\n");
		return;
	}

	Binomial *current = malloc(m * sizeof(Binomial));
	int count = 0;

	printf("\nFinding all coprime combinations of size %zu...\n", m);
	FindCoprimeCombinationsRecursive(coefficients, n, current, m, 0, 0, &count);
	printf("Total coprime combinations found: %d\n", count);

	free(current);
}

int main()
{
	int maxN = 100;
	int maxK = 0;
	int minDifference = 3;
	Binomial *coefficients = NULL;
	for(int n = 4; n <= maxN; n++)
	{
		maxK = n / 2;
		for(int k = 3; k < maxK; k++)
		{
			if(n - k >= 3)
			{
				Binomial b = CreateBinomial();
				FillBinomial(b, n, k);
				//PrintBinomial(b);
				arrput(coefficients, b);
			}
		}
	}
	FindCoprimeCombinations(coefficients, arrlen(coefficients), 4);
	for(size_t i = 0; i < arrlen(coefficients); i++){DestroyBinomial(coefficients[i]);}
	arrfree(coefficients);
}