00015 {
00016 static double c0 = .833333333333333e-01;
00017 static double c1 = -.277777777760991e-02;
00018 static double c2 = .793650666825390e-03;
00019 static double c3 = -.595202931351870e-03;
00020 static double c4 = .837308034031215e-03;
00021 static double c5 = -.165322962780713e-02;
00022 static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;
00023 /*
00024      ..
00025      .. Executable Statements ..
00026 */
00027     if(*a <= *b) goto S10;
00028     h = *b/ *a;
00029     c = 1.0e0/(1.0e0+h);
00030     x = h/(1.0e0+h);
00031     d = *a+(*b-0.5e0);
00032     goto S20;
00033 S10:
00034     h = *a/ *b;
00035     c = h/(1.0e0+h);
00036     x = 1.0e0/(1.0e0+h);
00037     d = *b+(*a-0.5e0);
00038 S20:
00039 /*
00040                 SET SN = (1 - X**N)/(1 - X)
00041 */
00042     x2 = x*x;
00043     s3 = 1.0e0+(x+x2);
00044     s5 = 1.0e0+(x+x2*s3);
00045     s7 = 1.0e0+(x+x2*s5);
00046     s9 = 1.0e0+(x+x2*s7);
00047     s11 = 1.0e0+(x+x2*s9);
00048 /*
00049                 SET W = DEL(B) - DEL(A + B)
00050 */
00051     t = pow(1.0e0/ *b,2.0);
00052     w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
00053     w *= (c/ *b);
00054 /*
00055                     COMBINE THE RESULTS
00056 */
00057     T1 = *a/ *b;
00058     u = d*alnrel(&T1);
00059     v = *a*(log(*b)-1.0e0);
00060     if(u <= v) goto S30;
00061     algdiv = w-v-u;
00062     return algdiv;
00063 S30:
00064     algdiv = w-u-v;
00065     return algdiv;
00066 } /* END */

00036 {
00037 #define hln2pi 0.91893853320467274178e0
00038 static double coef[5] = {
00039     0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
00040     -0.594997310889e-3,0.8065880899e-3
00041 };
00042 static double scoefd[4] = {
00043     0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
00044     0.1000000000000000000e1
00045 };
00046 static double scoefn[9] = {
00047     0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
00048     0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
00049     0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
00050 };
00051 static int K1 = 9;
00052 static int K3 = 4;
00053 static int K5 = 5;
00054 static double alngam,offset,prod,xx;
00055 static int i,n;
00056 static double T2,T4,T6;
00057 /*
00058      ..
00059      .. Executable Statements ..
00060 */
00061     if(!(*x <= 6.0e0)) goto S70;
00062     prod = 1.0e0;
00063     xx = *x;
00064     if(!(*x > 3.0e0)) goto S30;
00065 S10:
00066     if(!(xx > 3.0e0)) goto S20;
00067     xx -= 1.0e0;
00068     prod *= xx;
00069     goto S10;
00070 S30:
00071 S20:
00072     if(!(*x < 2.0e0)) goto S60;
00073 S40:
00074     if(!(xx < 2.0e0)) goto S50;
00075     prod /= xx;
00076     xx += 1.0e0;
00077     goto S40;
00078 S60:
00079 S50:
00080     T2 = xx-2.0e0;
00081     T4 = xx-2.0e0;
00082     alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
00083 /*
00084      COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
00085 */
00086     alngam *= prod;
00087     alngam = log(alngam);
00088     goto S110;
00089 S70:
00090     offset = hln2pi;
00091 /*
00092      IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
00093 */
00094     n = fifidint(12.0e0-*x);
00095     if(!(n > 0)) goto S90;
00096     prod = 1.0e0;
00097     for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));
00098     offset -= log(prod);
00099     xx = *x+(double)n;
00100     goto S100;
00101 S90:
00102     xx = *x;
00103 S100:
00104 /*
00105      COMPUTE POWER SERIES
00106 */
00107     T6 = 1.0e0/pow(xx,2.0);
00108     alngam = devlpl(coef,&K5,&T6)/xx;
00109     alngam += (offset+(xx-0.5e0)*log(xx)-xx);
00110 S110:
00111     return alngam;
00112 #undef hln2pi
00113 } /* END */

00008 {
00009 static double p1 = -.129418923021993e+01;
00010 static double p2 = .405303492862024e+00;
00011 static double p3 = -.178874546012214e-01;
00012 static double q1 = -.162752256355323e+01;
00013 static double q2 = .747811014037616e+00;
00014 static double q3 = -.845104217945565e-01;
00015 static double alnrel,t,t2,w,x;
00016 /*
00017      ..
00018      .. Executable Statements ..
00019 */
00020     if(fabs(*a) > 0.375e0) goto S10;
00021     t = *a/(*a+2.0e0);
00022     t2 = t*t;
00023     w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);
00024     alnrel = 2.0e0*t*w;
00025     return alnrel;
00026 S10:
00027     x = 1.e0+*a;
00028     alnrel = log(x);
00029     return alnrel;
00030 } /* END */

00010 {
00011 static double g = .577215664901533e0;
00012 static double apser,aj,bx,c,j,s,t,tol;
00013 /*
00014      ..
00015      .. Executable Statements ..
00016 */
00017     bx = *b**x;
00018     t = *x-bx;
00019     if(*b**eps > 2.e-2) goto S10;
00020     c = log(*x)+psi(b)+g+t;
00021     goto S20;
00022 S10:
00023     c = log(bx)+g+t;
00024 S20:
00025     tol = 5.0e0**eps*fabs(c);
00026     j = 1.0e0;
00027     s = 0.0e0;
00028 S30:
00029     j += 1.0e0;
00030     t *= (*x-bx/j);
00031     aj = t/j;
00032     s += aj;
00033     if(fabs(aj) > tol) goto S30;
00034     apser = -(*a*(c+s));
00035     return apser;
00036 } /* END */

00011 {
00012 static double e0 = 1.12837916709551e0;
00013 static double e1 = .353553390593274e0;
00014 static int num = 20;
00015 /*
00016 ------------------------
00017      ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
00018             ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
00019             THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
00020 ------------------------
00021      E0 = 2/SQRT(PI)
00022      E1 = 2**(-3/2)
00023 ------------------------
00024 */
00025 static int K3 = 1;
00026 static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
00027     z2,zn,znm1;
00028 static int i,im1,imj,j,m,mm1,mmj,n,np1;
00029 static double a0[21],b0[21],c[21],d[21],T1,T2;
00030 /*
00031      ..
00032      .. Executable Statements ..
00033 */
00034     basym = 0.0e0;
00035     if(*a >= *b) goto S10;
00036     h = *a/ *b;
00037     r0 = 1.0e0/(1.0e0+h);
00038     r1 = (*b-*a)/ *b;
00039     w0 = 1.0e0/sqrt(*a*(1.0e0+h));
00040     goto S20;
00041 S10:
00042     h = *b/ *a;
00043     r0 = 1.0e0/(1.0e0+h);
00044     r1 = (*b-*a)/ *a;
00045     w0 = 1.0e0/sqrt(*b*(1.0e0+h));
00046 S20:
00047     T1 = -(*lambda/ *a);
00048     T2 = *lambda/ *b;
00049     f = *a*rlog1(&T1)+*b*rlog1(&T2);
00050     t = exp(-f);
00051     if(t == 0.0e0) return basym;
00052     z0 = sqrt(f);
00053     z = 0.5e0*(z0/e1);
00054     z2 = f+f;
00055     a0[0] = 2.0e0/3.0e0*r1;
00056     c[0] = -(0.5e0*a0[0]);
00057     d[0] = -c[0];
00058     j0 = 0.5e0/e0*erfc1(&K3,&z0);
00059     j1 = e1;
00060     sum = j0+d[0]*w0*j1;
00061     s = 1.0e0;
00062     h2 = h*h;
00063     hn = 1.0e0;
00064     w = w0;
00065     znm1 = z;
00066     zn = z2;
00067     for(n=2; n<=num; n+=2) {
00068         hn = h2*hn;
00069         a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
00070         np1 = n+1;
00071         s += hn;
00072         a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
00073         for(i=n; i<=np1; i++) {
00074             r = -(0.5e0*((double)i+1.0e0));
00075             b0[0] = r*a0[0];
00076             for(m=2; m<=i; m++) {
00077                 bsum = 0.0e0;
00078                 mm1 = m-1;
00079                 for(j=1; j<=mm1; j++) {
00080                     mmj = m-j;
00081                     bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
00082                 }
00083                 b0[m-1] = r*a0[m-1]+bsum/(double)m;
00084             }
00085             c[i-1] = b0[i-1]/((double)i+1.0e0);
00086             dsum = 0.0e0;
00087             im1 = i-1;
00088             for(j=1; j<=im1; j++) {
00089                 imj = i-j;
00090                 dsum += (d[imj-1]*c[j-1]);
00091             }
00092             d[i-1] = -(dsum+c[i-1]);
00093         }
00094         j0 = e1*znm1+((double)n-1.0e0)*j0;
00095         j1 = e1*zn+(double)n*j1;
00096         znm1 = z2*znm1;
00097         zn = z2*zn;
00098         w = w0*w;
00099         t0 = d[n-1]*w*j0;
00100         w = w0*w;
00101         t1 = d[np1-1]*w*j1;
00102         sum += (t0+t1);
00103         if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
00104     }
00105 S80:
00106     u = exp(-bcorr(a,b));
00107     basym = e0*t*u*sum;
00108     return basym;
00109 } /* END */

00012 {
00013 static double c0 = .833333333333333e-01;
00014 static double c1 = -.277777777760991e-02;
00015 static double c2 = .793650666825390e-03;
00016 static double c3 = -.595202931351870e-03;
00017 static double c4 = .837308034031215e-03;
00018 static double c5 = -.165322962780713e-02;
00019 static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;
00020 /*
00021      ..
00022      .. Executable Statements ..
00023 */
00024     a = fifdmin1(*a0,*b0);
00025     b = fifdmax1(*a0,*b0);
00026     h = a/b;
00027     c = h/(1.0e0+h);
00028     x = 1.0e0/(1.0e0+h);
00029     x2 = x*x;
00030 /*
00031                 SET SN = (1 - X**N)/(1 - X)
00032 */
00033     s3 = 1.0e0+(x+x2);
00034     s5 = 1.0e0+(x+x2*s3);
00035     s7 = 1.0e0+(x+x2*s5);
00036     s9 = 1.0e0+(x+x2*s7);
00037     s11 = 1.0e0+(x+x2*s9);
00038 /*
00039                 SET W = DEL(B) - DEL(A + B)
00040 */
00041     t = pow(1.0e0/b,2.0);
00042     w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
00043     w *= (c/b);
00044 /*
00045                    COMPUTE  DEL(A) + W
00046 */
00047     t = pow(1.0e0/a,2.0);
00048     bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w;
00049     return bcorr;
00050 } /* END */

00010 {
00011 static double e = .918938533204673e0;
00012 static double betaln,a,b,c,h,u,v,w,z;
00013 static int i,n;
00014 static double T1;
00015 /*
00016      ..
00017      .. Executable Statements ..
00018 */
00019     a = fifdmin1(*a0,*b0);
00020     b = fifdmax1(*a0,*b0);
00021     if(a >= 8.0e0) goto S100;
00022     if(a >= 1.0e0) goto S20;
00023 /*
00024 -----------------------------------------------------------------------
00025                    PROCEDURE WHEN A .LT. 1
00026 -----------------------------------------------------------------------
00027 */
00028     if(b >= 8.0e0) goto S10;
00029     T1 = a+b;
00030     betaln = gamln(&a)+(gamln(&b)-gamln(&T1));
00031     return betaln;
00032 S10:
00033     betaln = gamln(&a)+algdiv(&a,&b);
00034     return betaln;
00035 S20:
00036 /*
00037 -----------------------------------------------------------------------
00038                 PROCEDURE WHEN 1 .LE. A .LT. 8
00039 -----------------------------------------------------------------------
00040 */
00041     if(a > 2.0e0) goto S40;
00042     if(b > 2.0e0) goto S30;
00043     betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b);
00044     return betaln;
00045 S30:
00046     w = 0.0e0;
00047     if(b < 8.0e0) goto S60;
00048     betaln = gamln(&a)+algdiv(&a,&b);
00049     return betaln;
00050 S40:
00051 /*
00052                 REDUCTION OF A WHEN B .LE. 1000
00053 */
00054     if(b > 1000.0e0) goto S80;
00055     n = a-1.0e0;
00056     w = 1.0e0;
00057     for(i=1; i<=n; i++) {
00058         a -= 1.0e0;
00059         h = a/b;
00060         w *= (h/(1.0e0+h));
00061     }
00062     w = log(w);
00063     if(b < 8.0e0) goto S60;
00064     betaln = w+gamln(&a)+algdiv(&a,&b);
00065     return betaln;
00066 S60:
00067 /*
00068                  REDUCTION OF B WHEN B .LT. 8
00069 */
00070     n = b-1.0e0;
00071     z = 1.0e0;
00072     for(i=1; i<=n; i++) {
00073         b -= 1.0e0;
00074         z *= (b/(a+b));
00075     }
00076     betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
00077     return betaln;
00078 S80:
00079 /*
00080                 REDUCTION OF A WHEN B .GT. 1000
00081 */
00082     n = a-1.0e0;
00083     w = 1.0e0;
00084     for(i=1; i<=n; i++) {
00085         a -= 1.0e0;
00086         w *= (a/(1.0e0+a/b));
00087     }
00088     betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
00089     return betaln;
00090 S100:
00091 /*
00092 -----------------------------------------------------------------------
00093                    PROCEDURE WHEN A .GE. 8
00094 -----------------------------------------------------------------------
00095 */
00096     w = bcorr(&a,&b);
00097     h = a/b;
00098     c = h/(1.0e0+h);
00099     u = -((a-0.5e0)*log(c));
00100     v = b*alnrel(&h);
00101     if(u <= v) goto S110;
00102     betaln = -(0.5e0*log(b))+e+w-v-u;
00103     return betaln;
00104 S110:
00105     betaln = -(0.5e0*log(b))+e+w-u-v;
00106     return betaln;
00107 } /* END */

00010 {
00011 static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;
00012 /*
00013      ..
00014      .. Executable Statements ..
00015 */
00016     bfrac = brcomp(a,b,x,y);
00017     if(bfrac == 0.0e0) return bfrac;
00018     c = 1.0e0+*lambda;
00019     c0 = *b/ *a;
00020     c1 = 1.0e0+1.0e0/ *a;
00021     yp1 = *y+1.0e0;
00022     n = 0.0e0;
00023     p = 1.0e0;
00024     s = *a+1.0e0;
00025     an = 0.0e0;
00026     bn = anp1 = 1.0e0;
00027     bnp1 = c/c1;
00028     r = c1/c;
00029 S10:
00030 /*
00031         CONTINUED FRACTION CALCULATION
00032 */
00033     n += 1.0e0;
00034     t = n/ *a;
00035     w = n*(*b-n)**x;
00036     e = *a/s;
00037     alpha = p*(p+c0)*e*e*(w**x);
00038     e = (1.0e0+t)/(c1+t+t);
00039     beta = n+w/s+e*(c+n*yp1);
00040     p = 1.0e0+t;
00041     s += 2.0e0;
00042 /*
00043         UPDATE AN, BN, ANP1, AND BNP1
00044 */
00045     t = alpha*an+beta*anp1;
00046     an = anp1;
00047     anp1 = t;
00048     t = alpha*bn+beta*bnp1;
00049     bn = bnp1;
00050     bnp1 = t;
00051     r0 = r;
00052     r = anp1/bnp1;
00053     if(fabs(r-r0) <= *eps*r) goto S20;
00054 /*
00055         RESCALE AN, BN, ANP1, AND BNP1
00056 */
00057     an /= bnp1;
00058     bn /= bnp1;
00059     anp1 = r;
00060     bnp1 = 1.0e0;
00061     goto S10;
00062 S20:
00063 /*
00064                  TERMINATION
00065 */
00066     bfrac *= r;
00067     return bfrac;
00068 } /* END */

00012 {
00013 static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
00014 static int i,n,nm1;
00015 static double c[30],d[30],T1;
00016 /*
00017      ..
00018      .. Executable Statements ..
00019 */
00020     bm1 = *b-0.5e0-0.5e0;
00021     nu = *a+0.5e0*bm1;
00022     if(*y > 0.375e0) goto S10;
00023     T1 = -*y;
00024     lnx = alnrel(&T1);
00025     goto S20;
00026 S10:
00027     lnx = log(*x);
00028 S20:
00029     z = -(nu*lnx);
00030     if(*b*z == 0.0e0) goto S70;
00031 /*
00032                  COMPUTATION OF THE EXPANSION
00033                  SET R = EXP(-Z)*Z**B/GAMMA(B)
00034 */
00035     r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
00036     r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
00037     u = algdiv(b,a)+*b*log(nu);
00038     u = r*exp(-u);
00039     if(u == 0.0e0) goto S70;
00040     grat1(b,&z,&r,&p,&q,eps);
00041     v = 0.25e0*pow(1.0e0/nu,2.0);
00042     t2 = 0.25e0*lnx*lnx;
00043     l = *w/u;
00044     j = q/r;
00045     sum = j;
00046     t = cn = 1.0e0;
00047     n2 = 0.0e0;
00048     for(n=1; n<=30; n++) {
00049         bp2n = *b+n2;
00050         j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
00051         n2 += 2.0e0;
00052         t *= t2;
00053         cn /= (n2*(n2+1.0e0));
00054         c[n-1] = cn;
00055         s = 0.0e0;
00056         if(n == 1) goto S40;
00057         nm1 = n-1;
00058         coef = *b-(double)n;
00059         for(i=1; i<=nm1; i++) {
00060             s += (coef*c[i-1]*d[n-i-1]);
00061             coef += *b;
00062         }
00063 S40:
00064         d[n-1] = bm1*cn+s/(double)n;
00065         dj = d[n-1]*j;
00066         sum += dj;
00067         if(sum <= 0.0e0) goto S70;
00068         if(fabs(dj) <= *eps*(sum+l)) goto S60;
00069     }
00070 S60:
00071 /*
00072                     ADD THE RESULTS TO W
00073 */
00074     *ierr = 0;
00075     *w += (u*sum);
00076     return;
00077 S70:
00078 /*
00079                THE EXPANSION CANNOT BE COMPUTED
00080 */
00081     *ierr = 1;
00082     return;
00083 } /* END */

00009 {
00010 static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
00011 static int i,m;
00012 /*
00013      ..
00014      .. Executable Statements ..
00015 */
00016     bpser = 0.0e0;
00017     if(*x == 0.0e0) return bpser;
00018 /*
00019 -----------------------------------------------------------------------
00020             COMPUTE THE FACTOR X**A/(A*BETA(A,B))
00021 -----------------------------------------------------------------------
00022 */
00023     a0 = fifdmin1(*a,*b);
00024     if(a0 < 1.0e0) goto S10;
00025     z = *a*log(*x)-betaln(a,b);
00026     bpser = exp(z)/ *a;
00027     goto S100;
00028 S10:
00029     b0 = fifdmax1(*a,*b);
00030     if(b0 >= 8.0e0) goto S90;
00031     if(b0 > 1.0e0) goto S40;
00032 /*
00033             PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
00034 */
00035     bpser = pow(*x,*a);
00036     if(bpser == 0.0e0) return bpser;
00037     apb = *a+*b;
00038     if(apb > 1.0e0) goto S20;
00039     z = 1.0e0+gam1(&apb);
00040     goto S30;
00041 S20:
00042     u = *a+*b-1.e0;
00043     z = (1.0e0+gam1(&u))/apb;
00044 S30:
00045     c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
00046     bpser *= (c*(*b/apb));
00047     goto S100;
00048 S40:
00049 /*
00050          PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
00051 */
00052     u = gamln1(&a0);
00053     m = b0-1.0e0;
00054     if(m < 1) goto S60;
00055     c = 1.0e0;
00056     for(i=1; i<=m; i++) {
00057         b0 -= 1.0e0;
00058         c *= (b0/(a0+b0));
00059     }
00060     u = log(c)+u;
00061 S60:
00062     z = *a*log(*x)-u;
00063     b0 -= 1.0e0;
00064     apb = a0+b0;
00065     if(apb > 1.0e0) goto S70;
00066     t = 1.0e0+gam1(&apb);
00067     goto S80;
00068 S70:
00069     u = a0+b0-1.e0;
00070     t = (1.0e0+gam1(&u))/apb;
00071 S80:
00072     bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
00073     goto S100;
00074 S90:
00075 /*
00076             PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
00077 */
00078     u = gamln1(&a0)+algdiv(&a0,&b0);
00079     z = *a*log(*x)-u;
00080     bpser = a0/ *a*exp(z);
00081 S100:
00082     if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
00083 /*
00084 -----------------------------------------------------------------------
00085                      COMPUTE THE SERIES
00086 -----------------------------------------------------------------------
00087 */
00088     sum = n = 0.0e0;
00089     c = 1.0e0;
00090     tol = *eps/ *a;
00091 S110:
00092     n += 1.0e0;
00093     c *= ((0.5e0+(0.5e0-*b/n))**x);
00094     w = c/(*a+n);
00095     sum += w;
00096     if(fabs(w) > tol) goto S110;
00097     bpser *= (1.0e0+*a*sum);
00098     return bpser;
00099 } /* END */

