/****************************************************************** * 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" * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *****************************************************************/ #include #include double GAMMA(double); double BT, P, Q; double BETA(double P, double 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 for computing â(x) ! ========================================== */ double GP,GQ,GPQ,PPQ; GP = GAMMA(P); GQ = GAMMA(Q); PPQ=P+Q; GPQ = GAMMA(PPQ); return (GP*GQ/GPQ); } double GAMMA(double X) { /* ================================================== ! Purpose: Compute gamma function gamma(x) ! Input : x --- Argument of gamma(x) ! ( x is not equal to 0,-1,-2) ! Output: gamma(x) ! ================================================== */ double GA, GR,PI, R, Z; double G[27]; int K, M, M1; PI=3.141592653589793; if (X==int(X)) if (X>0.0) { GA=1.0; M1=X-1; for (K=2; K<=M1; K++) GA=GA*K; } else GA=1e+100; else { if (fabs(X)>1.0) { Z=fabs(X); M=int(Z); R=1.0; for (K=1; K<=M; K++) R=R*(Z-K); Z=Z-M; } else Z=X; G[1]=1.0; G[2]=0.5772156649015329; G[3]=-0.6558780715202538; G[4]=-0.420026350340952e-1; G[5]=0.1665386113822915; G[6]=-0.421977345555443e-1; G[7]=-0.96219715278770e-2; G[8]=0.72189432466630e-2; G[9]=-0.11651675918591e-2; G[10]=-0.2152416741149e-3; G[11]=0.1280502823882e-3; G[12]=-0.201348547807e-4; G[13]=-0.12504934821e-5; G[14]=0.11330272320e-5; G[15]=-0.2056338417e-6; G[16]=0.61160950e-8; G[17]=0.50020075e-8; G[18]=-0.11812746e-8; G[19]=0.1043427e-9; G[20]=0.77823e-11; G[21]=-0.36968e-11; G[22]=0.51e-12; G[23]=-0.206e-13; G[24]=-0.54e-14; G[25]=0.14e-14; G[26]=0.1e-15; GR=G[26]; for (K=25; K>0; K--) GR=GR*Z+G[K]; GA=1.0/(GR*Z); if (fabs(X)>1.0) { GA=GA*R; if (X < 0.0) GA=-PI/(X*GA*sin(PI*X)); } } return GA; } void main() { printf("\n p q B(p,q)\n"); printf(" ---------------------------------\n"); P=1.5; Q=2.0; BT=BETA(P,Q); printf(" %5.1f %5.1f %e\n", P, Q, BT); P=2.5; Q=2.0; BT=BETA(P,Q); printf(" %5.1f %5.1f %e\n", P, Q, BT); P=1.5; Q=3.0; BT=BETA(P,Q); printf(" %5.1f %5.1f %e\n\n", P, Q, BT); } // end of file mbeta.cpp