FUNCTION gammad(x, p) RESULT(fn_val)
! ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3
! Computation of the Incomplete Gamma Integral
! Auxiliary functions required: ALOGAM = logarithm of the gamma
! function, and ALNORM = algorithm AS66
! ELF90-compatible version by Alan Miller
! Latest revision - 27 October 2000
! N.B. Argument IFAULT has been removed
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
REAL (dp), INTENT(IN) :: x, p
REAL (dp) :: fn_val
! Local variables
REAL (dp) :: pn1, pn2, pn3, pn4, pn5, pn6, arg, c, rn, a, b, an
REAL (dp), PARAMETER :: zero = 0.d0, one = 1.d0, two = 2.d0, &
oflo = 1.d+37, three = 3.d0, nine = 9.d0, &
tol = 1.d-14, xbig = 1.d+8, plimit = 1000.d0, &
elimit = -88.d0
! EXTERNAL alogam, alnorm
fn_val = zero
! Check that we have valid values for X and P
IF (p <= zero .OR. x < zero) THEN
WRITE(*, *) 'AS239: Either p <= 0 or x < 0'
RETURN
END IF
IF (x == zero) RETURN
! Use a normal approximation if P > PLIMIT
IF (p > plimit) THEN
pn1 = three * SQRT(p) * ((x / p) ** (one / three) + one /(nine * p) - one)
fn_val = alnorm(pn1, .false.)
RETURN
END IF
! If X is extremely large compared to P then set fn_val = 1
IF (x > xbig) THEN
fn_val = one
RETURN
END IF
IF (x <= one .OR. x < p) THEN
! Use Pearson's series expansion.
! (Note that P is not large enough to force overflow in ALOGAM).
! No need to test IFAULT on exit since P > 0.
arg = p * LOG(x) - x - alogam(p + one, ifault)
c = one
fn_val = one
a = p
40 a = a + one
c = c * x / a
fn_val = fn_val + c
IF (c > tol) GO TO 40
arg = arg + LOG(fn_val)
fn_val = zero
IF (arg >= elimit) fn_val = EXP(arg)
ELSE
! Use a continued fraction expansion
arg = p * LOG(x) - x - alogam(p, ifault)
a = one - p
b = a + x + one
c = zero
pn1 = one
pn2 = x
pn3 = x + one
pn4 = x * b
fn_val = pn3 / pn4
60 a = a + one
b = b + two
c = c + one
an = a * c
pn5 = b * pn3 - an * pn1
pn6 = b * pn4 - an * pn2
IF (ABS(pn6) > zero) THEN
rn = pn5 / pn6
IF (ABS(fn_val - rn) <= MIN(tol, tol * rn)) GO TO 80
fn_val = rn
END IF
pn1 = pn3
pn2 = pn4
pn3 = pn5
pn4 = pn6
IF (ABS(pn5) >= oflo) THEN
! Re-scale terms in continued fraction if terms are large
pn1 = pn1 / oflo
pn2 = pn2 / oflo
pn3 = pn3 / oflo
pn4 = pn4 / oflo
END IF
GO TO 60
80 arg = arg + LOG(fn_val)
fn_val = one
IF (arg >= elimit) fn_val = one - EXP(arg)
END IF
RETURN
END FUNCTION gammad