00001 #include "cdflib.h"
00002 double dlnbet(double *a0,double *b0)
00003 /*
00004 **********************************************************************
00005  
00006      double dlnbet(a0,b0)
00007           Double precision LN of the complete BETa
00008  
00009  
00010                               Function
00011  
00012  
00013      Returns the natural log of the complete beta function,
00014      i.e.,
00015  
00016                   ln( Gamma(a)*Gamma(b) / Gamma(a+b)
00017  
00018  
00019                               Arguments
00020  
00021  
00022    A,B --> The (symmetric) arguments to the complete beta
00023                   DOUBLE PRECISION A, B
00024  
00025  
00026                               Method
00027  
00028  
00029      Renames BETALN from:
00030      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00031      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00032      Trans. Math.  Softw. 18 (1993), 360-373.
00033  
00034 **********************************************************************
00035 -----------------------------------------------------------------------
00036      EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
00037 -----------------------------------------------------------------------
00038      E = 0.5*LN(2*PI)
00039 --------------------------
00040 */
00041 {
00042 static double e = .918938533204673e0;
00043 static double dlnbet,a,b,c,h,u,v,w,z;
00044 static int i,n;
00045 static double T1;
00046 /*
00047      ..
00048      .. Executable Statements ..
00049 */
00050     a = fifdmin1(*a0,*b0);
00051     b = fifdmax1(*a0,*b0);
00052     if(a >= 8.0e0) goto S100;
00053     if(a >= 1.0e0) goto S20;
00054 /*
00055 -----------------------------------------------------------------------
00056                    PROCEDURE WHEN A .LT. 1
00057 -----------------------------------------------------------------------
00058 */
00059     if(b >= 8.0e0) goto S10;
00060     T1 = a+b;
00061     dlnbet = gamln(&a)+(gamln(&b)-gamln(&T1));
00062     return dlnbet;
00063 S10:
00064     dlnbet = gamln(&a)+algdiv(&a,&b);
00065     return dlnbet;
00066 S20:
00067 /*
00068 -----------------------------------------------------------------------
00069                 PROCEDURE WHEN 1 .LE. A .LT. 8
00070 -----------------------------------------------------------------------
00071 */
00072     if(a > 2.0e0) goto S40;
00073     if(b > 2.0e0) goto S30;
00074     dlnbet = gamln(&a)+gamln(&b)-gsumln(&a,&b);
00075     return dlnbet;
00076 S30:
00077     w = 0.0e0;
00078     if(b < 8.0e0) goto S60;
00079     dlnbet = gamln(&a)+algdiv(&a,&b);
00080     return dlnbet;
00081 S40:
00082 /*
00083                 REDUCTION OF A WHEN B .LE. 1000
00084 */
00085     if(b > 1000.0e0) goto S80;
00086     n = a-1.0e0;
00087     w = 1.0e0;
00088     for(i=1; i<=n; i++) {
00089         a -= 1.0e0;
00090         h = a/b;
00091         w *= (h/(1.0e0+h));
00092     }
00093     w = log(w);
00094     if(b < 8.0e0) goto S60;
00095     dlnbet = w+gamln(&a)+algdiv(&a,&b);
00096     return dlnbet;
00097 S60:
00098 /*
00099                  REDUCTION OF B WHEN B .LT. 8
00100 */
00101     n = b-1.0e0;
00102     z = 1.0e0;
00103     for(i=1; i<=n; i++) {
00104         b -= 1.0e0;
00105         z *= (b/(a+b));
00106     }
00107     dlnbet = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
00108     return dlnbet;
00109 S80:
00110 /*
00111                 REDUCTION OF A WHEN B .GT. 1000
00112 */
00113     n = a-1.0e0;
00114     w = 1.0e0;
00115     for(i=1; i<=n; i++) {
00116         a -= 1.0e0;
00117         w *= (a/(1.0e0+a/b));
00118     }
00119     dlnbet = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
00120     return dlnbet;
00121 S100:
00122 /*
00123 -----------------------------------------------------------------------
00124                    PROCEDURE WHEN A .GE. 8
00125 -----------------------------------------------------------------------
00126 */
00127     w = bcorr(&a,&b);
00128     h = a/b;
00129     c = h/(1.0e0+h);
00130     u = -((a-0.5e0)*log(c));
00131     v = b*alnrel(&h);
00132     if(u <= v) goto S110;
00133     dlnbet = -(0.5e0*log(b))+e+w-v-u;
00134     return dlnbet;
00135 S110:
00136     dlnbet = -(0.5e0*log(b))+e+w-u-v;
00137     return dlnbet;
00138 } /* END */