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