00001 #include "cdflib.h"
00002 void gratio(double *a,double *x,double *ans,double *qans,int *ind)
00003 /*
00004  ----------------------------------------------------------------------
00005         EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
00006                       P(A,X) AND Q(A,X)
00007  
00008                         ----------
00009  
00010      IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
00011      ARE NOT BOTH 0.
00012  
00013      ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
00014      P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
00015      IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
00016      POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
00017      IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
00018      6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
00019      IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.
00020  
00021      ERROR RETURN ...
00022         ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
00023      WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
00024      P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
00025      X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.
00026  ----------------------------------------------------------------------
00027      WRITTEN BY ALFRED H. MORRIS, JR.
00028         NAVAL SURFACE WEAPONS CENTER
00029         DAHLGREN, VIRGINIA
00030      --------------------
00031 */
00032 {
00033 static double alog10 = 2.30258509299405e0;
00034 static double d10 = -.185185185185185e-02;
00035 static double d20 = .413359788359788e-02;
00036 static double d30 = .649434156378601e-03;
00037 static double d40 = -.861888290916712e-03;
00038 static double d50 = -.336798553366358e-03;
00039 static double d60 = .531307936463992e-03;
00040 static double d70 = .344367606892378e-03;
00041 static double rt2pin = .398942280401433e0;
00042 static double rtpi = 1.77245385090552e0;
00043 static double third = .333333333333333e0;
00044 static double acc0[3] = {
00045     5.e-15,5.e-7,5.e-4
00046 };
00047 static double big[3] = {
00048     20.0e0,14.0e0,10.0e0
00049 };
00050 static double d0[13] = {
00051     .833333333333333e-01,-.148148148148148e-01,.115740740740741e-02,
00052     .352733686067019e-03,-.178755144032922e-03,.391926317852244e-04,
00053     -.218544851067999e-05,-.185406221071516e-05,.829671134095309e-06,
00054     -.176659527368261e-06,.670785354340150e-08,.102618097842403e-07,
00055     -.438203601845335e-08
00056 };
00057 static double d1[12] = {
00058     -.347222222222222e-02,.264550264550265e-02,-.990226337448560e-03,
00059     .205761316872428e-03,-.401877572016461e-06,-.180985503344900e-04,
00060     .764916091608111e-05,-.161209008945634e-05,.464712780280743e-08,
00061     .137863344691572e-06,-.575254560351770e-07,.119516285997781e-07
00062 };
00063 static double d2[10] = {
00064     -.268132716049383e-02,.771604938271605e-03,.200938786008230e-05,
00065     -.107366532263652e-03,.529234488291201e-04,-.127606351886187e-04,
00066     .342357873409614e-07,.137219573090629e-05,-.629899213838006e-06,
00067     .142806142060642e-06
00068 };
00069 static double d3[8] = {
00070     .229472093621399e-03,-.469189494395256e-03,.267720632062839e-03,
00071     -.756180167188398e-04,-.239650511386730e-06,.110826541153473e-04,
00072     -.567495282699160e-05,.142309007324359e-05
00073 };
00074 static double d4[6] = {
00075     .784039221720067e-03,-.299072480303190e-03,-.146384525788434e-05,
00076     .664149821546512e-04,-.396836504717943e-04,.113757269706784e-04
00077 };
00078 static double d5[4] = {
00079     -.697281375836586e-04,.277275324495939e-03,-.199325705161888e-03,
00080     .679778047793721e-04
00081 };
00082 static double d6[2] = {
00083     -.592166437353694e-03,.270878209671804e-03
00084 };
00085 static double e00[3] = {
00086     .25e-3,.25e-1,.14e0
00087 };
00088 static double x00[3] = {
00089     31.0e0,17.0e0,9.7e0
00090 };
00091 static int K1 = 1;
00092 static int K2 = 0;
00093 static double a2n,a2nm1,acc,am0,amn,an,an0,apn,b2n,b2nm1,c,c0,c1,c2,c3,c4,c5,c6,
00094     cma,e,e0,g,h,j,l,r,rta,rtx,s,sum,t,t1,tol,twoa,u,w,x0,y,z;
00095 static int i,iop,m,max,n;
00096 static double wk[20],T3;
00097 static int T4,T5;
00098 static double T6,T7;
00099 /*
00100      ..
00101      .. Executable Statements ..
00102 */
00103 /*
00104      --------------------
00105      ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
00106             FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
00107 */
00108     e = spmpar(&K1);
00109     if(*a < 0.0e0 || *x < 0.0e0) goto S430;
00110     if(*a == 0.0e0 && *x == 0.0e0) goto S430;
00111     if(*a**x == 0.0e0) goto S420;
00112     iop = *ind+1;
00113     if(iop != 1 && iop != 2) iop = 3;
00114     acc = fifdmax1(acc0[iop-1],e);
00115     e0 = e00[iop-1];
00116     x0 = x00[iop-1];
00117 /*
00118             SELECT THE APPROPRIATE ALGORITHM
00119 */
00120     if(*a >= 1.0e0) goto S10;
00121     if(*a == 0.5e0) goto S390;
00122     if(*x < 1.1e0) goto S160;
00123     t1 = *a*log(*x)-*x;
00124     u = *a*exp(t1);
00125     if(u == 0.0e0) goto S380;
00126     r = u*(1.0e0+gam1(a));
00127     goto S250;
00128 S10:
00129     if(*a >= big[iop-1]) goto S30;
00130     if(*a > *x || *x >= x0) goto S20;
00131     twoa = *a+*a;
00132     m = fifidint(twoa);
00133     if(twoa != (double)m) goto S20;
00134     i = m/2;
00135     if(*a == (double)i) goto S210;
00136     goto S220;
00137 S20:
00138     t1 = *a*log(*x)-*x;
00139     r = exp(t1)/Xgamm(a);
00140     goto S40;
00141 S30:
00142     l = *x/ *a;
00143     if(l == 0.0e0) goto S370;
00144     s = 0.5e0+(0.5e0-l);
00145     z = rlog(&l);
00146     if(z >= 700.0e0/ *a) goto S410;
00147     y = *a*z;
00148     rta = sqrt(*a);
00149     if(fabs(s) <= e0/rta) goto S330;
00150     if(fabs(s) <= 0.4e0) goto S270;
00151     t = pow(1.0e0/ *a,2.0);
00152     t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
00153     t1 -= y;
00154     r = rt2pin*rta*exp(t1);
00155 S40:
00156     if(r == 0.0e0) goto S420;
00157     if(*x <= fifdmax1(*a,alog10)) goto S50;
00158     if(*x < x0) goto S250;
00159     goto S100;
00160 S50:
00161 /*
00162                  TAYLOR SERIES FOR P/R
00163 */
00164     apn = *a+1.0e0;
00165     t = *x/apn;
00166     wk[0] = t;
00167     for(n=2; n<=20; n++) {
00168         apn += 1.0e0;
00169         t *= (*x/apn);
00170         if(t <= 1.e-3) goto S70;
00171         wk[n-1] = t;
00172     }
00173     n = 20;
00174 S70:
00175     sum = t;
00176     tol = 0.5e0*acc;
00177 S80:
00178     apn += 1.0e0;
00179     t *= (*x/apn);
00180     sum += t;
00181     if(t > tol) goto S80;
00182     max = n-1;
00183     for(m=1; m<=max; m++) {
00184         n -= 1;
00185         sum += wk[n-1];
00186     }
00187     *ans = r/ *a*(1.0e0+sum);
00188     *qans = 0.5e0+(0.5e0-*ans);
00189     return;
00190 S100:
00191 /*
00192                  ASYMPTOTIC EXPANSION
00193 */
00194     amn = *a-1.0e0;
00195     t = amn/ *x;
00196     wk[0] = t;
00197     for(n=2; n<=20; n++) {
00198         amn -= 1.0e0;
00199         t *= (amn/ *x);
00200         if(fabs(t) <= 1.e-3) goto S120;
00201         wk[n-1] = t;
00202     }
00203     n = 20;
00204 S120:
00205     sum = t;
00206 S130:
00207     if(fabs(t) <= acc) goto S140;
00208     amn -= 1.0e0;
00209     t *= (amn/ *x);
00210     sum += t;
00211     goto S130;
00212 S140:
00213     max = n-1;
00214     for(m=1; m<=max; m++) {
00215         n -= 1;
00216         sum += wk[n-1];
00217     }
00218     *qans = r/ *x*(1.0e0+sum);
00219     *ans = 0.5e0+(0.5e0-*qans);
00220     return;
00221 S160:
00222 /*
00223              TAYLOR SERIES FOR P(A,X)/X**A
00224 */
00225     an = 3.0e0;
00226     c = *x;
00227     sum = *x/(*a+3.0e0);
00228     tol = 3.0e0*acc/(*a+1.0e0);
00229 S170:
00230     an += 1.0e0;
00231     c = -(c*(*x/an));
00232     t = c/(*a+an);
00233     sum += t;
00234     if(fabs(t) > tol) goto S170;
00235     j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
00236     z = *a*log(*x);
00237     h = gam1(a);
00238     g = 1.0e0+h;
00239     if(*x < 0.25e0) goto S180;
00240     if(*a < *x/2.59e0) goto S200;
00241     goto S190;
00242 S180:
00243     if(z > -.13394e0) goto S200;
00244 S190:
00245     w = exp(z);
00246     *ans = w*g*(0.5e0+(0.5e0-j));
00247     *qans = 0.5e0+(0.5e0-*ans);
00248     return;
00249 S200:
00250     l = rexp(&z);
00251     w = 0.5e0+(0.5e0+l);
00252     *qans = (w*j-l)*g-h;
00253     if(*qans < 0.0e0) goto S380;
00254     *ans = 0.5e0+(0.5e0-*qans);
00255     return;
00256 S210:
00257 /*
00258              FINITE SUMS FOR Q WHEN A .GE. 1
00259                  AND 2*A IS AN INTEGER
00260 */
00261     sum = exp(-*x);
00262     t = sum;
00263     n = 1;
00264     c = 0.0e0;
00265     goto S230;
00266 S220:
00267     rtx = sqrt(*x);
00268     sum = erfc1(&K2,&rtx);
00269     t = exp(-*x)/(rtpi*rtx);
00270     n = 0;
00271     c = -0.5e0;
00272 S230:
00273     if(n == i) goto S240;
00274     n += 1;
00275     c += 1.0e0;
00276     t = *x*t/c;
00277     sum += t;
00278     goto S230;
00279 S240:
00280     *qans = sum;
00281     *ans = 0.5e0+(0.5e0-*qans);
00282     return;
00283 S250:
00284 /*
00285               CONTINUED FRACTION EXPANSION
00286 */
00287     tol = fifdmax1(5.0e0*e,acc);
00288     a2nm1 = a2n = 1.0e0;
00289     b2nm1 = *x;
00290     b2n = *x+(1.0e0-*a);
00291     c = 1.0e0;
00292 S260:
00293     a2nm1 = *x*a2n+c*a2nm1;
00294     b2nm1 = *x*b2n+c*b2nm1;
00295     am0 = a2nm1/b2nm1;
00296     c += 1.0e0;
00297     cma = c-*a;
00298     a2n = a2nm1+cma*a2n;
00299     b2n = b2nm1+cma*b2n;
00300     an0 = a2n/b2n;
00301     if(fabs(an0-am0) >= tol*an0) goto S260;
00302     *qans = r*an0;
00303     *ans = 0.5e0+(0.5e0-*qans);
00304     return;
00305 S270:
00306 /*
00307                 GENERAL TEMME EXPANSION
00308 */
00309     if(fabs(s) <= 2.0e0*e && *a*e*e > 3.28e-3) goto S430;
00310     c = exp(-y);
00311     T3 = sqrt(y);
00312     w = 0.5e0*erfc1(&K1,&T3);
00313     u = 1.0e0/ *a;
00314     z = sqrt(z+z);
00315     if(l < 1.0e0) z = -z;
00316     T4 = iop-2;
00317     if(T4 < 0) goto S280;
00318     else if(T4 == 0) goto S290;
00319     else  goto S300;
00320 S280:
00321     if(fabs(s) <= 1.e-3) goto S340;
00322     c0 = ((((((((((((d0[12]*z+d0[11])*z+d0[10])*z+d0[9])*z+d0[8])*z+d0[7])*z+d0[
00323       6])*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
00324     c1 = (((((((((((d1[11]*z+d1[10])*z+d1[9])*z+d1[8])*z+d1[7])*z+d1[6])*z+d1[5]
00325       )*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
00326     c2 = (((((((((d2[9]*z+d2[8])*z+d2[7])*z+d2[6])*z+d2[5])*z+d2[4])*z+d2[3])*z+
00327       d2[2])*z+d2[1])*z+d2[0])*z+d20;
00328     c3 = (((((((d3[7]*z+d3[6])*z+d3[5])*z+d3[4])*z+d3[3])*z+d3[2])*z+d3[1])*z+
00329       d3[0])*z+d30;
00330     c4 = (((((d4[5]*z+d4[4])*z+d4[3])*z+d4[2])*z+d4[1])*z+d4[0])*z+d40;
00331     c5 = (((d5[3]*z+d5[2])*z+d5[1])*z+d5[0])*z+d50;
00332     c6 = (d6[1]*z+d6[0])*z+d60;
00333     t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
00334     goto S310;
00335 S290:
00336     c0 = (((((d0[5]*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
00337     c1 = (((d1[3]*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
00338     c2 = d2[0]*z+d20;
00339     t = (c2*u+c1)*u+c0;
00340     goto S310;
00341 S300:
00342     t = ((d0[2]*z+d0[1])*z+d0[0])*z-third;
00343 S310:
00344     if(l < 1.0e0) goto S320;
00345     *qans = c*(w+rt2pin*t/rta);
00346     *ans = 0.5e0+(0.5e0-*qans);
00347     return;
00348 S320:
00349     *ans = c*(w-rt2pin*t/rta);
00350     *qans = 0.5e0+(0.5e0-*ans);
00351     return;
00352 S330:
00353 /*
00354                TEMME EXPANSION FOR L = 1
00355 */
00356     if(*a*e*e > 3.28e-3) goto S430;
00357     c = 0.5e0+(0.5e0-y);
00358     w = (0.5e0-sqrt(y)*(0.5e0+(0.5e0-y/3.0e0))/rtpi)/c;
00359     u = 1.0e0/ *a;
00360     z = sqrt(z+z);
00361     if(l < 1.0e0) z = -z;
00362     T5 = iop-2;
00363     if(T5 < 0) goto S340;
00364     else if(T5 == 0) goto S350;
00365     else  goto S360;
00366 S340:
00367     c0 = ((((((d0[6]*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-
00368       third;
00369     c1 = (((((d1[5]*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
00370     c2 = ((((d2[4]*z+d2[3])*z+d2[2])*z+d2[1])*z+d2[0])*z+d20;
00371     c3 = (((d3[3]*z+d3[2])*z+d3[1])*z+d3[0])*z+d30;
00372     c4 = (d4[1]*z+d4[0])*z+d40;
00373     c5 = (d5[1]*z+d5[0])*z+d50;
00374     c6 = d6[0]*z+d60;
00375     t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
00376     goto S310;
00377 S350:
00378     c0 = (d0[1]*z+d0[0])*z-third;
00379     c1 = d1[0]*z+d10;
00380     t = (d20*u+c1)*u+c0;
00381     goto S310;
00382 S360:
00383     t = d0[0]*z-third;
00384     goto S310;
00385 S370:
00386 /*
00387                      SPECIAL CASES
00388 */
00389     *ans = 0.0e0;
00390     *qans = 1.0e0;
00391     return;
00392 S380:
00393     *ans = 1.0e0;
00394     *qans = 0.0e0;
00395     return;
00396 S390:
00397     if(*x >= 0.25e0) goto S400;
00398     T6 = sqrt(*x);
00399     *ans = erf1(&T6);
00400     *qans = 0.5e0+(0.5e0-*ans);
00401     return;
00402 S400:
00403     T7 = sqrt(*x);
00404     *qans = erfc1(&K2,&T7);
00405     *ans = 0.5e0+(0.5e0-*qans);
00406     return;
00407 S410:
00408     if(fabs(s) <= 2.0e0*e) goto S430;
00409 S420:
00410     if(*x <= *a) goto S370;
00411     goto S380;
00412 S430:
00413 /*
00414                      ERROR RETURN
00415 */
00416     *ans = 2.0e0;
00417     return;
00418 } /* END */