{*****************************************************************
*      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" *
*                                                                *
*                              TPW Release By J-P Moreau, Paris. *
*                                      (www.jpmoreau.fr)         *
*****************************************************************}
Program MBeta;

Uses WinCrt1;

Var
    BT, P, Q: double;

Function GAMMA(X: double): double; Forward;

Function BETA(P, Q: double): double;

{       ==========================================
!       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)
!       ========================================== }
Var     GP,GQ,GPQ,PPQ: double;
Begin
        GP := GAMMA(P);
        GQ := GAMMA(Q);
        PPQ:=P+Q;
        GPQ := GAMMA(PPQ);
        BETA := (GP*GQ/GPQ)
End;


Function GAMMA(X: double): double;

{       ==================================================
!       Purpose: Compute gamma function gamma(x)
!       Input :  x  --- Argument of gamma(x)
!                       ( x is not equal to 0,-1,-2)
!       Output:  gamma(x)
!       ================================================== }
Var
        GA, GR, R, Z: double;
        G: Array[1..26] of double;
        K, M, M1: integer;
Begin
        if X=int(X) then
	   if (X>0.0) then
           begin
              GA:=1.0;
              M1:=Round(X-1);
              for K:=2 to M1 do  GA:=GA*K
           end
           else
              GA:=1E+100
	 else
         begin
	    if abs(X)>1.0 then
            begin
              Z:=abs(X);
              M:=Round(Z);
              R:=1.0;
              for K:=1 to M do  R:=R*(Z-K);
              Z:=Z-M
            end
            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 Downto 1 do  GR:=GR*Z+G[K];
           GA:=1.0/(GR*Z);
           if abs(X)>1.0 then
           begin
              GA:=GA*R;
              if X < 0.0 then GA:=-PI/(X*GA*sin(PI*X))
           end
        end;
        GAMMA := GA
End;


BEGIN
      writeln;
      writeln('    p       q                B(p,q)        ');
      writeln('  -----------------------------------------');
      P:=1.5;
      Q:=2.0;
      BT:=BETA(P,Q);
      writeln('  ',P:5:1,'   ',Q:5:1,'     ',BT);
      P:=2.5;
      Q:=2.0;
      BT:=BETA(P,Q);
      writeln('  ',P:5:1,'   ',Q:5:1,'     ',BT);
      P:=1.5;
      Q:=3.0;
      BT:=BETA(P,Q);
      writeln('  ',P:5:1,'   ',Q:5:1,'     ',BT);
      Readkey;
      DonewinCrt;
END.

{ end of file mbeta.pas }