00038 {
00039 static int K1 = 1;
00040 static double a0,b0,eps,lambda,t,x0,y0,z;
00041 static int ierr1,ind,n;
00042 static double T2,T3,T4,T5;
00043 /*
00044      ..
00045      .. Executable Statements ..
00046 */
00047 /*
00048      ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
00049             FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
00050 */
00051     eps = spmpar(&K1);
00052     *w = *w1 = 0.0e0;
00053     if(*a < 0.0e0 || *b < 0.0e0) goto S270;
00054     if(*a == 0.0e0 && *b == 0.0e0) goto S280;
00055     if(*x < 0.0e0 || *x > 1.0e0) goto S290;
00056     if(*y < 0.0e0 || *y > 1.0e0) goto S300;
00057     z = *x+*y-0.5e0-0.5e0;
00058     if(fabs(z) > 3.0e0*eps) goto S310;
00059     *ierr = 0;
00060     if(*x == 0.0e0) goto S210;
00061     if(*y == 0.0e0) goto S230;
00062     if(*a == 0.0e0) goto S240;
00063     if(*b == 0.0e0) goto S220;
00064     eps = fifdmax1(eps,1.e-15);
00065     if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
00066     ind = 0;
00067     a0 = *a;
00068     b0 = *b;
00069     x0 = *x;
00070     y0 = *y;
00071     if(fifdmin1(a0,b0) > 1.0e0) goto S40;
00072 /*
00073              PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
00074 */
00075     if(*x <= 0.5e0) goto S10;
00076     ind = 1;
00077     a0 = *b;
00078     b0 = *a;
00079     x0 = *y;
00080     y0 = *x;
00081 S10:
00082     if(b0 < fifdmin1(eps,eps*a0)) goto S90;
00083     if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
00084     if(fifdmax1(a0,b0) > 1.0e0) goto S20;
00085     if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
00086     if(pow(x0,a0) <= 0.9e0) goto S110;
00087     if(x0 >= 0.3e0) goto S120;
00088     n = 20;
00089     goto S140;
00090 S20:
00091     if(b0 <= 1.0e0) goto S110;
00092     if(x0 >= 0.3e0) goto S120;
00093     if(x0 >= 0.1e0) goto S30;
00094     if(pow(x0*b0,a0) <= 0.7e0) goto S110;
00095 S30:
00096     if(b0 > 15.0e0) goto S150;
00097     n = 20;
00098     goto S140;
00099 S40:
00100 /*
00101              PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
00102 */
00103     if(*a > *b) goto S50;
00104     lambda = *a-(*a+*b)**x;
00105     goto S60;
00106 S50:
00107     lambda = (*a+*b)**y-*b;
00108 S60:
00109     if(lambda >= 0.0e0) goto S70;
00110     ind = 1;
00111     a0 = *b;
00112     b0 = *a;
00113     x0 = *y;
00114     y0 = *x;
00115     lambda = fabs(lambda);
00116 S70:
00117     if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
00118     if(b0 < 40.0e0) goto S160;
00119     if(a0 > b0) goto S80;
00120     if(a0 <= 100.0e0) goto S130;
00121     if(lambda > 0.03e0*a0) goto S130;
00122     goto S200;
00123 S80:
00124     if(b0 <= 100.0e0) goto S130;
00125     if(lambda > 0.03e0*b0) goto S130;
00126     goto S200;
00127 S90:
00128 /*
00129             EVALUATION OF THE APPROPRIATE ALGORITHM
00130 */
00131     *w = fpser(&a0,&b0,&x0,&eps);
00132     *w1 = 0.5e0+(0.5e0-*w);
00133     goto S250;
00134 S100:
00135     *w1 = apser(&a0,&b0,&x0,&eps);
00136     *w = 0.5e0+(0.5e0-*w1);
00137     goto S250;
00138 S110:
00139     *w = bpser(&a0,&b0,&x0,&eps);
00140     *w1 = 0.5e0+(0.5e0-*w);
00141     goto S250;
00142 S120:
00143     *w1 = bpser(&b0,&a0,&y0,&eps);
00144     *w = 0.5e0+(0.5e0-*w1);
00145     goto S250;
00146 S130:
00147     T2 = 15.0e0*eps;
00148     *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
00149     *w1 = 0.5e0+(0.5e0-*w);
00150     goto S250;
00151 S140:
00152     *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
00153     b0 += (double)n;
00154 S150:
00155     T3 = 15.0e0*eps;
00156     bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
00157     *w = 0.5e0+(0.5e0-*w1);
00158     goto S250;
00159 S160:
00160     n = b0;
00161     b0 -= (double)n;
00162     if(b0 != 0.0e0) goto S170;
00163     n -= 1;
00164     b0 = 1.0e0;
00165 S170:
00166     *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
00167     if(x0 > 0.7e0) goto S180;
00168     *w += bpser(&a0,&b0,&x0,&eps);
00169     *w1 = 0.5e0+(0.5e0-*w);
00170     goto S250;
00171 S180:
00172     if(a0 > 15.0e0) goto S190;
00173     n = 20;
00174     *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
00175     a0 += (double)n;
00176 S190:
00177     T4 = 15.0e0*eps;
00178     bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
00179     *w1 = 0.5e0+(0.5e0-*w);
00180     goto S250;
00181 S200:
00182     T5 = 100.0e0*eps;
00183     *w = basym(&a0,&b0,&lambda,&T5);
00184     *w1 = 0.5e0+(0.5e0-*w);
00185     goto S250;
00186 S210:
00187 /*
00188                TERMINATION OF THE PROCEDURE
00189 */
00190     if(*a == 0.0e0) goto S320;
00191 S220:
00192     *w = 0.0e0;
00193     *w1 = 1.0e0;
00194     return;
00195 S230:
00196     if(*b == 0.0e0) goto S330;
00197 S240:
00198     *w = 1.0e0;
00199     *w1 = 0.0e0;
00200     return;
00201 S250:
00202     if(ind == 0) return;
00203     t = *w;
00204     *w = *w1;
00205     *w1 = t;
00206     return;
00207 S260:
00208 /*
00209            PROCEDURE FOR A AND B .LT. 1.E-3*EPS
00210 */
00211     *w = *b/(*a+*b);
00212     *w1 = *a/(*a+*b);
00213     return;
00214 S270:
00215 /*
00216                        ERROR RETURN
00217 */
00218     *ierr = 1;
00219     return;
00220 S280:
00221     *ierr = 2;
00222     return;
00223 S290:
00224     *ierr = 3;
00225     return;
00226 S300:
00227     *ierr = 4;
00228     return;
00229 S310:
00230     *ierr = 5;
00231     return;
00232 S320:
00233     *ierr = 6;
00234     return;
00235 S330:
00236     *ierr = 7;
00237     return;
00238 } /* END */

00008 {
00009 static double Const = .398942280401433e0;
00010 static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
00011 static int i,n;
00012 /*
00013 -----------------
00014      CONST = 1/SQRT(2*PI)
00015 -----------------
00016 */
00017 static double T1,T2,T3,T4;
00018 /*
00019      ..
00020      .. Executable Statements ..
00021 */
00022     a0 = fifdmin1(*a,*b);
00023     if(a0 >= 8.0e0) goto S130;
00024     if(*x > 0.375e0) goto S10;
00025     lnx = log(*x);
00026     T1 = -*x;
00027     lny = alnrel(&T1);
00028     goto S30;
00029 S10:
00030     if(*y > 0.375e0) goto S20;
00031     T2 = -*y;
00032     lnx = alnrel(&T2);
00033     lny = log(*y);
00034     goto S30;
00035 S20:
00036     lnx = log(*x);
00037     lny = log(*y);
00038 S30:
00039     z = *a*lnx+*b*lny;
00040     if(a0 < 1.0e0) goto S40;
00041     z -= betaln(a,b);
00042     brcmp1 = esum(mu,&z);
00043     return brcmp1;
00044 S40:
00045 /*
00046 -----------------------------------------------------------------------
00047               PROCEDURE FOR A .LT. 1 OR B .LT. 1
00048 -----------------------------------------------------------------------
00049 */
00050     b0 = fifdmax1(*a,*b);
00051     if(b0 >= 8.0e0) goto S120;
00052     if(b0 > 1.0e0) goto S70;
00053 /*
00054                    ALGORITHM FOR B0 .LE. 1
00055 */
00056     brcmp1 = esum(mu,&z);
00057     if(brcmp1 == 0.0e0) return brcmp1;
00058     apb = *a+*b;
00059     if(apb > 1.0e0) goto S50;
00060     z = 1.0e0+gam1(&apb);
00061     goto S60;
00062 S50:
00063     u = *a+*b-1.e0;
00064     z = (1.0e0+gam1(&u))/apb;
00065 S60:
00066     c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
00067     brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
00068     return brcmp1;
00069 S70:
00070 /*
00071                 ALGORITHM FOR 1 .LT. B0 .LT. 8
00072 */
00073     u = gamln1(&a0);
00074     n = b0-1.0e0;
00075     if(n < 1) goto S90;
00076     c = 1.0e0;
00077     for(i=1; i<=n; i++) {
00078         b0 -= 1.0e0;
00079         c *= (b0/(a0+b0));
00080     }
00081     u = log(c)+u;
00082 S90:
00083     z -= u;
00084     b0 -= 1.0e0;
00085     apb = a0+b0;
00086     if(apb > 1.0e0) goto S100;
00087     t = 1.0e0+gam1(&apb);
00088     goto S110;
00089 S100:
00090     u = a0+b0-1.e0;
00091     t = (1.0e0+gam1(&u))/apb;
00092 S110:
00093     brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
00094     return brcmp1;
00095 S120:
00096 /*
00097                    ALGORITHM FOR B0 .GE. 8
00098 */
00099     u = gamln1(&a0)+algdiv(&a0,&b0);
00100     T3 = z-u;
00101     brcmp1 = a0*esum(mu,&T3);
00102     return brcmp1;
00103 S130:
00104 /*
00105 -----------------------------------------------------------------------
00106               PROCEDURE FOR A .GE. 8 AND B .GE. 8
00107 -----------------------------------------------------------------------
00108 */
00109     if(*a > *b) goto S140;
00110     h = *a/ *b;
00111     x0 = h/(1.0e0+h);
00112     y0 = 1.0e0/(1.0e0+h);
00113     lambda = *a-(*a+*b)**x;
00114     goto S150;
00115 S140:
00116     h = *b/ *a;
00117     x0 = 1.0e0/(1.0e0+h);
00118     y0 = h/(1.0e0+h);
00119     lambda = (*a+*b)**y-*b;
00120 S150:
00121     e = -(lambda/ *a);
00122     if(fabs(e) > 0.6e0) goto S160;
00123     u = rlog1(&e);
00124     goto S170;
00125 S160:
00126     u = e-log(*x/x0);
00127 S170:
00128     e = lambda/ *b;
00129     if(fabs(e) > 0.6e0) goto S180;
00130     v = rlog1(&e);
00131     goto S190;
00132 S180:
00133     v = e-log(*y/y0);
00134 S190:
00135     T4 = -(*a*u+*b*v);
00136     z = esum(mu,&T4);
00137     brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
00138     return brcmp1;
00139 } /* END */

00008 {
00009 static double Const = .398942280401433e0;
00010 static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
00011 static int i,n;
00012 /*
00013 -----------------
00014      CONST = 1/SQRT(2*PI)
00015 -----------------
00016 */
00017 static double T1,T2;
00018 /*
00019      ..
00020      .. Executable Statements ..
00021 */
00022     brcomp = 0.0e0;
00023     if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
00024     a0 = fifdmin1(*a,*b);
00025     if(a0 >= 8.0e0) goto S130;
00026     if(*x > 0.375e0) goto S10;
00027     lnx = log(*x);
00028     T1 = -*x;
00029     lny = alnrel(&T1);
00030     goto S30;
00031 S10:
00032     if(*y > 0.375e0) goto S20;
00033     T2 = -*y;
00034     lnx = alnrel(&T2);
00035     lny = log(*y);
00036     goto S30;
00037 S20:
00038     lnx = log(*x);
00039     lny = log(*y);
00040 S30:
00041     z = *a*lnx+*b*lny;
00042     if(a0 < 1.0e0) goto S40;
00043     z -= betaln(a,b);
00044     brcomp = exp(z);
00045     return brcomp;
00046 S40:
00047 /*
00048 -----------------------------------------------------------------------
00049               PROCEDURE FOR A .LT. 1 OR B .LT. 1
00050 -----------------------------------------------------------------------
00051 */
00052     b0 = fifdmax1(*a,*b);
00053     if(b0 >= 8.0e0) goto S120;
00054     if(b0 > 1.0e0) goto S70;
00055 /*
00056                    ALGORITHM FOR B0 .LE. 1
00057 */
00058     brcomp = exp(z);
00059     if(brcomp == 0.0e0) return brcomp;
00060     apb = *a+*b;
00061     if(apb > 1.0e0) goto S50;
00062     z = 1.0e0+gam1(&apb);
00063     goto S60;
00064 S50:
00065     u = *a+*b-1.e0;
00066     z = (1.0e0+gam1(&u))/apb;
00067 S60:
00068     c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
00069     brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
00070     return brcomp;
00071 S70:
00072 /*
00073                 ALGORITHM FOR 1 .LT. B0 .LT. 8
00074 */
00075     u = gamln1(&a0);
00076     n = b0-1.0e0;
00077     if(n < 1) goto S90;
00078     c = 1.0e0;
00079     for(i=1; i<=n; i++) {
00080         b0 -= 1.0e0;
00081         c *= (b0/(a0+b0));
00082     }
00083     u = log(c)+u;
00084 S90:
00085     z -= u;
00086     b0 -= 1.0e0;
00087     apb = a0+b0;
00088     if(apb > 1.0e0) goto S100;
00089     t = 1.0e0+gam1(&apb);
00090     goto S110;
00091 S100:
00092     u = a0+b0-1.e0;
00093     t = (1.0e0+gam1(&u))/apb;
00094 S110:
00095     brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
00096     return brcomp;
00097 S120:
00098 /*
00099                    ALGORITHM FOR B0 .GE. 8
00100 */
00101     u = gamln1(&a0)+algdiv(&a0,&b0);
00102     brcomp = a0*exp(z-u);
00103     return brcomp;
00104 S130:
00105 /*
00106 -----------------------------------------------------------------------
00107               PROCEDURE FOR A .GE. 8 AND B .GE. 8
00108 -----------------------------------------------------------------------
00109 */
00110     if(*a > *b) goto S140;
00111     h = *a/ *b;
00112     x0 = h/(1.0e0+h);
00113     y0 = 1.0e0/(1.0e0+h);
00114     lambda = *a-(*a+*b)**x;
00115     goto S150;
00116 S140:
00117     h = *b/ *a;
00118     x0 = 1.0e0/(1.0e0+h);
00119     y0 = h/(1.0e0+h);
00120     lambda = (*a+*b)**y-*b;
00121 S150:
00122     e = -(lambda/ *a);
00123     if(fabs(e) > 0.6e0) goto S160;
00124     u = rlog1(&e);
00125     goto S170;
00126 S160:
00127     u = e-log(*x/x0);
00128 S170:
00129     e = lambda/ *b;
00130     if(fabs(e) > 0.6e0) goto S180;
00131     v = rlog1(&e);
00132     goto S190;
00133 S180:
00134     v = e-log(*y/y0);
00135 S190:
00136     z = exp(-(*a*u+*b*v));
00137     brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
00138     return brcomp;
00139 } /* END */

00009 {
00010 static int K1 = 1;
00011 static int K2 = 0;
00012 static double bup,ap1,apb,d,l,r,t,w;
00013 static int i,k,kp1,mu,nm1;
00014 /*
00015      ..
00016      .. Executable Statements ..
00017 */
00018 /*
00019           OBTAIN THE SCALING FACTOR EXP(-MU) AND
00020              EXP(MU)*(X**A*Y**B/BETA(A,B))/A
00021 */
00022     apb = *a+*b;
00023     ap1 = *a+1.0e0;
00024     mu = 0;
00025     d = 1.0e0;
00026     if(*n == 1 || *a < 1.0e0) goto S10;
00027     if(apb < 1.1e0*ap1) goto S10;
00028     mu = fabs(exparg(&K1));
00029     k = exparg(&K2);
00030     if(k < mu) mu = k;
00031     t = mu;
00032     d = exp(-t);
00033 S10:
00034     bup = brcmp1(&mu,a,b,x,y)/ *a;
00035     if(*n == 1 || bup == 0.0e0) return bup;
00036     nm1 = *n-1;
00037     w = d;
00038 /*
00039           LET K BE THE INDEX OF THE MAXIMUM TERM
00040 */
00041     k = 0;
00042     if(*b <= 1.0e0) goto S50;
00043     if(*y > 1.e-4) goto S20;
00044     k = nm1;
00045     goto S30;
00046 S20:
00047     r = (*b-1.0e0)**x/ *y-*a;
00048     if(r < 1.0e0) goto S50;
00049     k = t = nm1;
00050     if(r < t) k = r;
00051 S30:
00052 /*
00053           ADD THE INCREASING TERMS OF THE SERIES
00054 */
00055     for(i=1; i<=k; i++) {
00056         l = i-1;
00057         d = (apb+l)/(ap1+l)**x*d;
00058         w += d;
00059     }
00060     if(k == nm1) goto S70;
00061 S50:
00062 /*
00063           ADD THE REMAINING TERMS OF THE SERIES
00064 */
00065     kp1 = k+1;
00066     for(i=kp1; i<=nm1; i++) {
00067         l = i-1;
00068         d = (apb+l)/(ap1+l)**x*d;
00069         w += d;
00070         if(d <= *eps*w) goto S70;
00071     }
00072 S70:
00073 /*
00074                TERMINATE THE PROCEDURE
00075 */
00076     bup *= w;
00077     return bup;
00078 } /* END */

