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

//gamma function using Lanczos approximation formula
//output result in log base e
//use exp() to convert back
//has a nice side effect: can store large values in small [power of e] form
double logGamma(double x)
{
    double tmp = (x-0.5) * log(x+4.5) - (x+4.5);
	double ser = 1.0 + 76.18009173     / (x+0) - 86.50532033    / (x+1)
	                 + 24.01409822     / (x+2) -  1.231739516   / (x+3)
	                 +  0.00120858003  / (x+4) -  0.00000536382 / (x+5);
	return tmp + log(ser * sqrt(2*M_PI) );	
}

//result from logGamma() are actually (n-1)!
double combination(int n, int r)
{
	return exp(logGamma(n+1)-( logGamma(r+1) + logGamma(n-r+1) ));
}

//primitive hamming weight counter
int hWeight(int x)
{
	int count, y;
	for (count=0, y=x; y; count++)
		y &= y-1; 
	return count;
}

//-------------------------------------------------------------------------------------
//recursively find the previous group's "hamming weight member count" and sum them
int rCummGroupCount(int bitsize, int hw)
{
	if (hw <= 0 || hw == bitsize) 
		return 1;
	else
		return round(combination(bitsize, hw)) + rCummGroupCount(bitsize,hw-1);
}
//-------------------------------------------------------------------------------------

int main(int argc, char* argv[])
{
	int bitsize = 4, integer = 14;
	int hw = hWeight(integer);
	int groupStartIndex = rCummGroupCount(bitsize,hw-1);
	printf("bitsize: %d\n", bitsize);
	printf("integer: %d  hamming weight: %d\n", integer, hw);
	printf("group start index: %d\n", groupStartIndex);
}