************80
! HYGFZ computes the hypergeometric function F(a,b,c,x) for complex argument.
Licensing:
This routine is copyrighted by Shanjie Zhang and Jianming Jin. However,
they give permission to incorporate this routine into a user program
provided that the copyright is acknowledged.
Modified:
03 August 2012
Author:
Shanjie Zhang, Jianming Jin
Reference:
Shanjie Zhang, Jianming Jin,
Computation of Special Functions,
Wiley, 1996,
ISBN: 0-471-11963-6,
LC: QA351.C45.
Parameters:
Input, real ( kind = 8 ) A, B, C, parameters.
Input, complex ( kind = 8 ) Z, the argument.
Output, complex ( kind = 8 ) ZHF, the value of F(a,b,c,z).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=8) | :: | a | ||||
| real(kind=8) | :: | b | ||||
| real(kind=8) | :: | c | ||||
| complex(kind=8) | :: | z | ||||
| complex(kind=8) | :: | zhf |
subroutine hygfz ( a, b, c, z, zhf ) !*****************************************************************************80 ! !! HYGFZ computes the hypergeometric function F(a,b,c,x) for complex argument. ! ! Licensing: ! ! This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, ! they give permission to incorporate this routine into a user program ! provided that the copyright is acknowledged. ! ! Modified: ! ! 03 August 2012 ! ! Author: ! ! Shanjie Zhang, Jianming Jin ! ! Reference: ! ! Shanjie Zhang, Jianming Jin, ! Computation of Special Functions, ! Wiley, 1996, ! ISBN: 0-471-11963-6, ! LC: QA351.C45. ! ! Parameters: ! ! Input, real ( kind = 8 ) A, B, C, parameters. ! ! Input, complex ( kind = 8 ) Z, the argument. ! ! Output, complex ( kind = 8 ) ZHF, the value of F(a,b,c,z). ! implicit none real ( kind = 8 ) a real ( kind = 8 ) a0 real ( kind = 8 ) aa real ( kind = 8 ) b real ( kind = 8 ) bb real ( kind = 8 ) c real ( kind = 8 ) ca real ( kind = 8 ) cb real ( kind = 8 ) el real ( kind = 8 ) eps real ( kind = 8 ) g0 real ( kind = 8 ) g1 real ( kind = 8 ) g2 real ( kind = 8 ) g3 real ( kind = 8 ) ga real ( kind = 8 ) gab real ( kind = 8 ) gabc real ( kind = 8 ) gam real ( kind = 8 ) gb real ( kind = 8 ) gba real ( kind = 8 ) gbm real ( kind = 8 ) gc real ( kind = 8 ) gca real ( kind = 8 ) gcab real ( kind = 8 ) gcb real ( kind = 8 ) gcbk real ( kind = 8 ) gm integer ( kind = 4 ) j integer ( kind = 4 ) k logical l0 logical l1 logical l2 logical l3 logical l4 logical l5 logical l6 integer ( kind = 4 ) m integer ( kind = 4 ) mab integer ( kind = 4 ) mcab integer ( kind = 4 ) nca integer ( kind = 4 ) ncb integer ( kind = 4 ) nm real ( kind = 8 ) pa real ( kind = 8 ) pac real ( kind = 8 ) pb real ( kind = 8 ) pca real ( kind = 8 ) pi real ( kind = 8 ) rk1 real ( kind = 8 ) rk2 real ( kind = 8 ) rm real ( kind = 8 ) sj1 real ( kind = 8 ) sj2 real ( kind = 8 ) sm real ( kind = 8 ) sp real ( kind = 8 ) sp0 real ( kind = 8 ) sq real ( kind = 8 ) t0 real ( kind = 8 ) w0 real ( kind = 8 ) ws real ( kind = 8 ) x real ( kind = 8 ) y complex ( kind = 8 ) z complex ( kind = 8 ) z00 complex ( kind = 8 ) z1 complex ( kind = 8 ) zc0 complex ( kind = 8 ) zc1 complex ( kind = 8 ) zf0 complex ( kind = 8 ) zf1 complex ( kind = 8 ) zhf complex ( kind = 8 ) zp complex ( kind = 8 ) zp0 complex ( kind = 8 ) zr complex ( kind = 8 ) zr0 complex ( kind = 8 ) zr1 complex ( kind = 8 ) zw x = real ( z, kind = 8 ) y = imag ( z ) eps = 1.0D-15 l0 = c == int ( c ) .and. c < 0.0D+00 l1 = abs ( 1.0D+00 - x ) < eps .and. y == 0.0D+00 .and. & c - a - b <= 0.0D+00 l2 = abs ( z + 1.0D+00 ) < eps .and. & abs ( c - a + b - 1.0D+00 ) < eps l3 = a == int ( a ) .and. a < 0.0D+00 l4 = b == int ( b ) .and. b < 0.0D+00 l5 = c - a == int ( c - a ) .and. c - a <= 0.0D+00 l6 = c - b == int ( c - b ) .and. c - b <= 0.0D+00 aa = a bb = b a0 = abs ( z ) if ( 0.95D+00 < a0 ) then eps = 1.0D-08 end if pi = 3.141592653589793D+00 el = 0.5772156649015329D+00 if ( l0 .or. l1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HYGFZ - Fatal error!' write ( *, '(a)' ) ' The hypergeometric series is divergent.' stop end if if ( a0 == 0.0D+00 .or. a == 0.0D+00 .or. b == 0.0D+00 ) then zhf = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) else if ( z == 1.0D+00.and. 0.0D+00 < c - a - b ) then call gammaf ( c, gc ) call gammaf ( c - a - b, gcab ) call gammaf ( c - a, gca ) call gammaf ( c - b, gcb ) zhf = gc * gcab / ( gca * gcb ) else if ( l2 ) then g0 = sqrt ( pi ) * 2.0D+00 ** ( - a ) call gammaf ( c, g1 ) call gammaf ( 1.0D+00 + a / 2.0D+00 - b, g2 ) call gammaf ( 0.5D+00 + 0.5D+00 * a, g3 ) zhf = g0 * g1 / ( g2 * g3 ) else if ( l3 .or. l4 ) then if ( l3 ) then nm = int ( abs ( a ) ) end if if ( l4 ) then nm = int ( abs ( b ) ) end if zhf = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) zr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) do k = 1, nm zr = zr * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( c + k - 1.0D+00 ) ) * z zhf = zhf + zr end do else if ( l5 .or. l6 ) then if ( l5 ) then nm = int ( abs ( c - a ) ) end if if ( l6 ) then nm = int ( abs ( c - b ) ) end if zhf = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) zr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) do k = 1, nm zr = zr * ( c - a + k - 1.0D+00 ) * ( c - b + k - 1.0D+00 ) & / ( k * ( c + k - 1.0D+00 ) ) * z zhf = zhf + zr end do zhf = ( 1.0D+00 - z ) ** ( c - a - b ) * zhf else if ( a0 <= 1.0D+00 ) then if ( x < 0.0D+00 ) then z1 = z / ( z - 1.0D+00 ) if ( a < c .and. b < a .and. 0.0D+00 < b ) then a = bb b = aa end if zc0 = 1.0D+00 / ( ( 1.0D+00 - z ) ** a ) zhf = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) zr0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) do k = 1, 500 zr0 = zr0 * ( a + k - 1.0D+00 ) * ( c - b + k - 1.0D+00 ) & / ( k * ( c + k - 1.0D+00 ) ) * z1 zhf = zhf + zr0 if ( abs ( zhf - zw ) < abs ( zhf ) * eps ) then exit end if zw = zhf end do zhf = zc0 * zhf else if ( 0.90D+00 <= a0 ) then gm = 0.0D+00 mcab = int ( c - a - b + eps * sign ( 1.0D+00, c - a - b ) ) if ( abs ( c - a - b - mcab ) < eps ) then m = int ( c - a - b ) call gammaf ( a, ga ) call gammaf ( b, gb ) call gammaf ( c, gc ) call gammaf ( a + m, gam ) call gammaf ( b + m, gbm ) call psi ( a, pa ) call psi ( b, pb ) if ( m /= 0 ) then gm = 1.0D+00 end if do j = 1, abs ( m ) - 1 gm = gm * j end do rm = 1.0D+00 do j = 1, abs ( m ) rm = rm * j end do zf0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) zr0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) zr1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) sp0 = 0.0D+00 sp = 0.0D+00 if ( 0 <= m ) then zc0 = gm * gc / ( gam * gbm ) zc1 = - gc * ( z - 1.0D+00 ) ** m / ( ga * gb * rm ) do k = 1, m - 1 zr0 = zr0 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( k - m ) ) * ( 1.0D+00 - z ) zf0 = zf0 + zr0 end do do k = 1, m sp0 = sp0 + 1.0D+00 / ( a + k - 1.0D+00 ) & + 1.0D+00 / ( b + k - 1.0D+00 ) - 1.0D+00 / k end do zf1 = pa + pb + sp0 + 2.0D+00 * el + log ( 1.0D+00 - z ) do k = 1, 500 sp = sp + ( 1.0D+00 - a ) & / ( k * ( a + k - 1.0D+00 ) ) + ( 1.0D+00 - b ) & / ( k * ( b + k - 1.0D+00 ) ) sm = 0.0D+00 do j = 1, m sm = sm + ( 1.0D+00 - a ) / ( ( j + k ) & * ( a + j + k - 1.0D+00 ) ) & + 1.0D+00 / ( b + j + k - 1.0D+00 ) end do zp = pa + pb + 2.0D+00 * el + sp + sm + log ( 1.0D+00 - z ) zr1 = zr1 * ( a + m + k - 1.0D+00 ) & * ( b + m + k - 1.0D+00 ) / ( k * ( m + k ) ) & * ( 1.0D+00 - z ) zf1 = zf1 + zr1 * zp if ( abs ( zf1 - zw ) < abs ( zf1 ) * eps ) then exit end if zw = zf1 end do zhf = zf0 * zc0 + zf1 * zc1 else if ( m < 0 ) then m = - m zc0 = gm * gc / ( ga * gb * ( 1.0D+00 - z ) ** m ) zc1 = - ( - 1.0D+00 ) ** m * gc / ( gam * gbm * rm ) do k = 1, m - 1 zr0 = zr0 * ( a - m + k - 1.0D+00 ) & * ( b - m + k - 1.0D+00 ) / ( k * ( k - m ) ) & * ( 1.0D+00 - z ) zf0 = zf0 + zr0 end do do k = 1, m sp0 = sp0 + 1.0D+00 / k end do zf1 = pa + pb - sp0 + 2.0D+00 * el + log ( 1.0D+00 - z ) do k = 1, 500 sp = sp + ( 1.0D+00 - a ) / ( k * ( a + k - 1.0D+00 ) ) & + ( 1.0D+00 - b ) / ( k * ( b + k - 1.0D+00 ) ) sm = 0.0D+00 do j = 1, m sm = sm + 1.0D+00 / ( j + k ) end do zp = pa + pb + 2.0D+00 * el + sp - sm + log ( 1.0D+00 - z ) zr1 = zr1 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( m + k ) ) * ( 1.0D+00 - z ) zf1 = zf1 + zr1 * zp if ( abs ( zf1 - zw ) < abs ( zf1 ) * eps ) then exit end if zw = zf1 end do zhf = zf0 * zc0 + zf1 * zc1 end if else call gammaf ( a, ga ) call gammaf ( b, gb ) call gammaf ( c, gc ) call gammaf ( c - a, gca ) call gammaf ( c - b, gcb ) call gammaf ( c - a - b, gcab ) call gammaf ( a + b - c, gabc ) zc0 = gc * gcab / ( gca * gcb ) zc1 = gc * gabc / ( ga * gb ) * ( 1.0D+00 - z ) ** ( c - a - b ) zhf = cmplx ( 0.0D+00, 0.0D+00, kind = 8 ) zr0 = zc0 zr1 = zc1 do k = 1, 500 zr0 = zr0 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( a + b - c + k ) ) * ( 1.0D+00 - z ) zr1 = zr1 * ( c - a + k - 1.0D+00 ) & * ( c - b + k - 1.0D+00 ) / ( k * ( c - a - b + k ) ) & * ( 1.0D+00 - z ) zhf = zhf + zr0 + zr1 if ( abs ( zhf - zw ) < abs ( zhf ) * eps ) then exit end if zw = zhf end do zhf = zhf + zc0 + zc1 end if else z00 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) if ( c - a < a .and. c - b < b ) then z00 = ( 1.0D+00 - z ) ** ( c - a - b ) a = c - a b = c - b end if zhf = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) zr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 ) do k = 1, 1500 zr = zr * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( c + k - 1.0D+00 ) ) * z zhf = zhf + zr if ( abs ( zhf - zw ) <= abs ( zhf ) * eps ) then exit end if zw = zhf end do zhf = z00 * zhf end if else if ( 1.0D+00 < a0 ) then mab = int ( a - b + eps * sign ( 1.0D+00, a - b ) ) if ( abs ( a - b - mab ) < eps .and. a0 <= 1.1D+00 ) then b = b + eps end if if ( eps < abs ( a - b - mab ) ) then call gammaf ( a, ga ) call gammaf ( b, gb ) call gammaf ( c, gc ) call gammaf ( a - b, gab ) call gammaf ( b - a, gba ) call gammaf ( c - a, gca ) call gammaf ( c - b, gcb ) zc0 = gc * gba / ( gca * gb * ( - z ) ** a ) zc1 = gc * gab / ( gcb * ga * ( - z ) ** b ) zr0 = zc0 zr1 = zc1 zhf = cmplx ( 0.0D+00, 0.0D+00 ) do k = 1, 500 zr0 = zr0 * ( a + k - 1.0D+00 ) * ( a - c + k ) & / ( ( a - b + k ) * k * z ) zr1 = zr1 * ( b + k - 1.0D+00 ) * ( b - c + k ) & / ( ( b - a + k ) * k * z ) zhf = zhf + zr0 + zr1 if ( abs ( ( zhf - zw ) / zhf ) <= eps ) then exit end if zw = zhf end do zhf = zhf + zc0 + zc1 else if ( a - b < 0.0D+00 ) then a = bb b = aa end if ca = c - a cb = c - b nca = int ( ca + eps * sign ( 1.0D+00, ca ) ) ncb = int ( cb + eps * sign ( 1.0D+00, cb ) ) if ( abs ( ca - nca ) < eps .or. abs ( cb - ncb ) < eps ) then c = c + eps end if call gammaf ( a, ga ) call gammaf ( c, gc ) call gammaf ( 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 gammaf ( a - b, gm ) zf0 = gm / gcb * zc0 zr = zc0 do k = 1, mab - 1 zr = zr * ( b + k - 1.0D+00 ) / ( k * z ) t0 = a - b - k call gammaf ( t0, g0 ) call gammaf ( c - b - k, gcbk ) zf0 = zf0 + zr * g0 / gcbk end do if ( mab == 0 ) then zf0 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 ) end if zc1 = gc / ( ga * gcb * ( - z ) ** a ) sp = -2.0D+00 * el - pa - pca do j = 1, mab sp = sp + 1.0D+00 / j end do zp0 = sp + log ( - z ) sq = 1.0D+00 do j = 1, mab sq = sq * ( b + j - 1.0D+00 ) * ( b - c + j ) / j end do zf1 = ( sq * zp0 ) * zc1 zr = zc1 rk1 = 1.0D+00 sj1 = 0.0D+00 do k = 1, 10000 zr = zr / z rk1 = rk1 * ( b + k - 1.0D+00 ) * ( b - c + k ) / ( k * k ) rk2 = rk1 do j = k + 1, k + mab rk2 = rk2 * ( b + j - 1.0D+00 ) * ( b - c + j ) / j end do sj1 = sj1 + ( a - 1.0D+00 ) / ( k * ( a + k - 1.0D+00 ) ) & + ( a - c - 1.0D+00 ) / ( k * ( a - c + k - 1.0D+00 ) ) sj2 = sj1 do j = k + 1, k + mab sj2 = sj2 + 1.0D+00 / j end do zp = -2.0D+00 * el - pa - pac + sj2 - 1.0D+00 / ( 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 ) then exit end if w0 = ws end do zhf = zf0 + zf1 end if end if a = aa b = bb if ( 150 < k ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HYGFZ - Warning!' write ( *, '(a)' ) ' The solution returned may have low accuracy.' end if return end subroutine hygfz