How to Calculate the Chi-Squared P-Value

/*
    Implementation of the Chi-Square Distribution &
	Incomplete Gamma Function in C

	Written By Jacob Wells
	July 31, 2012
    Based on the formulas found here:

    Wikipedia - Incomplete Gamma Function -> Evaluation formulae -> Connection with Kummer's confluent hypergeometric function
    http://en.wikipedia.org/wiki/Regularized_Gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function

    Wikipedia - Chi-squared Distribution -> Cumulative distribution function
    http://en.wikipedia.org/wiki/Chi-squared_distribution#Cumulative_distribution_function

    These functions are placed in the Public Domain, and may be used by anyone, anywhere, for any reason, absolutely free of charge.

*/


#include <stdio.h>
#include <math.h>
#include "chisqr.h"
#include "gamma.h"


static double igf(double S, double Z);


double chisqr(int Dof, double Cv)
{
    printf("Dof:  %i\n", Dof);
    printf("Cv:  %f\n", Cv);
    if(Cv < 0 || Dof < 1)
    {
        return 0.0;
    }
	double K = ((double)Dof) * 0.5;
	double X = Cv * 0.5;
	if(Dof == 2)
	{
		return exp(-1.0 * X);
	}

	double PValue = igf(K, X);
    if(isnan(PValue) || isinf(PValue) || PValue <= 1e-8)
    {
        return 1e-14;
    }
	PValue /= approx_gamma(K);
    //PValue /= gamma(K);  // First verion in gamma.c
	//PValue /= tgamma(K);  // <math.h> implementation of the gamma function (slightly more accurate & probably faster)

	return (1.0 - PValue);

}

/*
	Incomplete Gamma Function

	I'll admit, I'm still not sure how many Iterations
	are needed to get keep the value always accurate.
	I had some problems getting accurate readings either
	over or under 200 iterations.  I would suggest doing
	some trials to make sure it works as you need.  As
	it is now, it works just fine.
*/
static double igf(double S, double Z)
{
	if(Z < 0.0)
	{
		return 0.0;
	}
	long double Sc = (1.0 / S);
	Sc *= powl(Z, S);
	Sc *= expl(-Z);

	long double Sum = 1.0;
	long double Nom = 1.0;
	long double Denom = 1.0;

	for(int I = 0; I < 200; I++) // 200
	{
		Nom *= Z;
		S++;
		Denom *= S;
		Sum += (Nom / Denom);
	}

	return Sum * Sc;
}