00025                           : 1..4
00026                iwhich = 1 : Calculate P and Q from X,Y,A and B
00027                iwhich = 2 : Calculate X and Y from P,Q,A and B
00028                iwhich = 3 : Calculate A from P,Q,X,Y and B
00029                iwhich = 4 : Calculate B from P,Q,X,Y and A
00030 
00031      P <--> The integral from 0 to X of the chi-square
00032             distribution.
00033             Input range: [0, 1].
00034 
00035      Q <--> 1-P.
00036             Input range: [0, 1].
00037             P + Q = 1.0.
00038 
00039      X <--> Upper limit of integration of beta density.
00040             Input range: [0,1].
00041             Search range: [0,1]
00042 
00043      Y <--> 1-X.
00044             Input range: [0,1].
00045             Search range: [0,1]
00046             X + Y = 1.0.
00047 
00048      A <--> The first parameter of the beta density.
00049             Input range: (0, +infinity).
00050             Search range: [1D-300,1D300]
00051 
00052      B <--> The second parameter of the beta density.
00053             Input range: (0, +infinity).
00054             Search range: [1D-300,1D300]
00055 
00056      STATUS <-- 0 if calculation completed correctly
00057                -I if input parameter number I is out of range
00058                 1 if answer appears to be lower than lowest
00059                   search bound
00060                 2 if answer appears to be higher than greatest
00061                   search bound
00062                 3 if P + Q .ne. 1
00063                 4 if X + Y .ne. 1
00064 
00065      BOUND <-- Undefined if STATUS is 0
00066 
00067                Bound exceeded by parameter number I if STATUS
00068                is negative.
00069 
00070                Lower search bound if STATUS is 1.
00071 
00072                Upper search bound if STATUS is 2.
00073 
00074 
00075                               Method
00076 
00077 
00078      Cumulative distribution function  (P)  is calculated directly by
00079      code associated with the following reference.
00080 
00081      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00082      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00083      Trans. Math.  Softw. 18 (1993), 360-373.
00084 
00085      Computation of other parameters involve a seach for a value that
00086      produces  the desired  value  of P.   The search relies  on  the
00087      monotinicity of P with the other parameter.
00088 
00089 
00090                               Note
00091 
00092 
00093      The beta density is proportional to
00094                t^(A-1) * (1-t)^(B-1)
00095 
00096 **********************************************************************/
00097 {
00098 #define tol (1.0e-8)
00099 #define atol (1.0e-50)
00100 #define zero (1.0e-300)
00101 #define inf 1.0e300
00102 #define one 1.0e0
00103 static int K1 = 1;
00104 static double K2 = 0.0e0;
00105 static double K3 = 1.0e0;
00106 static double K8 = 0.5e0;
00107 static double K9 = 5.0e0;
00108 static double fx,xhi,xlo,cum,ccum,xy,pq;
00109 static unsigned long qhi,qleft,qporq;
00110 static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
00111 /*
00112      ..
00113      .. Executable Statements ..
00114 */
00115 /*
00116      Check arguments
00117 */
00118     if(!(*which < 1 || *which > 4)) goto S30;
00119     if(!(*which < 1)) goto S10;
00120     *bound = 1.0e0;
00121     goto S20;
00122 S10:
00123     *bound = 4.0e0;
00124 S20:
00125     *status = -1;
00126     return;
00127 S30:
00128     if(*which == 1) goto S70;
00129 /*
00130      P
00131 */
00132     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00133     if(!(*p < 0.0e0)) goto S40;
00134     *bound = 0.0e0;
00135     goto S50;
00136 S40:
00137     *bound = 1.0e0;
00138 S50:
00139     *status = -2;
00140     return;
00141 S70:
00142 S60:
00143     if(*which == 1) goto S110;
00144 /*
00145      Q
00146 */
00147     if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
00148     if(!(*q < 0.0e0)) goto S80;
00149     *bound = 0.0e0;
00150     goto S90;
00151 S80:
00152     *bound = 1.0e0;
00153 S90:
00154     *status = -3;
00155     return;
00156 S110:
00157 S100:
00158     if(*which == 2) goto S150;
00159 /*
00160      X
00161 */
00162     if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
00163     if(!(*x < 0.0e0)) goto S120;
00164     *bound = 0.0e0;
00165     goto S130;
00166 S120:
00167     *bound = 1.0e0;
00168 S130:
00169     *status = -4;
00170     return;
00171 S150:
00172 S140:
00173     if(*which == 2) goto S190;
00174 /*
00175      Y
00176 */
00177     if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
00178     if(!(*y < 0.0e0)) goto S160;
00179     *bound = 0.0e0;
00180     goto S170;
00181 S160:
00182     *bound = 1.0e0;
00183 S170:
00184     *status = -5;
00185     return;
00186 S190:
00187 S180:
00188     if(*which == 3) goto S210;
00189 /*
00190      A
00191 */
00192     if(!(*a <= 0.0e0)) goto S200;
00193     *bound = 0.0e0;
00194     *status = -6;
00195     return;
00196 S210:
00197 S200:
00198     if(*which == 4) goto S230;
00199 /*
00200      B
00201 */
00202     if(!(*b <= 0.0e0)) goto S220;
00203     *bound = 0.0e0;
00204     *status = -7;
00205     return;
00206 S230:
00207 S220:
00208     if(*which == 1) goto S270;
00209 /*
00210      P + Q
00211 */
00212     pq = *p+*q;
00213     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
00214     if(!(pq < 0.0e0)) goto S240;
00215     *bound = 0.0e0;
00216     goto S250;
00217 S240:
00218     *bound = 1.0e0;
00219 S250:
00220     *status = 3;
00221     return;
00222 S270:
00223 S260:
00224     if(*which == 2) goto S310;
00225 /*
00226      X + Y
00227 */
00228     xy = *x+*y;
00229     if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
00230     if(!(xy < 0.0e0)) goto S280;
00231     *bound = 0.0e0;
00232     goto S290;
00233 S280:
00234     *bound = 1.0e0;
00235 S290:
00236     *status = 4;
00237     return;
00238 S310:
00239 S300:
00240     if(!(*which == 1)) qporq = *p <= *q;
00241 /*
00242      Select the minimum of P or Q
00243      Calculate ANSWERS
00244 */
00245     if(1 == *which) {
00246 /*
00247      Calculating P and Q
00248 */
00249         cumbet(x,y,a,b,p,q);
00250         *status = 0;
00251     }
00252     else if(2 == *which) {
00253 /*
00254      Calculating X and Y
00255 */
00256         T4 = atol;
00257         T5 = tol;
00258         dstzr(&K2,&K3,&T4,&T5);
00259         if(!qporq) goto S340;
00260         *status = 0;
00261         dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
00262         *y = one-*x;
00263 S320:
00264         if(!(*status == 1)) goto S330;
00265         cumbet(x,y,a,b,&cum,&ccum);
00266         fx = cum-*p;
00267         dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
00268         *y = one-*x;
00269         goto S320;
00270 S330:
00271         goto S370;
00272 S340:
00273         *status = 0;
00274         dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
00275         *x = one-*y;
00276 S350:
00277         if(!(*status == 1)) goto S360;
00278         cumbet(x,y,a,b,&cum,&ccum);
00279         fx = ccum-*q;
00280         dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
00281         *x = one-*y;
00282         goto S350;
00283 S370:
00284 S360:
00285         if(!(*status == -1)) goto S400;
00286         if(!qleft) goto S380;
00287         *status = 1;
00288         *bound = 0.0e0;
00289         goto S390;
00290 S380:
00291         *status = 2;
00292         *bound = 1.0e0;
00293 S400:
00294 S390:
00295         ;
00296     }
00297     else if(3 == *which) {
00298 /*
00299      Computing A
00300 */
00301         *a = 5.0e0;
00302         T6 = zero;
00303         T7 = inf;
00304         T10 = atol;
00305         T11 = tol;
00306         dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
00307         *status = 0;
00308         dinvr(status,a,&fx,&qleft,&qhi);
00309 S410:
00310         if(!(*status == 1)) goto S440;
00311         cumbet(x,y,a,b,&cum,&ccum);
00312         if(!qporq) goto S420;
00313         fx = cum-*p;
00314         goto S430;
00315 S420:
00316         fx = ccum-*q;
00317 S430:
00318         dinvr(status,a,&fx,&qleft,&qhi);
00319         goto S410;
00320 S440:
00321         if(!(*status == -1)) goto S470;
00322         if(!qleft) goto S450;
00323         *status = 1;
00324         *bound = zero;
00325         goto S460;
00326 S450:
00327         *status = 2;
00328         *bound = inf;
00329 S470:
00330 S460:
00331         ;
00332     }
00333     else if(4 == *which) {
00334 /*
00335      Computing B
00336 */
00337         *b = 5.0e0;
00338         T12 = zero;
00339         T13 = inf;
00340         T14 = atol;
00341         T15 = tol;
00342         dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
00343         *status = 0;
00344         dinvr(status,b,&fx,&qleft,&qhi);
00345 S480:
00346         if(!(*status == 1)) goto S510;
00347         cumbet(x,y,a,b,&cum,&ccum);
00348         if(!qporq) goto S490;
00349         fx = cum-*p;
00350         goto S500;
00351 S490:
00352         fx = ccum-*q;
00353 S500:
00354         dinvr(status,b,&fx,&qleft,&qhi);
00355         goto S480;
00356 S510:
00357         if(!(*status == -1)) goto S540;
00358         if(!qleft) goto S520;
00359         *status = 1;
00360         *bound = zero;
00361         goto S530;
00362 S520:
00363         *status = 2;
00364         *bound = inf;
00365 S530:
00366         ;
00367     }
00368 S540:
00369     return;
00370 #undef tol
00371 #undef atol
00372 #undef zero
00373 #undef inf
00374 #undef one
00375 } /* END */

00025                           : 1..4
00026                iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
00027                iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
00028                iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
00029                iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
00030 
00031      P <--> The cumulation from 0 to S of the binomial distribution.
00032             (Probablility of S or fewer successes in XN trials each
00033             with probability of success PR.)
00034             Input range: [0,1].
00035 
00036      Q <--> 1-P.
00037             Input range: [0, 1].
00038             P + Q = 1.0.
00039 
00040      S <--> The number of successes observed.
00041             Input range: [0, XN]
00042             Search range: [0, XN]
00043 
00044      XN  <--> The number of binomial trials.
00045               Input range: (0, +infinity).
00046               Search range: [1E-300, 1E300]
00047 
00048      PR  <--> The probability of success in each binomial trial.
00049               Input range: [0,1].
00050               Search range: [0,1]
00051 
00052      OMPR  <--> 1-PR
00053               Input range: [0,1].
00054               Search range: [0,1]
00055               PR + OMPR = 1.0
00056 
00057      STATUS <-- 0 if calculation completed correctly
00058                -I if input parameter number I is out of range
00059                 1 if answer appears to be lower than lowest
00060                   search bound
00061                 2 if answer appears to be higher than greatest
00062                   search bound
00063                 3 if P + Q .ne. 1
00064                 4 if PR + OMPR .ne. 1
00065 
00066      BOUND <-- Undefined if STATUS is 0
00067 
00068                Bound exceeded by parameter number I if STATUS
00069                is negative.
00070 
00071                Lower search bound if STATUS is 1.
00072 
00073                Upper search bound if STATUS is 2.
00074 
00075 
00076                               Method
00077 
00078 
00079      Formula  26.5.24    of   Abramowitz  and    Stegun,  Handbook   of
00080      Mathematical   Functions (1966) is   used  to reduce the  binomial
00081      distribution  to  the  cumulative incomplete    beta distribution.
00082 
00083      Computation of other parameters involve a seach for a value that
00084      produces  the desired  value  of P.   The search relies  on  the
00085      monotinicity of P with the other parameter.
00086 
00087 
00088 **********************************************************************/
00089 {
00090 #define atol (1.0e-50)
00091 #define tol (1.0e-8)
00092 #define zero (1.0e-300)
00093 #define inf 1.0e300
00094 #define one 1.0e0
00095 static int K1 = 1;
00096 static double K2 = 0.0e0;
00097 static double K3 = 0.5e0;
00098 static double K4 = 5.0e0;
00099 static double K11 = 1.0e0;
00100 static double fx,xhi,xlo,cum,ccum,pq,prompr;
00101 static unsigned long qhi,qleft,qporq;
00102 static double T5,T6,T7,T8,T9,T10,T12,T13;
00103 /*
00104      ..
00105      .. Executable Statements ..
00106 */
00107 /*
00108      Check arguments
00109 */
00110     if(!(*which < 1 && *which > 4)) goto S30;
00111     if(!(*which < 1)) goto S10;
00112     *bound = 1.0e0;
00113     goto S20;
00114 S10:
00115     *bound = 4.0e0;
00116 S20:
00117     *status = -1;
00118     return;
00119 S30:
00120     if(*which == 1) goto S70;
00121 /*
00122      P
00123 */
00124     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00125     if(!(*p < 0.0e0)) goto S40;
00126     *bound = 0.0e0;
00127     goto S50;
00128 S40:
00129     *bound = 1.0e0;
00130 S50:
00131     *status = -2;
00132     return;
00133 S70:
00134 S60:
00135     if(*which == 1) goto S110;
00136 /*
00137      Q
00138 */
00139     if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
00140     if(!(*q < 0.0e0)) goto S80;
00141     *bound = 0.0e0;
00142     goto S90;
00143 S80:
00144     *bound = 1.0e0;
00145 S90:
00146     *status = -3;
00147     return;
00148 S110:
00149 S100:
00150     if(*which == 3) goto S130;
00151 /*
00152      XN
00153 */
00154     if(!(*xn <= 0.0e0)) goto S120;
00155     *bound = 0.0e0;
00156     *status = -5;
00157     return;
00158 S130:
00159 S120:
00160     if(*which == 2) goto S170;
00161 /*
00162      S
00163 */
00164     if(!(*s < 0.0e0 || *which != 3 && *s > *xn)) goto S160;
00165     if(!(*s < 0.0e0)) goto S140;
00166     *bound = 0.0e0;
00167     goto S150;
00168 S140:
00169     *bound = *xn;
00170 S150:
00171     *status = -4;
00172     return;
00173 S170:
00174 S160:
00175     if(*which == 4) goto S210;
00176 /*
00177      PR
00178 */
00179     if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S200;
00180     if(!(*pr < 0.0e0)) goto S180;
00181     *bound = 0.0e0;
00182     goto S190;
00183 S180:
00184     *bound = 1.0e0;
00185 S190:
00186     *status = -6;
00187     return;
00188 S210:
00189 S200:
00190     if(*which == 4) goto S250;
00191 /*
00192      OMPR
00193 */
00194     if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S240;
00195     if(!(*ompr < 0.0e0)) goto S220;
00196     *bound = 0.0e0;
00197     goto S230;
00198 S220:
00199     *bound = 1.0e0;
00200 S230:
00201     *status = -7;
00202     return;
00203 S250:
00204 S240:
00205     if(*which == 1) goto S290;
00206 /*
00207      P + Q
00208 */
00209     pq = *p+*q;
00210     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S280;
00211     if(!(pq < 0.0e0)) goto S260;
00212     *bound = 0.0e0;
00213     goto S270;
00214 S260:
00215     *bound = 1.0e0;
00216 S270:
00217     *status = 3;
00218     return;
00219 S290:
00220 S280:
00221     if(*which == 4) goto S330;
00222 /*
00223      PR + OMPR
00224 */
00225     prompr = *pr+*ompr;
00226     if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S320;
00227     if(!(prompr < 0.0e0)) goto S300;
00228     *bound = 0.0e0;
00229     goto S310;
00230 S300:
00231     *bound = 1.0e0;
00232 S310:
00233     *status = 4;
00234     return;
00235 S330:
00236 S320:
00237     if(!(*which == 1)) qporq = *p <= *q;
00238 /*
00239      Select the minimum of P or Q
00240      Calculate ANSWERS
00241 */
00242     if(1 == *which) {
00243 /*
00244      Calculating P
00245 */
00246         cumbin(s,xn,pr,ompr,p,q);
00247         *status = 0;
00248     }
00249     else if(2 == *which) {
00250 /*
00251      Calculating S
00252 */
00253         *s = 5.0e0;
00254         T5 = atol;
00255         T6 = tol;
00256         dstinv(&K2,xn,&K3,&K3,&K4,&T5,&T6);
00257         *status = 0;
00258         dinvr(status,s,&fx,&qleft,&qhi);
00259 S340:
00260         if(!(*status == 1)) goto S370;
00261         cumbin(s,xn,pr,ompr,&cum,&ccum);
00262         if(!qporq) goto S350;
00263         fx = cum-*p;
00264         goto S360;
00265 S350:
00266         fx = ccum-*q;
00267 S360:
00268         dinvr(status,s,&fx,&qleft,&qhi);
00269         goto S340;
00270 S370:
00271         if(!(*status == -1)) goto S400;
00272         if(!qleft) goto S380;
00273         *status = 1;
00274         *bound = 0.0e0;
00275         goto S390;
00276 S380:
00277         *status = 2;
00278         *bound = *xn;
00279 S400:
00280 S390:
00281         ;
00282     }
00283     else if(3 == *which) {
00284 /*
00285      Calculating XN
00286 */
00287         *xn = 5.0e0;
00288         T7 = zero;
00289         T8 = inf;
00290         T9 = atol;
00291         T10 = tol;
00292         dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
00293         *status = 0;
00294         dinvr(status,xn,&fx,&qleft,&qhi);
00295 S410:
00296         if(!(*status == 1)) goto S440;
00297         cumbin(s,xn,pr,ompr,&cum,&ccum);
00298         if(!qporq) goto S420;
00299         fx = cum-*p;
00300         goto S430;
00301 S420:
00302         fx = ccum-*q;
00303 S430:
00304         dinvr(status,xn,&fx,&qleft,&qhi);
00305         goto S410;
00306 S440:
00307         if(!(*status == -1)) goto S470;
00308         if(!qleft) goto S450;
00309         *status = 1;
00310         *bound = zero;
00311         goto S460;
00312 S450:
00313         *status = 2;
00314         *bound = inf;
00315 S470:
00316 S460:
00317         ;
00318     }
00319     else if(4 == *which) {
00320 /*
00321      Calculating PR and OMPR
00322 */
00323         T12 = atol;
00324         T13 = tol;
00325         dstzr(&K2,&K11,&T12,&T13);
00326         if(!qporq) goto S500;
00327         *status = 0;
00328         dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
00329         *ompr = one-*pr;
00330 S480:
00331         if(!(*status == 1)) goto S490;
00332         cumbin(s,xn,pr,ompr,&cum,&ccum);
00333         fx = cum-*p;
00334         dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
00335         *ompr = one-*pr;
00336         goto S480;
00337 S490:
00338         goto S530;
00339 S500:
00340         *status = 0;
00341         dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
00342         *pr = one-*ompr;
00343 S510:
00344         if(!(*status == 1)) goto S520;
00345         cumbin(s,xn,pr,ompr,&cum,&ccum);
00346         fx = ccum-*q;
00347         dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
00348         *pr = one-*ompr;
00349         goto S510;
00350 S530:
00351 S520:
00352         if(!(*status == -1)) goto S560;
00353         if(!qleft) goto S540;
00354         *status = 1;
00355         *bound = 0.0e0;
00356         goto S550;
00357 S540:
00358         *status = 2;
00359         *bound = 1.0e0;
00360 S550:
00361         ;
00362     }
00363 S560:
00364     return;
00365 #undef atol
00366 #undef tol
00367 #undef zero
00368 #undef inf
00369 #undef one
00370 } /* END */

