/* This code was taken from http://www.fourmilab.ch/random/ where it states that:

   This software is in the public domain. Permission to use, copy, modify, and distribute 
   this software and its documentation for any purpose and without fee is hereby granted, 
   without any conditions or restrictions. This software is provided “as is” without 
   express or implied warranty. */

/*

    Compute probability of measured Chi Square value.
    
    This code was developed by Gary Perlman of the Wang
    Institute (full citation below) and has been minimally
    modified for use in this program.
    
*/

#include <math.h>

/*HEADER
	Module:       z.c
	Purpose:      compute approximations to normal z distribution probabilities
	Programmer:   Gary Perlman
	Organization: Wang Institute, Tyngsboro, MA 01879
	Copyright:    none
	Tabstops:     4
*/

#define	Z_MAX          6.0            /* maximum meaningful z value */

/*FUNCTION poz: probability of normal z value */
/*ALGORITHM
	Adapted from a polynomial approximation in:
		Ibbetson D, Algorithm 209
		Collected Algorithms of the CACM 1963 p. 616
	Note:
		This routine has six digit accuracy, so it is only useful for absolute
		z values < 6.  For z values >= to 6.0, poz() returns 0.0.
*/
static double        /*VAR returns cumulative probability from -oo to z */
poz(const double z)  /*VAR normal z value */
{
    double y, x, w;

    if (z == 0.0) {
    	x = 0.0;
    } else {
	y = 0.5 * fabs(z);
	if (y >= (Z_MAX * 0.5)) {
    	    x = 1.0;
	} else if (y < 1.0) {
	   w = y * y;
	   x = ((((((((0.000124818987 * w
		   -0.001075204047) * w +0.005198775019) * w
		   -0.019198292004) * w +0.059054035642) * w
		   -0.151968751364) * w +0.319152932694) * w
		   -0.531923007300) * w +0.797884560593) * y * 2.0;
	} else {
	    y -= 2.0;
	    x = (((((((((((((-0.000045255659 * y
		    +0.000152529290) * y -0.000019538132) * y
		    -0.000676904986) * y +0.001390604284) * y
		    -0.000794620820) * y -0.002034254874) * y
		    +0.006549791214) * y -0.010557625006) * y
		    +0.011630447319) * y -0.009279453341) * y
		    +0.005353579108) * y -0.002141268741) * y
		    +0.000535310849) * y +0.999936657524;
    	}
    }
    return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
}

/*
	Module:       chisq.c
	Purpose:      compute approximations to chisquare distribution probabilities
	Contents:     pochisq()
	Uses:         poz() in z.c (Algorithm 209)
	Programmer:   Gary Perlman
	Organization: Wang Institute, Tyngsboro, MA 01879
	Copyright:    none
	Tabstops:     4
*/

#define	LOG_SQRT_PI     0.5723649429247000870717135 /* log (sqrt (pi)) */
#define	I_SQRT_PI       0.5641895835477562869480795 /* 1 / sqrt (pi) */
#define	BIGX           20.0         /* max value to represent exp (x) */
#define	ex(x)             (((x) < -BIGX) ? 0.0 : exp(x))

/*FUNCTION pochisq: probability of chi sqaure value */
/*ALGORITHM Compute probability of chi square value.
	Adapted from:
		Hill, I. D. and Pike, M. C.  Algorithm 299
		Collected Algorithms for the CACM 1967 p. 243
	Updated for rounding errors based on remark in
		ACM TOMS June 1985, page 185
*/

double pochisq(
    	const double ax,    /* obtained chi-square value */
     	const int df	    /* degrees of freedom */
     	)
{
    double x = ax;
    double a, y, s;
    double e, c, z;
    int even;	    	    /* true if df is an even number */

    if (x <= 0.0 || df < 1) {
    	return 1.0;
    }

    a = 0.5 * x;
    even = (2 * (df / 2)) == df;
    y = 0.0;
    if (df > 1) {
    	y = ex(-a);
    }
    s = (even ? y : (2.0 * poz(-sqrt(x))));
    if (df > 2) {
	x = 0.5 * (df - 1.0);
	z = (even ? 1.0 : 0.5);
	if (a > BIGX) {
    	    e = (even ? 0.0 : LOG_SQRT_PI);
    	    c = log(a);
    	    while (z <= x) {
		e = log(z) + e;
		s += ex(c * z - a - e);
		z += 1.0;
    	    }
    	    return (s);
    	} else {
	    e = (even ? 1.0 : (I_SQRT_PI / sqrt(a)));
	    c = 0.0;
	    while (z <= x) {
		    e = e * (a / z);
		    c = c + e;
		    z += 1.0;
    	    }
	    return (c * y + s);
    	}
    } else {
    	return s;
    }
}