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

#include "calculus.h"


/*
	type could be 'i' (int), 'l' (long), 'f' (float) or 'd' (double)
*/

double *Moments(   void * untyped_array , int n, char type )
{

	double sum[3]={0,};
	double *moment;
	double val=0;
	int i,j;
	
	double epsilon=1e-6;
	
	int *iarray=NULL;
	long *larray=NULL;
	float *farray=NULL;
	double *darray=NULL;
	
	
	switch(type){
		
		case 'i':
			iarray=(int *)untyped_array;
			break;
		
		case 'l':
			larray=(long *)untyped_array;
			break;
		
		case 'f':
			farray=(float *)untyped_array;
			break;
		
		case 'd':
			darray=(double *)untyped_array;
			break;

		default:
			fprintf(stderr, "Moments: type of array incorrect, please use i, l, f or d\n");
			exit(4);
	}

	
	moment = (double *)calloc( 3, sizeof(double) );
	if(!moment)fprintf(stderr, "Moments: cannot allocate moment, bye\n"), exit(3);
	
	for(i=0;i<n;i++){
	
		switch(type){
		
			case 'i':
				val = iarray[i];
				break;
		
			case 'l':
				val = larray[i];
				break;
		
			case 'f':
				val = farray[i];
				break;
		
			case 'd':
				val = darray[i];
				break;
		}
	
		for(j=0;j<3;j++)
			sum[j] += (double )pow( val, j+1.0 );
	}
	

	/*
		mean
	*/
	moment[0] = sum[0] / n;
	
	/*
		variance
	*/
	
	if( fabs( sum[0]*sum[0] - n*sum[1] ) < epsilon )
		moment[1] = 0;
	else
		moment[1] = (sum[1]/n - pow(moment[0], 2.0));

	/*
		skewness
	*/
	if(moment[1] == 0)
		moment[2] = 0;
	else{
		moment[2] = sum[2]/n - 3.0*moment[0]*moment[1] - pow(moment[0],3);
		moment[2] /= pow( moment[1], 1.5 );
	}


	/*
		Correct for bias
	*/
	moment[1] *= ( n / (n-1.0) );
	moment[2] *= sqrt( n*(n-1.0) ) / (n-2.0); 
	

	return moment;
}

/*
	All to compute n choose n1
*/
long double LOG_factoriel_step(int n, int steps)
{

	if (steps > n)
		printf("step should be smaller than n, bye\n"), exit(2);
	if (steps == 0)
		return 0;
	return logl(n) + LOG_factoriel_step(n - 1, steps - 1);

}
long double LOG_factoriel(int n)
{

	if (n <= 1)
		return 0;
	return logl(n) + LOG_factoriel(n - 1);
}
long double log_N_choose_n(int n, int n1)
{

	long double     f = 0;
	int             n2 = n - n1;

	if (n1 > n2)
		f = LOG_factoriel_step(n, n - n1) - LOG_factoriel(n2);
	else
		f = LOG_factoriel_step(n, n - n2) - LOG_factoriel(n1);

	return f;
}

long N_choose_n(int n, int n1)
{
	return lround(exp(log_N_choose_n(n, n1)));
}