00001 #include "cdflib.h"
00002 void gaminv(double *a,double *x,double *x0,double *p,double *q,
00003             int *ierr)
00004 /*
00005  ----------------------------------------------------------------------
00006             INVERSE INCOMPLETE GAMMA RATIO FUNCTION
00007  
00008      GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1.
00009      THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER
00010      ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X
00011      TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE
00012      PARTICULAR COMPUTER ARITHMETIC BEING USED.
00013  
00014                       ------------
00015  
00016      X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0,
00017      AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT
00018      NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN
00019      A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE
00020      IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X.
00021  
00022      X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER
00023      DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET
00024      X0 .LE. 0.
00025  
00026      IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
00027      WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING
00028      VALUES ...
00029  
00030        IERR =  0    THE SOLUTION WAS OBTAINED. ITERATION WAS
00031                     NOT USED.
00032        IERR.GT.0    THE SOLUTION WAS OBTAINED. IERR ITERATIONS
00033                     WERE PERFORMED.
00034        IERR = -2    (INPUT ERROR) A .LE. 0
00035        IERR = -3    NO SOLUTION WAS OBTAINED. THE RATIO Q/A
00036                     IS TOO LARGE.
00037        IERR = -4    (INPUT ERROR) P + Q .NE. 1
00038        IERR = -6    20 ITERATIONS WERE PERFORMED. THE MOST
00039                     RECENT VALUE OBTAINED FOR X IS GIVEN.
00040                     THIS CANNOT OCCUR IF X0 .LE. 0.
00041        IERR = -7    ITERATION FAILED. NO VALUE IS GIVEN FOR X.
00042                     THIS MAY OCCUR WHEN X IS APPROXIMATELY 0.
00043        IERR = -8    A VALUE FOR X HAS BEEN OBTAINED, BUT THE
00044                     ROUTINE IS NOT CERTAIN OF ITS ACCURACY.
00045                     ITERATION CANNOT BE PERFORMED IN THIS
00046                     CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY
00047                     WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS
00048                     POSITIVE THEN THIS CAN OCCUR WHEN A IS
00049                     EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY
00050                     LARGE (SAY A .GE. 1.E20).
00051  ----------------------------------------------------------------------
00052      WRITTEN BY ALFRED H. MORRIS, JR.
00053         NAVAL SURFACE WEAPONS CENTER
00054         DAHLGREN, VIRGINIA
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      .. Executable Statements ..
00095 */
00096 /*
00097      ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
00098             E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
00099             XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
00100             LARGEST POSITIVE NUMBER.
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         SELECTION OF THE INITIAL APPROXIMATION XN OF X
00122                        WHEN A .LT. 1
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         SELECTION OF THE INITIAL APPROXIMATION XN OF X
00184                        WHEN A .GT. 1
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                  SCHRODER ITERATION USING P
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                  SCHRODER ITERATION USING Q
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                        SPECIAL CASES
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                        ERROR RETURN
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 } /* END */