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