cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c CDF of the noncentral chisquare                                  c
c x = the value at which the cdf is to be computed, x > 0          c
c a = degrees of freedom, a > 0                                    c
c d = noncentrality parameter, d > 0                               c
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


        implicit double precision(a-h,o-z)
        print *,'Enter the observed value x  = '
        read *, x
        print *,'Enter the value of df       = '
        read *, df
        print *,'Enter the nc parameter      = '
        read *, d
        print *, 'P(X <= x)  = ',chncdf(x, df, d)
        end
c                             
        double precision function chncdf(xx, df, elambda)
        implicit double precision (a-h,o-z)
        data one, half, zero /1.0d0, 0.5d0, 0.0d0/

        if(xx .le. zero) then
          chncdf = zero
          return
        end if

        x = half*xx
        del = half*elambda
        k = int(del)
        a = half*df + k
        gamkf = gamf(x, a)
        gamkb = gamkf

        if(del .eq. 0.0) then
          chncdf = gamkf
          return
        end if

        poikf = poipro(k, del)
        poikb = poikf
        xtermf = dexp((a-one)*dlog(x)-x-alng(a))
        xtermb = xtermf*x/a
        sum = poikf * gamkf
        remain = one - poikf
1       i = i + 1
          xtermf = xtermf*x/(a+i-one)
          gamkf = gamkf - xtermf
          poikf = poikf * del/(k+i)
          termf = poikf * gamkf
          sum = sum + termf
          error = remain * gamkf
          remain = remain - poikf
c
c Do forward and backward computation k times or until convergence 
c
        if (i .gt. k) then
          if(error .le. 1.0d-12 .or. i .gt. 1000) goto 2
          goto 1
        else
          xtermb = xtermb * (a-i+one)/x
          gamkb = gamkb + xtermb
          poikb = poikb * (k-i+one)/del
          termb = gamkb * poikb
          sum = sum + termb
          remain = remain - poikb
          if(remain .le. 1.0d-12 .or. i .gt. 1000)goto 2
          goto 1
        end if
2       chncdf = sum
        end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        double precision function alng(x)
        implicit double precision (a-h, o-z)
        double precision b(8)
        logical indx
        data b/8.33333333333333d-2, 3.33333333333333d-2,
     +         2.52380952380952d-1, 5.25606469002695d-1,
     +         1.01152306812684d0,  1.51747364915329d0,
     +         2.26948897420496d0,  3.00991738325940d0/ 
        if(x .lt. 8.0d0) then
        xx = x + 8.0d0
        indx = .true.
        else
        indx = .false.
        xx = x
        end if
c
        fterm = (xx-0.5d0)*dlog(xx) - xx + 9.1893853320467d-1
        sum = b(1)/(xx+b(2)/(xx+b(3)/(xx+b(4)/(xx+b(5)/(xx+b(6)
     +/(xx+b(7)/(xx+b(8))))))))
        alng = sum + fterm
        if(indx)  alng = alng-dlog(x+7.0d0)-dlog(x+6.0d0)-dlog
     +        (x+5.0d0)-dlog(x+4.0d0)-dlog(x+3.0d0)-dlog(x+2.0d0)
     +        -dlog(x+1.0d0)-dlog(x)
        end

                               
c                                  
        double precision  function poipro(k,el)
        implicit double precision (a-h, o-z)
        ek = k
        poipro = dexp(-el+ek*dlog(el)-alng(ek+1.0d0))
        end
c
        double precision function gamf(x,a)
        implicit double precision (a-h,o-z)
        com  = dexp(a*dlog(x)-alng(a+1.0d0)-x)
        term = 1.0d0
        sum  = 1.0d0
        one  = 1.0d0
1       term = term*x/(a+one)
        sum  = sum + term
        if (term .le. 1.0d-12) goto 2
        one = one + 1.0d0
        goto 1
2       gamf  = com*sum
        if (gamf .ge. 1.0d0) gamf = 1.0d0
        end
c