mhygfz.f90

数值计算和数值分析在Fortran下的特殊函数库,是数值计算的必备

MODULE hygfz_func
 
! From the book "Computation of Special Functions"
!      by Shanjie Zhang and Jianming Jin
!   Copyright 1996 by John Wiley & Sons, Inc.
! The authors state:
!   "However, we give permission to the reader who purchases this book
!    to incorporate any of these programs into his or her programs
!    provided that the copyright is acknowledged."
 
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
REAL (dp), PARAMETER, PRIVATE  :: pi = 3.141592653589793_dp,  &
                                  el = .5772156649015329_dp

CONTAINS
 

SUBROUTINE hygfz(a, b, c, z, zhf)

!       ======================================================
!       Purpose: Compute the hypergeometric function for a
!                complex argument, F(a,b,c,z)
!       Input :  a --- Parameter
!                b --- Parameter
!                c --- Parameter,  c <> 0,-1,-2,...
!                z --- Complex argument
!       Output:  ZHF --- F(a,b,c,z)
!       Routines called:
!            (1) GAMMA for computing gamma function
!            (2) PSI for computing psi function
!       ======================================================

REAL (dp), INTENT(IN OUT)  :: a
REAL (dp), INTENT(IN OUT)  :: b
REAL (dp), INTENT(IN OUT)  :: c
COMPLEX (dp), INTENT(IN)   :: z
COMPLEX (dp), INTENT(OUT)  :: zhf

LOGICAL       :: l0, l1, l2, l3, l4, l5, l6
COMPLEX (dp)  :: z00, z1, zc0, zc1, zf0, zf1, zp, zp0, zr, zr0, zr1, zw
REAL (dp)     :: a0, aa, bb, ca, cb, eps, g0, g1, g2, g3,  &
                 ga, gab, gabc, gam, gb, gba, gbm, gc, gca, gcab, gcb, gcbk,  &
                 gm, pa, pac, pb, pca, rk1, rk2, rm, sj1, sj2, sm, sp, sp0,  &
                 sq, t0, w0, ws, x, y
INTEGER       :: j, k, m, mab, mcab, nca, ncb, nm

x = REAL(z, KIND=dp)
y = AIMAG(z)
eps = 1.0D-15
l0 = c == INT(c) .AND. c < 0.0_dp
l1 = ABS(1.0_dp-x) < eps .AND. y == 0.0_dp .AND. c - a - b  <= 0.0_dp
l2 = ABS(z+1.0_dp) < eps .AND. ABS(c-a+b-1.0_dp) < eps
l3 = a == INT(a) .AND. a < 0.0_dp
l4 = b == INT(b) .AND. b < 0.0_dp
l5 = c - a == INT(c-a) .AND. c - a <= 0.0_dp
l6 = c - b == INT(c-b) .AND. c - b <= 0.0_dp
aa = a
bb = b
a0 = ABS(z)
IF (a0 > 0.95_dp) eps = 1.0D-8
IF (l0 .OR. l1) THEN
  WRITE (*,*) 'The hypergeometric series is divergent'
  RETURN
END IF
IF (a0 == 0.0_dp .OR. a == 0.0_dp .OR. b == 0.0_dp) THEN
  zhf = (1.0_dp,0.0_dp)
ELSE IF (z == 1.0_dp .AND. c-a-b > 0.0_dp) THEN
  CALL gamma(c, gc)
  CALL gamma(c-a-b, gcab)
  CALL gamma(c-a, gca)
  CALL gamma(c-b, gcb)
  zhf = gc * gcab / (gca*gcb)
ELSE IF (l2) THEN
  g0 = SQRT(pi) * 2.0_dp ** (-a)
  CALL gamma(c, g1)
  CALL gamma(1.0_dp+a/2.0_dp-b, g2)
  CALL gamma(0.5_dp+0.5_dp*a, g3)
  zhf = g0 * g1 / (g2*g3)
ELSE IF (l3 .OR. l4) THEN
  IF (l3) nm = INT(ABS(a))
  IF (l4) nm = INT(ABS(b))
  zhf = (1.0_dp,0.0_dp)
  zr = (1.0_dp,0.0_dp)
  DO  k = 1, nm
    zr = zr * (a+k-1) * (b+k-1) / (k*(c+k-1)) * z
    zhf = zhf + zr
  END DO
ELSE IF (l5 .OR. l6) THEN
  IF (l5) nm = INT(ABS(c-a))
  IF (l6) nm = INT(ABS(c-b))
  zhf = (1.0_dp,0.0_dp)
  zr = (1.0_dp,0.0_dp)
  DO  k = 1, nm
    zr = zr * (c-a+k-1) * (c-b+k-1) / (k*(c+k-1)) * z
    zhf = zhf + zr
  END DO
  zhf = (1.0_dp-z) ** (c-a-b) * zhf
