f_dst_gammai_core Function

elemental function f_dst_gammai_core(a, x) result(p)

! ==== Description
!! Incomplete gamma function. Needed by gamma and chi-squared cdf.
!! Uses Fortran 2008+ intrinsics.

! ==== Declarations
  real(wp), intent(in)   :: a, x
  real(wp)               :: p
  real(wp)               :: sum, term, lngamma_a
  real(wp)               :: ap, del, b, c, d, h
  real(wp)   , parameter :: eps = 1.0e-12_wp     !! convergence threshold
  real(wp)   , parameter :: fpmin = 1.0e-30_wp   !! small number to prevent division by zero
  integer(i4), parameter :: i_max = c_bisect_i          !! max. number of iterations
  integer(i4)            :: i

! ==== Instructions
  if (x .lt. 0.0_wp .or. a .le. 0.0_wp) then
     p = 0.0_wp
     return
  endif

  if (x .eq. 0.0_wp) then
     p = 0.0_wp
     return
  endif

  lngamma_a = log(gamma(a))

  if (x .lt. a + 1.0_wp) then
     ! use series expansion
     sum = 1.0_wp / a
     term = sum
     ap = a
     do
        ap = ap + 1.0_wp
        term = term * x / ap
        sum = sum + term
        if (abs(term) .lt. abs(sum) * eps) exit
     enddo
     p = sum * exp(-x + a * log(x) - lngamma_a)
  else
     ! use continued fraction expansion
     b = x + 1.0_wp - a
     c = 1.0_wp / fpmin
     d = 1.0_wp / b
     h = d
     i = 1
     do
       i = i + 1
       if (mod(i, 2) .eq. 0) then
          del = (i / 2) * (a - (i / 2))
       else
          del = -((i - 1) / 2) * ((i - 1) / 2)
       endif
       b = b + 2.0_wp
       d = del * d + b
       if (abs(d) .lt. fpmin) d = fpmin
       c = b + del / c
       if (abs(c) .lt. fpmin) c = fpmin
       d = 1.0_wp / d
       h = h * d * c
       if (abs(d * c - 1.0_wp) .lt. eps) exit
    enddo
    p = 1.0_wp - exp(-x + a * log(x) - lngamma_a) * h
  endif

end function f_dst_gammai_core