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 */