ELSE IF (a0 <= 1.0_dp) THEN
  IF (x < 0.0_dp) THEN
    z1 = z / (z-1.0_dp)
    IF (c > a .AND. b < a .AND. b > 0.0) THEN
      a = bb
      b = aa
    END IF
    zc0 = 1.0_dp / ((1.0_dp-z)**a)
    zhf = (1.0_dp,0.0_dp)
    zr0 = (1.0_dp,0.0_dp)
    DO  k = 1, 500
      zr0 = zr0 * (a+k-1) * (c-b+k-1) / (k*(c+k-1)) * z1
      zhf = zhf + zr0
      IF (ABS(zhf-zw) < ABS(zhf)*eps) EXIT
      zw = zhf
    END DO
    zhf = zc0 * zhf
  ELSE IF (a0 >= 0.90_dp) THEN
    gm = 0.0_dp
    mcab = INT(c-a-b+eps*SIGN(1.0_dp, c-a-b))
    IF (ABS(c-a-b-mcab) < eps) THEN
      m = INT(c-a-b)
      CALL gamma(a, ga)
      CALL gamma(b, gb)
      CALL gamma(c, gc)
      CALL gamma(a+m, gam)
      CALL gamma(b+m, gbm)
      CALL psi(a, pa)
      CALL psi(b, pb)
      IF (m /= 0) gm = 1.0_dp
      DO  j = 1, ABS(m) - 1
        gm = gm * j
      END DO
      rm = 1.0_dp
      DO  j = 1, ABS(m)
        rm = rm * j
      END DO
      zf0 = (1.0_dp,0.0_dp)
      zr0 = (1.0_dp,0.0_dp)
      zr1 = (1.0_dp,0.0_dp)
      sp0 = 0._dp
      sp = 0.0_dp
      IF (m >= 0) THEN
        zc0 = gm * gc / (gam*gbm)
        zc1 = -gc * (z-1.0_dp) ** m / (ga*gb*rm)
        DO  k = 1, m - 1
          zr0 = zr0 * (a+k-1) * (b+k-1) / (k*(k-m)) * (1._dp - z)
          zf0 = zf0 + zr0
        END DO
        DO  k = 1, m
          sp0 = sp0 + 1.0_dp / (a+k-1) + 1.0 / (b+k-1) - 1._dp / k
        END DO
        zf1 = pa + pb + sp0 + 2.0_dp * el + LOG(1.0_dp-z)
        DO  k = 1, 500
          sp = sp + (1.0_dp-a) / (k*(a+k-1)) + (1.0_dp-b) / (k*( b+k-1))
          sm = 0.0_dp
          DO  j = 1, m
            sm = sm + (1.0_dp-a) / ((j+k)*(a+j+k-1)) + 1.0_dp / (b+j+k-1)
          END DO
          zp = pa + pb + 2.0_dp * el + sp + sm + LOG(1.0_dp-z)
          zr1 = zr1 * (a+m+k-1) * (b+m+k-1) / (k*(m+k)) * (1.0_dp-z)
          zf1 = zf1 + zr1 * zp
          IF (ABS(zf1-zw) < ABS(zf1)*eps) EXIT
          zw = zf1
        END DO
        zhf = zf0 * zc0 + zf1 * zc1
      ELSE IF (m < 0) THEN
        m = -m
        zc0 = gm * gc / (ga*gb*(1.0_dp-z)**m)
        zc1 = -(-1) ** m * gc / (gam*gbm*rm)
        DO  k = 1, m - 1
          zr0 = zr0 * (a-m+k-1) * (b-m+k-1) / (k*(k-m)) * (1.0_dp-z)
          zf0 = zf0 + zr0
        END DO
        DO  k = 1, m
          sp0 = sp0 + 1.0_dp / k
        END DO
        zf1 = pa + pb - sp0 + 2.0_dp * el + LOG(1.0_dp-z)
        DO  k = 1, 500
          sp = sp + (1.0_dp-a) / (k*(a+k-1)) + (1.0_dp-b) / (k*( b+k-1))
          sm = 0.0_dp
          DO  j = 1, m
            sm = sm + 1.0_dp / (j+k)
          END DO
          zp = pa + pb + 2.0_dp * el + sp - sm + LOG(1.0_dp-z)
          zr1 = zr1 * (a+k-1._dp) * (b+k-1) / (k*(m+k)) * (1._dp - z)
          zf1 = zf1 + zr1 * zp
          IF (ABS(zf1-zw) < ABS(zf1)*eps) EXIT
          zw = zf1
        END DO
        zhf = zf0 * zc0 + zf1 * zc1
      END IF
    ELSE
      CALL gamma(a, ga)
      CALL gamma(b, gb)
      CALL gamma(c, gc)
      CALL gamma(c-a, gca)
      CALL gamma(c-b, gcb)
      CALL gamma(c-a-b, gcab)
      CALL gamma(a+b-c, gabc)
      zc0 = gc * gcab / (gca*gcb)
      zc1 = gc * gabc / (ga*gb) * (1.0_dp-z) ** (c-a-b)
      zhf = (0.0_dp,0.0_dp)
      zr0 = zc0
      zr1 = zc1
      DO  k = 1, 500
        zr0 = zr0 * (a+k-1) * (b+k-1) / (k*(a+b-c+k)) * ( 1._dp-z)
        zr1 = zr1 * (c-a+k-1) * (c-b+k-1) / (k*(c-a-b+k)) * (1.0_dp-z)
        zhf = zhf + zr0 + zr1
        IF (ABS(zhf-zw) < ABS(zhf)*eps) EXIT
        zw = zhf
      END DO
      zhf = zhf + zc0 + zc1
    END IF
  ELSE
    z00 = (1.0_dp,0.0_dp)
    IF (c-a < a .AND. c-b < b) THEN
      z00 = (1.0_dp-z) ** (c-a-b)
      a = c - a
      b = c - b
    END IF
    zhf = (1.0_dp,0._dp)
    zr = (1.0_dp,0.0_dp)
    DO  k = 1, 1500
      zr = zr * (a+k-1) * (b+k-1) / (k*(c+k-1)) * z
      zhf = zhf + zr
      IF (ABS(zhf-zw) <= ABS(zhf)*eps) EXIT
      zw = zhf
    END DO
    zhf = z00 * zhf
  END IF