00025                           : 1..3
00026                iwhich = 1 : Calculate P and Q from X and DF
00027                iwhich = 2 : Calculate X from P,Q and DF
00028                iwhich = 3 : Calculate DF from P,Q and X
00029 
00030      P <--> The integral from 0 to X of the chi-square
00031             distribution.
00032             Input range: [0, 1].
00033 
00034      Q <--> 1-P.
00035             Input range: (0, 1].
00036             P + Q = 1.0.
00037 
00038      X <--> Upper limit of integration of the non-central
00039             chi-square distribution.
00040             Input range: [0, +infinity).
00041             Search range: [0,1E300]
00042 
00043      DF <--> Degrees of freedom of the
00044              chi-square distribution.
00045              Input range: (0, +infinity).
00046              Search range: [ 1E-300, 1E300]
00047 
00048      STATUS <-- 0 if calculation completed correctly
00049                -I if input parameter number I is out of range
00050                 1 if answer appears to be lower than lowest
00051                   search bound
00052                 2 if answer appears to be higher than greatest
00053                   search bound
00054                 3 if P + Q .ne. 1
00055                10 indicates error returned from cumgam.  See
00056                   references in cdfgam
00057 
00058      BOUND <-- Undefined if STATUS is 0
00059 
00060                Bound exceeded by parameter number I if STATUS
00061                is negative.
00062 
00063                Lower search bound if STATUS is 1.
00064 
00065                Upper search bound if STATUS is 2.
00066 
00067 
00068                               Method
00069 
00070 
00071      Formula    26.4.19   of Abramowitz  and     Stegun, Handbook  of
00072      Mathematical Functions   (1966) is used   to reduce the chisqure
00073      distribution to the incomplete distribution.
00074 
00075      Computation of other parameters involve a seach for a value that
00076      produces  the desired  value  of P.   The search relies  on  the
00077      monotinicity of P with the other parameter.
00078 
00079 **********************************************************************/
00080 {
00081 #define tol (1.0e-8)
00082 #define atol (1.0e-50)
00083 #define zero (1.0e-300)
00084 #define inf 1.0e300
00085 static int K1 = 1;
00086 static double K2 = 0.0e0;
00087 static double K4 = 0.5e0;
00088 static double K5 = 5.0e0;
00089 static double fx,cum,ccum,pq,porq;
00090 static unsigned long qhi,qleft,qporq;
00091 static double T3,T6,T7,T8,T9,T10,T11;
00092 /*
00093      ..
00094      .. Executable Statements ..
00095 */
00096 /*
00097      Check arguments
00098 */
00099     if(!(*which < 1 || *which > 3)) goto S30;
00100     if(!(*which < 1)) goto S10;
00101     *bound = 1.0e0;
00102     goto S20;
00103 S10:
00104     *bound = 3.0e0;
00105 S20:
00106     *status = -1;
00107     return;
00108 S30:
00109     if(*which == 1) goto S70;
00110 /*
00111      P
00112 */
00113     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00114     if(!(*p < 0.0e0)) goto S40;
00115     *bound = 0.0e0;
00116     goto S50;
00117 S40:
00118     *bound = 1.0e0;
00119 S50:
00120     *status = -2;
00121     return;
00122 S70:
00123 S60:
00124     if(*which == 1) goto S110;
00125 /*
00126      Q
00127 */
00128     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00129     if(!(*q <= 0.0e0)) goto S80;
00130     *bound = 0.0e0;
00131     goto S90;
00132 S80:
00133     *bound = 1.0e0;
00134 S90:
00135     *status = -3;
00136     return;
00137 S110:
00138 S100:
00139     if(*which == 2) goto S130;
00140 /*
00141      X
00142 */
00143     if(!(*x < 0.0e0)) goto S120;
00144     *bound = 0.0e0;
00145     *status = -4;
00146     return;
00147 S130:
00148 S120:
00149     if(*which == 3) goto S150;
00150 /*
00151      DF
00152 */
00153     if(!(*df <= 0.0e0)) goto S140;
00154     *bound = 0.0e0;
00155     *status = -5;
00156     return;
00157 S150:
00158 S140:
00159     if(*which == 1) goto S190;
00160 /*
00161      P + Q
00162 */
00163     pq = *p+*q;
00164     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
00165     if(!(pq < 0.0e0)) goto S160;
00166     *bound = 0.0e0;
00167     goto S170;
00168 S160:
00169     *bound = 1.0e0;
00170 S170:
00171     *status = 3;
00172     return;
00173 S190:
00174 S180:
00175     if(*which == 1) goto S220;
00176 /*
00177      Select the minimum of P or Q
00178 */
00179     qporq = *p <= *q;
00180     if(!qporq) goto S200;
00181     porq = *p;
00182     goto S210;
00183 S200:
00184     porq = *q;
00185 S220:
00186 S210:
00187 /*
00188      Calculate ANSWERS
00189 */
00190     if(1 == *which) {
00191 /*
00192      Calculating P and Q
00193 */
00194         *status = 0;
00195         cumchi(x,df,p,q);
00196         if(porq > 1.5e0) {
00197             *status = 10;
00198             return;
00199         }
00200     }
00201     else if(2 == *which) {
00202 /*
00203      Calculating X
00204 */
00205         *x = 5.0e0;
00206         T3 = inf;
00207         T6 = atol;
00208         T7 = tol;
00209         dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00210         *status = 0;
00211         dinvr(status,x,&fx,&qleft,&qhi);
00212 S230:
00213         if(!(*status == 1)) goto S270;
00214         cumchi(x,df,&cum,&ccum);
00215         if(!qporq) goto S240;
00216         fx = cum-*p;
00217         goto S250;
00218 S240:
00219         fx = ccum-*q;
00220 S250:
00221         if(!(fx+porq > 1.5e0)) goto S260;
00222         *status = 10;
00223         return;
00224 S260:
00225         dinvr(status,x,&fx,&qleft,&qhi);
00226         goto S230;
00227 S270:
00228         if(!(*status == -1)) goto S300;
00229         if(!qleft) goto S280;
00230         *status = 1;
00231         *bound = 0.0e0;
00232         goto S290;
00233 S280:
00234         *status = 2;
00235         *bound = inf;
00236 S300:
00237 S290:
00238         ;
00239     }
00240     else if(3 == *which) {
00241 /*
00242      Calculating DF
00243 */
00244         *df = 5.0e0;
00245         T8 = zero;
00246         T9 = inf;
00247         T10 = atol;
00248         T11 = tol;
00249         dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
00250         *status = 0;
00251         dinvr(status,df,&fx,&qleft,&qhi);
00252 S310:
00253         if(!(*status == 1)) goto S350;
00254         cumchi(x,df,&cum,&ccum);
00255         if(!qporq) goto S320;
00256         fx = cum-*p;
00257         goto S330;
00258 S320:
00259         fx = ccum-*q;
00260 S330:
00261         if(!(fx+porq > 1.5e0)) goto S340;
00262         *status = 10;
00263         return;
00264 S340:
00265         dinvr(status,df,&fx,&qleft,&qhi);
00266         goto S310;
00267 S350:
00268         if(!(*status == -1)) goto S380;
00269         if(!qleft) goto S360;
00270         *status = 1;
00271         *bound = zero;
00272         goto S370;
00273 S360:
00274         *status = 2;
00275         *bound = inf;
00276 S370:
00277         ;
00278     }
00279 S380:
00280     return;
00281 #undef tol
00282 #undef atol
00283 #undef zero
00284 #undef inf
00285 } /* END */

00025                           : 1..4
00026                iwhich = 1 : Calculate P and Q from X and DF
00027                iwhich = 2 : Calculate X from P,DF and PNONC
00028                iwhich = 3 : Calculate DF from P,X and PNONC
00029                iwhich = 3 : Calculate PNONC from P,X and DF
00030 
00031      P <--> The integral from 0 to X of the non-central chi-square
00032             distribution.
00033             Input range: [0, 1-1E-16).
00034 
00035      Q <--> 1-P.
00036             Q is not used by this subroutine and is only included
00037             for similarity with other cdf* routines.
00038 
00039      X <--> Upper limit of integration of the non-central
00040             chi-square distribution.
00041             Input range: [0, +infinity).
00042             Search range: [0,1E300]
00043 
00044      DF <--> Degrees of freedom of the non-central
00045              chi-square distribution.
00046              Input range: (0, +infinity).
00047              Search range: [ 1E-300, 1E300]
00048 
00049      PNONC <--> Non-centrality parameter of the non-central
00050                 chi-square distribution.
00051                 Input range: [0, +infinity).
00052                 Search range: [0,1E4]
00053 
00054      STATUS <-- 0 if calculation completed correctly
00055                -I if input parameter number I is out of range
00056                 1 if answer appears to be lower than lowest
00057                   search bound
00058                 2 if answer appears to be higher than greatest
00059                   search bound
00060 
00061      BOUND <-- Undefined if STATUS is 0
00062 
00063                Bound exceeded by parameter number I if STATUS
00064                is negative.
00065 
00066                Lower search bound if STATUS is 1.
00067 
00068                Upper search bound if STATUS is 2.
00069 
00070 
00071                               Method
00072 
00073 
00074      Formula  26.4.25   of   Abramowitz   and   Stegun,  Handbook  of
00075      Mathematical  Functions (1966) is used to compute the cumulative
00076      distribution function.
00077 
00078      Computation of other parameters involve a seach for a value that
00079      produces  the desired  value  of P.   The search relies  on  the
00080      monotinicity of P with the other parameter.
00081 
00082 
00083                             WARNING
00084 
00085      The computation time  required for this  routine is proportional
00086      to the noncentrality  parameter  (PNONC).  Very large  values of
00087      this parameter can consume immense  computer resources.  This is
00088      why the search range is bounded by 10,000.
00089 
00090 **********************************************************************/
00091 {
00092 #define tent4 1.0e4
00093 #define tol (1.0e-8)
00094 #define atol (1.0e-50)
00095 #define zero (1.0e-300)
00096 #define one (1.0e0-1.0e-16)
00097 #define inf 1.0e300
00098 static double K1 = 0.0e0;
00099 static double K3 = 0.5e0;
00100 static double K4 = 5.0e0;
00101 static double fx,cum,ccum;
00102 static unsigned long qhi,qleft;
00103 static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;
00104 /*
00105      ..
00106      .. Executable Statements ..
00107 */
00108 /*
00109      Check arguments
00110 */
00111     if(!(*which < 1 || *which > 4)) goto S30;
00112     if(!(*which < 1)) goto S10;
00113     *bound = 1.0e0;
00114     goto S20;
00115 S10:
00116     *bound = 4.0e0;
00117 S20:
00118     *status = -1;
00119     return;
00120 S30:
00121     if(*which == 1) goto S70;
00122 /*
00123      P
00124 */
00125     if(!(*p < 0.0e0 || *p > one)) goto S60;
00126     if(!(*p < 0.0e0)) goto S40;
00127     *bound = 0.0e0;
00128     goto S50;
00129 S40:
00130     *bound = one;
00131 S50:
00132     *status = -2;
00133     return;
00134 S70:
00135 S60:
00136     if(*which == 2) goto S90;
00137 /*
00138      X
00139 */
00140     if(!(*x < 0.0e0)) goto S80;
00141     *bound = 0.0e0;
00142     *status = -4;
00143     return;
00144 S90:
00145 S80:
00146     if(*which == 3) goto S110;
00147 /*
00148      DF
00149 */
00150     if(!(*df <= 0.0e0)) goto S100;
00151     *bound = 0.0e0;
00152     *status = -5;
00153     return;
00154 S110:
00155 S100:
00156     if(*which == 4) goto S130;
00157 /*
00158      PNONC
00159 */
00160     if(!(*pnonc < 0.0e0)) goto S120;
00161     *bound = 0.0e0;
00162     *status = -6;
00163     return;
00164 S130:
00165 S120:
00166 /*
00167      Calculate ANSWERS
00168 */
00169     if(1 == *which) {
00170 /*
00171      Calculating P and Q
00172 */
00173         cumchn(x,df,pnonc,p,q);
00174         *status = 0;
00175     }
00176     else if(2 == *which) {
00177 /*
00178      Calculating X
00179 */
00180         *x = 5.0e0;
00181         T2 = inf;
00182         T5 = atol;
00183         T6 = tol;
00184         dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
00185         *status = 0;
00186         dinvr(status,x,&fx,&qleft,&qhi);
00187 S140:
00188         if(!(*status == 1)) goto S150;
00189         cumchn(x,df,pnonc,&cum,&ccum);
00190         fx = cum-*p;
00191         dinvr(status,x,&fx,&qleft,&qhi);
00192         goto S140;
00193 S150:
00194         if(!(*status == -1)) goto S180;
00195         if(!qleft) goto S160;
00196         *status = 1;
00197         *bound = 0.0e0;
00198         goto S170;
00199 S160:
00200         *status = 2;
00201         *bound = inf;
00202 S180:
00203 S170:
00204         ;
00205     }
00206     else if(3 == *which) {
00207 /*
00208      Calculating DF
00209 */
00210         *df = 5.0e0;
00211         T7 = zero;
00212         T8 = inf;
00213         T9 = atol;
00214         T10 = tol;
00215         dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
00216         *status = 0;
00217         dinvr(status,df,&fx,&qleft,&qhi);
00218 S190:
00219         if(!(*status == 1)) goto S200;
00220         cumchn(x,df,pnonc,&cum,&ccum);
00221         fx = cum-*p;
00222         dinvr(status,df,&fx,&qleft,&qhi);
00223         goto S190;
00224 S200:
00225         if(!(*status == -1)) goto S230;
00226         if(!qleft) goto S210;
00227         *status = 1;
00228         *bound = zero;
00229         goto S220;
00230 S210:
00231         *status = 2;
00232         *bound = inf;
00233 S230:
00234 S220:
00235         ;
00236     }
00237     else if(4 == *which) {
00238 /*
00239      Calculating PNONC
00240 */
00241         *pnonc = 5.0e0;
00242         T11 = tent4;
00243         T12 = atol;
00244         T13 = tol;
00245         dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);
00246         *status = 0;
00247         dinvr(status,pnonc,&fx,&qleft,&qhi);
00248 S240:
00249         if(!(*status == 1)) goto S250;
00250         cumchn(x,df,pnonc,&cum,&ccum);
00251         fx = cum-*p;
00252         dinvr(status,pnonc,&fx,&qleft,&qhi);
00253         goto S240;
00254 S250:
00255         if(!(*status == -1)) goto S280;
00256         if(!qleft) goto S260;
00257         *status = 1;
00258         *bound = zero;
00259         goto S270;
00260 S260:
00261         *status = 2;
00262         *bound = tent4;
00263 S270:
00264         ;
00265     }
00266 S280:
00267     return;
00268 #undef tent4
00269 #undef tol
00270 #undef atol
00271 #undef zero
00272 #undef one
00273 #undef inf
00274 } /* END */

00025                           : 1..4
00026                iwhich = 1 : Calculate P and Q from F,DFN and DFD
00027                iwhich = 2 : Calculate F from P,Q,DFN and DFD
00028                iwhich = 3 : Calculate DFN from P,Q,F and DFD
00029                iwhich = 4 : Calculate DFD from P,Q,F and DFN
00030 
00031        P <--> The integral from 0 to F of the f-density.
00032               Input range: [0,1].
00033 
00034        Q <--> 1-P.
00035               Input range: (0, 1].
00036               P + Q = 1.0.
00037 
00038        F <--> Upper limit of integration of the f-density.
00039               Input range: [0, +infinity).
00040               Search range: [0,1E300]
00041 
00042      DFN < --> Degrees of freedom of the numerator sum of squares.
00043                Input range: (0, +infinity).
00044                Search range: [ 1E-300, 1E300]
00045 
00046      DFD < --> Degrees of freedom of the denominator sum of squares.
00047                Input range: (0, +infinity).
00048                Search range: [ 1E-300, 1E300]
00049 
00050      STATUS <-- 0 if calculation completed correctly
00051                -I if input parameter number I is out of range
00052                 1 if answer appears to be lower than lowest
00053                   search bound
00054                 2 if answer appears to be higher than greatest
00055                   search bound
00056                 3 if P + Q .ne. 1
00057 
00058      BOUND <-- Undefined if STATUS is 0
00059 
00060                Bound exceeded by parameter number I if STATUS
00061                is negative.
00062 
00063                Lower search bound if STATUS is 1.
00064 
00065                Upper search bound if STATUS is 2.
00066 
00067 
00068                               Method
00069 
00070 
00071      Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of
00072      Mathematical  Functions (1966) is used to reduce the computation
00073      of the  cumulative  distribution function for the  F  variate to
00074      that of an incomplete beta.
00075 
00076      Computation of other parameters involve a seach for a value that
00077      produces  the desired  value  of P.   The search relies  on  the
00078      monotinicity of P with the other parameter.
00079 
00080                               WARNING
00081 
00082      The value of the  cumulative  F distribution is  not necessarily
00083      monotone in  either degrees of freedom.  There  thus may  be two
00084      values  that  provide a given CDF  value.   This routine assumes
00085      monotonicity and will find an arbitrary one of the two values.
00086 
00087 **********************************************************************/
00088 {
00089 #define tol (1.0e-8)
00090 #define atol (1.0e-50)
00091 #define zero (1.0e-300)
00092 #define inf 1.0e300
00093 static int K1 = 1;
00094 static double K2 = 0.0e0;
00095 static double K4 = 0.5e0;
00096 static double K5 = 5.0e0;
00097 static double pq,fx,cum,ccum;
00098 static unsigned long qhi,qleft,qporq;
00099 static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
00100 /*
00101      ..
00102      .. Executable Statements ..
00103 */
00104 /*
00105      Check arguments
00106 */
00107     if(!(*which < 1 || *which > 4)) goto S30;
00108     if(!(*which < 1)) goto S10;
00109     *bound = 1.0e0;
00110     goto S20;
00111 S10:
00112     *bound = 4.0e0;
00113 S20:
00114     *status = -1;
00115     return;
00116 S30:
00117     if(*which == 1) goto S70;
00118 /*
00119      P
00120 */
00121     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00122     if(!(*p < 0.0e0)) goto S40;
00123     *bound = 0.0e0;
00124     goto S50;
00125 S40:
00126     *bound = 1.0e0;
00127 S50:
00128     *status = -2;
00129     return;
00130 S70:
00131 S60:
00132     if(*which == 1) goto S110;
00133 /*
00134      Q
00135 */
00136     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00137     if(!(*q <= 0.0e0)) goto S80;
00138     *bound = 0.0e0;
00139     goto S90;
00140 S80:
00141     *bound = 1.0e0;
00142 S90:
00143     *status = -3;
00144     return;
00145 S110:
00146 S100:
00147     if(*which == 2) goto S130;
00148 /*
00149      F
00150 */
00151     if(!(*f < 0.0e0)) goto S120;
00152     *bound = 0.0e0;
00153     *status = -4;
00154     return;
00155 S130:
00156 S120:
00157     if(*which == 3) goto S150;
00158 /*
00159      DFN
00160 */
00161     if(!(*dfn <= 0.0e0)) goto S140;
00162     *bound = 0.0e0;
00163     *status = -5;
00164     return;
00165 S150:
00166 S140:
00167     if(*which == 4) goto S170;
00168 /*
00169      DFD
00170 */
00171     if(!(*dfd <= 0.0e0)) goto S160;
00172     *bound = 0.0e0;
00173     *status = -6;
00174     return;
00175 S170:
00176 S160:
00177     if(*which == 1) goto S210;
00178 /*
00179      P + Q
00180 */
00181     pq = *p+*q;
00182     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
00183     if(!(pq < 0.0e0)) goto S180;
00184     *bound = 0.0e0;
00185     goto S190;
00186 S180:
00187     *bound = 1.0e0;
00188 S190:
00189     *status = 3;
00190     return;
00191 S210:
00192 S200:
00193     if(!(*which == 1)) qporq = *p <= *q;
00194 /*
00195      Select the minimum of P or Q
00196      Calculate ANSWERS
00197 */
00198     if(1 == *which) {
00199 /*
00200      Calculating P
00201 */
00202         cumf(f,dfn,dfd,p,q);
00203         *status = 0;
00204     }
00205     else if(2 == *which) {
00206 /*
00207      Calculating F
00208 */
00209         *f = 5.0e0;
00210         T3 = inf;
00211         T6 = atol;
00212         T7 = tol;
00213         dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00214         *status = 0;
00215         dinvr(status,f,&fx,&qleft,&qhi);
00216 S220:
00217         if(!(*status == 1)) goto S250;
00218         cumf(f,dfn,dfd,&cum,&ccum);
00219         if(!qporq) goto S230;
00220         fx = cum-*p;
00221         goto S240;
00222 S230:
00223         fx = ccum-*q;
00224 S240:
00225         dinvr(status,f,&fx,&qleft,&qhi);
00226         goto S220;
00227 S250:
00228         if(!(*status == -1)) goto S280;
00229         if(!qleft) goto S260;
00230         *status = 1;
00231         *bound = 0.0e0;
00232         goto S270;
00233 S260:
00234         *status = 2;
00235         *bound = inf;
00236 S280:
00237 S270:
00238         ;
00239     }
00240     else if(3 == *which) {
00241 /*
00242      Calculating DFN
00243 */
00244         *dfn = 5.0e0;
00245         T8 = zero;
00246         T9 = inf;
00247         T10 = atol;
00248         T11 = tol;
00249         dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
00250         *status = 0;
00251         dinvr(status,dfn,&fx,&qleft,&qhi);
00252 S290:
00253         if(!(*status == 1)) goto S320;
00254         cumf(f,dfn,dfd,&cum,&ccum);
00255         if(!qporq) goto S300;
00256         fx = cum-*p;
00257         goto S310;
00258 S300:
00259         fx = ccum-*q;
00260 S310:
00261         dinvr(status,dfn,&fx,&qleft,&qhi);
00262         goto S290;
00263 S320:
00264         if(!(*status == -1)) goto S350;
00265         if(!qleft) goto S330;
00266         *status = 1;
00267         *bound = zero;
00268         goto S340;
00269 S330:
00270         *status = 2;
00271         *bound = inf;
00272 S350:
00273 S340:
00274         ;
00275     }
00276     else if(4 == *which) {
00277 /*
00278      Calculating DFD
00279 */
00280         *dfd = 5.0e0;
00281         T12 = zero;
00282         T13 = inf;
00283         T14 = atol;
00284         T15 = tol;
00285         dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15);
00286         *status = 0;
00287         dinvr(status,dfd,&fx,&qleft,&qhi);
00288 S360:
00289         if(!(*status == 1)) goto S390;
00290         cumf(f,dfn,dfd,&cum,&ccum);
00291         if(!qporq) goto S370;
00292         fx = cum-*p;
00293         goto S380;
00294 S370:
00295         fx = ccum-*q;
00296 S380:
00297         dinvr(status,dfd,&fx,&qleft,&qhi);
00298         goto S360;
00299 S390:
00300         if(!(*status == -1)) goto S420;
00301         if(!qleft) goto S400;
00302         *status = 1;
00303         *bound = zero;
00304         goto S410;
00305 S400:
00306         *status = 2;
00307         *bound = inf;
00308 S410:
00309         ;
00310     }
00311 S420:
00312     return;
00313 #undef tol
00314 #undef atol
00315 #undef zero
00316 #undef inf
00317 } /* END */

