How to Calculate the Chi-Squared P-Value in BCB

// original code from http://www.codeproject.com/Articles/432194/How-to-Calculate-the-Chi-Squared-P-Value
// this code work on C++ Builder XE.
// Just create a console project using C++ and choose option 'Use VCL'
// delete the generate Unit1.cpp (or something) file and add this one


//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

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

#define A 15

//---------------------------------------------------------------------------

double approx_gamma(double Z)
{
    const double RECIP_E = 0.36787944117144232159552377016147;  // RECIP_E = (E^-1) = (1.0 / E)
    const double TWOPI = 6.283185307179586476925286766559;  // TWOPI = 2.0 * PI

    double D = 1.0 / (10.0 * Z);
    D = 1.0 / ((12 * Z) - D);
    D = (D + Z) * RECIP_E;
    D = pow(D, Z);
    D *= sqrt(TWOPI / Z);

    return D;
}

double mygamma(double N)
{
    /*
        The constant SQRT2PI is defined as sqrt(2.0 * PI);
        For speed the constant is already defined in decimal
        form.  However, if you wish to ensure that you achieve
        maximum precision on your own machine, you can calculate
        it yourself using (sqrt(atan(1.0) * 8.0))
    */

    //const long double SQRT2PI = sqrtl(atanl(1.0) * 8.0);
    const long double SQRT2PI = 2.5066282746310005024157652848110452530069867406099383;
    long double F = 1.0;
    long double Ck;
    long double Sum = SQRT2PI;
    long double Z = (long double)N;
    long double Sc = powl((Z + A), (Z + 0.5));
    int K;

    Sc *= expl(-1.0 * (Z + A));
    Sc /= Z;

    for(K = 1; K < A; K++)
    {
        Z++;
        Ck = powl(A - K, K - 0.5);
        Ck *= expl(A - K);
        Ck /= F;

        Sum += (Ck / Z);

        F *= (-1.0 * K);
    }

    return (double)(Sum * Sc);
}

static double igf(double S, double Z)
{
    long double Sc = (1.0 / S);
    long double Sum = 1.0;
    long double Nom = 1.0;
    long double Denom = 1.0;
    int I;

    if(Z < 0.0)
    {
        return 0.0;
    }

    Sc *= powl(Z, S);
    Sc *= expl(-Z);


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

    return Sum * Sc;
}

double chisqr(int Dof, double Cv)
{
    double K = ((double)Dof) * 0.5;
    double X = Cv * 0.5;
    double PValue;

    printf("Dof:  %i\n", Dof);
    printf("Cv:  %f\n", Cv);

    if(Cv < 0 || Dof < 1)
    {
        return 0.0;
    }

    if(Dof == 2)
    {
        return exp(-1.0 * X);
    }

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

    printf("%.28f",(1.0 - PValue));

    return 0;
}

#pragma argsused
int _tmain(int argc, _TCHAR* argv[])
{
    chisqr(255, 290.285192);

    getchar();
    return 0;
}
//---------------------------------------------------------------------------