!****************************************************
!*  Program to demonstrate Chi_Ssquare_Cumul()      *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
!*                                                  *
!*             F90 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                               *
!*                                                  *
!****************************************************
PROGRAM CHISQ

  real*8  e,m1,x,x1,x2,xx2,x3,y
  integer m

  print *, ' '
  write(*,'(" How many degrees of freedom: ")',advance='no')
  read *, m
  print *,'What is the range (X1,X2):'
  write(*,'("   X1: ")',advance='no')
  read *, x1
  write(*,'("   X2: ")',advance='no')
  read *, x2
  write(*,'(" What is the table step size: ")',advance='no')
  read *, x3
  write(*,'(" Summation truncation error bound: ")',advance='no')
  read *, e
  print *, ' '
  print *, '   X       Chi-Square CDF '
  print *, '--------------------------'
  x=x1; xx2=x2
  do while (x <= xx2)  
    call Chi_Square_Cumul(e,m,x,x2,y)
    write(*,10)  x, y
    x = x + x3
  end do
  print *, ' '
  stop
10 format('  ',F4.0,'        ',F7.5)
end

!***************************************************
!* 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.                     *
!* *************************************************
Subroutine LN_FactX(x,y)  
  real*8 x, x1, y
  x1 = 1.d0 / (x * x)
  y = (x + 0.5d0) * dlog(x) - x * (1.d0 - x1 / 12 + x1 * x1 / 360.d0 - x1 * x1 * x1 / 1260.d0 + x1 * x1 * x1 * x1 / 1680.d0)
  y = y + 0.918938533205d0
  return
end

!******************************************************
!* 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.                                              *
!******************************************************
Subroutine Chi_Square(m,x,y)
  real*8 c, m1, x, y
  !Save X
  m1 = x
  !Perform calculation
  x =  (dfloat(m) / 2.d0 - 1.d0)
  !Call LN(X!) subroutine
  call LN_FactX(x,y)
  x = m1
  c = -x / 2 + (dfloat(m) / 2.d0 - 1.d0) * dlog(x) - (dfloat(m) / 2.d0) * dlog(2.d0) - y
  y = dexp(c)
  return
end

!*****************************************************************
!*       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.         *
!*****************************************************************
Subroutine Chi_Square_Cumul(e,m,x,x2,y)  
  ! Labels: 100, 200
  integer m, m2
  real*8 e,x,x2,y,y1
  y1 = 1; x2 = x; m2 = m + 2
  x2 = x2 / m2
100 y1 = y1 + x2
  if (x2 < e) goto 200
  m2 = m2 + 2;
  ! This form is used to avoid overflow
  x2 = x2 * (x / m2);
  ! Loop to continue sum
  goto 100
  ! Obtain y, the probability density function
200 call Chi_Square(m,x,y)
  y = (y1 * y * 2) * (x / m)
  return
end

! End of file chi-sq.f90