00025                           : 1..5
00026                iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
00027                iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
00028                iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
00029                iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
00030                iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
00031 
00032        P <--> The integral from 0 to F of the non-central f-density.
00033               Input range: [0,1-1E-16).
00034 
00035        Q <--> 1-P.
00036               Q is not used by this subroutine and is only included
00037               for similarity with other cdf* routines.
00038 
00039        F <--> Upper limit of integration of the non-central f-density.
00040               Input range: [0, +infinity).
00041               Search range: [0,1E300]
00042 
00043      DFN < --> Degrees of freedom of the numerator sum of squares.
00044                Input range: (0, +infinity).
00045                Search range: [ 1E-300, 1E300]
00046 
00047      DFD < --> Degrees of freedom of the denominator sum of squares.
00048                Must be in range: (0, +infinity).
00049                Input range: (0, +infinity).
00050                Search range: [ 1E-300, 1E300]
00051 
00052      PNONC <-> The non-centrality parameter
00053                Input range: [0,infinity)
00054                Search range: [0,1E4]
00055 
00056      STATUS <-- 0 if calculation completed correctly
00057                -I if input parameter number I is out of range
00058                 1 if answer appears to be lower than lowest
00059                   search bound
00060                 2 if answer appears to be higher than greatest
00061                   search bound
00062                 3 if P + Q .ne. 1
00063 
00064      BOUND <-- Undefined if STATUS is 0
00065 
00066                Bound exceeded by parameter number I if STATUS
00067                is negative.
00068 
00069                Lower search bound if STATUS is 1.
00070 
00071                Upper search bound if STATUS is 2.
00072 
00073 
00074                               Method
00075 
00076 
00077      Formula  26.6.20   of   Abramowitz   and   Stegun,  Handbook  of
00078      Mathematical  Functions (1966) is used to compute the cumulative
00079      distribution function.
00080 
00081      Computation of other parameters involve a seach for a value that
00082      produces  the desired  value  of P.   The search relies  on  the
00083      monotinicity of P with the other parameter.
00084 
00085                             WARNING
00086 
00087      The computation time  required for this  routine is proportional
00088      to the noncentrality  parameter  (PNONC).  Very large  values of
00089      this parameter can consume immense  computer resources.  This is
00090      why the search range is bounded by 10,000.
00091 
00092                               WARNING
00093 
00094      The  value  of the  cumulative  noncentral F distribution is not
00095      necessarily monotone in either degrees  of freedom.  There  thus
00096      may be two values that provide a given  CDF value.  This routine
00097      assumes monotonicity  and will find  an arbitrary one of the two
00098      values.
00099 
00100 **********************************************************************/
00101 {
00102 #define tent4 1.0e4
00103 #define tol (1.0e-8)
00104 #define atol (1.0e-50)
00105 #define zero (1.0e-300)
00106 #define one (1.0e0-1.0e-16)
00107 #define inf 1.0e300
00108 static double K1 = 0.0e0;
00109 static double K3 = 0.5e0;
00110 static double K4 = 5.0e0;
00111 static double fx,cum,ccum;
00112 static unsigned long qhi,qleft;
00113 static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17;
00114 /*
00115      ..
00116      .. Executable Statements ..
00117 */
00118 /*
00119      Check arguments
00120 */
00121     if(!(*which < 1 || *which > 5)) goto S30;
00122     if(!(*which < 1)) goto S10;
00123     *bound = 1.0e0;
00124     goto S20;
00125 S10:
00126     *bound = 5.0e0;
00127 S20:
00128     *status = -1;
00129     return;
00130 S30:
00131     if(*which == 1) goto S70;
00132 /*
00133      P
00134 */
00135     if(!(*p < 0.0e0 || *p > one)) goto S60;
00136     if(!(*p < 0.0e0)) goto S40;
00137     *bound = 0.0e0;
00138     goto S50;
00139 S40:
00140     *bound = one;
00141 S50:
00142     *status = -2;
00143     return;
00144 S70:
00145 S60:
00146     if(*which == 2) goto S90;
00147 /*
00148      F
00149 */
00150     if(!(*f < 0.0e0)) goto S80;
00151     *bound = 0.0e0;
00152     *status = -4;
00153     return;
00154 S90:
00155 S80:
00156     if(*which == 3) goto S110;
00157 /*
00158      DFN
00159 */
00160     if(!(*dfn <= 0.0e0)) goto S100;
00161     *bound = 0.0e0;
00162     *status = -5;
00163     return;
00164 S110:
00165 S100:
00166     if(*which == 4) goto S130;
00167 /*
00168      DFD
00169 */
00170     if(!(*dfd <= 0.0e0)) goto S120;
00171     *bound = 0.0e0;
00172     *status = -6;
00173     return;
00174 S130:
00175 S120:
00176     if(*which == 5) goto S150;
00177 /*
00178      PHONC
00179 */
00180     if(!(*phonc < 0.0e0)) goto S140;
00181     *bound = 0.0e0;
00182     *status = -7;
00183     return;
00184 S150:
00185 S140:
00186 /*
00187      Calculate ANSWERS
00188 */
00189     if(1 == *which) {
00190 /*
00191      Calculating P
00192 */
00193         cumfnc(f,dfn,dfd,phonc,p,q);
00194         *status = 0;
00195     }
00196     else if(2 == *which) {
00197 /*
00198      Calculating F
00199 */
00200         *f = 5.0e0;
00201         T2 = inf;
00202         T5 = atol;
00203         T6 = tol;
00204         dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
00205         *status = 0;
00206         dinvr(status,f,&fx,&qleft,&qhi);
00207 S160:
00208         if(!(*status == 1)) goto S170;
00209         cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
00210         fx = cum-*p;
00211         dinvr(status,f,&fx,&qleft,&qhi);
00212         goto S160;
00213 S170:
00214         if(!(*status == -1)) goto S200;
00215         if(!qleft) goto S180;
00216         *status = 1;
00217         *bound = 0.0e0;
00218         goto S190;
00219 S180:
00220         *status = 2;
00221         *bound = inf;
00222 S200:
00223 S190:
00224         ;
00225     }
00226     else if(3 == *which) {
00227 /*
00228      Calculating DFN
00229 */
00230         *dfn = 5.0e0;
00231         T7 = zero;
00232         T8 = inf;
00233         T9 = atol;
00234         T10 = tol;
00235         dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
00236         *status = 0;
00237         dinvr(status,dfn,&fx,&qleft,&qhi);
00238 S210:
00239         if(!(*status == 1)) goto S220;
00240         cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
00241         fx = cum-*p;
00242         dinvr(status,dfn,&fx,&qleft,&qhi);
00243         goto S210;
00244 S220:
00245         if(!(*status == -1)) goto S250;
00246         if(!qleft) goto S230;
00247         *status = 1;
00248         *bound = zero;
00249         goto S240;
00250 S230:
00251         *status = 2;
00252         *bound = inf;
00253 S250:
00254 S240:
00255         ;
00256     }
00257     else if(4 == *which) {
00258 /*
00259      Calculating DFD
00260 */
00261         *dfd = 5.0e0;
00262         T11 = zero;
00263         T12 = inf;
00264         T13 = atol;
00265         T14 = tol;
00266         dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
00267         *status = 0;
00268         dinvr(status,dfd,&fx,&qleft,&qhi);
00269 S260:
00270         if(!(*status == 1)) goto S270;
00271         cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
00272         fx = cum-*p;
00273         dinvr(status,dfd,&fx,&qleft,&qhi);
00274         goto S260;
00275 S270:
00276         if(!(*status == -1)) goto S300;
00277         if(!qleft) goto S280;
00278         *status = 1;
00279         *bound = zero;
00280         goto S290;
00281 S280:
00282         *status = 2;
00283         *bound = inf;
00284 S300:
00285 S290:
00286         ;
00287     }
00288     else if(5 == *which) {
00289 /*
00290      Calculating PHONC
00291 */
00292         *phonc = 5.0e0;
00293         T15 = tent4;
00294         T16 = atol;
00295         T17 = tol;
00296         dstinv(&K1,&T15,&K3,&K3,&K4,&T16,&T17);
00297         *status = 0;
00298         dinvr(status,phonc,&fx,&qleft,&qhi);
00299 S310:
00300         if(!(*status == 1)) goto S320;
00301         cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
00302         fx = cum-*p;
00303         dinvr(status,phonc,&fx,&qleft,&qhi);
00304         goto S310;
00305 S320:
00306         if(!(*status == -1)) goto S350;
00307         if(!qleft) goto S330;
00308         *status = 1;
00309         *bound = 0.0e0;
00310         goto S340;
00311 S330:
00312         *status = 2;
00313         *bound = tent4;
00314 S340:
00315         ;
00316     }
00317 S350:
00318     return;
00319 #undef tent4
00320 #undef tol
00321 #undef atol
00322 #undef zero
00323 #undef one
00324 #undef inf
00325 } /* END */

00025                           : 1..4
00026                iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE
00027                iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE
00028                iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE
00029                iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE
00030 
00031      P <--> The integral from 0 to X of the gamma density.
00032             Input range: [0,1].
00033 
00034      Q <--> 1-P.
00035             Input range: (0, 1].
00036             P + Q = 1.0.
00037 
00038      X <--> The upper limit of integration of the gamma density.
00039             Input range: [0, +infinity).
00040             Search range: [0,1E300]
00041 
00042      SHAPE <--> The shape parameter of the gamma density.
00043                 Input range: (0, +infinity).
00044                 Search range: [1E-300,1E300]
00045 
00046      SCALE <--> The scale parameter of the gamma density.
00047                 Input range: (0, +infinity).
00048                 Search range: (1E-300,1E300]
00049 
00050      STATUS <-- 0 if calculation completed correctly
00051                -I if input parameter number I is out of range
00052                 1 if answer appears to be lower than lowest
00053                   search bound
00054                 2 if answer appears to be higher than greatest
00055                   search bound
00056                 3 if P + Q .ne. 1
00057                 10 if the gamma or inverse gamma routine cannot
00058                    compute the answer.  Usually happens only for
00059                    X and SHAPE very large (gt 1E10 or more)
00060 
00061      BOUND <-- Undefined if STATUS is 0
00062 
00063                Bound exceeded by parameter number I if STATUS
00064                is negative.
00065 
00066                Lower search bound if STATUS is 1.
00067 
00068                Upper search bound if STATUS is 2.
00069 
00070 
00071                               Method
00072 
00073 
00074      Cumulative distribution function (P) is calculated directly by
00075      the code associated with:
00076 
00077      DiDinato, A. R. and Morris, A. H. Computation of the  incomplete
00078      gamma function  ratios  and their  inverse.   ACM  Trans.  Math.
00079      Softw. 12 (1986), 377-393.
00080 
00081      Computation of other parameters involve a seach for a value that
00082      produces  the desired  value  of P.   The search relies  on  the
00083      monotinicity of P with the other parameter.
00084 
00085 
00086                               Note
00087 
00088 
00089 
00090      The gamma density is proportional to
00091        T**(SHAPE - 1) * EXP(- SCALE * T)
00092 
00093 **********************************************************************/
00094 {
00095 #define tol (1.0e-8)
00096 #define atol (1.0e-50)
00097 #define zero (1.0e-300)
00098 #define inf 1.0e300
00099 static int K1 = 1;
00100 static double K5 = 0.5e0;
00101 static double K6 = 5.0e0;
00102 static double xx,fx,xscale,cum,ccum,pq,porq;
00103 static int ierr;
00104 static unsigned long qhi,qleft,qporq;
00105 static double T2,T3,T4,T7,T8,T9;
00106 /*
00107      ..
00108      .. Executable Statements ..
00109 */
00110 /*
00111      Check arguments
00112 */
00113     if(!(*which < 1 || *which > 4)) goto S30;
00114     if(!(*which < 1)) goto S10;
00115     *bound = 1.0e0;
00116     goto S20;
00117 S10:
00118     *bound = 4.0e0;
00119 S20:
00120     *status = -1;
00121     return;
00122 S30:
00123     if(*which == 1) goto S70;
00124 /*
00125      P
00126 */
00127     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00128     if(!(*p < 0.0e0)) goto S40;
00129     *bound = 0.0e0;
00130     goto S50;
00131 S40:
00132     *bound = 1.0e0;
00133 S50:
00134     *status = -2;
00135     return;
00136 S70:
00137 S60:
00138     if(*which == 1) goto S110;
00139 /*
00140      Q
00141 */
00142     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00143     if(!(*q <= 0.0e0)) goto S80;
00144     *bound = 0.0e0;
00145     goto S90;
00146 S80:
00147     *bound = 1.0e0;
00148 S90:
00149     *status = -3;
00150     return;
00151 S110:
00152 S100:
00153     if(*which == 2) goto S130;
00154 /*
00155      X
00156 */
00157     if(!(*x < 0.0e0)) goto S120;
00158     *bound = 0.0e0;
00159     *status = -4;
00160     return;
00161 S130:
00162 S120:
00163     if(*which == 3) goto S150;
00164 /*
00165      SHAPE
00166 */
00167     if(!(*shape <= 0.0e0)) goto S140;
00168     *bound = 0.0e0;
00169     *status = -5;
00170     return;
00171 S150:
00172 S140:
00173     if(*which == 4) goto S170;
00174 /*
00175      SCALE
00176 */
00177     if(!(*scale <= 0.0e0)) goto S160;
00178     *bound = 0.0e0;
00179     *status = -6;
00180     return;
00181 S170:
00182 S160:
00183     if(*which == 1) goto S210;
00184 /*
00185      P + Q
00186 */
00187     pq = *p+*q;
00188     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
00189     if(!(pq < 0.0e0)) goto S180;
00190     *bound = 0.0e0;
00191     goto S190;
00192 S180:
00193     *bound = 1.0e0;
00194 S190:
00195     *status = 3;
00196     return;
00197 S210:
00198 S200:
00199     if(*which == 1) goto S240;
00200 /*
00201      Select the minimum of P or Q
00202 */
00203     qporq = *p <= *q;
00204     if(!qporq) goto S220;
00205     porq = *p;
00206     goto S230;
00207 S220:
00208     porq = *q;
00209 S240:
00210 S230:
00211 /*
00212      Calculate ANSWERS
00213 */
00214     if(1 == *which) {
00215 /*
00216      Calculating P
00217 */
00218         *status = 0;
00219         xscale = *x**scale;
00220         cumgam(&xscale,shape,p,q);
00221         if(porq > 1.5e0) *status = 10;
00222     }
00223     else if(2 == *which) {
00224 /*
00225      Computing X
00226 */
00227         T2 = -1.0e0;
00228         gaminv(shape,&xx,&T2,p,q,&ierr);
00229         if(ierr < 0.0e0) {
00230             *status = 10;
00231             return;
00232         }
00233         else  {
00234             *x = xx/ *scale;
00235             *status = 0;
00236         }
00237     }
00238     else if(3 == *which) {
00239 /*
00240      Computing SHAPE
00241 */
00242         *shape = 5.0e0;
00243         xscale = *x**scale;
00244         T3 = zero;
00245         T4 = inf;
00246         T7 = atol;
00247         T8 = tol;
00248         dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
00249         *status = 0;
00250         dinvr(status,shape,&fx,&qleft,&qhi);
00251 S250:
00252         if(!(*status == 1)) goto S290;
00253         cumgam(&xscale,shape,&cum,&ccum);
00254         if(!qporq) goto S260;
00255         fx = cum-*p;
00256         goto S270;
00257 S260:
00258         fx = ccum-*q;
00259 S270:
00260         if(!(qporq && cum > 1.5e0 || !qporq && ccum > 1.5e0)) goto S280;
00261         *status = 10;
00262         return;
00263 S280:
00264         dinvr(status,shape,&fx,&qleft,&qhi);
00265         goto S250;
00266 S290:
00267         if(!(*status == -1)) goto S320;
00268         if(!qleft) goto S300;
00269         *status = 1;
00270         *bound = zero;
00271         goto S310;
00272 S300:
00273         *status = 2;
00274         *bound = inf;
00275 S320:
00276 S310:
00277         ;
00278     }
00279     else if(4 == *which) {
00280 /*
00281      Computing SCALE
00282 */
00283         T9 = -1.0e0;
00284         gaminv(shape,&xx,&T9,p,q,&ierr);
00285         if(ierr < 0.0e0) {
00286             *status = 10;
00287             return;
00288         }
00289         else  {
00290             *scale = xx/ *x;
00291             *status = 0;
00292         }
00293     }
00294     return;
00295 #undef tol
00296 #undef atol
00297 #undef zero
00298 #undef inf
00299 } /* END */

