betaf Subroutine

subroutine betaf ( p, q, bt )

  !*****************************************************************************80
  !
  !! BETA computes the Beta function B(p,q).
  !
  !  Licensing:
  !
  !    The original FORTRAN77 version of this routine is copyrighted by 
  !    Shanjie Zhang and Jianming Jin.  However, they give permission to 
  !    incorporate this routine into a user program that the copyright 
  !    is acknowledged.
  !
  !  Modified:
  !
  !    12 March 2012
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
  !    FORTRAN90 version by John Burkardt.
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) P, Q, the parameters.
  !    0 < P, 0 < Q.
  !
  !    Output, real ( kind = 8 ) BT, the value of B(P,Q).
  !
  implicit none

  real ( kind = 8 ) bt
  real ( kind = 8 ) gp
  real ( kind = 8 ) gpq
  real ( kind = 8 ) gq
  real ( kind = 8 ) p
  real ( kind = 8 ) ppq
  real ( kind = 8 ) q

  call gammaf ( p, gp )
  call gammaf ( q, gq )
  ppq = p + q
  call gammaf ( ppq, gpq )
  bt = gp * gq / gpq

  return
end subroutine betaf