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