MODULE Complex_Incomplete_Gamma

! Code converted using TO_F90 by Alan Miller
! Date: 2002-03-31  Time: 00:10:14

! --- Written By Eric Kostlan & Dmitry Gokhman
! --- March  1986
! --- For documentation, see:
! http://www.math.utsa.edu/~gokhman/papers/igf.html

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

PRIVATE
PUBLIC  :: cdig


CONTAINS


FUNCTION cdig(alpha, x) RESULT(fn_val)
 
COMPLEX (dp), INTENT(IN)  :: alpha
COMPLEX (dp), INTENT(IN)  :: x
COMPLEX (dp)              :: fn_val

COMPLEX (dp)  :: p, q
INTEGER       :: i, ilim
REAL (dp), PARAMETER     :: zero = 0.0_dp, xlim = 1.0_dp
COMPLEX (dp), PARAMETER  :: cone = (1.0_dp, 0.0_dp)
REAL (dp), PARAMETER     :: re = (0.36787944117144232_dp, 0.0_dp)
INTEGER, PARAMETER       :: ibuf = 34

! --- If x is near the negative real axis, then shift to x=1.
IF (dnrm(x) < xlim .OR. REAL(x, KIND=dp) < zero .AND. ABS(AIMAG(x)) < xlim) THEN
  fn_val = re / cdh(alpha, cone)
  ilim = REAL(x/re, KIND=dp)
  DO  i = 0, ibuf - ilim
    CALL term(alpha, x, i, p, q)
    fn_val = fn_val + p * q
  END DO
ELSE
  fn_val = EXP(-x + alpha*LOG(x)) / cdh(alpha, x)
END IF
RETURN
END FUNCTION cdig



FUNCTION cdh(alpha, x) RESULT(fn_val)
! --- Written By Eric Kostlan & Dmitry Gokhman
! --- March  1986

COMPLEX (dp), INTENT(IN)  :: alpha
COMPLEX (dp), INTENT(IN)  :: x
COMPLEX (dp)              :: fn_val

COMPLEX (dp)  :: term, sum, cn, alpha1
INTEGER       :: i, n
REAL (dp), PARAMETER  :: one = 1.0_dp

! --- If Re(alpha-x) is too big, shift alpha.
n = REAL(alpha-x, KIND=dp)
IF (n > 0) THEN
  cn = n
  alpha1 = alpha - cn
  term = one / x
  sum = term
  DO  i = 1, n - 1
    cn = n - i
    term = term * (alpha1 + cn) / x
    sum = term + sum
  END DO
  sum = sum + term * alpha1 / cdhs(alpha1, x)
  fn_val = one / sum
ELSE
  fn_val = cdhs(alpha, x)
END IF
RETURN
END FUNCTION cdh



FUNCTION cdhs(alpha, x) RESULT(fn_val)
! --- Written By Eric Kostlan & Dmitry Gokhman
! --- March  1986

COMPLEX (dp), INTENT(IN)  :: alpha
COMPLEX (dp), INTENT(IN)  :: x
COMPLEX (dp)              :: fn_val

COMPLEX (dp)  :: p0, q0, p1, q1, r0, r1, ci, factor
INTEGER       :: i
REAL (dp), PARAMETER  :: zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp
REAL (dp), PARAMETER  :: tol1 = 1.0D+10, tol2 = 1.0D-10, error = 5.D-18
INTEGER, PARAMETER    :: ilim = 100000

q0 = one
q1 = one
p0 = x
p1 = x + one - alpha
DO  i = 1, ilim
  ci = i
  IF (p0 /= zero .AND. q0 /= zero .AND. q1 /= zero) THEN
    r0 = p0 / q0
    r1 = p1 / q1
    IF (dnrm(r0-r1) <= dnrm(r1)*error) THEN
      fn_val = r1
      RETURN
    END IF
! --------- Occasionally renormalize the sequences to avoid over(under)flow.
    IF (dnrm(p0) > tol1 .OR. dnrm(p0) < tol2 .OR. dnrm(q0) > tol1  &
          .OR. dnrm(q0) < tol2) THEN
      factor = p0 * q0
      p0 = p0 / factor
      q0 = q0 / factor
      p1 = p1 / factor
      q1 = q1 / factor
    END IF
  END IF
  p0 = x * p1 + ci * p0
  q0 = x * q1 + ci * q0
  p1 = p0 + (ci+one-alpha) * p1
  q1 = q0 + (ci+one-alpha) * q1
END DO
! --- If the peripheral routines are written correctly,
! --- the following four statements should never be executed.
WRITE(*, *) 'cdhs:  *** Warning: i >', ilim
WRITE(*, *) 'cdhs:  *** r0,r1= ', r0, r1
fn_val = half * (r0+r1)
RETURN
END FUNCTION cdhs



SUBROUTINE term(alpha, x, i, p, q)
! --- Calculate p*q = -1**i(1 - x**(alpha+i))/(alpha+i)i ! carefully.

COMPLEX (dp), INTENT(IN)   :: alpha
COMPLEX (dp), INTENT(IN)   :: x
INTEGER, INTENT(IN)        :: i
COMPLEX (dp), INTENT(OUT)  :: p
COMPLEX (dp), INTENT(OUT)  :: q

COMPLEX (dp)  :: ci, cdlx, alphai
REAL (dp), PARAMETER  :: zero = 0.0_dp, one = 1.0_dp, two = 2.0_dp
REAL (dp), PARAMETER  :: tol = 3.0D-7, xlim = 39.0_dp

IF (i == 0) q = one
ci = i
alphai = alpha + ci
IF (x == zero) THEN
  p = one / alphai
  IF (i /= 0) q = -q / ci
  RETURN
END IF
cdlx = LOG(x)

! --- If (1 - x**alphai) = -x**alphai on the computer,
! --- then change the inductive scheme to avoid overflow.
IF (REAL(alphai*cdlx, KIND=dp) > xlim .AND. i /= 0) THEN
  p = p * (alphai - one) / alphai
  q = -q * x / ci
  RETURN
END IF
IF (dnrm(alphai) > tol) THEN
  p = (one - x**alphai) / alphai
ELSE
  p = -cdlx * (one + cdlx*alphai/two)
END IF
IF (i /= 0) q = -q / ci
RETURN
END SUBROUTINE term



FUNCTION dnrm(z) RESULT(fn_val)

COMPLEX (dp), INTENT(IN)  :: z
REAL (dp)                 :: fn_val

dnrm = ABS(REAL(z, KIND=dp)) + ABS(AIMAG(z))
RETURN
END FUNCTION dnrm

END MODULE Complex_Incomplete_Gamma