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

/*LINTLIBRARY*/

#include <math.h>
#include "z.c"
/*static char sccsfid[] = "@(#) chisq.c 5.2 (|stat) 12/08/86";*/

#define        CHI_EPSILON     0.000001    /* accuracy of critchi approximation */
#define        CHI_MAX     99999.0         /* maximum chi square value */

#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))

double pochisq ();
double critchi ();

/*#ifdef CHISQTEST
double Prob[] = { .10, .05, .01, .005, .001, -1.0 };
main ()
       {
       int     df;
       int     p;
       printf ("%-4s ", "df");
       for (p = 0; Prob[p] > 0.0; p++)
               printf ("%8.3f ", Prob[p]);
       putchar ('\n');
       for (df = 1; df < 30; df++)
               {
               printf ("%4d ", df);
               for (p = 0; Prob[p] > 0.0; p++)
                       printf ("%8.3f ", critchi (Prob[p], df));
               putchar ('\n');
               }
       }
#endif CHISQTEST*/

/*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 (x, df)
double x;       /* obtained chi-square value */
int    df;      /* degrees of freedom */
       {
       double  a, y, s;
       double  e, c, z;
       double  poz ();   /* computes probability of normal z score */
       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;
       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);
       }

/*FUNCTION critchi: compute critical chi square value to produce given p */
double

critchi (p, df)
double p;
int    df;
       {
       double  minchisq = 0.0;
       double  maxchisq = CHI_MAX;
       double  chisqval;
       
       if (p <= 0.0)
               return (maxchisq);
       else if (p >= 1.0)
               return (0.0);
       
       chisqval = df / sqrt (p);    /* fair first value */
       while (maxchisq - minchisq > CHI_EPSILON)
               {
               if (pochisq (chisqval, df) < p)
                       maxchisq = chisqval;
               else
                       minchisq = chisqval;
               chisqval = (maxchisq + minchisq) * 0.5;
               }
       return (chisqval);
       }