subroutine hyper_2f1_real ( a, b, c, x, hf ) c*********************************************************************72 c cc hyper_2f1_real() computes the hypergeometric function 2F1(a,b,c;x). c c Licensing: c c This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, c they give permission to incorporate this routine into a user program c provided that the copyright is acknowledged. c c Modified: c c 04 April 2012 c c Author: c c Shanjie Zhang, Jianming Jin c c Reference: c c Shanjie Zhang, Jianming Jin, c Computation of Special Functions, c Wiley, 1996, c ISBN: 0-471-11963-6, c LC: QA351.C45. c c Input: c c double precision A, B, C, parameters. c c double precision X, the argument. c c Output: c c double precision HF, the function value. c implicit double precision (a-h,o-z) double precision a double precision b double precision c double precision c0 double precision el double precision eps double precision g0 double precision g1 double precision g2 double precision g3 double precision gc double precision hf integer k logical l0 logical l1 logical l2 logical l3 logical l4 logical l5 double precision pi double precision r double precision x pi = 3.141592653589793D+00 el = 0.5772156649015329D+00 l0 = c .eq. int ( c ) .and. c .lt. 0.0D+00 l1 = 1.0D+00 - x .lt. 1.0D-15 .and. c - a - b .le. 0.0D+00 l2 = a .eq. int ( a ) .and. a .lt. 0.0D+00 l3 = b .eq. int ( b ) .and. b .lt. 0.0D+00 l4 = c - a .eq. int ( c - a ) .and. c - a .le. 0.0D+00 l5 = c - b .eq. int ( c - b ) .and. c - b .le. 0.0D+00 if ( l0 .or. l1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'hyper_2f1_real(): Fatal error!' write ( *, '(a)' ) ' The hypergeometric series is divergent.' stop end if if ( 0.95D+00 .lt. x ) then eps = 1.0D-08 else eps = 1.0D-15 end if if ( x .eq. 0.0D+00 .or.a .eq. 0.0D+00 .or.b .eq. 0.0D+00 ) then hf = 1.0D+00 return else if ( 1.0D+00 - x .eq. eps .and. 0.0D+00 .lt. c - a - b ) then gc = gamma ( c ) gcab = gamma ( c - a - b ) gca = gamma ( c - a ) gcb = gamma ( c - b ) hf = gc * gcab / ( gca * gcb ) return else if ( 1.0D+00 + x .le. eps .and. & abs ( c - a + b - 1.0D+00 ) .le. eps ) then g0 = dsqrt ( pi ) * 2.0D+00 ** ( - a ) g1 = gamma ( c ) g2 = gamma ( 1.0D+00 + a / 2.0D+00 - b ) g3 = gamma ( 0.5D+00 + 0.5D+00 * a ) hf = g0 * g1 / ( g2 * g3 ) return else if ( l2 .or. l3 ) then if ( l2 ) then nm = int ( abs ( a ) ) end if if ( l3 ) then nm = int ( abs ( b ) ) end if hf = 1.0D+00 r = 1.0D+00 do k = 1, nm r = r * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00) & / ( k * ( c + k - 1.0D+00 ) ) * x hf = hf + r end do return else if ( l4 .or. l5 ) then if ( l4 ) then nm = int ( abs ( c - a ) ) end if if ( l5 ) then nm = int ( abs ( c - b ) ) end if hf = 1.0D+00 r = 1.0D+00 do k = 1, nm r = r * ( c - a + k - 1.0D+00 ) * ( c - b + k - 1.0D+00 ) & / ( k * ( c + k - 1.0D+00 ) ) * x hf = hf + r end do hf = ( 1.0D+00 - x ) ** ( c - a - b ) * hf return end if aa = a bb = b x1 = x if ( x .lt. 0.0D+00 ) then x = x / ( x - 1.0D+00 ) if ( a .lt. c .and. b .lt. a .and. 0.0D+00 .lt. b ) then a = bb b = aa end if b = c - b end if if ( 0.75D+00 .le. x ) then gm = 0.0D+00 if ( abs ( c - a - b - int ( c - a - b ) ) .lt. 1.0D-15 ) then m = int ( c - a - b ) ga = gamma ( a ) gb = gamma ( b ) gc = gamma ( c ) gam = gamma ( a + m ) gbm = gamma ( b + m ) call psi ( a, pa ) call psi ( b, pb ) if ( m .ne. 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 f0 = 1.0D+00 r0 = 1.0D+00 r1 = 1.0D+00 sp0 = 0.0D+00 sp = 0.0D+00 if ( 0 .le. m ) then c0 = gm * gc / ( gam * gbm ) c1 = - gc * ( x - 1.0D+00 ) ** m / ( ga * gb * rm ) do k = 1, m - 1 r0 = r0 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( k - m ) ) * ( 1.0D+00 - x ) f0 = f0 + r0 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 f1 = pa + pb + sp0 + 2.0D+00 * el + log ( 1.0D+00 - x ) do k = 1, 250 sp = sp + ( 1.0D+00 - a ) / ( k * ( a + k - 1.0D+00 ) ) & + ( 1.0 - b ) / ( k * ( b + k - 1.0 ) ) sm = 0.0D+00 do j = 1, m sm = sm + ( 1.0D+00 -a ) / ( ( j + k ) & * ( a + j + k - 1.0 ) ) + 1.0 / & ( b + j + k - 1.0D+00 ) end do rp = pa + pb + 2.0D+00 * el + sp + sm & + log ( 1.0D+00 - x ) r1 = r1 * ( a + m + k - 1.0D+00 ) * ( b + m + k - 1.0 ) & / ( k * ( m + k ) ) * ( 1.0D+00 - x ) f1 = f1 + r1 * rp if ( abs ( f1 - hw ) .lt. abs ( f1 ) * eps ) then go to 10 end if hw = f1 end do 10 continue hf = f0 * c0 + f1 * c1 else if ( m .lt. 0 ) then m = - m c0 = gm * gc / ( ga * gb * ( 1.0D+00 - x ) ** m ) c1 = - ( -1.0 ) ** m * gc / ( gam * gbm * rm ) do k = 1, m - 1 r0 = r0 * ( a - m + k - 1.0D+00 ) & * ( b - m + k - 1.0D+00 ) & / ( k * ( k - m ) ) * ( 1.0D+00 - x ) f0 = f0 + r0 end do do k = 1, m sp0 = sp0 + 1.0D+00 / k end do f1 = pa + pb - sp0 + 2.0D+00 * el + log ( 1.0D+00 - x ) do k = 1, 250 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 rp = pa+pb+2.0D+00*el+sp-sm+log(1.0D+00-x) r1 = r1*(a+k-1.0D+00)*(b+k-1.0D+00) & /(k*(m+k))*(1.0D+00-x) f1 = f1+r1*rp if (abs(f1-hw) .lt. abs(f1)*eps) then go to 20 end if hw = f1 end do 20 continue hf = f0 * c0 + f1 * c1 end if else ga = gamma ( a ) gb = gamma ( b ) gc = gamma ( c ) gca = gamma ( c - a ) gcb = gamma ( c - b ) gcab = gamma ( c - a - b ) gabc = gamma ( a + b - c ) c0 = gc * gcab / ( gca * gcb ) c1 = gc * gabc / ( ga * gb ) & * ( 1.0D+00 - x ) ** ( c - a - b ) hf = 0.0D+00 r0 = c0 r1 = c1 do k = 1, 250 r0 = r0 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & /(k*(a+b-c+k))*(1.0D+00-x) r1 = r1*(c-a+k-1.0D+00)*(c-b+k-1.0D+00)/(k*(c-a-b+k)) & *(1.0D+00-x) hf = hf + r0 + r1 if ( abs ( hf - hw ) .lt. abs ( hf ) * eps ) then go to 30 end if hw = hf end do 30 continue hf = hf + c0 + c1 end if else a0 = 1.0D+00 if ( c .gt. a .and. c .lt. 2.0D+00 * a .and. & c .gt. b .and. c .lt. 2.0D+00 * b ) then a0 = ( 1.0D+00 - x ) ** ( c - a - b ) a = c - a b = c - b end if hf = 1.0D+00 r = 1.0D+00 do k = 1, 250 r = r * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) & / ( k * ( c + k - 1.0D+00 ) ) * x hf = hf + r if ( abs ( hf - hw ) .le. abs ( hf ) * eps ) then go to 40 end if hw = hf end do 40 continue hf = a0 * hf end if if ( x1 .lt. 0.0D+00 ) then x = x1 c0 = 1.0D+00 / ( 1.0D+00 - x ) ** aa hf = c0 * hf end if a = aa b = bb if ( 120 .lt. k ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'hyper_2f1_real(): Warning!' write ( *, '(a)' ) ' Requested accuracy not achieved.' end if return end subroutine hyper_2f1_complex ( a, b, c, z, zhf ) c*********************************************************************72 c cc hyper_2f1_complex() computes the hypergeometric function 2F1(a,b,c;z) for complex argument. c c Licensing: c c This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, c they give permission to incorporate this routine into a user program c provided that the copyright is acknowledged. c c Modified: c c 03 August 2012 c c Author: c c Shanjie Zhang, Jianming Jin c c Reference: c c Shanjie Zhang, Jianming Jin, c Computation of Special Functions, c Wiley, 1996, c ISBN: 0-471-11963-6, c LC: QA351.C45. c c Input: c c double precision A, B, C, parameters. c c double complex Z, the argument. c c Output: c c Output, double complex ZHF, the value of F(a,b,c,z). c implicit none double precision a double precision a0 double precision aa double precision b double precision bb double precision c double precision ca double precision cb double precision el double precision eps double precision g0 double precision g1 double precision g2 double precision g3 double precision ga double precision gab double precision gabc double precision gam double precision gb double precision gba double precision gbm double precision gc double precision gca double precision gcab double precision gcb double precision gcbk double precision gm integer j integer k logical l0 logical l1 logical l2 logical l3 logical l4 logical l5 logical l6 integer m integer mab integer mcab integer nca integer ncb integer nm double precision pa double precision pac double precision pb double precision pca double precision pi double precision rk1 double precision rk2 double precision rm double precision sj1 double precision sj2 double precision sm double precision sp double precision sp0 double precision sq double precision t0 double precision w0 double precision ws double precision x double precision y double complex z double complex z00 double complex z1 double complex zc0 double complex zc1 double complex zf0 double complex zf1 double complex zhf double complex zp double complex zp0 double complex zr double complex zr0 double complex zr1 double complex zw x = real ( z ) y = dimag ( z ) eps = 1.0D-15 l0 = c .eq. int ( c ) .and. c .lt. 0.0D+00 l1 = abs ( 1.0D+00 - x ) .lt. eps .and. y .eq. 0.0D+00 .and. & c - a - b .le. 0.0D+00 l2 = cdabs ( z + 1.0D+00 ) .lt. eps .and. & abs ( c - a + b - 1.0D+00 ) .lt. eps l3 = a .eq. int ( a ) .and. a .lt. 0.0D+00 l4 = b .eq. int ( b ) .and. b .lt. 0.0D+00 l5 = c - a .eq. int ( c - a ) .and. c - a .le. 0.0D+00 l6 = c - b .eq. int ( c - b ) .and. c - b .le. 0.0D+00 aa = a bb = b a0 = cdabs ( z ) if ( 0.95D+00 .lt. 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)' ) 'hyper_2f1_complex(): Fatal error!' write ( *, '(a)' ) ' The hypergeometric series is divergent.' stop end if if ( a0 .eq. 0.0D+00 .or. a .eq. 0.0D+00 .or. & b .eq. 0.0D+00 ) then zhf = dcmplx ( 1.0D+00, 0.0D+00 ) else if ( z .eq. 1.0D+00.and. 0.0D+00 .lt. c - a - b ) then gc = gamma ( c ) gcab = gamma ( c - a - b ) gca = gamma ( c - a ) gcb = gamma ( c - b ) zhf = gc * gcab / ( gca * gcb ) else if ( l2 ) then g0 = dsqrt ( pi ) * 2.0D+00 ** ( - a ) g1 = gamma ( c ) g2 = gamma ( 1.0D+00 + a / 2.0D+00 - b ) g3 = gamma ( 0.5D+00 + 0.5D+00 * a ) 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 = dcmplx ( 1.0D+00, 0.0D+00 ) zr = dcmplx ( 1.0D+00, 0.0D+00 ) 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 = dcmplx ( 1.0D+00, 0.0D+00 ) zr = dcmplx ( 1.0D+00, 0.0D+00 ) 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 .le. 1.0D+00 ) then if ( x .lt. 0.0D+00 ) then z1 = z / ( z - 1.0D+00 ) if ( a .lt. c .and. b .lt. a .and. 0.0D+00 .lt. b ) then a = bb b = aa end if zc0 = 1.0D+00 / ( ( 1.0D+00 - z ) ** a ) zhf = dcmplx ( 1.0D+00, 0.0D+00 ) zr0 = dcmplx ( 1.0D+00, 0.0D+00 ) 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 ( cdabs ( zhf - zw ) .lt. cdabs ( zhf ) * eps ) then go to 10 end if zw = zhf end do 10 continue zhf = zc0 * zhf else if ( 0.90D+00 .le. a0 ) then gm = 0.0D+00 mcab = int ( c - a - b + eps * dsign ( 1.0D+00, c - a - b ) ) if ( abs ( c - a - b - mcab ) .lt. eps ) then m = int ( c - a - b ) ga = gamma ( a ) gb = gamma ( b ) gc = gamma ( c ) gam = gamma ( a + m ) gbm = gamma ( b + m ) call psi ( a, pa ) call psi ( b, pb ) if ( m .ne. 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 = dcmplx ( 1.0D+00, 0.0D+00 ) zr0 = dcmplx ( 1.0D+00, 0.0D+00 ) zr1 = dcmplx ( 1.0D+00, 0.0D+00 ) sp0 = 0.0D+00 sp = 0.0D+00 if ( 0 .le. 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 & + cdlog ( 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 & + cdlog ( 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 ( cdabs ( zf1 - zw ) .lt. cdabs ( zf1 ) * eps ) then go to 20 end if zw = zf1 end do 20 continue zhf = zf0 * zc0 + zf1 * zc1 else if ( m .lt. 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 + cdlog ( 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 & + cdlog ( 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 ( cdabs ( zf1 - zw ) .lt. cdabs ( zf1 ) * eps ) then go to 30 end if zw = zf1 end do 30 continue zhf = zf0 * zc0 + zf1 * zc1 end if else ga = gamma ( a ) gb = gamma ( b ) gc = gamma ( c ) gca = gamma ( c - a ) gcb = gamma ( c - b ) gcab = gamma ( c - a - b ) gabc = gamma ( a + b - c ) zc0 = gc * gcab / ( gca * gcb ) zc1 = gc * gabc / ( ga * gb ) & * ( 1.0D+00 - z ) ** ( c - a - b ) zhf = dcmplx ( 0.0D+00, 0.0D+00 ) 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 ( cdabs ( zhf - zw ) .lt. cdabs ( zhf ) * eps ) then go to 40 end if zw = zhf end do 40 continue zhf = zhf + zc0 + zc1 end if else z00 = dcmplx ( 1.0D+00, 0.0D+00 ) if ( c - a .lt. a .and. c - b .lt. b ) then z00 = ( 1.0D+00 - z ) ** ( c - a - b ) a = c - a b = c - b end if zhf = dcmplx ( 1.0D+00, 0.0D+00 ) zr = dcmplx ( 1.0D+00, 0.0D+00 ) 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 ( cdabs ( zhf - zw ) .le. cdabs ( zhf ) * eps ) then go to 50 end if zw = zhf end do 50 continue zhf = z00 * zhf end if else if ( 1.0D+00 .lt. a0 ) then mab = int ( a - b + eps * dsign ( 1.0D+00, a - b ) ) if ( abs ( a - b - mab ) .lt. eps .and. a0 .le. 1.1D+00 ) then b = b + eps end if if ( eps .lt. abs ( a - b - mab ) ) then ga = gamma ( a ) gb = gamma ( b ) gc = gamma ( c ) gab = gamma ( a - b ) gba = gamma ( b - a ) gca = gamma ( c - a ) gcb = gamma ( c - b ) zc0 = gc * gba / ( gca * gb * ( - z ) ** a ) zc1 = gc * gab / ( gcb * ga * ( - z ) ** b ) zr0 = zc0 zr1 = zc1 zhf = dcmplx ( 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 ( cdabs ( ( zhf - zw ) / zhf ) .le. eps ) then go to 60 end if zw = zhf end do 60 continue zhf = zhf + zc0 + zc1 else if ( a - b .lt. 0.0D+00 ) then a = bb b = aa end if ca = c - a cb = c - b nca = int ( ca + eps * dsign ( 1.0D+00, ca ) ) ncb = int ( cb + eps * dsign ( 1.0D+00, cb ) ) if ( abs ( ca - nca ) .lt. eps .or. & abs ( cb - ncb ) .lt. eps ) then c = c + eps end if ga = gamma ( a ) gc = gamma ( c ) gcb = gamma ( c - b ) call psi ( a, pa ) call psi ( c - a, pca ) call psi ( a - c, pac ) mab = int ( a - b + eps ) zc0 = gc / ( ga * ( - z ) ** b ) gm = gamma ( a - b ) zf0 = gm / gcb * zc0 zr = zc0 do k = 1, mab - 1 zr = zr * ( b + k - 1.0D+00 ) / ( k * z ) t0 = a - b - k g0 = gamma ( t0 ) gcbk = gamma ( c - b - k ) zf0 = zf0 + zr * g0 / gcbk end do if ( mab .eq. 0 ) then zf0 = dcmplx ( 0.0D+00, 0.0D+00 ) 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 + cdlog ( - 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 / dtan ( pi * ( k + a - c ) ) + cdlog ( - z ) zf1 = zf1 + rk2 * zr * zp ws = cdabs ( zf1 ) if ( abs ( ( ws - w0 ) / ws ) .lt. eps ) then go to 70 end if w0 = ws end do 70 continue zhf = zf0 + zf1 end if end if a = aa b = bb if ( 150 .lt. k ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'hyper_2f1_complex(): Warning!' write ( *, '(a)' ) & ' The solution returned may have low accuracy.' end if return end subroutine psi ( x, ps ) c*********************************************************************72 c cc psi() computes the Psi function. c c Licensing: c c This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, c they give permission to incorporate this routine into a user program c provided that the copyright is acknowledged. c c Modified: c c 22 July 2012 c c Author: c c Shanjie Zhang, Jianming Jin c c Reference: c c Shanjie Zhang, Jianming Jin, c Computation of Special Functions, c Wiley, 1996, c ISBN: 0-471-11963-6, c LC: QA351.C45. c c Input: c c double precision X, the argument. c c Output: c c double precision PS, the value of the function. c implicit none double precision a1 double precision a2 double precision a3 double precision a4 double precision a5 double precision a6 double precision a7 double precision a8 double precision el integer k integer n double precision pi double precision ps double precision s double precision x double precision x2 double precision xa xa = abs ( x ) pi = 3.141592653589793D+00 el = 0.5772156649015329D+00 s = 0.0D+00 if ( x .eq. int ( x ) .and. x .le. 0.0D+00 ) then ps = 1.0D+300 return else if ( xa .eq. int ( xa ) ) then n = int ( xa ) do k = 1, n - 1 s = s + 1.0D+00 / k end do ps = - el + s else if ( xa + 0.5D+00 .eq. int ( xa + 0.5D+00 ) ) then n = int ( xa - 0.5D+00 ) do k = 1, n s = s + 1.0D+00 / ( 2.0D+00 * k - 1.0D+00 ) end do ps = - el + 2.0D+00 * s - 1.386294361119891D+00 else if ( xa .lt. 10.0D+00 ) then n = 10 - int ( xa ) do k = 0, n - 1 s = s + 1.0D+00 / ( xa + k ) end do xa = xa + n end if x2 = 1.0D+00 / ( xa * xa ) a1 = -0.8333333333333D-01 a2 = 0.83333333333333333D-02 a3 = -0.39682539682539683D-02 a4 = 0.41666666666666667D-02 a5 = -0.75757575757575758D-02 a6 = 0.21092796092796093D-01 a7 = -0.83333333333333333D-01 a8 = 0.4432598039215686D+00 ps = log ( xa ) - 0.5D+00 / xa + x2 * ((((((( & a8 * x2 & + a7 ) * x2 & + a6 ) * x2 & + a5 ) * x2 & + a4 ) * x2 & + a3 ) * x2 & + a2 ) * x2 & + a1 ) ps = ps - s end if if ( x .lt. 0.0D+00 ) then ps = ps - pi * cos ( pi * x ) / sin ( pi * x ) - 1.0D+00 / x end if return end