Doxygen Source Code Documentation
cdf_11.c File Reference
#include "cdflib.h"
Go to the source code of this file.
Functions | |
double | brcmp1 (int *mu, double *a, double *b, double *x, double *y) |
Function Documentation
|
Definition at line 2 of file cdf_11.c. References a, algdiv(), alnrel(), bcorr(), betaln(), brcmp1(), c, esum(), fifdmax1(), fifdmin1(), gam1(), gamln1(), i, rlog1(), v, x0, and y0. Referenced by brcmp1(), and bup().
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 */ |