00034                           : 1..4
00035                iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
00036                iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
00037                iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
00038                iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
00039 
00040      P <--> The cumulation from 0 to S of the  negative
00041             binomial distribution.
00042             Input range: [0,1].
00043 
00044      Q <--> 1-P.
00045             Input range: (0, 1].
00046             P + Q = 1.0.
00047 
00048      S <--> The upper limit of cumulation of the binomial distribution.
00049             There are F or fewer failures before the XNth success.
00050             Input range: [0, +infinity).
00051             Search range: [0, 1E300]
00052 
00053      XN  <--> The number of successes.
00054               Input range: [0, +infinity).
00055               Search range: [0, 1E300]
00056 
00057      PR  <--> The probability of success in each binomial trial.
00058               Input range: [0,1].
00059               Search range: [0,1].
00060 
00061      OMPR  <--> 1-PR
00062               Input range: [0,1].
00063               Search range: [0,1]
00064               PR + OMPR = 1.0
00065 
00066      STATUS <-- 0 if calculation completed correctly
00067                -I if input parameter number I is out of range
00068                 1 if answer appears to be lower than lowest
00069                   search bound
00070                 2 if answer appears to be higher than greatest
00071                   search bound
00072                 3 if P + Q .ne. 1
00073                 4 if PR + OMPR .ne. 1
00074 
00075      BOUND <-- Undefined if STATUS is 0
00076 
00077                Bound exceeded by parameter number I if STATUS
00078                is negative.
00079 
00080                Lower search bound if STATUS is 1.
00081 
00082                Upper search bound if STATUS is 2.
00083 
00084 
00085                               Method
00086 
00087 
00088      Formula   26.5.26   of   Abramowitz  and  Stegun,  Handbook   of
00089      Mathematical Functions (1966) is used  to  reduce calculation of
00090      the cumulative distribution  function to that of  an  incomplete
00091      beta.
00092 
00093      Computation of other parameters involve a seach for a value that
00094      produces  the desired  value  of P.   The search relies  on  the
00095      monotinicity of P with the other parameter.
00096 
00097 **********************************************************************/
00098 {
00099 #define tol (1.0e-8)
00100 #define atol (1.0e-50)
00101 #define inf 1.0e300
00102 #define one 1.0e0
00103 static int K1 = 1;
00104 static double K2 = 0.0e0;
00105 static double K4 = 0.5e0;
00106 static double K5 = 5.0e0;
00107 static double K11 = 1.0e0;
00108 static double fx,xhi,xlo,pq,prompr,cum,ccum;
00109 static unsigned long qhi,qleft,qporq;
00110 static double T3,T6,T7,T8,T9,T10,T12,T13;
00111 /*
00112      ..
00113      .. Executable Statements ..
00114 */
00115 /*
00116      Check arguments
00117 */
00118     if(!(*which < 1 || *which > 4)) goto S30;
00119     if(!(*which < 1)) goto S10;
00120     *bound = 1.0e0;
00121     goto S20;
00122 S10:
00123     *bound = 4.0e0;
00124 S20:
00125     *status = -1;
00126     return;
00127 S30:
00128     if(*which == 1) goto S70;
00129 /*
00130      P
00131 */
00132     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00133     if(!(*p < 0.0e0)) goto S40;
00134     *bound = 0.0e0;
00135     goto S50;
00136 S40:
00137     *bound = 1.0e0;
00138 S50:
00139     *status = -2;
00140     return;
00141 S70:
00142 S60:
00143     if(*which == 1) goto S110;
00144 /*
00145      Q
00146 */
00147     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00148     if(!(*q <= 0.0e0)) goto S80;
00149     *bound = 0.0e0;
00150     goto S90;
00151 S80:
00152     *bound = 1.0e0;
00153 S90:
00154     *status = -3;
00155     return;
00156 S110:
00157 S100:
00158     if(*which == 2) goto S130;
00159 /*
00160      S
00161 */
00162     if(!(*s < 0.0e0)) goto S120;
00163     *bound = 0.0e0;
00164     *status = -4;
00165     return;
00166 S130:
00167 S120:
00168     if(*which == 3) goto S150;
00169 /*
00170      XN
00171 */
00172     if(!(*xn < 0.0e0)) goto S140;
00173     *bound = 0.0e0;
00174     *status = -5;
00175     return;
00176 S150:
00177 S140:
00178     if(*which == 4) goto S190;
00179 /*
00180      PR
00181 */
00182     if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
00183     if(!(*pr < 0.0e0)) goto S160;
00184     *bound = 0.0e0;
00185     goto S170;
00186 S160:
00187     *bound = 1.0e0;
00188 S170:
00189     *status = -6;
00190     return;
00191 S190:
00192 S180:
00193     if(*which == 4) goto S230;
00194 /*
00195      OMPR
00196 */
00197     if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
00198     if(!(*ompr < 0.0e0)) goto S200;
00199     *bound = 0.0e0;
00200     goto S210;
00201 S200:
00202     *bound = 1.0e0;
00203 S210:
00204     *status = -7;
00205     return;
00206 S230:
00207 S220:
00208     if(*which == 1) goto S270;
00209 /*
00210      P + Q
00211 */
00212     pq = *p+*q;
00213     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
00214     if(!(pq < 0.0e0)) goto S240;
00215     *bound = 0.0e0;
00216     goto S250;
00217 S240:
00218     *bound = 1.0e0;
00219 S250:
00220     *status = 3;
00221     return;
00222 S270:
00223 S260:
00224     if(*which == 4) goto S310;
00225 /*
00226      PR + OMPR
00227 */
00228     prompr = *pr+*ompr;
00229     if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
00230     if(!(prompr < 0.0e0)) goto S280;
00231     *bound = 0.0e0;
00232     goto S290;
00233 S280:
00234     *bound = 1.0e0;
00235 S290:
00236     *status = 4;
00237     return;
00238 S310:
00239 S300:
00240     if(!(*which == 1)) qporq = *p <= *q;
00241 /*
00242      Select the minimum of P or Q
00243      Calculate ANSWERS
00244 */
00245     if(1 == *which) {
00246 /*
00247      Calculating P
00248 */
00249         cumnbn(s,xn,pr,ompr,p,q);
00250         *status = 0;
00251     }
00252     else if(2 == *which) {
00253 /*
00254      Calculating S
00255 */
00256         *s = 5.0e0;
00257         T3 = inf;
00258         T6 = atol;
00259         T7 = tol;
00260         dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00261         *status = 0;
00262         dinvr(status,s,&fx,&qleft,&qhi);
00263 S320:
00264         if(!(*status == 1)) goto S350;
00265         cumnbn(s,xn,pr,ompr,&cum,&ccum);
00266         if(!qporq) goto S330;
00267         fx = cum-*p;
00268         goto S340;
00269 S330:
00270         fx = ccum-*q;
00271 S340:
00272         dinvr(status,s,&fx,&qleft,&qhi);
00273         goto S320;
00274 S350:
00275         if(!(*status == -1)) goto S380;
00276         if(!qleft) goto S360;
00277         *status = 1;
00278         *bound = 0.0e0;
00279         goto S370;
00280 S360:
00281         *status = 2;
00282         *bound = inf;
00283 S380:
00284 S370:
00285         ;
00286     }
00287     else if(3 == *which) {
00288 /*
00289      Calculating XN
00290 */
00291         *xn = 5.0e0;
00292         T8 = inf;
00293         T9 = atol;
00294         T10 = tol;
00295         dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
00296         *status = 0;
00297         dinvr(status,xn,&fx,&qleft,&qhi);
00298 S390:
00299         if(!(*status == 1)) goto S420;
00300         cumnbn(s,xn,pr,ompr,&cum,&ccum);
00301         if(!qporq) goto S400;
00302         fx = cum-*p;
00303         goto S410;
00304 S400:
00305         fx = ccum-*q;
00306 S410:
00307         dinvr(status,xn,&fx,&qleft,&qhi);
00308         goto S390;
00309 S420:
00310         if(!(*status == -1)) goto S450;
00311         if(!qleft) goto S430;
00312         *status = 1;
00313         *bound = 0.0e0;
00314         goto S440;
00315 S430:
00316         *status = 2;
00317         *bound = inf;
00318 S450:
00319 S440:
00320         ;
00321     }
00322     else if(4 == *which) {
00323 /*
00324      Calculating PR and OMPR
00325 */
00326         T12 = atol;
00327         T13 = tol;
00328         dstzr(&K2,&K11,&T12,&T13);
00329         if(!qporq) goto S480;
00330         *status = 0;
00331         dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
00332         *ompr = one-*pr;
00333 S460:
00334         if(!(*status == 1)) goto S470;
00335         cumnbn(s,xn,pr,ompr,&cum,&ccum);
00336         fx = cum-*p;
00337         dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
00338         *ompr = one-*pr;
00339         goto S460;
00340 S470:
00341         goto S510;
00342 S480:
00343         *status = 0;
00344         dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
00345         *pr = one-*ompr;
00346 S490:
00347         if(!(*status == 1)) goto S500;
00348         cumnbn(s,xn,pr,ompr,&cum,&ccum);
00349         fx = ccum-*q;
00350         dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
00351         *pr = one-*ompr;
00352         goto S490;
00353 S510:
00354 S500:
00355         if(!(*status == -1)) goto S540;
00356         if(!qleft) goto S520;
00357         *status = 1;
00358         *bound = 0.0e0;
00359         goto S530;
00360 S520:
00361         *status = 2;
00362         *bound = 1.0e0;
00363 S530:
00364         ;
00365     }
00366 S540:
00367     return;
00368 #undef tol
00369 #undef atol
00370 #undef inf
00371 #undef one
00372 } /* END */

00025                 : 1..4
00026                iwhich = 1 : Calculate P and Q from X,MEAN and SD
00027                iwhich = 2 : Calculate X from P,Q,MEAN and SD
00028                iwhich = 3 : Calculate MEAN from P,Q,X and SD
00029                iwhich = 4 : Calculate SD from P,Q,X and MEAN
00030 
00031      P <--> The integral from -infinity to X of the normal density.
00032             Input range: (0,1].
00033 
00034      Q <--> 1-P.
00035             Input range: (0, 1].
00036             P + Q = 1.0.
00037 
00038      X < --> Upper limit of integration of the normal-density.
00039              Input range: ( -infinity, +infinity)
00040 
00041      MEAN <--> The mean of the normal density.
00042                Input range: (-infinity, +infinity)
00043 
00044      SD <--> Standard Deviation of the normal density.
00045              Input range: (0, +infinity).
00046 
00047      STATUS <-- 0 if calculation completed correctly
00048                -I if input parameter number I is out of range
00049                 1 if answer appears to be lower than lowest
00050                   search bound
00051                 2 if answer appears to be higher than greatest
00052                   search bound
00053                 3 if P + Q .ne. 1
00054 
00055      BOUND <-- Undefined if STATUS is 0
00056 
00057                Bound exceeded by parameter number I if STATUS
00058                is negative.
00059 
00060                Lower search bound if STATUS is 1.
00061 
00062                Upper search bound if STATUS is 2.
00063 
00064 
00065                               Method
00066 
00067 
00068 
00069 
00070      A slightly modified version of ANORM from
00071 
00072      Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
00073      Package of Special Function Routines and Test Drivers"
00074      acm Transactions on Mathematical Software. 19, 22-32.
00075 
00076      is used to calulate the  cumulative standard normal distribution.
00077 
00078      The rational functions from pages  90-95  of Kennedy and Gentle,
00079      Statistical  Computing,  Marcel  Dekker, NY,  1980 are  used  as
00080      starting values to Newton's Iterations which compute the inverse
00081      standard normal.  Therefore no  searches  are necessary for  any
00082      parameter.
00083 
00084      For X < -15, the asymptotic expansion for the normal is used  as
00085      the starting value in finding the inverse standard normal.
00086      This is formula 26.2.12 of Abramowitz and Stegun.
00087 
00088 
00089                               Note
00090 
00091 
00092       The normal density is proportional to
00093       exp( - 0.5 * (( X - MEAN)/SD)**2)
00094 
00095 **********************************************************************/
00096 {
00097 static int K1 = 1;
00098 static double z,pq;
00099 /*
00100      ..
00101      .. Executable Statements ..
00102 */
00103 /*
00104      Check arguments
00105 */
00106     *status = 0;
00107     if(!(*which < 1 || *which > 4)) goto S30;
00108     if(!(*which < 1)) goto S10;
00109     *bound = 1.0e0;
00110     goto S20;
00111 S10:
00112     *bound = 4.0e0;
00113 S20:
00114     *status = -1;
00115     return;
00116 S30:
00117     if(*which == 1) goto S70;
00118 /*
00119      P
00120 */
00121     if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
00122     if(!(*p <= 0.0e0)) goto S40;
00123     *bound = 0.0e0;
00124     goto S50;
00125 S40:
00126     *bound = 1.0e0;
00127 S50:
00128     *status = -2;
00129     return;
00130 S70:
00131 S60:
00132     if(*which == 1) goto S110;
00133 /*
00134      Q
00135 */
00136     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00137     if(!(*q <= 0.0e0)) goto S80;
00138     *bound = 0.0e0;
00139     goto S90;
00140 S80:
00141     *bound = 1.0e0;
00142 S90:
00143     *status = -3;
00144     return;
00145 S110:
00146 S100:
00147     if(*which == 1) goto S150;
00148 /*
00149      P + Q
00150 */
00151     pq = *p+*q;
00152     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S140;
00153     if(!(pq < 0.0e0)) goto S120;
00154     *bound = 0.0e0;
00155     goto S130;
00156 S120:
00157     *bound = 1.0e0;
00158 S130:
00159     *status = 3;
00160     return;
00161 S150:
00162 S140:
00163     if(*which == 4) goto S170;
00164 /*
00165      SD
00166 */
00167     if(!(*sd <= 0.0e0)) goto S160;
00168     *bound = 0.0e0;
00169     *status = -6;
00170     return;
00171 S170:
00172 S160:
00173 /*
00174      Calculate ANSWERS
00175 */
00176     if(1 == *which) {
00177 /*
00178      Computing P
00179 */
00180         z = (*x-*mean)/ *sd;
00181         cumnor(&z,p,q);
00182     }
00183     else if(2 == *which) {
00184 /*
00185      Computing X
00186 */
00187         z = dinvnr(p,q);
00188         *x = *sd*z+*mean;
00189     }
00190     else if(3 == *which) {
00191 /*
00192      Computing the MEAN
00193 */
00194         z = dinvnr(p,q);
00195         *mean = *x-*sd*z;
00196     }
00197     else if(4 == *which) {
00198 /*
00199      Computing SD
00200 */
00201         z = dinvnr(p,q);
00202         *sd = (*x-*mean)/z;
00203     }
00204     return;
00205 } /* END */

00025                           : 1..3
00026                iwhich = 1 : Calculate P and Q from S and XLAM
00027                iwhich = 2 : Calculate A from P,Q and XLAM
00028                iwhich = 3 : Calculate XLAM from P,Q and S
00029 
00030         P <--> The cumulation from 0 to S of the poisson density.
00031                Input range: [0,1].
00032 
00033         Q <--> 1-P.
00034                Input range: (0, 1].
00035                P + Q = 1.0.
00036 
00037         S <--> Upper limit of cumulation of the Poisson.
00038                Input range: [0, +infinity).
00039                Search range: [0,1E300]
00040 
00041      XLAM <--> Mean of the Poisson distribution.
00042                Input range: [0, +infinity).
00043                Search range: [0,1E300]
00044 
00045      STATUS <-- 0 if calculation completed correctly
00046                -I if input parameter number I is out of range
00047                 1 if answer appears to be lower than lowest
00048                   search bound
00049                 2 if answer appears to be higher than greatest
00050                   search bound
00051                 3 if P + Q .ne. 1
00052 
00053      BOUND <-- Undefined if STATUS is 0
00054 
00055                Bound exceeded by parameter number I if STATUS
00056                is negative.
00057 
00058                Lower search bound if STATUS is 1.
00059 
00060                Upper search bound if STATUS is 2.
00061 
00062 
00063                               Method
00064 
00065 
00066      Formula   26.4.21  of   Abramowitz  and   Stegun,   Handbook  of
00067      Mathematical Functions (1966) is used  to reduce the computation
00068      of  the cumulative distribution function to that  of computing a
00069      chi-square, hence an incomplete gamma function.
00070 
00071      Cumulative  distribution function  (P) is  calculated  directly.
00072      Computation of other parameters involve a seach for a value that
00073      produces  the desired value of  P.   The  search relies  on  the
00074      monotinicity of P with the other parameter.
00075 
00076 **********************************************************************/
00077 {
00078 #define tol (1.0e-8)
00079 #define atol (1.0e-50)
00080 #define inf 1.0e300
00081 static int K1 = 1;
00082 static double K2 = 0.0e0;
00083 static double K4 = 0.5e0;
00084 static double K5 = 5.0e0;
00085 static double fx,cum,ccum,pq;
00086 static unsigned long qhi,qleft,qporq;
00087 static double T3,T6,T7,T8,T9,T10;
00088 /*
00089      ..
00090      .. Executable Statements ..
00091 */
00092 /*
00093      Check arguments
00094 */
00095     if(!(*which < 1 || *which > 3)) goto S30;
00096     if(!(*which < 1)) goto S10;
00097     *bound = 1.0e0;
00098     goto S20;
00099 S10:
00100     *bound = 3.0e0;
00101 S20:
00102     *status = -1;
00103     return;
00104 S30:
00105     if(*which == 1) goto S70;
00106 /*
00107      P
00108 */
00109     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00110     if(!(*p < 0.0e0)) goto S40;
00111     *bound = 0.0e0;
00112     goto S50;
00113 S40:
00114     *bound = 1.0e0;
00115 S50:
00116     *status = -2;
00117     return;
00118 S70:
00119 S60:
00120     if(*which == 1) goto S110;
00121 /*
00122      Q
00123 */
00124     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00125     if(!(*q <= 0.0e0)) goto S80;
00126     *bound = 0.0e0;
00127     goto S90;
00128 S80:
00129     *bound = 1.0e0;
00130 S90:
00131     *status = -3;
00132     return;
00133 S110:
00134 S100:
00135     if(*which == 2) goto S130;
00136 /*
00137      S
00138 */
00139     if(!(*s < 0.0e0)) goto S120;
00140     *bound = 0.0e0;
00141     *status = -4;
00142     return;
00143 S130:
00144 S120:
00145     if(*which == 3) goto S150;
00146 /*
00147      XLAM
00148 */
00149     if(!(*xlam < 0.0e0)) goto S140;
00150     *bound = 0.0e0;
00151     *status = -5;
00152     return;
00153 S150:
00154 S140:
00155     if(*which == 1) goto S190;
00156 /*
00157      P + Q
00158 */
00159     pq = *p+*q;
00160     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
00161     if(!(pq < 0.0e0)) goto S160;
00162     *bound = 0.0e0;
00163     goto S170;
00164 S160:
00165     *bound = 1.0e0;
00166 S170:
00167     *status = 3;
00168     return;
00169 S190:
00170 S180:
00171     if(!(*which == 1)) qporq = *p <= *q;
00172 /*
00173      Select the minimum of P or Q
00174      Calculate ANSWERS
00175 */
00176     if(1 == *which) {
00177 /*
00178      Calculating P
00179 */
00180         cumpoi(s,xlam,p,q);
00181         *status = 0;
00182     }
00183     else if(2 == *which) {
00184 /*
00185      Calculating S
00186 */
00187         *s = 5.0e0;
00188         T3 = inf;
00189         T6 = atol;
00190         T7 = tol;
00191         dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00192         *status = 0;
00193         dinvr(status,s,&fx,&qleft,&qhi);
00194 S200:
00195         if(!(*status == 1)) goto S230;
00196         cumpoi(s,xlam,&cum,&ccum);
00197         if(!qporq) goto S210;
00198         fx = cum-*p;
00199         goto S220;
00200 S210:
00201         fx = ccum-*q;
00202 S220:
00203         dinvr(status,s,&fx,&qleft,&qhi);
00204         goto S200;
00205 S230:
00206         if(!(*status == -1)) goto S260;
00207         if(!qleft) goto S240;
00208         *status = 1;
00209         *bound = 0.0e0;
00210         goto S250;
00211 S240:
00212         *status = 2;
00213         *bound = inf;
00214 S260:
00215 S250:
00216         ;
00217     }
00218     else if(3 == *which) {
00219 /*
00220      Calculating XLAM
00221 */
00222         *xlam = 5.0e0;
00223         T8 = inf;
00224         T9 = atol;
00225         T10 = tol;
00226         dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
00227         *status = 0;
00228         dinvr(status,xlam,&fx,&qleft,&qhi);
00229 S270:
00230         if(!(*status == 1)) goto S300;
00231         cumpoi(s,xlam,&cum,&ccum);
00232         if(!qporq) goto S280;
00233         fx = cum-*p;
00234         goto S290;
00235 S280:
00236         fx = ccum-*q;
00237 S290:
00238         dinvr(status,xlam,&fx,&qleft,&qhi);
00239         goto S270;
00240 S300:
00241         if(!(*status == -1)) goto S330;
00242         if(!qleft) goto S310;
00243         *status = 1;
00244         *bound = 0.0e0;
00245         goto S320;
00246 S310:
00247         *status = 2;
00248         *bound = inf;
00249 S320:
00250         ;
00251     }
00252 S330:
00253     return;
00254 #undef tol
00255 #undef atol
00256 #undef inf
00257 } /* END */

