Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_59.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 void gaminv(double *a,double *x,double *x0,double *p,double *q,
00003 int *ierr)
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 {
00058 static double a0 = 3.31125922108741e0;
00059 static double a1 = 11.6616720288968e0;
00060 static double a2 = 4.28342155967104e0;
00061 static double a3 = .213623493715853e0;
00062 static double b1 = 6.61053765625462e0;
00063 static double b2 = 6.40691597760039e0;
00064 static double b3 = 1.27364489782223e0;
00065 static double b4 = .036117081018842e0;
00066 static double c = .577215664901533e0;
00067 static double ln10 = 2.302585e0;
00068 static double tol = 1.e-5;
00069 static double amin[2] = {
00070 500.0e0,100.0e0
00071 };
00072 static double bmin[2] = {
00073 1.e-28,1.e-13
00074 };
00075 static double dmin[2] = {
00076 1.e-06,1.e-04
00077 };
00078 static double emin[2] = {
00079 2.e-03,6.e-03
00080 };
00081 static double eps0[2] = {
00082 1.e-10,1.e-08
00083 };
00084 static int K1 = 1;
00085 static int K2 = 2;
00086 static int K3 = 3;
00087 static int K8 = 0;
00088 static double am1,amax,ap1,ap2,ap3,apn,b,c1,c2,c3,c4,c5,d,e,e2,eps,g,h,pn,qg,qn,
00089 r,rta,s,s2,sum,t,u,w,xmax,xmin,xn,y,z;
00090 static int iop;
00091 static double T4,T5,T6,T7,T9;
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 e = spmpar(&K1);
00103 xmin = spmpar(&K2);
00104 xmax = spmpar(&K3);
00105 *x = 0.0e0;
00106 if(*a <= 0.0e0) goto S300;
00107 t = *p+*q-1.e0;
00108 if(fabs(t) > e) goto S320;
00109 *ierr = 0;
00110 if(*p == 0.0e0) return;
00111 if(*q == 0.0e0) goto S270;
00112 if(*a == 1.0e0) goto S280;
00113 e2 = 2.0e0*e;
00114 amax = 0.4e-10/(e*e);
00115 iop = 1;
00116 if(e > 1.e-10) iop = 2;
00117 eps = eps0[iop-1];
00118 xn = *x0;
00119 if(*x0 > 0.0e0) goto S160;
00120
00121
00122
00123
00124 if(*a > 1.0e0) goto S80;
00125 T4 = *a+1.0e0;
00126 g = Xgamm(&T4);
00127 qg = *q*g;
00128 if(qg == 0.0e0) goto S360;
00129 b = qg/ *a;
00130 if(qg > 0.6e0**a) goto S40;
00131 if(*a >= 0.30e0 || b < 0.35e0) goto S10;
00132 t = exp(-(b+c));
00133 u = t*exp(t);
00134 xn = t*exp(u);
00135 goto S160;
00136 S10:
00137 if(b >= 0.45e0) goto S40;
00138 if(b == 0.0e0) goto S360;
00139 y = -log(b);
00140 s = 0.5e0+(0.5e0-*a);
00141 z = log(y);
00142 t = y-s*z;
00143 if(b < 0.15e0) goto S20;
00144 xn = y-s*log(t)-log(1.0e0+s/(t+1.0e0));
00145 goto S220;
00146 S20:
00147 if(b <= 0.01e0) goto S30;
00148 u = ((t+2.0e0*(3.0e0-*a))*t+(2.0e0-*a)*(3.0e0-*a))/((t+(5.0e0-*a))*t+2.0e0);
00149 xn = y-s*log(t)-log(u);
00150 goto S220;
00151 S30:
00152 c1 = -(s*z);
00153 c2 = -(s*(1.0e0+c1));
00154 c3 = s*((0.5e0*c1+(2.0e0-*a))*c1+(2.5e0-1.5e0**a));
00155 c4 = -(s*(((c1/3.0e0+(2.5e0-1.5e0**a))*c1+((*a-6.0e0)**a+7.0e0))*c1+(
00156 (11.0e0**a-46.0)**a+47.0e0)/6.0e0));
00157 c5 = -(s*((((-(c1/4.0e0)+(11.0e0**a-17.0e0)/6.0e0)*c1+((-(3.0e0**a)+13.0e0)*
00158 *a-13.0e0))*c1+0.5e0*(((2.0e0**a-25.0e0)**a+72.0e0)**a-61.0e0))*c1+((
00159 (25.0e0**a-195.0e0)**a+477.0e0)**a-379.0e0)/12.0e0));
00160 xn = (((c5/y+c4)/y+c3)/y+c2)/y+c1+y;
00161 if(*a > 1.0e0) goto S220;
00162 if(b > bmin[iop-1]) goto S220;
00163 *x = xn;
00164 return;
00165 S40:
00166 if(b**q > 1.e-8) goto S50;
00167 xn = exp(-(*q/ *a+c));
00168 goto S70;
00169 S50:
00170 if(*p <= 0.9e0) goto S60;
00171 T5 = -*q;
00172 xn = exp((alnrel(&T5)+gamln1(a))/ *a);
00173 goto S70;
00174 S60:
00175 xn = exp(log(*p*g)/ *a);
00176 S70:
00177 if(xn == 0.0e0) goto S310;
00178 t = 0.5e0+(0.5e0-xn/(*a+1.0e0));
00179 xn /= t;
00180 goto S160;
00181 S80:
00182
00183
00184
00185
00186 if(*q <= 0.5e0) goto S90;
00187 w = log(*p);
00188 goto S100;
00189 S90:
00190 w = log(*q);
00191 S100:
00192 t = sqrt(-(2.0e0*w));
00193 s = t-(((a3*t+a2)*t+a1)*t+a0)/((((b4*t+b3)*t+b2)*t+b1)*t+1.0e0);
00194 if(*q > 0.5e0) s = -s;
00195 rta = sqrt(*a);
00196 s2 = s*s;
00197 xn = *a+s*rta+(s2-1.0e0)/3.0e0+s*(s2-7.0e0)/(36.0e0*rta)-((3.0e0*s2+7.0e0)*
00198 s2-16.0e0)/(810.0e0**a)+s*((9.0e0*s2+256.0e0)*s2-433.0e0)/(38880.0e0**a*
00199 rta);
00200 xn = fifdmax1(xn,0.0e0);
00201 if(*a < amin[iop-1]) goto S110;
00202 *x = xn;
00203 d = 0.5e0+(0.5e0-*x/ *a);
00204 if(fabs(d) <= dmin[iop-1]) return;
00205 S110:
00206 if(*p <= 0.5e0) goto S130;
00207 if(xn < 3.0e0**a) goto S220;
00208 y = -(w+gamln(a));
00209 d = fifdmax1(2.0e0,*a*(*a-1.0e0));
00210 if(y < ln10*d) goto S120;
00211 s = 1.0e0-*a;
00212 z = log(y);
00213 goto S30;
00214 S120:
00215 t = *a-1.0e0;
00216 T6 = -(t/(xn+1.0e0));
00217 xn = y+t*log(xn)-alnrel(&T6);
00218 T7 = -(t/(xn+1.0e0));
00219 xn = y+t*log(xn)-alnrel(&T7);
00220 goto S220;
00221 S130:
00222 ap1 = *a+1.0e0;
00223 if(xn > 0.70e0*ap1) goto S170;
00224 w += gamln(&ap1);
00225 if(xn > 0.15e0*ap1) goto S140;
00226 ap2 = *a+2.0e0;
00227 ap3 = *a+3.0e0;
00228 *x = exp((w+*x)/ *a);
00229 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
00230 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
00231 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2*(1.0e0+*x/ap3))))/ *a);
00232 xn = *x;
00233 if(xn > 1.e-2*ap1) goto S140;
00234 if(xn <= emin[iop-1]*ap1) return;
00235 goto S170;
00236 S140:
00237 apn = ap1;
00238 t = xn/apn;
00239 sum = 1.0e0+t;
00240 S150:
00241 apn += 1.0e0;
00242 t *= (xn/apn);
00243 sum += t;
00244 if(t > 1.e-4) goto S150;
00245 t = w-log(sum);
00246 xn = exp((xn+t)/ *a);
00247 xn *= (1.0e0-(*a*log(xn)-xn-t)/(*a-xn));
00248 goto S170;
00249 S160:
00250
00251
00252
00253 if(*p > 0.5e0) goto S220;
00254 S170:
00255 if(*p <= 1.e10*xmin) goto S350;
00256 am1 = *a-0.5e0-0.5e0;
00257 S180:
00258 if(*a <= amax) goto S190;
00259 d = 0.5e0+(0.5e0-xn/ *a);
00260 if(fabs(d) <= e2) goto S350;
00261 S190:
00262 if(*ierr >= 20) goto S330;
00263 *ierr += 1;
00264 gratio(a,&xn,&pn,&qn,&K8);
00265 if(pn == 0.0e0 || qn == 0.0e0) goto S350;
00266 r = rcomp(a,&xn);
00267 if(r == 0.0e0) goto S350;
00268 t = (pn-*p)/r;
00269 w = 0.5e0*(am1-xn);
00270 if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S200;
00271 *x = xn*(1.0e0-t);
00272 if(*x <= 0.0e0) goto S340;
00273 d = fabs(t);
00274 goto S210;
00275 S200:
00276 h = t*(1.0e0+w*t);
00277 *x = xn*(1.0e0-h);
00278 if(*x <= 0.0e0) goto S340;
00279 if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
00280 d = fabs(h);
00281 S210:
00282 xn = *x;
00283 if(d > tol) goto S180;
00284 if(d <= eps) return;
00285 if(fabs(*p-pn) <= tol**p) return;
00286 goto S180;
00287 S220:
00288
00289
00290
00291 if(*q <= 1.e10*xmin) goto S350;
00292 am1 = *a-0.5e0-0.5e0;
00293 S230:
00294 if(*a <= amax) goto S240;
00295 d = 0.5e0+(0.5e0-xn/ *a);
00296 if(fabs(d) <= e2) goto S350;
00297 S240:
00298 if(*ierr >= 20) goto S330;
00299 *ierr += 1;
00300 gratio(a,&xn,&pn,&qn,&K8);
00301 if(pn == 0.0e0 || qn == 0.0e0) goto S350;
00302 r = rcomp(a,&xn);
00303 if(r == 0.0e0) goto S350;
00304 t = (*q-qn)/r;
00305 w = 0.5e0*(am1-xn);
00306 if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S250;
00307 *x = xn*(1.0e0-t);
00308 if(*x <= 0.0e0) goto S340;
00309 d = fabs(t);
00310 goto S260;
00311 S250:
00312 h = t*(1.0e0+w*t);
00313 *x = xn*(1.0e0-h);
00314 if(*x <= 0.0e0) goto S340;
00315 if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
00316 d = fabs(h);
00317 S260:
00318 xn = *x;
00319 if(d > tol) goto S230;
00320 if(d <= eps) return;
00321 if(fabs(*q-qn) <= tol**q) return;
00322 goto S230;
00323 S270:
00324
00325
00326
00327 *x = xmax;
00328 return;
00329 S280:
00330 if(*q < 0.9e0) goto S290;
00331 T9 = -*p;
00332 *x = -alnrel(&T9);
00333 return;
00334 S290:
00335 *x = -log(*q);
00336 return;
00337 S300:
00338
00339
00340
00341 *ierr = -2;
00342 return;
00343 S310:
00344 *ierr = -3;
00345 return;
00346 S320:
00347 *ierr = -4;
00348 return;
00349 S330:
00350 *ierr = -6;
00351 return;
00352 S340:
00353 *ierr = -7;
00354 return;
00355 S350:
00356 *x = xn;
00357 *ierr = -8;
00358 return;
00359 S360:
00360 *x = xmax;
00361 *ierr = -8;
00362 return;
00363 }