Doxygen Source Code Documentation
cdf_10.c File Reference
#include "cdflib.h"
Go to the source code of this file.
Functions | |
void | bratio (double *a, double *b, double *x, double *y, double *w, double *w1, int *ierr) |
Function Documentation
|
Definition at line 2 of file cdf_10.c. References a, apser(), basym(), bfrac(), bgrat(), bpser(), bup(), fifdmax1(), fifdmin1(), fpser(), ind, spmpar(), x0, and y0. Referenced by cumbet(), cumf(), and cumfnc().
00038 { 00039 static int K1 = 1; 00040 static double a0,b0,eps,lambda,t,x0,y0,z; 00041 static int ierr1,ind,n; 00042 static double T2,T3,T4,T5; 00043 /* 00044 .. 00045 .. Executable Statements .. 00046 */ 00047 /* 00048 ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST 00049 FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0 00050 */ 00051 eps = spmpar(&K1); 00052 *w = *w1 = 0.0e0; 00053 if(*a < 0.0e0 || *b < 0.0e0) goto S270; 00054 if(*a == 0.0e0 && *b == 0.0e0) goto S280; 00055 if(*x < 0.0e0 || *x > 1.0e0) goto S290; 00056 if(*y < 0.0e0 || *y > 1.0e0) goto S300; 00057 z = *x+*y-0.5e0-0.5e0; 00058 if(fabs(z) > 3.0e0*eps) goto S310; 00059 *ierr = 0; 00060 if(*x == 0.0e0) goto S210; 00061 if(*y == 0.0e0) goto S230; 00062 if(*a == 0.0e0) goto S240; 00063 if(*b == 0.0e0) goto S220; 00064 eps = fifdmax1(eps,1.e-15); 00065 if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260; 00066 ind = 0; 00067 a0 = *a; 00068 b0 = *b; 00069 x0 = *x; 00070 y0 = *y; 00071 if(fifdmin1(a0,b0) > 1.0e0) goto S40; 00072 /* 00073 PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1 00074 */ 00075 if(*x <= 0.5e0) goto S10; 00076 ind = 1; 00077 a0 = *b; 00078 b0 = *a; 00079 x0 = *y; 00080 y0 = *x; 00081 S10: 00082 if(b0 < fifdmin1(eps,eps*a0)) goto S90; 00083 if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100; 00084 if(fifdmax1(a0,b0) > 1.0e0) goto S20; 00085 if(a0 >= fifdmin1(0.2e0,b0)) goto S110; 00086 if(pow(x0,a0) <= 0.9e0) goto S110; 00087 if(x0 >= 0.3e0) goto S120; 00088 n = 20; 00089 goto S140; 00090 S20: 00091 if(b0 <= 1.0e0) goto S110; 00092 if(x0 >= 0.3e0) goto S120; 00093 if(x0 >= 0.1e0) goto S30; 00094 if(pow(x0*b0,a0) <= 0.7e0) goto S110; 00095 S30: 00096 if(b0 > 15.0e0) goto S150; 00097 n = 20; 00098 goto S140; 00099 S40: 00100 /* 00101 PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1 00102 */ 00103 if(*a > *b) goto S50; 00104 lambda = *a-(*a+*b)**x; 00105 goto S60; 00106 S50: 00107 lambda = (*a+*b)**y-*b; 00108 S60: 00109 if(lambda >= 0.0e0) goto S70; 00110 ind = 1; 00111 a0 = *b; 00112 b0 = *a; 00113 x0 = *y; 00114 y0 = *x; 00115 lambda = fabs(lambda); 00116 S70: 00117 if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110; 00118 if(b0 < 40.0e0) goto S160; 00119 if(a0 > b0) goto S80; 00120 if(a0 <= 100.0e0) goto S130; 00121 if(lambda > 0.03e0*a0) goto S130; 00122 goto S200; 00123 S80: 00124 if(b0 <= 100.0e0) goto S130; 00125 if(lambda > 0.03e0*b0) goto S130; 00126 goto S200; 00127 S90: 00128 /* 00129 EVALUATION OF THE APPROPRIATE ALGORITHM 00130 */ 00131 *w = fpser(&a0,&b0,&x0,&eps); 00132 *w1 = 0.5e0+(0.5e0-*w); 00133 goto S250; 00134 S100: 00135 *w1 = apser(&a0,&b0,&x0,&eps); 00136 *w = 0.5e0+(0.5e0-*w1); 00137 goto S250; 00138 S110: 00139 *w = bpser(&a0,&b0,&x0,&eps); 00140 *w1 = 0.5e0+(0.5e0-*w); 00141 goto S250; 00142 S120: 00143 *w1 = bpser(&b0,&a0,&y0,&eps); 00144 *w = 0.5e0+(0.5e0-*w1); 00145 goto S250; 00146 S130: 00147 T2 = 15.0e0*eps; 00148 *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2); 00149 *w1 = 0.5e0+(0.5e0-*w); 00150 goto S250; 00151 S140: 00152 *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps); 00153 b0 += (double)n; 00154 S150: 00155 T3 = 15.0e0*eps; 00156 bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1); 00157 *w = 0.5e0+(0.5e0-*w1); 00158 goto S250; 00159 S160: 00160 n = b0; 00161 b0 -= (double)n; 00162 if(b0 != 0.0e0) goto S170; 00163 n -= 1; 00164 b0 = 1.0e0; 00165 S170: 00166 *w = bup(&b0,&a0,&y0,&x0,&n,&eps); 00167 if(x0 > 0.7e0) goto S180; 00168 *w += bpser(&a0,&b0,&x0,&eps); 00169 *w1 = 0.5e0+(0.5e0-*w); 00170 goto S250; 00171 S180: 00172 if(a0 > 15.0e0) goto S190; 00173 n = 20; 00174 *w += bup(&a0,&b0,&x0,&y0,&n,&eps); 00175 a0 += (double)n; 00176 S190: 00177 T4 = 15.0e0*eps; 00178 bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1); 00179 *w1 = 0.5e0+(0.5e0-*w); 00180 goto S250; 00181 S200: 00182 T5 = 100.0e0*eps; 00183 *w = basym(&a0,&b0,&lambda,&T5); 00184 *w1 = 0.5e0+(0.5e0-*w); 00185 goto S250; 00186 S210: 00187 /* 00188 TERMINATION OF THE PROCEDURE 00189 */ 00190 if(*a == 0.0e0) goto S320; 00191 S220: 00192 *w = 0.0e0; 00193 *w1 = 1.0e0; 00194 return; 00195 S230: 00196 if(*b == 0.0e0) goto S330; 00197 S240: 00198 *w = 1.0e0; 00199 *w1 = 0.0e0; 00200 return; 00201 S250: 00202 if(ind == 0) return; 00203 t = *w; 00204 *w = *w1; 00205 *w1 = t; 00206 return; 00207 S260: 00208 /* 00209 PROCEDURE FOR A AND B .LT. 1.E-3*EPS 00210 */ 00211 *w = *b/(*a+*b); 00212 *w1 = *a/(*a+*b); 00213 return; 00214 S270: 00215 /* 00216 ERROR RETURN 00217 */ 00218 *ierr = 1; 00219 return; 00220 S280: 00221 *ierr = 2; 00222 return; 00223 S290: 00224 *ierr = 3; 00225 return; 00226 S300: 00227 *ierr = 4; 00228 return; 00229 S310: 00230 *ierr = 5; 00231 return; 00232 S320: 00233 *ierr = 6; 00234 return; 00235 S330: 00236 *ierr = 7; 00237 return; 00238 } /* END */ |