DECLARE FUNCTION BETA# (P#, Q#)
DECLARE FUNCTION GAMMA# (X#)
'******************************************************************
'*      Purpose: This program computes the beta function          *
'*               B(p,q) for p > 0 and q > 0 using                 *
'*               subroutine BETA                                  *
'*      Input :  p  --- Parameter  ( p > 0 )                      *
'*               q  --- Parameter  ( q > 0 )                      *
'*      Output:  BT --- B(p,q)                                    *
'*      Examples:                                                 *
'*                p       q           B(p,q)                      *
'*              ---------------------------------                 *
'*               1.5     2.0     .2666666667D+00                  *
'*               2.5     2.0     .1142857143D+00                  *
'*               1.5     3.0     .1523809524D+00                  *
'* -------------------------------------------------------------- *
'* REFERENCE: "Fortran Routines for Computation of Special        *
'*             Functions jin.ece.uiuc.edu/routines/routines.html" *
'*                                                                *
'*                            Basic Release By J-P Moreau, Paris. *
'*                                     (www.jpmoreau.fr)          *
'******************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

      PRINT
      PRINT "    p       q                B(p,q)        "
      PRINT "  -----------------------------------------"
      P = 1.5
      Q = 2#
      BT = BETA(P, Q)
      PRINT "  "; P; "   "; Q; "     "; BT
      P = 2.5
      Q = 2#
      BT = BETA(P, Q)
      PRINT "  "; P; "   "; Q; "     "; BT
      P = 1.5
      Q = 3#
      BT = BETA(P, Q)
      PRINT "  "; P; "   "; Q; "     "; BT
      PRINT
END

' end of file mbeta.bas

FUNCTION BETA (P, Q)

'   ==========================================
'    Purpose: Compute the beta function B(p,q)
'       Input :  p  --- Parameter  ( p > 0 )
'                q  --- Parameter  ( q > 0 )
'       Output:  B(p,q)
'    Routine called: GAMMA(x)
'   ==========================================
	GP = GAMMA(P)
	GQ = GAMMA(Q)
	PPQ = P + Q
	GPQ = GAMMA(PPQ)
	BETA = (GP * GQ / GPQ)
END FUNCTION

'*******************************************
'*           FUNCTION  GAMMA(X)            *
'* --------------------------------------- *
'* Returns the value of Gamma(x) in double *
'* precision as EXP(LN(GAMMA(X))) for X>0. *
'* (input: x  output: y)                   *
'*******************************************
Function Gamma(x)
  dim cof(6)
  cof(1)=76.18009173#   : cof(2)=-86.50532033#
  cof(3)=24.01409822#   : cof(4)=-1.231739516#
  cof(5)=0.00120858003# : cof(6)=-0.00000536382#
  stp=2.50662827465#
  half=0.5# : one=1# : fpf=5.5#
  xx=x-one
  tmp=xx+fpf
  tmp=(xx+half)*LOG(tmp)-tmp
  ser=one
  for j=1 to 6
    xx=xx+one
    ser=ser+cof(j)/xx
  next j
  Gamma = EXP(tmp+LOG(stp*ser))
End Function