Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_11.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 double brcmp1(int *mu,double *a,double *b,double *x,double *y)
00003
00004
00005
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
00015
00016
00017 static double T1,T2,T3,T4;
00018
00019
00020
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
00048
00049
00050 b0 = fifdmax1(*a,*b);
00051 if(b0 >= 8.0e0) goto S120;
00052 if(b0 > 1.0e0) goto S70;
00053
00054
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
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
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
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 }