mhygfx.f90

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

MODULE hygfx_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."

! Latest revision - 16 January 2002
! Corrections by Alan Miller (amiller @ bigpond.net.au)
! In 4 places HW was used without being initialized.

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
REAL (dp), PARAMETER, PRIVATE  :: pi = 3.141592653589793D0,  &
                                  el = .5772156649015329D0

CONTAINS
 

SUBROUTINE hygfx(a, b, c, x, hf)

!       ====================================================
!       Purpose: Compute hypergeometric function F(a,b,c,x)
!       Input :  a --- Parameter
!                b --- Parameter
!                c --- Parameter, c <> 0,-1,-2,...
!                x --- Argument   ( x < 1 )
!       Output:  HF --- F(a,b,c,x)
!       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)      :: c
REAL (dp), INTENT(IN OUT)  :: x
REAL (dp), INTENT(OUT)     :: hf

LOGICAL    :: l0, l1, l2, l3, l4, l5
REAL (dp)  :: a0, aa, bb, c0, c1, eps, f0, f1, g0, g1, g2, g3,   &
              ga, gabc, gam, gb, gbm, gc, gca, gcb, gcab, gm, hw,  &
              pa, pb, r, r0, r1, rm, rp, sm, sp, sp0, x1
INTEGER    :: j, k, m, nm

l0 = c == INT(c) .AND. c < 0.0
l1 = 1.0D0 - x < 1.0D-15 .AND. c - a - b <= 0.0
l2 = a == INT(a) .AND. a < 0.0
l3 = b == INT(b) .AND. b < 0.0
l4 = c - a == INT(c-a) .AND. c - a <= 0.0
l5 = c - b == INT(c-b) .AND. c - b <= 0.0
IF (l0 .OR. l1) THEN
  WRITE (*,*) 'The hypergeometric series is divergent'
  RETURN
END IF
eps = 1.0D-15
IF (x > 0.95) eps = 1.0D-8
IF (x == 0.0 .OR. a == 0.0.OR.b == 0.0) THEN
  hf = 1.0D0
  RETURN
ELSE IF (1.0D0-x == eps .AND. c-a-b > 0.0) THEN
  CALL gamma(c,gc)
  CALL gamma(c-a-b,gcab)
  CALL gamma(c-a,gca)
  CALL gamma(c-b,gcb)
  hf = gc * gcab / (gca*gcb)
  RETURN
ELSE IF (1.0D0+x <= eps .AND. ABS(c-a+b-1.0) <= eps) THEN
  g0 = SQRT(pi) * 2.0D0 ** (-a)
  CALL gamma(c,g1)
  CALL gamma(1.0D0+a/2.0-b,g2)
  CALL gamma(0.5D0+0.5*a,g3)
  hf = g0 * g1 / (g2*g3)
  RETURN
ELSE IF (l2 .OR. l3) THEN
  IF (l2) nm = INT(ABS(a))
  IF (l3) nm = INT(ABS(b))
  hf = 1.0D0
  r = 1.0D0
  DO  k = 1, nm
    r = r * (a+k-1.0D0) * (b+k-1.0D0) / (k*(c+k-1.0D0)) * x
    hf = hf + r
  END DO
  RETURN
ELSE IF (l4 .OR. l5) THEN
  IF (l4) nm = INT(ABS(c-a))
  IF (l5) nm = INT(ABS(c-b))
  hf = 1.0D0
  r = 1.0D0
  DO  k = 1, nm
    r = r * (c-a+k-1.0D0) * (c-b+k-1.0D0) / (k*(c+k-1.0D0)) * x
    hf = hf + r
  END DO
  hf = (1.0D0-x) ** (c-a-b) * hf
  RETURN
END IF
aa = a
bb = b
x1 = x
IF (x < 0.0D0) THEN
  x = x / (x-1.0D0)
  IF (c > a .AND. b < a.AND.b > 0.0) THEN
    a = bb
    b = aa
  END IF
  b = c - b
