incob Subroutine

subroutine incob ( a, b, x, bix )

  !*****************************************************************************80
  !
  !! INCOB computes the incomplete beta function Ix(a,b).
  !
  !  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:
  !
  !    22 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 ) BIX, the function value.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) bix
  real ( kind = 8 ) bt
  real ( kind = 8 ) dk(51)
  real ( kind = 8 ) fk(51)
  integer ( kind = 4 ) k
  real ( kind = 8 ) s0
  real ( kind = 8 ) t1
  real ( kind = 8 ) t2
  real ( kind = 8 ) ta
  real ( kind = 8 ) tb
  real ( kind = 8 ) x

  s0 = ( a + 1.0D+00 ) / ( a + b + 2.0D+00 )
  call betaf ( a, b, bt )

  if ( x <= s0 ) then

     do k = 1, 20
        dk(2*k) = k * ( b - k ) * x / &
             ( a + 2.0D+00 * k - 1.0D+00 ) / ( a + 2.0D+00 * k )
     end do

     do k = 0, 20
        dk(2*k+1) = - ( a + k ) * ( a + b + k ) * x &
             / ( a + 2.0D+00 * k ) / ( a + 2.0D+00 * k + 1.0D+00 )
     end do

     t1 = 0.0D+00
     do k = 20, 1, -1
        t1 = dk(k) / ( 1.0D+00 + t1 )
     end do
     ta = 1.0D+00 / ( 1.0D+00 + t1 )
     bix = x ** a * ( 1.0D+00 - x ) ** b / ( a * bt ) * ta

  else

     do k = 1, 20
        fk(2*k) = k * ( a - k ) * ( 1.0D+00 - x ) &
             / ( b + 2.0D+00 * k - 1.0D+00 ) / ( b + 2.0D+00 * k )
     end do

     do k = 0,20
        fk(2*k+1) = - ( b + k ) * ( a + b + k ) * ( 1.0D+00 - x ) &
             / ( b + 2.0D+00 * k ) / ( b + 2.0D+00 * k + 1.0D+00 )
     end do

     t2 = 0.0D+00
     do k = 20, 1, -1
        t2 = fk(k) / ( 1.0D+00 + t2 )
     end do
     tb = 1.0D+00 / ( 1.0D+00 + t2 )
     bix = 1.0D+00 - x ** a * ( 1.0D+00 - x ) ** b / ( b * bt ) * tb

  end if

  return
end subroutine incob