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