chgubi Subroutine

subroutine chgubi ( a, b, x, hu, id )

  !*****************************************************************************80
  !
  !! CHGUBI: confluent hypergeometric function with integer argument B.
  !
  !  Discussion:
  !
  !    This procedure computes the confluent hypergeometric function
  !    U(a,b,x) with integer ( kind = 4 ) b ( b = ±1,±2,... )
  !
  !  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:
  !
  !    31 July 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, parameters.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) HU, the value of U(a,b,x).
  !
  !    Output, integer ( kind = 4 ) ID, the estimated number of significant
  !    digits.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) a1
  real ( kind = 8 ) a2
  real ( kind = 8 ) b
  real ( kind = 8 ) da1
  real ( kind = 8 ) da2
  real ( kind = 8 ) db1
  real ( kind = 8 ) db2
  real ( kind = 8 ) el
  real ( kind = 8 ) ga
  real ( kind = 8 ) ga1
  real ( kind = 8 ) h0
  real ( kind = 8 ) hm1
  real ( kind = 8 ) hm2
  real ( kind = 8 ) hm3
  real ( kind = 8 ) hmax
  real ( kind = 8 ) hmin
  real ( kind = 8 ) hu
  real ( kind = 8 ) hu1
  real ( kind = 8 ) hu2
  real ( kind = 8 ) hw
  integer ( kind = 4 ) id
  integer ( kind = 4 ) id1
  integer ( kind = 4 ) id2
  integer ( kind = 4 ) j 
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  real ( kind = 8 ) ps
  real ( kind = 8 ) r
  real ( kind = 8 ) rn
  real ( kind = 8 ) rn1
  real ( kind = 8 ) s0
  real ( kind = 8 ) s1
  real ( kind = 8 ) s2
  real ( kind = 8 ) sa
  real ( kind = 8 ) sb
  real ( kind = 8 ) ua
  real ( kind = 8 ) ub
  real ( kind = 8 ) x

  id = -100
  el = 0.5772156649015329D+00
  n = abs ( b - 1 )
  rn1 = 1.0D+00
  rn = 1.0D+00
  do j = 1, n
     rn = rn * j
     if ( j == n - 1 ) then
        rn1 = rn
     end if
  end do

  call psi ( a, ps )
  call gammaf ( a, ga )

  if ( 0.0D+00 < b ) then
     a0 = a
     a1 = a - n
     a2 = a1
     call gammaf ( a1, ga1 )
     ua = ( - 1 ) ** ( n - 1 ) / ( rn * ga1 )
     ub = rn1 / ga * x ** ( - n )
  else
     a0 = a + n
     a1 = a0
     a2 = a
     call gammaf ( a1, ga1 )
     ua = ( - 1 ) ** ( n - 1 ) / ( rn * ga ) * x ** n
     ub = rn1 / ga1
  end if

  hm1 = 1.0D+00
  r = 1.0D+00
  hmax = 0.0D+00
  hmin = 1.0D+300

  do k = 1, 150
     r = r * ( a0 + k - 1.0D+00 ) * x / ( ( n + k ) * k )
     hm1 = hm1 + r
     hu1 = abs ( hm1 )
     hmax = max ( hmax, hu1 )
     hmin = min ( hmin, hu1 )
     if ( abs ( hm1 - h0 ) < abs ( hm1 ) * 1.0D-15 ) then
        exit
     end if
     h0 = hm1
  end do

  da1 = log10 ( hmax )
  if ( hmin /= 0.0D+00 ) then
     da2 = log10 ( hmin )
  end if
  id = 15 - abs ( da1 - da2 )
  hm1 = hm1 * log ( x )
  s0 = 0.0D+00
  do m = 1, n
     if ( 0.0D+00 <= b ) then
        s0 = s0 - 1.0D+00 / m
     else
        s0 = s0 + ( 1.0D+00 - a ) / ( m * ( a + m - 1.0D+00 ) )
     end if
  end do
  hm2 = ps + 2.0D+00 * el + s0
  r = 1.0D+00
  hmax = 0.0D+00
  hmin = 1.0D+300
  do k = 1, 150
     s1 = 0.0D+00
     s2 = 0.0D+00
     if ( 0.0D+00 < b ) then
        do m = 1, k
           s1 = s1 - ( m + 2.0D+00 * a - 2.0D+00 ) / ( m * ( m + a - 1.0D+00 ) )
        end do
        do m = 1, n
           s2 = s2 + 1.0D+00 / ( k + m )
        end do
     else
        do m = 1, k + n
           s1 = s1 + ( 1.0D+00 - a ) / ( m * ( m + a - 1.0D+00 ) )
        end do
        do m = 1, k
           s2 = s2 + 1.0D+00 / m
        end do
     end if
     hw = 2.0D+00 * el + ps + s1 - s2
     r = r * ( a0 + k - 1.0D+00 ) * x / ( ( n + k ) * k )
     hm2 = hm2 + r * hw
     hu2 = abs ( hm2 )
     hmax = max ( hmax, hu2 )
     hmin = min ( hmin, hu2 )

     if ( abs ( ( hm2 - h0 ) / hm2 ) < 1.0D-15 ) then
        exit
     end if

     h0 = hm2

  end do

  db1 = log10 ( hmax )
  if ( hmin /= 0.0D+00 ) then
     db2 = log10 ( hmin )
  end if
  id1 = 15 - abs ( db1 - db2 )
  id = min ( id, id1 )

  if ( n == 0 ) then
     hm3 = 0.0D+00
  else
     hm3 = 1.0D+00
  end if

  r = 1.0D+00
  do k = 1, n - 1
     r = r * ( a2 + k - 1.0D+00 ) / ( ( k - n ) * k ) * x
     hm3 = hm3 + r
  end do

  sa = ua * ( hm1 + hm2 )
  sb = ub * hm3
  hu = sa + sb

  if ( sa /= 0.0D+00 ) then
     id1 = int ( log10 ( abs ( sa ) ) )
  end if

  if ( hu /= 0.0D+00 ) then
     id2 = int ( log10 ( abs ( hu ) ) )
  end if

  if ( sa * sb < 0.0D+00 ) then
     id = id - abs ( id1 - id2 )
  end if

  return
end subroutine chgubi