php/archived-pecl-math-stats
1
0
Fork 0
mirror of https://github.com/php/pecl-math-stats.git synced 2026-03-23 23:22:18 +01:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
archived-pecl-math-stats/dcdflib.c
y-uti 7f9445ff2d Fix bug in the initilization of S in cdfbin function 
Closes #12
2016-02-04 00:01:55 +09:00

9086 lines
221 KiB
C
Raw Permalink Blame History

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cdflib.h"
/*
* 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 = *xn;
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) ftnstop(" SMALL, X, BIG not monotone in INVR");
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.
**********************************************************************
*/
{
E0000(1,NULL,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
zstpmu);
}
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;
}
/************************************************************************
FTNSTOP:
Prints msg to standard error and then exits
************************************************************************/
void ftnstop(char* msg)
/* msg - error message */
{
if (msg != NULL) fprintf(stderr,"%s\n",msg);
exit(EXIT_FAILURE); /* EXIT_FAILURE from stdlib.h, or use an int */
}
Reference in New Issue View Git Blame Copy Permalink