ELSE IF (a0 > 1.0_dp) THEN
  mab = INT(a-b+eps*SIGN(1.0_dp, a-b))
  IF (ABS(a-b-mab) < eps .AND. a0 <= 1.1_dp) b = b + eps
  IF (ABS(a-b-mab) > eps) THEN
    CALL gamma(a, ga)
    CALL gamma(b, gb)
    CALL gamma(c, gc)
    CALL gamma(a-b, gab)
    CALL gamma(b-a, gba)
    CALL gamma(c-a, gca)
    CALL gamma(c-b, gcb)
    zc0 = gc * gba / (gca*gb*(-z)**a)
    zc1 = gc * gab / (gcb*ga*(-z)**b)
    zr0 = zc0
    zr1 = zc1
    zhf = (0.0_dp,0.0_dp)
    DO  k = 1, 500
      zr0 = zr0 * (a+k-1) * (a-c+k) / ((a-b+k)*k*z)
      zr1 = zr1 * (b+k-1) * (b-c+k) / ((b-a+k)*k*z)
      zhf = zhf + zr0 + zr1
      IF (ABS((zhf-zw)/zhf) <= eps) EXIT
      zw = zhf
    END DO
    zhf = zhf + zc0 + zc1
  ELSE
    IF (a-b < 0.0_dp) THEN
      a = bb
      b = aa
    END IF
    ca = c - a
    cb = c - b
    nca = INT(ca+eps*SIGN(1.0_dp, ca))
    ncb = INT(cb+eps*SIGN(1.0_dp, cb))
    IF (ABS(ca-nca) < eps .OR. ABS(cb-ncb) < eps) c = c + eps
    CALL gamma(a, ga)
    CALL gamma(c, gc)
    CALL gamma(c-b, gcb)
    CALL psi(a, pa)
    CALL psi(c-a, pca)
    CALL psi(a-c, pac)
    mab = INT(a-b+eps)
    zc0 = gc / (ga*(-z)**b)
    CALL gamma(a-b, gm)
    zf0 = gm / gcb * zc0
    zr = zc0
    DO  k = 1, mab - 1
      zr = zr * (b+k-1) / (k*z)
      t0 = a - b - k
      CALL gamma(t0,g0)
      CALL gamma(c-b-k,gcbk)
      zf0 = zf0 + zr * g0 / gcbk
    END DO
    IF (mab == 0) zf0 = (0.0_dp,0.0_dp)
    zc1 = gc / (ga*gcb*(-z)**a)
    sp = -2.0_dp * el - pa - pca
    DO  j = 1, mab
      sp = sp + 1.0_dp / j
    END DO
    zp0 = sp + LOG(-z)
    sq = 1.0_dp
    DO  j = 1, mab
      sq = sq * (b+j-1) * (b-c+j) / j
    END DO
    zf1 = (sq*zp0) * zc1
    zr = zc1
    rk1 = 1.0_dp
    sj1 = 0.0_dp
    DO  k = 1, 10000
      zr = zr / z
      rk1 = rk1 * (b+k-1) * (b-c+k) / (k*k)
      rk2 = rk1
      DO  j = k + 1, k + mab
        rk2 = rk2 * (b+j-1) * (b-c+j) / j
      END DO
      sj1 = sj1 + (a-1.0_dp) / (k*(a+k-1)) + (a-c-1.0_dp) / (k*(a-c+k-1))
      sj2 = sj1
      DO  j = k + 1, k + mab
        sj2 = sj2 + 1.0_dp / j
      END DO
      zp = -2.0_dp * el - pa - pac + sj2 - 1.0_dp / (k+a-c) - pi /  &
           TAN(pi*(k+a-c)) + LOG(-z)
      zf1 = zf1 + rk2 * zr * zp
      ws = ABS(zf1)
      IF (ABS((ws-w0)/ws) < eps) EXIT
      w0 = ws
    END DO
    zhf = zf0 + zf1
  END IF
END IF
a = aa
b = bb
IF (k > 150) WRITE (*,5000)
RETURN
5000 FORMAT (' Warning! You should check the accuracy')
END SUBROUTINE hygfz



SUBROUTINE gamma(x, ga)

!       ==================================================
!       Purpose: Compute gamma function