00001 #include "cdflib.h"
00002 double brcmp1(int *mu,double *a,double *b,double *x,double *y)
00003 /*
00004 -----------------------------------------------------------------------
00005           EVALUATION OF  EXP(MU) * (X**A*Y**B/BETA(A,B))
00006 -----------------------------------------------------------------------
00007 */
00008 {
00009 static double Const = .398942280401433e0;
00010 static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
00011 static int i,n;
00012 /*
00013 -----------------
00014      CONST = 1/SQRT(2*PI)
00015 -----------------
00016 */
00017 static double T1,T2,T3,T4;
00018 /*
00019      ..
00020      .. Executable Statements ..
00021 */
00022     a0 = fifdmin1(*a,*b);
00023     if(a0 >= 8.0e0) goto S130;
00024     if(*x > 0.375e0) goto S10;
00025     lnx = log(*x);
00026     T1 = -*x;
00027     lny = alnrel(&T1);
00028     goto S30;
00029 S10:
00030     if(*y > 0.375e0) goto S20;
00031     T2 = -*y;
00032     lnx = alnrel(&T2);
00033     lny = log(*y);
00034     goto S30;
00035 S20:
00036     lnx = log(*x);
00037     lny = log(*y);
00038 S30:
00039     z = *a*lnx+*b*lny;
00040     if(a0 < 1.0e0) goto S40;
00041     z -= betaln(a,b);
00042     brcmp1 = esum(mu,&z);
00043     return brcmp1;
00044 S40:
00045 /*
00046 -----------------------------------------------------------------------
00047               PROCEDURE FOR A .LT. 1 OR B .LT. 1
00048 -----------------------------------------------------------------------
00049 */
00050     b0 = fifdmax1(*a,*b);
00051     if(b0 >= 8.0e0) goto S120;
00052     if(b0 > 1.0e0) goto S70;
00053 /*
00054                    ALGORITHM FOR B0 .LE. 1
00055 */
00056     brcmp1 = esum(mu,&z);
00057     if(brcmp1 == 0.0e0) return brcmp1;
00058     apb = *a+*b;
00059     if(apb > 1.0e0) goto S50;
00060     z = 1.0e0+gam1(&apb);
00061     goto S60;
00062 S50:
00063     u = *a+*b-1.e0;
00064     z = (1.0e0+gam1(&u))/apb;
00065 S60:
00066     c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
00067     brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
00068     return brcmp1;
00069 S70:
00070 /*
00071                 ALGORITHM FOR 1 .LT. B0 .LT. 8
00072 */
00073     u = gamln1(&a0);
00074     n = b0-1.0e0;
00075     if(n < 1) goto S90;
00076     c = 1.0e0;
00077     for(i=1; i<=n; i++) {
00078         b0 -= 1.0e0;
00079         c *= (b0/(a0+b0));
00080     }
00081     u = log(c)+u;
00082 S90:
00083     z -= u;
00084     b0 -= 1.0e0;
00085     apb = a0+b0;
00086     if(apb > 1.0e0) goto S100;
00087     t = 1.0e0+gam1(&apb);
00088     goto S110;
00089 S100:
00090     u = a0+b0-1.e0;
00091     t = (1.0e0+gam1(&u))/apb;
00092 S110:
00093     brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
00094     return brcmp1;
00095 S120:
00096 /*
00097                    ALGORITHM FOR B0 .GE. 8
00098 */
00099     u = gamln1(&a0)+algdiv(&a0,&b0);
00100     T3 = z-u;
00101     brcmp1 = a0*esum(mu,&T3);
00102     return brcmp1;
00103 S130:
00104 /*
00105 -----------------------------------------------------------------------
00106               PROCEDURE FOR A .GE. 8 AND B .GE. 8
00107 -----------------------------------------------------------------------
00108 */
00109     if(*a > *b) goto S140;
00110     h = *a/ *b;
00111     x0 = h/(1.0e0+h);
00112     y0 = 1.0e0/(1.0e0+h);
00113     lambda = *a-(*a+*b)**x;
00114     goto S150;
00115 S140:
00116     h = *b/ *a;
00117     x0 = 1.0e0/(1.0e0+h);
00118     y0 = h/(1.0e0+h);
00119     lambda = (*a+*b)**y-*b;
00120 S150:
00121     e = -(lambda/ *a);
00122     if(fabs(e) > 0.6e0) goto S160;
00123     u = rlog1(&e);
00124     goto S170;
00125 S160:
00126     u = e-log(*x/x0);
00127 S170:
00128     e = lambda/ *b;
00129     if(fabs(e) > 0.6e0) goto S180;
00130     v = rlog1(&e);
00131     goto S190;
00132 S180:
00133     v = e-log(*y/y0);
00134 S190:
00135     T4 = -(*a*u+*b*v);
00136     z = esum(mu,&T4);
00137     brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
00138     return brcmp1;
00139 } /* END */