/***************************************************
*    Program to demonstrate Chi_Square_Cumul()     *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*             C++ Version By J.-P. Moreau, Paris.  *
*                     (www.jpmoreau.fr)            *
* ------------------------------------------------ *
* SAMPLE RUN:                                      *
* How many degrees of freedom: 100                 *
* What is the range (X1,X2):                       *
*   X1: 50                                         *
*   X2: 150                                        *
* What is the table step size: 5                   *
* Summation truncation error bound: 1e-6           *
*                                                  *
*   X      Chi-Square CDF                          *
* ------------------------                         *
*   50       0.00001                               *
*   55       0.00007                               *
*   60       0.00052                               *
*   65       0.00261                               *
*   70       0.00985                               *
*   75       0.02918                               *
*   80       0.07034                               *
*   85       0.14206                               *
*   90       0.24680                               *
*   95       0.37742                               *
*  100       0.51881                               *
*  105       0.65350                               *
*  110       0.76780                               *
*  115       0.85507                               *
*  120       0.91559                               *
*  125       0.95401                               *
*  130       0.97649                               *
*  135       0.98869                               *
*  140       0.99486                               *
*  145       0.99779                               *
*  150       0.99910                               *
*                                                  *
***************************************************/
#include <stdio.h>
#include <math.h>

double e,m1,x,x1,x2,xx2,x3,y;
int    m;


/**************************************************
* Series approximation subroutine LN(X!)          *
* Accuracy better then 6 places for x>=3          *
* Accuracy better than 12 places for x>10         *
* Advantage is that very large values of the      *
* argument can be used without fear of over flow. *
* x is the input, y is the output.                *
* ----------------------------------------------- *
* Reference: CRC Math Tables.                     *
* ************************************************/
void LN_FactX()  {
  double x1;
  x1 = 1.0 / (x * x);
  y = (x + 0.5) * log(x) - x * (1 - x1 / 12 + x1 * x1 / 360 - x1 * x1 * x1 / 1260 + x1 * x1 * x1 * x1 / 1680);
  y += 0.918938533205;
}


/*****************************************************
* Chi-square function subroutine. This program takes *
* a given degree of freedom, m and value, x, and     *
* calculates the chi-square density distribution     *
* function value, y. Subroutine used: LN(X!).        *
* -------------------------------------------------- *
* Reference: Texas Instruments SR-51 owners Manual,  *
* 1974.                                              *
*****************************************************/
void Chi_Square()  {
  double c;
  //Save X
  m1 = x;
  //Perform calculation
  x = (double) (m / 2 - 1);
  //Call LN(X!) subroutine
  LN_FactX();
  x = m1;
  c = -x / 2 + (m / 2 - 1) * log(x) - (m / 2) * log(2) - y;
  y = exp(c);
}

/****************************************************************
*        Chi-square cumulative distribution subroutine          *
* ------------------------------------------------------------- *
* The program is fairly accurate and calls upon the chi-square  *
* probability density function subroutine. The input parameter  *
* is m, the number of degrees of freedom. Also required is the  *
* ordinate value, x. The subroutine returns y, the cummulative  *
* distribution integral from 0 to x. This program also requires *
* an accuracy parameter, e, to determine the level of summation.*
* Calls Chi_Square() subroutine.                                *
* ------------------------------------------------------------- *
* Reference: Hewlett-Packard statistics programs, 1974.         *
****************************************************************/
void Chi_Square_Cumul()  {
  // Labels: e100, e200
  int m2; double y1;
  y1 = 1; x2 = x; m2 = m + 2;
  x2 = x2 / m2;
  e100: y1 = y1 + x2;
  if (x2 < e) goto e200;
  m2 = m2 + 2;
  // This form is used to avoid overflow
  x2 = x2 * (x / m2);
  // Loop to continue sum
  goto e100;
  // Obtain y, the probability density function
  e200: Chi_Square();
  y = (y1 * y * 2) * (x / m);
}


void main()  {
  printf("\n");
  printf(" How many degrees of freedom: "); scanf("%d",&m);
  printf(" What is the range (X1,X2):\n");
  printf("   X1: "); scanf("%lf",&x1);
  printf("   X2: "); scanf("%lf",&x2);
  printf(" What is the table step size: "); scanf("%lf",&x3);
  printf(" Summation truncation error bound: "); scanf("%lf",&e);
  printf("\n   X       Chi-Square CDF ");
  printf("\n -------------------------\n");
  x=x1; xx2=x2;
  while (x <= xx2)  {
    Chi_Square_Cumul();
    printf("  %3.0f        %7.5f\n",x,y);
    x +=  x3;
  }
  printf("\n");
}


// End of file chi-sq.cpp