END IF
IF (x >= 0.75D0) THEN
  gm = 0.0D0
  IF (ABS(c-a-b-INT(c-a-b)) < 1.0D-15) 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.0D0
    DO  j = 1, ABS(m) - 1
      gm = gm * j
    END DO
    rm = 1.0D0
    DO  j = 1, ABS(m)
      rm = rm * j
    END DO
    f0 = 1.0D0
    r0 = 1.0D0
    r1 = 1.0D0
    sp0 = 0.d0
    sp = 0.0D0
    IF (m >= 0) THEN
      c0 = gm * gc / (gam*gbm)
      c1 = -gc * (x-1.0D0) ** m / (ga*gb*rm)
      DO  k = 1, m - 1
        r0 = r0 * (a+k-1.0D0) * (b+k-1.0) / (k*(k-m)) * (1.0-x)
        f0 = f0 + r0
      END DO
      DO  k = 1, m
        sp0 = sp0 + 1.0D0 / (a+k-1.0) + 1.0 / (b+k-1.0) - 1.0 / k
      END DO
      f1 = pa + pb + sp0 + 2.0D0 * el + LOG(1.0D0-x)
      hw = f1
      DO  k = 1, 250
        sp = sp + (1.0D0-a) / (k*(a+k-1.0)) + (1.0-b) / (k*(b+k- 1.0))
        sm = 0.0D0
        DO  j = 1, m
          sm = sm + (1.0D0-a) / ((j+k)*(a+j+k-1.0)) + 1.0 / (b+j+k -1.0)
        END DO
        rp = pa + pb + 2.0D0 * el + sp + sm + LOG(1.0D0-x)
        r1 = r1 * (a+m+k-1.0D0) * (b+m+k-1.0) / (k*(m+k)) * (1.0-x )
        f1 = f1 + r1 * rp
        IF (ABS(f1-hw) < ABS(f1)*eps) EXIT
        hw = f1
      END DO
      hf = f0 * c0 + f1 * c1
    ELSE IF (m < 0) THEN
      m = -m
      c0 = gm * gc / (ga*gb*(1.0D0-x)**m)
      c1 = -(-1) ** m * gc / (gam*gbm*rm)
      DO  k = 1, m - 1
        r0 = r0 * (a-m+k-1.0D0) * (b-m+k-1.0) / (k*(k-m)) * (1.0-x )
        f0 = f0 + r0
      END DO
      DO  k = 1, m
        sp0 = sp0 + 1.0D0 / k
      END DO
      f1 = pa + pb - sp0 + 2.0D0 * el + LOG(1.0D0-x)
      hw = f1
      DO  k = 1, 250
        sp = sp + (1.0D0-a) / (k*(a+k-1.0)) + (1.0-b) / (k*(b+k- 1.0))
        sm = 0.0D0
        DO  j = 1, m
          sm = sm + 1.0D0 / (j+k)
        END DO
        rp = pa + pb + 2.0D0 * el + sp - sm + LOG(1.0D0-x)
        r1 = r1 * (a+k-1.0D0) * (b+k-1.0) / (k*(m+k)) * (1.0-x)
        f1 = f1 + r1 * rp
        IF (ABS(f1-hw) < ABS(f1)*eps) EXIT
        hw = f1
      END DO
      hf = f0 * c0 + f1 * c1
    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)
    c0 = gc * gcab / (gca*gcb)
    c1 = gc * gabc / (ga*gb) * (1.0D0-x) ** (c-a-b)
    hf = 0.0D0
    hw = hf
    r0 = c0
    r1 = c1
    DO  k = 1, 250
      r0 = r0 * (a+k-1.0D0) * (b+k-1.0) / (k*(a+b-c+k)) * (1.0-x)
      r1 = r1 * (c-a+k-1.0D0) * (c-b+k-1.0) / (k*(c-a-b+k)) * (1.0 -x)
      hf = hf + r0 + r1
      IF (ABS(hf-hw) < ABS(hf)*eps) EXIT
      hw = hf
    END DO
    hf = hf + c0 + c1
  END IF
ELSE
  a0 = 1.0D0
  IF (c > a .AND. c < 2.0D0*a .AND. c > b .AND. c < 2.0D0*b) THEN
    a0 = (1.0D0-x) ** (c-a-b)
    a = c - a
    b = c - b
  END IF
  hf = 1.0D0
  hw = hf
  r = 1.0D0
  DO  k = 1, 250
    r = r * (a+k-1.0D0) * (b+k-1.0D0) / (k*(c+k-1.0D0)) * x
    hf = hf + r
    IF (ABS(hf-hw) <= ABS(hf)*eps) EXIT
    hw = hf
  END DO
  hf = a0 * hf
END IF
IF (x1 < 0.0D0) THEN
  x = x1
  c0 = 1.0D0 / (1.0D0-x) ** aa
  hf = c0 * hf
END IF
a = aa
b = bb
IF (k > 120) WRITE (*,5000)
RETURN
5000 FORMAT (' Warning! You should check the accuracy')
END SUBROUTINE hygfx



SUBROUTINE gamma(x, ga)

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