00025                           : 1..3
00026                iwhich = 1 : Calculate P and Q from T and DF
00027                iwhich = 2 : Calculate T from P,Q and DF
00028                iwhich = 3 : Calculate DF from P,Q and T
00029 
00030         P <--> The integral from -infinity to t of the t-density.
00031                Input range: (0,1].
00032 
00033         Q <--> 1-P.
00034                Input range: (0, 1].
00035                P + Q = 1.0.
00036 
00037         T <--> Upper limit of integration of the t-density.
00038                Input range: ( -infinity, +infinity).
00039                Search range: [ -1E300, 1E300 ]
00040 
00041         DF <--> Degrees of freedom of the t-distribution.
00042                 Input range: (0 , +infinity).
00043                 Search range: [1e-300, 1E10]
00044 
00045      STATUS <-- 0 if calculation completed correctly
00046                -I if input parameter number I is out of range
00047                 1 if answer appears to be lower than lowest
00048                   search bound
00049                 2 if answer appears to be higher than greatest
00050                   search bound
00051                 3 if P + Q .ne. 1
00052 
00053      BOUND <-- Undefined if STATUS is 0
00054 
00055                Bound exceeded by parameter number I if STATUS
00056                is negative.
00057 
00058                Lower search bound if STATUS is 1.
00059 
00060                Upper search bound if STATUS is 2.
00061 
00062 
00063                               Method
00064 
00065 
00066      Formula  26.5.27  of   Abramowitz   and  Stegun,   Handbook   of
00067      Mathematical Functions  (1966) is used to reduce the computation
00068      of the cumulative distribution function to that of an incomplete
00069      beta.
00070 
00071      Computation of other parameters involve a seach for a value that
00072      produces  the desired  value  of P.   The search relies  on  the
00073      monotinicity of P with the other parameter.
00074 
00075 **********************************************************************/
00076 {
00077 #define tol (1.0e-8)
00078 #define atol (1.0e-50)
00079 #define zero (1.0e-300)
00080 #define inf 1.0e300
00081 #define maxdf 1.0e10
00082 static int K1 = 1;
00083 static double K4 = 0.5e0;
00084 static double K5 = 5.0e0;
00085 static double fx,cum,ccum,pq;
00086 static unsigned long qhi,qleft,qporq;
00087 static double T2,T3,T6,T7,T8,T9,T10,T11;
00088 /*
00089      ..
00090      .. Executable Statements ..
00091 */
00092 /*
00093      Check arguments
00094 */
00095     if(!(*which < 1 || *which > 3)) goto S30;
00096     if(!(*which < 1)) goto S10;
00097     *bound = 1.0e0;
00098     goto S20;
00099 S10:
00100     *bound = 3.0e0;
00101 S20:
00102     *status = -1;
00103     return;
00104 S30:
00105     if(*which == 1) goto S70;
00106 /*
00107      P
00108 */
00109     if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
00110     if(!(*p <= 0.0e0)) goto S40;
00111     *bound = 0.0e0;
00112     goto S50;
00113 S40:
00114     *bound = 1.0e0;
00115 S50:
00116     *status = -2;
00117     return;
00118 S70:
00119 S60:
00120     if(*which == 1) goto S110;
00121 /*
00122      Q
00123 */
00124     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00125     if(!(*q <= 0.0e0)) goto S80;
00126     *bound = 0.0e0;
00127     goto S90;
00128 S80:
00129     *bound = 1.0e0;
00130 S90:
00131     *status = -3;
00132     return;
00133 S110:
00134 S100:
00135     if(*which == 3) goto S130;
00136 /*
00137      DF
00138 */
00139     if(!(*df <= 0.0e0)) goto S120;
00140     *bound = 0.0e0;
00141     *status = -5;
00142     return;
00143 S130:
00144 S120:
00145     if(*which == 1) goto S170;
00146 /*
00147      P + Q
00148 */
00149     pq = *p+*q;
00150     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S160;
00151     if(!(pq < 0.0e0)) goto S140;
00152     *bound = 0.0e0;
00153     goto S150;
00154 S140:
00155     *bound = 1.0e0;
00156 S150:
00157     *status = 3;
00158     return;
00159 S170:
00160 S160:
00161     if(!(*which == 1)) qporq = *p <= *q;
00162 /*
00163      Select the minimum of P or Q
00164      Calculate ANSWERS
00165 */
00166     if(1 == *which) {
00167 /*
00168      Computing P and Q
00169 */
00170         cumt(t,df,p,q);
00171         *status = 0;
00172     }
00173     else if(2 == *which) {
00174 /*
00175      Computing T
00176      .. Get initial approximation for T
00177 */
00178         *t = dt1(p,q,df);
00179         T2 = -inf;
00180         T3 = inf;
00181         T6 = atol;
00182         T7 = tol;
00183         dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
00184         *status = 0;
00185         dinvr(status,t,&fx,&qleft,&qhi);
00186 S180:
00187         if(!(*status == 1)) goto S210;
00188         cumt(t,df,&cum,&ccum);
00189         if(!qporq) goto S190;
00190         fx = cum-*p;
00191         goto S200;
00192 S190:
00193         fx = ccum-*q;
00194 S200:
00195         dinvr(status,t,&fx,&qleft,&qhi);
00196         goto S180;
00197 S210:
00198         if(!(*status == -1)) goto S240;
00199         if(!qleft) goto S220;
00200         *status = 1;
00201         *bound = -inf;
00202         goto S230;
00203 S220:
00204         *status = 2;
00205         *bound = inf;
00206 S240:
00207 S230:
00208         ;
00209     }
00210     else if(3 == *which) {
00211 /*
00212      Computing DF
00213 */
00214         *df = 5.0e0;
00215         T8 = zero;
00216         T9 = maxdf;
00217         T10 = atol;
00218         T11 = tol;
00219         dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
00220         *status = 0;
00221         dinvr(status,df,&fx,&qleft,&qhi);
00222 S250:
00223         if(!(*status == 1)) goto S280;
00224         cumt(t,df,&cum,&ccum);
00225         if(!qporq) goto S260;
00226         fx = cum-*p;
00227         goto S270;
00228 S260:
00229         fx = ccum-*q;
00230 S270:
00231         dinvr(status,df,&fx,&qleft,&qhi);
00232         goto S250;
00233 S280:
00234         if(!(*status == -1)) goto S310;
00235         if(!qleft) goto S290;
00236         *status = 1;
00237         *bound = zero;
00238         goto S300;
00239 S290:
00240         *status = 2;
00241         *bound = maxdf;
00242 S300:
00243         ;
00244     }
00245 S310:
00246     return;
00247 #undef tol
00248 #undef atol
00249 #undef zero
00250 #undef inf
00251 #undef maxdf
00252 } /* END */

00056 {
00057 static int ierr;
00058 /*
00059      ..
00060      .. Executable Statements ..
00061 */
00062     if(!(*x <= 0.0e0)) goto S10;
00063     *cum = 0.0e0;
00064     *ccum = 1.0e0;
00065     return;
00066 S10:
00067     if(!(*y <= 0.0e0)) goto S20;
00068     *cum = 1.0e0;
00069     *ccum = 0.0e0;
00070     return;
00071 S20:
00072     bratio(a,b,x,y,cum,ccum,&ierr);
00073 /*
00074      Call bratio routine
00075 */
00076     return;
00077 } /* END */

00051 {
00052 static double T1,T2;
00053 /*
00054      ..
00055      .. Executable Statements ..
00056 */
00057     if(!(*s < *xn)) goto S10;
00058     T1 = *s+1.0e0;
00059     T2 = *xn-*s;
00060     cumbet(pr,ompr,&T1,&T2,ccum,cum);
00061     goto S20;
00062 S10:
00063     *cum = 1.0e0;
00064     *ccum = 0.0e0;
00065 S20:
00066     return;
00067 } /* END */

00041 {
00042 static double a,xx;
00043 /*
00044      ..
00045      .. Executable Statements ..
00046 */
00047     a = *df*0.5e0;
00048     xx = *x*0.5e0;
00049     cumgam(&xx,&a,cum,ccum);
00050     return;
00051 } /* END */

00074 {
00075 #define dg(i) (*df+2.0e0*(double)(i))
00076 #define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps*sum)
00077 #define qtired(i) (int)((i) > ntired)
00078 static double eps = 1.0e-5;
00079 static int ntired = 1000;
00080 static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
00081     sumadj,term,wt,xnonc;
00082 static int i,icent,iterb,iterf;
00083 static double T1,T2,T3;
00084 /*
00085      ..
00086      .. Executable Statements ..
00087 */
00088     if(!(*x <= 0.0e0)) goto S10;
00089     *cum = 0.0e0;
00090     *ccum = 1.0e0;
00091     return;
00092 S10:
00093     if(!(*pnonc <= 1.0e-10)) goto S20;
00094 /*
00095      When non-centrality parameter is (essentially) zero,
00096      use cumulative chi-square distribution
00097 */
00098     cumchi(x,df,cum,ccum);
00099     return;
00100 S20:
00101     xnonc = *pnonc/2.0e0;
00102 /*
00103 **********************************************************************
00104      The following code calcualtes the weight, chi-square, and
00105      adjustment term for the central term in the infinite series.
00106      The central term is the one in which the poisson weight is
00107      greatest.  The adjustment term is the amount that must
00108      be subtracted from the chi-square to move up two degrees
00109      of freedom.
00110 **********************************************************************
00111 */
00112     icent = fifidint(xnonc);
00113     if(icent == 0) icent = 1;
00114     chid2 = *x/2.0e0;
00115 /*
00116      Calculate central weight term
00117 */
00118     T1 = (double)(icent+1);
00119     lfact = alngam(&T1);
00120     lcntwt = -xnonc+(double)icent*log(xnonc)-lfact;
00121     centwt = exp(lcntwt);
00122 /*
00123      Calculate central chi-square
00124 */
00125     T2 = dg(icent);
00126     cumchi(x,&T2,&pcent,ccum);
00127 /*
00128      Calculate central adjustment term
00129 */
00130     dfd2 = dg(icent)/2.0e0;
00131     T3 = 1.0e0+dfd2;
00132     lfact = alngam(&T3);
00133     lcntaj = dfd2*log(chid2)-chid2-lfact;
00134     centaj = exp(lcntaj);
00135     sum = centwt*pcent;
00136 /*
00137 **********************************************************************
00138      Sum backwards from the central term towards zero.
00139      Quit whenever either
00140      (1) the zero term is reached, or
00141      (2) the term gets small relative to the sum, or
00142      (3) More than NTIRED terms are totaled.
00143 **********************************************************************
00144 */
00145     iterb = 0;
00146     sumadj = 0.0e0;
00147     adj = centaj;
00148     wt = centwt;
00149     i = icent;
00150     goto S40;
00151 S30:
00152     if(qtired(iterb) || qsmall(term) || i == 0) goto S50;
00153 S40:
00154     dfd2 = dg(i)/2.0e0;
00155 /*
00156      Adjust chi-square for two fewer degrees of freedom.
00157      The adjusted value ends up in PTERM.
00158 */
00159     adj = adj*dfd2/chid2;
00160     sumadj += adj;
00161     pterm = pcent+sumadj;
00162 /*
00163      Adjust poisson weight for J decreased by one
00164 */
00165     wt *= ((double)i/xnonc);
00166     term = wt*pterm;
00167     sum += term;
00168     i -= 1;
00169     iterb += 1;
00170     goto S30;
00171 S50:
00172     iterf = 0;
00173 /*
00174 **********************************************************************
00175      Now sum forward from the central term towards infinity.
00176      Quit when either
00177      (1) the term gets small relative to the sum, or
00178      (2) More than NTIRED terms are totaled.
00179 **********************************************************************
00180 */
00181     sumadj = adj = centaj;
00182     wt = centwt;
00183     i = icent;
00184     goto S70;
00185 S60:
00186     if(qtired(iterf) || qsmall(term)) goto S80;
00187 S70:
00188 /*
00189      Update weights for next higher J
00190 */
00191     wt *= (xnonc/(double)(i+1));
00192 /*
00193      Calculate PTERM and add term to sum
00194 */
00195     pterm = pcent-sumadj;
00196     term = wt*pterm;
00197     sum += term;
00198 /*
00199      Update adjustment term for DF for next iteration
00200 */
00201     i += 1;
00202     dfd2 = dg(i)/2.0e0;
00203     adj = adj*chid2/dfd2;
00204     sumadj += adj;
00205     iterf += 1;
00206     goto S60;
00207 S80:
00208     *cum = sum;
00209     *ccum = 0.5e0+(0.5e0-*cum);
00210     return;
00211 #undef dg
00212 #undef qsmall
00213 #undef qtired
00214 } /* END */

00050 {
00051 #define half 0.5e0
00052 #define done 1.0e0
00053 static double dsum,prod,xx,yy;
00054 static int ierr;
00055 static double T1,T2;
00056 /*
00057      ..
00058      .. Executable Statements ..
00059 */
00060     if(!(*f <= 0.0e0)) goto S10;
00061     *cum = 0.0e0;
00062     *ccum = 1.0e0;
00063     return;
00064 S10:
00065     prod = *dfn**f;
00066 /*
00067      XX is such that the incomplete beta with parameters
00068      DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM
00069      YY is 1 - XX
00070      Calculate the smaller of XX and YY accurately
00071 */
00072     dsum = *dfd+prod;
00073     xx = *dfd/dsum;
00074     if(xx > half) {
00075         yy = prod/dsum;
00076         xx = done-yy;
00077     }
00078     else  yy = done-xx;
00079     T1 = *dfd*half;
00080     T2 = *dfn*half;
00081     bratio(&T1,&T2,&xx,&yy,ccum,cum,&ierr);
00082     return;
00083 #undef half
00084 #undef done
00085 } /* END */

00065 {
00066 #define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
00067 #define half 0.5e0
00068 #define done 1.0e0
00069 static double eps = 1.0e-4;
00070 static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
00071     upterm,xmult,xnonc;
00072 static int i,icent,ierr;
00073 static double T1,T2,T3,T4,T5,T6;
00074 /*
00075      ..
00076      .. Executable Statements ..
00077 */
00078     if(!(*f <= 0.0e0)) goto S10;
00079     *cum = 0.0e0;
00080     *ccum = 1.0e0;
00081     return;
00082 S10:
00083     if(!(*pnonc < 1.0e-10)) goto S20;
00084 /*
00085      Handle case in which the non-centrality parameter is
00086      (essentially) zero.
00087 */
00088     cumf(f,dfn,dfd,cum,ccum);
00089     return;
00090 S20:
00091     xnonc = *pnonc/2.0e0;
00092 /*
00093      Calculate the central term of the poisson weighting factor.
00094 */
00095     icent = xnonc;
00096     if(icent == 0) icent = 1;
00097 /*
00098      Compute central weight term
00099 */
00100     T1 = (double)(icent+1);
00101     centwt = exp(-xnonc+(double)icent*log(xnonc)-alngam(&T1));
00102 /*
00103      Compute central incomplete beta term
00104      Assure that minimum of arg to beta and 1 - arg is computed
00105           accurately.
00106 */
00107     prod = *dfn**f;
00108     dsum = *dfd+prod;
00109     yy = *dfd/dsum;
00110     if(yy > half) {
00111         xx = prod/dsum;
00112         yy = done-xx;
00113     }
00114     else  xx = done-yy;
00115     T2 = *dfn*half+(double)icent;
00116     T3 = *dfd*half;
00117     bratio(&T2,&T3,&xx,&yy,&betdn,&dummy,&ierr);
00118     adn = *dfn/2.0e0+(double)icent;
00119     aup = adn;
00120     b = *dfd/2.0e0;
00121     betup = betdn;
00122     sum = centwt*betdn;
00123 /*
00124      Now sum terms backward from icent until convergence or all done
00125 */
00126     xmult = centwt;
00127     i = icent;
00128     T4 = adn+b;
00129     T5 = adn+1.0e0;
00130     dnterm = exp(alngam(&T4)-alngam(&T5)-alngam(&b)+adn*log(xx)+b*log(yy));
00131 S30:
00132     if(qsmall(xmult*betdn) || i <= 0) goto S40;
00133     xmult *= ((double)i/xnonc);
00134     i -= 1;
00135     adn -= 1.0;
00136     dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
00137     betdn += dnterm;
00138     sum += (xmult*betdn);
00139     goto S30;
00140 S40:
00141     i = icent+1;
00142 /*
00143      Now sum forwards until convergence
00144 */
00145     xmult = centwt;
00146     if(aup-1.0+b == 0) upterm = exp(-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+
00147       b*log(yy));
00148     else  {
00149         T6 = aup-1.0+b;
00150         upterm = exp(alngam(&T6)-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+b*
00151           log(yy));
00152     }
00153     goto S60;
00154 S50:
00155     if(qsmall(xmult*betup)) goto S70;
00156 S60:
00157     xmult *= (xnonc/(double)i);
00158     i += 1;
00159     aup += 1.0;
00160     upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
00161     betup -= upterm;
00162     sum += (xmult*betup);
00163     goto S50;
00164 S70:
00165     *cum = sum;
00166     *ccum = 0.5e0+(0.5e0-*cum);
00167     return;
00168 #undef qsmall
00169 #undef half
00170 #undef done
00171 } /* END */

00044 {
00045 static int K1 = 0;
00046 /*
00047      ..
00048      .. Executable Statements ..
00049 */
00050     if(!(*x <= 0.0e0)) goto S10;
00051     *cum = 0.0e0;
00052     *ccum = 1.0e0;
00053     return;
00054 S10:
00055     gratio(a,x,cum,ccum,&K1);
00056 /*
00057      Call gratio routine
00058 */
00059     return;
00060 } /* END */

00057 {
00058 static double T1;
00059 /*
00060      ..
00061      .. Executable Statements ..
00062 */
00063     T1 = *s+1.e0;
00064     cumbet(pr,ompr,xn,&T1,cum,ccum);
00065     return;
00066 } /* END */

00025                                     :
00026 
00027      Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
00028      Package of Special Function Routines and Test Drivers"
00029      acm Transactions on Mathematical Software. 19, 22-32.
00030 
00031      with slight modifications to return ccum and to deal with
00032      machine constants.
00033  
00034 **********************************************************************
00035   Original Comments:
00036 ------------------------------------------------------------------
00037  
00038  This function evaluates the normal distribution function:
00039  
00040                               / x
00041                      1       |       -t*t/2
00042           P(x) = ----------- |      e       dt
00043                  sqrt(2 pi)  |
00044                              /-oo
00045  
00046    The main computation evaluates near-minimax approximations
00047    derived from those in "Rational Chebyshev approximations for
00048    the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
00049    This transportable program uses rational functions that
00050    theoretically approximate the normal distribution function to
00051    at least 18 significant decimal digits.  The accuracy achieved
00052    depends on the arithmetic system, the compiler, the intrinsic
00053    functions, and proper selection of the machine-dependent
00054    constants.
00055  
00056 *******************************************************************
00057 *******************************************************************
00058  
00059  Explanation of machine-dependent constants.
00060  
00061    MIN   = smallest machine representable number.
00062  
00063    EPS   = argument below which anorm(x) may be represented by
00064            0.5  and above which  x*x  will not underflow.
00065            A conservative value is the largest machine number X
00066            such that   1.0 + X = 1.0   to machine precision.
00067 *******************************************************************
00068 *******************************************************************
00069  
00070  Error returns
00071  
00072   The program returns  ANORM = 0     for  ARG .LE. XLOW.
00073  
00074  
00075  Intrinsic functions required are:
00076  
00077      ABS, AINT, EXP
00078  
00079  
00080   Author: W. J. Cody
00081           Mathematics and Computer Science Division
00082           Argonne National Laboratory
00083           Argonne, IL 60439
00084  
00085   Latest modification: March 15, 1992
00086  
00087 ------------------------------------------------------------------
00088 */
00089 {
00090 static double a[5] = {
00091     2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
00092     1.8154981253343561249e04,6.5682337918207449113e-2
00093 };
00094 static double b[4] = {
00095     4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
00096     4.5507789335026729956e04
00097 };
00098 static double c[9] = {
00099     3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
00100     5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
00101     1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
00102 };
00103 static double d[8] = {
00104     2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
00105     6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
00106     3.8912003286093271411e04,1.9685429676859990727e04
00107 };
00108 static double half = 0.5e0;
00109 static double p[6] = {
00110     2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
00111     1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
00112 };
00113 static double one = 1.0e0;
00114 static double q[5] = {
00115     1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
00116     3.78239633202758244e-3,7.29751555083966205e-5
00117 };
00118 static double sixten = 1.60e0;
00119 static double sqrpi = 3.9894228040143267794e-1;
00120 static double thrsh = 0.66291e0;
00121 static double root32 = 5.656854248e0;
00122 static double zero = 0.0e0;
00123 static int K1 = 1;
00124 static int K2 = 2;
00125 static int i;
00126 static double del,eps,temp,x,xden,xnum,y,xsq,min;
00127 /*
00128 ------------------------------------------------------------------
00129   Machine dependent constants
00130 ------------------------------------------------------------------
00131 */
00132     eps = spmpar(&K1)*0.5e0;
00133     min = spmpar(&K2);
00134     x = *arg;
00135     y = fabs(x);
00136     if(y <= thrsh) {
00137 /*
00138 ------------------------------------------------------------------
00139   Evaluate  anorm  for  |X| <= 0.66291
00140 ------------------------------------------------------------------
00141 */
00142         xsq = zero;
00143         if(y > eps) xsq = x*x;
00144         xnum = a[4]*xsq;
00145         xden = xsq;
00146         for(i=0; i<3; i++) {
00147             xnum = (xnum+a[i])*xsq;
00148             xden = (xden+b[i])*xsq;
00149         }
00150         *result = x*(xnum+a[3])/(xden+b[3]);
00151         temp = *result;
00152         *result = half+temp;
00153         *ccum = half-temp;
00154     }
00155 /*
00156 ------------------------------------------------------------------
00157   Evaluate  anorm  for 0.66291 <= |X| <= sqrt(32)
00158 ------------------------------------------------------------------
00159 */
00160     else if(y <= root32) {
00161         xnum = c[8]*y;
00162         xden = y;
00163         for(i=0; i<7; i++) {
00164             xnum = (xnum+c[i])*y;
00165             xden = (xden+d[i])*y;
00166         }
00167         *result = (xnum+c[7])/(xden+d[7]);
00168         xsq = fifdint(y*sixten)/sixten;
00169         del = (y-xsq)*(y+xsq);
00170         *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
00171         *ccum = one-*result;
00172         if(x > zero) {
00173             temp = *result;
00174             *result = *ccum;
00175             *ccum = temp;
00176         }
00177     }
00178 /*
00179 ------------------------------------------------------------------
00180   Evaluate  anorm  for |X| > sqrt(32)
00181 ------------------------------------------------------------------
00182 */
00183     else  {
00184         *result = zero;
00185         xsq = one/(x*x);
00186         xnum = p[5]*xsq;
00187         xden = xsq;
00188         for(i=0; i<4; i++) {
00189             xnum = (xnum+p[i])*xsq;
00190             xden = (xden+q[i])*xsq;
00191         }
00192         *result = xsq*(xnum+p[4])/(xden+q[4]);
00193         *result = (sqrpi-*result)/y;
00194         xsq = fifdint(x*sixten)/sixten;
00195         del = (x-xsq)*(x+xsq);
00196         *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
00197         *ccum = one-*result;
00198         if(x > zero) {
00199             temp = *result;
00200             *result = *ccum;
00201             *ccum = temp;
00202         }
00203     }
00204     if(*result < min) *result = 0.0e0;
00205 /*
00206 ------------------------------------------------------------------
00207   Fix up for negative argument, erf, etc.
00208 ------------------------------------------------------------------
00209 ----------Last card of ANORM ----------
00210 */
00211     if(*ccum < min) *ccum = 0.0e0;
} /* END */

