+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include "cdflib.h"
-#include <assert.h>
-
-static void E0000(int,int*,double*,double*,unsigned long*,
- unsigned long*,double*,double*,double*,
- double*,double*,double*,double*);
-static void E0001(int,int*,double*,double*,double*,double*,
- unsigned long*,unsigned long*,double*,double*,
- double*,double*);
-
-/*
- * A comment about ints and longs - whether ints or longs are used should
- * make no difference, but where double r-values are assigned to ints the
- * r-value is cast converted to a long, which is then assigned to the int
- * to be compatible with the operation of fifidint.
- */
-/*
------------------------------------------------------------------------
-
- COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
-
- --------
-
- IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
- LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
-
------------------------------------------------------------------------
-*/
-double algdiv(double *a,double *b)
-{
-static double c0 = .833333333333333e-01;
-static double c1 = -.277777777760991e-02;
-static double c2 = .793650666825390e-03;
-static double c3 = -.595202931351870e-03;
-static double c4 = .837308034031215e-03;
-static double c5 = -.165322962780713e-02;
-static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*a <= *b) goto S10;
- h = *b/ *a;
- c = 1.0e0/(1.0e0+h);
- x = h/(1.0e0+h);
- d = *a+(*b-0.5e0);
- goto S20;
-S10:
- h = *a/ *b;
- c = h/(1.0e0+h);
- x = 1.0e0/(1.0e0+h);
- d = *b+(*a-0.5e0);
-S20:
-/*
- SET SN = (1 - X**N)/(1 - X)
-*/
- x2 = x*x;
- s3 = 1.0e0+(x+x2);
- s5 = 1.0e0+(x+x2*s3);
- s7 = 1.0e0+(x+x2*s5);
- s9 = 1.0e0+(x+x2*s7);
- s11 = 1.0e0+(x+x2*s9);
-/*
- SET W = DEL(B) - DEL(A + B)
-*/
- t = pow(1.0e0/ *b,2.0);
- w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
- w *= (c/ *b);
-/*
- COMBINE THE RESULTS
-*/
- T1 = *a/ *b;
- u = d*alnrel(&T1);
- v = *a*(log(*b)-1.0e0);
- if(u <= v) goto S30;
- algdiv = w-v-u;
- return algdiv;
-S30:
- algdiv = w-u-v;
- return algdiv;
-}
-double alngam(double *x)
-/*
-**********************************************************************
-
- double alngam(double *x)
- double precision LN of the GAMma function
-
-
- Function
-
-
- Returns the natural logarithm of GAMMA(X).
-
-
- Arguments
-
-
- X --> value at which scaled log gamma is to be returned
- X is DOUBLE PRECISION
-
-
- Method
-
-
- If X .le. 6.0, then use recursion to get X below 3
- then apply rational approximation number 5236 of
- Hart et al, Computer Approximations, John Wiley and
- Sons, NY, 1968.
-
- If X .gt. 6.0, then use recursion to get X to at least 12 and
- then use formula 5423 of the same source.
-
-**********************************************************************
-*/
-{
-#define hln2pi 0.91893853320467274178e0
-static double coef[5] = {
- 0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
- -0.594997310889e-3,0.8065880899e-3
-};
-static double scoefd[4] = {
- 0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
- 0.1000000000000000000e1
-};
-static double scoefn[9] = {
- 0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
- 0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
- 0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
-};
-static int K1 = 9;
-static int K3 = 4;
-static int K5 = 5;
-static double alngam,offset,prod,xx;
-static int i,n;
-static double T2,T4,T6;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*x <= 6.0e0)) goto S70;
- prod = 1.0e0;
- xx = *x;
- if(!(*x > 3.0e0)) goto S30;
-S10:
- if(!(xx > 3.0e0)) goto S20;
- xx -= 1.0e0;
- prod *= xx;
- goto S10;
-S30:
-S20:
- if(!(*x < 2.0e0)) goto S60;
-S40:
- if(!(xx < 2.0e0)) goto S50;
- prod /= xx;
- xx += 1.0e0;
- goto S40;
-S60:
-S50:
- T2 = xx-2.0e0;
- T4 = xx-2.0e0;
- alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
-/*
- COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
-*/
- alngam *= prod;
- alngam = log(alngam);
- goto S110;
-S70:
- offset = hln2pi;
-/*
- IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
-*/
- n = fifidint(12.0e0-*x);
- if(!(n > 0)) goto S90;
- prod = 1.0e0;
- for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));
- offset -= log(prod);
- xx = *x+(double)n;
- goto S100;
-S90:
- xx = *x;
-S100:
-/*
- COMPUTE POWER SERIES
-*/
- T6 = 1.0e0/pow(xx,2.0);
- alngam = devlpl(coef,&K5,&T6)/xx;
- alngam += (offset+(xx-0.5e0)*log(xx)-xx);
-S110:
- return alngam;
-#undef hln2pi
-}
-double alnrel(double *a)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE FUNCTION LN(1 + A)
------------------------------------------------------------------------
-*/
-{
-static double p1 = -.129418923021993e+01;
-static double p2 = .405303492862024e+00;
-static double p3 = -.178874546012214e-01;
-static double q1 = -.162752256355323e+01;
-static double q2 = .747811014037616e+00;
-static double q3 = -.845104217945565e-01;
-static double alnrel,t,t2,w,x;
-/*
- ..
- .. Executable Statements ..
-*/
- if(fabs(*a) > 0.375e0) goto S10;
- t = *a/(*a+2.0e0);
- t2 = t*t;
- w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);
- alnrel = 2.0e0*t*w;
- return alnrel;
-S10:
- x = 1.e0+*a;
- alnrel = log(x);
- return alnrel;
-}
-double apser(double *a,double *b,double *x,double *eps)
-/*
------------------------------------------------------------------------
- APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
- A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
- A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
------------------------------------------------------------------------
-*/
-{
-static double g = .577215664901533e0;
-static double apser,aj,bx,c,j,s,t,tol;
-/*
- ..
- .. Executable Statements ..
-*/
- bx = *b**x;
- t = *x-bx;
- if(*b**eps > 2.e-2) goto S10;
- c = log(*x)+psi(b)+g+t;
- goto S20;
-S10:
- c = log(bx)+g+t;
-S20:
- tol = 5.0e0**eps*fabs(c);
- j = 1.0e0;
- s = 0.0e0;
-S30:
- j += 1.0e0;
- t *= (*x-bx/j);
- aj = t/j;
- s += aj;
- if(fabs(aj) > tol) goto S30;
- apser = -(*a*(c+s));
- return apser;
-}
-double basym(double *a,double *b,double *lambda,double *eps)
-/*
------------------------------------------------------------------------
- ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
- LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED.
- IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
- A AND B ARE GREATER THAN OR EQUAL TO 15.
------------------------------------------------------------------------
-*/
-{
-static double e0 = 1.12837916709551e0;
-static double e1 = .353553390593274e0;
-static int num = 20;
-/*
-------------------------
- ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
- ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
- THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
-------------------------
- E0 = 2/SQRT(PI)
- E1 = 2**(-3/2)
-------------------------
-*/
-static int K3 = 1;
-static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
- z2,zn,znm1;
-static int i,im1,imj,j,m,mm1,mmj,n,np1;
-static double a0[21],b0[21],c[21],d[21],T1,T2;
-/*
- ..
- .. Executable Statements ..
-*/
- basym = 0.0e0;
- if(*a >= *b) goto S10;
- h = *a/ *b;
- r0 = 1.0e0/(1.0e0+h);
- r1 = (*b-*a)/ *b;
- w0 = 1.0e0/sqrt(*a*(1.0e0+h));
- goto S20;
-S10:
- h = *b/ *a;
- r0 = 1.0e0/(1.0e0+h);
- r1 = (*b-*a)/ *a;
- w0 = 1.0e0/sqrt(*b*(1.0e0+h));
-S20:
- T1 = -(*lambda/ *a);
- T2 = *lambda/ *b;
- f = *a*rlog1(&T1)+*b*rlog1(&T2);
- t = exp(-f);
- if(t == 0.0e0) return basym;
- z0 = sqrt(f);
- z = 0.5e0*(z0/e1);
- z2 = f+f;
- a0[0] = 2.0e0/3.0e0*r1;
- c[0] = -(0.5e0*a0[0]);
- d[0] = -c[0];
- j0 = 0.5e0/e0*erfc1(&K3,&z0);
- j1 = e1;
- sum = j0+d[0]*w0*j1;
- s = 1.0e0;
- h2 = h*h;
- hn = 1.0e0;
- w = w0;
- znm1 = z;
- zn = z2;
- for(n=2; n<=num; n+=2) {
- hn = h2*hn;
- a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
- np1 = n+1;
- s += hn;
- a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
- for(i=n; i<=np1; i++) {
- r = -(0.5e0*((double)i+1.0e0));
- b0[0] = r*a0[0];
- for(m=2; m<=i; m++) {
- bsum = 0.0e0;
- mm1 = m-1;
- for(j=1; j<=mm1; j++) {
- mmj = m-j;
- bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
- }
- b0[m-1] = r*a0[m-1]+bsum/(double)m;
- }
- c[i-1] = b0[i-1]/((double)i+1.0e0);
- dsum = 0.0e0;
- im1 = i-1;
- for(j=1; j<=im1; j++) {
- imj = i-j;
- dsum += (d[imj-1]*c[j-1]);
- }
- d[i-1] = -(dsum+c[i-1]);
- }
- j0 = e1*znm1+((double)n-1.0e0)*j0;
- j1 = e1*zn+(double)n*j1;
- znm1 = z2*znm1;
- zn = z2*zn;
- w = w0*w;
- t0 = d[n-1]*w*j0;
- w = w0*w;
- t1 = d[np1-1]*w*j1;
- sum += (t0+t1);
- if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
- }
-S80:
- u = exp(-bcorr(a,b));
- basym = e0*t*u*sum;
- return basym;
-}
-double bcorr(double *a0,double *b0)
-/*
------------------------------------------------------------------------
-
- EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE
- LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
- IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8.
-
------------------------------------------------------------------------
-*/
-{
-static double c0 = .833333333333333e-01;
-static double c1 = -.277777777760991e-02;
-static double c2 = .793650666825390e-03;
-static double c3 = -.595202931351870e-03;
-static double c4 = .837308034031215e-03;
-static double c5 = -.165322962780713e-02;
-static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;
-/*
- ..
- .. Executable Statements ..
-*/
- a = fifdmin1(*a0,*b0);
- b = fifdmax1(*a0,*b0);
- h = a/b;
- c = h/(1.0e0+h);
- x = 1.0e0/(1.0e0+h);
- x2 = x*x;
-/*
- SET SN = (1 - X**N)/(1 - X)
-*/
- s3 = 1.0e0+(x+x2);
- s5 = 1.0e0+(x+x2*s3);
- s7 = 1.0e0+(x+x2*s5);
- s9 = 1.0e0+(x+x2*s7);
- s11 = 1.0e0+(x+x2*s9);
-/*
- SET W = DEL(B) - DEL(A + B)
-*/
- t = pow(1.0e0/b,2.0);
- w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
- w *= (c/b);
-/*
- COMPUTE DEL(A) + W
-*/
- t = pow(1.0e0/a,2.0);
- bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w;
- return bcorr;
-}
-double betaln(double *a0,double *b0)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
------------------------------------------------------------------------
- E = 0.5*LN(2*PI)
---------------------------
-*/
-{
-static double e = .918938533204673e0;
-static double betaln,a,b,c,h,u,v,w,z;
-static int i,n;
-static double T1;
-/*
- ..
- .. Executable Statements ..
-*/
- a = fifdmin1(*a0,*b0);
- b = fifdmax1(*a0,*b0);
- if(a >= 8.0e0) goto S100;
- if(a >= 1.0e0) goto S20;
-/*
------------------------------------------------------------------------
- PROCEDURE WHEN A .LT. 1
------------------------------------------------------------------------
-*/
- if(b >= 8.0e0) goto S10;
- T1 = a+b;
- betaln = gamln(&a)+(gamln(&b)-gamln(&T1));
- return betaln;
-S10:
- betaln = gamln(&a)+algdiv(&a,&b);
- return betaln;
-S20:
-/*
------------------------------------------------------------------------
- PROCEDURE WHEN 1 .LE. A .LT. 8
------------------------------------------------------------------------
-*/
- if(a > 2.0e0) goto S40;
- if(b > 2.0e0) goto S30;
- betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b);
- return betaln;
-S30:
- w = 0.0e0;
- if(b < 8.0e0) goto S60;
- betaln = gamln(&a)+algdiv(&a,&b);
- return betaln;
-S40:
-/*
- REDUCTION OF A WHEN B .LE. 1000
-*/
- if(b > 1000.0e0) goto S80;
- n = (long)(a - 1.0e0);
- w = 1.0e0;
- for(i=1; i<=n; i++) {
- a -= 1.0e0;
- h = a/b;
- w *= (h/(1.0e0+h));
- }
- w = log(w);
- if(b < 8.0e0) goto S60;
- betaln = w+gamln(&a)+algdiv(&a,&b);
- return betaln;
-S60:
-/*
- REDUCTION OF B WHEN B .LT. 8
-*/
- n = (long)(b - 1.0e0);
- z = 1.0e0;
- for(i=1; i<=n; i++) {
- b -= 1.0e0;
- z *= (b/(a+b));
- }
- betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
- return betaln;
-S80:
-/*
- REDUCTION OF A WHEN B .GT. 1000
-*/
- n = (long)(a - 1.0e0);
- w = 1.0e0;
- for(i=1; i<=n; i++) {
- a -= 1.0e0;
- w *= (a/(1.0e0+a/b));
- }
- betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
- return betaln;
-S100:
-/*
------------------------------------------------------------------------
- PROCEDURE WHEN A .GE. 8
------------------------------------------------------------------------
-*/
- w = bcorr(&a,&b);
- h = a/b;
- c = h/(1.0e0+h);
- u = -((a-0.5e0)*log(c));
- v = b*alnrel(&h);
- if(u <= v) goto S110;
- betaln = -(0.5e0*log(b))+e+w-v-u;
- return betaln;
-S110:
- betaln = -(0.5e0*log(b))+e+w-u-v;
- return betaln;
-}
-double bfrac(double *a,double *b,double *x,double *y,double *lambda,
- double *eps)
-/*
------------------------------------------------------------------------
- CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
- IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B.
------------------------------------------------------------------------
-*/
-{
-static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;
-/*
- ..
- .. Executable Statements ..
-*/
- bfrac = brcomp(a,b,x,y);
- if(bfrac == 0.0e0) return bfrac;
- c = 1.0e0+*lambda;
- c0 = *b/ *a;
- c1 = 1.0e0+1.0e0/ *a;
- yp1 = *y+1.0e0;
- n = 0.0e0;
- p = 1.0e0;
- s = *a+1.0e0;
- an = 0.0e0;
- bn = anp1 = 1.0e0;
- bnp1 = c/c1;
- r = c1/c;
-S10:
-/*
- CONTINUED FRACTION CALCULATION
-*/
- n += 1.0e0;
- t = n/ *a;
- w = n*(*b-n)**x;
- e = *a/s;
- alpha = p*(p+c0)*e*e*(w**x);
- e = (1.0e0+t)/(c1+t+t);
- beta = n+w/s+e*(c+n*yp1);
- p = 1.0e0+t;
- s += 2.0e0;
-/*
- UPDATE AN, BN, ANP1, AND BNP1
-*/
- t = alpha*an+beta*anp1;
- an = anp1;
- anp1 = t;
- t = alpha*bn+beta*bnp1;
- bn = bnp1;
- bnp1 = t;
- r0 = r;
- r = anp1/bnp1;
- if(fabs(r-r0) <= *eps*r) goto S20;
-/*
- RESCALE AN, BN, ANP1, AND BNP1
-*/
- an /= bnp1;
- bn /= bnp1;
- anp1 = r;
- bnp1 = 1.0e0;
- goto S10;
-S20:
-/*
- TERMINATION
-*/
- bfrac *= r;
- return bfrac;
-}
-void bgrat(double *a,double *b,double *x,double *y,double *w,
- double *eps,int *ierr)
-/*
------------------------------------------------------------------------
- ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
- THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
- THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED.
- IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
------------------------------------------------------------------------
-*/
-{
-static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
-static int i,n,nm1;
-static double c[30],d[30],T1;
-/*
- ..
- .. Executable Statements ..
-*/
- bm1 = *b-0.5e0-0.5e0;
- nu = *a+0.5e0*bm1;
- if(*y > 0.375e0) goto S10;
- T1 = -*y;
- lnx = alnrel(&T1);
- goto S20;
-S10:
- lnx = log(*x);
-S20:
- z = -(nu*lnx);
- if(*b*z == 0.0e0) goto S70;
-/*
- COMPUTATION OF THE EXPANSION
- SET R = EXP(-Z)*Z**B/GAMMA(B)
-*/
- r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
- r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
- u = algdiv(b,a)+*b*log(nu);
- u = r*exp(-u);
- if(u == 0.0e0) goto S70;
- grat1(b,&z,&r,&p,&q,eps);
- v = 0.25e0*pow(1.0e0/nu,2.0);
- t2 = 0.25e0*lnx*lnx;
- l = *w/u;
- j = q/r;
- sum = j;
- t = cn = 1.0e0;
- n2 = 0.0e0;
- for(n=1; n<=30; n++) {
- bp2n = *b+n2;
- j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
- n2 += 2.0e0;
- t *= t2;
- cn /= (n2*(n2+1.0e0));
- c[n-1] = cn;
- s = 0.0e0;
- if(n == 1) goto S40;
- nm1 = n-1;
- coef = *b-(double)n;
- for(i=1; i<=nm1; i++) {
- s += (coef*c[i-1]*d[n-i-1]);
- coef += *b;
- }
-S40:
- d[n-1] = bm1*cn+s/(double)n;
- dj = d[n-1]*j;
- sum += dj;
- if(sum <= 0.0e0) goto S70;
- if(fabs(dj) <= *eps*(sum+l)) goto S60;
- }
-S60:
-/*
- ADD THE RESULTS TO W
-*/
- *ierr = 0;
- *w += (u*sum);
- return;
-S70:
-/*
- THE EXPANSION CANNOT BE COMPUTED
-*/
- *ierr = 1;
- return;
-}
-double bpser(double *a,double *b,double *x,double *eps)
-/*
------------------------------------------------------------------------
- POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
- OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED.
------------------------------------------------------------------------
-*/
-{
-static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
-static int i,m;
-/*
- ..
- .. Executable Statements ..
-*/
- bpser = 0.0e0;
- if(*x == 0.0e0) return bpser;
-/*
------------------------------------------------------------------------
- COMPUTE THE FACTOR X**A/(A*BETA(A,B))
------------------------------------------------------------------------
-*/
- a0 = fifdmin1(*a,*b);
- if(a0 < 1.0e0) goto S10;
- z = *a*log(*x)-betaln(a,b);
- bpser = exp(z)/ *a;
- goto S100;
-S10:
- b0 = fifdmax1(*a,*b);
- if(b0 >= 8.0e0) goto S90;
- if(b0 > 1.0e0) goto S40;
-/*
- PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
-*/
- bpser = pow(*x,*a);
- if(bpser == 0.0e0) return bpser;
- apb = *a+*b;
- if(apb > 1.0e0) goto S20;
- z = 1.0e0+gam1(&apb);
- goto S30;
-S20:
- u = *a+*b-1.e0;
- z = (1.0e0+gam1(&u))/apb;
-S30:
- c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
- bpser *= (c*(*b/apb));
- goto S100;
-S40:
-/*
- PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
-*/
- u = gamln1(&a0);
- m = (long)(b0 - 1.0e0);
- if(m < 1) goto S60;
- c = 1.0e0;
- for(i=1; i<=m; i++) {
- b0 -= 1.0e0;
- c *= (b0/(a0+b0));
- }
- u = log(c)+u;
-S60:
- z = *a*log(*x)-u;
- b0 -= 1.0e0;
- apb = a0+b0;
- if(apb > 1.0e0) goto S70;
- t = 1.0e0+gam1(&apb);
- goto S80;
-S70:
- u = a0+b0-1.e0;
- t = (1.0e0+gam1(&u))/apb;
-S80:
- bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
- goto S100;
-S90:
-/*
- PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
-*/
- u = gamln1(&a0)+algdiv(&a0,&b0);
- z = *a*log(*x)-u;
- bpser = a0/ *a*exp(z);
-S100:
- if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
-/*
------------------------------------------------------------------------
- COMPUTE THE SERIES
------------------------------------------------------------------------
-*/
- sum = n = 0.0e0;
- c = 1.0e0;
- tol = *eps/ *a;
-S110:
- n += 1.0e0;
- c *= ((0.5e0+(0.5e0-*b/n))**x);
- w = c/(*a+n);
- sum += w;
- if(fabs(w) > tol) goto S110;
- bpser *= (1.0e0+*a*sum);
- return bpser;
-}
-void bratio(double *a,double *b,double *x,double *y,double *w,
- double *w1,int *ierr)
-/*
------------------------------------------------------------------------
-
- EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
-
- --------------------
-
- IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
- AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES
-
- W = IX(A,B)
- W1 = 1 - IX(A,B)
-
- IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
- IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
- W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
- THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
- ONE OF THE FOLLOWING VALUES ...
-
- IERR = 1 IF A OR B IS NEGATIVE
- IERR = 2 IF A = B = 0
- IERR = 3 IF X .LT. 0 OR X .GT. 1
- IERR = 4 IF Y .LT. 0 OR Y .GT. 1
- IERR = 5 IF X + Y .NE. 1
- IERR = 6 IF X = A = 0
- IERR = 7 IF Y = B = 0
-
---------------------
- WRITTEN BY ALFRED H. MORRIS, JR.
- NAVAL SURFACE WARFARE CENTER
- DAHLGREN, VIRGINIA
- REVISED ... NOV 1991
------------------------------------------------------------------------
-*/
-{
-static int K1 = 1;
-static double a0,b0,eps,lambda,t,x0,y0,z;
-static int ierr1,ind,n;
-static double T2,T3,T4,T5;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
- FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
-*/
- eps = spmpar(&K1);
- *w = *w1 = 0.0e0;
- if(*a < 0.0e0 || *b < 0.0e0) goto S270;
- if(*a == 0.0e0 && *b == 0.0e0) goto S280;
- if(*x < 0.0e0 || *x > 1.0e0) goto S290;
- if(*y < 0.0e0 || *y > 1.0e0) goto S300;
- z = *x+*y-0.5e0-0.5e0;
- if(fabs(z) > 3.0e0*eps) goto S310;
- *ierr = 0;
- if(*x == 0.0e0) goto S210;
- if(*y == 0.0e0) goto S230;
- if(*a == 0.0e0) goto S240;
- if(*b == 0.0e0) goto S220;
- eps = fifdmax1(eps,1.e-15);
- if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
- ind = 0;
- a0 = *a;
- b0 = *b;
- x0 = *x;
- y0 = *y;
- if(fifdmin1(a0,b0) > 1.0e0) goto S40;
-/*
- PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
-*/
- if(*x <= 0.5e0) goto S10;
- ind = 1;
- a0 = *b;
- b0 = *a;
- x0 = *y;
- y0 = *x;
-S10:
- if(b0 < fifdmin1(eps,eps*a0)) goto S90;
- if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
- if(fifdmax1(a0,b0) > 1.0e0) goto S20;
- if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
- if(pow(x0,a0) <= 0.9e0) goto S110;
- if(x0 >= 0.3e0) goto S120;
- n = 20;
- goto S140;
-S20:
- if(b0 <= 1.0e0) goto S110;
- if(x0 >= 0.3e0) goto S120;
- if(x0 >= 0.1e0) goto S30;
- if(pow(x0*b0,a0) <= 0.7e0) goto S110;
-S30:
- if(b0 > 15.0e0) goto S150;
- n = 20;
- goto S140;
-S40:
-/*
- PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
-*/
- if(*a > *b) goto S50;
- lambda = *a-(*a+*b)**x;
- goto S60;
-S50:
- lambda = (*a+*b)**y-*b;
-S60:
- if(lambda >= 0.0e0) goto S70;
- ind = 1;
- a0 = *b;
- b0 = *a;
- x0 = *y;
- y0 = *x;
- lambda = fabs(lambda);
-S70:
- if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
- if(b0 < 40.0e0) goto S160;
- if(a0 > b0) goto S80;
- if(a0 <= 100.0e0) goto S130;
- if(lambda > 0.03e0*a0) goto S130;
- goto S200;
-S80:
- if(b0 <= 100.0e0) goto S130;
- if(lambda > 0.03e0*b0) goto S130;
- goto S200;
-S90:
-/*
- EVALUATION OF THE APPROPRIATE ALGORITHM
-*/
- *w = fpser(&a0,&b0,&x0,&eps);
- *w1 = 0.5e0+(0.5e0-*w);
- goto S250;
-S100:
- *w1 = apser(&a0,&b0,&x0,&eps);
- *w = 0.5e0+(0.5e0-*w1);
- goto S250;
-S110:
- *w = bpser(&a0,&b0,&x0,&eps);
- *w1 = 0.5e0+(0.5e0-*w);
- goto S250;
-S120:
- *w1 = bpser(&b0,&a0,&y0,&eps);
- *w = 0.5e0+(0.5e0-*w1);
- goto S250;
-S130:
- T2 = 15.0e0*eps;
- *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
- *w1 = 0.5e0+(0.5e0-*w);
- goto S250;
-S140:
- *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
- b0 += (double)n;
-S150:
- T3 = 15.0e0*eps;
- bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
- *w = 0.5e0+(0.5e0-*w1);
- goto S250;
-S160:
- n = (long)(b0);
- b0 -= (double)n;
- if(b0 != 0.0e0) goto S170;
- n -= 1;
- b0 = 1.0e0;
-S170:
- *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
- if(x0 > 0.7e0) goto S180;
- *w += bpser(&a0,&b0,&x0,&eps);
- *w1 = 0.5e0+(0.5e0-*w);
- goto S250;
-S180:
- if(a0 > 15.0e0) goto S190;
- n = 20;
- *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
- a0 += (double)n;
-S190:
- T4 = 15.0e0*eps;
- bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
- *w1 = 0.5e0+(0.5e0-*w);
- goto S250;
-S200:
- T5 = 100.0e0*eps;
- *w = basym(&a0,&b0,&lambda,&T5);
- *w1 = 0.5e0+(0.5e0-*w);
- goto S250;
-S210:
-/*
- TERMINATION OF THE PROCEDURE
-*/
- if(*a == 0.0e0) goto S320;
-S220:
- *w = 0.0e0;
- *w1 = 1.0e0;
- return;
-S230:
- if(*b == 0.0e0) goto S330;
-S240:
- *w = 1.0e0;
- *w1 = 0.0e0;
- return;
-S250:
- if(ind == 0) return;
- t = *w;
- *w = *w1;
- *w1 = t;
- return;
-S260:
-/*
- PROCEDURE FOR A AND B .LT. 1.E-3*EPS
-*/
- *w = *b/(*a+*b);
- *w1 = *a/(*a+*b);
- return;
-S270:
-/*
- ERROR RETURN
-*/
- *ierr = 1;
- return;
-S280:
- *ierr = 2;
- return;
-S290:
- *ierr = 3;
- return;
-S300:
- *ierr = 4;
- return;
-S310:
- *ierr = 5;
- return;
-S320:
- *ierr = 6;
- return;
-S330:
- *ierr = 7;
- return;
-}
-double brcmp1(int *mu,double *a,double *b,double *x,double *y)
-/*
------------------------------------------------------------------------
- EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B))
------------------------------------------------------------------------
-*/
-{
-static double Const = .398942280401433e0;
-static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
-static int i,n;
-/*
------------------
- CONST = 1/SQRT(2*PI)
------------------
-*/
-static double T1,T2,T3,T4;
-/*
- ..
- .. Executable Statements ..
-*/
- a0 = fifdmin1(*a,*b);
- if(a0 >= 8.0e0) goto S130;
- if(*x > 0.375e0) goto S10;
- lnx = log(*x);
- T1 = -*x;
- lny = alnrel(&T1);
- goto S30;
-S10:
- if(*y > 0.375e0) goto S20;
- T2 = -*y;
- lnx = alnrel(&T2);
- lny = log(*y);
- goto S30;
-S20:
- lnx = log(*x);
- lny = log(*y);
-S30:
- z = *a*lnx+*b*lny;
- if(a0 < 1.0e0) goto S40;
- z -= betaln(a,b);
- brcmp1 = esum(mu,&z);
- return brcmp1;
-S40:
-/*
------------------------------------------------------------------------
- PROCEDURE FOR A .LT. 1 OR B .LT. 1
------------------------------------------------------------------------
-*/
- b0 = fifdmax1(*a,*b);
- if(b0 >= 8.0e0) goto S120;
- if(b0 > 1.0e0) goto S70;
-/*
- ALGORITHM FOR B0 .LE. 1
-*/
- brcmp1 = esum(mu,&z);
- if(brcmp1 == 0.0e0) return brcmp1;
- apb = *a+*b;
- if(apb > 1.0e0) goto S50;
- z = 1.0e0+gam1(&apb);
- goto S60;
-S50:
- u = *a+*b-1.e0;
- z = (1.0e0+gam1(&u))/apb;
-S60:
- c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
- brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
- return brcmp1;
-S70:
-/*
- ALGORITHM FOR 1 .LT. B0 .LT. 8
-*/
- u = gamln1(&a0);
- n = (long)(b0 - 1.0e0);
- if(n < 1) goto S90;
- c = 1.0e0;
- for(i=1; i<=n; i++) {
- b0 -= 1.0e0;
- c *= (b0/(a0+b0));
- }
- u = log(c)+u;
-S90:
- z -= u;
- b0 -= 1.0e0;
- apb = a0+b0;
- if(apb > 1.0e0) goto S100;
- t = 1.0e0+gam1(&apb);
- goto S110;
-S100:
- u = a0+b0-1.e0;
- t = (1.0e0+gam1(&u))/apb;
-S110:
- brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
- return brcmp1;
-S120:
-/*
- ALGORITHM FOR B0 .GE. 8
-*/
- u = gamln1(&a0)+algdiv(&a0,&b0);
- T3 = z-u;
- brcmp1 = a0*esum(mu,&T3);
- return brcmp1;
-S130:
-/*
------------------------------------------------------------------------
- PROCEDURE FOR A .GE. 8 AND B .GE. 8
------------------------------------------------------------------------
-*/
- if(*a > *b) goto S140;
- h = *a/ *b;
- x0 = h/(1.0e0+h);
- y0 = 1.0e0/(1.0e0+h);
- lambda = *a-(*a+*b)**x;
- goto S150;
-S140:
- h = *b/ *a;
- x0 = 1.0e0/(1.0e0+h);
- y0 = h/(1.0e0+h);
- lambda = (*a+*b)**y-*b;
-S150:
- e = -(lambda/ *a);
- if(fabs(e) > 0.6e0) goto S160;
- u = rlog1(&e);
- goto S170;
-S160:
- u = e-log(*x/x0);
-S170:
- e = lambda/ *b;
- if(fabs(e) > 0.6e0) goto S180;
- v = rlog1(&e);
- goto S190;
-S180:
- v = e-log(*y/y0);
-S190:
- T4 = -(*a*u+*b*v);
- z = esum(mu,&T4);
- brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
- return brcmp1;
-}
-double brcomp(double *a,double *b,double *x,double *y)
-/*
------------------------------------------------------------------------
- EVALUATION OF X**A*Y**B/BETA(A,B)
------------------------------------------------------------------------
-*/
-{
-static double Const = .398942280401433e0;
-static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
-static int i,n;
-/*
------------------
- CONST = 1/SQRT(2*PI)
------------------
-*/
-static double T1,T2;
-/*
- ..
- .. Executable Statements ..
-*/
- brcomp = 0.0e0;
- if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
- a0 = fifdmin1(*a,*b);
- if(a0 >= 8.0e0) goto S130;
- if(*x > 0.375e0) goto S10;
- lnx = log(*x);
- T1 = -*x;
- lny = alnrel(&T1);
- goto S30;
-S10:
- if(*y > 0.375e0) goto S20;
- T2 = -*y;
- lnx = alnrel(&T2);
- lny = log(*y);
- goto S30;
-S20:
- lnx = log(*x);
- lny = log(*y);
-S30:
- z = *a*lnx+*b*lny;
- if(a0 < 1.0e0) goto S40;
- z -= betaln(a,b);
- brcomp = exp(z);
- return brcomp;
-S40:
-/*
------------------------------------------------------------------------
- PROCEDURE FOR A .LT. 1 OR B .LT. 1
------------------------------------------------------------------------
-*/
- b0 = fifdmax1(*a,*b);
- if(b0 >= 8.0e0) goto S120;
- if(b0 > 1.0e0) goto S70;
-/*
- ALGORITHM FOR B0 .LE. 1
-*/
- brcomp = exp(z);
- if(brcomp == 0.0e0) return brcomp;
- apb = *a+*b;
- if(apb > 1.0e0) goto S50;
- z = 1.0e0+gam1(&apb);
- goto S60;
-S50:
- u = *a+*b-1.e0;
- z = (1.0e0+gam1(&u))/apb;
-S60:
- c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
- brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
- return brcomp;
-S70:
-/*
- ALGORITHM FOR 1 .LT. B0 .LT. 8
-*/
- u = gamln1(&a0);
- n = (long)(b0 - 1.0e0);
- if(n < 1) goto S90;
- c = 1.0e0;
- for(i=1; i<=n; i++) {
- b0 -= 1.0e0;
- c *= (b0/(a0+b0));
- }
- u = log(c)+u;
-S90:
- z -= u;
- b0 -= 1.0e0;
- apb = a0+b0;
- if(apb > 1.0e0) goto S100;
- t = 1.0e0+gam1(&apb);
- goto S110;
-S100:
- u = a0+b0-1.e0;
- t = (1.0e0+gam1(&u))/apb;
-S110:
- brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
- return brcomp;
-S120:
-/*
- ALGORITHM FOR B0 .GE. 8
-*/
- u = gamln1(&a0)+algdiv(&a0,&b0);
- brcomp = a0*exp(z-u);
- return brcomp;
-S130:
-/*
------------------------------------------------------------------------
- PROCEDURE FOR A .GE. 8 AND B .GE. 8
------------------------------------------------------------------------
-*/
- if(*a > *b) goto S140;
- h = *a/ *b;
- x0 = h/(1.0e0+h);
- y0 = 1.0e0/(1.0e0+h);
- lambda = *a-(*a+*b)**x;
- goto S150;
-S140:
- h = *b/ *a;
- x0 = 1.0e0/(1.0e0+h);
- y0 = h/(1.0e0+h);
- lambda = (*a+*b)**y-*b;
-S150:
- e = -(lambda/ *a);
- if(fabs(e) > 0.6e0) goto S160;
- u = rlog1(&e);
- goto S170;
-S160:
- u = e-log(*x/x0);
-S170:
- e = lambda/ *b;
- if(fabs(e) > 0.6e0) goto S180;
- v = rlog1(&e);
- goto S190;
-S180:
- v = e-log(*y/y0);
-S190:
- z = exp(-(*a*u+*b*v));
- brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
- return brcomp;
-}
-double bup(double *a,double *b,double *x,double *y,int *n,double *eps)
-/*
------------------------------------------------------------------------
- EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
- EPS IS THE TOLERANCE USED.
------------------------------------------------------------------------
-*/
-{
-static int K1 = 1;
-static int K2 = 0;
-static double bup,ap1,apb,d,l,r,t,w;
-static int i,k,kp1,mu,nm1;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- OBTAIN THE SCALING FACTOR EXP(-MU) AND
- EXP(MU)*(X**A*Y**B/BETA(A,B))/A
-*/
- apb = *a+*b;
- ap1 = *a+1.0e0;
- mu = 0;
- d = 1.0e0;
- if(*n == 1 || *a < 1.0e0) goto S10;
- if(apb < 1.1e0*ap1) goto S10;
- mu = (long)(fabs(exparg(&K1)));
- k = (long)(exparg(&K2));
- if(k < mu) mu = k;
- t = mu;
- d = exp(-t);
-S10:
- bup = brcmp1(&mu,a,b,x,y)/ *a;
- if(*n == 1 || bup == 0.0e0) return bup;
- nm1 = *n-1;
- w = d;
-/*
- LET K BE THE INDEX OF THE MAXIMUM TERM
-*/
- k = 0;
- if(*b <= 1.0e0) goto S50;
- if(*y > 1.e-4) goto S20;
- k = nm1;
- goto S30;
-S20:
- r = (*b-1.0e0)**x/ *y-*a;
- if(r < 1.0e0) goto S50;
- t = nm1;
- k = (long)(t);
- if(r < t) k = (long)(r);
-S30:
-/*
- ADD THE INCREASING TERMS OF THE SERIES
-*/
- for(i=1; i<=k; i++) {
- l = i-1;
- d = (apb+l)/(ap1+l)**x*d;
- w += d;
- }
- if(k == nm1) goto S70;
-S50:
-/*
- ADD THE REMAINING TERMS OF THE SERIES
-*/
- kp1 = k+1;
- for(i=kp1; i<=nm1; i++) {
- l = i-1;
- d = (apb+l)/(ap1+l)**x*d;
- w += d;
- if(d <= *eps*w) goto S70;
- }
-S70:
-/*
- TERMINATE THE PROCEDURE
-*/
- bup *= w;
- return bup;
-}
-void cdfbet(int *which,double *p,double *q,double *x,double *y,
- double *a,double *b,int *status,double *bound)
-/**********************************************************************
-
- void cdfbet(int *which,double *p,double *q,double *x,double *y,
- double *a,double *b,int *status,double *bound)
-
- Cumulative Distribution Function
- BETa Distribution
-
-
- Function
-
-
- Calculates any one parameter of the beta distribution given
- values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next four argument
- values is to be calculated from the others.
- Legal range: 1..4
- iwhich = 1 : Calculate P and Q from X,Y,A and B
- iwhich = 2 : Calculate X and Y from P,Q,A and B
- iwhich = 3 : Calculate A from P,Q,X,Y and B
- iwhich = 4 : Calculate B from P,Q,X,Y and A
-
- P <--> The integral from 0 to X of the chi-square
- distribution.
- Input range: [0, 1].
-
- Q <--> 1-P.
- Input range: [0, 1].
- P + Q = 1.0.
-
- X <--> Upper limit of integration of beta density.
- Input range: [0,1].
- Search range: [0,1]
-
- Y <--> 1-X.
- Input range: [0,1].
- Search range: [0,1]
- X + Y = 1.0.
-
- A <--> The first parameter of the beta density.
- Input range: (0, +infinity).
- Search range: [1D-100,1D100]
-
- B <--> The second parameter of the beta density.
- Input range: (0, +infinity).
- Search range: [1D-100,1D100]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
- 4 if X + Y .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Cumulative distribution function (P) is calculated directly by
- code associated with the following reference.
-
- DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant
- Digit Computation of the Incomplete Beta Function Ratios. ACM
- Trans. Math. Softw. 18 (1993), 360-373.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-
- Note
-
-
- The beta density is proportional to
- t^(A-1) * (1-t)^(B-1)
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define inf 1.0e100
-#define one 1.0e0
-static int K1 = 1;
-static double K2 = 0.0e0;
-static double K3 = 1.0e0;
-static double K8 = 0.5e0;
-static double K9 = 5.0e0;
-static double fx,xhi,xlo,cum,ccum,xy,pq;
-static unsigned long qhi,qleft,qporq;
-static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q < 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 2) goto S150;
-/*
- X
-*/
- if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
- if(!(*x < 0.0e0)) goto S120;
- *bound = 0.0e0;
- goto S130;
-S120:
- *bound = 1.0e0;
-S130:
- *status = -4;
- return;
-S150:
-S140:
- if(*which == 2) goto S190;
-/*
- Y
-*/
- if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
- if(!(*y < 0.0e0)) goto S160;
- *bound = 0.0e0;
- goto S170;
-S160:
- *bound = 1.0e0;
-S170:
- *status = -5;
- return;
-S190:
-S180:
- if(*which == 3) goto S210;
-/*
- A
-*/
- if(!(*a <= 0.0e0)) goto S200;
- *bound = 0.0e0;
- *status = -6;
- return;
-S210:
-S200:
- if(*which == 4) goto S230;
-/*
- B
-*/
- if(!(*b <= 0.0e0)) goto S220;
- *bound = 0.0e0;
- *status = -7;
- return;
-S230:
-S220:
- if(*which == 1) goto S270;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
- if(!(pq < 0.0e0)) goto S240;
- *bound = 0.0e0;
- goto S250;
-S240:
- *bound = 1.0e0;
-S250:
- *status = 3;
- return;
-S270:
-S260:
- if(*which == 2) goto S310;
-/*
- X + Y
-*/
- xy = *x+*y;
- if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
- if(!(xy < 0.0e0)) goto S280;
- *bound = 0.0e0;
- goto S290;
-S280:
- *bound = 1.0e0;
-S290:
- *status = 4;
- return;
-S310:
-S300:
- if(!(*which == 1)) qporq = *p <= *q;
-/*
- Select the minimum of P or Q
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P and Q
-*/
- cumbet(x,y,a,b,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating X and Y
-*/
- T4 = atol;
- T5 = tol;
- dstzr(&K2,&K3,&T4,&T5);
- if(!qporq) goto S340;
- *status = 0;
- dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
- *y = one-*x;
-S320:
- if(!(*status == 1)) goto S330;
- cumbet(x,y,a,b,&cum,&ccum);
- fx = cum-*p;
- dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
- *y = one-*x;
- goto S320;
-S330:
- goto S370;
-S340:
- *status = 0;
- dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
- *x = one-*y;
-S350:
- if(!(*status == 1)) goto S360;
- cumbet(x,y,a,b,&cum,&ccum);
- fx = ccum-*q;
- dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
- *x = one-*y;
- goto S350;
-S370:
-S360:
- if(!(*status == -1)) goto S400;
- if(!qleft) goto S380;
- *status = 1;
- *bound = 0.0e0;
- goto S390;
-S380:
- *status = 2;
- *bound = 1.0e0;
-S400:
-S390:
- ;
- }
- else if(3 == *which) {
-/*
- Computing A
-*/
- *a = 5.0e0;
- T6 = zero;
- T7 = inf;
- T10 = atol;
- T11 = tol;
- dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
- *status = 0;
- dinvr(status,a,&fx,&qleft,&qhi);
-S410:
- if(!(*status == 1)) goto S440;
- cumbet(x,y,a,b,&cum,&ccum);
- if(!qporq) goto S420;
- fx = cum-*p;
- goto S430;
-S420:
- fx = ccum-*q;
-S430:
- dinvr(status,a,&fx,&qleft,&qhi);
- goto S410;
-S440:
- if(!(*status == -1)) goto S470;
- if(!qleft) goto S450;
- *status = 1;
- *bound = zero;
- goto S460;
-S450:
- *status = 2;
- *bound = inf;
-S470:
-S460:
- ;
- }
- else if(4 == *which) {
-/*
- Computing B
-*/
- *b = 5.0e0;
- T12 = zero;
- T13 = inf;
- T14 = atol;
- T15 = tol;
- dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
- *status = 0;
- dinvr(status,b,&fx,&qleft,&qhi);
-S480:
- if(!(*status == 1)) goto S510;
- cumbet(x,y,a,b,&cum,&ccum);
- if(!qporq) goto S490;
- fx = cum-*p;
- goto S500;
-S490:
- fx = ccum-*q;
-S500:
- dinvr(status,b,&fx,&qleft,&qhi);
- goto S480;
-S510:
- if(!(*status == -1)) goto S540;
- if(!qleft) goto S520;
- *status = 1;
- *bound = zero;
- goto S530;
-S520:
- *status = 2;
- *bound = inf;
-S530:
- ;
- }
-S540:
- return;
-#undef tol
-#undef atol
-#undef zero
-#undef inf
-#undef one
-}
-void cdfbin(int *which,double *p,double *q,double *s,double *xn,
- double *pr,double *ompr,int *status,double *bound)
-/**********************************************************************
-
- void cdfbin(int *which,double *p,double *q,double *s,double *xn,
- double *pr,double *ompr,int *status,double *bound)
-
- Cumulative Distribution Function
- BINomial distribution
-
-
- Function
-
-
- Calculates any one parameter of the binomial
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next four argument
- values is to be calculated from the others.
- Legal range: 1..4
- iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
- iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
- iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
- iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
-
- P <--> The cumulation from 0 to S of the binomial distribution.
- (Probablility of S or fewer successes in XN trials each
- with probability of success PR.)
- Input range: [0,1].
-
- Q <--> 1-P.
- Input range: [0, 1].
- P + Q = 1.0.
-
- S <--> The number of successes observed.
- Input range: [0, XN]
- Search range: [0, XN]
-
- XN <--> The number of binomial trials.
- Input range: (0, +infinity).
- Search range: [1E-100, 1E100]
-
- PR <--> The probability of success in each binomial trial.
- Input range: [0,1].
- Search range: [0,1]
-
- OMPR <--> 1-PR
- Input range: [0,1].
- Search range: [0,1]
- PR + OMPR = 1.0
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
- 4 if PR + OMPR .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.5.24 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the binomial
- distribution to the cumulative incomplete beta distribution.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-
-**********************************************************************/
-{
-#define atol 1.0e-50
-#define tol 1.0e-8
-#define zero 1.0e-100
-#define inf 1.0e100
-#define one 1.0e0
-static int K1 = 1;
-static double K2 = 0.0e0;
-static double K3 = 0.5e0;
-static double K4 = 5.0e0;
-static double K11 = 1.0e0;
-static double fx,xhi,xlo,cum,ccum,pq,prompr;
-static unsigned long qhi,qleft,qporq;
-static double T5,T6,T7,T8,T9,T10,T12,T13;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 && *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q < 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 3) goto S130;
-/*
- XN
-*/
- if(!(*xn <= 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -5;
- return;
-S130:
-S120:
- if(*which == 2) goto S170;
-/*
- S
-*/
- if(!(*s < 0.0e0 || (*which != 3 && *s > *xn))) goto S160;
- if(!(*s < 0.0e0)) goto S140;
- *bound = 0.0e0;
- goto S150;
-S140:
- *bound = *xn;
-S150:
- *status = -4;
- return;
-S170:
-S160:
- if(*which == 4) goto S210;
-/*
- PR
-*/
- if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S200;
- if(!(*pr < 0.0e0)) goto S180;
- *bound = 0.0e0;
- goto S190;
-S180:
- *bound = 1.0e0;
-S190:
- *status = -6;
- return;
-S210:
-S200:
- if(*which == 4) goto S250;
-/*
- OMPR
-*/
- if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S240;
- if(!(*ompr < 0.0e0)) goto S220;
- *bound = 0.0e0;
- goto S230;
-S220:
- *bound = 1.0e0;
-S230:
- *status = -7;
- return;
-S250:
-S240:
- if(*which == 1) goto S290;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S280;
- if(!(pq < 0.0e0)) goto S260;
- *bound = 0.0e0;
- goto S270;
-S260:
- *bound = 1.0e0;
-S270:
- *status = 3;
- return;
-S290:
-S280:
- if(*which == 4) goto S330;
-/*
- PR + OMPR
-*/
- prompr = *pr+*ompr;
- if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S320;
- if(!(prompr < 0.0e0)) goto S300;
- *bound = 0.0e0;
- goto S310;
-S300:
- *bound = 1.0e0;
-S310:
- *status = 4;
- return;
-S330:
-S320:
- if(!(*which == 1)) qporq = *p <= *q;
-/*
- Select the minimum of P or Q
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P
-*/
- cumbin(s,xn,pr,ompr,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating S
-*/
- *s = 5.0e0;
- T5 = atol;
- T6 = tol;
- dstinv(&K2,xn,&K3,&K3,&K4,&T5,&T6);
- *status = 0;
- dinvr(status,s,&fx,&qleft,&qhi);
-S340:
- if(!(*status == 1)) goto S370;
- cumbin(s,xn,pr,ompr,&cum,&ccum);
- if(!qporq) goto S350;
- fx = cum-*p;
- goto S360;
-S350:
- fx = ccum-*q;
-S360:
- dinvr(status,s,&fx,&qleft,&qhi);
- goto S340;
-S370:
- if(!(*status == -1)) goto S400;
- if(!qleft) goto S380;
- *status = 1;
- *bound = 0.0e0;
- goto S390;
-S380:
- *status = 2;
- *bound = *xn;
-S400:
-S390:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating XN
-*/
- *xn = 5.0e0;
- T7 = zero;
- T8 = inf;
- T9 = atol;
- T10 = tol;
- dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
- *status = 0;
- dinvr(status,xn,&fx,&qleft,&qhi);
-S410:
- if(!(*status == 1)) goto S440;
- cumbin(s,xn,pr,ompr,&cum,&ccum);
- if(!qporq) goto S420;
- fx = cum-*p;
- goto S430;
-S420:
- fx = ccum-*q;
-S430:
- dinvr(status,xn,&fx,&qleft,&qhi);
- goto S410;
-S440:
- if(!(*status == -1)) goto S470;
- if(!qleft) goto S450;
- *status = 1;
- *bound = zero;
- goto S460;
-S450:
- *status = 2;
- *bound = inf;
-S470:
-S460:
- ;
- }
- else if(4 == *which) {
-/*
- Calculating PR and OMPR
-*/
- T12 = atol;
- T13 = tol;
- dstzr(&K2,&K11,&T12,&T13);
- if(!qporq) goto S500;
- *status = 0;
- dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
- *ompr = one-*pr;
-S480:
- if(!(*status == 1)) goto S490;
- cumbin(s,xn,pr,ompr,&cum,&ccum);
- fx = cum-*p;
- dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
- *ompr = one-*pr;
- goto S480;
-S490:
- goto S530;
-S500:
- *status = 0;
- dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
- *pr = one-*ompr;
-S510:
- if(!(*status == 1)) goto S520;
- cumbin(s,xn,pr,ompr,&cum,&ccum);
- fx = ccum-*q;
- dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
- *pr = one-*ompr;
- goto S510;
-S530:
-S520:
- if(!(*status == -1)) goto S560;
- if(!qleft) goto S540;
- *status = 1;
- *bound = 0.0e0;
- goto S550;
-S540:
- *status = 2;
- *bound = 1.0e0;
-S550:
- ;
- }
-S560:
- return;
-#undef atol
-#undef tol
-#undef zero
-#undef inf
-#undef one
-}
-void cdfchi(int *which,double *p,double *q,double *x,double *df,
- int *status,double *bound)
-/**********************************************************************
-
- void cdfchi(int *which,double *p,double *q,double *x,double *df,
- int *status,double *bound)
-
- Cumulative Distribution Function
- CHI-Square distribution
-
-
- Function
-
-
- Calculates any one parameter of the chi-square
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next three argument
- values is to be calculated from the others.
- Legal range: 1..3
- iwhich = 1 : Calculate P and Q from X and DF
- iwhich = 2 : Calculate X from P,Q and DF
- iwhich = 3 : Calculate DF from P,Q and X
-
- P <--> The integral from 0 to X of the chi-square
- distribution.
- Input range: [0, 1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- X <--> Upper limit of integration of the non-central
- chi-square distribution.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- DF <--> Degrees of freedom of the
- chi-square distribution.
- Input range: (0, +infinity).
- Search range: [ 1E-100, 1E100]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
- 10 indicates error returned from cumgam. See
- references in cdfgam
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.4.19 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the chisqure
- distribution to the incomplete distribution.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define inf 1.0e100
-static int K1 = 1;
-static double K2 = 0.0e0;
-static double K4 = 0.5e0;
-static double K5 = 5.0e0;
-static double fx,cum,ccum,pq,porq;
-static unsigned long qhi,qleft,qporq;
-static double T3,T6,T7,T8,T9,T10,T11;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 3)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 3.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 2) goto S130;
-/*
- X
-*/
- if(!(*x < 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -4;
- return;
-S130:
-S120:
- if(*which == 3) goto S150;
-/*
- DF
-*/
- if(!(*df <= 0.0e0)) goto S140;
- *bound = 0.0e0;
- *status = -5;
- return;
-S150:
-S140:
- if(*which == 1) goto S190;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
- if(!(pq < 0.0e0)) goto S160;
- *bound = 0.0e0;
- goto S170;
-S160:
- *bound = 1.0e0;
-S170:
- *status = 3;
- return;
-S190:
-S180:
- if(*which == 1) goto S220;
-/*
- Select the minimum of P or Q
-*/
- qporq = *p <= *q;
- if(!qporq) goto S200;
- porq = *p;
- goto S210;
-S200:
- porq = *q;
-S220:
-S210:
-/*
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P and Q
-*/
- *status = 0;
- cumchi(x,df,p,q);
- if(porq > 1.5e0) {
- *status = 10;
- return;
- }
- }
- else if(2 == *which) {
-/*
- Calculating X
-*/
- *x = 5.0e0;
- T3 = inf;
- T6 = atol;
- T7 = tol;
- dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
- *status = 0;
- dinvr(status,x,&fx,&qleft,&qhi);
-S230:
- if(!(*status == 1)) goto S270;
- cumchi(x,df,&cum,&ccum);
- if(!qporq) goto S240;
- fx = cum-*p;
- goto S250;
-S240:
- fx = ccum-*q;
-S250:
- if(!(fx+porq > 1.5e0)) goto S260;
- *status = 10;
- return;
-S260:
- dinvr(status,x,&fx,&qleft,&qhi);
- goto S230;
-S270:
- if(!(*status == -1)) goto S300;
- if(!qleft) goto S280;
- *status = 1;
- *bound = 0.0e0;
- goto S290;
-S280:
- *status = 2;
- *bound = inf;
-S300:
-S290:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating DF
-*/
- *df = 5.0e0;
- T8 = zero;
- T9 = inf;
- T10 = atol;
- T11 = tol;
- dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
- *status = 0;
- dinvr(status,df,&fx,&qleft,&qhi);
-S310:
- if(!(*status == 1)) goto S350;
- cumchi(x,df,&cum,&ccum);
- if(!qporq) goto S320;
- fx = cum-*p;
- goto S330;
-S320:
- fx = ccum-*q;
-S330:
- if(!(fx+porq > 1.5e0)) goto S340;
- *status = 10;
- return;
-S340:
- dinvr(status,df,&fx,&qleft,&qhi);
- goto S310;
-S350:
- if(!(*status == -1)) goto S380;
- if(!qleft) goto S360;
- *status = 1;
- *bound = zero;
- goto S370;
-S360:
- *status = 2;
- *bound = inf;
-S370:
- ;
- }
-S380:
- return;
-#undef tol
-#undef atol
-#undef zero
-#undef inf
-}
-void cdfchn(int *which,double *p,double *q,double *x,double *df,
- double *pnonc,int *status,double *bound)
-/**********************************************************************
-
- void cdfchn(int *which,double *p,double *q,double *x,double *df,
- double *pnonc,int *status,double *bound)
-
- Cumulative Distribution Function
- Non-central Chi-Square
-
-
- Function
-
-
- Calculates any one parameter of the non-central chi-square
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next three argument
- values is to be calculated from the others.
- Input range: 1..4
- iwhich = 1 : Calculate P and Q from X and DF
- iwhich = 2 : Calculate X from P,DF and PNONC
- iwhich = 3 : Calculate DF from P,X and PNONC
- iwhich = 3 : Calculate PNONC from P,X and DF
-
- P <--> The integral from 0 to X of the non-central chi-square
- distribution.
- Input range: [0, 1-1E-16).
-
- Q <--> 1-P.
- Q is not used by this subroutine and is only included
- for similarity with other cdf* routines.
-
- X <--> Upper limit of integration of the non-central
- chi-square distribution.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- DF <--> Degrees of freedom of the non-central
- chi-square distribution.
- Input range: (0, +infinity).
- Search range: [ 1E-100, 1E100]
-
- PNONC <--> Non-centrality parameter of the non-central
- chi-square distribution.
- Input range: [0, +infinity).
- Search range: [0,1E4]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.4.25 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to compute the cumulative
- distribution function.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-
- WARNING
-
- The computation time required for this routine is proportional
- to the noncentrality parameter (PNONC). Very large values of
- this parameter can consume immense computer resources. This is
- why the search range is bounded by 10,000.
-
-**********************************************************************/
-{
-#define tent4 1.0e4
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define one ( 1.0e0 - 1.0e-16 )
-#define inf 1.0e100
-static double K1 = 0.0e0;
-static double K3 = 0.5e0;
-static double K4 = 5.0e0;
-static double fx,cum,ccum;
-static unsigned long qhi,qleft;
-static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > one)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = one;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 2) goto S90;
-/*
- X
-*/
- if(!(*x < 0.0e0)) goto S80;
- *bound = 0.0e0;
- *status = -4;
- return;
-S90:
-S80:
- if(*which == 3) goto S110;
-/*
- DF
-*/
- if(!(*df <= 0.0e0)) goto S100;
- *bound = 0.0e0;
- *status = -5;
- return;
-S110:
-S100:
- if(*which == 4) goto S130;
-/*
- PNONC
-*/
- if(!(*pnonc < 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -6;
- return;
-S130:
-S120:
-/*
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P and Q
-*/
- cumchn(x,df,pnonc,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating X
-*/
- *x = 5.0e0;
- T2 = inf;
- T5 = atol;
- T6 = tol;
- dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
- *status = 0;
- dinvr(status,x,&fx,&qleft,&qhi);
-S140:
- if(!(*status == 1)) goto S150;
- cumchn(x,df,pnonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,x,&fx,&qleft,&qhi);
- goto S140;
-S150:
- if(!(*status == -1)) goto S180;
- if(!qleft) goto S160;
- *status = 1;
- *bound = 0.0e0;
- goto S170;
-S160:
- *status = 2;
- *bound = inf;
-S180:
-S170:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating DF
-*/
- *df = 5.0e0;
- T7 = zero;
- T8 = inf;
- T9 = atol;
- T10 = tol;
- dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
- *status = 0;
- dinvr(status,df,&fx,&qleft,&qhi);
-S190:
- if(!(*status == 1)) goto S200;
- cumchn(x,df,pnonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,df,&fx,&qleft,&qhi);
- goto S190;
-S200:
- if(!(*status == -1)) goto S230;
- if(!qleft) goto S210;
- *status = 1;
- *bound = zero;
- goto S220;
-S210:
- *status = 2;
- *bound = inf;
-S230:
-S220:
- ;
- }
- else if(4 == *which) {
-/*
- Calculating PNONC
-*/
- *pnonc = 5.0e0;
- T11 = tent4;
- T12 = atol;
- T13 = tol;
- dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);
- *status = 0;
- dinvr(status,pnonc,&fx,&qleft,&qhi);
-S240:
- if(!(*status == 1)) goto S250;
- cumchn(x,df,pnonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,pnonc,&fx,&qleft,&qhi);
- goto S240;
-S250:
- if(!(*status == -1)) goto S280;
- if(!qleft) goto S260;
- *status = 1;
- *bound = zero;
- goto S270;
-S260:
- *status = 2;
- *bound = tent4;
-S270:
- ;
- }
-S280:
- return;
-#undef tent4
-#undef tol
-#undef atol
-#undef zero
-#undef one
-#undef inf
-}
-void cdff(int *which,double *p,double *q,double *f,double *dfn,
- double *dfd,int *status,double *bound)
-/**********************************************************************
-
- void cdff(int *which,double *p,double *q,double *f,double *dfn,
- double *dfd,int *status,double *bound)
-
- Cumulative Distribution Function
- F distribution
-
-
- Function
-
-
- Calculates any one parameter of the F distribution
- given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next four argument
- values is to be calculated from the others.
- Legal range: 1..4
- iwhich = 1 : Calculate P and Q from F,DFN and DFD
- iwhich = 2 : Calculate F from P,Q,DFN and DFD
- iwhich = 3 : Calculate DFN from P,Q,F and DFD
- iwhich = 4 : Calculate DFD from P,Q,F and DFN
-
- P <--> The integral from 0 to F of the f-density.
- Input range: [0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- F <--> Upper limit of integration of the f-density.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- DFN < --> Degrees of freedom of the numerator sum of squares.
- Input range: (0, +infinity).
- Search range: [ 1E-100, 1E100]
-
- DFD < --> Degrees of freedom of the denominator sum of squares.
- Input range: (0, +infinity).
- Search range: [ 1E-100, 1E100]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.6.2 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the computation
- of the cumulative distribution function for the F variate to
- that of an incomplete beta.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
- WARNING
-
- The value of the cumulative F distribution is not necessarily
- monotone in either degrees of freedom. There thus may be two
- values that provide a given CDF value. This routine assumes
- monotonicity and will find an arbitrary one of the two values.
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define inf 1.0e100
-static int K1 = 1;
-static double K2 = 0.0e0;
-static double K4 = 0.5e0;
-static double K5 = 5.0e0;
-static double pq,fx,cum,ccum;
-static unsigned long qhi,qleft,qporq;
-static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 2) goto S130;
-/*
- F
-*/
- if(!(*f < 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -4;
- return;
-S130:
-S120:
- if(*which == 3) goto S150;
-/*
- DFN
-*/
- if(!(*dfn <= 0.0e0)) goto S140;
- *bound = 0.0e0;
- *status = -5;
- return;
-S150:
-S140:
- if(*which == 4) goto S170;
-/*
- DFD
-*/
- if(!(*dfd <= 0.0e0)) goto S160;
- *bound = 0.0e0;
- *status = -6;
- return;
-S170:
-S160:
- if(*which == 1) goto S210;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
- if(!(pq < 0.0e0)) goto S180;
- *bound = 0.0e0;
- goto S190;
-S180:
- *bound = 1.0e0;
-S190:
- *status = 3;
- return;
-S210:
-S200:
- if(!(*which == 1)) qporq = *p <= *q;
-/*
- Select the minimum of P or Q
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P
-*/
- cumf(f,dfn,dfd,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating F
-*/
- *f = 5.0e0;
- T3 = inf;
- T6 = atol;
- T7 = tol;
- dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
- *status = 0;
- dinvr(status,f,&fx,&qleft,&qhi);
-S220:
- if(!(*status == 1)) goto S250;
- cumf(f,dfn,dfd,&cum,&ccum);
- if(!qporq) goto S230;
- fx = cum-*p;
- goto S240;
-S230:
- fx = ccum-*q;
-S240:
- dinvr(status,f,&fx,&qleft,&qhi);
- goto S220;
-S250:
- if(!(*status == -1)) goto S280;
- if(!qleft) goto S260;
- *status = 1;
- *bound = 0.0e0;
- goto S270;
-S260:
- *status = 2;
- *bound = inf;
-S280:
-S270:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating DFN
-*/
- *dfn = 5.0e0;
- T8 = zero;
- T9 = inf;
- T10 = atol;
- T11 = tol;
- dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
- *status = 0;
- dinvr(status,dfn,&fx,&qleft,&qhi);
-S290:
- if(!(*status == 1)) goto S320;
- cumf(f,dfn,dfd,&cum,&ccum);
- if(!qporq) goto S300;
- fx = cum-*p;
- goto S310;
-S300:
- fx = ccum-*q;
-S310:
- dinvr(status,dfn,&fx,&qleft,&qhi);
- goto S290;
-S320:
- if(!(*status == -1)) goto S350;
- if(!qleft) goto S330;
- *status = 1;
- *bound = zero;
- goto S340;
-S330:
- *status = 2;
- *bound = inf;
-S350:
-S340:
- ;
- }
- else if(4 == *which) {
-/*
- Calculating DFD
-*/
- *dfd = 5.0e0;
- T12 = zero;
- T13 = inf;
- T14 = atol;
- T15 = tol;
- dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15);
- *status = 0;
- dinvr(status,dfd,&fx,&qleft,&qhi);
-S360:
- if(!(*status == 1)) goto S390;
- cumf(f,dfn,dfd,&cum,&ccum);
- if(!qporq) goto S370;
- fx = cum-*p;
- goto S380;
-S370:
- fx = ccum-*q;
-S380:
- dinvr(status,dfd,&fx,&qleft,&qhi);
- goto S360;
-S390:
- if(!(*status == -1)) goto S420;
- if(!qleft) goto S400;
- *status = 1;
- *bound = zero;
- goto S410;
-S400:
- *status = 2;
- *bound = inf;
-S410:
- ;
- }
-S420:
- return;
-#undef tol
-#undef atol
-#undef zero
-#undef inf
-}
-void cdffnc(int *which,double *p,double *q,double *f,double *dfn,
- double *dfd,double *phonc,int *status,double *bound)
-/**********************************************************************
-
- void cdffnc(int *which,double *p,double *q,double *f,double *dfn,
- double *dfd,double *phonc,int *status,double *bound)
-
- Cumulative Distribution Function
- Non-central F distribution
-
-
- Function
-
-
- Calculates any one parameter of the Non-central F
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next five argument
- values is to be calculated from the others.
- Legal range: 1..5
- iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
- iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
- iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
- iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
- iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
-
- P <--> The integral from 0 to F of the non-central f-density.
- Input range: [0,1-1E-16).
-
- Q <--> 1-P.
- Q is not used by this subroutine and is only included
- for similarity with other cdf* routines.
-
- F <--> Upper limit of integration of the non-central f-density.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- DFN < --> Degrees of freedom of the numerator sum of squares.
- Input range: (0, +infinity).
- Search range: [ 1E-100, 1E100]
-
- DFD < --> Degrees of freedom of the denominator sum of squares.
- Must be in range: (0, +infinity).
- Input range: (0, +infinity).
- Search range: [ 1E-100, 1E100]
-
- PNONC <-> The non-centrality parameter
- Input range: [0,infinity)
- Search range: [0,1E4]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.6.20 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to compute the cumulative
- distribution function.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
- WARNING
-
- The computation time required for this routine is proportional
- to the noncentrality parameter (PNONC). Very large values of
- this parameter can consume immense computer resources. This is
- why the search range is bounded by 10,000.
-
- WARNING
-
- The value of the cumulative noncentral F distribution is not
- necessarily monotone in either degrees of freedom. There thus
- may be two values that provide a given CDF value. This routine
- assumes monotonicity and will find an arbitrary one of the two
- values.
-
-**********************************************************************/
-{
-#define tent4 1.0e4
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define one ( 1.0e0 - 1.0e-16 )
-#define inf 1.0e100
-static double K1 = 0.0e0;
-static double K3 = 0.5e0;
-static double K4 = 5.0e0;
-static double fx,cum,ccum;
-static unsigned long qhi,qleft;
-static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 5)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 5.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > one)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = one;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 2) goto S90;
-/*
- F
-*/
- if(!(*f < 0.0e0)) goto S80;
- *bound = 0.0e0;
- *status = -4;
- return;
-S90:
-S80:
- if(*which == 3) goto S110;
-/*
- DFN
-*/
- if(!(*dfn <= 0.0e0)) goto S100;
- *bound = 0.0e0;
- *status = -5;
- return;
-S110:
-S100:
- if(*which == 4) goto S130;
-/*
- DFD
-*/
- if(!(*dfd <= 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -6;
- return;
-S130:
-S120:
- if(*which == 5) goto S150;
-/*
- PHONC
-*/
- if(!(*phonc < 0.0e0)) goto S140;
- *bound = 0.0e0;
- *status = -7;
- return;
-S150:
-S140:
-/*
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P
-*/
- cumfnc(f,dfn,dfd,phonc,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating F
-*/
- *f = 5.0e0;
- T2 = inf;
- T5 = atol;
- T6 = tol;
- dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
- *status = 0;
- dinvr(status,f,&fx,&qleft,&qhi);
-S160:
- if(!(*status == 1)) goto S170;
- cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,f,&fx,&qleft,&qhi);
- goto S160;
-S170:
- if(!(*status == -1)) goto S200;
- if(!qleft) goto S180;
- *status = 1;
- *bound = 0.0e0;
- goto S190;
-S180:
- *status = 2;
- *bound = inf;
-S200:
-S190:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating DFN
-*/
- *dfn = 5.0e0;
- T7 = zero;
- T8 = inf;
- T9 = atol;
- T10 = tol;
- dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
- *status = 0;
- dinvr(status,dfn,&fx,&qleft,&qhi);
-S210:
- if(!(*status == 1)) goto S220;
- cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,dfn,&fx,&qleft,&qhi);
- goto S210;
-S220:
- if(!(*status == -1)) goto S250;
- if(!qleft) goto S230;
- *status = 1;
- *bound = zero;
- goto S240;
-S230:
- *status = 2;
- *bound = inf;
-S250:
-S240:
- ;
- }
- else if(4 == *which) {
-/*
- Calculating DFD
-*/
- *dfd = 5.0e0;
- T11 = zero;
- T12 = inf;
- T13 = atol;
- T14 = tol;
- dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
- *status = 0;
- dinvr(status,dfd,&fx,&qleft,&qhi);
-S260:
- if(!(*status == 1)) goto S270;
- cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,dfd,&fx,&qleft,&qhi);
- goto S260;
-S270:
- if(!(*status == -1)) goto S300;
- if(!qleft) goto S280;
- *status = 1;
- *bound = zero;
- goto S290;
-S280:
- *status = 2;
- *bound = inf;
-S300:
-S290:
- ;
- }
- else if(5 == *which) {
-/*
- Calculating PHONC
-*/
- *phonc = 5.0e0;
- T15 = tent4;
- T16 = atol;
- T17 = tol;
- dstinv(&K1,&T15,&K3,&K3,&K4,&T16,&T17);
- *status = 0;
- dinvr(status,phonc,&fx,&qleft,&qhi);
-S310:
- if(!(*status == 1)) goto S320;
- cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
- fx = cum-*p;
- dinvr(status,phonc,&fx,&qleft,&qhi);
- goto S310;
-S320:
- if(!(*status == -1)) goto S350;
- if(!qleft) goto S330;
- *status = 1;
- *bound = 0.0e0;
- goto S340;
-S330:
- *status = 2;
- *bound = tent4;
-S340:
- ;
- }
-S350:
- return;
-#undef tent4
-#undef tol
-#undef atol
-#undef zero
-#undef one
-#undef inf
-}
-void cdfgam(int *which,double *p,double *q,double *x,double *shape,
- double *scale,int *status,double *bound)
-/**********************************************************************
-
- void cdfgam(int *which,double *p,double *q,double *x,double *shape,
- double *scale,int *status,double *bound)
-
- Cumulative Distribution Function
- GAMma Distribution
-
-
- Function
-
-
- Calculates any one parameter of the gamma
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next four argument
- values is to be calculated from the others.
- Legal range: 1..4
- iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE
- iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE
- iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE
- iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE
-
- P <--> The integral from 0 to X of the gamma density.
- Input range: [0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- X <--> The upper limit of integration of the gamma density.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- SHAPE <--> The shape parameter of the gamma density.
- Input range: (0, +infinity).
- Search range: [1E-100,1E100]
-
- SCALE <--> The scale parameter of the gamma density.
- Input range: (0, +infinity).
- Search range: (1E-100,1E100]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
- 10 if the gamma or inverse gamma routine cannot
- compute the answer. Usually happens only for
- X and SHAPE very large (gt 1E10 or more)
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Cumulative distribution function (P) is calculated directly by
- the code associated with:
-
- DiDinato, A. R. and Morris, A. H. Computation of the incomplete
- gamma function ratios and their inverse. ACM Trans. Math.
- Softw. 12 (1986), 377-393.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-
- Note
-
-
-
- The gamma density is proportional to
- T**(SHAPE - 1) * EXP(- SCALE * T)
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define inf 1.0e100
-static int K1 = 1;
-static double K5 = 0.5e0;
-static double K6 = 5.0e0;
-static double xx,fx,xscale,cum,ccum,pq,porq;
-static int ierr;
-static unsigned long qhi,qleft,qporq;
-static double T2,T3,T4,T7,T8,T9;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 2) goto S130;
-/*
- X
-*/
- if(!(*x < 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -4;
- return;
-S130:
-S120:
- if(*which == 3) goto S150;
-/*
- SHAPE
-*/
- if(!(*shape <= 0.0e0)) goto S140;
- *bound = 0.0e0;
- *status = -5;
- return;
-S150:
-S140:
- if(*which == 4) goto S170;
-/*
- SCALE
-*/
- if(!(*scale <= 0.0e0)) goto S160;
- *bound = 0.0e0;
- *status = -6;
- return;
-S170:
-S160:
- if(*which == 1) goto S210;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
- if(!(pq < 0.0e0)) goto S180;
- *bound = 0.0e0;
- goto S190;
-S180:
- *bound = 1.0e0;
-S190:
- *status = 3;
- return;
-S210:
-S200:
- if(*which == 1) goto S240;
-/*
- Select the minimum of P or Q
-*/
- qporq = *p <= *q;
- if(!qporq) goto S220;
- porq = *p;
- goto S230;
-S220:
- porq = *q;
-S240:
-S230:
-/*
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P
-*/
- *status = 0;
- xscale = *x**scale;
- cumgam(&xscale,shape,p,q);
- if(porq > 1.5e0) *status = 10;
- }
- else if(2 == *which) {
-/*
- Computing X
-*/
- T2 = -1.0e0;
- gaminv(shape,&xx,&T2,p,q,&ierr);
- if(ierr < 0.0e0) {
- *status = 10;
- return;
- }
- else {
- *x = xx/ *scale;
- *status = 0;
- }
- }
- else if(3 == *which) {
-/*
- Computing SHAPE
-*/
- *shape = 5.0e0;
- xscale = *x**scale;
- T3 = zero;
- T4 = inf;
- T7 = atol;
- T8 = tol;
- dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
- *status = 0;
- dinvr(status,shape,&fx,&qleft,&qhi);
-S250:
- if(!(*status == 1)) goto S290;
- cumgam(&xscale,shape,&cum,&ccum);
- if(!qporq) goto S260;
- fx = cum-*p;
- goto S270;
-S260:
- fx = ccum-*q;
-S270:
- if(!((qporq && cum > 1.5e0) || (!qporq && ccum > 1.5e0))) goto S280;
- *status = 10;
- return;
-S280:
- dinvr(status,shape,&fx,&qleft,&qhi);
- goto S250;
-S290:
- if(!(*status == -1)) goto S320;
- if(!qleft) goto S300;
- *status = 1;
- *bound = zero;
- goto S310;
-S300:
- *status = 2;
- *bound = inf;
-S320:
-S310:
- ;
- }
- else if(4 == *which) {
-/*
- Computing SCALE
-*/
- T9 = -1.0e0;
- gaminv(shape,&xx,&T9,p,q,&ierr);
- if(ierr < 0.0e0) {
- *status = 10;
- return;
- }
- else {
- *scale = xx/ *x;
- *status = 0;
- }
- }
- return;
-#undef tol
-#undef atol
-#undef zero
-#undef inf
-}
-void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
- double *pr,double *ompr,int *status,double *bound)
-/**********************************************************************
-
- void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
- double *pr,double *ompr,int *status,double *bound)
-
- Cumulative Distribution Function
- Negative BiNomial distribution
-
-
- Function
-
-
- Calculates any one parameter of the negative binomial
- distribution given values for the others.
-
- The cumulative negative binomial distribution returns the
- probability that there will be F or fewer failures before the
- XNth success in binomial trials each of which has probability of
- success PR.
-
- The individual term of the negative binomial is the probability of
- S failures before XN successes and is
- Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next four argument
- values is to be calculated from the others.
- Legal range: 1..4
- iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
- iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
- iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
- iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
-
- P <--> The cumulation from 0 to S of the negative
- binomial distribution.
- Input range: [0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- S <--> The upper limit of cumulation of the binomial distribution.
- There are F or fewer failures before the XNth success.
- Input range: [0, +infinity).
- Search range: [0, 1E100]
-
- XN <--> The number of successes.
- Input range: [0, +infinity).
- Search range: [0, 1E100]
-
- PR <--> The probability of success in each binomial trial.
- Input range: [0,1].
- Search range: [0,1].
-
- OMPR <--> 1-PR
- Input range: [0,1].
- Search range: [0,1]
- PR + OMPR = 1.0
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
- 4 if PR + OMPR .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.5.26 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce calculation of
- the cumulative distribution function to that of an incomplete
- beta.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define inf 1.0e100
-#define one 1.0e0
-static int K1 = 1;
-static double K2 = 0.0e0;
-static double K4 = 0.5e0;
-static double K5 = 5.0e0;
-static double K11 = 1.0e0;
-static double fx,xhi,xlo,pq,prompr,cum,ccum;
-static unsigned long qhi,qleft,qporq;
-static double T3,T6,T7,T8,T9,T10,T12,T13;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 2) goto S130;
-/*
- S
-*/
- if(!(*s < 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -4;
- return;
-S130:
-S120:
- if(*which == 3) goto S150;
-/*
- XN
-*/
- if(!(*xn < 0.0e0)) goto S140;
- *bound = 0.0e0;
- *status = -5;
- return;
-S150:
-S140:
- if(*which == 4) goto S190;
-/*
- PR
-*/
- if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
- if(!(*pr < 0.0e0)) goto S160;
- *bound = 0.0e0;
- goto S170;
-S160:
- *bound = 1.0e0;
-S170:
- *status = -6;
- return;
-S190:
-S180:
- if(*which == 4) goto S230;
-/*
- OMPR
-*/
- if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
- if(!(*ompr < 0.0e0)) goto S200;
- *bound = 0.0e0;
- goto S210;
-S200:
- *bound = 1.0e0;
-S210:
- *status = -7;
- return;
-S230:
-S220:
- if(*which == 1) goto S270;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
- if(!(pq < 0.0e0)) goto S240;
- *bound = 0.0e0;
- goto S250;
-S240:
- *bound = 1.0e0;
-S250:
- *status = 3;
- return;
-S270:
-S260:
- if(*which == 4) goto S310;
-/*
- PR + OMPR
-*/
- prompr = *pr+*ompr;
- if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
- if(!(prompr < 0.0e0)) goto S280;
- *bound = 0.0e0;
- goto S290;
-S280:
- *bound = 1.0e0;
-S290:
- *status = 4;
- return;
-S310:
-S300:
- if(!(*which == 1)) qporq = *p <= *q;
-/*
- Select the minimum of P or Q
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P
-*/
- cumnbn(s,xn,pr,ompr,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating S
-*/
- *s = 5.0e0;
- T3 = inf;
- T6 = atol;
- T7 = tol;
- dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
- *status = 0;
- dinvr(status,s,&fx,&qleft,&qhi);
-S320:
- if(!(*status == 1)) goto S350;
- cumnbn(s,xn,pr,ompr,&cum,&ccum);
- if(!qporq) goto S330;
- fx = cum-*p;
- goto S340;
-S330:
- fx = ccum-*q;
-S340:
- dinvr(status,s,&fx,&qleft,&qhi);
- goto S320;
-S350:
- if(!(*status == -1)) goto S380;
- if(!qleft) goto S360;
- *status = 1;
- *bound = 0.0e0;
- goto S370;
-S360:
- *status = 2;
- *bound = inf;
-S380:
-S370:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating XN
-*/
- *xn = 5.0e0;
- T8 = inf;
- T9 = atol;
- T10 = tol;
- dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
- *status = 0;
- dinvr(status,xn,&fx,&qleft,&qhi);
-S390:
- if(!(*status == 1)) goto S420;
- cumnbn(s,xn,pr,ompr,&cum,&ccum);
- if(!qporq) goto S400;
- fx = cum-*p;
- goto S410;
-S400:
- fx = ccum-*q;
-S410:
- dinvr(status,xn,&fx,&qleft,&qhi);
- goto S390;
-S420:
- if(!(*status == -1)) goto S450;
- if(!qleft) goto S430;
- *status = 1;
- *bound = 0.0e0;
- goto S440;
-S430:
- *status = 2;
- *bound = inf;
-S450:
-S440:
- ;
- }
- else if(4 == *which) {
-/*
- Calculating PR and OMPR
-*/
- T12 = atol;
- T13 = tol;
- dstzr(&K2,&K11,&T12,&T13);
- if(!qporq) goto S480;
- *status = 0;
- dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
- *ompr = one-*pr;
-S460:
- if(!(*status == 1)) goto S470;
- cumnbn(s,xn,pr,ompr,&cum,&ccum);
- fx = cum-*p;
- dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
- *ompr = one-*pr;
- goto S460;
-S470:
- goto S510;
-S480:
- *status = 0;
- dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
- *pr = one-*ompr;
-S490:
- if(!(*status == 1)) goto S500;
- cumnbn(s,xn,pr,ompr,&cum,&ccum);
- fx = ccum-*q;
- dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
- *pr = one-*ompr;
- goto S490;
-S510:
-S500:
- if(!(*status == -1)) goto S540;
- if(!qleft) goto S520;
- *status = 1;
- *bound = 0.0e0;
- goto S530;
-S520:
- *status = 2;
- *bound = 1.0e0;
-S530:
- ;
- }
-S540:
- return;
-#undef tol
-#undef atol
-#undef inf
-#undef one
-}
-void cdfnor(int *which,double *p,double *q,double *x,double *mean,
- double *sd,int *status,double *bound)
-/**********************************************************************
-
- void cdfnor(int *which,double *p,double *q,double *x,double *mean,
- double *sd,int *status,double *bound)
-
- Cumulative Distribution Function
- NORmal distribution
-
-
- Function
-
-
- Calculates any one parameter of the normal
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which of the next parameter
- values is to be calculated using values of the others.
- Legal range: 1..4
- iwhich = 1 : Calculate P and Q from X,MEAN and SD
- iwhich = 2 : Calculate X from P,Q,MEAN and SD
- iwhich = 3 : Calculate MEAN from P,Q,X and SD
- iwhich = 4 : Calculate SD from P,Q,X and MEAN
-
- P <--> The integral from -infinity to X of the normal density.
- Input range: (0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- X < --> Upper limit of integration of the normal-density.
- Input range: ( -infinity, +infinity)
-
- MEAN <--> The mean of the normal density.
- Input range: (-infinity, +infinity)
-
- SD <--> Standard Deviation of the normal density.
- Input range: (0, +infinity).
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
-
-
- A slightly modified version of ANORM from
-
- Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
- Package of Special Function Routines and Test Drivers"
- acm Transactions on Mathematical Software. 19, 22-32.
-
- is used to calulate the cumulative standard normal distribution.
-
- The rational functions from pages 90-95 of Kennedy and Gentle,
- Statistical Computing, Marcel Dekker, NY, 1980 are used as
- starting values to Newton's Iterations which compute the inverse
- standard normal. Therefore no searches are necessary for any
- parameter.
-
- For X < -15, the asymptotic expansion for the normal is used as
- the starting value in finding the inverse standard normal.
- This is formula 26.2.12 of Abramowitz and Stegun.
-
-
- Note
-
-
- The normal density is proportional to
- exp( - 0.5 * (( X - MEAN)/SD)**2)
-
-**********************************************************************/
-{
-static int K1 = 1;
-static double z,pq;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- *status = 0;
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 4.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p <= 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 1) goto S150;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S140;
- if(!(pq < 0.0e0)) goto S120;
- *bound = 0.0e0;
- goto S130;
-S120:
- *bound = 1.0e0;
-S130:
- *status = 3;
- return;
-S150:
-S140:
- if(*which == 4) goto S170;
-/*
- SD
-*/
- if(!(*sd <= 0.0e0)) goto S160;
- *bound = 0.0e0;
- *status = -6;
- return;
-S170:
-S160:
-/*
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Computing P
-*/
- z = (*x-*mean)/ *sd;
- cumnor(&z,p,q);
- }
- else if(2 == *which) {
-/*
- Computing X
-*/
- z = dinvnr(p,q);
- *x = *sd*z+*mean;
- }
- else if(3 == *which) {
-/*
- Computing the MEAN
-*/
- z = dinvnr(p,q);
- *mean = *x-*sd*z;
- }
- else if(4 == *which) {
-/*
- Computing SD
-*/
- z = dinvnr(p,q);
- *sd = (*x-*mean)/z;
- }
- return;
-}
-void cdfpoi(int *which,double *p,double *q,double *s,double *xlam,
- int *status,double *bound)
-/**********************************************************************
-
- void cdfpoi(int *which,double *p,double *q,double *s,double *xlam,
- int *status,double *bound)
-
- Cumulative Distribution Function
- POIsson distribution
-
-
- Function
-
-
- Calculates any one parameter of the Poisson
- distribution given values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which argument
- value is to be calculated from the others.
- Legal range: 1..3
- iwhich = 1 : Calculate P and Q from S and XLAM
- iwhich = 2 : Calculate A from P,Q and XLAM
- iwhich = 3 : Calculate XLAM from P,Q and S
-
- P <--> The cumulation from 0 to S of the poisson density.
- Input range: [0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- S <--> Upper limit of cumulation of the Poisson.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- XLAM <--> Mean of the Poisson distribution.
- Input range: [0, +infinity).
- Search range: [0,1E100]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.4.21 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the computation
- of the cumulative distribution function to that of computing a
- chi-square, hence an incomplete gamma function.
-
- Cumulative distribution function (P) is calculated directly.
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define inf 1.0e100
-static int K1 = 1;
-static double K2 = 0.0e0;
-static double K4 = 0.5e0;
-static double K5 = 5.0e0;
-static double fx,cum,ccum,pq;
-static unsigned long qhi,qleft,qporq;
-static double T3,T6,T7,T8,T9,T10;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 3)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 3.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 2) goto S130;
-/*
- S
-*/
- if(!(*s < 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -4;
- return;
-S130:
-S120:
- if(*which == 3) goto S150;
-/*
- XLAM
-*/
- if(!(*xlam < 0.0e0)) goto S140;
- *bound = 0.0e0;
- *status = -5;
- return;
-S150:
-S140:
- if(*which == 1) goto S190;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
- if(!(pq < 0.0e0)) goto S160;
- *bound = 0.0e0;
- goto S170;
-S160:
- *bound = 1.0e0;
-S170:
- *status = 3;
- return;
-S190:
-S180:
- if(!(*which == 1)) qporq = *p <= *q;
-/*
- Select the minimum of P or Q
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Calculating P
-*/
- cumpoi(s,xlam,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Calculating S
-*/
- *s = 5.0e0;
- T3 = inf;
- T6 = atol;
- T7 = tol;
- dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
- *status = 0;
- dinvr(status,s,&fx,&qleft,&qhi);
-S200:
- if(!(*status == 1)) goto S230;
- cumpoi(s,xlam,&cum,&ccum);
- if(!qporq) goto S210;
- fx = cum-*p;
- goto S220;
-S210:
- fx = ccum-*q;
-S220:
- dinvr(status,s,&fx,&qleft,&qhi);
- goto S200;
-S230:
- if(!(*status == -1)) goto S260;
- if(!qleft) goto S240;
- *status = 1;
- *bound = 0.0e0;
- goto S250;
-S240:
- *status = 2;
- *bound = inf;
-S260:
-S250:
- ;
- }
- else if(3 == *which) {
-/*
- Calculating XLAM
-*/
- *xlam = 5.0e0;
- T8 = inf;
- T9 = atol;
- T10 = tol;
- dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
- *status = 0;
- dinvr(status,xlam,&fx,&qleft,&qhi);
-S270:
- if(!(*status == 1)) goto S300;
- cumpoi(s,xlam,&cum,&ccum);
- if(!qporq) goto S280;
- fx = cum-*p;
- goto S290;
-S280:
- fx = ccum-*q;
-S290:
- dinvr(status,xlam,&fx,&qleft,&qhi);
- goto S270;
-S300:
- if(!(*status == -1)) goto S330;
- if(!qleft) goto S310;
- *status = 1;
- *bound = 0.0e0;
- goto S320;
-S310:
- *status = 2;
- *bound = inf;
-S320:
- ;
- }
-S330:
- return;
-#undef tol
-#undef atol
-#undef inf
-}
-void cdft(int *which,double *p,double *q,double *t,double *df,
- int *status,double *bound)
-/**********************************************************************
-
- void cdft(int *which,double *p,double *q,double *t,double *df,
- int *status,double *bound)
-
- Cumulative Distribution Function
- T distribution
-
-
- Function
-
-
- Calculates any one parameter of the t distribution given
- values for the others.
-
-
- Arguments
-
-
- WHICH --> Integer indicating which argument
- values is to be calculated from the others.
- Legal range: 1..3
- iwhich = 1 : Calculate P and Q from T and DF
- iwhich = 2 : Calculate T from P,Q and DF
- iwhich = 3 : Calculate DF from P,Q and T
-
- P <--> The integral from -infinity to t of the t-density.
- Input range: (0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- T <--> Upper limit of integration of the t-density.
- Input range: ( -infinity, +infinity).
- Search range: [ -1E100, 1E100 ]
-
- DF <--> Degrees of freedom of the t-distribution.
- Input range: (0 , +infinity).
- Search range: [1e-100, 1E10]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
-
- Method
-
-
- Formula 26.5.27 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the computation
- of the cumulative distribution function to that of an incomplete
- beta.
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-**********************************************************************/
-{
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define inf 1.0e100
-#define rtinf 1.0e100
-#define maxdf 1.0e10
-static int K1 = 1;
-static double K4 = 0.5e0;
-static double K5 = 5.0e0;
-static double fx,cum,ccum,pq;
-static unsigned long qhi,qleft,qporq;
-static double T2,T3,T6,T7,T8,T9,T10,T11;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Check arguments
-*/
- if(!(*which < 1 || *which > 3)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 3.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
-/*
- P
-*/
- if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
- if(!(*p <= 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = 1.0e0;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 1) goto S110;
-/*
- Q
-*/
- if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
- if(!(*q <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- goto S90;
-S80:
- *bound = 1.0e0;
-S90:
- *status = -3;
- return;
-S110:
-S100:
- if(*which == 3) goto S130;
-/*
- DF
-*/
- if(!(*df <= 0.0e0)) goto S120;
- *bound = 0.0e0;
- *status = -5;
- return;
-S130:
-S120:
- if(*which == 1) goto S170;
-/*
- P + Q
-*/
- pq = *p+*q;
- if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S160;
- if(!(pq < 0.0e0)) goto S140;
- *bound = 0.0e0;
- goto S150;
-S140:
- *bound = 1.0e0;
-S150:
- *status = 3;
- return;
-S170:
-S160:
- if(!(*which == 1)) qporq = *p <= *q;
-/*
- Select the minimum of P or Q
- Calculate ANSWERS
-*/
- if(1 == *which) {
-/*
- Computing P and Q
-*/
- cumt(t,df,p,q);
- *status = 0;
- }
- else if(2 == *which) {
-/*
- Computing T
- .. Get initial approximation for T
-*/
- *t = dt1(p,q,df);
- T2 = -rtinf;
- T3 = rtinf;
- T6 = atol;
- T7 = tol;
- dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
- *status = 0;
- dinvr(status,t,&fx,&qleft,&qhi);
-S180:
- if(!(*status == 1)) goto S210;
- cumt(t,df,&cum,&ccum);
- if(!qporq) goto S190;
- fx = cum-*p;
- goto S200;
-S190:
- fx = ccum-*q;
-S200:
- dinvr(status,t,&fx,&qleft,&qhi);
- goto S180;
-S210:
- if(!(*status == -1)) goto S240;
- if(!qleft) goto S220;
- *status = 1;
- *bound = -rtinf;
- goto S230;
-S220:
- *status = 2;
- *bound = rtinf;
-S240:
-S230:
- ;
- }
- else if(3 == *which) {
-/*
- Computing DF
-*/
- *df = 5.0e0;
- T8 = zero;
- T9 = maxdf;
- T10 = atol;
- T11 = tol;
- dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
- *status = 0;
- dinvr(status,df,&fx,&qleft,&qhi);
-S250:
- if(!(*status == 1)) goto S280;
- cumt(t,df,&cum,&ccum);
- if(!qporq) goto S260;
- fx = cum-*p;
- goto S270;
-S260:
- fx = ccum-*q;
-S270:
- dinvr(status,df,&fx,&qleft,&qhi);
- goto S250;
-S280:
- if(!(*status == -1)) goto S310;
- if(!qleft) goto S290;
- *status = 1;
- *bound = zero;
- goto S300;
-S290:
- *status = 2;
- *bound = maxdf;
-S300:
- ;
- }
-S310:
- return;
-#undef tol
-#undef atol
-#undef zero
-#undef inf
-#undef rtinf
-#undef maxdf
-}
-void cdftnc(int *which,double *p,double *q,double *t,double *df,
- double *pnonc,int *status,double *bound)
-/**********************************************************************
-
- void cdftnc(int *which,double *p,double *q,double *t,double *df,
- double *pnonc,int *status,double *bound)
-
- Cumulative Distribution Function
- Non-Central T distribution
-
- Function
-
- Calculates any one parameter of the noncentral t distribution give
- values for the others.
-
- Arguments
-
- WHICH --> Integer indicating which argument
- values is to be calculated from the others.
- Legal range: 1..3
- iwhich = 1 : Calculate P and Q from T,DF,PNONC
- iwhich = 2 : Calculate T from P,Q,DF,PNONC
- iwhich = 3 : Calculate DF from P,Q,T
- iwhich = 4 : Calculate PNONC from P,Q,DF,T
-
- P <--> The integral from -infinity to t of the noncentral t-den
- Input range: (0,1].
-
- Q <--> 1-P.
- Input range: (0, 1].
- P + Q = 1.0.
-
- T <--> Upper limit of integration of the noncentral t-density.
- Input range: ( -infinity, +infinity).
- Search range: [ -1E100, 1E100 ]
-
- DF <--> Degrees of freedom of the noncentral t-distribution.
- Input range: (0 , +infinity).
- Search range: [1e-100, 1E10]
-
- PNONC <--> Noncentrality parameter of the noncentral t-distributio
- Input range: [-infinity , +infinity).
- Search range: [-1e4, 1E4]
-
- STATUS <-- 0 if calculation completed correctly
- -I if input parameter number I is out of range
- 1 if answer appears to be lower than lowest
- search bound
- 2 if answer appears to be higher than greatest
- search bound
- 3 if P + Q .ne. 1
-
- BOUND <-- Undefined if STATUS is 0
-
- Bound exceeded by parameter number I if STATUS
- is negative.
-
- Lower search bound if STATUS is 1.
-
- Upper search bound if STATUS is 2.
-
- Method
-
- Upper tail of the cumulative noncentral t is calculated usin
- formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuou
- Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
-
- Computation of other parameters involve a seach for a value that
- produces the desired value of P. The search relies on the
- monotinicity of P with the other parameter.
-
-**********************************************************************/
-{
-#define tent4 1.0e4
-#define tol 1.0e-8
-#define atol 1.0e-50
-#define zero 1.0e-100
-#define one ( 1.0e0 - 1.0e-16 )
-#define inf 1.0e100
-static double K3 = 0.5e0;
-static double K4 = 5.0e0;
-static double ccum,cum,fx;
-static unsigned long qhi,qleft;
-static double T1,T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*which < 1 || *which > 4)) goto S30;
- if(!(*which < 1)) goto S10;
- *bound = 1.0e0;
- goto S20;
-S10:
- *bound = 5.0e0;
-S20:
- *status = -1;
- return;
-S30:
- if(*which == 1) goto S70;
- if(!(*p < 0.0e0 || *p > one)) goto S60;
- if(!(*p < 0.0e0)) goto S40;
- *bound = 0.0e0;
- goto S50;
-S40:
- *bound = one;
-S50:
- *status = -2;
- return;
-S70:
-S60:
- if(*which == 3) goto S90;
- if(!(*df <= 0.0e0)) goto S80;
- *bound = 0.0e0;
- *status = -5;
- return;
-S90:
-S80:
- if(*which == 4) goto S100;
-S100:
- if(1 == *which) {
- cumtnc(t,df,pnonc,p,q);
- *status = 0;
- }
- else if(2 == *which) {
- *t = 5.0e0;
- T1 = -inf;
- T2 = inf;
- T5 = atol;
- T6 = tol;
- dstinv(&T1,&T2,&K3,&K3,&K4,&T5,&T6);
- *status = 0;
- dinvr(status,t,&fx,&qleft,&qhi);
-S110:
- if(!(*status == 1)) goto S120;
- cumtnc(t,df,pnonc,&cum,&ccum);
- fx = cum - *p;
- dinvr(status,t,&fx,&qleft,&qhi);
- goto S110;
-S120:
- if(!(*status == -1)) goto S150;
- if(!qleft) goto S130;
- *status = 1;
- *bound = -inf;
- goto S140;
-S130:
- *status = 2;
- *bound = inf;
-S150:
-S140:
- ;
- }
- else if(3 == *which) {
- *df = 5.0e0;
- T7 = zero;
- T8 = tent4;
- T9 = atol;
- T10 = tol;
- dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
- *status = 0;
- dinvr(status,df,&fx,&qleft,&qhi);
-S160:
- if(!(*status == 1)) goto S170;
- cumtnc(t,df,pnonc,&cum,&ccum);
- fx = cum - *p;
- dinvr(status,df,&fx,&qleft,&qhi);
- goto S160;
-S170:
- if(!(*status == -1)) goto S200;
- if(!qleft) goto S180;
- *status = 1;
- *bound = zero;
- goto S190;
-S180:
- *status = 2;
- *bound = inf;
-S200:
-S190:
- ;
- }
- else if(4 == *which) {
- *pnonc = 5.0e0;
- T11 = -tent4;
- T12 = tent4;
- T13 = atol;
- T14 = tol;
- dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
- *status = 0;
- dinvr(status,pnonc,&fx,&qleft,&qhi);
-S210:
- if(!(*status == 1)) goto S220;
- cumtnc(t,df,pnonc,&cum,&ccum);
- fx = cum - *p;
- dinvr(status,pnonc,&fx,&qleft,&qhi);
- goto S210;
-S220:
- if(!(*status == -1)) goto S250;
- if(!qleft) goto S230;
- *status = 1;
- *bound = 0.0e0;
- goto S240;
-S230:
- *status = 2;
- *bound = tent4;
-S240:
- ;
- }
-S250:
- return;
-#undef tent4
-#undef tol
-#undef atol
-#undef zero
-#undef one
-#undef inf
-}
-void cumbet(double *x,double *y,double *a,double *b,double *cum,
- double *ccum)
-/*
-**********************************************************************
-
- void cumbet(double *x,double *y,double *a,double *b,double *cum,
- double *ccum)
-
- Double precision cUMulative incomplete BETa distribution
-
-
- Function
-
-
- Calculates the cdf to X of the incomplete beta distribution
- with parameters a and b. This is the integral from 0 to x
- of (1/B(a,b))*f(t)) where f(t) = t**(a-1) * (1-t)**(b-1)
-
-
- Arguments
-
-
- X --> Upper limit of integration.
- X is DOUBLE PRECISION
-
- Y --> 1 - X.
- Y is DOUBLE PRECISION
-
- A --> First parameter of the beta distribution.
- A is DOUBLE PRECISION
-
- B --> Second parameter of the beta distribution.
- B is DOUBLE PRECISION
-
- CUM <-- Cumulative incomplete beta distribution.
- CUM is DOUBLE PRECISION
-
- CCUM <-- Compliment of Cumulative incomplete beta distribution.
- CCUM is DOUBLE PRECISION
-
-
- Method
-
-
- Calls the routine BRATIO.
-
- References
-
- Didonato, Armido R. and Morris, Alfred H. Jr. (1992) Algorithim
- 708 Significant Digit Computation of the Incomplete Beta Function
- Ratios. ACM ToMS, Vol.18, No. 3, Sept. 1992, 360-373.
-
-**********************************************************************
-*/
-{
-static int ierr;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*x <= 0.0e0)) goto S10;
- *cum = 0.0e0;
- *ccum = 1.0e0;
- return;
-S10:
- if(!(*y <= 0.0e0)) goto S20;
- *cum = 1.0e0;
- *ccum = 0.0e0;
- return;
-S20:
- bratio(a,b,x,y,cum,ccum,&ierr);
-/*
- Call bratio routine
-*/
- return;
-}
-void cumbin(double *s,double *xn,double *pr,double *ompr,
- double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumbin(double *s,double *xn,double *pr,double *ompr,
- double *cum,double *ccum)
-
- CUmulative BINomial distribution
-
-
- Function
-
-
- Returns the probability of 0 to S successes in XN binomial
- trials, each of which has a probability of success, PBIN.
-
-
- Arguments
-
-
- S --> The upper limit of cumulation of the binomial distribution.
- S is DOUBLE PRECISION
-
- XN --> The number of binomial trials.
- XN is DOUBLE PRECISIO
-
- PBIN --> The probability of success in each binomial trial.
- PBIN is DOUBLE PRECIS
-
- OMPR --> 1 - PBIN
- OMPR is DOUBLE PRECIS
-
- CUM <-- Cumulative binomial distribution.
- CUM is DOUBLE PRECISI
-
- CCUM <-- Compliment of Cumulative binomial distribution.
- CCUM is DOUBLE PRECIS
-
-
- Method
-
-
- Formula 26.5.24 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the binomial
- distribution to the cumulative beta distribution.
-
-**********************************************************************
-*/
-{
-static double T1,T2;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*s < *xn)) goto S10;
- T1 = *s+1.0e0;
- T2 = *xn-*s;
- cumbet(pr,ompr,&T1,&T2,ccum,cum);
- goto S20;
-S10:
- *cum = 1.0e0;
- *ccum = 0.0e0;
-S20:
- return;
-}
-void cumchi(double *x,double *df,double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumchi(double *x,double *df,double *cum,double *ccum)
- CUMulative of the CHi-square distribution
-
-
- Function
-
-
- Calculates the cumulative chi-square distribution.
-
-
- Arguments
-
-
- X --> Upper limit of integration of the
- chi-square distribution.
- X is DOUBLE PRECISION
-
- DF --> Degrees of freedom of the
- chi-square distribution.
- DF is DOUBLE PRECISION
-
- CUM <-- Cumulative chi-square distribution.
- CUM is DOUBLE PRECISIO
-
- CCUM <-- Compliment of Cumulative chi-square distribution.
- CCUM is DOUBLE PRECISI
-
-
- Method
-
-
- Calls incomplete gamma function (CUMGAM)
-
-**********************************************************************
-*/
-{
-static double a,xx;
-/*
- ..
- .. Executable Statements ..
-*/
- a = *df*0.5e0;
- xx = *x*0.5e0;
- cumgam(&xx,&a,cum,ccum);
- return;
-}
-void cumchn(double *x,double *df,double *pnonc,double *cum,
- double *ccum)
-/**********************************************************************
-
- void cumchn(double *x,double *df,double *pnonc,double *cum,
- double *ccum)
-
- CUMulative of the Non-central CHi-square distribution
-
- Function
-
- Calculates the cumulative non-central chi-square
- distribution, i.e., the probability that a random variable
- which follows the non-central chi-square distribution, with
- non-centrality parameter PNONC and continuous degrees of
- freedom DF, is less than or equal to X.
-
- Arguments
-
- X --> Upper limit of integration of the non-central
- chi-square distribution.
-
- DF --> Degrees of freedom of the non-central
- chi-square distribution.
-
- PNONC --> Non-centrality parameter of the non-central
- chi-square distribution.
-
- CUM <-- Cumulative non-central chi-square distribution.
-
- CCUM <-- Compliment of Cumulative non-central chi-square distribut
-
-
- Method
-
- Uses formula 26.4.25 of Abramowitz and Stegun, Handbook of
- Mathematical Functions, US NBS (1966) to calculate the
- non-central chi-square.
-
- Variables
-
- EPS --- Convergence criterion. The sum stops when a
- term is less than EPS*SUM.
-
- CCUM <-- Compliment of Cumulative non-central
- chi-square distribution.
-
-**********************************************************************/
-{
-#define dg(i) (*df + 2.0e0 * (double)(i))
-#define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps * sum)
-static double eps = 1.0e-5;
-static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
- sumadj,term,wt,xnonc;
-static int i,icent;
-static double T1,T2,T3;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*x <= 0.0e0)) goto S10;
- *cum = 0.0e0;
- *ccum = 1.0e0;
- return;
-S10:
- if(!(*pnonc <= 1.0e-10 )) goto S20;
-/*
- When non-centrality parameter is (essentially) zero,
- use cumulative chi-square distribution
-*/
- cumchi(x,df,cum,ccum);
- return;
-S20:
- xnonc = *pnonc / 2.0e0;
-/*
-***********************************************************************
- The following code calcualtes the weight, chi-square, and
- adjustment term for the central term in the infinite series.
- The central term is the one in which the poisson weight is
- greatest. The adjustment term is the amount that must
- be subtracted from the chi-square to move up two degrees
- of freedom.
-***********************************************************************
-*/
- icent = fifidint(xnonc);
- if(icent == 0) icent = 1;
- chid2 = *x / 2.0e0;
-/*
- Calculate central weight term
-*/
- T1 = (double)(icent + 1);
- lfact = alngam(&T1);
- lcntwt = -xnonc + (double)icent * log(xnonc) - lfact;
- centwt = exp(lcntwt);
-/*
- Calculate central chi-square
-*/
- T2 = dg(icent);
- cumchi(x,&T2,&pcent,ccum);
-/*
- Calculate central adjustment term
-*/
- dfd2 = dg(icent) / 2.0e0;
- T3 = 1.0e0 + dfd2;
- lfact = alngam(&T3);
- lcntaj = dfd2 * log(chid2) - chid2 - lfact;
- centaj = exp(lcntaj);
- sum = centwt * pcent;
-/*
-***********************************************************************
- Sum backwards from the central term towards zero.
- Quit whenever either
- (1) the zero term is reached, or
- (2) the term gets small relative to the sum
-***********************************************************************
-*/
- sumadj = 0.0e0;
- adj = centaj;
- wt = centwt;
- i = icent;
- goto S40;
-S30:
- if(qsmall(term) || i == 0) goto S50;
-S40:
- dfd2 = dg(i) / 2.0e0;
-/*
- Adjust chi-square for two fewer degrees of freedom.
- The adjusted value ends up in PTERM.
-*/
- adj = adj * dfd2 / chid2;
- sumadj += adj;
- pterm = pcent + sumadj;
-/*
- Adjust poisson weight for J decreased by one
-*/
- wt *= ((double)i / xnonc);
- term = wt * pterm;
- sum += term;
- i -= 1;
- goto S30;
-S50:
-/*
-***********************************************************************
- Now sum forward from the central term towards infinity.
- Quit when either
- (1) the term gets small relative to the sum, or
-***********************************************************************
-*/
- sumadj = adj = centaj;
- wt = centwt;
- i = icent;
- goto S70;
-S60:
- if(qsmall(term)) goto S80;
-S70:
-/*
- Update weights for next higher J
-*/
- wt *= (xnonc / (double)(i + 1));
-/*
- Calculate PTERM and add term to sum
-*/
- pterm = pcent - sumadj;
- term = wt * pterm;
- sum += term;
-/*
- Update adjustment term for DF for next iteration
-*/
- i += 1;
- dfd2 = dg(i) / 2.0e0;
- adj = adj * chid2 / dfd2;
- sumadj += adj;
- goto S60;
-S80:
- *cum = sum;
- *ccum = 0.5e0 + (0.5e0 - *cum);
- return;
-#undef dg
-#undef qsmall
-}
-void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
- CUMulative F distribution
-
-
- Function
-
-
- Computes the integral from 0 to F of the f-density with DFN
- and DFD degrees of freedom.
-
-
- Arguments
-
-
- F --> Upper limit of integration of the f-density.
- F is DOUBLE PRECISION
-
- DFN --> Degrees of freedom of the numerator sum of squares.
- DFN is DOUBLE PRECISI
-
- DFD --> Degrees of freedom of the denominator sum of squares.
- DFD is DOUBLE PRECISI
-
- CUM <-- Cumulative f distribution.
- CUM is DOUBLE PRECISI
-
- CCUM <-- Compliment of Cumulative f distribution.
- CCUM is DOUBLE PRECIS
-
-
- Method
-
-
- Formula 26.5.28 of Abramowitz and Stegun is used to reduce
- the cumulative F to a cumulative beta distribution.
-
-
- Note
-
-
- If F is less than or equal to 0, 0 is returned.
-
-**********************************************************************
-*/
-{
-#define half 0.5e0
-#define done 1.0e0
-static double dsum,prod,xx,yy;
-static int ierr;
-static double T1,T2;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*f <= 0.0e0)) goto S10;
- *cum = 0.0e0;
- *ccum = 1.0e0;
- return;
-S10:
- prod = *dfn**f;
-/*
- XX is such that the incomplete beta with parameters
- DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM
- YY is 1 - XX
- Calculate the smaller of XX and YY accurately
-*/
- dsum = *dfd+prod;
- xx = *dfd/dsum;
- if(xx > half) {
- yy = prod/dsum;
- xx = done-yy;
- }
- else yy = done-xx;
- T1 = *dfd*half;
- T2 = *dfn*half;
- bratio(&T1,&T2,&xx,&yy,ccum,cum,&ierr);
- return;
-#undef half
-#undef done
-}
-void cumfnc(double *f,double *dfn,double *dfd,double *pnonc,
- double *cum,double *ccum)
-/*
-**********************************************************************
-
- F -NON- -C-ENTRAL F DISTRIBUTION
-
-
-
- Function
-
-
- COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
- DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC
-
-
- Arguments
-
-
- X --> UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION
-
- DFN --> DEGREES OF FREEDOM OF NUMERATOR
-
- DFD --> DEGREES OF FREEDOM OF DENOMINATOR
-
- PNONC --> NONCENTRALITY PARAMETER.
-
- CUM <-- CUMULATIVE NONCENTRAL F DISTRIBUTION
-
- CCUM <-- COMPLIMENT OF CUMMULATIVE
-
-
- Method
-
-
- USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
- SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
- (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
- THE CONVERGENCE CRITERION IS MET.
-
- FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
- BY FORMULA 26.5.16.
-
-
- REFERENCE
-
-
- HANDBOOD OF MATHEMATICAL FUNCTIONS
- EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
- NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
- MARCH 1965
- P 947, EQUATIONS 26.6.17, 26.6.18
-
-
- Note
-
-
- THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
- TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20). EPS IS
- SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED.
-
-**********************************************************************
-*/
-{
-#define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
-#define half 0.5e0
-#define done 1.0e0
-static double eps = 1.0e-4;
-static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
- upterm,xmult,xnonc;
-static int i,icent,ierr;
-static double T1,T2,T3,T4,T5,T6;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*f <= 0.0e0)) goto S10;
- *cum = 0.0e0;
- *ccum = 1.0e0;
- return;
-S10:
- if(!(*pnonc < 1.0e-10)) goto S20;
-/*
- Handle case in which the non-centrality parameter is
- (essentially) zero.
-*/
- cumf(f,dfn,dfd,cum,ccum);
- return;
-S20:
- xnonc = *pnonc/2.0e0;
-/*
- Calculate the central term of the poisson weighting factor.
-*/
- icent = (long)(xnonc);
- if(icent == 0) icent = 1;
-/*
- Compute central weight term
-*/
- T1 = (double)(icent+1);
- centwt = exp(-xnonc+(double)icent*log(xnonc)-alngam(&T1));
-/*
- Compute central incomplete beta term
- Assure that minimum of arg to beta and 1 - arg is computed
- accurately.
-*/
- prod = *dfn**f;
- dsum = *dfd+prod;
- yy = *dfd/dsum;
- if(yy > half) {
- xx = prod/dsum;
- yy = done-xx;
- }
- else xx = done-yy;
- T2 = *dfn*half+(double)icent;
- T3 = *dfd*half;
- bratio(&T2,&T3,&xx,&yy,&betdn,&dummy,&ierr);
- adn = *dfn/2.0e0+(double)icent;
- aup = adn;
- b = *dfd/2.0e0;
- betup = betdn;
- sum = centwt*betdn;
-/*
- Now sum terms backward from icent until convergence or all done
-*/
- xmult = centwt;
- i = icent;
- T4 = adn+b;
- T5 = adn+1.0e0;
- dnterm = exp(alngam(&T4)-alngam(&T5)-alngam(&b)+adn*log(xx)+b*log(yy));
-S30:
- if(qsmall(xmult*betdn) || i <= 0) goto S40;
- xmult *= ((double)i/xnonc);
- i -= 1;
- adn -= 1.0;
- dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
- betdn += dnterm;
- sum += (xmult*betdn);
- goto S30;
-S40:
- i = icent+1;
-/*
- Now sum forwards until convergence
-*/
- xmult = centwt;
- if(aup-1.0+b == 0) upterm = exp(-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+
- b*log(yy));
- else {
- T6 = aup-1.0+b;
- upterm = exp(alngam(&T6)-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+b*
- log(yy));
- }
- goto S60;
-S50:
- if(qsmall(xmult*betup)) goto S70;
-S60:
- xmult *= (xnonc/(double)i);
- i += 1;
- aup += 1.0;
- upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
- betup -= upterm;
- sum += (xmult*betup);
- goto S50;
-S70:
- *cum = sum;
- *ccum = 0.5e0+(0.5e0-*cum);
- return;
-#undef qsmall
-#undef half
-#undef done
-}
-void cumgam(double *x,double *a,double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumgam(double *x,double *a,double *cum,double *ccum)
- Double precision cUMulative incomplete GAMma distribution
-
-
- Function
-
-
- Computes the cumulative of the incomplete gamma
- distribution, i.e., the integral from 0 to X of
- (1/GAM(A))*EXP(-T)*T**(A-1) DT
- where GAM(A) is the complete gamma function of A, i.e.,
- GAM(A) = integral from 0 to infinity of
- EXP(-T)*T**(A-1) DT
-
-
- Arguments
-
-
- X --> The upper limit of integration of the incomplete gamma.
- X is DOUBLE PRECISION
-
- A --> The shape parameter of the incomplete gamma.
- A is DOUBLE PRECISION
-
- CUM <-- Cumulative incomplete gamma distribution.
- CUM is DOUBLE PRECISION
-
- CCUM <-- Compliment of Cumulative incomplete gamma distribution.
- CCUM is DOUBLE PRECISIO
-
-
- Method
-
-
- Calls the routine GRATIO.
-
-**********************************************************************
-*/
-{
-static int K1 = 0;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*x <= 0.0e0)) goto S10;
- *cum = 0.0e0;
- *ccum = 1.0e0;
- return;
-S10:
- gratio(a,x,cum,ccum,&K1);
-/*
- Call gratio routine
-*/
- return;
-}
-void cumnbn(double *s,double *xn,double *pr,double *ompr,
- double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumnbn(double *s,double *xn,double *pr,double *ompr,
- double *cum,double *ccum)
-
- CUmulative Negative BINomial distribution
-
-
- Function
-
-
- Returns the probability that it there will be S or fewer failures
- before there are XN successes, with each binomial trial having
- a probability of success PR.
-
- Prob(# failures = S | XN successes, PR) =
- ( XN + S - 1 )
- ( ) * PR^XN * (1-PR)^S
- ( S )
-
-
- Arguments
-
-
- S --> The number of failures
- S is DOUBLE PRECISION
-
- XN --> The number of successes
- XN is DOUBLE PRECISIO
-
- PR --> The probability of success in each binomial trial.
- PR is DOUBLE PRECISIO
-
- OMPR --> 1 - PR
- OMPR is DOUBLE PRECIS
-
- CUM <-- Cumulative negative binomial distribution.
- CUM is DOUBLE PRECISI
-
- CCUM <-- Compliment of Cumulative negative binomial distribution.
- CCUM is DOUBLE PRECIS
-
-
- Method
-
-
- Formula 26.5.26 of Abramowitz and Stegun, Handbook of
- Mathematical Functions (1966) is used to reduce the negative
- binomial distribution to the cumulative beta distribution.
-
-**********************************************************************
-*/
-{
-static double T1;
-/*
- ..
- .. Executable Statements ..
-*/
- T1 = *s+1.e0;
- cumbet(pr,ompr,xn,&T1,cum,ccum);
- return;
-}
-void cumnor(double *arg,double *result,double *ccum)
-/*
-**********************************************************************
-
- void cumnor(double *arg,double *result,double *ccum)
-
-
- Function
-
-
- Computes the cumulative of the normal distribution, i.e.,
- the integral from -infinity to x of
- (1/sqrt(2*pi)) exp(-u*u/2) du
-
- X --> Upper limit of integration.
- X is DOUBLE PRECISION
-
- RESULT <-- Cumulative normal distribution.
- RESULT is DOUBLE PRECISION
-
- CCUM <-- Compliment of Cumulative normal distribution.
- CCUM is DOUBLE PRECISION
-
- Renaming of function ANORM from:
-
- Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
- Package of Special Function Routines and Test Drivers"
- acm Transactions on Mathematical Software. 19, 22-32.
-
- with slight modifications to return ccum and to deal with
- machine constants.
-
-**********************************************************************
- Original Comments:
-------------------------------------------------------------------
-
- This function evaluates the normal distribution function:
-
- / x
- 1 | -t*t/2
- P(x) = ----------- | e dt
- sqrt(2 pi) |
- /-oo
-
- The main computation evaluates near-minimax approximations
- derived from those in "Rational Chebyshev approximations for
- the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
- This transportable program uses rational functions that
- theoretically approximate the normal distribution function to
- at least 18 significant decimal digits. The accuracy achieved
- depends on the arithmetic system, the compiler, the intrinsic
- functions, and proper selection of the machine-dependent
- constants.
-
-*******************************************************************
-*******************************************************************
-
- Explanation of machine-dependent constants.
-
- MIN = smallest machine representable number.
-
- EPS = argument below which anorm(x) may be represented by
- 0.5 and above which x*x will not underflow.
- A conservative value is the largest machine number X
- such that 1.0 + X = 1.0 to machine precision.
-*******************************************************************
-*******************************************************************
-
- Error returns
-
- The program returns ANORM = 0 for ARG .LE. XLOW.
-
-
- Intrinsic functions required are:
-
- ABS, AINT, EXP
-
-
- Author: W. J. Cody
- Mathematics and Computer Science Division
- Argonne National Laboratory
- Argonne, IL 60439
-
- Latest modification: March 15, 1992
-
-------------------------------------------------------------------
-*/
-{
-static double a[5] = {
- 2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
- 1.8154981253343561249e04,6.5682337918207449113e-2
-};
-static double b[4] = {
- 4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
- 4.5507789335026729956e04
-};
-static double c[9] = {
- 3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
- 5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
- 1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
-};
-static double d[8] = {
- 2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
- 6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
- 3.8912003286093271411e04,1.9685429676859990727e04
-};
-static double half = 0.5e0;
-static double p[6] = {
- 2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
- 1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
-};
-static double one = 1.0e0;
-static double q[5] = {
- 1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
- 3.78239633202758244e-3,7.29751555083966205e-5
-};
-static double sixten = 1.60e0;
-static double sqrpi = 3.9894228040143267794e-1;
-static double thrsh = 0.66291e0;
-static double root32 = 5.656854248e0;
-static double zero = 0.0e0;
-static int K1 = 1;
-static int K2 = 2;
-static int i;
-static double del,eps,temp,x,xden,xnum,y,xsq,min;
-/*
-------------------------------------------------------------------
- Machine dependent constants
-------------------------------------------------------------------
-*/
- eps = spmpar(&K1)*0.5e0;
- min = spmpar(&K2);
- x = *arg;
- y = fabs(x);
- if(y <= thrsh) {
-/*
-------------------------------------------------------------------
- Evaluate anorm for |X| <= 0.66291
-------------------------------------------------------------------
-*/
- xsq = zero;
- if(y > eps) xsq = x*x;
- xnum = a[4]*xsq;
- xden = xsq;
- for(i=0; i<3; i++) {
- xnum = (xnum+a[i])*xsq;
- xden = (xden+b[i])*xsq;
- }
- *result = x*(xnum+a[3])/(xden+b[3]);
- temp = *result;
- *result = half+temp;
- *ccum = half-temp;
- }
-/*
-------------------------------------------------------------------
- Evaluate anorm for 0.66291 <= |X| <= sqrt(32)
-------------------------------------------------------------------
-*/
- else if(y <= root32) {
- xnum = c[8]*y;
- xden = y;
- for(i=0; i<7; i++) {
- xnum = (xnum+c[i])*y;
- xden = (xden+d[i])*y;
- }
- *result = (xnum+c[7])/(xden+d[7]);
- xsq = fifdint(y*sixten)/sixten;
- del = (y-xsq)*(y+xsq);
- *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
- *ccum = one-*result;
- if(x > zero) {
- temp = *result;
- *result = *ccum;
- *ccum = temp;
- }
- }
-/*
-------------------------------------------------------------------
- Evaluate anorm for |X| > sqrt(32)
-------------------------------------------------------------------
-*/
- else {
- *result = zero;
- xsq = one/(x*x);
- xnum = p[5]*xsq;
- xden = xsq;
- for(i=0; i<4; i++) {
- xnum = (xnum+p[i])*xsq;
- xden = (xden+q[i])*xsq;
- }
- *result = xsq*(xnum+p[4])/(xden+q[4]);
- *result = (sqrpi-*result)/y;
- xsq = fifdint(x*sixten)/sixten;
- del = (x-xsq)*(x+xsq);
- *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
- *ccum = one-*result;
- if(x > zero) {
- temp = *result;
- *result = *ccum;
- *ccum = temp;
- }
- }
- if(*result < min) *result = 0.0e0;
-/*
-------------------------------------------------------------------
- Fix up for negative argument, erf, etc.
-------------------------------------------------------------------
-----------Last card of ANORM ----------
-*/
- if(*ccum < min) *ccum = 0.0e0;
-}
-void cumpoi(double *s,double *xlam,double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumpoi(double *s,double *xlam,double *cum,double *ccum)
- CUMulative POIsson distribution
-
-
- Function
-
-
- Returns the probability of S or fewer events in a Poisson
- distribution with mean XLAM.
-
-
- Arguments
-
-
- S --> Upper limit of cumulation of the Poisson.
- S is DOUBLE PRECISION
-
- XLAM --> Mean of the Poisson distribution.
- XLAM is DOUBLE PRECIS
-
- CUM <-- Cumulative poisson distribution.
- CUM is DOUBLE PRECISION
-
- CCUM <-- Compliment of Cumulative poisson distribution.
- CCUM is DOUBLE PRECIS
-
-
- Method
-
-
- Uses formula 26.4.21 of Abramowitz and Stegun, Handbook of
- Mathematical Functions to reduce the cumulative Poisson to
- the cumulative chi-square distribution.
-
-**********************************************************************
-*/
-{
-static double chi,df;
-/*
- ..
- .. Executable Statements ..
-*/
- df = 2.0e0*(*s+1.0e0);
- chi = 2.0e0**xlam;
- cumchi(&chi,&df,ccum,cum);
- return;
-}
-void cumt(double *t,double *df,double *cum,double *ccum)
-/*
-**********************************************************************
-
- void cumt(double *t,double *df,double *cum,double *ccum)
- CUMulative T-distribution
-
-
- Function
-
-
- Computes the integral from -infinity to T of the t-density.
-
-
- Arguments
-
-
- T --> Upper limit of integration of the t-density.
- T is DOUBLE PRECISION
-
- DF --> Degrees of freedom of the t-distribution.
- DF is DOUBLE PRECISIO
-
- CUM <-- Cumulative t-distribution.
- CCUM is DOUBLE PRECIS
-
- CCUM <-- Compliment of Cumulative t-distribution.
- CCUM is DOUBLE PRECIS
-
-
- Method
-
-
- Formula 26.5.27 of Abramowitz and Stegun, Handbook of
- Mathematical Functions is used to reduce the t-distribution
- to an incomplete beta.
-
-**********************************************************************
-*/
-{
-static double K2 = 0.5e0;
-static double xx,a,oma,tt,yy,dfptt,T1;
-/*
- ..
- .. Executable Statements ..
-*/
- tt = *t**t;
- dfptt = *df+tt;
- xx = *df/dfptt;
- yy = tt/dfptt;
- T1 = 0.5e0**df;
- cumbet(&xx,&yy,&T1,&K2,&a,&oma);
- if(!(*t <= 0.0e0)) goto S10;
- *cum = 0.5e0*a;
- *ccum = oma+*cum;
- goto S20;
-S10:
- *ccum = 0.5e0*a;
- *cum = oma+*ccum;
-S20:
- return;
-}
-void cumtnc(double *t,double *df,double *pnonc,double *cum,
- double *ccum)
-/**********************************************************************
-
- void cumtnc(double *t,double *df,double *pnonc,double *cum,
- double *ccum)
-
- CUMulative Non-Central T-distribution
-
-
- Function
-
-
- Computes the integral from -infinity to T of the non-central
- t-density.
-
-
- Arguments
-
-
- T --> Upper limit of integration of the non-central t-density.
-
- DF --> Degrees of freedom of the non-central t-distribution.
-
- PNONC --> Non-centrality parameter of the non-central t distibutio
-
- CUM <-- Cumulative t-distribution.
-
- CCUM <-- Compliment of Cumulative t-distribution.
-
-
- Method
-
- Upper tail of the cumulative noncentral t using
- formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuous
- Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
-
- This implementation starts the calculation at i = lambda,
- which is near the largest Di. It then sums forward and backward.
-**********************************************************************/
-{
-#define one 1.0e0
-#define zero 0.0e0
-#define half 0.5e0
-#define two 2.0e0
-#define onep5 1.5e0
-#define conv 1.0e-7
-#define tiny 1.0e-10
-static double alghdf,b,bb,bbcent,bcent,cent,d,dcent,dpnonc,dum1,dum2,e,ecent,
- halfdf,lambda,lnomx,lnx,omx,pnonc2,s,scent,ss,sscent,t2,term,tt,twoi,x,xi,
- xlnd,xlne;
-static int ierr;
-static unsigned long qrevs;
-static double T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- Case pnonc essentially zero
-*/
- if(fabs(*pnonc) <= tiny) {
- cumt(t,df,cum,ccum);
- return;
- }
- qrevs = *t < zero;
- if(qrevs) {
- tt = -*t;
- dpnonc = -*pnonc;
- }
- else {
- tt = *t;
- dpnonc = *pnonc;
- }
- pnonc2 = dpnonc * dpnonc;
- t2 = tt * tt;
- if(fabs(tt) <= tiny) {
- T1 = -*pnonc;
- cumnor(&T1,cum,ccum);
- return;
- }
- lambda = half * pnonc2;
- x = *df / (*df + t2);
- omx = one - x;
- lnx = log(x);
- lnomx = log(omx);
- halfdf = half * *df;
- alghdf = gamln(&halfdf);
-/*
- ******************** Case i = lambda
-*/
- cent = fifidint(lambda);
- if(cent < one) cent = one;
-/*
- Compute d=T(2i) in log space and offset by exp(-lambda)
-*/
- T2 = cent + one;
- xlnd = cent * log(lambda) - gamln(&T2) - lambda;
- dcent = exp(xlnd);
-/*
- Compute e=t(2i+1) in log space offset by exp(-lambda)
-*/
- T3 = cent + onep5;
- xlne = (cent + half) * log(lambda) - gamln(&T3) - lambda;
- ecent = exp(xlne);
- if(dpnonc < zero) ecent = -ecent;
-/*
- Compute bcent=B(2*cent)
-*/
- T4 = cent + half;
- bratio(&halfdf,&T4,&x,&omx,&bcent,&dum1,&ierr);
-/*
- compute bbcent=B(2*cent+1)
-*/
- T5 = cent + one;
- bratio(&halfdf,&T5,&x,&omx,&bbcent,&dum2,&ierr);
-/*
- Case bcent and bbcent are essentially zero
- Thus t is effectively infinite
-*/
- if(bcent + bbcent < tiny) {
- if(qrevs) {
- *cum = zero;
- *ccum = one;
- }
- else {
- *cum = one;
- *ccum = zero;
- }
- return;
- }
-/*
- Case bcent and bbcent are essentially one
- Thus t is effectively zero
-*/
- if(dum1 + dum2 < tiny) {
- T6 = -*pnonc;
- cumnor(&T6,cum,ccum);
- return;
- }
-/*
- First term in ccum is D*B + E*BB
-*/
- *ccum = dcent * bcent + ecent * bbcent;
-/*
- compute s(cent) = B(2*(cent+1)) - B(2*cent))
-*/
- T7 = halfdf + cent + half;
- T8 = cent + onep5;
- scent = gamln(&T7) - gamln(&T8) - alghdf + halfdf * lnx + (cent + half) *
- lnomx;
- scent = exp(scent);
-/*
- compute ss(cent) = B(2*cent+3) - B(2*cent+1)
-*/
- T9 = halfdf + cent + one;
- T10 = cent + two;
- sscent = gamln(&T9) - gamln(&T10) - alghdf + halfdf * lnx + (cent + one) *
- lnomx;
- sscent = exp(sscent);
-/*
- ******************** Sum Forward
-*/
- xi = cent + one;
- twoi = two * xi;
- d = dcent;
- e = ecent;
- b = bcent;
- bb = bbcent;
- s = scent;
- ss = sscent;
-S10:
- b += s;
- bb += ss;
- d = lambda / xi * d;
- e = lambda / (xi + half) * e;
- term = d * b + e * bb;
- *ccum += term;
- s = s * omx * (*df + twoi - one) / (twoi + one);
- ss = ss * omx * (*df + twoi) / (twoi + two);
- xi += one;
- twoi = two * xi;
- if(fabs(term) > conv * *ccum) goto S10;
-/*
- ******************** Sum Backward
-*/
- xi = cent;
- twoi = two * xi;
- d = dcent;
- e = ecent;
- b = bcent;
- bb = bbcent;
- s = scent * (one + twoi) / ((*df + twoi - one) * omx);
- ss = sscent * (two + twoi) / ((*df + twoi) * omx);
-S20:
- b -= s;
- bb -= ss;
- d *= (xi / lambda);
- e *= ((xi + half) / lambda);
- term = d * b + e * bb;
- *ccum += term;
- xi -= one;
- if(xi < half) goto S30;
- twoi = two * xi;
- s = s * (one + twoi) / ((*df + twoi - one) * omx);
- ss = ss * (two + twoi) / ((*df + twoi) * omx);
- if(fabs(term) > conv * *ccum) goto S20;
-S30:
- if(qrevs) {
- *cum = half * *ccum;
- *ccum = one - *cum;
- }
- else {
- *ccum = half * *ccum;
- *cum = one - *ccum;
- }
-/*
- Due to roundoff error the answer may not lie between zero and one
- Force it to do so
-*/
- *cum = fifdmax1(fifdmin1(*cum,one),zero);
- *ccum = fifdmax1(fifdmin1(*ccum,one),zero);
- return;
-#undef one
-#undef zero
-#undef half
-#undef two
-#undef onep5
-#undef conv
-#undef tiny
-}
-double devlpl(double a[],int *n,double *x)
-/*
-**********************************************************************
-
- double devlpl(double a[],int *n,double *x)
- Double precision EVALuate a PoLynomial at X
-
-
- Function
-
-
- returns
- A(1) + A(2)*X + ... + A(N)*X**(N-1)
-
-
- Arguments
-
-
- A --> Array of coefficients of the polynomial.
- A is DOUBLE PRECISION(N)
-
- N --> Length of A, also degree of polynomial - 1.
- N is INTEGER
-
- X --> Point at which the polynomial is to be evaluated.
- X is DOUBLE PRECISION
-
-**********************************************************************
-*/
-{
-static double devlpl,term;
-static int i;
-/*
- ..
- .. Executable Statements ..
-*/
- term = a[*n-1];
- for(i= *n-1-1; i>=0; i--) term = a[i]+term**x;
- devlpl = term;
- return devlpl;
-}
-double dinvnr(double *p,double *q)
-/*
-**********************************************************************
-
- double dinvnr(double *p,double *q)
- Double precision NoRmal distribution INVerse
-
-
- Function
-
-
- Returns X such that CUMNOR(X) = P, i.e., the integral from -
- infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
-
-
- Arguments
-
-
- P --> The probability whose normal deviate is sought.
- P is DOUBLE PRECISION
-
- Q --> 1-P
- P is DOUBLE PRECISION
-
-
- Method
-
-
- The rational function on page 95 of Kennedy and Gentle,
- Statistical Computing, Marcel Dekker, NY , 1980 is used as a start
- value for the Newton method of finding roots.
-
-
- Note
-
-
- If P or Q .lt. machine EPS returns +/- DINVNR(EPS)
-
-**********************************************************************
-*/
-{
-#define maxit 100
-#define eps 1.0e-13
-#define r2pi 0.3989422804014326e0
-#define nhalf -0.5e0
-#define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
-static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
-static int i;
-static unsigned long qporq;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- FIND MINIMUM OF P AND Q
-*/
- qporq = *p <= *q;
- if(!qporq) goto S10;
- pp = *p;
- goto S20;
-S10:
- pp = *q;
-S20:
-/*
- INITIALIZATION STEP
-*/
- strtx = stvaln(&pp);
- xcur = strtx;
-/*
- NEWTON INTERATIONS
-*/
- for(i=1; i<=maxit; i++) {
- cumnor(&xcur,&cum,&ccum);
- dx = (cum-pp)/dennor(xcur);
- xcur -= dx;
- if(fabs(dx/xcur) < eps) goto S40;
- }
- dinvnr = strtx;
-/*
- IF WE GET HERE, NEWTON HAS FAILED
-*/
- if(!qporq) dinvnr = -dinvnr;
- return dinvnr;
-S40:
-/*
- IF WE GET HERE, NEWTON HAS SUCCEDED
-*/
- dinvnr = xcur;
- if(!qporq) dinvnr = -dinvnr;
- return dinvnr;
-#undef maxit
-#undef eps
-#undef r2pi
-#undef nhalf
-#undef dennor
-}
-/* DEFINE DINVR */
-static void E0000(int IENTRY,int *status,double *x,double *fx,
- unsigned long *qleft,unsigned long *qhi,double *zabsst,
- double *zabsto,double *zbig,double *zrelst,
- double *zrelto,double *zsmall,double *zstpmu)
-{
-#define qxmon(zx,zy,zz) (int)((zx) <= (zy) && (zy) <= (zz))
-static double absstp,abstol,big,fbig,fsmall,relstp,reltol,small,step,stpmul,xhi,
- xlb,xlo,xsave,xub,yy;
-static int i99999;
-static unsigned long qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup;
- switch(IENTRY){case 0: goto DINVR; case 1: goto DSTINV;}
-DINVR:
- if(*status > 0) goto S310;
- qcond = !qxmon(small,*x,big);
- if(qcond)
- {
- *status = -999;
- return;
- }
-
- xsave = *x;
-/*
- See that SMALL and BIG bound the zero and set QINCR
-*/
- *x = small;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 1;
- goto S300;
-S10:
- fsmall = *fx;
- *x = big;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 2;
- goto S300;
-S20:
- fbig = *fx;
- qincr = fbig > fsmall;
- if(!qincr) goto S50;
- if(fsmall <= 0.0e0) goto S30;
- *status = -1;
- *qleft = *qhi = 1;
- return;
-S30:
- if(fbig >= 0.0e0) goto S40;
- *status = -1;
- *qleft = *qhi = 0;
- return;
-S40:
- goto S80;
-S50:
- if(fsmall >= 0.0e0) goto S60;
- *status = -1;
- *qleft = 1;
- *qhi = 0;
- return;
-S60:
- if(fbig <= 0.0e0) goto S70;
- *status = -1;
- *qleft = 0;
- *qhi = 1;
- return;
-S80:
-S70:
- *x = xsave;
- step = fifdmax1(absstp,relstp*fabs(*x));
-/*
- YY = F(X) - Y
- GET-FUNCTION-VALUE
-*/
- i99999 = 3;
- goto S300;
-S90:
- yy = *fx;
- if(!(yy == 0.0e0)) goto S100;
- *status = 0;
- qok = 1;
- return;
-S100:
- qup = (qincr && yy < 0.0e0) || (!qincr && yy > 0.0e0);
-/*
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- HANDLE CASE IN WHICH WE MUST STEP HIGHER
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-*/
- if(!qup) goto S170;
- xlb = xsave;
- xub = fifdmin1(xlb+step,big);
- goto S120;
-S110:
- if(qcond) goto S150;
-S120:
-/*
- YY = F(XUB) - Y
-*/
- *x = xub;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 4;
- goto S300;
-S130:
- yy = *fx;
- qbdd = (qincr && yy >= 0.0e0) || (!qincr && yy <= 0.0e0);
- qlim = xub >= big;
- qcond = qbdd || qlim;
- if(qcond) goto S140;
- step = stpmul*step;
- xlb = xub;
- xub = fifdmin1(xlb+step,big);
-S140:
- goto S110;
-S150:
- if(!(qlim && !qbdd)) goto S160;
- *status = -1;
- *qleft = 0;
- *qhi = !qincr;
- *x = big;
- return;
-S160:
- goto S240;
-S170:
-/*
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- HANDLE CASE IN WHICH WE MUST STEP LOWER
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-*/
- xub = xsave;
- xlb = fifdmax1(xub-step,small);
- goto S190;
-S180:
- if(qcond) goto S220;
-S190:
-/*
- YY = F(XLB) - Y
-*/
- *x = xlb;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 5;
- goto S300;
-S200:
- yy = *fx;
- qbdd = (qincr && yy <= 0.0e0) || (!qincr && yy >= 0.0e0);
- qlim = xlb <= small;
- qcond = qbdd || qlim;
- if(qcond) goto S210;
- step = stpmul*step;
- xub = xlb;
- xlb = fifdmax1(xub-step,small);
-S210:
- goto S180;
-S220:
- if(!(qlim && !qbdd)) goto S230;
- *status = -1;
- *qleft = 1;
- *qhi = qincr;
- *x = small;
- return;
-S240:
-S230:
- dstzr(&xlb,&xub,&abstol,&reltol);
-/*
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
-++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-*/
- *status = 0;
- goto S260;
-S250:
- if(!(*status == 1)) goto S290;
-S260:
- dzror(status,x,fx,&xlo,&xhi,&qdum1,&qdum2);
- if(!(*status == 1)) goto S280;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 6;
- goto S300;
-S280:
-S270:
- goto S250;
-S290:
- *x = xlo;
- *status = 0;
- return;
-DSTINV:
- small = *zsmall;
- big = *zbig;
- absstp = *zabsst;
- relstp = *zrelst;
- stpmul = *zstpmu;
- abstol = *zabsto;
- reltol = *zrelto;
- return;
-S300:
-/*
- TO GET-FUNCTION-VALUE
-*/
- *status = 1;
- return;
-S310:
- switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S90;case
- 4: goto S130;case 5: goto S200;case 6: goto S270;default: break;}
-#undef qxmon
-}
-void dinvr(int *status,double *x,double *fx,
- unsigned long *qleft,unsigned long *qhi)
-/*
-**********************************************************************
-
- void dinvr(int *status,double *x,double *fx,
- unsigned long *qleft,unsigned long *qhi)
-
- Double precision
- bounds the zero of the function and invokes zror
- Reverse Communication
-
-
- Function
-
-
- Bounds the function and invokes ZROR to perform the zero
- finding. STINVR must have been called before this routine
- in order to set its parameters.
-
-
- Arguments
-
-
- STATUS <--> At the beginning of a zero finding problem, STATUS
- should be set to 0 and INVR invoked. (The value
- of parameters other than X will be ignored on this cal
-
- When INVR needs the function evaluated, it will set
- STATUS to 1 and return. The value of the function
- should be set in FX and INVR again called without
- changing any of its other parameters.
-
- When INVR has finished without error, it will return
- with STATUS 0. In that case X is approximately a root
- of F(X).
-
- If INVR cannot bound the function, it returns status
- -1 and sets QLEFT and QHI.
- INTEGER STATUS
-
- X <-- The value of X at which F(X) is to be evaluated.
- DOUBLE PRECISION X
-
- FX --> The value of F(X) calculated when INVR returns with
- STATUS = 1.
- DOUBLE PRECISION FX
-
- QLEFT <-- Defined only if QMFINV returns .FALSE. In that
- case it is .TRUE. If the stepping search terminated
- unsucessfully at SMALL. If it is .FALSE. the search
- terminated unsucessfully at BIG.
- QLEFT is LOGICAL
-
- QHI <-- Defined only if QMFINV returns .FALSE. In that
- case it is .TRUE. if F(X) .GT. Y at the termination
- of the search and .FALSE. if F(X) .LT. Y at the
- termination of the search.
- QHI is LOGICAL
-
-**********************************************************************
-*/
-{
- E0000(0,status,x,fx,qleft,qhi,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
-}
-void dstinv(double *zsmall,double *zbig,double *zabsst,
- double *zrelst,double *zstpmu,double *zabsto,
- double *zrelto)
-/*
-**********************************************************************
- void dstinv(double *zsmall,double *zbig,double *zabsst,
- double *zrelst,double *zstpmu,double *zabsto,
- double *zrelto)
-
- Double Precision - SeT INverse finder - Reverse Communication
- Function
- Concise Description - Given a monotone function F finds X
- such that F(X) = Y. Uses Reverse communication -- see invr.
- This routine sets quantities needed by INVR.
- More Precise Description of INVR -
- F must be a monotone function, the results of QMFINV are
- otherwise undefined. QINCR must be .TRUE. if F is non-
- decreasing and .FALSE. if F is non-increasing.
- QMFINV will return .TRUE. if and only if F(SMALL) and
- F(BIG) bracket Y, i. e.,
- QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
- QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
- if QMFINV returns .TRUE., then the X returned satisfies
- the following condition. let
- TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
- then if QINCR is .TRUE.,
- F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
- and if QINCR is .FALSE.
- F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
- Arguments
- SMALL --> The left endpoint of the interval to be
- searched for a solution.
- SMALL is DOUBLE PRECISION
- BIG --> The right endpoint of the interval to be
- searched for a solution.
- BIG is DOUBLE PRECISION
- ABSSTP, RELSTP --> The initial step size in the search
- is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
- ABSSTP is DOUBLE PRECISION
- RELSTP is DOUBLE PRECISION
- STPMUL --> When a step doesn't bound the zero, the step
- size is multiplied by STPMUL and another step
- taken. A popular value is 2.0
- DOUBLE PRECISION STPMUL
- ABSTOL, RELTOL --> Two numbers that determine the accuracy
- of the solution. See function for a precise definition.
- ABSTOL is DOUBLE PRECISION
- RELTOL is DOUBLE PRECISION
- Method
- Compares F(X) with Y for the input value of X then uses QINCR
- to determine whether to step left or right to bound the
- desired x. the initial step size is
- MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
- Iteratively steps right or left until it bounds X.
- At each step which doesn't bound X, the step size is doubled.
- The routine is careful never to step beyond SMALL or BIG. If
- it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
- after setting QLEFT and QHI.
- If X is successfully bounded then Algorithm R of the paper
- 'Two Efficient Algorithms with Guaranteed Convergence for
- Finding a Zero of a Function' by J. C. P. Bus and
- T. J. Dekker in ACM Transactions on Mathematical
- Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
- to find the zero of the function F(X)-Y. This is routine
- QRZERO.
-**********************************************************************
-*/
-{
- int status=0;
- E0000(1,&status,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
- zstpmu);
- assert(status == 0 );
-}
-double dt1(double *p,double *q,double *df)
-/*
-**********************************************************************
-
- double dt1(double *p,double *q,double *df)
- Double precision Initalize Approximation to
- INVerse of the cumulative T distribution
-
-
- Function
-
-
- Returns the inverse of the T distribution function, i.e.,
- the integral from 0 to INVT of the T density is P. This is an
- initial approximation
-
-
- Arguments
-
-
- P --> The p-value whose inverse from the T distribution is
- desired.
- P is DOUBLE PRECISION
-
- Q --> 1-P.
- Q is DOUBLE PRECISION
-
- DF --> Degrees of freedom of the T distribution.
- DF is DOUBLE PRECISION
-
-**********************************************************************
-*/
-{
-static double coef[4][5] = {
- {1.0e0,1.0e0,0.0e0,0.0e0,0.0e0},
- {3.0e0,16.0e0,5.0e0,0.0e0,0.0e0},
- {-15.0e0,17.0e0,19.0e0,3.0e0,0.0e0},
- {-945.0e0,-1920.0e0,1482.0e0,776.0e0,79.0e0}
-};
-static double denom[4] = {
- 4.0e0,96.0e0,384.0e0,92160.0e0
-};
-static int ideg[4] = {
- 2,3,4,5
-};
-static double dt1,denpow,sum,term,x,xp,xx;
-static int i;
-/*
- ..
- .. Executable Statements ..
-*/
- x = fabs(dinvnr(p,q));
- xx = x*x;
- sum = x;
- denpow = 1.0e0;
- for(i=0; i<4; i++) {
- term = devlpl(&coef[i][0],&ideg[i],&xx)*x;
- denpow *= *df;
- sum += (term/(denpow*denom[i]));
- }
- if(!(*p >= 0.5e0)) goto S20;
- xp = sum;
- goto S30;
-S20:
- xp = -sum;
-S30:
- dt1 = xp;
- return dt1;
-}
-/* DEFINE DZROR */
-static void E0001(int IENTRY,int *status,double *x,double *fx,
- double *xlo,double *xhi,unsigned long *qleft,
- unsigned long *qhi,double *zabstl,double *zreltl,
- double *zxhi,double *zxlo)
-{
-#define ftol(zx) (0.5e0*fifdmax1(abstol,reltol*fabs((zx))))
-static double a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q,reltol,tol,w,xxhi,xxlo;
-static int ext,i99999;
-static unsigned long first,qrzero;
- switch(IENTRY){case 0: goto DZROR; case 1: goto DSTZR;}
-DZROR:
- if(*status > 0) goto S280;
- *xlo = xxlo;
- *xhi = xxhi;
- b = *x = *xlo;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 1;
- goto S270;
-S10:
- fb = *fx;
- *xlo = *xhi;
- a = *x = *xlo;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 2;
- goto S270;
-S20:
-/*
- Check that F(ZXLO) < 0 < F(ZXHI) or
- F(ZXLO) > 0 > F(ZXHI)
-*/
- if(!(fb < 0.0e0)) goto S40;
- if(!(*fx < 0.0e0)) goto S30;
- *status = -1;
- *qleft = *fx < fb;
- *qhi = 0;
- return;
-S40:
-S30:
- if(!(fb > 0.0e0)) goto S60;
- if(!(*fx > 0.0e0)) goto S50;
- *status = -1;
- *qleft = *fx > fb;
- *qhi = 1;
- return;
-S60:
-S50:
- fa = *fx;
- first = 1;
-S70:
- c = a;
- fc = fa;
- ext = 0;
-S80:
- if(!(fabs(fc) < fabs(fb))) goto S100;
- if(!(c != a)) goto S90;
- d = a;
- fd = fa;
-S90:
- a = b;
- fa = fb;
- *xlo = c;
- b = *xlo;
- fb = fc;
- c = a;
- fc = fa;
-S100:
- tol = ftol(*xlo);
- m = (c+b)*.5e0;
- mb = m-b;
- if(!(fabs(mb) > tol)) goto S240;
- if(!(ext > 3)) goto S110;
- w = mb;
- goto S190;
-S110:
- tol = fifdsign(tol,mb);
- p = (b-a)*fb;
- if(!first) goto S120;
- q = fa-fb;
- first = 0;
- goto S130;
-S120:
- fdb = (fd-fb)/(d-b);
- fda = (fd-fa)/(d-a);
- p = fda*p;
- q = fdb*fa-fda*fb;
-S130:
- if(!(p < 0.0e0)) goto S140;
- p = -p;
- q = -q;
-S140:
- if(ext == 3) p *= 2.0e0;
- if(!(p*1.0e0 == 0.0e0 || p <= q*tol)) goto S150;
- w = tol;
- goto S180;
-S150:
- if(!(p < mb*q)) goto S160;
- w = p/q;
- goto S170;
-S160:
- w = mb;
-S190:
-S180:
-S170:
- d = a;
- fd = fa;
- a = b;
- fa = fb;
- b += w;
- *xlo = b;
- *x = *xlo;
-/*
- GET-FUNCTION-VALUE
-*/
- i99999 = 3;
- goto S270;
-S200:
- fb = *fx;
- if(!(fc*fb >= 0.0e0)) goto S210;
- goto S70;
-S210:
- if(!(w == mb)) goto S220;
- ext = 0;
- goto S230;
-S220:
- ext += 1;
-S230:
- goto S80;
-S240:
- *xhi = c;
- qrzero = (fc >= 0.0e0 && fb <= 0.0e0) || (fc < 0.0e0 && fb >= 0.0e0);
- if(!qrzero) goto S250;
- *status = 0;
- goto S260;
-S250:
- *status = -1;
-S260:
- return;
-DSTZR:
- xxlo = *zxlo;
- xxhi = *zxhi;
- abstol = *zabstl;
- reltol = *zreltl;
- return;
-S270:
-/*
- TO GET-FUNCTION-VALUE
-*/
- *status = 1;
- return;
-S280:
- switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S200;
- default: break;}
-#undef ftol
-}
-void dzror(int *status,double *x,double *fx,double *xlo,
- double *xhi,unsigned long *qleft,unsigned long *qhi)
-/*
-**********************************************************************
-
- void dzror(int *status,double *x,double *fx,double *xlo,
- double *xhi,unsigned long *qleft,unsigned long *qhi)
-
- Double precision ZeRo of a function -- Reverse Communication
-
-
- Function
-
-
- Performs the zero finding. STZROR must have been called before
- this routine in order to set its parameters.
-
-
- Arguments
-
-
- STATUS <--> At the beginning of a zero finding problem, STATUS
- should be set to 0 and ZROR invoked. (The value
- of other parameters will be ignored on this call.)
-
- When ZROR needs the function evaluated, it will set
- STATUS to 1 and return. The value of the function
- should be set in FX and ZROR again called without
- changing any of its other parameters.
-
- When ZROR has finished without error, it will return
- with STATUS 0. In that case (XLO,XHI) bound the answe
-
- If ZROR finds an error (which implies that F(XLO)-Y an
- F(XHI)-Y have the same sign, it returns STATUS -1. In
- this case, XLO and XHI are undefined.
- INTEGER STATUS
-
- X <-- The value of X at which F(X) is to be evaluated.
- DOUBLE PRECISION X
-
- FX --> The value of F(X) calculated when ZROR returns with
- STATUS = 1.
- DOUBLE PRECISION FX
-
- XLO <-- When ZROR returns with STATUS = 0, XLO bounds the
- inverval in X containing the solution below.
- DOUBLE PRECISION XLO
-
- XHI <-- When ZROR returns with STATUS = 0, XHI bounds the
- inverval in X containing the solution above.
- DOUBLE PRECISION XHI
-
- QLEFT <-- .TRUE. if the stepping search terminated unsucessfully
- at XLO. If it is .FALSE. the search terminated
- unsucessfully at XHI.
- QLEFT is LOGICAL
-
- QHI <-- .TRUE. if F(X) .GT. Y at the termination of the
- search and .FALSE. if F(X) .LT. Y at the
- termination of the search.
- QHI is LOGICAL
-
-**********************************************************************
-*/
-{
- E0001(0,status,x,fx,xlo,xhi,qleft,qhi,NULL,NULL,NULL,NULL);
-}
-void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl)
-/*
-**********************************************************************
- void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl)
- Double precision SeT ZeRo finder - Reverse communication version
- Function
- Sets quantities needed by ZROR. The function of ZROR
- and the quantities set is given here.
- Concise Description - Given a function F
- find XLO such that F(XLO) = 0.
- More Precise Description -
- Input condition. F is a double precision function of a single
- double precision argument and XLO and XHI are such that
- F(XLO)*F(XHI) .LE. 0.0
- If the input condition is met, QRZERO returns .TRUE.
- and output values of XLO and XHI satisfy the following
- F(XLO)*F(XHI) .LE. 0.
- ABS(F(XLO) .LE. ABS(F(XHI)
- ABS(XLO-XHI) .LE. TOL(X)
- where
- TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
- If this algorithm does not find XLO and XHI satisfying
- these conditions then QRZERO returns .FALSE. This
- implies that the input condition was not met.
- Arguments
- XLO --> The left endpoint of the interval to be
- searched for a solution.
- XLO is DOUBLE PRECISION
- XHI --> The right endpoint of the interval to be
- for a solution.
- XHI is DOUBLE PRECISION
- ABSTOL, RELTOL --> Two numbers that determine the accuracy
- of the solution. See function for a
- precise definition.
- ABSTOL is DOUBLE PRECISION
- RELTOL is DOUBLE PRECISION
- Method
- Algorithm R of the paper 'Two Efficient Algorithms with
- Guaranteed Convergence for Finding a Zero of a Function'
- by J. C. P. Bus and T. J. Dekker in ACM Transactions on
- Mathematical Software, Volume 1, no. 4 page 330
- (Dec. '75) is employed to find the zero of F(X)-Y.
-**********************************************************************
-*/
-{
- E0001(1,NULL,NULL,NULL,NULL,NULL,NULL,NULL,zabstl,zreltl,zxhi,zxlo);
-}
-double erf1(double *x)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE REAL ERROR FUNCTION
------------------------------------------------------------------------
-*/
-{
-static double c = .564189583547756e0;
-static double a[5] = {
- .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
- .479137145607681e-01,.128379167095513e+00
-};
-static double b[3] = {
- .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
-};
-static double p[8] = {
- -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
- 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
- 4.51918953711873e+02,3.00459261020162e+02
-};
-static double q[8] = {
- 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
- 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
- 7.90950925327898e+02,3.00459260956983e+02
-};
-static double r[5] = {
- 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
- 4.65807828718470e+00,2.82094791773523e-01
-};
-static double s[4] = {
- 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
- 1.80124575948747e+01
-};
-static double erf1,ax,bot,t,top,x2;
-/*
- ..
- .. Executable Statements ..
-*/
- ax = fabs(*x);
- if(ax > 0.5e0) goto S10;
- t = *x**x;
- top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
- bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
- erf1 = *x*(top/bot);
- return erf1;
-S10:
- if(ax > 4.0e0) goto S20;
- top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
- 7];
- bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
- 7];
- erf1 = 0.5e0+(0.5e0-exp(-(*x**x))*top/bot);
- if(*x < 0.0e0) erf1 = -erf1;
- return erf1;
-S20:
- if(ax >= 5.8e0) goto S30;
- x2 = *x**x;
- t = 1.0e0/x2;
- top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
- bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
- erf1 = (c-top/(x2*bot))/ax;
- erf1 = 0.5e0+(0.5e0-exp(-x2)*erf1);
- if(*x < 0.0e0) erf1 = -erf1;
- return erf1;
-S30:
- erf1 = fifdsign(1.0e0,*x);
- return erf1;
-}
-double erfc1(int *ind,double *x)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
-
- ERFC1(IND,X) = ERFC(X) IF IND = 0
- ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE
------------------------------------------------------------------------
-*/
-{
-static double c = .564189583547756e0;
-static double a[5] = {
- .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
- .479137145607681e-01,.128379167095513e+00
-};
-static double b[3] = {
- .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
-};
-static double p[8] = {
- -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
- 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
- 4.51918953711873e+02,3.00459261020162e+02
-};
-static double q[8] = {
- 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
- 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
- 7.90950925327898e+02,3.00459260956983e+02
-};
-static double r[5] = {
- 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
- 4.65807828718470e+00,2.82094791773523e-01
-};
-static double s[4] = {
- 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
- 1.80124575948747e+01
-};
-static int K1 = 1;
-static double erfc1,ax,bot,e,t,top,w;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- ABS(X) .LE. 0.5
-*/
- ax = fabs(*x);
- if(ax > 0.5e0) goto S10;
- t = *x**x;
- top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
- bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
- erfc1 = 0.5e0+(0.5e0-*x*(top/bot));
- if(*ind != 0) erfc1 = exp(t)*erfc1;
- return erfc1;
-S10:
-/*
- 0.5 .LT. ABS(X) .LE. 4
-*/
- if(ax > 4.0e0) goto S20;
- top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
- 7];
- bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
- 7];
- erfc1 = top/bot;
- goto S40;
-S20:
-/*
- ABS(X) .GT. 4
-*/
- if(*x <= -5.6e0) goto S60;
- if(*ind != 0) goto S30;
- if(*x > 100.0e0) goto S70;
- if(*x**x > -exparg(&K1)) goto S70;
-S30:
- t = pow(1.0e0/ *x,2.0);
- top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
- bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
- erfc1 = (c-t*top/bot)/ax;
-S40:
-/*
- FINAL ASSEMBLY
-*/
- if(*ind == 0) goto S50;
- if(*x < 0.0e0) erfc1 = 2.0e0*exp(*x**x)-erfc1;
- return erfc1;
-S50:
- w = *x**x;
- t = w;
- e = w-t;
- erfc1 = (0.5e0+(0.5e0-e))*exp(-t)*erfc1;
- if(*x < 0.0e0) erfc1 = 2.0e0-erfc1;
- return erfc1;
-S60:
-/*
- LIMIT VALUE FOR LARGE NEGATIVE X
-*/
- erfc1 = 2.0e0;
- if(*ind != 0) erfc1 = 2.0e0*exp(*x**x);
- return erfc1;
-S70:
-/*
- LIMIT VALUE FOR LARGE POSITIVE X
- WHEN IND = 0
-*/
- erfc1 = 0.0e0;
- return erfc1;
-}
-double esum(int *mu,double *x)
-/*
------------------------------------------------------------------------
- EVALUATION OF EXP(MU + X)
------------------------------------------------------------------------
-*/
-{
-static double esum,w;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*x > 0.0e0) goto S10;
- if(*mu < 0) goto S20;
- w = (double)*mu+*x;
- if(w > 0.0e0) goto S20;
- esum = exp(w);
- return esum;
-S10:
- if(*mu > 0) goto S20;
- w = (double)*mu+*x;
- if(w < 0.0e0) goto S20;
- esum = exp(w);
- return esum;
-S20:
- w = *mu;
- esum = exp(w)*exp(*x);
- return esum;
-}
-double exparg(int *l)
-/*
---------------------------------------------------------------------
- IF L = 0 THEN EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
- EXP(W) CAN BE COMPUTED.
-
- IF L IS NONZERO THEN EXPARG(L) = THE LARGEST NEGATIVE W FOR
- WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.
-
- NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.
---------------------------------------------------------------------
-*/
-{
-static int K1 = 4;
-static int K2 = 9;
-static int K3 = 10;
-static double exparg,lnb;
-static int b,m;
-/*
- ..
- .. Executable Statements ..
-*/
- b = ipmpar(&K1);
- if(b != 2) goto S10;
- lnb = .69314718055995e0;
- goto S40;
-S10:
- if(b != 8) goto S20;
- lnb = 2.0794415416798e0;
- goto S40;
-S20:
- if(b != 16) goto S30;
- lnb = 2.7725887222398e0;
- goto S40;
-S30:
- lnb = log((double)b);
-S40:
- if(*l == 0) goto S50;
- m = ipmpar(&K2)-1;
- exparg = 0.99999e0*((double)m*lnb);
- return exparg;
-S50:
- m = ipmpar(&K3);
- exparg = 0.99999e0*((double)m*lnb);
- return exparg;
-}
-double fpser(double *a,double *b,double *x,double *eps)
-/*
------------------------------------------------------------------------
-
- EVALUATION OF I (A,B)
- X
-
- FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
-
------------------------------------------------------------------------
-
- SET FPSER = X**A
-*/
-{
-static int K1 = 1;
-static double fpser,an,c,s,t,tol;
-/*
- ..
- .. Executable Statements ..
-*/
- fpser = 1.0e0;
- if(*a <= 1.e-3**eps) goto S10;
- fpser = 0.0e0;
- t = *a*log(*x);
- if(t < exparg(&K1)) return fpser;
- fpser = exp(t);
-S10:
-/*
- NOTE THAT 1/B(A,B) = B
-*/
- fpser = *b/ *a*fpser;
- tol = *eps/ *a;
- an = *a+1.0e0;
- t = *x;
- s = t/an;
-S20:
- an += 1.0e0;
- t = *x*t;
- c = t/an;
- s += c;
- if(fabs(c) > tol) goto S20;
- fpser *= (1.0e0+*a*s);
- return fpser;
-}
-double gam1(double *a)
-/*
- ------------------------------------------------------------------
- COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5
- ------------------------------------------------------------------
-*/
-{
-static double s1 = .273076135303957e+00;
-static double s2 = .559398236957378e-01;
-static double p[7] = {
- .577215664901533e+00,-.409078193005776e+00,-.230975380857675e+00,
- .597275330452234e-01,.766968181649490e-02,-.514889771323592e-02,
- .589597428611429e-03
-};
-static double q[5] = {
- .100000000000000e+01,.427569613095214e+00,.158451672430138e+00,
- .261132021441447e-01,.423244297896961e-02
-};
-static double r[9] = {
- -.422784335098468e+00,-.771330383816272e+00,-.244757765222226e+00,
- .118378989872749e+00,.930357293360349e-03,-.118290993445146e-01,
- .223047661158249e-02,.266505979058923e-03,-.132674909766242e-03
-};
-static double gam1,bot,d,t,top,w,T1;
-/*
- ..
- .. Executable Statements ..
-*/
- t = *a;
- d = *a-0.5e0;
- if(d > 0.0e0) t = d-0.5e0;
- T1 = t;
- if(T1 < 0) goto S40;
- else if(T1 == 0) goto S10;
- else goto S20;
-S10:
- gam1 = 0.0e0;
- return gam1;
-S20:
- top = (((((p[6]*t+p[5])*t+p[4])*t+p[3])*t+p[2])*t+p[1])*t+p[0];
- bot = (((q[4]*t+q[3])*t+q[2])*t+q[1])*t+1.0e0;
- w = top/bot;
- if(d > 0.0e0) goto S30;
- gam1 = *a*w;
- return gam1;
-S30:
- gam1 = t/ *a*(w-0.5e0-0.5e0);
- return gam1;
-S40:
- top = (((((((r[8]*t+r[7])*t+r[6])*t+r[5])*t+r[4])*t+r[3])*t+r[2])*t+r[1])*t+
- r[0];
- bot = (s2*t+s1)*t+1.0e0;
- w = top/bot;
- if(d > 0.0e0) goto S50;
- gam1 = *a*(w+0.5e0+0.5e0);
- return gam1;
-S50:
- gam1 = t*w/ *a;
- return gam1;
-}
-void gaminv(double *a,double *x,double *x0,double *p,double *q,
- int *ierr)
-/*
- ----------------------------------------------------------------------
- INVERSE INCOMPLETE GAMMA RATIO FUNCTION
-
- GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1.
- THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER
- ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X
- TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE
- PARTICULAR COMPUTER ARITHMETIC BEING USED.
-
- ------------
-
- X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0,
- AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT
- NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN
- A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE
- IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X.
-
- X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER
- DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET
- X0 .LE. 0.
-
- IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
- WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING
- VALUES ...
-
- IERR = 0 THE SOLUTION WAS OBTAINED. ITERATION WAS
- NOT USED.
- IERR.GT.0 THE SOLUTION WAS OBTAINED. IERR ITERATIONS
- WERE PERFORMED.
- IERR = -2 (INPUT ERROR) A .LE. 0
- IERR = -3 NO SOLUTION WAS OBTAINED. THE RATIO Q/A
- IS TOO LARGE.
- IERR = -4 (INPUT ERROR) P + Q .NE. 1
- IERR = -6 20 ITERATIONS WERE PERFORMED. THE MOST
- RECENT VALUE OBTAINED FOR X IS GIVEN.
- THIS CANNOT OCCUR IF X0 .LE. 0.
- IERR = -7 ITERATION FAILED. NO VALUE IS GIVEN FOR X.
- THIS MAY OCCUR WHEN X IS APPROXIMATELY 0.
- IERR = -8 A VALUE FOR X HAS BEEN OBTAINED, BUT THE
- ROUTINE IS NOT CERTAIN OF ITS ACCURACY.
- ITERATION CANNOT BE PERFORMED IN THIS
- CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY
- WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS
- POSITIVE THEN THIS CAN OCCUR WHEN A IS
- EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY
- LARGE (SAY A .GE. 1.E20).
- ----------------------------------------------------------------------
- WRITTEN BY ALFRED H. MORRIS, JR.
- NAVAL SURFACE WEAPONS CENTER
- DAHLGREN, VIRGINIA
- -------------------
-*/
-{
-static double a0 = 3.31125922108741e0;
-static double a1 = 11.6616720288968e0;
-static double a2 = 4.28342155967104e0;
-static double a3 = .213623493715853e0;
-static double b1 = 6.61053765625462e0;
-static double b2 = 6.40691597760039e0;
-static double b3 = 1.27364489782223e0;
-static double b4 = .036117081018842e0;
-static double c = .577215664901533e0;
-static double ln10 = 2.302585e0;
-static double tol = 1.e-5;
-static double amin[2] = {
- 500.0e0,100.0e0
-};
-static double bmin[2] = {
- 1.e-28,1.e-13
-};
-static double dmin[2] = {
- 1.e-06,1.e-04
-};
-static double emin[2] = {
- 2.e-03,6.e-03
-};
-static double eps0[2] = {
- 1.e-10,1.e-08
-};
-static int K1 = 1;
-static int K2 = 2;
-static int K3 = 3;
-static int K8 = 0;
-static double am1,amax,ap1,ap2,ap3,apn,b,c1,c2,c3,c4,c5,d,e,e2,eps,g,h,pn,qg,qn,
- r,rta,s,s2,sum,t,u,w,xmax,xmin,xn,y,z;
-static int iop;
-static double T4,T5,T6,T7,T9;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
- E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
- XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
- LARGEST POSITIVE NUMBER.
-*/
- e = spmpar(&K1);
- xmin = spmpar(&K2);
- xmax = spmpar(&K3);
- *x = 0.0e0;
- if(*a <= 0.0e0) goto S300;
- t = *p+*q-1.e0;
- if(fabs(t) > e) goto S320;
- *ierr = 0;
- if(*p == 0.0e0) return;
- if(*q == 0.0e0) goto S270;
- if(*a == 1.0e0) goto S280;
- e2 = 2.0e0*e;
- amax = 0.4e-10/(e*e);
- iop = 1;
- if(e > 1.e-10) iop = 2;
- eps = eps0[iop-1];
- xn = *x0;
- if(*x0 > 0.0e0) goto S160;
-/*
- SELECTION OF THE INITIAL APPROXIMATION XN OF X
- WHEN A .LT. 1
-*/
- if(*a > 1.0e0) goto S80;
- T4 = *a+1.0e0;
- g = Xgamm(&T4);
- qg = *q*g;
- if(qg == 0.0e0) goto S360;
- b = qg/ *a;
- if(qg > 0.6e0**a) goto S40;
- if(*a >= 0.30e0 || b < 0.35e0) goto S10;
- t = exp(-(b+c));
- u = t*exp(t);
- xn = t*exp(u);
- goto S160;
-S10:
- if(b >= 0.45e0) goto S40;
- if(b == 0.0e0) goto S360;
- y = -log(b);
- s = 0.5e0+(0.5e0-*a);
- z = log(y);
- t = y-s*z;
- if(b < 0.15e0) goto S20;
- xn = y-s*log(t)-log(1.0e0+s/(t+1.0e0));
- goto S220;
-S20:
- if(b <= 0.01e0) goto S30;
- u = ((t+2.0e0*(3.0e0-*a))*t+(2.0e0-*a)*(3.0e0-*a))/((t+(5.0e0-*a))*t+2.0e0);
- xn = y-s*log(t)-log(u);
- goto S220;
-S30:
- c1 = -(s*z);
- c2 = -(s*(1.0e0+c1));
- c3 = s*((0.5e0*c1+(2.0e0-*a))*c1+(2.5e0-1.5e0**a));
- c4 = -(s*(((c1/3.0e0+(2.5e0-1.5e0**a))*c1+((*a-6.0e0)**a+7.0e0))*c1+(
- (11.0e0**a-46.0)**a+47.0e0)/6.0e0));
- c5 = -(s*((((-(c1/4.0e0)+(11.0e0**a-17.0e0)/6.0e0)*c1+((-(3.0e0**a)+13.0e0)*
- *a-13.0e0))*c1+0.5e0*(((2.0e0**a-25.0e0)**a+72.0e0)**a-61.0e0))*c1+((
- (25.0e0**a-195.0e0)**a+477.0e0)**a-379.0e0)/12.0e0));
- xn = (((c5/y+c4)/y+c3)/y+c2)/y+c1+y;
- if(*a > 1.0e0) goto S220;
- if(b > bmin[iop-1]) goto S220;
- *x = xn;
- return;
-S40:
- if(b**q > 1.e-8) goto S50;
- xn = exp(-(*q/ *a+c));
- goto S70;
-S50:
- if(*p <= 0.9e0) goto S60;
- T5 = -*q;
- xn = exp((alnrel(&T5)+gamln1(a))/ *a);
- goto S70;
-S60:
- xn = exp(log(*p*g)/ *a);
-S70:
- if(xn == 0.0e0) goto S310;
- t = 0.5e0+(0.5e0-xn/(*a+1.0e0));
- xn /= t;
- goto S160;
-S80:
-/*
- SELECTION OF THE INITIAL APPROXIMATION XN OF X
- WHEN A .GT. 1
-*/
- if(*q <= 0.5e0) goto S90;
- w = log(*p);
- goto S100;
-S90:
- w = log(*q);
-S100:
- t = sqrt(-(2.0e0*w));
- s = t-(((a3*t+a2)*t+a1)*t+a0)/((((b4*t+b3)*t+b2)*t+b1)*t+1.0e0);
- if(*q > 0.5e0) s = -s;
- rta = sqrt(*a);
- s2 = s*s;
- xn = *a+s*rta+(s2-1.0e0)/3.0e0+s*(s2-7.0e0)/(36.0e0*rta)-((3.0e0*s2+7.0e0)*
- s2-16.0e0)/(810.0e0**a)+s*((9.0e0*s2+256.0e0)*s2-433.0e0)/(38880.0e0**a*
- rta);
- xn = fifdmax1(xn,0.0e0);
- if(*a < amin[iop-1]) goto S110;
- *x = xn;
- d = 0.5e0+(0.5e0-*x/ *a);
- if(fabs(d) <= dmin[iop-1]) return;
-S110:
- if(*p <= 0.5e0) goto S130;
- if(xn < 3.0e0**a) goto S220;
- y = -(w+gamln(a));
- d = fifdmax1(2.0e0,*a*(*a-1.0e0));
- if(y < ln10*d) goto S120;
- s = 1.0e0-*a;
- z = log(y);
- goto S30;
-S120:
- t = *a-1.0e0;
- T6 = -(t/(xn+1.0e0));
- xn = y+t*log(xn)-alnrel(&T6);
- T7 = -(t/(xn+1.0e0));
- xn = y+t*log(xn)-alnrel(&T7);
- goto S220;
-S130:
- ap1 = *a+1.0e0;
- if(xn > 0.70e0*ap1) goto S170;
- w += gamln(&ap1);
- if(xn > 0.15e0*ap1) goto S140;
- ap2 = *a+2.0e0;
- ap3 = *a+3.0e0;
- *x = exp((w+*x)/ *a);
- *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
- *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
- *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2*(1.0e0+*x/ap3))))/ *a);
- xn = *x;
- if(xn > 1.e-2*ap1) goto S140;
- if(xn <= emin[iop-1]*ap1) return;
- goto S170;
-S140:
- apn = ap1;
- t = xn/apn;
- sum = 1.0e0+t;
-S150:
- apn += 1.0e0;
- t *= (xn/apn);
- sum += t;
- if(t > 1.e-4) goto S150;
- t = w-log(sum);
- xn = exp((xn+t)/ *a);
- xn *= (1.0e0-(*a*log(xn)-xn-t)/(*a-xn));
- goto S170;
-S160:
-/*
- SCHRODER ITERATION USING P
-*/
- if(*p > 0.5e0) goto S220;
-S170:
- if(*p <= 1.e10*xmin) goto S350;
- am1 = *a-0.5e0-0.5e0;
-S180:
- if(*a <= amax) goto S190;
- d = 0.5e0+(0.5e0-xn/ *a);
- if(fabs(d) <= e2) goto S350;
-S190:
- if(*ierr >= 20) goto S330;
- *ierr += 1;
- gratio(a,&xn,&pn,&qn,&K8);
- if(pn == 0.0e0 || qn == 0.0e0) goto S350;
- r = rcomp(a,&xn);
- if(r == 0.0e0) goto S350;
- t = (pn-*p)/r;
- w = 0.5e0*(am1-xn);
- if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S200;
- *x = xn*(1.0e0-t);
- if(*x <= 0.0e0) goto S340;
- d = fabs(t);
- goto S210;
-S200:
- h = t*(1.0e0+w*t);
- *x = xn*(1.0e0-h);
- if(*x <= 0.0e0) goto S340;
- if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
- d = fabs(h);
-S210:
- xn = *x;
- if(d > tol) goto S180;
- if(d <= eps) return;
- if(fabs(*p-pn) <= tol**p) return;
- goto S180;
-S220:
-/*
- SCHRODER ITERATION USING Q
-*/
- if(*q <= 1.e10*xmin) goto S350;
- am1 = *a-0.5e0-0.5e0;
-S230:
- if(*a <= amax) goto S240;
- d = 0.5e0+(0.5e0-xn/ *a);
- if(fabs(d) <= e2) goto S350;
-S240:
- if(*ierr >= 20) goto S330;
- *ierr += 1;
- gratio(a,&xn,&pn,&qn,&K8);
- if(pn == 0.0e0 || qn == 0.0e0) goto S350;
- r = rcomp(a,&xn);
- if(r == 0.0e0) goto S350;
- t = (*q-qn)/r;
- w = 0.5e0*(am1-xn);
- if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S250;
- *x = xn*(1.0e0-t);
- if(*x <= 0.0e0) goto S340;
- d = fabs(t);
- goto S260;
-S250:
- h = t*(1.0e0+w*t);
- *x = xn*(1.0e0-h);
- if(*x <= 0.0e0) goto S340;
- if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
- d = fabs(h);
-S260:
- xn = *x;
- if(d > tol) goto S230;
- if(d <= eps) return;
- if(fabs(*q-qn) <= tol**q) return;
- goto S230;
-S270:
-/*
- SPECIAL CASES
-*/
- *x = xmax;
- return;
-S280:
- if(*q < 0.9e0) goto S290;
- T9 = -*p;
- *x = -alnrel(&T9);
- return;
-S290:
- *x = -log(*q);
- return;
-S300:
-/*
- ERROR RETURN
-*/
- *ierr = -2;
- return;
-S310:
- *ierr = -3;
- return;
-S320:
- *ierr = -4;
- return;
-S330:
- *ierr = -6;
- return;
-S340:
- *ierr = -7;
- return;
-S350:
- *x = xn;
- *ierr = -8;
- return;
-S360:
- *x = xmax;
- *ierr = -8;
- return;
-}
-double gamln(double *a)
-/*
------------------------------------------------------------------------
- EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
------------------------------------------------------------------------
- WRITTEN BY ALFRED H. MORRIS
- NAVAL SURFACE WARFARE CENTER
- DAHLGREN, VIRGINIA
---------------------------
- D = 0.5*(LN(2*PI) - 1)
---------------------------
-*/
-{
-static double c0 = .833333333333333e-01;
-static double c1 = -.277777777760991e-02;
-static double c2 = .793650666825390e-03;
-static double c3 = -.595202931351870e-03;
-static double c4 = .837308034031215e-03;
-static double c5 = -.165322962780713e-02;
-static double d = .418938533204673e0;
-static double gamln,t,w;
-static int i,n;
-static double T1;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*a > 0.8e0) goto S10;
- gamln = gamln1(a)-log(*a);
- return gamln;
-S10:
- if(*a > 2.25e0) goto S20;
- t = *a-0.5e0-0.5e0;
- gamln = gamln1(&t);
- return gamln;
-S20:
- if(*a >= 10.0e0) goto S40;
- n = (long)(*a - 1.25e0);
- t = *a;
- w = 1.0e0;
- for(i=1; i<=n; i++) {
- t -= 1.0e0;
- w = t*w;
- }
- T1 = t-1.0e0;
- gamln = gamln1(&T1)+log(w);
- return gamln;
-S40:
- t = pow(1.0e0/ *a,2.0);
- w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a;
- gamln = d+w+(*a-0.5e0)*(log(*a)-1.0e0);
- return gamln;
-}
-double gamln1(double *a)
-/*
------------------------------------------------------------------------
- EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
------------------------------------------------------------------------
-*/
-{
-static double p0 = .577215664901533e+00;
-static double p1 = .844203922187225e+00;
-static double p2 = -.168860593646662e+00;
-static double p3 = -.780427615533591e+00;
-static double p4 = -.402055799310489e+00;
-static double p5 = -.673562214325671e-01;
-static double p6 = -.271935708322958e-02;
-static double q1 = .288743195473681e+01;
-static double q2 = .312755088914843e+01;
-static double q3 = .156875193295039e+01;
-static double q4 = .361951990101499e+00;
-static double q5 = .325038868253937e-01;
-static double q6 = .667465618796164e-03;
-static double r0 = .422784335098467e+00;
-static double r1 = .848044614534529e+00;
-static double r2 = .565221050691933e+00;
-static double r3 = .156513060486551e+00;
-static double r4 = .170502484022650e-01;
-static double r5 = .497958207639485e-03;
-static double s1 = .124313399877507e+01;
-static double s2 = .548042109832463e+00;
-static double s3 = .101552187439830e+00;
-static double s4 = .713309612391000e-02;
-static double s5 = .116165475989616e-03;
-static double gamln1,w,x;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*a >= 0.6e0) goto S10;
- w = ((((((p6**a+p5)**a+p4)**a+p3)**a+p2)**a+p1)**a+p0)/((((((q6**a+q5)**a+
- q4)**a+q3)**a+q2)**a+q1)**a+1.0e0);
- gamln1 = -(*a*w);
- return gamln1;
-S10:
- x = *a-0.5e0-0.5e0;
- w = (((((r5*x+r4)*x+r3)*x+r2)*x+r1)*x+r0)/(((((s5*x+s4)*x+s3)*x+s2)*x+s1)*x
- +1.0e0);
- gamln1 = x*w;
- return gamln1;
-}
-double Xgamm(double *a)
-/*
------------------------------------------------------------------------
-
- EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS
-
- -----------
-
- GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
- BE COMPUTED.
-
------------------------------------------------------------------------
- WRITTEN BY ALFRED H. MORRIS, JR.
- NAVAL SURFACE WEAPONS CENTER
- DAHLGREN, VIRGINIA
------------------------------------------------------------------------
-*/
-{
-static double d = .41893853320467274178e0;
-static double pi = 3.1415926535898e0;
-static double r1 = .820756370353826e-03;
-static double r2 = -.595156336428591e-03;
-static double r3 = .793650663183693e-03;
-static double r4 = -.277777777770481e-02;
-static double r5 = .833333333333333e-01;
-static double p[7] = {
- .539637273585445e-03,.261939260042690e-02,.204493667594920e-01,
- .730981088720487e-01,.279648642639792e+00,.553413866010467e+00,1.0e0
-};
-static double q[7] = {
- -.832979206704073e-03,.470059485860584e-02,.225211131035340e-01,
- -.170458969313360e+00,-.567902761974940e-01,.113062953091122e+01,1.0e0
-};
-static int K2 = 3;
-static int K3 = 0;
-static double Xgamm,bot,g,lnx,s,t,top,w,x,z;
-static int i,j,m,n,T1;
-/*
- ..
- .. Executable Statements ..
-*/
- Xgamm = 0.0e0;
- x = *a;
- if(fabs(*a) >= 15.0e0) goto S110;
-/*
------------------------------------------------------------------------
- EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15
------------------------------------------------------------------------
-*/
- t = 1.0e0;
- m = fifidint(*a)-1;
-/*
- LET T BE THE PRODUCT OF A-J WHEN A .GE. 2
-*/
- T1 = m;
- if(T1 < 0) goto S40;
- else if(T1 == 0) goto S30;
- else goto S10;
-S10:
- for(j=1; j<=m; j++) {
- x -= 1.0e0;
- t = x*t;
- }
-S30:
- x -= 1.0e0;
- goto S80;
-S40:
-/*
- LET T BE THE PRODUCT OF A+J WHEN A .LT. 1
-*/
- t = *a;
- if(*a > 0.0e0) goto S70;
- m = -m-1;
- if(m == 0) goto S60;
- for(j=1; j<=m; j++) {
- x += 1.0e0;
- t = x*t;
- }
-S60:
- x += (0.5e0+0.5e0);
- t = x*t;
- if(t == 0.0e0) return Xgamm;
-S70:
-/*
- THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
- CODE MAY BE OMITTED IF DESIRED.
-*/
- if(fabs(t) >= 1.e-30) goto S80;
- if(fabs(t)*spmpar(&K2) <= 1.0001e0) return Xgamm;
- Xgamm = 1.0e0/t;
- return Xgamm;
-S80:
-/*
- COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1
-*/
- top = p[0];
- bot = q[0];
- for(i=1; i<7; i++) {
- top = p[i]+x*top;
- bot = q[i]+x*bot;
- }
- Xgamm = top/bot;
-/*
- TERMINATION
-*/
- if(*a < 1.0e0) goto S100;
- Xgamm *= t;
- return Xgamm;
-S100:
- Xgamm /= t;
- return Xgamm;
-S110:
-/*
------------------------------------------------------------------------
- EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15
------------------------------------------------------------------------
-*/
- if(fabs(*a) >= 1.e3) return Xgamm;
- if(*a > 0.0e0) goto S120;
- x = -*a;
- n = (long)(x);
- t = x-(double)n;
- if(t > 0.9e0) t = 1.0e0-t;
- s = sin(pi*t)/pi;
- if(fifmod(n,2) == 0) s = -s;
- if(s == 0.0e0) return Xgamm;
-S120:
-/*
- COMPUTE THE MODIFIED ASYMPTOTIC SUM
-*/
- t = 1.0e0/(x*x);
- g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x;
-/*
- ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X)
- BUT LESS ACCURACY WILL NORMALLY BE OBTAINED.
-*/
- lnx = log(x);
-/*
- FINAL ASSEMBLY
-*/
- z = x;
- g = d+g+(z-0.5e0)*(lnx-1.e0);
- w = g;
- t = g-w;
- if(w > 0.99999e0*exparg(&K3)) return Xgamm;
- Xgamm = exp(w)*(1.0e0+t);
- if(*a < 0.0e0) Xgamm = 1.0e0/(Xgamm*s)/x;
- return Xgamm;
-}
-void grat1(double *a,double *x,double *r,double *p,double *q,
- double *eps)
-{
-static int K2 = 0;
-static double a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,t,tol,w,z,T1,T3;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
------------------------------------------------------------------------
- EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
- P(A,X) AND Q(A,X)
- IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED.
- THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A).
------------------------------------------------------------------------
-*/
- if(*a**x == 0.0e0) goto S120;
- if(*a == 0.5e0) goto S100;
- if(*x < 1.1e0) goto S10;
- goto S60;
-S10:
-/*
- TAYLOR SERIES FOR P(A,X)/X**A
-*/
- an = 3.0e0;
- c = *x;
- sum = *x/(*a+3.0e0);
- tol = 0.1e0**eps/(*a+1.0e0);
-S20:
- an += 1.0e0;
- c = -(c*(*x/an));
- t = c/(*a+an);
- sum += t;
- if(fabs(t) > tol) goto S20;
- j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
- z = *a*log(*x);
- h = gam1(a);
- g = 1.0e0+h;
- if(*x < 0.25e0) goto S30;
- if(*a < *x/2.59e0) goto S50;
- goto S40;
-S30:
- if(z > -.13394e0) goto S50;
-S40:
- w = exp(z);
- *p = w*g*(0.5e0+(0.5e0-j));
- *q = 0.5e0+(0.5e0-*p);
- return;
-S50:
- l = rexp(&z);
- w = 0.5e0+(0.5e0+l);
- *q = (w*j-l)*g-h;
- if(*q < 0.0e0) goto S90;
- *p = 0.5e0+(0.5e0-*q);
- return;
-S60:
-/*
- CONTINUED FRACTION EXPANSION
-*/
- a2nm1 = a2n = 1.0e0;
- b2nm1 = *x;
- b2n = *x+(1.0e0-*a);
- c = 1.0e0;
-S70:
- a2nm1 = *x*a2n+c*a2nm1;
- b2nm1 = *x*b2n+c*b2nm1;
- am0 = a2nm1/b2nm1;
- c += 1.0e0;
- cma = c-*a;
- a2n = a2nm1+cma*a2n;
- b2n = b2nm1+cma*b2n;
- an0 = a2n/b2n;
- if(fabs(an0-am0) >= *eps*an0) goto S70;
- *q = *r*an0;
- *p = 0.5e0+(0.5e0-*q);
- return;
-S80:
-/*
- SPECIAL CASES
-*/
- *p = 0.0e0;
- *q = 1.0e0;
- return;
-S90:
- *p = 1.0e0;
- *q = 0.0e0;
- return;
-S100:
- if(*x >= 0.25e0) goto S110;
- T1 = sqrt(*x);
- *p = erf1(&T1);
- *q = 0.5e0+(0.5e0-*p);
- return;
-S110:
- T3 = sqrt(*x);
- *q = erfc1(&K2,&T3);
- *p = 0.5e0+(0.5e0-*q);
- return;
-S120:
- if(*x <= *a) goto S80;
- goto S90;
-}
-void gratio(double *a,double *x,double *ans,double *qans,int *ind)
-/*
- ----------------------------------------------------------------------
- EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
- P(A,X) AND Q(A,X)
-
- ----------
-
- IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
- ARE NOT BOTH 0.
-
- ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
- P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
- IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
- POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
- IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
- 6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
- IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.
-
- ERROR RETURN ...
- ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
- WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
- P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
- X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.
- ----------------------------------------------------------------------
- WRITTEN BY ALFRED H. MORRIS, JR.
- NAVAL SURFACE WEAPONS CENTER
- DAHLGREN, VIRGINIA
- --------------------
-*/
-{
-static double alog10 = 2.30258509299405e0;
-static double d10 = -.185185185185185e-02;
-static double d20 = .413359788359788e-02;
-static double d30 = .649434156378601e-03;
-static double d40 = -.861888290916712e-03;
-static double d50 = -.336798553366358e-03;
-static double d60 = .531307936463992e-03;
-static double d70 = .344367606892378e-03;
-static double rt2pin = .398942280401433e0;
-static double rtpi = 1.77245385090552e0;
-static double third = .333333333333333e0;
-static double acc0[3] = {
- 5.e-15,5.e-7,5.e-4
-};
-static double big[3] = {
- 20.0e0,14.0e0,10.0e0
-};
-static double d0[13] = {
- .833333333333333e-01,-.148148148148148e-01,.115740740740741e-02,
- .352733686067019e-03,-.178755144032922e-03,.391926317852244e-04,
- -.218544851067999e-05,-.185406221071516e-05,.829671134095309e-06,
- -.176659527368261e-06,.670785354340150e-08,.102618097842403e-07,
- -.438203601845335e-08
-};
-static double d1[12] = {
- -.347222222222222e-02,.264550264550265e-02,-.990226337448560e-03,
- .205761316872428e-03,-.401877572016461e-06,-.180985503344900e-04,
- .764916091608111e-05,-.161209008945634e-05,.464712780280743e-08,
- .137863344691572e-06,-.575254560351770e-07,.119516285997781e-07
-};
-static double d2[10] = {
- -.268132716049383e-02,.771604938271605e-03,.200938786008230e-05,
- -.107366532263652e-03,.529234488291201e-04,-.127606351886187e-04,
- .342357873409614e-07,.137219573090629e-05,-.629899213838006e-06,
- .142806142060642e-06
-};
-static double d3[8] = {
- .229472093621399e-03,-.469189494395256e-03,.267720632062839e-03,
- -.756180167188398e-04,-.239650511386730e-06,.110826541153473e-04,
- -.567495282699160e-05,.142309007324359e-05
-};
-static double d4[6] = {
- .784039221720067e-03,-.299072480303190e-03,-.146384525788434e-05,
- .664149821546512e-04,-.396836504717943e-04,.113757269706784e-04
-};
-static double d5[4] = {
- -.697281375836586e-04,.277275324495939e-03,-.199325705161888e-03,
- .679778047793721e-04
-};
-static double d6[2] = {
- -.592166437353694e-03,.270878209671804e-03
-};
-static double e00[3] = {
- .25e-3,.25e-1,.14e0
-};
-static double x00[3] = {
- 31.0e0,17.0e0,9.7e0
-};
-static int K1 = 1;
-static int K2 = 0;
-static double a2n,a2nm1,acc,am0,amn,an,an0,apn,b2n,b2nm1,c,c0,c1,c2,c3,c4,c5,c6,
- cma,e,e0,g,h,j,l,r,rta,rtx,s,sum,t,t1,tol,twoa,u,w,x0,y,z;
-static int i,iop,m,max,n;
-static double wk[20],T3;
-static int T4,T5;
-static double T6,T7;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
- --------------------
- ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
- FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
-*/
- e = spmpar(&K1);
- if(*a < 0.0e0 || *x < 0.0e0) goto S430;
- if(*a == 0.0e0 && *x == 0.0e0) goto S430;
- if(*a**x == 0.0e0) goto S420;
- iop = *ind+1;
- if(iop != 1 && iop != 2) iop = 3;
- acc = fifdmax1(acc0[iop-1],e);
- e0 = e00[iop-1];
- x0 = x00[iop-1];
-/*
- SELECT THE APPROPRIATE ALGORITHM
-*/
- if(*a >= 1.0e0) goto S10;
- if(*a == 0.5e0) goto S390;
- if(*x < 1.1e0) goto S160;
- t1 = *a*log(*x)-*x;
- u = *a*exp(t1);
- if(u == 0.0e0) goto S380;
- r = u*(1.0e0+gam1(a));
- goto S250;
-S10:
- if(*a >= big[iop-1]) goto S30;
- if(*a > *x || *x >= x0) goto S20;
- twoa = *a+*a;
- m = fifidint(twoa);
- if(twoa != (double)m) goto S20;
- i = m/2;
- if(*a == (double)i) goto S210;
- goto S220;
-S20:
- t1 = *a*log(*x)-*x;
- r = exp(t1)/Xgamm(a);
- goto S40;
-S30:
- l = *x/ *a;
- if(l == 0.0e0) goto S370;
- s = 0.5e0+(0.5e0-l);
- z = rlog(&l);
- if(z >= 700.0e0/ *a) goto S410;
- y = *a*z;
- rta = sqrt(*a);
- if(fabs(s) <= e0/rta) goto S330;
- if(fabs(s) <= 0.4e0) goto S270;
- t = pow(1.0e0/ *a,2.0);
- t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
- t1 -= y;
- r = rt2pin*rta*exp(t1);
-S40:
- if(r == 0.0e0) goto S420;
- if(*x <= fifdmax1(*a,alog10)) goto S50;
- if(*x < x0) goto S250;
- goto S100;
-S50:
-/*
- TAYLOR SERIES FOR P/R
-*/
- apn = *a+1.0e0;
- t = *x/apn;
- wk[0] = t;
- for(n=2; n<=20; n++) {
- apn += 1.0e0;
- t *= (*x/apn);
- if(t <= 1.e-3) goto S70;
- wk[n-1] = t;
- }
- n = 20;
-S70:
- sum = t;
- tol = 0.5e0*acc;
-S80:
- apn += 1.0e0;
- t *= (*x/apn);
- sum += t;
- if(t > tol) goto S80;
- max = n-1;
- for(m=1; m<=max; m++) {
- n -= 1;
- sum += wk[n-1];
- }
- *ans = r/ *a*(1.0e0+sum);
- *qans = 0.5e0+(0.5e0-*ans);
- return;
-S100:
-/*
- ASYMPTOTIC EXPANSION
-*/
- amn = *a-1.0e0;
- t = amn/ *x;
- wk[0] = t;
- for(n=2; n<=20; n++) {
- amn -= 1.0e0;
- t *= (amn/ *x);
- if(fabs(t) <= 1.e-3) goto S120;
- wk[n-1] = t;
- }
- n = 20;
-S120:
- sum = t;
-S130:
- if(fabs(t) <= acc) goto S140;
- amn -= 1.0e0;
- t *= (amn/ *x);
- sum += t;
- goto S130;
-S140:
- max = n-1;
- for(m=1; m<=max; m++) {
- n -= 1;
- sum += wk[n-1];
- }
- *qans = r/ *x*(1.0e0+sum);
- *ans = 0.5e0+(0.5e0-*qans);
- return;
-S160:
-/*
- TAYLOR SERIES FOR P(A,X)/X**A
-*/
- an = 3.0e0;
- c = *x;
- sum = *x/(*a+3.0e0);
- tol = 3.0e0*acc/(*a+1.0e0);
-S170:
- an += 1.0e0;
- c = -(c*(*x/an));
- t = c/(*a+an);
- sum += t;
- if(fabs(t) > tol) goto S170;
- j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
- z = *a*log(*x);
- h = gam1(a);
- g = 1.0e0+h;
- if(*x < 0.25e0) goto S180;
- if(*a < *x/2.59e0) goto S200;
- goto S190;
-S180:
- if(z > -.13394e0) goto S200;
-S190:
- w = exp(z);
- *ans = w*g*(0.5e0+(0.5e0-j));
- *qans = 0.5e0+(0.5e0-*ans);
- return;
-S200:
- l = rexp(&z);
- w = 0.5e0+(0.5e0+l);
- *qans = (w*j-l)*g-h;
- if(*qans < 0.0e0) goto S380;
- *ans = 0.5e0+(0.5e0-*qans);
- return;
-S210:
-/*
- FINITE SUMS FOR Q WHEN A .GE. 1
- AND 2*A IS AN INTEGER
-*/
- sum = exp(-*x);
- t = sum;
- n = 1;
- c = 0.0e0;
- goto S230;
-S220:
- rtx = sqrt(*x);
- sum = erfc1(&K2,&rtx);
- t = exp(-*x)/(rtpi*rtx);
- n = 0;
- c = -0.5e0;
-S230:
- if(n == i) goto S240;
- n += 1;
- c += 1.0e0;
- t = *x*t/c;
- sum += t;
- goto S230;
-S240:
- *qans = sum;
- *ans = 0.5e0+(0.5e0-*qans);
- return;
-S250:
-/*
- CONTINUED FRACTION EXPANSION
-*/
- tol = fifdmax1(5.0e0*e,acc);
- a2nm1 = a2n = 1.0e0;
- b2nm1 = *x;
- b2n = *x+(1.0e0-*a);
- c = 1.0e0;
-S260:
- a2nm1 = *x*a2n+c*a2nm1;
- b2nm1 = *x*b2n+c*b2nm1;
- am0 = a2nm1/b2nm1;
- c += 1.0e0;
- cma = c-*a;
- a2n = a2nm1+cma*a2n;
- b2n = b2nm1+cma*b2n;
- an0 = a2n/b2n;
- if(fabs(an0-am0) >= tol*an0) goto S260;
- *qans = r*an0;
- *ans = 0.5e0+(0.5e0-*qans);
- return;
-S270:
-/*
- GENERAL TEMME EXPANSION
-*/
- if(fabs(s) <= 2.0e0*e && *a*e*e > 3.28e-3) goto S430;
- c = exp(-y);
- T3 = sqrt(y);
- w = 0.5e0*erfc1(&K1,&T3);
- u = 1.0e0/ *a;
- z = sqrt(z+z);
- if(l < 1.0e0) z = -z;
- T4 = iop-2;
- if(T4 < 0) goto S280;
- else if(T4 == 0) goto S290;
- else goto S300;
-S280:
- if(fabs(s) <= 1.e-3) goto S340;
- c0 = ((((((((((((d0[12]*z+d0[11])*z+d0[10])*z+d0[9])*z+d0[8])*z+d0[7])*z+d0[
- 6])*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
- c1 = (((((((((((d1[11]*z+d1[10])*z+d1[9])*z+d1[8])*z+d1[7])*z+d1[6])*z+d1[5]
- )*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
- c2 = (((((((((d2[9]*z+d2[8])*z+d2[7])*z+d2[6])*z+d2[5])*z+d2[4])*z+d2[3])*z+
- d2[2])*z+d2[1])*z+d2[0])*z+d20;
- c3 = (((((((d3[7]*z+d3[6])*z+d3[5])*z+d3[4])*z+d3[3])*z+d3[2])*z+d3[1])*z+
- d3[0])*z+d30;
- c4 = (((((d4[5]*z+d4[4])*z+d4[3])*z+d4[2])*z+d4[1])*z+d4[0])*z+d40;
- c5 = (((d5[3]*z+d5[2])*z+d5[1])*z+d5[0])*z+d50;
- c6 = (d6[1]*z+d6[0])*z+d60;
- t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
- goto S310;
-S290:
- c0 = (((((d0[5]*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
- c1 = (((d1[3]*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
- c2 = d2[0]*z+d20;
- t = (c2*u+c1)*u+c0;
- goto S310;
-S300:
- t = ((d0[2]*z+d0[1])*z+d0[0])*z-third;
-S310:
- if(l < 1.0e0) goto S320;
- *qans = c*(w+rt2pin*t/rta);
- *ans = 0.5e0+(0.5e0-*qans);
- return;
-S320:
- *ans = c*(w-rt2pin*t/rta);
- *qans = 0.5e0+(0.5e0-*ans);
- return;
-S330:
-/*
- TEMME EXPANSION FOR L = 1
-*/
- if(*a*e*e > 3.28e-3) goto S430;
- c = 0.5e0+(0.5e0-y);
- w = (0.5e0-sqrt(y)*(0.5e0+(0.5e0-y/3.0e0))/rtpi)/c;
- u = 1.0e0/ *a;
- z = sqrt(z+z);
- if(l < 1.0e0) z = -z;
- T5 = iop-2;
- if(T5 < 0) goto S340;
- else if(T5 == 0) goto S350;
- else goto S360;
-S340:
- c0 = ((((((d0[6]*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-
- third;
- c1 = (((((d1[5]*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
- c2 = ((((d2[4]*z+d2[3])*z+d2[2])*z+d2[1])*z+d2[0])*z+d20;
- c3 = (((d3[3]*z+d3[2])*z+d3[1])*z+d3[0])*z+d30;
- c4 = (d4[1]*z+d4[0])*z+d40;
- c5 = (d5[1]*z+d5[0])*z+d50;
- c6 = d6[0]*z+d60;
- t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
- goto S310;
-S350:
- c0 = (d0[1]*z+d0[0])*z-third;
- c1 = d1[0]*z+d10;
- t = (d20*u+c1)*u+c0;
- goto S310;
-S360:
- t = d0[0]*z-third;
- goto S310;
-S370:
-/*
- SPECIAL CASES
-*/
- *ans = 0.0e0;
- *qans = 1.0e0;
- return;
-S380:
- *ans = 1.0e0;
- *qans = 0.0e0;
- return;
-S390:
- if(*x >= 0.25e0) goto S400;
- T6 = sqrt(*x);
- *ans = erf1(&T6);
- *qans = 0.5e0+(0.5e0-*ans);
- return;
-S400:
- T7 = sqrt(*x);
- *qans = erfc1(&K2,&T7);
- *ans = 0.5e0+(0.5e0-*qans);
- return;
-S410:
- if(fabs(s) <= 2.0e0*e) goto S430;
-S420:
- if(*x <= *a) goto S370;
- goto S380;
-S430:
-/*
- ERROR RETURN
-*/
- *ans = 2.0e0;
- return;
-}
-double gsumln(double *a,double *b)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
- FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2
------------------------------------------------------------------------
-*/
-{
-static double gsumln,x,T1,T2;
-/*
- ..
- .. Executable Statements ..
-*/
- x = *a+*b-2.e0;
- if(x > 0.25e0) goto S10;
- T1 = 1.0e0+x;
- gsumln = gamln1(&T1);
- return gsumln;
-S10:
- if(x > 1.25e0) goto S20;
- gsumln = gamln1(&x)+alnrel(&x);
- return gsumln;
-S20:
- T2 = x-1.0e0;
- gsumln = gamln1(&T2)+log(x*(1.0e0+x));
- return gsumln;
-}
-double psi(double *xx)
-/*
----------------------------------------------------------------------
-
- EVALUATION OF THE DIGAMMA FUNCTION
-
- -----------
-
- PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
- BE COMPUTED.
-
- THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
- APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
- CODY, STRECOK AND THACHER.
-
----------------------------------------------------------------------
- PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
- PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
- A.H. MORRIS (NSWC).
----------------------------------------------------------------------
-*/
-{
-static double dx0 = 1.461632144968362341262659542325721325e0;
-static double piov4 = .785398163397448e0;
-static double p1[7] = {
- .895385022981970e-02,.477762828042627e+01,.142441585084029e+03,
- .118645200713425e+04,.363351846806499e+04,.413810161269013e+04,
- .130560269827897e+04
-};
-static double p2[4] = {
- -.212940445131011e+01,-.701677227766759e+01,-.448616543918019e+01,
- -.648157123766197e+00
-};
-static double q1[6] = {
- .448452573429826e+02,.520752771467162e+03,.221000799247830e+04,
- .364127349079381e+04,.190831076596300e+04,.691091682714533e-05
-};
-static double q2[4] = {
- .322703493791143e+02,.892920700481861e+02,.546117738103215e+02,
- .777788548522962e+01
-};
-static int K1 = 3;
-static int K2 = 1;
-static double psi,aug,den,sgn,upper,w,x,xmax1,xmx0,xsmall,z;
-static int i,m,n,nq;
-/*
- ..
- .. Executable Statements ..
-*/
-/*
----------------------------------------------------------------------
- MACHINE DEPENDENT CONSTANTS ...
- XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
- WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED
- AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
- ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
- PSI MAY BE REPRESENTED AS ALOG(X).
- XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
- MAY BE REPRESENTED BY 1/X.
----------------------------------------------------------------------
-*/
- xmax1 = ipmpar(&K1);
- xmax1 = fifdmin1(xmax1,1.0e0/spmpar(&K2));
- xsmall = 1.e-9;
- x = *xx;
- aug = 0.0e0;
- if(x >= 0.5e0) goto S50;
-/*
----------------------------------------------------------------------
- X .LT. 0.5, USE REFLECTION FORMULA
- PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
----------------------------------------------------------------------
-*/
- if(fabs(x) > xsmall) goto S10;
- if(x == 0.0e0) goto S100;
-/*
----------------------------------------------------------------------
- 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE
- FOR PI*COTAN(PI*X)
----------------------------------------------------------------------
-*/
- aug = -(1.0e0/x);
- goto S40;
-S10:
-/*
----------------------------------------------------------------------
- REDUCTION OF ARGUMENT FOR COTAN
----------------------------------------------------------------------
-*/
- w = -x;
- sgn = piov4;
- if(w > 0.0e0) goto S20;
- w = -w;
- sgn = -sgn;
-S20:
-/*
----------------------------------------------------------------------
- MAKE AN ERROR EXIT IF X .LE. -XMAX1
----------------------------------------------------------------------
-*/
- if(w >= xmax1) goto S100;
- nq = fifidint(w);
- w -= (double)nq;
- nq = fifidint(w*4.0e0);
- w = 4.0e0*(w-(double)nq*.25e0);
-/*
----------------------------------------------------------------------
- W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X.
- ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
- QUADRANT AND DETERMINE SIGN
----------------------------------------------------------------------
-*/
- n = nq/2;
- if(n+n != nq) w = 1.0e0-w;
- z = piov4*w;
- m = n/2;
- if(m+m != n) sgn = -sgn;
-/*
----------------------------------------------------------------------
- DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X)
----------------------------------------------------------------------
-*/
- n = (nq+1)/2;
- m = n/2;
- m += m;
- if(m != n) goto S30;
-/*
----------------------------------------------------------------------
- CHECK FOR SINGULARITY
----------------------------------------------------------------------
-*/
- if(z == 0.0e0) goto S100;
-/*
----------------------------------------------------------------------
- USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
- SIN/COS AS A SUBSTITUTE FOR TAN
----------------------------------------------------------------------
-*/
- aug = sgn*(cos(z)/sin(z)*4.0e0);
- goto S40;
-S30:
- aug = sgn*(sin(z)/cos(z)*4.0e0);
-S40:
- x = 1.0e0-x;
-S50:
- if(x > 3.0e0) goto S70;
-/*
----------------------------------------------------------------------
- 0.5 .LE. X .LE. 3.0
----------------------------------------------------------------------
-*/
- den = x;
- upper = p1[0]*x;
- for(i=1; i<=5; i++) {
- den = (den+q1[i-1])*x;
- upper = (upper+p1[i+1-1])*x;
- }
- den = (upper+p1[6])/(den+q1[5]);
- xmx0 = x-dx0;
- psi = den*xmx0+aug;
- return psi;
-S70:
-/*
----------------------------------------------------------------------
- IF X .GE. XMAX1, PSI = LN(X)
----------------------------------------------------------------------
-*/
- if(x >= xmax1) goto S90;
-/*
----------------------------------------------------------------------
- 3.0 .LT. X .LT. XMAX1
----------------------------------------------------------------------
-*/
- w = 1.0e0/(x*x);
- den = w;
- upper = p2[0]*w;
- for(i=1; i<=3; i++) {
- den = (den+q2[i-1])*w;
- upper = (upper+p2[i+1-1])*w;
- }
- aug = upper/(den+q2[3])-0.5e0/x+aug;
-S90:
- psi = aug+log(x);
- return psi;
-S100:
-/*
----------------------------------------------------------------------
- ERROR RETURN
----------------------------------------------------------------------
-*/
- psi = 0.0e0;
- return psi;
-}
-double rcomp(double *a,double *x)
-/*
- -------------------
- EVALUATION OF EXP(-X)*X**A/GAMMA(A)
- -------------------
- RT2PIN = 1/SQRT(2*PI)
- -------------------
-*/
-{
-static double rt2pin = .398942280401433e0;
-static double rcomp,t,t1,u;
-/*
- ..
- .. Executable Statements ..
-*/
- rcomp = 0.0e0;
- if(*a >= 20.0e0) goto S20;
- t = *a*log(*x)-*x;
- if(*a >= 1.0e0) goto S10;
- rcomp = *a*exp(t)*(1.0e0+gam1(a));
- return rcomp;
-S10:
- rcomp = exp(t)/Xgamm(a);
- return rcomp;
-S20:
- u = *x/ *a;
- if(u == 0.0e0) return rcomp;
- t = pow(1.0e0/ *a,2.0);
- t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
- t1 -= (*a*rlog(&u));
- rcomp = rt2pin*sqrt(*a)*exp(t1);
- return rcomp;
-}
-double rexp(double *x)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE FUNCTION EXP(X) - 1
------------------------------------------------------------------------
-*/
-{
-static double p1 = .914041914819518e-09;
-static double p2 = .238082361044469e-01;
-static double q1 = -.499999999085958e+00;
-static double q2 = .107141568980644e+00;
-static double q3 = -.119041179760821e-01;
-static double q4 = .595130811860248e-03;
-static double rexp,w;
-/*
- ..
- .. Executable Statements ..
-*/
- if(fabs(*x) > 0.15e0) goto S10;
- rexp = *x*(((p2**x+p1)**x+1.0e0)/((((q4**x+q3)**x+q2)**x+q1)**x+1.0e0));
- return rexp;
-S10:
- w = exp(*x);
- if(*x > 0.0e0) goto S20;
- rexp = w-0.5e0-0.5e0;
- return rexp;
-S20:
- rexp = w*(0.5e0+(0.5e0-1.0e0/w));
- return rexp;
-}
-double rlog(double *x)
-/*
- -------------------
- COMPUTATION OF X - 1 - LN(X)
- -------------------
-*/
-{
-static double a = .566749439387324e-01;
-static double b = .456512608815524e-01;
-static double p0 = .333333333333333e+00;
-static double p1 = -.224696413112536e+00;
-static double p2 = .620886815375787e-02;
-static double q1 = -.127408923933623e+01;
-static double q2 = .354508718369557e+00;
-static double rlog,r,t,u,w,w1;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*x < 0.61e0 || *x > 1.57e0) goto S40;
- if(*x < 0.82e0) goto S10;
- if(*x > 1.18e0) goto S20;
-/*
- ARGUMENT REDUCTION
-*/
- u = *x-0.5e0-0.5e0;
- w1 = 0.0e0;
- goto S30;
-S10:
- u = *x-0.7e0;
- u /= 0.7e0;
- w1 = a-u*0.3e0;
- goto S30;
-S20:
- u = 0.75e0**x-1.e0;
- w1 = b+u/3.0e0;
-S30:
-/*
- SERIES EXPANSION
-*/
- r = u/(u+2.0e0);
- t = r*r;
- w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
- rlog = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
- return rlog;
-S40:
- r = *x-0.5e0-0.5e0;
- rlog = r-log(*x);
- return rlog;
-}
-double rlog1(double *x)
-/*
------------------------------------------------------------------------
- EVALUATION OF THE FUNCTION X - LN(1 + X)
------------------------------------------------------------------------
-*/
-{
-static double a = .566749439387324e-01;
-static double b = .456512608815524e-01;
-static double p0 = .333333333333333e+00;
-static double p1 = -.224696413112536e+00;
-static double p2 = .620886815375787e-02;
-static double q1 = -.127408923933623e+01;
-static double q2 = .354508718369557e+00;
-static double rlog1,h,r,t,w,w1;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*x < -0.39e0 || *x > 0.57e0) goto S40;
- if(*x < -0.18e0) goto S10;
- if(*x > 0.18e0) goto S20;
-/*
- ARGUMENT REDUCTION
-*/
- h = *x;
- w1 = 0.0e0;
- goto S30;
-S10:
- h = *x+0.3e0;
- h /= 0.7e0;
- w1 = a-h*0.3e0;
- goto S30;
-S20:
- h = 0.75e0**x-0.25e0;
- w1 = b+h/3.0e0;
-S30:
-/*
- SERIES EXPANSION
-*/
- r = h/(h+2.0e0);
- t = r*r;
- w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
- rlog1 = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
- return rlog1;
-S40:
- w = *x+0.5e0+0.5e0;
- rlog1 = *x-log(w);
- return rlog1;
-}
-double spmpar(int *i)
-/*
------------------------------------------------------------------------
-
- SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
- THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
- I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
- SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
- ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
-
- SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
-
- SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
-
- SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
-
------------------------------------------------------------------------
- WRITTEN BY
- ALFRED H. MORRIS, JR.
- NAVAL SURFACE WARFARE CENTER
- DAHLGREN VIRGINIA
------------------------------------------------------------------------
------------------------------------------------------------------------
- MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE
- CONSTANTS FOR THE COMPUTER BEING USED. THIS MODIFICATION WAS
- MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION
------------------------------------------------------------------------
-*/
-{
-static int K1 = 4;
-static int K2 = 8;
-static int K3 = 9;
-static int K4 = 10;
-static double spmpar,b,binv,bm1,one,w,z;
-static int emax,emin,ibeta,m;
-/*
- ..
- .. Executable Statements ..
-*/
- if(*i > 1) goto S10;
- b = ipmpar(&K1);
- m = ipmpar(&K2);
- spmpar = pow(b,(double)(1-m));
- return spmpar;
-S10:
- if(*i > 2) goto S20;
- b = ipmpar(&K1);
- emin = ipmpar(&K3);
- one = 1.0;
- binv = one/b;
- w = pow(b,(double)(emin+2));
- spmpar = w*binv*binv*binv;
- return spmpar;
-S20:
- ibeta = ipmpar(&K1);
- m = ipmpar(&K2);
- emax = ipmpar(&K4);
- b = ibeta;
- bm1 = ibeta-1;
- one = 1.0;
- z = pow(b,(double)(m-1));
- w = ((z-one)*b+bm1)/(b*z);
- z = pow(b,(double)(emax-2));
- spmpar = w*z*b*b;
- return spmpar;
-}
-double stvaln(double *p)
-/*
-**********************************************************************
-
- double stvaln(double *p)
- STarting VALue for Neton-Raphon
- calculation of Normal distribution Inverse
-
-
- Function
-
-
- Returns X such that CUMNOR(X) = P, i.e., the integral from -
- infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
-
-
- Arguments
-
-
- P --> The probability whose normal deviate is sought.
- P is DOUBLE PRECISION
-
-
- Method
-
-
- The rational function on page 95 of Kennedy and Gentle,
- Statistical Computing, Marcel Dekker, NY , 1980.
-
-**********************************************************************
-*/
-{
-static double xden[5] = {
- 0.993484626060e-1,0.588581570495e0,0.531103462366e0,0.103537752850e0,
- 0.38560700634e-2
-};
-static double xnum[5] = {
- -0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1,
- -0.453642210148e-4
-};
-static int K1 = 5;
-static double stvaln,sign,y,z;
-/*
- ..
- .. Executable Statements ..
-*/
- if(!(*p <= 0.5e0)) goto S10;
- sign = -1.0e0;
- z = *p;
- goto S20;
-S10:
- sign = 1.0e0;
- z = 1.0e0-*p;
-S20:
- y = sqrt(-(2.0e0*log(z)));
- stvaln = y+devlpl(xnum,&K1,&y)/devlpl(xden,&K1,&y);
- stvaln = sign*stvaln;
- return stvaln;
-}
-/************************************************************************
-FIFDINT:
-Truncates a double precision number to an integer and returns the
-value in a double.
-************************************************************************/
-double fifdint(double a)
-/* a - number to be truncated */
-{
- long temp;
- temp = (long)(a);
- return (double)(temp);
-}
-/************************************************************************
-FIFDMAX1:
-returns the maximum of two numbers a and b
-************************************************************************/
-double fifdmax1(double a,double b)
-/* a - first number */
-/* b - second number */
-{
- if (a < b) return b;
- else return a;
-}
-/************************************************************************
-FIFDMIN1:
-returns the minimum of two numbers a and b
-************************************************************************/
-double fifdmin1(double a,double b)
-/* a - first number */
-/* b - second number */
-{
- if (a < b) return a;
- else return b;
-}
-/************************************************************************
-FIFDSIGN:
-transfers the sign of the variable "sign" to the variable "mag"
-************************************************************************/
-double fifdsign(double mag,double sign)
-/* mag - magnitude */
-/* sign - sign to be transfered */
-{
- if (mag < 0) mag = -mag;
- if (sign < 0) mag = -mag;
- return mag;
-
-}
-/************************************************************************
-FIFIDINT:
-Truncates a double precision number to a long integer
-************************************************************************/
-long fifidint(double a)
-/* a - number to be truncated */
-{
- return (long)(a);
-}
-/************************************************************************
-FIFMOD:
-returns the modulo of a and b
-************************************************************************/
-long fifmod(long a,long b)
-/* a - numerator */
-/* b - denominator */
-{
- return a % b;
-}
projects / pspp-builds.git / commitdiff