00042 {
00043 static double chi,df;
00044 /*
00045      ..
00046      .. Executable Statements ..
00047 */
00048     df = 2.0e0*(*s+1.0e0);
00049     chi = 2.0e0**xlam;
00050     cumchi(&chi,&df,ccum,cum);
00051     return;
00052 } /* END */

00041 {
00042 static double K2 = 0.5e0;
00043 static double xx,a,oma,tt,yy,dfptt,T1;
00044 /*
00045      ..
00046      .. Executable Statements ..
00047 */
00048     tt = *t**t;
00049     dfptt = *df+tt;
00050     xx = *df/dfptt;
00051     yy = tt/dfptt;
00052     T1 = 0.5e0**df;
00053     cumbet(&xx,&yy,&T1,&K2,&a,&oma);
00054     if(!(*t <= 0.0e0)) goto S10;
00055     *cum = 0.5e0*a;
00056     *ccum = oma+*cum;
00057     goto S20;
00058 S10:
00059     *ccum = 0.5e0*a;
00060     *cum = oma+*ccum;
00061 S20:
00062     return;
00063 } /* END */

00035 {
00036 static double dbetrm,T1,T2,T3;
00037 /*
00038      ..
00039      .. Executable Statements ..
00040 */
00041 /*
00042      Try to sum from smallest to largest
00043 */
00044     T1 = *a+*b;
00045     dbetrm = -dstrem(&T1);
00046     T2 = fifdmax1(*a,*b);
00047     dbetrm += dstrem(&T2);
00048     T3 = fifdmin1(*a,*b);
00049     dbetrm += dstrem(&T3);
00050     return dbetrm;
00051 } /* END */

00031 {
00032 static double devlpl,term;
00033 static int i;
00034 /*
00035      ..
00036      .. Executable Statements ..
00037 */
00038     term = a[*n-1];
00039     for(i= *n-1-1; i>=0; i--) term = a[i]+term**x;
00040     devlpl = term;
00041     return devlpl;
00042 } /* END */

00020                                            :
00021  
00022      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00023      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00024      Trans. Math.  Softw. 18 (1993), 360-373.
00025  
00026 **********************************************************************
00027 */
00028 {
00029 static double p1 = .914041914819518e-09;
00030 static double p2 = .238082361044469e-01;
00031 static double q1 = -.499999999085958e+00;
00032 static double q2 = .107141568980644e+00;
00033 static double q3 = -.119041179760821e-01;
00034 static double q4 = .595130811860248e-03;
00035 static double dexpm1,w;
00036 /*
00037      ..
00038      .. Executable Statements ..
00039 */
00040     if(fabs(*x) > 0.15e0) goto S10;
00041     dexpm1 = *x*(((p2**x+p1)**x+1.0e0)/((((q4**x+q3)**x+q2)**x+q1)**x+1.0e0));
00042     return dexpm1;
00043 S10:
00044     w = exp(*x);
00045     if(*x > 0.0e0) goto S20;
00046     dexpm1 = w-0.5e0-0.5e0;
00047     return dexpm1;
00048 S20:
00049     dexpm1 = w*(0.5e0+(0.5e0-1.0e0/w));
00050     return dexpm1;
} /* END */

00042 {
00043 #define maxit 100
00044 #define eps (1.0e-13)
00045 #define r2pi 0.3989422804014326e0
00046 #define nhalf (-0.5e0)
00047 #define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
00048 static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
00049 static int i;
00050 static unsigned long qporq;
00051 /*
00052      ..
00053      .. Executable Statements ..
00054 */
00055 /*
00056      FIND MINIMUM OF P AND Q
00057 */
00058     qporq = *p <= *q;
00059     if(!qporq) goto S10;
00060     pp = *p;
00061     goto S20;
00062 S10:
00063     pp = *q;
00064 S20:
00065 /*
00066      INITIALIZATION STEP
00067 */
00068     strtx = stvaln(&pp);
00069     xcur = strtx;
00070 /*
00071      NEWTON INTERATIONS
00072 */
00073     for(i=1; i<=maxit; i++) {
00074         cumnor(&xcur,&cum,&ccum);
00075         dx = (cum-pp)/dennor(xcur);
00076         xcur -= dx;
00077         if(fabs(dx/xcur) < eps) goto S40;
00078     }
00079     dinvnr = strtx;
00080 /*
00081      IF WE GET HERE, NEWTON HAS FAILED
00082 */
00083     if(!qporq) dinvnr = -dinvnr;
00084     return dinvnr;
00085 S40:
00086 /*
00087      IF WE GET HERE, NEWTON HAS SUCCEDED
00088 */
00089     dinvnr = xcur;
00090     if(!qporq) dinvnr = -dinvnr;
00091     return dinvnr;
00092 #undef maxit
00093 #undef eps
00094 #undef r2pi
00095 #undef nhalf
00096 #undef dennor
00097 } /* END */

00064 {
00065     E0000(0,status,x,fx,qleft,qhi,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
00066 } /* END */

00038 {
00039 #define dlsqpi 0.91893853320467274177e0
00040 static double coef[12] = {
00041     -1.0e0,3.0e0,-15.0e0,105.0e0,-945.0e0,10395.0e0,-135135.0e0,2027025.0e0,
00042     -34459425.0e0,654729075.0e0,-13749310575.e0,316234143225.0e0
00043 };
00044 static int K1 = 12;
00045 static double dlanor,approx,correc,xx,xx2,T2;
00046 /*
00047      ..
00048      .. Executable Statements ..
00049 */
00050     xx = fabs(*x);
00051     if(xx < 5.0e0) ftnstop(" Argument too small in DLANOR");
00052     approx = -dlsqpi-0.5e0*xx*xx-log(xx);
00053     xx2 = xx*xx;
00054     T2 = 1.0e0/xx2;
00055     correc = devlpl(coef,&K1,&T2)/xx2;
00056     correc = dln1px(&correc);
00057     dlanor = approx+correc;
00058     return dlanor;
00059 #undef dlsqpi
00060 } /* END */

00034 {
00035 static double dln1mx,T1;
00036 /*
00037      ..
00038      .. Executable Statements ..
00039 */
00040     T1 = -*x;
00041     dln1mx = dln1px(&T1);
00042     return dln1mx;
00043 } /* END */

00029                         :
00030      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00031      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00032      Trans. Math.  Softw. 18 (1993), 360-373.
00033  
00034 **********************************************************************
00035 -----------------------------------------------------------------------
00036             EVALUATION OF THE FUNCTION LN(1 + A)
00037 -----------------------------------------------------------------------
00038 */
00039 {
00040 static double p1 = -.129418923021993e+01;
00041 static double p2 = .405303492862024e+00;
00042 static double p3 = -.178874546012214e-01;
00043 static double q1 = -.162752256355323e+01;
00044 static double q2 = .747811014037616e+00;
00045 static double q3 = -.845104217945565e-01;
00046 static double dln1px,t,t2,w,x;
00047 /*
00048      ..
00049      .. Executable Statements ..
00050 */
00051     if(fabs(*a) > 0.375e0) goto S10;
00052     t = *a/(*a+2.0e0);
00053     t2 = t*t;
00054     w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);
00055     dln1px = 2.0e0*t*w;
00056     return dln1px;
00057 S10:
00058     x = 1.e0+*a;
00059     dln1px = log(x);
00060     return dln1px;
} /* END */

00029                         :
00030      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00031      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00032      Trans. Math.  Softw. 18 (1993), 360-373.
00033  
00034 **********************************************************************
00035 -----------------------------------------------------------------------
00036      EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
00037 -----------------------------------------------------------------------
00038      E = 0.5*LN(2*PI)
00039 --------------------------
00040 */
00041 {
00042 static double e = .918938533204673e0;
00043 static double dlnbet,a,b,c,h,u,v,w,z;
00044 static int i,n;
00045 static double T1;
00046 /*
00047      ..
00048      .. Executable Statements ..
00049 */
00050     a = fifdmin1(*a0,*b0);
00051     b = fifdmax1(*a0,*b0);
00052     if(a >= 8.0e0) goto S100;
00053     if(a >= 1.0e0) goto S20;
00054 /*
00055 -----------------------------------------------------------------------
00056                    PROCEDURE WHEN A .LT. 1
00057 -----------------------------------------------------------------------
00058 */
00059     if(b >= 8.0e0) goto S10;
00060     T1 = a+b;
00061     dlnbet = gamln(&a)+(gamln(&b)-gamln(&T1));
00062     return dlnbet;
00063 S10:
00064     dlnbet = gamln(&a)+algdiv(&a,&b);
00065     return dlnbet;
00066 S20:
00067 /*
00068 -----------------------------------------------------------------------
00069                 PROCEDURE WHEN 1 .LE. A .LT. 8
00070 -----------------------------------------------------------------------
00071 */
00072     if(a > 2.0e0) goto S40;
00073     if(b > 2.0e0) goto S30;
00074     dlnbet = gamln(&a)+gamln(&b)-gsumln(&a,&b);
00075     return dlnbet;
00076 S30:
00077     w = 0.0e0;
00078     if(b < 8.0e0) goto S60;
00079     dlnbet = gamln(&a)+algdiv(&a,&b);
00080     return dlnbet;
00081 S40:
00082 /*
00083                 REDUCTION OF A WHEN B .LE. 1000
00084 */
00085     if(b > 1000.0e0) goto S80;
00086     n = a-1.0e0;
00087     w = 1.0e0;
00088     for(i=1; i<=n; i++) {
00089         a -= 1.0e0;
00090         h = a/b;
00091         w *= (h/(1.0e0+h));
00092     }
00093     w = log(w);
00094     if(b < 8.0e0) goto S60;
00095     dlnbet = w+gamln(&a)+algdiv(&a,&b);
00096     return dlnbet;
00097 S60:
00098 /*
00099                  REDUCTION OF B WHEN B .LT. 8
00100 */
00101     n = b-1.0e0;
00102     z = 1.0e0;
00103     for(i=1; i<=n; i++) {
00104         b -= 1.0e0;
00105         z *= (b/(a+b));
00106     }
00107     dlnbet = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
00108     return dlnbet;
00109 S80:
00110 /*
00111                 REDUCTION OF A WHEN B .GT. 1000
00112 */
00113     n = a-1.0e0;
00114     w = 1.0e0;
00115     for(i=1; i<=n; i++) {
00116         a -= 1.0e0;
00117         w *= (a/(1.0e0+a/b));
00118     }
00119     dlnbet = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
00120     return dlnbet;
00121 S100:
00122 /*
00123 -----------------------------------------------------------------------
00124                    PROCEDURE WHEN A .GE. 8
00125 -----------------------------------------------------------------------
00126 */
00127     w = bcorr(&a,&b);
00128     h = a/b;
00129     c = h/(1.0e0+h);
00130     u = -((a-0.5e0)*log(c));
00131     v = b*alnrel(&h);
00132     if(u <= v) goto S110;
00133     dlnbet = -(0.5e0*log(b))+e+w-v-u;
00134     return dlnbet;
00135 S110:
00136     dlnbet = -(0.5e0*log(b))+e+w-u-v;
00137     return dlnbet;
} /* END */

00026                        :
00027      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00028      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00029      Trans. Math.  Softw. 18 (1993), 360-373.
00030  
00031 **********************************************************************
00032 -----------------------------------------------------------------------
00033             EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
00034 -----------------------------------------------------------------------
00035      WRITTEN BY ALFRED H. MORRIS
00036           NAVAL SURFACE WARFARE CENTER
00037           DAHLGREN, VIRGINIA
00038 --------------------------
00039      D = 0.5*(LN(2*PI) - 1)
00040 --------------------------
00041 */
00042 {
00043 static double c0 = .833333333333333e-01;
00044 static double c1 = -.277777777760991e-02;
00045 static double c2 = .793650666825390e-03;
00046 static double c3 = -.595202931351870e-03;
00047 static double c4 = .837308034031215e-03;
00048 static double c5 = -.165322962780713e-02;
00049 static double d = .418938533204673e0;
00050 static double dlngam,t,w;
00051 static int i,n;
00052 static double T1;
00053 /*
00054      ..
00055      .. Executable Statements ..
00056 */
00057     if(*a > 0.8e0) goto S10;
00058     dlngam = gamln1(a)-log(*a);
00059     return dlngam;
00060 S10:
00061     if(*a > 2.25e0) goto S20;
00062     t = *a-0.5e0-0.5e0;
00063     dlngam = gamln1(&t);
00064     return dlngam;
00065 S20:
00066     if(*a >= 10.0e0) goto S40;
00067     n = *a-1.25e0;
00068     t = *a;
00069     w = 1.0e0;
00070     for(i=1; i<=n; i++) {
00071         t -= 1.0e0;
00072         w = t*w;
00073     }
00074     T1 = t-1.0e0;
00075     dlngam = gamln1(&T1)+log(w);
00076     return dlngam;
00077 S40:
00078     t = pow(1.0e0/ *a,2.0);
00079     w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a;
00080     dlngam = d+w+(*a-0.5e0)*(log(*a)-1.0e0);
00081     return dlngam;
} /* END */

00069 {
00070     E0000(1,NULL,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
00071     zstpmu);
00072 } /* END */

00003 {
00004 /*
00005 **********************************************************************
00006      double dstrem(double *z)
00007              Double precision Sterling Remainder
00008                               Function
00009      Returns   Log(Gamma(Z))  -  Sterling(Z)  where   Sterling(Z)  is
00010      Sterling's Approximation to Log(Gamma(Z))
00011      Sterling(Z) = LOG( SQRT( 2*PI ) ) + ( Z-0.5 ) * LOG( Z ) - Z
00012                               Arguments
00013      Z --> Value at which Sterling remainder calculated
00014            Must be positive.
00015                   DOUBLE PRECISION Z
00016                               Method
00017      If Z >= 6 uses 9 terms of series in Bernoulli numbers
00018      (Values calculated using Maple)
00019      Otherwise computes difference explicitly
00020 **********************************************************************
00021 */
00022 #define hln2pi 0.91893853320467274178e0
00023 #define ncoef 10
00024 static double coef[ncoef] = {
00025     0.0e0,0.0833333333333333333333333333333e0,
00026     -0.00277777777777777777777777777778e0,0.000793650793650793650793650793651e0,
00027     -0.000595238095238095238095238095238e0,
00028     0.000841750841750841750841750841751e0,-0.00191752691752691752691752691753e0,
00029     0.00641025641025641025641025641026e0,-0.0295506535947712418300653594771e0,
00030     0.179644372368830573164938490016e0
00031 };
00032 static int K1 = 10;
00033 static double dstrem,sterl,T2;
00034 /*
00035      ..
00036      .. Executable Statements ..
00037 */
00038 /*
00039     For information, here are the next 11 coefficients of the
00040     remainder term in Sterling's formula
00041             -1.39243221690590111642743221691
00042             13.4028640441683919944789510007
00043             -156.848284626002017306365132452
00044             2193.10333333333333333333333333
00045             -36108.7712537249893571732652192
00046             691472.268851313067108395250776
00047             -0.152382215394074161922833649589D8
00048             0.382900751391414141414141414141D9
00049             -0.108822660357843910890151491655D11
00050             0.347320283765002252252252252252D12
00051             -0.123696021422692744542517103493D14
00052 */
00053     if(*z <= 0.0e0) ftnstop("Zero or negative argument in DSTREM");
00054     if(!(*z > 6.0e0)) goto S10;
00055     T2 = 1.0e0/pow(*z,2.0);
00056     dstrem = devlpl(coef,&K1,&T2)**z;
00057     goto S20;
00058 S10:
00059     sterl = hln2pi+(*z-0.5e0)*log(*z)-*z;
00060     dstrem = dlngam(z)-sterl;
00061 S20:
00062     return dstrem;
00063 #undef hln2pi
00064 #undef ncoef
00065 } /* END */

00046 {
00047     E0001(1,NULL,NULL,NULL,NULL,NULL,NULL,NULL,zabstl,zreltl,zxhi,zxlo);
00048 } /* END */

00034 {
00035 static double coef[4][5] = {
00036     1.0e0,1.0e0,0.0e0,0.0e0,0.0e0,3.0e0,16.0e0,5.0e0,0.0e0,0.0e0,-15.0e0,17.0e0,
00037     19.0e0,3.0e0,0.0e0,-945.0e0,-1920.0e0,1482.0e0,776.0e0,79.0e0
00038 };
00039 static double denom[4] = {
00040     4.0e0,96.0e0,384.0e0,92160.0e0
00041 };
00042 static int ideg[4] = {
00043     2,3,4,5
00044 };
00045 static double dt1,denpow,sum,term,x,xp,xx;
00046 static int i;
00047 /*
00048      ..
00049      .. Executable Statements ..
00050 */
00051     x = fabs(dinvnr(p,q));
00052     xx = x*x;
00053     sum = x;
00054     denpow = 1.0e0;
00055     for(i=0; i<4; i++) {
00056         term = devlpl(&coef[i][0],&ideg[i],&xx)*x;
00057         denpow *= *df;
00058         sum += (term/(denpow*denom[i]));
00059     }
00060     if(!(*p >= 0.5e0)) goto S20;
00061     xp = sum;
00062     goto S30;
00063 S20:
00064     xp = -sum;
00065 S30:
00066     dt1 = xp;
00067     return dt1;
00068 } /* END */