10
# include "cdflib.hpp"
12
//****************************************************************************80
14
double algdiv ( double *a, double *b )
16
//****************************************************************************80
20
// ALGDIV computes ln ( Gamma ( B ) / Gamma ( A + B ) ) when 8 <= B.
24
// In this algorithm, DEL(X) is the function defined by
26
// ln ( Gamma(X) ) = ( X - 0.5 ) * ln ( X ) - X + 0.5 * ln ( 2 * PI )
31
// Input, double *A, *B, define the arguments.
33
// Output, double ALGDIV, the value of ln(Gamma(B)/Gamma(A+B)).
38
static double c0 = 0.833333333333333e-01;
39
static double c1 = -0.277777777760991e-02;
40
static double c2 = 0.793650666825390e-03;
41
static double c3 = -0.595202931351870e-03;
42
static double c4 = 0.837308034031215e-03;
43
static double c5 = -0.165322962780713e-02;
62
c = 1.0e0 / ( 1.0e0 + h );
63
x = h / ( 1.0e0 + h );
64
d = *a + ( *b - 0.5e0 );
69
c = h / ( 1.0e0 + h );
70
x = 1.0e0 / ( 1.0e0 + h );
71
d = *b + ( *a - 0.5e0 );
74
// SET SN = (1 - X**N)/(1 - X)
77
s3 = 1.0e0 + ( x + x2 );
78
s5 = 1.0e0 + ( x + x2 * s3 );
79
s7 = 1.0e0 + ( x + x2 * s5 );
80
s9 = 1.0e0 + ( x + x2 * s7 );
81
s11 = 1.0e0 + ( x + x2 * s9 );
83
// SET W = DEL(B) - DEL(A + B)
85
t = pow ( 1.0e0 / *b, 2.0 );
96
// Combine the results.
99
u = d * alnrel ( &T1 );
100
v = *a * ( log ( *b ) - 1.0e0 );
112
//****************************************************************************80
114
double alnrel ( double *a )
116
//****************************************************************************80
120
// ALNREL evaluates the function ln ( 1 + A ).
128
// Armido DiDinato, Alfred Morris,
130
// Significant Digit Computation of the Incomplete Beta Function Ratios,
131
// ACM Transactions on Mathematical Software,
132
// Volume 18, 1993, pages 360-373.
136
// Input, double *A, the argument.
138
// Output, double ALNREL, the value of ln ( 1 + A ).
142
static double p1 = -0.129418923021993e+01;
143
static double p2 = 0.405303492862024e+00;
144
static double p3 = -0.178874546012214e-01;
145
static double q1 = -0.162752256355323e+01;
146
static double q2 = 0.747811014037616e+00;
147
static double q3 = -0.845104217945565e-01;
153
if ( fabs ( *a ) <= 0.375e0 )
155
t = *a / ( *a + 2.0e0 );
157
w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)
158
/ (((q3*t2+q2)*t2+q1)*t2+1.0e0);
159
alnrel = 2.0e0 * t * w;
168
//****************************************************************************80
170
double apser ( double *a, double *b, double *x, double *eps )
172
//****************************************************************************80
176
// APSER computes the incomplete beta ratio I(SUB(1-X))(B,A).
180
// APSER is used only for cases where
182
// A <= min ( EPS, EPS * B ),
188
// Input, double *A, *B, *X, the parameters of teh
189
// incomplete beta ratio.
191
// Input, double *EPS, a tolerance.
193
// Output, double APSER, the computed value of the
194
// incomplete beta ratio.
197
static double g = 0.577215664901533e0;
198
static double apser,aj,bx,c,j,s,t,tol;
202
if(*b**eps > 2.e-2) goto S10;
203
c = log(*x)+psi(b)+g+t;
208
tol = 5.0e0**eps*fabs(c);
216
if(fabs(aj) > tol) goto S30;
220
//****************************************************************************80
222
double bcorr ( double *a0, double *b0 )
224
//****************************************************************************80
228
// BCORR evaluates DEL(A0) + DEL(B0) - DEL(A0 + B0).
232
// The function DEL(A) is a remainder term that is used in the expression:
234
// ln ( Gamma ( A ) ) = ( A - 0.5 ) * ln ( A )
235
// - A + 0.5 * ln ( 2 * PI ) + DEL ( A ),
237
// or, in other words, DEL ( A ) is defined as:
239
// DEL ( A ) = ln ( Gamma ( A ) ) - ( A - 0.5 ) * ln ( A )
240
// + A + 0.5 * ln ( 2 * PI ).
244
// Input, double *A0, *B0, the arguments.
245
// It is assumed that 8 <= A0 and 8 <= B0.
247
// Output, double *BCORR, the value of the function.
250
static double c0 = 0.833333333333333e-01;
251
static double c1 = -0.277777777760991e-02;
252
static double c2 = 0.793650666825390e-03;
253
static double c3 = -0.595202931351870e-03;
254
static double c4 = 0.837308034031215e-03;
255
static double c5 = -0.165322962780713e-02;
256
static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;
258
a = fifdmin1 ( *a0, *b0 );
259
b = fifdmax1 ( *a0, *b0 );
261
c = h / ( 1.0e0 + h );
262
x = 1.0e0 / ( 1.0e0 + h );
265
// SET SN = (1 - X**N)/(1 - X)
267
s3 = 1.0e0 + ( x + x2 );
268
s5 = 1.0e0 + ( x + x2 * s3 );
269
s7 = 1.0e0 + ( x + x2 * s5 );
270
s9 = 1.0e0 + ( x + x2 * s7 );
271
s11 = 1.0e0 + ( x + x2 * s9 );
273
// SET W = DEL(B) - DEL(A + B)
275
t = pow ( 1.0e0 / b, 2.0 );
277
w = (((( c5 * s11 * t + c4
284
// COMPUTE DEL(A) + W
286
t = pow ( 1.0e0 / a, 2.0 );
288
bcorr = ((((( c5 * t + c4 )
295
//****************************************************************************80
297
double beta ( double a, double b )
299
//****************************************************************************80
303
// BETA evaluates the beta function.
315
// Input, double A, B, the arguments of the beta function.
317
// Output, double BETA, the value of the beta function.
320
return ( exp ( beta_log ( &a, &b ) ) );
322
//****************************************************************************80
324
double beta_asym ( double *a, double *b, double *lambda, double *eps )
326
//****************************************************************************80
330
// BETA_ASYM computes an asymptotic expansion for IX(A,B), for large A and B.
334
// Input, double *A, *B, the parameters of the function.
335
// A and B should be nonnegative. It is assumed that both A and B
336
// are greater than or equal to 15.
338
// Input, double *LAMBDA, the value of ( A + B ) * Y - B.
339
// It is assumed that 0 <= LAMBDA.
341
// Input, double *EPS, the tolerance.
344
static double e0 = 1.12837916709551e0;
345
static double e1 = .353553390593274e0;
348
// NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
349
// ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
350
// THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
356
static double bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
358
static int i,im1,imj,j,m,mm1,mmj,n,np1;
359
static double a0[21],b0[21],c[21],d[21],T1,T2;
362
if(*a >= *b) goto S10;
364
r0 = 1.0e0/(1.0e0+h);
366
w0 = 1.0e0/sqrt(*a*(1.0e0+h));
370
r0 = 1.0e0/(1.0e0+h);
372
w0 = 1.0e0/sqrt(*b*(1.0e0+h));
376
f = *a*rlog1(&T1)+*b*rlog1(&T2);
378
if(t == 0.0e0) return value;
382
a0[0] = 2.0e0/3.0e0*r1;
383
c[0] = -(0.5e0*a0[0]);
385
j0 = 0.5e0/e0 * error_fc ( &K3, &z0 );
394
for ( n = 2; n <= num; n += 2 )
397
a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
400
a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
401
for ( i = n; i <= np1; i++ )
403
r = -(0.5e0*((double)i+1.0e0));
405
for ( m = 2; m <= i; m++ )
409
for ( j = 1; j <= mm1; j++ )
412
bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
414
b0[m-1] = r*a0[m-1]+bsum/(double)m;
416
c[i-1] = b0[i-1]/((double)i+1.0e0);
419
for ( j = 1; j <= im1; j++ )
422
dsum += (d[imj-1]*c[j-1]);
424
d[i-1] = -(dsum+c[i-1]);
426
j0 = e1*znm1+((double)n-1.0e0)*j0;
427
j1 = e1*zn+(double)n*j1;
435
if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
438
u = exp(-bcorr(a,b));
442
//****************************************************************************80
444
double beta_frac ( double *a, double *b, double *x, double *y, double *lambda,
447
//****************************************************************************80
451
// BETA_FRAC evaluates a continued fraction expansion for IX(A,B).
455
// Input, double *A, *B, the parameters of the function.
456
// A and B should be nonnegative. It is assumed that both A and
457
// B are greater than 1.
459
// Input, double *X, *Y. X is the argument of the
460
// function, and should satisy 0 <= X <= 1. Y should equal 1 - X.
462
// Input, double *LAMBDA, the value of ( A + B ) * Y - B.
464
// Input, double *EPS, a tolerance.
466
// Output, double BETA_FRAC, the value of the continued
467
// fraction approximation for IX(A,B).
470
static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;
472
bfrac = beta_rcomp ( a, b, x, y );
474
if ( bfrac == 0.0e0 )
481
c1 = 1.0e0+1.0e0/ *a;
491
// CONTINUED FRACTION CALCULATION
498
alpha = p*(p+c0)*e*e*(w**x);
499
e = (1.0e0+t)/(c1+t+t);
500
beta = n+w/s+e*(c+n*yp1);
504
// UPDATE AN, BN, ANP1, AND BNP1
506
t = alpha*an+beta*anp1;
509
t = alpha*bn+beta*bnp1;
515
if ( fabs(r-r0) <= (*eps) * r )
520
// RESCALE AN, BN, ANP1, AND BNP1
534
//****************************************************************************80
536
void beta_grat ( double *a, double *b, double *x, double *y, double *w,
537
double *eps,int *ierr )
539
//****************************************************************************80
543
// BETA_GRAT evaluates an asymptotic expansion for IX(A,B).
547
// Input, double *A, *B, the parameters of the function.
548
// A and B should be nonnegative. It is assumed that 15 <= A
549
// and B <= 1, and that B is less than A.
551
// Input, double *X, *Y. X is the argument of the
552
// function, and should satisy 0 <= X <= 1. Y should equal 1 - X.
554
// Input/output, double *W, a quantity to which the
555
// result of the computation is to be added on output.
557
// Input, double *EPS, a tolerance.
559
// Output, int *IERR, an error flag, which is 0 if no error
563
static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
565
static double c[30],d[30],T1;
567
bm1 = *b-0.5e0-0.5e0;
569
if(*y > 0.375e0) goto S10;
577
if(*b*z == 0.0e0) goto S70;
579
// COMPUTATION OF THE EXPANSION
580
// SET R = EXP(-Z)*Z**B/GAMMA(B)
582
r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
583
r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
584
u = algdiv(b,a)+*b*log(nu);
586
if(u == 0.0e0) goto S70;
587
gamma_rat1 ( b, &z, &r, &p, &q, eps );
588
v = 0.25e0*pow(1.0e0/nu,2.0);
595
for ( n = 1; n <= 30; n++ )
598
j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
601
cn /= (n2*(n2+1.0e0));
607
for ( i = 1; i <= nm1; i++ )
609
s = s + (coef*c[i-1]*d[n-i-1]);
613
d[n-1] = bm1*cn+s/(double)n;
616
if(sum <= 0.0e0) goto S70;
617
if(fabs(dj) <= *eps*(sum+l)) goto S60;
621
// ADD THE RESULTS TO W
628
// THE EXPANSION CANNOT BE COMPUTED
633
//****************************************************************************80
635
void beta_inc ( double *a, double *b, double *x, double *y, double *w,
636
double *w1, int *ierr )
638
//****************************************************************************80
642
// BETA_INC evaluates the incomplete beta function IX(A,B).
646
// Alfred H Morris, Jr,
647
// Naval Surface Weapons Center,
648
// Dahlgren, Virginia.
652
// Input, double *A, *B, the parameters of the function.
653
// A and B should be nonnegative.
655
// Input, double *X, *Y. X is the argument of the
656
// function, and should satisy 0 <= X <= 1. Y should equal 1 - X.
658
// Output, double *W, *W1, the values of IX(A,B) and
661
// Output, int *IERR, the error flag.
662
// 0, no error was detected.
663
// 1, A or B is negative;
665
// 3, X < 0 or 1 < X;
666
// 4, Y < 0 or 1 < Y;
673
static double a0,b0,eps,lambda,t,x0,y0,z;
674
static int ierr1,ind,n;
675
static double T2,T3,T4,T5;
677
// EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
678
// NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
680
eps = dpmpar ( &K1 );
682
if(*a < 0.0e0 || *b < 0.0e0) goto S270;
683
if(*a == 0.0e0 && *b == 0.0e0) goto S280;
684
if(*x < 0.0e0 || *x > 1.0e0) goto S290;
685
if(*y < 0.0e0 || *y > 1.0e0) goto S300;
686
z = *x+*y-0.5e0-0.5e0;
687
if(fabs(z) > 3.0e0*eps) goto S310;
689
if(*x == 0.0e0) goto S210;
690
if(*y == 0.0e0) goto S230;
691
if(*a == 0.0e0) goto S240;
692
if(*b == 0.0e0) goto S220;
693
eps = fifdmax1(eps,1.e-15);
694
if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
700
if(fifdmin1(a0,b0) > 1.0e0) goto S40;
702
// PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
704
if(*x <= 0.5e0) goto S10;
711
if(b0 < fifdmin1(eps,eps*a0)) goto S90;
712
if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
713
if(fifdmax1(a0,b0) > 1.0e0) goto S20;
714
if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
715
if(pow(x0,a0) <= 0.9e0) goto S110;
716
if(x0 >= 0.3e0) goto S120;
720
if(b0 <= 1.0e0) goto S110;
721
if(x0 >= 0.3e0) goto S120;
722
if(x0 >= 0.1e0) goto S30;
723
if(pow(x0*b0,a0) <= 0.7e0) goto S110;
725
if(b0 > 15.0e0) goto S150;
730
// PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
732
if(*a > *b) goto S50;
733
lambda = *a-(*a+*b)**x;
736
lambda = (*a+*b)**y-*b;
738
if(lambda >= 0.0e0) goto S70;
744
lambda = fabs(lambda);
746
if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
747
if(b0 < 40.0e0) goto S160;
748
if(a0 > b0) goto S80;
749
if(a0 <= 100.0e0) goto S130;
750
if(lambda > 0.03e0*a0) goto S130;
753
if(b0 <= 100.0e0) goto S130;
754
if(lambda > 0.03e0*b0) goto S130;
758
// EVALUATION OF THE APPROPRIATE ALGORITHM
760
*w = fpser(&a0,&b0,&x0,&eps);
761
*w1 = 0.5e0+(0.5e0-*w);
764
*w1 = apser(&a0,&b0,&x0,&eps);
765
*w = 0.5e0+(0.5e0-*w1);
768
*w = beta_pser(&a0,&b0,&x0,&eps);
769
*w1 = 0.5e0+(0.5e0-*w);
772
*w1 = beta_pser(&b0,&a0,&y0,&eps);
773
*w = 0.5e0+(0.5e0-*w1);
777
*w = beta_frac ( &a0,&b0,&x0,&y0,&lambda,&T2 );
778
*w1 = 0.5e0+(0.5e0-*w);
781
*w1 = beta_up ( &b0, &a0, &y0, &x0, &n, &eps );
785
beta_grat (&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
786
*w = 0.5e0+(0.5e0-*w1);
791
if(b0 != 0.0e0) goto S170;
795
*w = beta_up ( &b0, &a0, &y0, &x0, &n, &eps );
796
if(x0 > 0.7e0) goto S180;
797
*w = *w + beta_pser(&a0,&b0,&x0,&eps);
798
*w1 = 0.5e0+(0.5e0-*w);
801
if(a0 > 15.0e0) goto S190;
803
*w = *w + beta_up ( &a0, &b0, &x0, &y0, &n, &eps );
807
beta_grat ( &a0, &b0, &x0, &y0, w, &T4, &ierr1 );
808
*w1 = 0.5e0+(0.5e0-*w);
812
*w = beta_asym ( &a0, &b0, &lambda, &T5 );
813
*w1 = 0.5e0+(0.5e0-*w);
817
// TERMINATION OF THE PROCEDURE
819
if(*a == 0.0e0) goto S320;
825
if(*b == 0.0e0) goto S330;
838
// PROCEDURE FOR A AND B .LT. 1.E-3*EPS
868
//****************************************************************************80
870
void beta_inc_values ( int *n_data, double *a, double *b, double *x,
873
//****************************************************************************80
877
// BETA_INC_VALUES returns some values of the incomplete Beta function.
881
// The incomplete Beta function may be written
883
// BETA_INC(A,B,X) = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT
884
// / Integral (0 to 1) T**(A-1) * (1-T)**(B-1) dT
888
// BETA_INC(A,B,0.0) = 0.0
889
// BETA_INC(A,B,1.0) = 1.0
891
// Note that in Mathematica, the expressions:
893
// BETA[A,B] = Integral (0 to 1) T**(A-1) * (1-T)**(B-1) dT
894
// BETA[X,A,B] = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT
896
// and thus, to evaluate the incomplete Beta function requires:
898
// BETA_INC(A,B,X) = BETA[X,A,B] / BETA[A,B]
910
// Milton Abramowitz and Irene Stegun,
911
// Handbook of Mathematical Functions,
912
// US Department of Commerce, 1964.
915
// Tables of the Incomplete Beta Function,
916
// Cambridge University Press, 1968.
920
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
921
// first call. On each call, the routine increments N_DATA by 1, and
922
// returns the corresponding data; when there is no more data, the
923
// output value of N_DATA will be 0 again.
925
// Output, double *A, *B, the parameters of the function.
927
// Output, double *X, the argument of the function.
929
// Output, double *FX, the value of the function.
934
double a_vec[N_MAX] = {
935
0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00,
936
1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00,
937
2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00,
938
2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00,
939
2.0E+00, 5.5E+00, 10.0E+00, 10.0E+00,
940
10.0E+00, 10.0E+00, 20.0E+00, 20.0E+00,
941
20.0E+00, 20.0E+00, 20.0E+00, 30.0E+00,
942
30.0E+00, 40.0E+00 };
943
double b_vec[N_MAX] = {
944
0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00,
945
0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00,
946
2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00,
947
2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00,
948
2.0E+00, 5.0E+00, 0.5E+00, 5.0E+00,
949
5.0E+00, 10.0E+00, 5.0E+00, 10.0E+00,
950
10.0E+00, 20.0E+00, 20.0E+00, 10.0E+00,
951
10.0E+00, 20.0E+00 };
952
double fx_vec[N_MAX] = {
953
0.0637686E+00, 0.2048328E+00, 1.0000000E+00, 0.0E+00,
954
0.0050126E+00, 0.0513167E+00, 0.2928932E+00, 0.5000000E+00,
955
0.028E+00, 0.104E+00, 0.216E+00, 0.352E+00,
956
0.500E+00, 0.648E+00, 0.784E+00, 0.896E+00,
957
0.972E+00, 0.4361909E+00, 0.1516409E+00, 0.0897827E+00,
958
1.0000000E+00, 0.5000000E+00, 0.4598773E+00, 0.2146816E+00,
959
0.9507365E+00, 0.5000000E+00, 0.8979414E+00, 0.2241297E+00,
960
0.7586405E+00, 0.7001783E+00 };
961
double x_vec[N_MAX] = {
962
0.01E+00, 0.10E+00, 1.00E+00, 0.0E+00,
963
0.01E+00, 0.10E+00, 0.50E+00, 0.50E+00,
964
0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00,
965
0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00,
966
0.9E+00, 0.50E+00, 0.90E+00, 0.50E+00,
967
1.00E+00, 0.50E+00, 0.80E+00, 0.60E+00,
968
0.80E+00, 0.50E+00, 0.60E+00, 0.70E+00,
969
0.80E+00, 0.70E+00 };
976
*n_data = *n_data + 1;
978
if ( N_MAX < *n_data )
988
*a = a_vec[*n_data-1];
989
*b = b_vec[*n_data-1];
990
*x = x_vec[*n_data-1];
991
*fx = fx_vec[*n_data-1];
996
//****************************************************************************80
998
double beta_log ( double *a0, double *b0 )
1000
//****************************************************************************80
1004
// BETA_LOG evaluates the logarithm of the beta function.
1008
// Armido DiDinato and Alfred Morris,
1010
// Significant Digit Computation of the Incomplete Beta Function Ratios,
1011
// ACM Transactions on Mathematical Software,
1012
// Volume 18, 1993, pages 360-373.
1016
// Input, double *A0, *B0, the parameters of the function.
1017
// A0 and B0 should be nonnegative.
1019
// Output, double *BETA_LOG, the value of the logarithm
1020
// of the Beta function.
1023
static double e = .918938533204673e0;
1024
static double value,a,b,c,h,u,v,w,z;
1028
a = fifdmin1(*a0,*b0);
1029
b = fifdmax1(*a0,*b0);
1030
if(a >= 8.0e0) goto S100;
1031
if(a >= 1.0e0) goto S20;
1033
// PROCEDURE WHEN A .LT. 1
1035
if(b >= 8.0e0) goto S10;
1037
value = gamma_log ( &a )+( gamma_log ( &b )- gamma_log ( &T1 ));
1040
value = gamma_log ( &a )+algdiv(&a,&b);
1044
// PROCEDURE WHEN 1 .LE. A .LT. 8
1046
if(a > 2.0e0) goto S40;
1047
if(b > 2.0e0) goto S30;
1048
value = gamma_log ( &a )+ gamma_log ( &b )-gsumln(&a,&b);
1052
if(b < 8.0e0) goto S60;
1053
value = gamma_log ( &a )+algdiv(&a,&b);
1057
// REDUCTION OF A WHEN B .LE. 1000
1059
if(b > 1000.0e0) goto S80;
1060
n = ( int ) ( a - 1.0e0 );
1062
for ( i = 1; i <= n; i++ )
1069
if(b < 8.0e0) goto S60;
1070
value = w+ gamma_log ( &a )+algdiv(&a,&b);
1074
// REDUCTION OF B WHEN B .LT. 8
1076
n = ( int ) ( b - 1.0e0 );
1078
for ( i = 1; i <= n; i++ )
1083
value = w+log(z)+( gamma_log ( &a )+( gamma_log ( &b )-gsumln(&a,&b)));
1087
// REDUCTION OF A WHEN B .GT. 1000
1089
n = ( int ) ( a - 1.0e0 );
1091
for ( i = 1; i <= n; i++ )
1094
w *= (a/(1.0e0+a/b));
1096
value = log(w)-(double)n*log(b)+( gamma_log ( &a )+algdiv(&a,&b));
1100
// PROCEDURE WHEN A .GE. 8
1105
u = -((a-0.5e0)*log(c));
1107
if(u <= v) goto S110;
1108
value = -(0.5e0*log(b))+e+w-v-u;
1111
value = -(0.5e0*log(b))+e+w-u-v;
1114
//****************************************************************************80
1116
double beta_pser ( double *a, double *b, double *x, double *eps )
1118
//****************************************************************************80
1122
// BETA_PSER uses a power series expansion to evaluate IX(A,B)(X).
1126
// BETA_PSER is used when B <= 1 or B*X <= 0.7.
1130
// Input, double *A, *B, the parameters.
1132
// Input, double *X, the point where the function
1133
// is to be evaluated.
1135
// Input, double *EPS, the tolerance.
1137
// Output, double BETA_PSER, the approximate value of IX(A,B)(X).
1140
static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
1144
if(*x == 0.0e0) return bpser;
1146
// COMPUTE THE FACTOR X**A/(A*BETA(A,B))
1148
a0 = fifdmin1(*a,*b);
1149
if(a0 < 1.0e0) goto S10;
1150
z = *a*log(*x)-beta_log(a,b);
1154
b0 = fifdmax1(*a,*b);
1155
if(b0 >= 8.0e0) goto S90;
1156
if(b0 > 1.0e0) goto S40;
1158
// PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
1161
if(bpser == 0.0e0) return bpser;
1163
if(apb > 1.0e0) goto S20;
1164
z = 1.0e0+gam1(&apb);
1168
z = (1.0e0+gam1(&u))/apb;
1170
c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1171
bpser *= (c*(*b/apb));
1175
// PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
1177
u = gamma_ln1 ( &a0 );
1178
m = ( int ) ( b0 - 1.0e0 );
1181
for ( i = 1; i <= m; i++ )
1191
if(apb > 1.0e0) goto S70;
1192
t = 1.0e0+gam1(&apb);
1196
t = (1.0e0+gam1(&u))/apb;
1198
bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
1202
// PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
1204
u = gamma_ln1 ( &a0 ) + algdiv ( &a0, &b0 );
1206
bpser = a0/ *a*exp(z);
1208
if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
1210
// COMPUTE THE SERIES
1217
c *= ((0.5e0+(0.5e0-*b/n))**x);
1220
if(fabs(w) > tol) goto S110;
1221
bpser *= (1.0e0+*a*sum);
1224
//****************************************************************************80
1226
double beta_rcomp ( double *a, double *b, double *x, double *y )
1228
//****************************************************************************80
1232
// BETA_RCOMP evaluates X**A * Y**B / Beta(A,B).
1236
// Input, double *A, *B, the parameters of the Beta function.
1237
// A and B should be nonnegative.
1239
// Input, double *X, *Y, define the numerator of the fraction.
1241
// Output, double BETA_RCOMP, the value of X**A * Y**B / Beta(A,B).
1244
static double Const = .398942280401433e0;
1245
static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
1248
// CONST = 1/SQRT(2*PI)
1250
static double T1,T2;
1253
if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
1254
a0 = fifdmin1(*a,*b);
1255
if(a0 >= 8.0e0) goto S130;
1256
if(*x > 0.375e0) goto S10;
1262
if(*y > 0.375e0) goto S20;
1272
if(a0 < 1.0e0) goto S40;
1278
// PROCEDURE FOR A .LT. 1 OR B .LT. 1
1280
b0 = fifdmax1(*a,*b);
1281
if(b0 >= 8.0e0) goto S120;
1282
if(b0 > 1.0e0) goto S70;
1284
// ALGORITHM FOR B0 .LE. 1
1287
if(brcomp == 0.0e0) return brcomp;
1289
if(apb > 1.0e0) goto S50;
1290
z = 1.0e0+gam1(&apb);
1294
z = (1.0e0+gam1(&u))/apb;
1296
c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1297
brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
1301
// ALGORITHM FOR 1 .LT. B0 .LT. 8
1303
u = gamma_ln1 ( &a0 );
1304
n = ( int ) ( b0 - 1.0e0 );
1307
for ( i = 1; i <= n; i++ )
1317
if(apb > 1.0e0) goto S100;
1318
t = 1.0e0+gam1(&apb);
1322
t = (1.0e0+gam1(&u))/apb;
1324
brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
1328
// ALGORITHM FOR B0 .GE. 8
1330
u = gamma_ln1 ( &a0 ) + algdiv ( &a0, &b0 );
1331
brcomp = a0*exp(z-u);
1335
// PROCEDURE FOR A .GE. 8 AND B .GE. 8
1337
if(*a > *b) goto S140;
1340
y0 = 1.0e0/(1.0e0+h);
1341
lambda = *a-(*a+*b)**x;
1345
x0 = 1.0e0/(1.0e0+h);
1347
lambda = (*a+*b)**y-*b;
1350
if(fabs(e) > 0.6e0) goto S160;
1357
if(fabs(e) > 0.6e0) goto S180;
1363
z = exp(-(*a*u+*b*v));
1364
brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
1367
//****************************************************************************80
1369
double beta_rcomp1 ( int *mu, double *a, double *b, double *x, double *y )
1371
//****************************************************************************80
1375
// BETA_RCOMP1 evaluates exp(MU) * X**A * Y**B / Beta(A,B).
1381
// Input, double A, B, the parameters of the Beta function.
1382
// A and B should be nonnegative.
1384
// Input, double X, Y, ?
1386
// Output, double BETA_RCOMP1, the value of
1387
// exp(MU) * X**A * Y**B / Beta(A,B).
1390
static double Const = .398942280401433e0;
1391
static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
1394
// CONST = 1/SQRT(2*PI)
1396
static double T1,T2,T3,T4;
1398
a0 = fifdmin1(*a,*b);
1399
if(a0 >= 8.0e0) goto S130;
1400
if(*x > 0.375e0) goto S10;
1406
if(*y > 0.375e0) goto S20;
1416
if(a0 < 1.0e0) goto S40;
1418
brcmp1 = esum(mu,&z);
1422
// PROCEDURE FOR A .LT. 1 OR B .LT. 1
1424
b0 = fifdmax1(*a,*b);
1425
if(b0 >= 8.0e0) goto S120;
1426
if(b0 > 1.0e0) goto S70;
1428
// ALGORITHM FOR B0 .LE. 1
1430
brcmp1 = esum(mu,&z);
1431
if(brcmp1 == 0.0e0) return brcmp1;
1433
if(apb > 1.0e0) goto S50;
1434
z = 1.0e0+gam1(&apb);
1438
z = (1.0e0+gam1(&u))/apb;
1440
c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1441
brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
1445
// ALGORITHM FOR 1 .LT. B0 .LT. 8
1447
u = gamma_ln1 ( &a0 );
1448
n = ( int ) ( b0 - 1.0e0 );
1451
for ( i = 1; i <= n; i++ )
1461
if(apb > 1.0e0) goto S100;
1462
t = 1.0e0+gam1(&apb);
1466
t = (1.0e0+gam1(&u))/apb;
1468
brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
1472
// ALGORITHM FOR B0 .GE. 8
1474
u = gamma_ln1 ( &a0 ) + algdiv ( &a0, &b0 );
1476
brcmp1 = a0*esum(mu,&T3);
1480
// PROCEDURE FOR A .GE. 8 AND B .GE. 8
1482
if(*a > *b) goto S140;
1485
y0 = 1.0e0/(1.0e0+h);
1486
lambda = *a-(*a+*b)**x;
1490
x0 = 1.0e0/(1.0e0+h);
1492
lambda = (*a+*b)**y-*b;
1495
if(fabs(e) > 0.6e0) goto S160;
1502
if(fabs(e) > 0.6e0) goto S180;
1510
brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
1513
//****************************************************************************80
1515
double beta_up ( double *a, double *b, double *x, double *y, int *n,
1518
//****************************************************************************80
1522
// BETA_UP evaluates IX(A,B) - IX(A+N,B) where N is a positive integer.
1526
// Input, double *A, *B, the parameters of the function.
1527
// A and B should be nonnegative.
1529
// Input, double *X, *Y, ?
1533
// Input, double *EPS, the tolerance.
1535
// Output, double BETA_UP, the value of IX(A,B) - IX(A+N,B).
1540
static double bup,ap1,apb,d,l,r,t,w;
1541
static int i,k,kp1,mu,nm1;
1543
// OBTAIN THE SCALING FACTOR EXP(-MU) AND
1544
// EXP(MU)*(X**A*Y**B/BETA(A,B))/A
1550
if(*n == 1 || *a < 1.0e0) goto S10;
1551
if(apb < 1.1e0*ap1) goto S10;
1552
mu = ( int ) fabs ( exparg(&K1) );
1553
k = ( int ) exparg ( &K2 );
1558
bup = beta_rcomp1 ( &mu, a, b, x, y ) / *a;
1559
if(*n == 1 || bup == 0.0e0) return bup;
1563
// LET K BE THE INDEX OF THE MAXIMUM TERM
1566
if(*b <= 1.0e0) goto S50;
1567
if(*y > 1.e-4) goto S20;
1571
r = (*b-1.0e0)**x/ *y-*a;
1572
if(r < 1.0e0) goto S50;
1575
if ( r < t ) k = ( int ) r;
1578
// ADD THE INCREASING TERMS OF THE SERIES
1580
for ( i = 1; i <= k; i++ )
1583
d = (apb+l)/(ap1+l)**x*d;
1586
if(k == nm1) goto S70;
1589
// ADD THE REMAINING TERMS OF THE SERIES
1592
for ( i = kp1; i <= nm1; i++ )
1595
d = (apb+l)/(ap1+l)**x*d;
1597
if(d <= *eps*w) goto S70;
1601
// TERMINATE THE PROCEDURE
1606
//****************************************************************************80
1608
void binomial_cdf_values ( int *n_data, int *a, double *b, int *x, double *fx )
1610
//****************************************************************************80
1614
// BINOMIAL_CDF_VALUES returns some values of the binomial CDF.
1618
// CDF(X)(A,B) is the probability of at most X successes in A trials,
1619
// given that the probability of success on a single trial is B.
1631
// Milton Abramowitz and Irene Stegun,
1632
// Handbook of Mathematical Functions,
1633
// US Department of Commerce, 1964.
1635
// Daniel Zwillinger,
1636
// CRC Standard Mathematical Tables and Formulae,
1637
// 30th Edition, CRC Press, 1996, pages 651-652.
1641
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
1642
// first call. On each call, the routine increments N_DATA by 1, and
1643
// returns the corresponding data; when there is no more data, the
1644
// output value of N_DATA will be 0 again.
1646
// Output, int *A, double *B, the parameters of the function.
1648
// Output, int *X, the argument of the function.
1650
// Output, double *FX, the value of the function.
1655
int a_vec[N_MAX] = {
1661
double b_vec[N_MAX] = {
1662
0.05E+00, 0.05E+00, 0.05E+00, 0.50E+00,
1663
0.50E+00, 0.25E+00, 0.25E+00, 0.25E+00,
1664
0.25E+00, 0.05E+00, 0.10E+00, 0.15E+00,
1665
0.20E+00, 0.25E+00, 0.30E+00, 0.40E+00,
1667
double fx_vec[N_MAX] = {
1668
0.9025E+00, 0.9975E+00, 1.0000E+00, 0.2500E+00,
1669
0.7500E+00, 0.3164E+00, 0.7383E+00, 0.9492E+00,
1670
0.9961E+00, 0.9999E+00, 0.9984E+00, 0.9901E+00,
1671
0.9672E+00, 0.9219E+00, 0.8497E+00, 0.6331E+00,
1673
int x_vec[N_MAX] = {
1685
*n_data = *n_data + 1;
1687
if ( N_MAX < *n_data )
1697
*a = a_vec[*n_data-1];
1698
*b = b_vec[*n_data-1];
1699
*x = x_vec[*n_data-1];
1700
*fx = fx_vec[*n_data-1];
1705
//****************************************************************************80
1707
void cdfbet ( int *which, double *p, double *q, double *x, double *y,
1708
double *a, double *b, int *status, double *bound )
1710
//****************************************************************************80
1714
// CDFBET evaluates the CDF of the Beta Distribution.
1718
// This routine calculates any one parameter of the beta distribution
1719
// given the others.
1721
// The value P of the cumulative distribution function is calculated
1722
// directly by code associated with the reference.
1724
// Computation of the other parameters involves a seach for a value that
1725
// produces the desired value of P. The search relies on the
1726
// monotonicity of P with respect to the other parameters.
1728
// The beta density is proportional to t^(A-1) * (1-t)^(B-1).
1736
// Armido DiDinato and Alfred Morris,
1738
// Significant Digit Computation of the Incomplete Beta Function Ratios,
1739
// ACM Transactions on Mathematical Software,
1740
// Volume 18, 1993, pages 360-373.
1744
// Input, int *WHICH, indicates which of the next four argument
1745
// values is to be calculated from the others.
1746
// 1: Calculate P and Q from X, Y, A and B;
1747
// 2: Calculate X and Y from P, Q, A and B;
1748
// 3: Calculate A from P, Q, X, Y and B;
1749
// 4: Calculate B from P, Q, X, Y and A.
1751
// Input/output, double *P, the integral from 0 to X of the
1752
// chi-square distribution. Input range: [0, 1].
1754
// Input/output, double *Q, equals 1-P. Input range: [0, 1].
1756
// Input/output, double *X, the upper limit of integration
1757
// of the beta density. If it is an input value, it should lie in
1758
// the range [0,1]. If it is an output value, it will be searched for
1759
// in the range [0,1].
1761
// Input/output, double *Y, equal to 1-X. If it is an input
1762
// value, it should lie in the range [0,1]. If it is an output value,
1763
// it will be searched for in the range [0,1].
1765
// Input/output, double *A, the first parameter of the beta
1766
// density. If it is an input value, it should lie in the range
1767
// (0, +infinity). If it is an output value, it will be searched
1768
// for in the range [1D-300,1D300].
1770
// Input/output, double *B, the second parameter of the beta
1771
// density. If it is an input value, it should lie in the range
1772
// (0, +infinity). If it is an output value, it will be searched
1773
// for in the range [1D-300,1D300].
1775
// Output, int *STATUS, reports the status of the computation.
1776
// 0, if the calculation completed correctly;
1777
// -I, if the input parameter number I is out of range;
1778
// +1, if the answer appears to be lower than lowest search bound;
1779
// +2, if the answer appears to be higher than greatest search bound;
1780
// +3, if P + Q /= 1;
1781
// +4, if X + Y /= 1.
1783
// Output, double *BOUND, is only defined if STATUS is nonzero.
1784
// If STATUS is negative, then this is the value exceeded by parameter I.
1785
// if STATUS is 1 or 2, this is the search bound that was exceeded.
1788
# define tol (1.0e-8)
1789
# define atol (1.0e-50)
1790
# define zero (1.0e-300)
1791
# define inf 1.0e300
1795
static double K2 = 0.0e0;
1796
static double K3 = 1.0e0;
1797
static double K8 = 0.5e0;
1798
static double K9 = 5.0e0;
1799
static double fx,xhi,xlo,cum,ccum,xy,pq;
1800
static unsigned long qhi,qleft,qporq;
1801
static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
1808
if(!(*which < 1 || *which > 4)) goto S30;
1809
if(!(*which < 1)) goto S10;
1818
if(*which == 1) goto S70;
1822
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
1823
if(!(*p < 0.0e0)) goto S40;
1833
if(*which == 1) goto S110;
1837
if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
1838
if(!(*q < 0.0e0)) goto S80;
1848
if(*which == 2) goto S150;
1852
if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
1853
if(!(*x < 0.0e0)) goto S120;
1863
if(*which == 2) goto S190;
1867
if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
1868
if(!(*y < 0.0e0)) goto S160;
1878
if(*which == 3) goto S210;
1882
if(!(*a <= 0.0e0)) goto S200;
1888
if(*which == 4) goto S230;
1892
if(!(*b <= 0.0e0)) goto S220;
1898
if(*which == 1) goto S270;
1903
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S260;
1904
if(!(pq < 0.0e0)) goto S240;
1914
if(*which == 2) goto S310;
1919
if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S300;
1920
if(!(xy < 0.0e0)) goto S280;
1930
if(!(*which == 1)) qporq = *p <= *q;
1932
// Select the minimum of P or Q
1933
// Calculate ANSWERS
1937
// Calculating P and Q
1939
cumbet(x,y,a,b,p,q);
1942
else if(2 == *which) {
1944
// Calculating X and Y
1948
dstzr(&K2,&K3,&T4,&T5);
1949
if(!qporq) goto S340;
1951
dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
1954
if(!(*status == 1)) goto S330;
1955
cumbet(x,y,a,b,&cum,&ccum);
1957
dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
1964
dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
1967
if(!(*status == 1)) goto S360;
1968
cumbet(x,y,a,b,&cum,&ccum);
1970
dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
1975
if(!(*status == -1)) goto S400;
1976
if(!qleft) goto S380;
1987
else if(3 == *which) {
1996
dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
1998
dinvr(status,a,&fx,&qleft,&qhi);
2000
if(!(*status == 1)) goto S440;
2001
cumbet(x,y,a,b,&cum,&ccum);
2002
if(!qporq) goto S420;
2008
dinvr(status,a,&fx,&qleft,&qhi);
2011
if(!(*status == -1)) goto S470;
2012
if(!qleft) goto S450;
2023
else if(4 == *which) {
2032
dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
2034
dinvr(status,b,&fx,&qleft,&qhi);
2036
if(!(*status == 1)) goto S510;
2037
cumbet(x,y,a,b,&cum,&ccum);
2038
if(!qporq) goto S490;
2044
dinvr(status,b,&fx,&qleft,&qhi);
2047
if(!(*status == -1)) goto S540;
2048
if(!qleft) goto S520;
2066
//****************************************************************************80
2068
void cdfbin ( int *which, double *p, double *q, double *s, double *xn,
2069
double *pr, double *ompr, int *status, double *bound )
2071
//****************************************************************************80
2075
// CDFBIN evaluates the CDF of the Binomial distribution.
2079
// This routine calculates any one parameter of the binomial distribution
2080
// given the others.
2082
// The value P of the cumulative distribution function is calculated
2085
// Computation of the other parameters involves a seach for a value that
2086
// produces the desired value of P. The search relies on the
2087
// monotonicity of P with respect to the other parameters.
2089
// P is the probablility of S or fewer successes in XN binomial trials,
2090
// each trial having an individual probability of success of PR.
2098
// Milton Abramowitz and Irene Stegun,
2099
// Handbook of Mathematical Functions
2100
// 1966, Formula 26.5.24.
2104
// Input, int *WHICH, indicates which of argument values is to
2105
// be calculated from the others.
2106
// 1: Calculate P and Q from S, XN, PR and OMPR;
2107
// 2: Calculate S from P, Q, XN, PR and OMPR;
2108
// 3: Calculate XN from P, Q, S, PR and OMPR;
2109
// 4: Calculate PR and OMPR from P, Q, S and XN.
2111
// Input/output, double *P, the cumulation, from 0 to S,
2112
// of the binomial distribution. If P is an input value, it should
2113
// lie in the range [0,1].
2115
// Input/output, double *Q, equal to 1-P. If Q is an input
2116
// value, it should lie in the range [0,1]. If Q is an output value,
2117
// it will lie in the range [0,1].
2119
// Input/output, double *S, the number of successes observed.
2120
// Whether this is an input or output value, it should lie in the
2123
// Input/output, double *XN, the number of binomial trials.
2124
// If this is an input value it should lie in the range: (0, +infinity).
2125
// If it is an output value it will be searched for in the
2126
// range [1.0D-300, 1.0D+300].
2128
// Input/output, double *PR, the probability of success in each
2129
// binomial trial. Whether this is an input or output value, it should
2130
// lie in the range: [0,1].
2132
// Input/output, double *OMPR, equal to 1-PR. Whether this is an
2133
// input or output value, it should lie in the range [0,1]. Also, it should
2134
// be the case that PR + OMPR = 1.
2136
// Output, int *STATUS, reports the status of the computation.
2137
// 0, if the calculation completed correctly;
2138
// -I, if the input parameter number I is out of range;
2139
// +1, if the answer appears to be lower than lowest search bound;
2140
// +2, if the answer appears to be higher than greatest search bound;
2141
// +3, if P + Q /= 1;
2142
// +4, if PR + OMPR /= 1.
2144
// Output, double *BOUND, is only defined if STATUS is nonzero.
2145
// If STATUS is negative, then this is the value exceeded by parameter I.
2146
// if STATUS is 1 or 2, this is the search bound that was exceeded.
2149
# define atol (1.0e-50)
2150
# define tol (1.0e-8)
2151
# define zero (1.0e-300)
2152
# define inf 1.0e300
2156
static double K2 = 0.0e0;
2157
static double K3 = 0.5e0;
2158
static double K4 = 5.0e0;
2159
static double K11 = 1.0e0;
2160
static double fx,xhi,xlo,cum,ccum,pq,prompr;
2161
static unsigned long qhi,qleft,qporq;
2162
static double T5,T6,T7,T8,T9,T10,T12,T13;
2169
if(!(*which < 1 && *which > 4)) goto S30;
2170
if(!(*which < 1)) goto S10;
2179
if(*which == 1) goto S70;
2183
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
2184
if(!(*p < 0.0e0)) goto S40;
2194
if(*which == 1) goto S110;
2198
if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
2199
if(!(*q < 0.0e0)) goto S80;
2209
if(*which == 3) goto S130;
2213
if(!(*xn <= 0.0e0)) goto S120;
2219
if(*which == 2) goto S170;
2223
if(!(*s < 0.0e0 || *which != 3 && *s > *xn)) goto S160;
2224
if(!(*s < 0.0e0)) goto S140;
2234
if(*which == 4) goto S210;
2238
if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S200;
2239
if(!(*pr < 0.0e0)) goto S180;
2249
if(*which == 4) goto S250;
2253
if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S240;
2254
if(!(*ompr < 0.0e0)) goto S220;
2264
if(*which == 1) goto S290;
2269
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S280;
2270
if(!(pq < 0.0e0)) goto S260;
2280
if(*which == 4) goto S330;
2285
if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S320;
2286
if(!(prompr < 0.0e0)) goto S300;
2296
if(!(*which == 1)) qporq = *p <= *q;
2298
// Select the minimum of P or Q
2299
// Calculate ANSWERS
2305
cumbin(s,xn,pr,ompr,p,q);
2308
else if(2 == *which) {
2315
dstinv(&K2,xn,&K3,&K3,&K4,&T5,&T6);
2317
dinvr(status,s,&fx,&qleft,&qhi);
2319
if(!(*status == 1)) goto S370;
2320
cumbin(s,xn,pr,ompr,&cum,&ccum);
2321
if(!qporq) goto S350;
2327
dinvr(status,s,&fx,&qleft,&qhi);
2330
if(!(*status == -1)) goto S400;
2331
if(!qleft) goto S380;
2342
else if(3 == *which) {
2351
dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
2353
dinvr(status,xn,&fx,&qleft,&qhi);
2355
if(!(*status == 1)) goto S440;
2356
cumbin(s,xn,pr,ompr,&cum,&ccum);
2357
if(!qporq) goto S420;
2363
dinvr(status,xn,&fx,&qleft,&qhi);
2366
if(!(*status == -1)) goto S470;
2367
if(!qleft) goto S450;
2378
else if(4 == *which) {
2380
// Calculating PR and OMPR
2384
dstzr(&K2,&K11,&T12,&T13);
2385
if(!qporq) goto S500;
2387
dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
2390
if(!(*status == 1)) goto S490;
2391
cumbin(s,xn,pr,ompr,&cum,&ccum);
2393
dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
2400
dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
2403
if(!(*status == 1)) goto S520;
2404
cumbin(s,xn,pr,ompr,&cum,&ccum);
2406
dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
2411
if(!(*status == -1)) goto S560;
2412
if(!qleft) goto S540;
2430
//****************************************************************************80
2432
void cdfchi ( int *which, double *p, double *q, double *x, double *df,
2433
int *status, double *bound )
2435
//****************************************************************************80
2439
// CDFCHI evaluates the CDF of the chi square distribution.
2443
// This routine calculates any one parameter of the chi square distribution
2444
// given the others.
2446
// The value P of the cumulative distribution function is calculated
2449
// Computation of the other parameters involves a seach for a value that
2450
// produces the desired value of P. The search relies on the
2451
// monotonicity of P with respect to the other parameters.
2453
// The CDF of the chi square distribution can be evaluated
2454
// within Mathematica by commands such as:
2456
// Needs["Statistics`ContinuousDistributions`"]
2457
// CDF [ ChiSquareDistribution [ DF ], X ]
2461
// Milton Abramowitz and Irene Stegun,
2462
// Handbook of Mathematical Functions
2463
// 1966, Formula 26.4.19.
2466
// The Mathematica Book,
2468
// Wolfram Media / Cambridge University Press, 1999.
2472
// Input, int *WHICH, indicates which argument is to be calculated
2474
// 1: Calculate P and Q from X and DF;
2475
// 2: Calculate X from P, Q and DF;
2476
// 3: Calculate DF from P, Q and X.
2478
// Input/output, double *P, the integral from 0 to X of
2479
// the chi-square distribution. If this is an input value, it should
2480
// lie in the range [0,1].
2482
// Input/output, double *Q, equal to 1-P. If Q is an input
2483
// value, it should lie in the range [0,1]. If Q is an output value,
2484
// it will lie in the range [0,1].
2486
// Input/output, double *X, the upper limit of integration
2487
// of the chi-square distribution. If this is an input
2488
// value, it should lie in the range: [0, +infinity). If it is an output
2489
// value, it will be searched for in the range: [0,1.0D+300].
2491
// Input/output, double *DF, the degrees of freedom of the
2492
// chi-square distribution. If this is an input value, it should lie
2493
// in the range: (0, +infinity). If it is an output value, it will be
2494
// searched for in the range: [ 1.0D-300, 1.0D+300].
2496
// Output, int *STATUS, reports the status of the computation.
2497
// 0, if the calculation completed correctly;
2498
// -I, if the input parameter number I is out of range;
2499
// +1, if the answer appears to be lower than lowest search bound;
2500
// +2, if the answer appears to be higher than greatest search bound;
2501
// +3, if P + Q /= 1;
2502
// +10, an error was returned from CUMGAM.
2504
// Output, double *BOUND, is only defined if STATUS is nonzero.
2505
// If STATUS is negative, then this is the value exceeded by parameter I.
2506
// if STATUS is 1 or 2, this is the search bound that was exceeded.
2509
# define tol (1.0e-8)
2510
# define atol (1.0e-50)
2511
# define zero (1.0e-300)
2512
# define inf 1.0e300
2515
static double K2 = 0.0e0;
2516
static double K4 = 0.5e0;
2517
static double K5 = 5.0e0;
2518
static double fx,cum,ccum,pq,porq;
2519
static unsigned long qhi,qleft,qporq;
2520
static double T3,T6,T7,T8,T9,T10,T11;
2527
if(!(*which < 1 || *which > 3)) goto S30;
2528
if(!(*which < 1)) goto S10;
2537
if(*which == 1) goto S70;
2541
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
2542
if(!(*p < 0.0e0)) goto S40;
2552
if(*which == 1) goto S110;
2556
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
2557
if(!(*q <= 0.0e0)) goto S80;
2567
if(*which == 2) goto S130;
2571
if(!(*x < 0.0e0)) goto S120;
2577
if(*which == 3) goto S150;
2581
if(!(*df <= 0.0e0)) goto S140;
2587
if(*which == 1) goto S190;
2592
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S180;
2593
if(!(pq < 0.0e0)) goto S160;
2603
if(*which == 1) goto S220;
2605
// Select the minimum of P or Q
2608
if(!qporq) goto S200;
2616
// Calculate ANSWERS
2620
// Calculating P and Q
2629
else if(2 == *which) {
2637
dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
2639
dinvr(status,x,&fx,&qleft,&qhi);
2641
if(!(*status == 1)) goto S270;
2642
cumchi(x,df,&cum,&ccum);
2643
if(!qporq) goto S240;
2649
if(!(fx+porq > 1.5e0)) goto S260;
2653
dinvr(status,x,&fx,&qleft,&qhi);
2656
if(!(*status == -1)) goto S300;
2657
if(!qleft) goto S280;
2668
else if(3 == *which) {
2677
dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
2679
dinvr(status,df,&fx,&qleft,&qhi);
2681
if(!(*status == 1)) goto S350;
2682
cumchi(x,df,&cum,&ccum);
2683
if(!qporq) goto S320;
2689
if(!(fx+porq > 1.5e0)) goto S340;
2693
dinvr(status,df,&fx,&qleft,&qhi);
2696
if(!(*status == -1)) goto S380;
2697
if(!qleft) goto S360;
2714
//****************************************************************************80
2716
void cdfchn ( int *which, double *p, double *q, double *x, double *df,
2717
double *pnonc, int *status, double *bound )
2719
//****************************************************************************80
2723
// CDFCHN evaluates the CDF of the Noncentral Chi-Square.
2727
// This routine calculates any one parameter of the noncentral chi-square
2728
// distribution given values for the others.
2730
// The value P of the cumulative distribution function is calculated
2733
// Computation of the other parameters involves a seach for a value that
2734
// produces the desired value of P. The search relies on the
2735
// monotonicity of P with respect to the other parameters.
2737
// The computation time required for this routine is proportional
2738
// to the noncentrality parameter (PNONC). Very large values of
2739
// this parameter can consume immense computer resources. This is
2740
// why the search range is bounded by 10,000.
2742
// The CDF of the noncentral chi square distribution can be evaluated
2743
// within Mathematica by commands such as:
2745
// Needs["Statistics`ContinuousDistributions`"]
2746
// CDF[ NoncentralChiSquareDistribution [ DF, LAMBDA ], X ]
2750
// Milton Abramowitz and Irene Stegun,
2751
// Handbook of Mathematical Functions
2752
// 1966, Formula 26.5.25.
2755
// The Mathematica Book,
2757
// Wolfram Media / Cambridge University Press, 1999.
2761
// Input, int *WHICH, indicates which argument is to be calculated
2763
// 1: Calculate P and Q from X, DF and PNONC;
2764
// 2: Calculate X from P, DF and PNONC;
2765
// 3: Calculate DF from P, X and PNONC;
2766
// 4: Calculate PNONC from P, X and DF.
2768
// Input/output, double *P, the integral from 0 to X of
2769
// the noncentral chi-square distribution. If this is an input
2770
// value, it should lie in the range: [0, 1.0-1.0D-16).
2772
// Input/output, double *Q, is generally not used by this
2773
// subroutine and is only included for similarity with other routines.
2774
// However, if P is to be computed, then a value will also be computed
2777
// Input, double *X, the upper limit of integration of the
2778
// noncentral chi-square distribution. If this is an input value, it
2779
// should lie in the range: [0, +infinity). If it is an output value,
2780
// it will be sought in the range: [0,1.0D+300].
2782
// Input/output, double *DF, the number of degrees of freedom
2783
// of the noncentral chi-square distribution. If this is an input value,
2784
// it should lie in the range: (0, +infinity). If it is an output value,
2785
// it will be searched for in the range: [ 1.0D-300, 1.0D+300].
2787
// Input/output, double *PNONC, the noncentrality parameter of
2788
// the noncentral chi-square distribution. If this is an input value, it
2789
// should lie in the range: [0, +infinity). If it is an output value,
2790
// it will be searched for in the range: [0,1.0D+4]
2792
// Output, int *STATUS, reports on the calculation.
2793
// 0, if calculation completed correctly;
2794
// -I, if input parameter number I is out of range;
2795
// 1, if the answer appears to be lower than the lowest search bound;
2796
// 2, if the answer appears to be higher than the greatest search bound.
2798
// Output, double *BOUND, is only defined if STATUS is nonzero.
2799
// If STATUS is negative, then this is the value exceeded by parameter I.
2800
// if STATUS is 1 or 2, this is the search bound that was exceeded.
2803
# define tent4 1.0e4
2804
# define tol (1.0e-8)
2805
# define atol (1.0e-50)
2806
# define zero (1.0e-300)
2807
# define one (1.0e0-1.0e-16)
2808
# define inf 1.0e300
2810
static double K1 = 0.0e0;
2811
static double K3 = 0.5e0;
2812
static double K4 = 5.0e0;
2813
static double fx,cum,ccum;
2814
static unsigned long qhi,qleft;
2815
static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;
2822
if(!(*which < 1 || *which > 4)) goto S30;
2823
if(!(*which < 1)) goto S10;
2832
if(*which == 1) goto S70;
2836
if(!(*p < 0.0e0 || *p > one)) goto S60;
2837
if(!(*p < 0.0e0)) goto S40;
2847
if(*which == 2) goto S90;
2851
if(!(*x < 0.0e0)) goto S80;
2857
if(*which == 3) goto S110;
2861
if(!(*df <= 0.0e0)) goto S100;
2867
if(*which == 4) goto S130;
2871
if(!(*pnonc < 0.0e0)) goto S120;
2878
// Calculate ANSWERS
2882
// Calculating P and Q
2884
cumchn(x,df,pnonc,p,q);
2887
else if(2 == *which) {
2895
dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
2897
dinvr(status,x,&fx,&qleft,&qhi);
2899
if(!(*status == 1)) goto S150;
2900
cumchn(x,df,pnonc,&cum,&ccum);
2902
dinvr(status,x,&fx,&qleft,&qhi);
2905
if(!(*status == -1)) goto S180;
2906
if(!qleft) goto S160;
2917
else if(3 == *which) {
2926
dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
2928
dinvr(status,df,&fx,&qleft,&qhi);
2930
if(!(*status == 1)) goto S200;
2931
cumchn(x,df,pnonc,&cum,&ccum);
2933
dinvr(status,df,&fx,&qleft,&qhi);
2936
if(!(*status == -1)) goto S230;
2937
if(!qleft) goto S210;
2948
else if(4 == *which) {
2950
// Calculating PNONC
2956
dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);
2958
dinvr(status,pnonc,&fx,&qleft,&qhi);
2960
if(!(*status == 1)) goto S250;
2961
cumchn(x,df,pnonc,&cum,&ccum);
2963
dinvr(status,pnonc,&fx,&qleft,&qhi);
2966
if(!(*status == -1)) goto S280;
2967
if(!qleft) goto S260;
2986
//****************************************************************************80
2988
void cdff ( int *which, double *p, double *q, double *f, double *dfn,
2989
double *dfd, int *status, double *bound )
2991
//****************************************************************************80
2995
// CDFF evaluates the CDF of the F distribution.
2999
// This routine calculates any one parameter of the F distribution
3000
// given the others.
3002
// The value P of the cumulative distribution function is calculated
3005
// Computation of the other parameters involves a seach for a value that
3006
// produces the desired value of P. The search relies on the
3007
// monotonicity of P with respect to the other parameters.
3009
// The value of the cumulative F distribution is not necessarily
3010
// monotone in either degree of freedom. There thus may be two
3011
// values that provide a given CDF value. This routine assumes
3012
// monotonicity and will find an arbitrary one of the two values.
3020
// Milton Abramowitz, Irene Stegun,
3021
// Handbook of Mathematical Functions
3022
// 1966, Formula 26.6.2.
3026
// Input, int *WHICH, indicates which argument is to be calculated
3028
// 1: Calculate P and Q from F, DFN and DFD;
3029
// 2: Calculate F from P, Q, DFN and DFD;
3030
// 3: Calculate DFN from P, Q, F and DFD;
3031
// 4: Calculate DFD from P, Q, F and DFN.
3033
// Input/output, double *P, the integral from 0 to F of
3034
// the F-density. If it is an input value, it should lie in the
3037
// Input/output, double *Q, equal to 1-P. If Q is an input
3038
// value, it should lie in the range [0,1]. If Q is an output value,
3039
// it will lie in the range [0,1].
3041
// Input/output, double *F, the upper limit of integration
3042
// of the F-density. If this is an input value, it should lie in the
3043
// range [0, +infinity). If it is an output value, it will be searched
3044
// for in the range [0,1.0D+300].
3046
// Input/output, double *DFN, the number of degrees of
3047
// freedom of the numerator sum of squares. If this is an input value,
3048
// it should lie in the range: (0, +infinity). If it is an output value,
3049
// it will be searched for in the range: [ 1.0D-300, 1.0D+300].
3051
// Input/output, double *DFD, the number of degrees of freedom
3052
// of the denominator sum of squares. If this is an input value, it should
3053
// lie in the range: (0, +infinity). If it is an output value, it will
3054
// be searched for in the range: [ 1.0D-300, 1.0D+300].
3056
// Output, int *STATUS, reports the status of the computation.
3057
// 0, if the calculation completed correctly;
3058
// -I, if the input parameter number I is out of range;
3059
// +1, if the answer appears to be lower than lowest search bound;
3060
// +2, if the answer appears to be higher than greatest search bound;
3061
// +3, if P + Q /= 1.
3063
// Output, double *BOUND, is only defined if STATUS is nonzero.
3064
// If STATUS is negative, then this is the value exceeded by parameter I.
3065
// if STATUS is 1 or 2, this is the search bound that was exceeded.
3068
# define tol (1.0e-8)
3069
# define atol (1.0e-50)
3070
# define zero (1.0e-300)
3071
# define inf 1.0e300
3074
static double K2 = 0.0e0;
3075
static double K4 = 0.5e0;
3076
static double K5 = 5.0e0;
3077
static double pq,fx,cum,ccum;
3078
static unsigned long qhi,qleft,qporq;
3079
static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
3086
if(!(*which < 1 || *which > 4)) goto S30;
3087
if(!(*which < 1)) goto S10;
3096
if(*which == 1) goto S70;
3100
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
3101
if(!(*p < 0.0e0)) goto S40;
3111
if(*which == 1) goto S110;
3115
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
3116
if(!(*q <= 0.0e0)) goto S80;
3126
if(*which == 2) goto S130;
3130
if(!(*f < 0.0e0)) goto S120;
3136
if(*which == 3) goto S150;
3140
if(!(*dfn <= 0.0e0)) goto S140;
3146
if(*which == 4) goto S170;
3150
if(!(*dfd <= 0.0e0)) goto S160;
3156
if(*which == 1) goto S210;
3161
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S200;
3162
if(!(pq < 0.0e0)) goto S180;
3172
if(!(*which == 1)) qporq = *p <= *q;
3174
// Select the minimum of P or Q
3175
// Calculate ANSWERS
3181
cumf(f,dfn,dfd,p,q);
3184
else if(2 == *which) {
3192
dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
3194
dinvr(status,f,&fx,&qleft,&qhi);
3196
if(!(*status == 1)) goto S250;
3197
cumf(f,dfn,dfd,&cum,&ccum);
3198
if(!qporq) goto S230;
3204
dinvr(status,f,&fx,&qleft,&qhi);
3207
if(!(*status == -1)) goto S280;
3208
if(!qleft) goto S260;
3222
// Note that, in the original calculation, the lower bound for DFN was 0.
3223
// Using DFN = 0 causes an error in CUMF when it calls BETA_INC.
3224
// The lower bound was set to the more reasonable value of 1.
3225
// JVB, 14 April 2007.
3227
else if ( 3 == *which )
3234
dstinv ( &T8, &T9, &K4, &K4, &K5, &T10, &T11 );
3240
dinvr ( status, dfn, &fx, &qleft, &qhi );
3242
while ( *status == 1 )
3244
cumf ( f, dfn, dfd, &cum, &ccum );
3254
dinvr ( status, dfn, &fx, &qleft, &qhi );
3257
if ( *status == -1 )
3274
// Note that, in the original calculation, the lower bound for DFD was 0.
3275
// Using DFD = 0 causes an error in CUMF when it calls BETA_INC.
3276
// The lower bound was set to the more reasonable value of 1.
3277
// JVB, 14 April 2007.
3280
else if ( 4 == *which )
3287
dstinv ( &T12, &T13, &K4, &K4, &K5, &T14, &T15 );
3292
dinvr ( status, dfd, &fx, &qleft, &qhi );
3294
while ( *status == 1 )
3296
cumf ( f, dfn, dfd, &cum, &ccum );
3306
dinvr ( status, dfd, &fx, &qleft, &qhi );
3309
if ( *status == -1 )
3330
//****************************************************************************80
3332
void cdffnc ( int *which, double *p, double *q, double *f, double *dfn,
3333
double *dfd, double *phonc, int *status, double *bound )
3335
//****************************************************************************80
3339
// CDFFNC evaluates the CDF of the Noncentral F distribution.
3343
// This routine originally used 1.0E+300 as the upper bound for the
3344
// interval in which many of the missing parameters are to be sought.
3345
// Since the underlying rootfinder routine needs to evaluate the
3346
// function at this point, it is no surprise that the program was
3347
// experiencing overflows. A less extravagant upper bound
3348
// is being tried for now!
3351
// This routine calculates any one parameter of the Noncentral F distribution
3352
// given the others.
3354
// The value P of the cumulative distribution function is calculated
3357
// Computation of the other parameters involves a seach for a value that
3358
// produces the desired value of P. The search relies on the
3359
// monotonicity of P with respect to the other parameters.
3361
// The computation time required for this routine is proportional
3362
// to the noncentrality parameter PNONC. Very large values of
3363
// this parameter can consume immense computer resources. This is
3364
// why the search range is bounded by 10,000.
3366
// The value of the cumulative noncentral F distribution is not
3367
// necessarily monotone in either degree of freedom. There thus
3368
// may be two values that provide a given CDF value. This routine
3369
// assumes monotonicity and will find an arbitrary one of the two
3372
// The CDF of the noncentral F distribution can be evaluated
3373
// within Mathematica by commands such as:
3375
// Needs["Statistics`ContinuousDistributions`"]
3376
// CDF [ NoncentralFRatioDistribution [ DFN, DFD, PNONC ], X ]
3384
// Milton Abramowitz and Irene Stegun,
3385
// Handbook of Mathematical Functions
3386
// 1966, Formula 26.6.20.
3389
// The Mathematica Book,
3391
// Wolfram Media / Cambridge University Press, 1999.
3395
// Input, int *WHICH, indicates which argument is to be calculated
3397
// 1: Calculate P and Q from F, DFN, DFD and PNONC;
3398
// 2: Calculate F from P, Q, DFN, DFD and PNONC;
3399
// 3: Calculate DFN from P, Q, F, DFD and PNONC;
3400
// 4: Calculate DFD from P, Q, F, DFN and PNONC;
3401
// 5: Calculate PNONC from P, Q, F, DFN and DFD.
3403
// Input/output, double *P, the integral from 0 to F of
3404
// the noncentral F-density. If P is an input value it should
3405
// lie in the range [0,1) (Not including 1!).
3407
// Dummy, double *Q, is not used by this subroutine,
3408
// and is only included for similarity with the other routines.
3409
// Its input value is not checked. If P is to be computed, the
3410
// Q is set to 1 - P.
3412
// Input/output, double *F, the upper limit of integration
3413
// of the noncentral F-density. If this is an input value, it should
3414
// lie in the range: [0, +infinity). If it is an output value, it
3415
// will be searched for in the range: [0,1.0D+30].
3417
// Input/output, double *DFN, the number of degrees of freedom
3418
// of the numerator sum of squares. If this is an input value, it should
3419
// lie in the range: (0, +infinity). If it is an output value, it will
3420
// be searched for in the range: [ 1.0, 1.0D+30].
3422
// Input/output, double *DFD, the number of degrees of freedom
3423
// of the denominator sum of squares. If this is an input value, it should
3424
// be in range: (0, +infinity). If it is an output value, it will be
3425
// searched for in the range [1.0, 1.0D+30].
3427
// Input/output, double *PNONC, the noncentrality parameter
3428
// If this is an input value, it should be nonnegative.
3429
// If it is an output value, it will be searched for in the range: [0,1.0D+4].
3431
// Output, int *STATUS, reports the status of the computation.
3432
// 0, if the calculation completed correctly;
3433
// -I, if the input parameter number I is out of range;
3434
// +1, if the answer appears to be lower than lowest search bound;
3435
// +2, if the answer appears to be higher than greatest search bound;
3436
// +3, if P + Q /= 1.
3438
// Output, double *BOUND, is only defined if STATUS is nonzero.
3439
// If STATUS is negative, then this is the value exceeded by parameter I.
3440
// if STATUS is 1 or 2, this is the search bound that was exceeded.
3443
# define tent4 1.0e4
3444
# define tol (1.0e-8)
3445
# define atol (1.0e-50)
3446
# define zero (1.0e-300)
3447
# define one (1.0e0-1.0e-16)
3448
# define inf 1.0e300
3450
static double K1 = 0.0e0;
3451
static double K3 = 0.5e0;
3452
static double K4 = 5.0e0;
3453
static double fx,cum,ccum;
3454
static unsigned long qhi,qleft;
3455
static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17;
3462
if(!(*which < 1 || *which > 5)) goto S30;
3463
if(!(*which < 1)) goto S10;
3472
if(*which == 1) goto S70;
3476
if(!(*p < 0.0e0 || *p > one)) goto S60;
3477
if(!(*p < 0.0e0)) goto S40;
3487
if(*which == 2) goto S90;
3491
if(!(*f < 0.0e0)) goto S80;
3497
if(*which == 3) goto S110;
3501
if(!(*dfn <= 0.0e0)) goto S100;
3507
if(*which == 4) goto S130;
3511
if(!(*dfd <= 0.0e0)) goto S120;
3517
if(*which == 5) goto S150;
3521
if(!(*phonc < 0.0e0)) goto S140;
3528
// Calculate ANSWERS
3534
cumfnc(f,dfn,dfd,phonc,p,q);
3537
else if(2 == *which) {
3545
dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
3547
dinvr(status,f,&fx,&qleft,&qhi);
3549
if(!(*status == 1)) goto S170;
3550
cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3552
dinvr(status,f,&fx,&qleft,&qhi);
3555
if(!(*status == -1)) goto S200;
3556
if(!qleft) goto S180;
3567
else if(3 == *which) {
3576
dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
3578
dinvr(status,dfn,&fx,&qleft,&qhi);
3580
if(!(*status == 1)) goto S220;
3581
cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3583
dinvr(status,dfn,&fx,&qleft,&qhi);
3586
if(!(*status == -1)) goto S250;
3587
if(!qleft) goto S230;
3598
else if(4 == *which) {
3607
dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
3609
dinvr(status,dfd,&fx,&qleft,&qhi);
3611
if(!(*status == 1)) goto S270;
3612
cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3614
dinvr(status,dfd,&fx,&qleft,&qhi);
3617
if(!(*status == -1)) goto S300;
3618
if(!qleft) goto S280;
3629
else if(5 == *which) {
3631
// Calculating PHONC
3637
dstinv(&K1,&T15,&K3,&K3,&K4,&T16,&T17);
3639
dinvr(status,phonc,&fx,&qleft,&qhi);
3641
if(!(*status == 1)) goto S320;
3642
cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3644
dinvr(status,phonc,&fx,&qleft,&qhi);
3647
if(!(*status == -1)) goto S350;
3648
if(!qleft) goto S330;
3667
//****************************************************************************80
3669
void cdfgam ( int *which, double *p, double *q, double *x, double *shape,
3670
double *scale, int *status, double *bound )
3672
//****************************************************************************80
3676
// CDFGAM evaluates the CDF of the Gamma Distribution.
3680
// This routine calculates any one parameter of the Gamma distribution
3681
// given the others.
3683
// The cumulative distribution function P is calculated directly.
3685
// Computation of the other parameters involves a seach for a value that
3686
// produces the desired value of P. The search relies on the
3687
// monotonicity of P with respect to the other parameters.
3689
// The gamma density is proportional to T**(SHAPE - 1) * EXP(- SCALE * T)
3693
// Armido DiDinato and Alfred Morris,
3694
// Computation of the incomplete gamma function ratios and their inverse,
3695
// ACM Transactions on Mathematical Software,
3696
// Volume 12, 1986, pages 377-393.
3700
// Input, int *WHICH, indicates which argument is to be calculated
3702
// 1: Calculate P and Q from X, SHAPE and SCALE;
3703
// 2: Calculate X from P, Q, SHAPE and SCALE;
3704
// 3: Calculate SHAPE from P, Q, X and SCALE;
3705
// 4: Calculate SCALE from P, Q, X and SHAPE.
3707
// Input/output, double *P, the integral from 0 to X of the
3708
// Gamma density. If this is an input value, it should lie in the
3711
// Input/output, double *Q, equal to 1-P. If Q is an input
3712
// value, it should lie in the range [0,1]. If Q is an output value,
3713
// it will lie in the range [0,1].
3715
// Input/output, double *X, the upper limit of integration of
3716
// the Gamma density. If this is an input value, it should lie in the
3717
// range: [0, +infinity). If it is an output value, it will lie in
3718
// the range: [0,1E300].
3720
// Input/output, double *SHAPE, the shape parameter of the
3721
// Gamma density. If this is an input value, it should lie in the range:
3722
// (0, +infinity). If it is an output value, it will be searched for
3723
// in the range: [1.0D-300,1.0D+300].
3725
// Input/output, double *SCALE, the scale parameter of the
3726
// Gamma density. If this is an input value, it should lie in the range
3727
// (0, +infinity). If it is an output value, it will be searched for
3728
// in the range: (1.0D-300,1.0D+300].
3730
// Output, int *STATUS, reports the status of the computation.
3731
// 0, if the calculation completed correctly;
3732
// -I, if the input parameter number I is out of range;
3733
// +1, if the answer appears to be lower than lowest search bound;
3734
// +2, if the answer appears to be higher than greatest search bound;
3735
// +3, if P + Q /= 1;
3736
// +10, if the Gamma or inverse Gamma routine cannot compute the answer.
3737
// This usually happens only for X and SHAPE very large (more than 1.0D+10.
3739
// Output, double *BOUND, is only defined if STATUS is nonzero.
3740
// If STATUS is negative, then this is the value exceeded by parameter I.
3741
// if STATUS is 1 or 2, this is the search bound that was exceeded.
3744
# define tol (1.0e-8)
3745
# define atol (1.0e-50)
3746
# define zero (1.0e-300)
3747
# define inf 1.0e300
3750
static double K5 = 0.5e0;
3751
static double K6 = 5.0e0;
3752
static double xx,fx,xscale,cum,ccum,pq,porq;
3754
static unsigned long qhi,qleft,qporq;
3755
static double T2,T3,T4,T7,T8,T9;
3762
if(!(*which < 1 || *which > 4)) goto S30;
3763
if(!(*which < 1)) goto S10;
3772
if(*which == 1) goto S70;
3776
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
3777
if(!(*p < 0.0e0)) goto S40;
3787
if(*which == 1) goto S110;
3791
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
3792
if(!(*q <= 0.0e0)) goto S80;
3802
if(*which == 2) goto S130;
3806
if(!(*x < 0.0e0)) goto S120;
3812
if(*which == 3) goto S150;
3816
if(!(*shape <= 0.0e0)) goto S140;
3822
if(*which == 4) goto S170;
3826
if(!(*scale <= 0.0e0)) goto S160;
3832
if(*which == 1) goto S210;
3837
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S200;
3838
if(!(pq < 0.0e0)) goto S180;
3848
if(*which == 1) goto S240;
3850
// Select the minimum of P or Q
3853
if(!qporq) goto S220;
3861
// Calculate ANSWERS
3869
cumgam(&xscale,shape,p,q);
3870
if(porq > 1.5e0) *status = 10;
3872
else if(2 == *which) {
3877
gamma_inc_inv ( shape, &xx, &T2, p, q, &ierr );
3887
else if(3 == *which) {
3897
dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
3899
dinvr(status,shape,&fx,&qleft,&qhi);
3901
if(!(*status == 1)) goto S290;
3902
cumgam(&xscale,shape,&cum,&ccum);
3903
if(!qporq) goto S260;
3909
if(!(qporq && cum > 1.5e0 || !qporq && ccum > 1.5e0)) goto S280;
3913
dinvr(status,shape,&fx,&qleft,&qhi);
3916
if(!(*status == -1)) goto S320;
3917
if(!qleft) goto S300;
3928
else if(4 == *which) {
3933
gamma_inc_inv ( shape, &xx, &T9, p, q, &ierr );
3949
//****************************************************************************80
3951
void cdfnbn ( int *which, double *p, double *q, double *s, double *xn,
3952
double *pr, double *ompr, int *status, double *bound )
3954
//****************************************************************************80
3958
// CDFNBN evaluates the CDF of the Negative Binomial distribution
3962
// This routine calculates any one parameter of the negative binomial
3963
// distribution given values for the others.
3965
// The cumulative negative binomial distribution returns the
3966
// probability that there will be F or fewer failures before the
3967
// S-th success in binomial trials each of which has probability of
3970
// The individual term of the negative binomial is the probability of
3971
// F failures before S successes and is
3972
// Choose( F, S+F-1 ) * PR^(S) * (1-PR)^F
3974
// Computation of other parameters involve a seach for a value that
3975
// produces the desired value of P. The search relies on the
3976
// monotonicity of P with respect to the other parameters.
3980
// Milton Abramowitz and Irene Stegun,
3981
// Handbook of Mathematical Functions
3982
// 1966, Formula 26.5.26.
3986
// Input, int WHICH, indicates which argument is to be calculated
3988
// 1: Calculate P and Q from F, S, PR and OMPR;
3989
// 2: Calculate F from P, Q, S, PR and OMPR;
3990
// 3: Calculate S from P, Q, F, PR and OMPR;
3991
// 4: Calculate PR and OMPR from P, Q, F and S.
3993
// Input/output, double P, the cumulation from 0 to F of
3994
// the negative binomial distribution. If P is an input value, it
3995
// should lie in the range [0,1].
3997
// Input/output, double Q, equal to 1-P. If Q is an input
3998
// value, it should lie in the range [0,1]. If Q is an output value,
3999
// it will lie in the range [0,1].
4001
// Input/output, double F, the upper limit of cumulation of
4002
// the binomial distribution. There are F or fewer failures before
4003
// the S-th success. If this is an input value, it may lie in the
4004
// range [0,+infinity), and if it is an output value, it will be searched
4005
// for in the range [0,1.0D+300].
4007
// Input/output, double S, the number of successes.
4008
// If this is an input value, it should lie in the range: [0, +infinity).
4009
// If it is an output value, it will be searched for in the range:
4012
// Input/output, double PR, the probability of success in each
4013
// binomial trial. Whether an input or output value, it should lie in the
4016
// Input/output, double OMPR, the value of (1-PR). Whether an
4017
// input or output value, it should lie in the range [0,1].
4019
// Output, int STATUS, reports the status of the computation.
4020
// 0, if the calculation completed correctly;
4021
// -I, if the input parameter number I is out of range;
4022
// +1, if the answer appears to be lower than lowest search bound;
4023
// +2, if the answer appears to be higher than greatest search bound;
4024
// +3, if P + Q /= 1;
4025
// +4, if PR + OMPR /= 1.
4027
// Output, double BOUND, is only defined if STATUS is nonzero.
4028
// If STATUS is negative, then this is the value exceeded by parameter I.
4029
// if STATUS is 1 or 2, this is the search bound that was exceeded.
4032
# define tol (1.0e-8)
4033
# define atol (1.0e-50)
4034
# define inf 1.0e300
4038
static double K2 = 0.0e0;
4039
static double K4 = 0.5e0;
4040
static double K5 = 5.0e0;
4041
static double K11 = 1.0e0;
4042
static double fx,xhi,xlo,pq,prompr,cum,ccum;
4043
static unsigned long qhi,qleft,qporq;
4044
static double T3,T6,T7,T8,T9,T10,T12,T13;
4051
if(!(*which < 1 || *which > 4)) goto S30;
4052
if(!(*which < 1)) goto S10;
4061
if(*which == 1) goto S70;
4065
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
4066
if(!(*p < 0.0e0)) goto S40;
4076
if(*which == 1) goto S110;
4080
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4081
if(!(*q <= 0.0e0)) goto S80;
4091
if(*which == 2) goto S130;
4095
if(!(*s < 0.0e0)) goto S120;
4101
if(*which == 3) goto S150;
4105
if(!(*xn < 0.0e0)) goto S140;
4111
if(*which == 4) goto S190;
4115
if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
4116
if(!(*pr < 0.0e0)) goto S160;
4126
if(*which == 4) goto S230;
4130
if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
4131
if(!(*ompr < 0.0e0)) goto S200;
4141
if(*which == 1) goto S270;
4146
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S260;
4147
if(!(pq < 0.0e0)) goto S240;
4157
if(*which == 4) goto S310;
4162
if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S300;
4163
if(!(prompr < 0.0e0)) goto S280;
4173
if(!(*which == 1)) qporq = *p <= *q;
4175
// Select the minimum of P or Q
4176
// Calculate ANSWERS
4182
cumnbn(s,xn,pr,ompr,p,q);
4185
else if(2 == *which) {
4193
dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
4195
dinvr(status,s,&fx,&qleft,&qhi);
4197
if(!(*status == 1)) goto S350;
4198
cumnbn(s,xn,pr,ompr,&cum,&ccum);
4199
if(!qporq) goto S330;
4205
dinvr(status,s,&fx,&qleft,&qhi);
4208
if(!(*status == -1)) goto S380;
4209
if(!qleft) goto S360;
4220
else if(3 == *which) {
4228
dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
4230
dinvr(status,xn,&fx,&qleft,&qhi);
4232
if(!(*status == 1)) goto S420;
4233
cumnbn(s,xn,pr,ompr,&cum,&ccum);
4234
if(!qporq) goto S400;
4240
dinvr(status,xn,&fx,&qleft,&qhi);
4243
if(!(*status == -1)) goto S450;
4244
if(!qleft) goto S430;
4255
else if(4 == *which) {
4257
// Calculating PR and OMPR
4261
dstzr(&K2,&K11,&T12,&T13);
4262
if(!qporq) goto S480;
4264
dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
4267
if(!(*status == 1)) goto S470;
4268
cumnbn(s,xn,pr,ompr,&cum,&ccum);
4270
dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
4277
dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
4280
if(!(*status == 1)) goto S500;
4281
cumnbn(s,xn,pr,ompr,&cum,&ccum);
4283
dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
4288
if(!(*status == -1)) goto S540;
4289
if(!qleft) goto S520;
4306
//****************************************************************************80
4308
void cdfnor ( int *which, double *p, double *q, double *x, double *mean,
4309
double *sd, int *status, double *bound )
4311
//****************************************************************************80
4315
// CDFNOR evaluates the CDF of the Normal distribution.
4319
// A slightly modified version of ANORM from SPECFUN
4320
// is used to calculate the cumulative standard normal distribution.
4322
// The rational functions from pages 90-95 of Kennedy and Gentle
4323
// are used as starting values to Newton's Iterations which
4324
// compute the inverse standard normal. Therefore no searches are
4325
// necessary for any parameter.
4327
// For X < -15, the asymptotic expansion for the normal is used as
4328
// the starting value in finding the inverse standard normal.
4330
// The normal density is proportional to
4331
// exp( - 0.5D+00 * (( X - MEAN)/SD)**2)
4335
// Milton Abramowitz and Irene Stegun,
4336
// Handbook of Mathematical Functions
4337
// 1966, Formula 26.2.12.
4340
// Algorithm 715: SPECFUN - A Portable FORTRAN Package of
4341
// Special Function Routines and Test Drivers,
4342
// ACM Transactions on Mathematical Software,
4343
// Volume 19, pages 22-32, 1993.
4345
// Kennedy and Gentle,
4346
// Statistical Computing,
4347
// Marcel Dekker, NY, 1980,
4352
// Input, int *WHICH, indicates which argument is to be calculated
4354
// 1: Calculate P and Q from X, MEAN and SD;
4355
// 2: Calculate X from P, Q, MEAN and SD;
4356
// 3: Calculate MEAN from P, Q, X and SD;
4357
// 4: Calculate SD from P, Q, X and MEAN.
4359
// Input/output, double *P, the integral from -infinity to X
4360
// of the Normal density. If this is an input or output value, it will
4361
// lie in the range [0,1].
4363
// Input/output, double *Q, equal to 1-P. If Q is an input
4364
// value, it should lie in the range [0,1]. If Q is an output value,
4365
// it will lie in the range [0,1].
4367
// Input/output, double *X, the upper limit of integration of
4368
// the Normal density.
4370
// Input/output, double *MEAN, the mean of the Normal density.
4372
// Input/output, double *SD, the standard deviation of the
4373
// Normal density. If this is an input value, it should lie in the
4374
// range (0,+infinity).
4376
// Output, int *STATUS, the status of the calculation.
4377
// 0, if calculation completed correctly;
4378
// -I, if input parameter number I is out of range;
4379
// 1, if answer appears to be lower than lowest search bound;
4380
// 2, if answer appears to be higher than greatest search bound;
4381
// 3, if P + Q /= 1.
4383
// Output, double *BOUND, is only defined if STATUS is nonzero.
4384
// If STATUS is negative, then this is the value exceeded by parameter I.
4385
// if STATUS is 1 or 2, this is the search bound that was exceeded.
4397
if(!(*which < 1 || *which > 4)) goto S30;
4398
if(!(*which < 1)) goto S10;
4407
if(*which == 1) goto S70;
4411
if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
4412
if(!(*p <= 0.0e0)) goto S40;
4422
if(*which == 1) goto S110;
4426
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4427
if(!(*q <= 0.0e0)) goto S80;
4437
if(*which == 1) goto S150;
4442
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S140;
4443
if(!(pq < 0.0e0)) goto S120;
4453
if(*which == 4) goto S170;
4457
if(!(*sd <= 0.0e0)) goto S160;
4464
// Calculate ANSWERS
4470
z = (*x-*mean)/ *sd;
4473
else if(2 == *which) {
4480
else if(3 == *which) {
4482
// Computing the MEAN
4487
else if(4 == *which) {
4496
//****************************************************************************80
4498
void cdfpoi ( int *which, double *p, double *q, double *s, double *xlam,
4499
int *status, double *bound )
4501
//****************************************************************************80
4505
// CDFPOI evaluates the CDF of the Poisson distribution.
4509
// This routine calculates any one parameter of the Poisson distribution
4510
// given the others.
4512
// The value P of the cumulative distribution function is calculated
4515
// Computation of other parameters involve a seach for a value that
4516
// produces the desired value of P. The search relies on the
4517
// monotonicity of P with respect to the other parameters.
4521
// Milton Abramowitz and Irene Stegun,
4522
// Handbook of Mathematical Functions
4523
// 1966, Formula 26.4.21.
4527
// Input, int *WHICH, indicates which argument is to be calculated
4529
// 1: Calculate P and Q from S and XLAM;
4530
// 2: Calculate A from P, Q and XLAM;
4531
// 3: Calculate XLAM from P, Q and S.
4533
// Input/output, double *P, the cumulation from 0 to S of the
4534
// Poisson density. Whether this is an input or output value, it will
4535
// lie in the range [0,1].
4537
// Input/output, double *Q, equal to 1-P. If Q is an input
4538
// value, it should lie in the range [0,1]. If Q is an output value,
4539
// it will lie in the range [0,1].
4541
// Input/output, double *S, the upper limit of cumulation of
4542
// the Poisson CDF. If this is an input value, it should lie in
4543
// the range: [0, +infinity). If it is an output value, it will be
4544
// searched for in the range: [0,1.0D+300].
4546
// Input/output, double *XLAM, the mean of the Poisson
4547
// distribution. If this is an input value, it should lie in the range
4548
// [0, +infinity). If it is an output value, it will be searched for
4549
// in the range: [0,1E300].
4551
// Output, int *STATUS, reports the status of the computation.
4552
// 0, if the calculation completed correctly;
4553
// -I, if the input parameter number I is out of range;
4554
// +1, if the answer appears to be lower than lowest search bound;
4555
// +2, if the answer appears to be higher than greatest search bound;
4556
// +3, if P + Q /= 1.
4558
// Output, double *BOUND, is only defined if STATUS is nonzero.
4559
// If STATUS is negative, then this is the value exceeded by parameter I.
4560
// if STATUS is 1 or 2, this is the search bound that was exceeded.
4563
# define tol (1.0e-8)
4564
# define atol (1.0e-50)
4565
# define inf 1.0e300
4568
static double K2 = 0.0e0;
4569
static double K4 = 0.5e0;
4570
static double K5 = 5.0e0;
4571
static double fx,cum,ccum,pq;
4572
static unsigned long qhi,qleft,qporq;
4573
static double T3,T6,T7,T8,T9,T10;
4580
if(!(*which < 1 || *which > 3)) goto S30;
4581
if(!(*which < 1)) goto S10;
4590
if(*which == 1) goto S70;
4594
if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
4595
if(!(*p < 0.0e0)) goto S40;
4605
if(*which == 1) goto S110;
4609
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4610
if(!(*q <= 0.0e0)) goto S80;
4620
if(*which == 2) goto S130;
4624
if(!(*s < 0.0e0)) goto S120;
4630
if(*which == 3) goto S150;
4634
if(!(*xlam < 0.0e0)) goto S140;
4640
if(*which == 1) goto S190;
4645
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S180;
4646
if(!(pq < 0.0e0)) goto S160;
4656
if(!(*which == 1)) qporq = *p <= *q;
4658
// Select the minimum of P or Q
4659
// Calculate ANSWERS
4668
else if(2 == *which) {
4676
dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
4678
dinvr(status,s,&fx,&qleft,&qhi);
4680
if(!(*status == 1)) goto S230;
4681
cumpoi(s,xlam,&cum,&ccum);
4682
if(!qporq) goto S210;
4688
dinvr(status,s,&fx,&qleft,&qhi);
4691
if(!(*status == -1)) goto S260;
4692
if(!qleft) goto S240;
4703
else if(3 == *which) {
4711
dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
4713
dinvr(status,xlam,&fx,&qleft,&qhi);
4715
if(!(*status == 1)) goto S300;
4716
cumpoi(s,xlam,&cum,&ccum);
4717
if(!qporq) goto S280;
4723
dinvr(status,xlam,&fx,&qleft,&qhi);
4726
if(!(*status == -1)) goto S330;
4727
if(!qleft) goto S310;
4743
//****************************************************************************80
4745
void cdft ( int *which, double *p, double *q, double *t, double *df,
4746
int *status, double *bound )
4748
//****************************************************************************80
4752
// CDFT evaluates the CDF of the T distribution.
4756
// This routine calculates any one parameter of the T distribution
4757
// given the others.
4759
// The value P of the cumulative distribution function is calculated
4762
// Computation of other parameters involve a seach for a value that
4763
// produces the desired value of P. The search relies on the
4764
// monotonicity of P with respect to the other parameters.
4766
// The original version of this routine allowed the search interval
4767
// to extend from -1.0E+300 to +1.0E+300, which is fine until you
4768
// try to evaluate a function at such a point!
4772
// Milton Abramowitz and Irene Stegun,
4773
// Handbook of Mathematical Functions
4774
// 1966, Formula 26.5.27.
4778
// Input, int *WHICH, indicates which argument is to be calculated
4780
// 1 : Calculate P and Q from T and DF;
4781
// 2 : Calculate T from P, Q and DF;
4782
// 3 : Calculate DF from P, Q and T.
4784
// Input/output, double *P, the integral from -infinity to T of
4785
// the T-density. Whether an input or output value, this will lie in the
4788
// Input/output, double *Q, equal to 1-P. If Q is an input
4789
// value, it should lie in the range [0,1]. If Q is an output value,
4790
// it will lie in the range [0,1].
4792
// Input/output, double *T, the upper limit of integration of
4793
// the T-density. If this is an input value, it may have any value.
4794
// It it is an output value, it will be searched for in the range
4795
// [ -1.0D+30, 1.0D+30 ].
4797
// Input/output, double *DF, the number of degrees of freedom
4798
// of the T distribution. If this is an input value, it should lie
4799
// in the range: (0 , +infinity). If it is an output value, it will be
4800
// searched for in the range: [1, 1.0D+10].
4802
// Output, int *STATUS, reports the status of the computation.
4803
// 0, if the calculation completed correctly;
4804
// -I, if the input parameter number I is out of range;
4805
// +1, if the answer appears to be lower than lowest search bound;
4806
// +2, if the answer appears to be higher than greatest search bound;
4807
// +3, if P + Q /= 1.
4809
// Output, double *BOUND, is only defined if STATUS is nonzero.
4810
// If STATUS is negative, then this is the value exceeded by parameter I.
4811
// if STATUS is 1 or 2, this is the search bound that was exceeded.
4814
# define tol (1.0e-8)
4815
# define atol (1.0e-50)
4816
# define zero (1.0e-300)
4818
# define maxdf 1.0e10
4821
static double K4 = 0.5e0;
4822
static double K5 = 5.0e0;
4823
static double fx,cum,ccum,pq;
4824
static unsigned long qhi,qleft,qporq;
4825
static double T2,T3,T6,T7,T8,T9,T10,T11;
4832
if(!(*which < 1 || *which > 3)) goto S30;
4833
if(!(*which < 1)) goto S10;
4842
if(*which == 1) goto S70;
4846
if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
4847
if(!(*p <= 0.0e0)) goto S40;
4857
if(*which == 1) goto S110;
4861
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4862
if(!(*q <= 0.0e0)) goto S80;
4872
if(*which == 3) goto S130;
4876
if(!(*df <= 0.0e0)) goto S120;
4882
if(*which == 1) goto S170;
4887
if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S160;
4888
if(!(pq < 0.0e0)) goto S140;
4898
if(!(*which == 1)) qporq = *p <= *q;
4900
// Select the minimum of P or Q
4901
// Calculate ANSWERS
4905
// Computing P and Q
4910
else if(2 == *which) {
4913
// .. Get initial approximation for T
4920
dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
4922
dinvr(status,t,&fx,&qleft,&qhi);
4924
if(!(*status == 1)) goto S210;
4925
cumt(t,df,&cum,&ccum);
4926
if(!qporq) goto S190;
4932
dinvr(status,t,&fx,&qleft,&qhi);
4935
if(!(*status == -1)) goto S240;
4936
if(!qleft) goto S220;
4947
else if(3 == *which) {
4956
dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
4958
dinvr(status,df,&fx,&qleft,&qhi);
4960
if(!(*status == 1)) goto S280;
4961
cumt(t,df,&cum,&ccum);
4962
if(!qporq) goto S260;
4968
dinvr(status,df,&fx,&qleft,&qhi);
4971
if(!(*status == -1)) goto S310;
4972
if(!qleft) goto S290;
4990
//****************************************************************************80
4992
void chi_noncentral_cdf_values ( int *n_data, double *x, double *lambda,
4993
int *df, double *cdf )
4995
//****************************************************************************80
4999
// CHI_NONCENTRAL_CDF_VALUES returns values of the noncentral chi CDF.
5003
// The CDF of the noncentral chi square distribution can be evaluated
5004
// within Mathematica by commands such as:
5006
// Needs["Statistics`ContinuousDistributions`"]
5007
// CDF [ NoncentralChiSquareDistribution [ DF, LAMBDA ], X ]
5020
// The Mathematica Book,
5022
// Wolfram Media / Cambridge University Press, 1999.
5026
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
5027
// first call. On each call, the routine increments N_DATA by 1, and
5028
// returns the corresponding data; when there is no more data, the
5029
// output value of N_DATA will be 0 again.
5031
// Output, double *X, the argument of the function.
5033
// Output, double *LAMBDA, the noncentrality parameter.
5035
// Output, int *DF, the number of degrees of freedom.
5037
// Output, double *CDF, the noncentral chi CDF.
5042
double cdf_vec[N_MAX] = {
5043
0.839944E+00, 0.695906E+00, 0.535088E+00,
5044
0.764784E+00, 0.620644E+00, 0.469167E+00,
5045
0.307088E+00, 0.220382E+00, 0.150025E+00,
5046
0.307116E-02, 0.176398E-02, 0.981679E-03,
5047
0.165175E-01, 0.202342E-03, 0.498448E-06,
5048
0.151325E-01, 0.209041E-02, 0.246502E-03,
5049
0.263684E-01, 0.185798E-01, 0.130574E-01,
5050
0.583804E-01, 0.424978E-01, 0.308214E-01,
5051
0.105788E+00, 0.794084E-01, 0.593201E-01 };
5052
int df_vec[N_MAX] = {
5062
double lambda_vec[N_MAX] = {
5063
0.5E+00, 0.5E+00, 0.5E+00,
5064
1.0E+00, 1.0E+00, 1.0E+00,
5065
5.0E+00, 5.0E+00, 5.0E+00,
5066
20.0E+00, 20.0E+00, 20.0E+00,
5067
30.0E+00, 30.0E+00, 30.0E+00,
5068
5.0E+00, 5.0E+00, 5.0E+00,
5069
2.0E+00, 3.0E+00, 4.0E+00,
5070
2.0E+00, 3.0E+00, 4.0E+00,
5071
2.0E+00, 3.0E+00, 4.0E+00 };
5072
double x_vec[N_MAX] = {
5073
3.000E+00, 3.000E+00, 3.000E+00,
5074
3.000E+00, 3.000E+00, 3.000E+00,
5075
3.000E+00, 3.000E+00, 3.000E+00,
5076
3.000E+00, 3.000E+00, 3.000E+00,
5077
60.000E+00, 60.000E+00, 60.000E+00,
5078
0.050E+00, 0.050E+00, 0.050E+00,
5079
4.000E+00, 4.000E+00, 4.000E+00,
5080
5.000E+00, 5.000E+00, 5.000E+00,
5081
6.000E+00, 6.000E+00, 6.000E+00 };
5088
*n_data = *n_data + 1;
5090
if ( N_MAX < *n_data )
5100
*x = x_vec[*n_data-1];
5101
*lambda = lambda_vec[*n_data-1];
5102
*df = df_vec[*n_data-1];
5103
*cdf = cdf_vec[*n_data-1];
5109
//****************************************************************************80
5111
void chi_square_cdf_values ( int *n_data, int *a, double *x, double *fx )
5113
//****************************************************************************80
5117
// CHI_SQUARE_CDF_VALUES returns some values of the Chi-Square CDF.
5121
// The value of CHI_CDF ( DF, X ) can be evaluated in Mathematica by
5124
// Needs["Statistics`ContinuousDistributions`"]
5125
// CDF[ChiSquareDistribution[DF], X ]
5137
// Milton Abramowitz and Irene Stegun,
5138
// Handbook of Mathematical Functions,
5139
// US Department of Commerce, 1964.
5142
// The Mathematica Book,
5144
// Wolfram Media / Cambridge University Press, 1999.
5148
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
5149
// first call. On each call, the routine increments N_DATA by 1, and
5150
// returns the corresponding data; when there is no more data, the
5151
// output value of N_DATA will be 0 again.
5153
// Output, int *A, the parameter of the function.
5155
// Output, double *X, the argument of the function.
5157
// Output, double *FX, the value of the function.
5162
int a_vec[N_MAX] = {
5169
double fx_vec[N_MAX] = {
5170
0.0796557E+00, 0.00498752E+00, 0.112463E+00, 0.00995017E+00,
5171
0.472911E+00, 0.181269E+00, 0.0597575E+00, 0.0175231E+00,
5172
0.682689E+00, 0.393469E+00, 0.198748E+00, 0.090204E+00,
5173
0.0374342E+00, 0.427593E+00, 0.608375E+00, 0.738536E+00,
5174
0.828203E+00, 0.88839E+00, 0.000172116E+00, 0.00365985E+00,
5176
double x_vec[N_MAX] = {
5177
0.01E+00, 0.01E+00, 0.02E+00, 0.02E+00,
5178
0.40E+00, 0.40E+00, 0.40E+00, 0.40E+00,
5179
1.00E+00, 1.00E+00, 1.00E+00, 1.00E+00,
5180
1.00E+00, 2.00E+00, 3.00E+00, 4.00E+00,
5181
5.00E+00, 6.00E+00, 1.00E+00, 2.00E+00,
5189
*n_data = *n_data + 1;
5191
if ( N_MAX < *n_data )
5200
*a = a_vec[*n_data-1];
5201
*x = x_vec[*n_data-1];
5202
*fx = fx_vec[*n_data-1];
5207
//****************************************************************************80
5209
void cumbet ( double *x, double *y, double *a, double *b, double *cum,
5212
//****************************************************************************80
5216
// CUMBET evaluates the cumulative incomplete beta distribution.
5220
// This routine calculates the CDF to X of the incomplete beta distribution
5221
// with parameters A and B. This is the integral from 0 to x
5222
// of (1/B(a,b))*f(t)) where f(t) = t**(a-1) * (1-t)**(b-1)
5230
// A R Didonato and Alfred Morris,
5232
// Significant Digit Computation of the Incomplete Beta Function Ratios.
5233
// ACM Transactions on Mathematical Software,
5234
// Volume 18, Number 3, September 1992, pages 360-373.
5238
// Input, double *X, the upper limit of integration.
5240
// Input, double *Y, the value of 1-X.
5242
// Input, double *A, *B, the parameters of the distribution.
5244
// Output, double *CUM, *CCUM, the values of the cumulative
5245
// density function and complementary cumulative density function.
5255
else if ( *y <= 0.0 )
5262
beta_inc ( a, b, x, y, cum, ccum, &ierr );
5266
//****************************************************************************80
5268
void cumbin ( double *s, double *xn, double *pr, double *ompr,
5269
double *cum, double *ccum )
5271
//****************************************************************************80
5275
// CUMBIN evaluates the cumulative binomial distribution.
5279
// This routine returns the probability of 0 to S successes in XN binomial
5280
// trials, each of which has a probability of success, PR.
5288
// Milton Abramowitz and Irene Stegun,
5289
// Handbook of Mathematical Functions
5290
// 1966, Formula 26.5.24.
5294
// Input, double *S, the upper limit of summation.
5296
// Input, double *XN, the number of trials.
5298
// Input, double *PR, the probability of success in one trial.
5300
// Input, double *OMPR, equals ( 1 - PR ).
5302
// Output, double *CUM, the cumulative binomial distribution.
5304
// Output, double *CCUM, the complement of the cumulative
5305
// binomial distribution.
5308
static double T1,T2;
5314
cumbet ( pr, ompr, &T1, &T2, ccum, cum );
5323
//****************************************************************************80
5325
void cumchi ( double *x, double *df, double *cum, double *ccum )
5327
//****************************************************************************80
5331
// CUMCHI evaluates the cumulative chi-square distribution.
5335
// Input, double *X, the upper limit of integration.
5337
// Input, double *DF, the degrees of freedom of the
5338
// chi-square distribution.
5340
// Output, double *CUM, the cumulative chi-square distribution.
5342
// Output, double *CCUM, the complement of the cumulative
5343
// chi-square distribution.
5351
cumgam ( &xx, &a, cum, ccum );
5354
//****************************************************************************80
5356
void cumchn ( double *x, double *df, double *pnonc, double *cum,
5359
//****************************************************************************80
5363
// CUMCHN evaluates the cumulative noncentral chi-square distribution.
5367
// Calculates the cumulative noncentral chi-square
5368
// distribution, i.e., the probability that a random variable
5369
// which follows the noncentral chi-square distribution, with
5370
// noncentrality parameter PNONC and continuous degrees of
5371
// freedom DF, is less than or equal to X.
5375
// Milton Abramowitz and Irene Stegun,
5376
// Handbook of Mathematical Functions
5377
// 1966, Formula 26.4.25.
5381
// Input, double *X, the upper limit of integration.
5383
// Input, double *DF, the number of degrees of freedom.
5385
// Input, double *PNONC, the noncentrality parameter of
5386
// the noncentral chi-square distribution.
5388
// Output, double *CUM, *CCUM, the CDF and complementary
5389
// CDF of the noncentral chi-square distribution.
5391
// Local Parameters:
5393
// Local, double EPS, the convergence criterion. The sum
5394
// stops when a term is less than EPS*SUM.
5396
// Local, int NTIRED, the maximum number of terms to be evaluated
5399
// Local, bool QCONV, is TRUE if convergence was achieved, that is,
5400
// the program did not stop on NTIRED criterion.
5403
# define dg(i) (*df+2.0e0*(double)(i))
5404
# define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps*sum)
5405
# define qtired(i) (int)((i) > ntired)
5407
static double eps = 1.0e-5;
5408
static int ntired = 1000;
5409
static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
5410
sumadj,term,wt,xnonc;
5411
static int i,icent,iterb,iterf;
5412
static double T1,T2,T3;
5414
if(!(*x <= 0.0e0)) goto S10;
5419
if(!(*pnonc <= 1.0e-10)) goto S20;
5421
// When non-centrality parameter is (essentially) zero,
5422
// use cumulative chi-square distribution
5424
cumchi(x,df,cum,ccum);
5427
xnonc = *pnonc/2.0e0;
5429
// The following code calculates the weight, chi-square, and
5430
// adjustment term for the central term in the infinite series.
5431
// The central term is the one in which the poisson weight is
5432
// greatest. The adjustment term is the amount that must
5433
// be subtracted from the chi-square to move up two degrees
5436
icent = fifidint(xnonc);
5437
if(icent == 0) icent = 1;
5440
// Calculate central weight term
5442
T1 = (double)(icent+1);
5443
lfact = gamma_log ( &T1 );
5444
lcntwt = -xnonc+(double)icent*log(xnonc)-lfact;
5445
centwt = exp(lcntwt);
5447
// Calculate central chi-square
5450
cumchi(x,&T2,&pcent,ccum);
5452
// Calculate central adjustment term
5454
dfd2 = dg(icent)/2.0e0;
5456
lfact = gamma_log ( &T3 );
5457
lcntaj = dfd2*log(chid2)-chid2-lfact;
5458
centaj = exp(lcntaj);
5461
// Sum backwards from the central term towards zero.
5462
// Quit whenever either
5463
// (1) the zero term is reached, or
5464
// (2) the term gets small relative to the sum, or
5465
// (3) More than NTIRED terms are totaled.
5474
if(qtired(iterb) || qsmall(term) || i == 0) goto S50;
5478
// Adjust chi-square for two fewer degrees of freedom.
5479
// The adjusted value ends up in PTERM.
5481
adj = adj*dfd2/chid2;
5482
sumadj = sumadj + adj;
5483
pterm = pcent+sumadj;
5485
// Adjust poisson weight for J decreased by one
5487
wt *= ((double)i/xnonc);
5496
// Now sum forward from the central term towards infinity.
5498
// (1) the term gets small relative to the sum, or
5499
// (2) More than NTIRED terms are totaled.
5501
sumadj = adj = centaj;
5506
if(qtired(iterf) || qsmall(term)) goto S80;
5509
// Update weights for next higher J
5511
wt *= (xnonc/(double)(i+1));
5513
// Calculate PTERM and add term to sum
5515
pterm = pcent-sumadj;
5519
// Update adjustment term for DF for next iteration
5523
adj = adj*chid2/dfd2;
5529
*ccum = 0.5e0+(0.5e0-*cum);
5535
//****************************************************************************80
5537
void cumf ( double *f, double *dfn, double *dfd, double *cum, double *ccum )
5539
//****************************************************************************80
5543
// CUMF evaluates the cumulative F distribution.
5547
// CUMF computes the integral from 0 to F of the F density with DFN
5548
// numerator and DFD denominator degrees of freedom.
5552
// Milton Abramowitz and Irene Stegun,
5553
// Handbook of Mathematical Functions
5554
// 1966, Formula 26.5.28.
5558
// Input, double *F, the upper limit of integration.
5560
// Input, double *DFN, *DFD, the number of degrees of
5561
// freedom for the numerator and denominator.
5563
// Output, double *CUM, *CCUM, the value of the F CDF and
5564
// the complementary F CDF.
5570
static double dsum,prod,xx,yy;
5572
static double T1,T2;
5574
if(!(*f <= 0.0e0)) goto S10;
5581
// XX is such that the incomplete beta with parameters
5582
// DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM
5584
// Calculate the smaller of XX and YY accurately
5601
beta_inc ( &T1, &T2, &xx, &yy, ccum, cum, &ierr );
5606
//****************************************************************************80
5608
void cumfnc ( double *f, double *dfn, double *dfd, double *pnonc,
5609
double *cum, double *ccum )
5611
//****************************************************************************80
5615
// CUMFNC evaluates the cumulative noncentral F distribution.
5619
// This routine computes the noncentral F distribution with DFN and DFD
5620
// degrees of freedom and noncentrality parameter PNONC.
5622
// The series is calculated backward and forward from J = LAMBDA/2
5623
// (this is the term with the largest Poisson weight) until
5624
// the convergence criterion is met.
5626
// The sum continues until a succeeding term is less than EPS
5627
// times the sum (or the sum is less than 1.0e-20). EPS is
5628
// set to 1.0e-4 in a data statement which can be changed.
5631
// The original version of this routine allowed the input values
5632
// of DFN and DFD to be negative (nonsensical) or zero (which
5633
// caused numerical overflow.) I have forced both these values
5634
// to be at least 1.
5642
// Milton Abramowitz and Irene Stegun,
5643
// Handbook of Mathematical Functions
5644
// 1966, Formula 26.5.16, 26.6.17, 26.6.18, 26.6.20.
5648
// Input, double *F, the upper limit of integration.
5650
// Input, double *DFN, *DFD, the number of degrees of freedom
5651
// in the numerator and denominator. Both DFN and DFD must be positive,
5652
// and normally would be integers. This routine requires that they
5653
// be no less than 1.
5655
// Input, double *PNONC, the noncentrality parameter.
5657
// Output, double *CUM, *CCUM, the noncentral F CDF and
5658
// complementary CDF.
5661
# define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
5665
static double eps = 1.0e-4;
5666
static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
5668
static int i,icent,ierr;
5669
static double T1,T2,T3,T4,T5,T6;
5671
if(!(*f <= 0.0e0)) goto S10;
5676
if(!(*pnonc < 1.0e-10)) goto S20;
5678
// Handle case in which the non-centrality parameter is
5679
// (essentially) zero.
5681
cumf(f,dfn,dfd,cum,ccum);
5684
xnonc = *pnonc/2.0e0;
5686
// Calculate the central term of the poisson weighting factor.
5688
icent = ( int ) xnonc;
5689
if(icent == 0) icent = 1;
5691
// Compute central weight term
5693
T1 = (double)(icent+1);
5694
centwt = exp(-xnonc+(double)icent*log(xnonc)- gamma_log ( &T1 ) );
5696
// Compute central incomplete beta term
5697
// Assure that minimum of arg to beta and 1 - arg is computed
5708
T2 = *dfn*half+(double)icent;
5710
beta_inc ( &T2, &T3, &xx, &yy, &betdn, &dummy, &ierr );
5711
adn = *dfn/2.0e0+(double)icent;
5717
// Now sum terms backward from icent until convergence or all done
5723
dnterm = exp( gamma_log ( &T4 ) - gamma_log ( &T5 )
5724
- gamma_log ( &b ) + adn * log ( xx ) + b * log(yy));
5726
if(qsmall(xmult*betdn) || i <= 0) goto S40;
5727
xmult *= ((double)i/xnonc);
5730
dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
5732
sum += (xmult*betdn);
5737
// Now sum forwards until convergence
5740
if(aup-1.0+b == 0) upterm = exp(-gamma_log ( &aup )
5741
- gamma_log ( &b ) + (aup-1.0)*log(xx)+
5745
upterm = exp( gamma_log ( &T6 ) - gamma_log ( &aup )
5746
- gamma_log ( &b ) + (aup-1.0)*log(xx)+b*
5751
if(qsmall(xmult*betup)) goto S70;
5753
xmult *= (xnonc/(double)i);
5756
upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
5758
sum += (xmult*betup);
5762
*ccum = 0.5e0+(0.5e0-*cum);
5768
//****************************************************************************80
5770
void cumgam ( double *x, double *a, double *cum, double *ccum )
5772
//****************************************************************************80
5776
// CUMGAM evaluates the cumulative incomplete gamma distribution.
5780
// This routine computes the cumulative distribution function of the
5781
// incomplete gamma distribution, i.e., the integral from 0 to X of
5783
// (1/GAM(A))*EXP(-T)*T**(A-1) DT
5785
// where GAM(A) is the complete gamma function of A, i.e.,
5787
// GAM(A) = integral from 0 to infinity of EXP(-T)*T**(A-1) DT
5791
// Input, double *X, the upper limit of integration.
5793
// Input, double *A, the shape parameter of the incomplete
5794
// Gamma distribution.
5796
// Output, double *CUM, *CCUM, the incomplete Gamma CDF and
5797
// complementary CDF.
5802
if(!(*x <= 0.0e0)) goto S10;
5807
gamma_inc ( a, x, cum, ccum, &K1 );
5809
// Call gratio routine
5813
//****************************************************************************80
5815
void cumnbn ( double *s, double *xn, double *pr, double *ompr,
5816
double *cum, double *ccum )
5818
//****************************************************************************80
5822
// CUMNBN evaluates the cumulative negative binomial distribution.
5826
// This routine returns the probability that there will be F or
5827
// fewer failures before there are S successes, with each binomial
5828
// trial having a probability of success PR.
5830
// Prob(# failures = F | S successes, PR) =
5832
// ( ) * PR^S * (1-PR)^F
5837
// Milton Abramowitz and Irene Stegun,
5838
// Handbook of Mathematical Functions
5839
// 1966, Formula 26.5.26.
5843
// Input, double *F, the number of failures.
5845
// Input, double *S, the number of successes.
5847
// Input, double *PR, *OMPR, the probability of success on
5848
// each binomial trial, and the value of (1-PR).
5850
// Output, double *CUM, *CCUM, the negative binomial CDF,
5851
// and the complementary CDF.
5857
cumbet(pr,ompr,xn,&T1,cum,ccum);
5860
//****************************************************************************80
5862
void cumnor ( double *arg, double *result, double *ccum )
5864
//****************************************************************************80
5868
// CUMNOR computes the cumulative normal distribution.
5872
// This function evaluates the normal distribution function:
5876
// P(x) = ----------- | e dt
5880
// This transportable program uses rational functions that
5881
// theoretically approximate the normal distribution function to
5882
// at least 18 significant decimal digits. The accuracy achieved
5883
// depends on the arithmetic system, the compiler, the intrinsic
5884
// functions, and proper selection of the machine-dependent
5890
// Mathematics and Computer Science Division
5891
// Argonne National Laboratory
5892
// Argonne, IL 60439
5897
// Rational Chebyshev approximations for the error function,
5898
// Mathematics of Computation,
5899
// 1969, pages 631-637.
5903
// SPECFUN - A Portable FORTRAN Package of Special Function Routines
5904
// and Test Drivers,
5905
// ACM Transactions on Mathematical Software,
5906
// Volume 19, 1993, pages 22-32.
5910
// Input, double *ARG, the upper limit of integration.
5912
// Output, double *CUM, *CCUM, the Normal density CDF and
5913
// complementary CDF.
5915
// Local Parameters:
5917
// Local, double EPS, the argument below which anorm(x)
5918
// may be represented by 0.5D+00 and above which x*x will not underflow.
5919
// A conservative value is the largest machine number X
5920
// such that 1.0D+00 + X = 1.0D+00 to machine precision.
5923
static double a[5] = {
5924
2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
5925
1.8154981253343561249e04,6.5682337918207449113e-2
5927
static double b[4] = {
5928
4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
5929
4.5507789335026729956e04
5931
static double c[9] = {
5932
3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
5933
5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
5934
1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
5936
static double d[8] = {
5937
2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
5938
6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
5939
3.8912003286093271411e04,1.9685429676859990727e04
5941
static double half = 0.5e0;
5942
static double p[6] = {
5943
2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
5944
1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
5946
static double one = 1.0e0;
5947
static double q[5] = {
5948
1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
5949
3.78239633202758244e-3,7.29751555083966205e-5
5951
static double sixten = 1.60e0;
5952
static double sqrpi = 3.9894228040143267794e-1;
5953
static double thrsh = 0.66291e0;
5954
static double root32 = 5.656854248e0;
5955
static double zero = 0.0e0;
5959
static double del,eps,temp,x,xden,xnum,y,xsq,min;
5961
// Machine dependent constants
5963
eps = dpmpar(&K1)*0.5e0;
5969
// Evaluate anorm for |X| <= 0.66291
5972
if(y > eps) xsq = x*x;
5975
for ( i = 0; i < 3; i++ )
5977
xnum = (xnum+a[i])*xsq;
5978
xden = (xden+b[i])*xsq;
5980
*result = x*(xnum+a[3])/(xden+b[3]);
5982
*result = half+temp;
5986
// Evaluate anorm for 0.66291 <= |X| <= sqrt(32)
5988
else if(y <= root32) {
5991
for ( i = 0; i < 7; i++ )
5993
xnum = (xnum+c[i])*y;
5994
xden = (xden+d[i])*y;
5996
*result = (xnum+c[7])/(xden+d[7]);
5997
xsq = fifdint(y*sixten)/sixten;
5998
del = (y-xsq)*(y+xsq);
5999
*result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
6000
*ccum = one-*result;
6008
// Evaluate anorm for |X| > sqrt(32)
6015
for ( i = 0; i < 4; i++ )
6017
xnum = (xnum+p[i])*xsq;
6018
xden = (xden+q[i])*xsq;
6020
*result = xsq*(xnum+p[4])/(xden+q[4]);
6021
*result = (sqrpi-*result)/y;
6022
xsq = fifdint(x*sixten)/sixten;
6023
del = (x-xsq)*(x+xsq);
6024
*result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
6025
*ccum = one-*result;
6032
if(*result < min) *result = 0.0e0;
6034
// Fix up for negative argument, erf, etc.
6036
if(*ccum < min) *ccum = 0.0e0;
6038
//****************************************************************************80
6040
void cumpoi ( double *s, double *xlam, double *cum, double *ccum )
6042
//****************************************************************************80
6046
// CUMPOI evaluates the cumulative Poisson distribution.
6050
// CUMPOI returns the probability of S or fewer events in a Poisson
6051
// distribution with mean XLAM.
6055
// Milton Abramowitz and Irene Stegun,
6056
// Handbook of Mathematical Functions,
6061
// Input, double *S, the upper limit of cumulation of the
6062
// Poisson density function.
6064
// Input, double *XLAM, the mean of the Poisson distribution.
6066
// Output, double *CUM, *CCUM, the Poisson density CDF and
6067
// complementary CDF.
6070
static double chi,df;
6072
df = 2.0e0*(*s+1.0e0);
6074
cumchi(&chi,&df,ccum,cum);
6077
//****************************************************************************80
6079
void cumt ( double *t, double *df, double *cum, double *ccum )
6081
//****************************************************************************80
6085
// CUMT evaluates the cumulative T distribution.
6089
// Milton Abramowitz and Irene Stegun,
6090
// Handbook of Mathematical Functions,
6095
// Input, double *T, the upper limit of integration.
6097
// Input, double *DF, the number of degrees of freedom of
6098
// the T distribution.
6100
// Output, double *CUM, *CCUM, the T distribution CDF and
6101
// complementary CDF.
6105
static double dfptt;
6106
static double K2 = 0.5e0;
6114
dfptt = ( *df ) + tt;
6117
T1 = 0.5e0 * ( *df );
6118
cumbet ( &xx, &yy, &T1, &K2, &a, &oma );
6123
*ccum = oma + ( *cum );
6128
*cum = oma + ( *ccum );
6132
//****************************************************************************80
6134
double dbetrm ( double *a, double *b )
6136
//****************************************************************************80
6140
// DBETRM computes the Sterling remainder for the complete beta function.
6144
// Log(Beta(A,B)) = Lgamma(A) + Lgamma(B) - Lgamma(A+B)
6145
// where Lgamma is the log of the (complete) gamma function
6147
// Let ZZ be approximation obtained if each log gamma is approximated
6148
// by Sterling's formula, i.e.,
6149
// Sterling(Z) = LOG( SQRT( 2*PI ) ) + ( Z-0.5D+00 ) * LOG( Z ) - Z
6151
// The Sterling remainder is Log(Beta(A,B)) - ZZ.
6155
// Input, double *A, *B, the parameters of the Beta function.
6157
// Output, double DBETRM, the Sterling remainder.
6160
static double dbetrm,T1,T2,T3;
6162
// Try to sum from smallest to largest
6165
dbetrm = -dstrem(&T1);
6166
T2 = fifdmax1(*a,*b);
6167
dbetrm += dstrem(&T2);
6168
T3 = fifdmin1(*a,*b);
6169
dbetrm += dstrem(&T3);
6172
//****************************************************************************80
6174
double dexpm1 ( double *x )
6176
//****************************************************************************80
6180
// DEXPM1 evaluates the function EXP(X) - 1.
6184
// Armido DiDinato and Alfred Morris,
6186
// Significant Digit Computation of the Incomplete Beta Function Ratios,
6187
// ACM Transactions on Mathematical Software,
6188
// Volume 18, 1993, pages 360-373.
6192
// Input, double *X, the value at which exp(X)-1 is desired.
6194
// Output, double DEXPM1, the value of exp(X)-1.
6197
static double p1 = .914041914819518e-09;
6198
static double p2 = .238082361044469e-01;
6199
static double q1 = -.499999999085958e+00;
6200
static double q2 = .107141568980644e+00;
6201
static double q3 = -.119041179760821e-01;
6202
static double q4 = .595130811860248e-03;
6203
static double dexpm1;
6206
if ( fabs(*x) <= 0.15e0 )
6219
else if ( *x <= 0.0e0 )
6222
dexpm1 = w-0.5e0-0.5e0;
6227
dexpm1 = w*(0.5e0+(0.5e0-1.0e0/w));
6232
//****************************************************************************80
6234
double dinvnr ( double *p, double *q )
6236
//****************************************************************************80
6240
// DINVNR computes the inverse of the normal distribution.
6244
// Returns X such that CUMNOR(X) = P, i.e., the integral from -
6245
// infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
6247
// The rational function on page 95 of Kennedy and Gentle is used as a start
6248
// value for the Newton method of finding roots.
6252
// Kennedy and Gentle,
6253
// Statistical Computing,
6254
// Marcel Dekker, NY, 1980,
6259
// Input, double *P, *Q, the probability, and the complementary
6262
// Output, double DINVNR, the argument X for which the
6263
// Normal CDF has the value P.
6267
# define eps (1.0e-13)
6268
# define r2pi 0.3989422804014326e0
6269
# define nhalf (-0.5e0)
6270
# define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
6272
static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
6274
static unsigned long qporq;
6277
// FIND MINIMUM OF P AND Q
6280
if(!qporq) goto S10;
6287
// INITIALIZATION STEP
6289
strtx = stvaln(&pp);
6292
// NEWTON INTERATIONS
6294
for ( i = 1; i <= maxit; i++ )
6296
cumnor(&xcur,&cum,&ccum);
6297
dx = (cum-pp)/dennor(xcur);
6299
if(fabs(dx/xcur) < eps) goto S40;
6303
// IF WE GET HERE, NEWTON HAS FAILED
6305
if(!qporq) dinvnr = -dinvnr;
6309
// IF WE GET HERE, NEWTON HAS SUCCEDED
6312
if(!qporq) dinvnr = -dinvnr;
6320
//****************************************************************************80
6322
void dinvr ( int *status, double *x, double *fx,
6323
unsigned long *qleft, unsigned long *qhi )
6325
//****************************************************************************80
6329
// DINVR bounds the zero of the function and invokes DZROR.
6333
// This routine seeks to find bounds on a root of the function and
6334
// invokes ZROR to perform the zero finding. STINVR must have been
6335
// called before this routine in order to set its parameters.
6339
// J C P Bus and T J Dekker,
6340
// Two Efficient Algorithms with Guaranteed Convergence for
6341
// Finding a Zero of a Function,
6342
// ACM Transactions on Mathematical Software,
6343
// Volume 1, Number 4, pages 330-345, 1975.
6347
// Input/output, integer STATUS. At the beginning of a zero finding
6348
// problem, STATUS should be set to 0 and INVR invoked. The value
6349
// of parameters other than X will be ignored on this call.
6350
// If INVR needs the function to be evaluated, it will set STATUS to 1
6351
// and return. The value of the function should be set in FX and INVR
6352
// again called without changing any of its other parameters.
6353
// If INVR finishes without error, it returns with STATUS 0, and X an
6354
// approximate root of F(X).
6355
// If INVR cannot bound the function, it returns a negative STATUS and
6356
// sets QLEFT and QHI.
6358
// Output, double precision X, the value at which F(X) is to be evaluated.
6360
// Input, double precision FX, the value of F(X) calculated by the user
6361
// on the previous call, when INVR returned with STATUS = 1.
6363
// Output, logical QLEFT, is defined only if QMFINV returns FALSE. In that
6364
// case, QLEFT is TRUE if the stepping search terminated unsucessfully
6365
// at SMALL, and FALSE if the search terminated unsucessfully at BIG.
6367
// Output, logical QHI, is defined only if QMFINV returns FALSE. In that
6368
// case, it is TRUE if Y < F(X) at the termination of the search and FALSE
6372
E0000(0,status,x,fx,qleft,qhi,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
6374
//****************************************************************************80
6376
double dlanor ( double *x )
6378
//****************************************************************************80
6382
// DLANOR evaluates the logarithm of the asymptotic Normal CDF.
6386
// This routine computes the logarithm of the cumulative normal distribution
6387
// from abs ( x ) to infinity for 5 <= abs ( X ).
6389
// The relative error at X = 5 is about 0.5D-5.
6393
// Milton Abramowitz and Irene Stegun,
6394
// Handbook of Mathematical Functions
6395
// 1966, Formula 26.2.12.
6399
// Input, double *X, the value at which the Normal CDF is to be
6400
// evaluated. It is assumed that 5 <= abs ( X ).
6402
// Output, double DLANOR, the logarithm of the asymptotic
6406
# define dlsqpi 0.91893853320467274177e0
6408
static double coef[12] = {
6409
-1.0e0,3.0e0,-15.0e0,105.0e0,-945.0e0,10395.0e0,-135135.0e0,2027025.0e0,
6410
-34459425.0e0,654729075.0e0,-13749310575.e0,316234143225.0e0
6413
static double dlanor,approx,correc,xx,xx2,T2;
6418
ftnstop(" Argument too small in DLANOR");
6420
approx = -dlsqpi-0.5e0*xx*xx-log(xx);
6423
correc = eval_pol ( coef, &K1, &T2 ) / xx2;
6424
correc = alnrel ( &correc );
6425
dlanor = approx+correc;
6429
//****************************************************************************80
6431
double dpmpar ( int *i )
6433
//****************************************************************************80
6437
// DPMPAR provides machine constants for double precision arithmetic.
6441
// DPMPAR PROVIDES THE double PRECISION MACHINE CONSTANTS FOR
6442
// THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
6443
// I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
6444
// double PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
6445
// ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
6447
// DPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
6449
// DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
6451
// DPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
6454
// ALFRED H. MORRIS, JR.
6455
// NAVAL SURFACE WARFARE CENTER
6456
// DAHLGREN VIRGINIA
6458
// MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE
6459
// CONSTANTS FOR THE COMPUTER BEING USED. THIS MODIFICATION WAS
6460
// MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION
6467
static double value,b,binv,bm1,one,w,z;
6468
static int emax,emin,ibeta,m;
6470
if(*i > 1) goto S10;
6473
value = pow(b,(double)(1-m));
6476
if(*i > 2) goto S20;
6481
w = pow(b,(double)(emin+2));
6482
value = w*binv*binv*binv;
6485
ibeta = ipmpar(&K1);
6491
z = pow(b,(double)(m-1));
6492
w = ((z-one)*b+bm1)/(b*z);
6493
z = pow(b,(double)(emax-2));
6497
//****************************************************************************80
6499
void dstinv ( double *zsmall, double *zbig, double *zabsst,
6500
double *zrelst, double *zstpmu, double *zabsto, double *zrelto )
6502
//****************************************************************************80
6506
// DSTINV seeks a value X such that F(X) = Y.
6510
// Double Precision - SeT INverse finder - Reverse Communication
6512
// Concise Description - Given a monotone function F finds X
6513
// such that F(X) = Y. Uses Reverse communication -- see invr.
6514
// This routine sets quantities needed by INVR.
6515
// More Precise Description of INVR -
6516
// F must be a monotone function, the results of QMFINV are
6517
// otherwise undefined. QINCR must be .TRUE. if F is non-
6518
// decreasing and .FALSE. if F is non-increasing.
6519
// QMFINV will return .TRUE. if and only if F(SMALL) and
6520
// F(BIG) bracket Y, i. e.,
6521
// QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
6522
// QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
6523
// if QMFINV returns .TRUE., then the X returned satisfies
6524
// the following condition. let
6525
// TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
6526
// then if QINCR is .TRUE.,
6527
// F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
6528
// and if QINCR is .FALSE.
6529
// F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
6531
// SMALL --> The left endpoint of the interval to be
6532
// searched for a solution.
6533
// SMALL is DOUBLE PRECISION
6534
// BIG --> The right endpoint of the interval to be
6535
// searched for a solution.
6536
// BIG is DOUBLE PRECISION
6537
// ABSSTP, RELSTP --> The initial step size in the search
6538
// is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
6539
// ABSSTP is DOUBLE PRECISION
6540
// RELSTP is DOUBLE PRECISION
6541
// STPMUL --> When a step doesn't bound the zero, the step
6542
// size is multiplied by STPMUL and another step
6543
// taken. A popular value is 2.0
6544
// DOUBLE PRECISION STPMUL
6545
// ABSTOL, RELTOL --> Two numbers that determine the accuracy
6546
// of the solution. See function for a precise definition.
6547
// ABSTOL is DOUBLE PRECISION
6548
// RELTOL is DOUBLE PRECISION
6550
// Compares F(X) with Y for the input value of X then uses QINCR
6551
// to determine whether to step left or right to bound the
6552
// desired x. the initial step size is
6553
// MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
6554
// Iteratively steps right or left until it bounds X.
6555
// At each step which doesn't bound X, the step size is doubled.
6556
// The routine is careful never to step beyond SMALL or BIG. If
6557
// it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
6558
// after setting QLEFT and QHI.
6559
// If X is successfully bounded then Algorithm R of the paper
6560
// 'Two Efficient Algorithms with Guaranteed Convergence for
6561
// Finding a Zero of a Function' by J. C. P. Bus and
6562
// T. J. Dekker in ACM Transactions on Mathematical
6563
// Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
6564
// to find the zero of the function F(X)-Y. This is routine
6568
E0000(1,NULL,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
6571
//****************************************************************************80
6573
double dstrem ( double *z )
6575
//****************************************************************************80
6579
// DSTREM computes the Sterling remainder ln ( Gamma ( Z ) ) - Sterling ( Z ).
6583
// This routine returns
6585
// ln ( Gamma ( Z ) ) - Sterling ( Z )
6587
// where Sterling(Z) is Sterling's approximation to ln ( Gamma ( Z ) ).
6589
// Sterling(Z) = ln ( sqrt ( 2 * PI ) ) + ( Z - 0.5 ) * ln ( Z ) - Z
6591
// If 6 <= Z, the routine uses 9 terms of a series in Bernoulli numbers,
6592
// with values calculated using Maple.
6594
// Otherwise, the difference is computed explicitly.
6602
// Input, double *Z, the value at which the Sterling
6603
// remainder is to be calculated. Z must be positive.
6605
// Output, double DSTREM, the Sterling remainder.
6608
# define hln2pi 0.91893853320467274178e0
6611
static double coef[ncoef] = {
6612
0.0e0,0.0833333333333333333333333333333e0,
6613
-0.00277777777777777777777777777778e0,0.000793650793650793650793650793651e0,
6614
-0.000595238095238095238095238095238e0,
6615
0.000841750841750841750841750841751e0,-0.00191752691752691752691752691753e0,
6616
0.00641025641025641025641025641026e0,-0.0295506535947712418300653594771e0,
6617
0.179644372368830573164938490016e0
6620
static double dstrem,sterl,T2;
6622
// For information, here are the next 11 coefficients of the
6623
// remainder term in Sterling's formula
6624
// -1.39243221690590111642743221691
6625
// 13.4028640441683919944789510007
6626
// -156.848284626002017306365132452
6627
// 2193.10333333333333333333333333
6628
// -36108.7712537249893571732652192
6629
// 691472.268851313067108395250776
6630
// -0.152382215394074161922833649589D8
6631
// 0.382900751391414141414141414141D9
6632
// -0.108822660357843910890151491655D11
6633
// 0.347320283765002252252252252252D12
6634
// -0.123696021422692744542517103493D14
6638
ftnstop ( "Zero or negative argument in DSTREM" );
6640
if(!(*z > 6.0e0)) goto S10;
6641
T2 = 1.0e0/pow(*z,2.0);
6642
dstrem = eval_pol ( coef, &K1, &T2 )**z;
6645
sterl = hln2pi+(*z-0.5e0)*log(*z)-*z;
6646
dstrem = gamma_log ( z ) - sterl;
6652
//****************************************************************************80
6654
void dstzr ( double *zxlo, double *zxhi, double *zabstl, double *zreltl )
6656
//****************************************************************************80
6660
// DSTXR sets quantities needed by the zero finder.
6664
// Double precision SeT ZeRo finder - Reverse communication version
6666
// Sets quantities needed by ZROR. The function of ZROR
6667
// and the quantities set is given here.
6668
// Concise Description - Given a function F
6669
// find XLO such that F(XLO) = 0.
6670
// More Precise Description -
6671
// Input condition. F is a double function of a single
6672
// double argument and XLO and XHI are such that
6673
// F(XLO)*F(XHI) .LE. 0.0
6674
// If the input condition is met, QRZERO returns .TRUE.
6675
// and output values of XLO and XHI satisfy the following
6676
// F(XLO)*F(XHI) .LE. 0.
6677
// ABS(F(XLO) .LE. ABS(F(XHI)
6678
// ABS(XLO-XHI) .LE. TOL(X)
6680
// TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
6681
// If this algorithm does not find XLO and XHI satisfying
6682
// these conditions then QRZERO returns .FALSE. This
6683
// implies that the input condition was not met.
6685
// XLO --> The left endpoint of the interval to be
6686
// searched for a solution.
6687
// XLO is DOUBLE PRECISION
6688
// XHI --> The right endpoint of the interval to be
6690
// XHI is DOUBLE PRECISION
6691
// ABSTOL, RELTOL --> Two numbers that determine the accuracy
6692
// of the solution. See function for a
6693
// precise definition.
6694
// ABSTOL is DOUBLE PRECISION
6695
// RELTOL is DOUBLE PRECISION
6697
// Algorithm R of the paper 'Two Efficient Algorithms with
6698
// Guaranteed Convergence for Finding a Zero of a Function'
6699
// by J. C. P. Bus and T. J. Dekker in ACM Transactions on
6700
// Mathematical Software, Volume 1, no. 4 page 330
6701
// (Dec. '75) is employed to find the zero of F(X)-Y.
6704
E0001(1,NULL,NULL,NULL,NULL,NULL,NULL,NULL,zabstl,zreltl,zxhi,zxlo);
6706
//****************************************************************************80
6708
double dt1 ( double *p, double *q, double *df )
6710
//****************************************************************************80
6714
// DT1 computes an approximate inverse of the cumulative T distribution.
6718
// Returns the inverse of the T distribution function, i.e.,
6719
// the integral from 0 to INVT of the T density is P. This is an
6720
// initial approximation.
6724
// Input, double *P, *Q, the value whose inverse from the
6725
// T distribution CDF is desired, and the value (1-P).
6727
// Input, double *DF, the number of degrees of freedom of the
6730
// Output, double DT1, the approximate value of X for which
6731
// the T density CDF with DF degrees of freedom has value P.
6734
static double coef[4][5] = {
6735
1.0e0,1.0e0,0.0e0,0.0e0,0.0e0,3.0e0,16.0e0,5.0e0,0.0e0,0.0e0,-15.0e0,17.0e0,
6736
19.0e0,3.0e0,0.0e0,-945.0e0,-1920.0e0,1482.0e0,776.0e0,79.0e0
6738
static double denom[4] = {
6739
4.0e0,96.0e0,384.0e0,92160.0e0
6741
static int ideg[4] = {
6744
static double dt1,denpow,sum,term,x,xp,xx;
6747
x = fabs(dinvnr(p,q));
6751
for ( i = 0; i < 4; i++ )
6753
term = eval_pol ( &coef[i][0], &ideg[i], &xx ) * x;
6755
sum += (term/(denpow*denom[i]));
6757
if(!(*p >= 0.5e0)) goto S20;
6766
//****************************************************************************80
6768
void dzror ( int *status, double *x, double *fx, double *xlo,
6769
double *xhi, unsigned long *qleft, unsigned long *qhi )
6771
//****************************************************************************80
6775
// DZROR seeks the zero of a function using reverse communication.
6779
// Performs the zero finding. STZROR must have been called before
6780
// this routine in order to set its parameters.
6786
// STATUS <--> At the beginning of a zero finding problem, STATUS
6787
// should be set to 0 and ZROR invoked. (The value
6788
// of other parameters will be ignored on this call.)
6790
// When ZROR needs the function evaluated, it will set
6791
// STATUS to 1 and return. The value of the function
6792
// should be set in FX and ZROR again called without
6793
// changing any of its other parameters.
6795
// When ZROR has finished without error, it will return
6796
// with STATUS 0. In that case (XLO,XHI) bound the answe
6798
// If ZROR finds an error (which implies that F(XLO)-Y an
6799
// F(XHI)-Y have the same sign, it returns STATUS -1. In
6800
// this case, XLO and XHI are undefined.
6803
// X <-- The value of X at which F(X) is to be evaluated.
6804
// DOUBLE PRECISION X
6806
// FX --> The value of F(X) calculated when ZROR returns with
6808
// DOUBLE PRECISION FX
6810
// XLO <-- When ZROR returns with STATUS = 0, XLO bounds the
6811
// inverval in X containing the solution below.
6812
// DOUBLE PRECISION XLO
6814
// XHI <-- When ZROR returns with STATUS = 0, XHI bounds the
6815
// inverval in X containing the solution above.
6816
// DOUBLE PRECISION XHI
6818
// QLEFT <-- .TRUE. if the stepping search terminated unsucessfully
6819
// at XLO. If it is .FALSE. the search terminated
6820
// unsucessfully at XHI.
6823
// QHI <-- .TRUE. if F(X) .GT. Y at the termination of the
6824
// search and .FALSE. if F(X) .LT. Y at the
6825
// termination of the search.
6830
E0001(0,status,x,fx,xlo,xhi,qleft,qhi,NULL,NULL,NULL,NULL);
6832
//****************************************************************************80
6834
static void E0000 ( int IENTRY, int *status, double *x, double *fx,
6835
unsigned long *qleft, unsigned long *qhi, double *zabsst,
6836
double *zabsto, double *zbig, double *zrelst,
6837
double *zrelto, double *zsmall, double *zstpmu )
6839
//****************************************************************************80
6843
// E0000 is a reverse-communication zero bounder.
6846
# define qxmon(zx,zy,zz) (int)((zx) <= (zy) && (zy) <= (zz))
6848
static double absstp;
6849
static double abstol;
6850
static double big,fbig,fsmall,relstp,reltol,small,step,stpmul,xhi,
6851
xlb,xlo,xsave,xub,yy;
6853
static unsigned long qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup;
6854
switch(IENTRY){case 0: goto DINVR; case 1: goto DSTINV;}
6856
if(*status > 0) goto S310;
6857
qcond = !qxmon(small,*x,big);
6860
ftnstop(" SMALL, X, BIG not monotone in INVR");
6864
// See that SMALL and BIG bound the zero and set QINCR
6868
// GET-FUNCTION-VALUE
6876
// GET-FUNCTION-VALUE
6882
qincr = fbig > fsmall;
6883
if(!qincr) goto S50;
6884
if(fsmall <= 0.0e0) goto S30;
6889
if(fbig >= 0.0e0) goto S40;
6896
if(fsmall >= 0.0e0) goto S60;
6902
if(fbig <= 0.0e0) goto S70;
6910
step = fifdmax1(absstp,relstp*fabs(*x));
6913
// GET-FUNCTION-VALUE
6919
if(!(yy == 0.0e0)) goto S100;
6924
qup = qincr && yy < 0.0e0 || !qincr && yy > 0.0e0;
6926
// HANDLE CASE IN WHICH WE MUST STEP HIGHER
6930
xub = fifdmin1(xlb+step,big);
6933
if(qcond) goto S150;
6940
// GET-FUNCTION-VALUE
6946
qbdd = qincr && yy >= 0.0e0 || !qincr && yy <= 0.0e0;
6948
qcond = qbdd || qlim;
6949
if(qcond) goto S140;
6952
xub = fifdmin1(xlb+step,big);
6956
if(!(qlim && !qbdd)) goto S160;
6966
// HANDLE CASE IN WHICH WE MUST STEP LOWER
6969
xlb = fifdmax1(xub-step,small);
6972
if(qcond) goto S220;
6979
// GET-FUNCTION-VALUE
6985
qbdd = qincr && yy <= 0.0e0 || !qincr && yy >= 0.0e0;
6986
qlim = xlb <= small;
6987
qcond = qbdd || qlim;
6988
if(qcond) goto S210;
6991
xlb = fifdmax1(xub-step,small);
6995
if(!(qlim && !qbdd)) goto S230;
7003
dstzr(&xlb,&xub,&abstol,&reltol);
7005
// IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
7010
if(!(*status == 1)) goto S290;
7012
dzror ( status, x, fx, &xlo, &xhi, &qdum1, &qdum2 );
7013
if(!(*status == 1)) goto S280;
7015
// GET-FUNCTION-VALUE
7037
// TO GET-FUNCTION-VALUE
7042
switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S90;case
7043
4: goto S130;case 5: goto S200;case 6: goto S270;default: break;}
7046
//****************************************************************************80
7048
static void E0001 ( int IENTRY, int *status, double *x, double *fx,
7049
double *xlo, double *xhi, unsigned long *qleft,
7050
unsigned long *qhi, double *zabstl, double *zreltl,
7051
double *zxhi, double *zxlo )
7053
//****************************************************************************80
7057
// E00001 is a reverse-communication zero finder.
7060
# define ftol(zx) (0.5e0*fifdmax1(abstol,reltol*fabs((zx))))
7062
static double a,abstol,b,c,d,fa,fb,fc,fd,fda;
7063
static double fdb,m,mb,p,q,reltol,tol,w,xxhi,xxlo;
7064
static int ext,i99999;
7065
static unsigned long first,qrzero;
7066
switch(IENTRY){case 0: goto DZROR; case 1: goto DSTZR;}
7068
if(*status > 0) goto S280;
7073
// GET-FUNCTION-VALUE
7082
// GET-FUNCTION-VALUE
7088
// Check that F(ZXLO) < 0 < F(ZXHI) or
7089
// F(ZXLO) > 0 > F(ZXHI)
7091
if(!(fb < 0.0e0)) goto S40;
7092
if(!(*fx < 0.0e0)) goto S30;
7099
if(!(fb > 0.0e0)) goto S60;
7100
if(!(*fx > 0.0e0)) goto S50;
7114
if(!(fabs(fc) < fabs(fb))) goto S100;
7115
if(!(c != a)) goto S90;
7130
if(!(fabs(mb) > tol)) goto S240;
7131
if(!(ext > 3)) goto S110;
7135
tol = fifdsign(tol,mb);
7137
if(!first) goto S120;
7142
fdb = (fd-fb)/(d-b);
7143
fda = (fd-fa)/(d-a);
7147
if(!(p < 0.0e0)) goto S140;
7151
if(ext == 3) p *= 2.0e0;
7152
if(!(p*1.0e0 == 0.0e0 || p <= q*tol)) goto S150;
7156
if(!(p < mb*q)) goto S160;
7172
// GET-FUNCTION-VALUE
7178
if(!(fc*fb >= 0.0e0)) goto S210;
7181
if(!(w == mb)) goto S220;
7190
qrzero = fc >= 0.0e0 && fb <= 0.0e0 || fc < 0.0e0 && fb >= 0.0e0;
7191
if(!qrzero) goto S250;
7206
// TO GET-FUNCTION-VALUE
7211
switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S200;
7215
//****************************************************************************80
7217
void erf_values ( int *n_data, double *x, double *fx )
7219
//****************************************************************************80
7223
// ERF_VALUES returns some values of the ERF or "error" function.
7227
// ERF(X) = ( 2 / sqrt ( PI ) * integral ( 0 <= T <= X ) exp ( - T^2 ) dT
7239
// Milton Abramowitz and Irene Stegun,
7240
// Handbook of Mathematical Functions,
7241
// US Department of Commerce, 1964.
7245
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
7246
// first call. On each call, the routine increments N_DATA by 1, and
7247
// returns the corresponding data; when there is no more data, the
7248
// output value of N_DATA will be 0 again.
7250
// Output, double *X, the argument of the function.
7252
// Output, double *FX, the value of the function.
7257
double fx_vec[N_MAX] = {
7258
0.0000000000E+00, 0.1124629160E+00, 0.2227025892E+00, 0.3286267595E+00,
7259
0.4283923550E+00, 0.5204998778E+00, 0.6038560908E+00, 0.6778011938E+00,
7260
0.7421009647E+00, 0.7969082124E+00, 0.8427007929E+00, 0.8802050696E+00,
7261
0.9103139782E+00, 0.9340079449E+00, 0.9522851198E+00, 0.9661051465E+00,
7262
0.9763483833E+00, 0.9837904586E+00, 0.9890905016E+00, 0.9927904292E+00,
7264
double x_vec[N_MAX] = {
7265
0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00,
7266
0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00,
7267
0.8E+00, 0.9E+00, 1.0E+00, 1.1E+00,
7268
1.2E+00, 1.3E+00, 1.4E+00, 1.5E+00,
7269
1.6E+00, 1.7E+00, 1.8E+00, 1.9E+00,
7277
*n_data = *n_data + 1;
7279
if ( N_MAX < *n_data )
7287
*x = x_vec[*n_data-1];
7288
*fx = fx_vec[*n_data-1];
7293
//****************************************************************************80
7295
double error_f ( double *x )
7297
//****************************************************************************80
7301
// ERROR_F evaluates the error function ERF.
7305
// Input, double *X, the argument.
7307
// Output, double ERROR_F, the value of the error function at X.
7310
static double c = .564189583547756e0;
7311
static double a[5] = {
7312
.771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
7313
.479137145607681e-01,.128379167095513e+00
7315
static double b[3] = {
7316
.301048631703895e-02,.538971687740286e-01,.375795757275549e+00
7318
static double p[8] = {
7319
-1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
7320
4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
7321
4.51918953711873e+02,3.00459261020162e+02
7323
static double q[8] = {
7324
1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
7325
2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7326
7.90950925327898e+02,3.00459260956983e+02
7328
static double r[5] = {
7329
2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
7330
4.65807828718470e+00,2.82094791773523e-01
7332
static double s[4] = {
7333
9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
7334
1.80124575948747e+01
7336
static double erf1,ax,bot,t,top,x2;
7339
if(ax > 0.5e0) goto S10;
7341
top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
7342
bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
7343
erf1 = *x*(top/bot);
7346
if(ax > 4.0e0) goto S20;
7347
top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7349
bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7351
erf1 = 0.5e0+(0.5e0-exp(-(*x**x))*top/bot);
7352
if(*x < 0.0e0) erf1 = -erf1;
7355
if(ax >= 5.8e0) goto S30;
7358
top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
7359
bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
7360
erf1 = (c-top/(x2*bot))/ax;
7361
erf1 = 0.5e0+(0.5e0-exp(-x2)*erf1);
7362
if(*x < 0.0e0) erf1 = -erf1;
7365
erf1 = fifdsign(1.0e0,*x);
7368
//****************************************************************************80
7370
double error_fc ( int *ind, double *x )
7372
//****************************************************************************80
7376
// ERROR_FC evaluates the complementary error function ERFC.
7384
// Input, int *IND, chooses the scaling.
7385
// If IND is nonzero, then the value returned has been multiplied by
7388
// Input, double *X, the argument of the function.
7390
// Output, double ERROR_FC, the value of the complementary
7394
static double c = .564189583547756e0;
7395
static double a[5] = {
7396
.771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
7397
.479137145607681e-01,.128379167095513e+00
7399
static double b[3] = {
7400
.301048631703895e-02,.538971687740286e-01,.375795757275549e+00
7402
static double p[8] = {
7403
-1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
7404
4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
7405
4.51918953711873e+02,3.00459261020162e+02
7407
static double q[8] = {
7408
1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
7409
2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7410
7.90950925327898e+02,3.00459260956983e+02
7412
static double r[5] = {
7413
2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
7414
4.65807828718470e+00,2.82094791773523e-01
7416
static double s[4] = {
7417
9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
7418
1.80124575948747e+01
7421
static double erfc1,ax,bot,e,t,top,w;
7427
if(ax > 0.5e0) goto S10;
7429
top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
7430
bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
7431
erfc1 = 0.5e0+(0.5e0-*x*(top/bot));
7432
if(*ind != 0) erfc1 = exp(t)*erfc1;
7436
// 0.5 .LT. ABS(X) .LE. 4
7438
if(ax > 4.0e0) goto S20;
7439
top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7441
bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7449
if(*x <= -5.6e0) goto S60;
7450
if(*ind != 0) goto S30;
7451
if(*x > 100.0e0) goto S70;
7452
if(*x**x > -exparg(&K1)) goto S70;
7454
t = pow(1.0e0/ *x,2.0);
7455
top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
7456
bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
7457
erfc1 = (c-t*top/bot)/ax;
7462
if(*ind == 0) goto S50;
7463
if(*x < 0.0e0) erfc1 = 2.0e0*exp(*x**x)-erfc1;
7469
erfc1 = (0.5e0+(0.5e0-e))*exp(-t)*erfc1;
7470
if(*x < 0.0e0) erfc1 = 2.0e0-erfc1;
7474
// LIMIT VALUE FOR LARGE NEGATIVE X
7477
if(*ind != 0) erfc1 = 2.0e0*exp(*x**x);
7481
// LIMIT VALUE FOR LARGE POSITIVE X
7487
//****************************************************************************80
7489
double esum ( int *mu, double *x )
7491
//****************************************************************************80
7495
// ESUM evaluates exp ( MU + X ).
7499
// Input, int *MU, part of the argument.
7501
// Input, double *X, part of the argument.
7503
// Output, double ESUM, the value of exp ( MU + X ).
7506
static double esum,w;
7508
if(*x > 0.0e0) goto S10;
7509
if(*mu < 0) goto S20;
7511
if(w > 0.0e0) goto S20;
7515
if(*mu > 0) goto S20;
7517
if(w < 0.0e0) goto S20;
7522
esum = exp(w)*exp(*x);
7525
//****************************************************************************80
7527
double eval_pol ( double a[], int *n, double *x )
7529
//****************************************************************************80
7533
// EVAL_POL evaluates a polynomial at X.
7537
// EVAL_POL = A(0) + A(1)*X + ... + A(N)*X**N
7545
// Input, double precision A(0:N), coefficients of the polynomial.
7547
// Input, int *N, length of A.
7549
// Input, double *X, the point at which the polynomial
7550
// is to be evaluated.
7552
// Output, double EVAL_POL, the value of the polynomial at X.
7555
static double devlpl,term;
7559
for ( i = *n-1-1; i >= 0; i-- )
7561
term = a[i]+term**x;
7567
//****************************************************************************80
7569
double exparg ( int *l )
7571
//****************************************************************************80
7575
// EXPARG returns the largest or smallest legal argument for EXP.
7579
// Only an approximate limit for the argument of EXP is desired.
7587
// Input, int *L, indicates which limit is desired.
7588
// If L = 0, then the largest positive argument for EXP is desired.
7589
// Otherwise, the largest negative argument for EXP for which the
7590
// result is nonzero is desired.
7592
// Output, double EXPARG, the desired value.
7598
static double exparg,lnb;
7602
if(b != 2) goto S10;
7603
lnb = .69314718055995e0;
7606
if(b != 8) goto S20;
7607
lnb = 2.0794415416798e0;
7610
if(b != 16) goto S30;
7611
lnb = 2.7725887222398e0;
7614
lnb = log((double)b);
7616
if(*l == 0) goto S50;
7618
exparg = 0.99999e0*((double)m*lnb);
7622
exparg = 0.99999e0*((double)m*lnb);
7625
//****************************************************************************80
7627
void f_cdf_values ( int *n_data, int *a, int *b, double *x, double *fx )
7629
//****************************************************************************80
7633
// F_CDF_VALUES returns some values of the F CDF test function.
7637
// The value of F_CDF ( DFN, DFD, X ) can be evaluated in Mathematica by
7640
// Needs["Statistics`ContinuousDistributions`"]
7641
// CDF[FRatioDistribution[ DFN, DFD ], X ]
7653
// Milton Abramowitz and Irene Stegun,
7654
// Handbook of Mathematical Functions,
7655
// US Department of Commerce, 1964.
7658
// The Mathematica Book,
7660
// Wolfram Media / Cambridge University Press, 1999.
7664
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
7665
// first call. On each call, the routine increments N_DATA by 1, and
7666
// returns the corresponding data; when there is no more data, the
7667
// output value of N_DATA will be 0 again.
7669
// Output, int *A, int *B, the parameters of the function.
7671
// Output, double *X, the argument of the function.
7673
// Output, double *FX, the value of the function.
7678
int a_vec[N_MAX] = {
7684
int b_vec[N_MAX] = {
7690
double fx_vec[N_MAX] = {
7691
0.500000E+00, 0.499971E+00, 0.499603E+00, 0.749699E+00,
7692
0.750466E+00, 0.751416E+00, 0.899987E+00, 0.899713E+00,
7693
0.900285E+00, 0.950025E+00, 0.950057E+00, 0.950193E+00,
7694
0.975013E+00, 0.990002E+00, 0.994998E+00, 0.999000E+00,
7695
0.568799E+00, 0.535145E+00, 0.514343E+00, 0.500000E+00 };
7696
double x_vec[N_MAX] = {
7697
1.00E+00, 0.528E+00, 1.89E+00, 1.69E+00,
7698
1.60E+00, 1.47E+00, 4.06E+00, 3.05E+00,
7699
2.09E+00, 6.61E+00, 3.71E+00, 3.00E+00,
7700
10.01E+00, 16.26E+00, 22.78E+00, 47.18E+00,
7701
1.00E+00, 1.00E+00, 1.00E+00, 1.00E+00 };
7708
*n_data = *n_data + 1;
7710
if ( N_MAX < *n_data )
7720
*a = a_vec[*n_data-1];
7721
*b = b_vec[*n_data-1];
7722
*x = x_vec[*n_data-1];
7723
*fx = fx_vec[*n_data-1];
7728
//****************************************************************************80
7730
void f_noncentral_cdf_values ( int *n_data, int *a, int *b, double *lambda,
7731
double *x, double *fx )
7733
//****************************************************************************80
7737
// F_NONCENTRAL_CDF_VALUES returns some values of the F CDF test function.
7741
// The value of NONCENTRAL_F_CDF ( DFN, DFD, LAMDA, X ) can be evaluated
7742
// in Mathematica by commands like:
7744
// Needs["Statistics`ContinuousDistributions`"]
7745
// CDF[NoncentralFRatioDistribution[ DFN, DFD, LAMBDA ], X ]
7757
// Milton Abramowitz and Irene Stegun,
7758
// Handbook of Mathematical Functions,
7759
// US Department of Commerce, 1964.
7762
// The Mathematica Book,
7764
// Wolfram Media / Cambridge University Press, 1999.
7768
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
7769
// first call. On each call, the routine increments N_DATA by 1, and
7770
// returns the corresponding data; when there is no more data, the
7771
// output value of N_DATA will be 0 again.
7773
// Output, int *A, int *B, double *LAMBDA, the
7774
// parameters of the function.
7776
// Output, double *X, the argument of the function.
7778
// Output, double *FX, the value of the function.
7783
int a_vec[N_MAX] = {
7790
int b_vec[N_MAX] = {
7797
double fx_vec[N_MAX] = {
7798
0.500000E+00, 0.636783E+00, 0.584092E+00, 0.323443E+00,
7799
0.450119E+00, 0.607888E+00, 0.705928E+00, 0.772178E+00,
7800
0.819105E+00, 0.317035E+00, 0.432722E+00, 0.450270E+00,
7801
0.426188E+00, 0.337744E+00, 0.422911E+00, 0.692767E+00,
7802
0.363217E+00, 0.421005E+00, 0.426667E+00, 0.446402E+00,
7803
0.844589E+00, 0.816368E+00 };
7804
double lambda_vec[N_MAX] = {
7805
0.00E+00, 0.000E+00, 0.25E+00, 1.00E+00,
7806
1.00E+00, 1.00E+00, 1.00E+00, 1.00E+00,
7807
1.00E+00, 2.00E+00, 1.00E+00, 1.00E+00,
7808
1.00E+00, 2.00E+00, 1.00E+00, 1.00E+00,
7809
0.00E+00, 1.00E+00, 1.00E+00, 1.00E+00,
7810
1.00E+00, 1.00E+00 };
7811
double x_vec[N_MAX] = {
7812
1.00E+00, 1.00E+00, 1.00E+00, 0.50E+00,
7813
1.00E+00, 2.00E+00, 3.00E+00, 4.00E+00,
7814
5.00E+00, 1.00E+00, 1.00E+00, 1.00E+00,
7815
1.00E+00, 1.00E+00, 1.00E+00, 2.00E+00,
7816
1.00E+00, 1.00E+00, 1.00E+00, 1.00E+00,
7817
2.00E+00, 2.00E+00 };
7824
*n_data = *n_data + 1;
7826
if ( N_MAX < *n_data )
7837
*a = a_vec[*n_data-1];
7838
*b = b_vec[*n_data-1];
7839
*lambda = lambda_vec[*n_data-1];
7840
*x = x_vec[*n_data-1];
7841
*fx = fx_vec[*n_data-1];
7847
//****************************************************************************80
7849
double fifdint ( double a )
7851
//****************************************************************************80
7855
// FIFDINT truncates a double number to an integer.
7859
// a - number to be truncated
7861
return (double) ((int) a);
7863
//****************************************************************************80
7865
double fifdmax1 ( double a, double b )
7867
//****************************************************************************80
7871
// FIFDMAX1 returns the maximum of two numbers a and b
7876
// b - second number
7888
//****************************************************************************80
7890
double fifdmin1 ( double a, double b )
7892
//****************************************************************************80
7896
// FIFDMIN1 returns the minimum of two numbers.
7901
// b - second number
7904
if (a < b) return a;
7907
//****************************************************************************80
7909
double fifdsign ( double mag, double sign )
7911
//****************************************************************************80
7915
// FIFDSIGN transfers the sign of the variable "sign" to the variable "mag"
7920
// sign - sign to be transfered
7923
if (mag < 0) mag = -mag;
7924
if (sign < 0) mag = -mag;
7928
//****************************************************************************80
7930
long fifidint ( double a )
7932
//****************************************************************************80
7936
// FIFIDINT truncates a double number to a long integer
7940
// a - number to be truncated
7952
//****************************************************************************80
7954
long fifmod ( long a, long b )
7956
//****************************************************************************80
7960
// FIFMOD returns the modulo of a and b
7970
//****************************************************************************80
7972
double fpser ( double *a, double *b, double *x, double *eps )
7974
//****************************************************************************80
7978
// FPSER evaluates IX(A,B)(X) for very small B.
7982
// This routine is appropriate for use when
7984
// B < min ( EPS, EPS * A )
7992
// Input, double *A, *B, parameters of the function.
7994
// Input, double *X, the point at which the function is to
7997
// Input, double *EPS, a tolerance.
7999
// Output, double FPSER, the value of IX(A,B)(X).
8003
static double fpser,an,c,s,t,tol;
8006
if(*a <= 1.e-3**eps) goto S10;
8009
if(t < exparg(&K1)) return fpser;
8013
// NOTE THAT 1/B(A,B) = B
8015
fpser = *b/ *a*fpser;
8025
if(fabs(c) > tol) goto S20;
8026
fpser *= (1.0e0+*a*s);
8029
//****************************************************************************80
8031
void ftnstop ( string msg )
8033
//****************************************************************************80
8037
// FTNSTOP prints a message to standard error and then exits.
8041
// Input, string MSG, the message to be printed.
8044
cerr << msg << "\n";
8048
//****************************************************************************80
8050
double gam1 ( double *a )
8052
//****************************************************************************80
8056
// GAM1 computes 1 / GAMMA(A+1) - 1 for -0.5D+00 <= A <= 1.5
8060
// Input, double *A, forms the argument of the Gamma function.
8062
// Output, double GAM1, the value of 1 / GAMMA ( A + 1 ) - 1.
8065
static double s1 = .273076135303957e+00;
8066
static double s2 = .559398236957378e-01;
8067
static double p[7] = {
8068
.577215664901533e+00,-.409078193005776e+00,-.230975380857675e+00,
8069
.597275330452234e-01,.766968181649490e-02,-.514889771323592e-02,
8070
.589597428611429e-03
8072
static double q[5] = {
8073
.100000000000000e+01,.427569613095214e+00,.158451672430138e+00,
8074
.261132021441447e-01,.423244297896961e-02
8076
static double r[9] = {
8077
-.422784335098468e+00,-.771330383816272e+00,-.244757765222226e+00,
8078
.118378989872749e+00,.930357293360349e-03,-.118290993445146e-01,
8079
.223047661158249e-02,.266505979058923e-03,-.132674909766242e-03
8081
static double gam1,bot,d,t,top,w,T1;
8085
if(d > 0.0e0) t = d-0.5e0;
8087
if(T1 < 0) goto S40;
8088
else if(T1 == 0) goto S10;
8094
top = (((((p[6]*t+p[5])*t+p[4])*t+p[3])*t+p[2])*t+p[1])*t+p[0];
8095
bot = (((q[4]*t+q[3])*t+q[2])*t+q[1])*t+1.0e0;
8097
if(d > 0.0e0) goto S30;
8101
gam1 = t/ *a*(w-0.5e0-0.5e0);
8104
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+
8106
bot = (s2*t+s1)*t+1.0e0;
8108
if(d > 0.0e0) goto S50;
8109
gam1 = *a*(w+0.5e0+0.5e0);
8115
//****************************************************************************80
8117
void gamma_inc ( double *a, double *x, double *ans, double *qans, int *ind )
8119
//****************************************************************************80
8123
// GAMMA_INC evaluates the incomplete gamma ratio functions P(A,X) and Q(A,X).
8127
// This is certified spaghetti code.
8131
// Alfred H Morris, Jr,
8132
// Naval Surface Weapons Center,
8133
// Dahlgren, Virginia.
8137
// Input, double *A, *X, the arguments of the incomplete
8138
// gamma ratio. A and X must be nonnegative. A and X cannot
8141
// Output, double *ANS, *QANS. On normal output,
8142
// ANS = P(A,X) and QANS = Q(A,X). However, ANS is set to 2 if
8143
// A or X is negative, or both are 0, or when the answer is
8144
// computationally indeterminate because A is extremely large
8145
// and X is very close to A.
8147
// Input, int *IND, indicates the accuracy request:
8148
// 0, as much accuracy as possible.
8149
// 1, to within 1 unit of the 6-th significant digit,
8150
// otherwise, to within 1 unit of the 3rd significant digit.
8153
static double alog10 = 2.30258509299405e0;
8154
static double d10 = -.185185185185185e-02;
8155
static double d20 = .413359788359788e-02;
8156
static double d30 = .649434156378601e-03;
8157
static double d40 = -.861888290916712e-03;
8158
static double d50 = -.336798553366358e-03;
8159
static double d60 = .531307936463992e-03;
8160
static double d70 = .344367606892378e-03;
8161
static double rt2pin = .398942280401433e0;
8162
static double rtpi = 1.77245385090552e0;
8163
static double third = .333333333333333e0;
8164
static double acc0[3] = {
8167
static double big[3] = {
8168
20.0e0,14.0e0,10.0e0
8170
static double d0[13] = {
8171
.833333333333333e-01,-.148148148148148e-01,.115740740740741e-02,
8172
.352733686067019e-03,-.178755144032922e-03,.391926317852244e-04,
8173
-.218544851067999e-05,-.185406221071516e-05,.829671134095309e-06,
8174
-.176659527368261e-06,.670785354340150e-08,.102618097842403e-07,
8175
-.438203601845335e-08
8177
static double d1[12] = {
8178
-.347222222222222e-02,.264550264550265e-02,-.990226337448560e-03,
8179
.205761316872428e-03,-.401877572016461e-06,-.180985503344900e-04,
8180
.764916091608111e-05,-.161209008945634e-05,.464712780280743e-08,
8181
.137863344691572e-06,-.575254560351770e-07,.119516285997781e-07
8183
static double d2[10] = {
8184
-.268132716049383e-02,.771604938271605e-03,.200938786008230e-05,
8185
-.107366532263652e-03,.529234488291201e-04,-.127606351886187e-04,
8186
.342357873409614e-07,.137219573090629e-05,-.629899213838006e-06,
8187
.142806142060642e-06
8189
static double d3[8] = {
8190
.229472093621399e-03,-.469189494395256e-03,.267720632062839e-03,
8191
-.756180167188398e-04,-.239650511386730e-06,.110826541153473e-04,
8192
-.567495282699160e-05,.142309007324359e-05
8194
static double d4[6] = {
8195
.784039221720067e-03,-.299072480303190e-03,-.146384525788434e-05,
8196
.664149821546512e-04,-.396836504717943e-04,.113757269706784e-04
8198
static double d5[4] = {
8199
-.697281375836586e-04,.277275324495939e-03,-.199325705161888e-03,
8200
.679778047793721e-04
8202
static double d6[2] = {
8203
-.592166437353694e-03,.270878209671804e-03
8205
static double e00[3] = {
8208
static double x00[3] = {
8213
static double a2n,a2nm1,acc,am0,amn,an,an0,apn,b2n,b2nm1,c,c0,c1,c2,c3,c4,c5,c6,
8214
cma,e,e0,g,h,j,l,r,rta,rtx,s,sum,t,t1,tol,twoa,u,w,x0,y,z;
8215
static int i,iop,m,max,n;
8216
static double wk[20],T3;
8218
static double T6,T7;
8221
// E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
8222
// NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
8225
if(*a < 0.0e0 || *x < 0.0e0) goto S430;
8226
if(*a == 0.0e0 && *x == 0.0e0) goto S430;
8227
if(*a**x == 0.0e0) goto S420;
8229
if(iop != 1 && iop != 2) iop = 3;
8230
acc = fifdmax1(acc0[iop-1],e);
8234
// SELECT THE APPROPRIATE ALGORITHM
8236
if(*a >= 1.0e0) goto S10;
8237
if(*a == 0.5e0) goto S390;
8238
if(*x < 1.1e0) goto S160;
8241
if(u == 0.0e0) goto S380;
8242
r = u*(1.0e0+gam1(a));
8245
if(*a >= big[iop-1]) goto S30;
8246
if(*a > *x || *x >= x0) goto S20;
8249
if(twoa != (double)m) goto S20;
8251
if(*a == (double)i) goto S210;
8255
r = exp(t1)/ gamma_x(a);
8259
if(l == 0.0e0) goto S370;
8260
s = 0.5e0+(0.5e0-l);
8262
if(z >= 700.0e0/ *a) goto S410;
8265
if(fabs(s) <= e0/rta) goto S330;
8266
if(fabs(s) <= 0.4e0) goto S270;
8267
t = pow(1.0e0/ *a,2.0);
8268
t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
8270
r = rt2pin*rta*exp(t1);
8272
if(r == 0.0e0) goto S420;
8273
if(*x <= fifdmax1(*a,alog10)) goto S50;
8274
if(*x < x0) goto S250;
8278
// TAYLOR SERIES FOR P/R
8283
for ( n = 2; n <= 20; n++ )
8287
if(t <= 1.e-3) goto S70;
8298
if(t > tol) goto S80;
8300
for ( m = 1; m <= max; m++ )
8305
*ans = r/ *a*(1.0e0+sum);
8306
*qans = 0.5e0+(0.5e0-*ans);
8310
// ASYMPTOTIC EXPANSION
8315
for ( n = 2; n <= 20; n++ )
8319
if(fabs(t) <= 1.e-3) goto S120;
8326
if(fabs(t) <= acc) goto S140;
8333
for ( m = 1; m <= max; m++ )
8338
*qans = r/ *x*(1.0e0+sum);
8339
*ans = 0.5e0+(0.5e0-*qans);
8343
// TAYLOR SERIES FOR P(A,X)/X**A
8347
sum = *x/(*a+3.0e0);
8348
tol = 3.0e0*acc/(*a+1.0e0);
8354
if(fabs(t) > tol) goto S170;
8355
j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
8359
if(*x < 0.25e0) goto S180;
8360
if(*a < *x/2.59e0) goto S200;
8363
if(z > -.13394e0) goto S200;
8366
*ans = w*g*(0.5e0+(0.5e0-j));
8367
*qans = 0.5e0+(0.5e0-*ans);
8371
w = 0.5e0+(0.5e0+l);
8372
*qans = (w*j-l)*g-h;
8373
if(*qans < 0.0e0) goto S380;
8374
*ans = 0.5e0+(0.5e0-*qans);
8378
// FINITE SUMS FOR Q WHEN A .GE. 1 AND 2*A IS AN INTEGER
8387
sum = error_fc ( &K2, &rtx );
8388
t = exp(-*x)/(rtpi*rtx);
8392
if(n == i) goto S240;
8400
*ans = 0.5e0+(0.5e0-*qans);
8404
// CONTINUED FRACTION EXPANSION
8406
tol = fifdmax1(5.0e0*e,acc);
8407
a2nm1 = a2n = 1.0e0;
8409
b2n = *x+(1.0e0-*a);
8412
a2nm1 = *x*a2n+c*a2nm1;
8413
b2nm1 = *x*b2n+c*b2nm1;
8417
a2n = a2nm1+cma*a2n;
8418
b2n = b2nm1+cma*b2n;
8420
if(fabs(an0-am0) >= tol*an0) goto S260;
8422
*ans = 0.5e0+(0.5e0-*qans);
8426
// GENERAL TEMME EXPANSION
8428
if(fabs(s) <= 2.0e0*e && *a*e*e > 3.28e-3) goto S430;
8431
w = 0.5e0 * error_fc ( &K1, &T3 );
8434
if(l < 1.0e0) z = -z;
8436
if(T4 < 0) goto S280;
8437
else if(T4 == 0) goto S290;
8440
if(fabs(s) <= 1.e-3) goto S340;
8441
c0 = ((((((((((((d0[12]*z+d0[11])*z+d0[10])*z+d0[9])*z+d0[8])*z+d0[7])*z+d0[
8442
6])*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
8443
c1 = (((((((((((d1[11]*z+d1[10])*z+d1[9])*z+d1[8])*z+d1[7])*z+d1[6])*z+d1[5]
8444
)*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8445
c2 = (((((((((d2[9]*z+d2[8])*z+d2[7])*z+d2[6])*z+d2[5])*z+d2[4])*z+d2[3])*z+
8446
d2[2])*z+d2[1])*z+d2[0])*z+d20;
8447
c3 = (((((((d3[7]*z+d3[6])*z+d3[5])*z+d3[4])*z+d3[3])*z+d3[2])*z+d3[1])*z+
8449
c4 = (((((d4[5]*z+d4[4])*z+d4[3])*z+d4[2])*z+d4[1])*z+d4[0])*z+d40;
8450
c5 = (((d5[3]*z+d5[2])*z+d5[1])*z+d5[0])*z+d50;
8451
c6 = (d6[1]*z+d6[0])*z+d60;
8452
t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
8455
c0 = (((((d0[5]*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
8456
c1 = (((d1[3]*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8461
t = ((d0[2]*z+d0[1])*z+d0[0])*z-third;
8463
if(l < 1.0e0) goto S320;
8464
*qans = c*(w+rt2pin*t/rta);
8465
*ans = 0.5e0+(0.5e0-*qans);
8468
*ans = c*(w-rt2pin*t/rta);
8469
*qans = 0.5e0+(0.5e0-*ans);
8473
// TEMME EXPANSION FOR L = 1
8475
if(*a*e*e > 3.28e-3) goto S430;
8476
c = 0.5e0+(0.5e0-y);
8477
w = (0.5e0-sqrt(y)*(0.5e0+(0.5e0-y/3.0e0))/rtpi)/c;
8480
if(l < 1.0e0) z = -z;
8482
if(T5 < 0) goto S340;
8483
else if(T5 == 0) goto S350;
8486
c0 = ((((((d0[6]*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-
8488
c1 = (((((d1[5]*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8489
c2 = ((((d2[4]*z+d2[3])*z+d2[2])*z+d2[1])*z+d2[0])*z+d20;
8490
c3 = (((d3[3]*z+d3[2])*z+d3[1])*z+d3[0])*z+d30;
8491
c4 = (d4[1]*z+d4[0])*z+d40;
8492
c5 = (d5[1]*z+d5[0])*z+d50;
8494
t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
8497
c0 = (d0[1]*z+d0[0])*z-third;
8499
t = (d20*u+c1)*u+c0;
8516
if(*x >= 0.25e0) goto S400;
8518
*ans = error_f ( &T6 );
8519
*qans = 0.5e0+(0.5e0-*ans);
8523
*qans = error_fc ( &K2, &T7 );
8524
*ans = 0.5e0+(0.5e0-*qans);
8527
if(fabs(s) <= 2.0e0*e) goto S430;
8529
if(*x <= *a) goto S370;
8538
//****************************************************************************80
8540
void gamma_inc_inv ( double *a, double *x, double *x0, double *p, double *q,
8543
//****************************************************************************80
8547
// GAMMA_INC_INV computes the inverse incomplete gamma ratio function.
8551
// The routine is given positive A, and nonnegative P and Q where P + Q = 1.
8552
// The value X is computed with the property that P(A,X) = P and Q(A,X) = Q.
8553
// Schroder iteration is employed. The routine attempts to compute X
8554
// to 10 significant digits if this is possible for the particular computer
8555
// arithmetic being used.
8559
// Alfred H Morris, Jr,
8560
// Naval Surface Weapons Center,
8561
// Dahlgren, Virginia.
8565
// Input, double *A, the parameter in the incomplete gamma
8566
// ratio. A must be positive.
8568
// Output, double *X, the computed point for which the
8569
// incomplete gamma functions have the values P and Q.
8571
// Input, double *X0, an optional initial approximation
8572
// for the solution X. If the user does not want to supply an
8573
// initial approximation, then X0 should be set to 0, or a negative
8576
// Input, double *P, *Q, the values of the incomplete gamma
8577
// functions, for which the corresponding argument is desired.
8579
// Output, int *IERR, error flag.
8580
// 0, the solution was obtained. Iteration was not used.
8581
// 0 < K, The solution was obtained. IERR iterations were performed.
8583
// -3, No solution was obtained. The ratio Q/A is too large.
8585
// -6, 20 iterations were performed. The most recent value obtained
8586
// for X is given. This cannot occur if X0 <= 0.
8587
// -7, Iteration failed. No value is given for X.
8588
// This may occur when X is approximately 0.
8589
// -8, A value for X has been obtained, but the routine is not certain
8590
// of its accuracy. Iteration cannot be performed in this
8591
// case. If X0 <= 0, this can occur only when P or Q is
8592
// approximately 0. If X0 is positive then this can occur when A is
8593
// exceedingly close to X and A is extremely large (say A .GE. 1.E20).
8596
static double a0 = 3.31125922108741e0;
8597
static double a1 = 11.6616720288968e0;
8598
static double a2 = 4.28342155967104e0;
8599
static double a3 = .213623493715853e0;
8600
static double b1 = 6.61053765625462e0;
8601
static double b2 = 6.40691597760039e0;
8602
static double b3 = 1.27364489782223e0;
8603
static double b4 = .036117081018842e0;
8604
static double c = .577215664901533e0;
8605
static double ln10 = 2.302585e0;
8606
static double tol = 1.e-5;
8607
static double amin[2] = {
8610
static double bmin[2] = {
8613
static double dmin[2] = {
8616
static double emin[2] = {
8619
static double eps0[2] = {
8626
static double am1,amax,ap1,ap2,ap3,apn,b,c1,c2,c3,c4,c5,d,e,e2,eps,g,h,pn,qg,qn,
8627
r,rta,s,s2,sum,t,u,w,xmax,xmin,xn,y,z;
8629
static double T4,T5,T6,T7,T9;
8632
// E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
8633
// E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
8634
// XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
8635
// LARGEST POSITIVE NUMBER.
8641
if(*a <= 0.0e0) goto S300;
8643
if(fabs(t) > e) goto S320;
8645
if(*p == 0.0e0) return;
8646
if(*q == 0.0e0) goto S270;
8647
if(*a == 1.0e0) goto S280;
8649
amax = 0.4e-10/(e*e);
8651
if(e > 1.e-10) iop = 2;
8654
if(*x0 > 0.0e0) goto S160;
8656
// SELECTION OF THE INITIAL APPROXIMATION XN OF X
8659
if(*a > 1.0e0) goto S80;
8663
if(qg == 0.0e0) goto S360;
8665
if(qg > 0.6e0**a) goto S40;
8666
if(*a >= 0.30e0 || b < 0.35e0) goto S10;
8672
if(b >= 0.45e0) goto S40;
8673
if(b == 0.0e0) goto S360;
8675
s = 0.5e0+(0.5e0-*a);
8678
if(b < 0.15e0) goto S20;
8679
xn = y-s*log(t)-log(1.0e0+s/(t+1.0e0));
8682
if(b <= 0.01e0) goto S30;
8683
u = ((t+2.0e0*(3.0e0-*a))*t+(2.0e0-*a)*(3.0e0-*a))/((t+(5.0e0-*a))*t+2.0e0);
8684
xn = y-s*log(t)-log(u);
8688
c2 = -(s*(1.0e0+c1));
8689
c3 = s*((0.5e0*c1+(2.0e0-*a))*c1+(2.5e0-1.5e0**a));
8690
c4 = -(s*(((c1/3.0e0+(2.5e0-1.5e0**a))*c1+((*a-6.0e0)**a+7.0e0))*c1+(
8691
(11.0e0**a-46.0)**a+47.0e0)/6.0e0));
8692
c5 = -(s*((((-(c1/4.0e0)+(11.0e0**a-17.0e0)/6.0e0)*c1+((-(3.0e0**a)+13.0e0)*
8693
*a-13.0e0))*c1+0.5e0*(((2.0e0**a-25.0e0)**a+72.0e0)**a-61.0e0))*c1+((
8694
(25.0e0**a-195.0e0)**a+477.0e0)**a-379.0e0)/12.0e0));
8695
xn = (((c5/y+c4)/y+c3)/y+c2)/y+c1+y;
8696
if(*a > 1.0e0) goto S220;
8697
if(b > bmin[iop-1]) goto S220;
8701
if(b**q > 1.e-8) goto S50;
8702
xn = exp(-(*q/ *a+c));
8705
if(*p <= 0.9e0) goto S60;
8707
xn = exp((alnrel(&T5)+ gamma_ln1 ( a ) ) / *a );
8710
xn = exp(log(*p*g)/ *a);
8712
if(xn == 0.0e0) goto S310;
8713
t = 0.5e0+(0.5e0-xn/(*a+1.0e0));
8718
// SELECTION OF THE INITIAL APPROXIMATION XN OF X
8721
if(*q <= 0.5e0) goto S90;
8727
t = sqrt(-(2.0e0*w));
8728
s = t-(((a3*t+a2)*t+a1)*t+a0)/((((b4*t+b3)*t+b2)*t+b1)*t+1.0e0);
8729
if(*q > 0.5e0) s = -s;
8732
xn = *a+s*rta+(s2-1.0e0)/3.0e0+s*(s2-7.0e0)/(36.0e0*rta)-((3.0e0*s2+7.0e0)*
8733
s2-16.0e0)/(810.0e0**a)+s*((9.0e0*s2+256.0e0)*s2-433.0e0)/(38880.0e0**a*
8735
xn = fifdmax1(xn,0.0e0);
8736
if(*a < amin[iop-1]) goto S110;
8738
d = 0.5e0+(0.5e0-*x/ *a);
8739
if(fabs(d) <= dmin[iop-1]) return;
8741
if(*p <= 0.5e0) goto S130;
8742
if(xn < 3.0e0**a) goto S220;
8743
y = -(w+ gamma_log ( a ) );
8744
d = fifdmax1(2.0e0,*a*(*a-1.0e0));
8745
if(y < ln10*d) goto S120;
8751
T6 = -(t/(xn+1.0e0));
8752
xn = y+t*log(xn)-alnrel(&T6);
8753
T7 = -(t/(xn+1.0e0));
8754
xn = y+t*log(xn)-alnrel(&T7);
8758
if(xn > 0.70e0*ap1) goto S170;
8759
w += gamma_log ( &ap1 );
8760
if(xn > 0.15e0*ap1) goto S140;
8763
*x = exp((w+*x)/ *a);
8764
*x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
8765
*x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
8766
*x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2*(1.0e0+*x/ap3))))/ *a);
8768
if(xn > 1.e-2*ap1) goto S140;
8769
if(xn <= emin[iop-1]*ap1) return;
8779
if(t > 1.e-4) goto S150;
8781
xn = exp((xn+t)/ *a);
8782
xn *= (1.0e0-(*a*log(xn)-xn-t)/(*a-xn));
8786
// SCHRODER ITERATION USING P
8788
if(*p > 0.5e0) goto S220;
8790
if(*p <= 1.e10*xmin) goto S350;
8791
am1 = *a-0.5e0-0.5e0;
8793
if(*a <= amax) goto S190;
8794
d = 0.5e0+(0.5e0-xn/ *a);
8795
if(fabs(d) <= e2) goto S350;
8797
if(*ierr >= 20) goto S330;
8799
gamma_inc ( a, &xn, &pn, &qn, &K8 );
8800
if(pn == 0.0e0 || qn == 0.0e0) goto S350;
8802
if(r == 0.0e0) goto S350;
8805
if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S200;
8807
if(*x <= 0.0e0) goto S340;
8813
if(*x <= 0.0e0) goto S340;
8814
if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
8818
if(d > tol) goto S180;
8819
if(d <= eps) return;
8820
if(fabs(*p-pn) <= tol**p) return;
8824
// SCHRODER ITERATION USING Q
8826
if(*q <= 1.e10*xmin) goto S350;
8827
am1 = *a-0.5e0-0.5e0;
8829
if(*a <= amax) goto S240;
8830
d = 0.5e0+(0.5e0-xn/ *a);
8831
if(fabs(d) <= e2) goto S350;
8833
if(*ierr >= 20) goto S330;
8835
gamma_inc ( a, &xn, &pn, &qn, &K8 );
8836
if(pn == 0.0e0 || qn == 0.0e0) goto S350;
8838
if(r == 0.0e0) goto S350;
8841
if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S250;
8843
if(*x <= 0.0e0) goto S340;
8849
if(*x <= 0.0e0) goto S340;
8850
if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
8854
if(d > tol) goto S230;
8855
if(d <= eps) return;
8856
if(fabs(*q-qn) <= tol**q) return;
8865
if(*q < 0.9e0) goto S290;
8899
//****************************************************************************80
8901
void gamma_inc_values ( int *n_data, double *a, double *x, double *fx )
8903
//****************************************************************************80
8907
// GAMMA_INC_VALUES returns some values of the incomplete Gamma function.
8911
// The (normalized) incomplete Gamma function P(A,X) is defined as:
8913
// PN(A,X) = 1/GAMMA(A) * Integral ( 0 <= T <= X ) T**(A-1) * exp(-T) dT.
8915
// With this definition, for all A and X,
8917
// 0 <= PN(A,X) <= 1
8921
// PN(A,INFINITY) = 1.0
8923
// Mathematica can compute this value as
8925
// 1 - GammaRegularized[A,X]
8937
// Milton Abramowitz and Irene Stegun,
8938
// Handbook of Mathematical Functions,
8939
// US Department of Commerce, 1964.
8943
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
8944
// first call. On each call, the routine increments N_DATA by 1, and
8945
// returns the corresponding data; when there is no more data, the
8946
// output value of N_DATA will be 0 again.
8948
// Output, double *A, the parameter of the function.
8950
// Output, double *X, the argument of the function.
8952
// Output, double *FX, the value of the function.
8957
double a_vec[N_MAX] = {
8958
0.1E+00, 0.1E+00, 0.1E+00, 0.5E+00,
8959
0.5E+00, 0.5E+00, 1.0E+00, 1.0E+00,
8960
1.0E+00, 1.1E+00, 1.1E+00, 1.1E+00,
8961
2.0E+00, 2.0E+00, 2.0E+00, 6.0E+00,
8962
6.0E+00, 11.0E+00, 26.0E+00, 41.0E+00 };
8963
double fx_vec[N_MAX] = {
8964
0.7420263E+00, 0.9119753E+00, 0.9898955E+00, 0.2931279E+00,
8965
0.7656418E+00, 0.9921661E+00, 0.0951626E+00, 0.6321206E+00,
8966
0.9932621E+00, 0.0757471E+00, 0.6076457E+00, 0.9933425E+00,
8967
0.0091054E+00, 0.4130643E+00, 0.9931450E+00, 0.0387318E+00,
8968
0.9825937E+00, 0.9404267E+00, 0.4863866E+00, 0.7359709E+00 };
8969
double x_vec[N_MAX] = {
8970
3.1622777E-02, 3.1622777E-01, 1.5811388E+00, 7.0710678E-02,
8971
7.0710678E-01, 3.5355339E+00, 0.1000000E+00, 1.0000000E+00,
8972
5.0000000E+00, 1.0488088E-01, 1.0488088E+00, 5.2440442E+00,
8973
1.4142136E-01, 1.4142136E+00, 7.0710678E+00, 2.4494897E+00,
8974
1.2247449E+01, 1.6583124E+01, 2.5495098E+01, 4.4821870E+01 };
8981
*n_data = *n_data + 1;
8983
if ( N_MAX < *n_data )
8992
*a = a_vec[*n_data-1];
8993
*x = x_vec[*n_data-1];
8994
*fx = fx_vec[*n_data-1];
8999
//****************************************************************************80
9001
double gamma_ln1 ( double *a )
9003
//****************************************************************************80
9007
// GAMMA_LN1 evaluates ln ( Gamma ( 1 + A ) ), for -0.2 <= A <= 1.25.
9011
// Input, double *A, defines the argument of the function.
9013
// Output, double GAMMA_LN1, the value of ln ( Gamma ( 1 + A ) ).
9016
static double p0 = .577215664901533e+00;
9017
static double p1 = .844203922187225e+00;
9018
static double p2 = -.168860593646662e+00;
9019
static double p3 = -.780427615533591e+00;
9020
static double p4 = -.402055799310489e+00;
9021
static double p5 = -.673562214325671e-01;
9022
static double p6 = -.271935708322958e-02;
9023
static double q1 = .288743195473681e+01;
9024
static double q2 = .312755088914843e+01;
9025
static double q3 = .156875193295039e+01;
9026
static double q4 = .361951990101499e+00;
9027
static double q5 = .325038868253937e-01;
9028
static double q6 = .667465618796164e-03;
9029
static double r0 = .422784335098467e+00;
9030
static double r1 = .848044614534529e+00;
9031
static double r2 = .565221050691933e+00;
9032
static double r3 = .156513060486551e+00;
9033
static double r4 = .170502484022650e-01;
9034
static double r5 = .497958207639485e-03;
9035
static double s1 = .124313399877507e+01;
9036
static double s2 = .548042109832463e+00;
9037
static double s3 = .101552187439830e+00;
9038
static double s4 = .713309612391000e-02;
9039
static double s5 = .116165475989616e-03;
9040
static double gamln1,w,x;
9042
if(*a >= 0.6e0) goto S10;
9043
w = ((((((p6**a+p5)**a+p4)**a+p3)**a+p2)**a+p1)**a+p0)/((((((q6**a+q5)**a+
9044
q4)**a+q3)**a+q2)**a+q1)**a+1.0e0);
9049
w = (((((r5*x+r4)*x+r3)*x+r2)*x+r1)*x+r0)/(((((s5*x+s4)*x+s3)*x+s2)*x+s1)*x
9054
//****************************************************************************80
9056
double gamma_log ( double *a )
9058
//****************************************************************************80
9062
// GAMMA_LOG evaluates ln ( Gamma ( A ) ) for positive A.
9066
// Alfred H Morris, Jr,
9067
// Naval Surface Weapons Center,
9068
// Dahlgren, Virginia.
9072
// Armido DiDinato and Alfred Morris,
9074
// Significant Digit Computation of the Incomplete Beta Function Ratios,
9075
// ACM Transactions on Mathematical Software,
9076
// Volume 18, 1993, pages 360-373.
9080
// Input, double *A, the argument of the function.
9081
// A should be positive.
9083
// Output, double GAMMA_LOG, the value of ln ( Gamma ( A ) ).
9086
static double c0 = .833333333333333e-01;
9087
static double c1 = -.277777777760991e-02;
9088
static double c2 = .793650666825390e-03;
9089
static double c3 = -.595202931351870e-03;
9090
static double c4 = .837308034031215e-03;
9091
static double c5 = -.165322962780713e-02;
9092
static double d = .418938533204673e0;
9093
static double gamln,t,w;
9097
if(*a > 0.8e0) goto S10;
9098
gamln = gamma_ln1 ( a ) - log ( *a );
9101
if(*a > 2.25e0) goto S20;
9103
gamln = gamma_ln1 ( &t );
9106
if(*a >= 10.0e0) goto S40;
9107
n = ( int ) ( *a - 1.25e0 );
9110
for ( i = 1; i <= n; i++ )
9116
gamln = gamma_ln1 ( &T1 ) + log ( w );
9119
t = pow(1.0e0/ *a,2.0);
9120
w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a;
9121
gamln = d+w+(*a-0.5e0)*(log(*a)-1.0e0);
9124
//****************************************************************************80
9126
void gamma_rat1 ( double *a, double *x, double *r, double *p, double *q,
9129
//****************************************************************************80
9133
// GAMMA_RAT1 evaluates the incomplete gamma ratio functions P(A,X) and Q(A,X).
9137
// Input, double *A, *X, the parameters of the functions.
9138
// It is assumed that A <= 1.
9140
// Input, double *R, the value exp(-X) * X**A / Gamma(A).
9142
// Output, double *P, *Q, the values of P(A,X) and Q(A,X).
9144
// Input, double *EPS, the tolerance.
9148
static double a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,t,tol,w,z,T1,T3;
9150
if(*a**x == 0.0e0) goto S120;
9151
if(*a == 0.5e0) goto S100;
9152
if(*x < 1.1e0) goto S10;
9156
// TAYLOR SERIES FOR P(A,X)/X**A
9160
sum = *x/(*a+3.0e0);
9161
tol = 0.1e0**eps/(*a+1.0e0);
9167
if(fabs(t) > tol) goto S20;
9168
j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
9172
if(*x < 0.25e0) goto S30;
9173
if(*a < *x/2.59e0) goto S50;
9176
if(z > -.13394e0) goto S50;
9179
*p = w*g*(0.5e0+(0.5e0-j));
9180
*q = 0.5e0+(0.5e0-*p);
9184
w = 0.5e0+(0.5e0+l);
9186
if(*q < 0.0e0) goto S90;
9187
*p = 0.5e0+(0.5e0-*q);
9191
// CONTINUED FRACTION EXPANSION
9193
a2nm1 = a2n = 1.0e0;
9195
b2n = *x+(1.0e0-*a);
9198
a2nm1 = *x*a2n+c*a2nm1;
9199
b2nm1 = *x*b2n+c*b2nm1;
9203
a2n = a2nm1+cma*a2n;
9204
b2n = b2nm1+cma*b2n;
9206
if(fabs(an0-am0) >= *eps*an0) goto S70;
9208
*p = 0.5e0+(0.5e0-*q);
9222
if(*x >= 0.25e0) goto S110;
9224
*p = error_f ( &T1 );
9225
*q = 0.5e0+(0.5e0-*p);
9229
*q = error_fc ( &K2, &T3 );
9230
*p = 0.5e0+(0.5e0-*q);
9233
if(*x <= *a) goto S80;
9236
//****************************************************************************80
9238
void gamma_values ( int *n_data, double *x, double *fx )
9240
//****************************************************************************80
9244
// GAMMA_VALUES returns some values of the Gamma function.
9248
// GAMMA(Z) = Integral ( 0 <= T < Infinity) T**(Z-1) EXP(-T) dT
9252
// GAMMA(X+1) = X*GAMMA(X)
9256
// 0 < X ( a software restriction).
9260
// GAMMA(0.5) = sqrt(PI)
9262
// For N a positive integer, GAMMA(N+1) = N!, the standard factorial.
9274
// Milton Abramowitz and Irene Stegun,
9275
// Handbook of Mathematical Functions,
9276
// US Department of Commerce, 1964.
9280
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
9281
// first call. On each call, the routine increments N_DATA by 1, and
9282
// returns the corresponding data; when there is no more data, the
9283
// output value of N_DATA will be 0 again.
9285
// Output, double *X, the argument of the function.
9287
// Output, double *FX, the value of the function.
9292
double fx_vec[N_MAX] = {
9293
4.590845E+00, 2.218160E+00, 1.489192E+00, 1.164230E+00,
9294
1.0000000000E+00, 0.9513507699E+00, 0.9181687424E+00, 0.8974706963E+00,
9295
0.8872638175E+00, 0.8862269255E+00, 0.8935153493E+00, 0.9086387329E+00,
9296
0.9313837710E+00, 0.9617658319E+00, 1.0000000000E+00, 3.6288000E+05,
9297
1.2164510E+17, 8.8417620E+30 };
9298
double x_vec[N_MAX] = {
9299
0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00,
9300
1.0E+00, 1.1E+00, 1.2E+00, 1.3E+00,
9301
1.4E+00, 1.5E+00, 1.6E+00, 1.7E+00,
9302
1.8E+00, 1.9E+00, 2.0E+00, 10.0E+00,
9303
20.0E+00, 30.0E+00 };
9310
*n_data = *n_data + 1;
9312
if ( N_MAX < *n_data )
9320
*x = x_vec[*n_data-1];
9321
*fx = fx_vec[*n_data-1];
9326
//****************************************************************************80
9328
double gamma_x ( double *a )
9330
//****************************************************************************80
9334
// GAMMA_X evaluates the gamma function.
9338
// This routine was renamed from "GAMMA" to avoid a conflict with the
9339
// C/C++ math library routine.
9343
// Alfred H Morris, Jr,
9344
// Naval Surface Weapons Center,
9345
// Dahlgren, Virginia.
9349
// Input, double *A, the argument of the Gamma function.
9351
// Output, double GAMMA_X, the value of the Gamma function.
9354
static double d = .41893853320467274178e0;
9355
static double pi = 3.1415926535898e0;
9356
static double r1 = .820756370353826e-03;
9357
static double r2 = -.595156336428591e-03;
9358
static double r3 = .793650663183693e-03;
9359
static double r4 = -.277777777770481e-02;
9360
static double r5 = .833333333333333e-01;
9361
static double p[7] = {
9362
.539637273585445e-03,.261939260042690e-02,.204493667594920e-01,
9363
.730981088720487e-01,.279648642639792e+00,.553413866010467e+00,1.0e0
9365
static double q[7] = {
9366
-.832979206704073e-03,.470059485860584e-02,.225211131035340e-01,
9367
-.170458969313360e+00,-.567902761974940e-01,.113062953091122e+01,1.0e0
9371
static double Xgamm,bot,g,lnx,s,t,top,w,x,z;
9372
static int i,j,m,n,T1;
9376
if(fabs(*a) >= 15.0e0) goto S110;
9378
// EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15
9383
// LET T BE THE PRODUCT OF A-J WHEN A .GE. 2
9386
if(T1 < 0) goto S40;
9387
else if(T1 == 0) goto S30;
9390
for ( j = 1; j <= m; j++ )
9400
// LET T BE THE PRODUCT OF A+J WHEN A .LT. 1
9403
if(*a > 0.0e0) goto S70;
9405
if(m == 0) goto S60;
9406
for ( j = 1; j <= m; j++ )
9414
if(t == 0.0e0) return Xgamm;
9417
// THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
9418
// CODE MAY BE OMITTED IF DESIRED.
9420
if(fabs(t) >= 1.e-30) goto S80;
9421
if(fabs(t)*dpmpar(&K2) <= 1.0001e0) return Xgamm;
9426
// COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1
9430
for ( i = 1; i < 7; i++ )
9439
if(*a < 1.0e0) goto S100;
9447
// EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15
9449
if(fabs(*a) >= 1.e3) return Xgamm;
9450
if(*a > 0.0e0) goto S120;
9454
if(t > 0.9e0) t = 1.0e0-t;
9456
if(fifmod(n,2) == 0) s = -s;
9457
if(s == 0.0e0) return Xgamm;
9460
// COMPUTE THE MODIFIED ASYMPTOTIC SUM
9463
g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x;
9465
// ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X)
9466
// BUT LESS ACCURACY WILL NORMALLY BE OBTAINED.
9473
g = d+g+(z-0.5e0)*(lnx-1.e0);
9476
if(w > 0.99999e0*exparg(&K3)) return Xgamm;
9477
Xgamm = exp(w)*(1.0e0+t);
9478
if(*a < 0.0e0) Xgamm = 1.0e0/(Xgamm*s)/x;
9481
//****************************************************************************80
9483
double gsumln ( double *a, double *b )
9485
//****************************************************************************80
9489
// GSUMLN evaluates the function ln(Gamma(A + B)).
9493
// GSUMLN is used for 1 <= A <= 2 and 1 <= B <= 2
9497
// Input, double *A, *B, values whose sum is the argument of
9498
// the Gamma function.
9500
// Output, double GSUMLN, the value of ln(Gamma(A+B)).
9503
static double gsumln,x,T1,T2;
9506
if(x > 0.25e0) goto S10;
9508
gsumln = gamma_ln1 ( &T1 );
9511
if(x > 1.25e0) goto S20;
9512
gsumln = gamma_ln1 ( &x ) + alnrel ( &x );
9516
gsumln = gamma_ln1 ( &T2 ) + log ( x * ( 1.0e0 + x ) );
9519
//****************************************************************************80
9521
int ipmpar ( int *i )
9523
//****************************************************************************80
9527
// IPMPAR returns integer machine constants.
9531
// Input arguments 1 through 3 are queries about integer arithmetic.
9532
// We assume integers are represented in the N-digit, base-A form
9534
// sign * ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) )
9536
// where 0 <= X(0:N-1) < A.
9540
// IPMPAR(1) = A, the base of integer arithmetic;
9541
// IPMPAR(2) = N, the number of base A digits;
9542
// IPMPAR(3) = A**N - 1, the largest magnitude.
9544
// It is assumed that the single and double precision floating
9545
// point arithmetics have the same base, say B, and that the
9546
// nonzero numbers are represented in the form
9548
// sign * (B**E) * (X(1)/B + ... + X(M)/B**M)
9550
// where X(1:M) is one of { 0, 1,..., B-1 }, and 1 <= X(1) and
9551
// EMIN <= E <= EMAX.
9553
// Input argument 4 is a query about the base of real arithmetic:
9555
// IPMPAR(4) = B, the base of single and double precision arithmetic.
9557
// Input arguments 5 through 7 are queries about single precision
9558
// floating point arithmetic:
9560
// IPMPAR(5) = M, the number of base B digits for single precision.
9561
// IPMPAR(6) = EMIN, the smallest exponent E for single precision.
9562
// IPMPAR(7) = EMAX, the largest exponent E for single precision.
9564
// Input arguments 8 through 10 are queries about double precision
9565
// floating point arithmetic:
9567
// IPMPAR(8) = M, the number of base B digits for double precision.
9568
// IPMPAR(9) = EMIN, the smallest exponent E for double precision.
9569
// IPMPAR(10) = EMAX, the largest exponent E for double precision.
9573
// Phyllis Fox, Andrew Hall, and Norman Schryer,
9575
// Framework for a Portable FORTRAN Subroutine Library,
9576
// ACM Transactions on Mathematical Software,
9577
// Volume 4, 1978, pages 176-188.
9581
// Input, int *I, the index of the desired constant.
9583
// Output, int IPMPAR, the value of the desired constant.
9586
static int imach[11];
9588
// MACHINE CONSTANTS FOR AMDAHL MACHINES.
9592
// imach[3] = 2147483647;
9601
// MACHINE CONSTANTS FOR THE AT&T 3B SERIES, AT&T
9602
// PC 7300, AND AT&T 6300.
9606
// imach[3] = 2147483647;
9612
// imach[9] = -1021;
9613
// imach[10] = 1024;
9615
// MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
9619
// imach[3] = 8589934591;
9628
// MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
9632
// imach[3] = 549755813887;
9641
// MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
9645
// imach[3] = 549755813887;
9651
// imach[9] = -32754;
9652
// imach[10] = 32780;
9654
// MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
9655
// 60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT
9656
// ARITHMETIC (NOS OPERATING SYSTEM).
9660
// imach[3] = 281474976710655;
9667
// imach[10] = 1070;
9669
// MACHINE CONSTANTS FOR THE CDC CYBER 995 64 BIT
9670
// ARITHMETIC (NOS/VE OPERATING SYSTEM).
9674
// imach[3] = 9223372036854775807;
9677
// imach[6] = -4096;
9680
// imach[9] = -4096;
9681
// imach[10] = 4095;
9683
// MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
9687
// imach[3] = 9223372036854775807;
9690
// imach[6] = -8189;
9693
// imach[9] = -8099;
9694
// imach[10] = 8190;
9696
// MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200.
9700
// imach[3] = 32767;
9709
// MACHINE CONSTANTS FOR THE HARRIS 220.
9713
// imach[3] = 8388607;
9722
// MACHINE CONSTANTS FOR THE HONEYWELL 600/6000
9723
// AND DPS 8/70 SERIES.
9727
// imach[3] = 34359738367;
9736
// MACHINE CONSTANTS FOR THE HP 2100
9737
// 3 WORD DOUBLE PRECISION OPTION WITH FTN4
9741
// imach[3] = 32767;
9750
// MACHINE CONSTANTS FOR THE HP 2100
9751
// 4 WORD DOUBLE PRECISION OPTION WITH FTN4
9755
// imach[3] = 32767;
9764
// MACHINE CONSTANTS FOR THE HP 9000.
9768
// imach[3] = 2147483647;
9774
// imach[9] = -1021;
9775
// imach[10] = 1024;
9777
// MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
9778
// THE ICL 2900, THE ITEL AS/6, THE XEROX SIGMA
9779
// 5/7/9 AND THE SEL SYSTEMS 85/86.
9783
// imach[3] = 2147483647;
9792
// MACHINE CONSTANTS FOR THE IBM PC.
9796
// imach[3] = 2147483647;
9802
// imach[9] = -1021;
9803
// imach[10] = 1024;
9805
// MACHINE CONSTANTS FOR THE MACINTOSH II - ABSOFT
9810
// imach[3] = 2147483647;
9816
// imach[9] = -1021;
9817
// imach[10] = 1024;
9819
// MACHINE CONSTANTS FOR THE MICROVAX - VMS FORTRAN.
9823
// imach[3] = 2147483647;
9832
// MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
9836
// imach[3] = 34359738367;
9845
// MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
9849
// imach[3] = 34359738367;
9858
// MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING
9859
// 32-BIT INTEGER ARITHMETIC.
9863
// imach[3] = 2147483647;
9872
// MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
9876
// imach[3] = 2147483647;
9882
// imach[9] = -1021;
9883
// imach[10] = 1024;
9885
// MACHINE CONSTANTS FOR THE SILICON GRAPHICS IRIS-4D
9886
// SERIES (MIPS R3000 PROCESSOR).
9890
// imach[3] = 2147483647;
9896
// imach[9] = -1021;
9897
// imach[10] = 1024;
9899
// MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
9900
// 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
9901
// PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300).
9905
imach[3] = 2147483647;
9914
// MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
9918
// imach[3] = 34359738367;
9924
// imach[9] = -1024;
9925
// imach[10] = 1023;
9927
// MACHINE CONSTANTS FOR THE VAX 11/780.
9931
// imach[3] = 2147483647;
9943
//****************************************************************************80
9945
void negative_binomial_cdf_values ( int *n_data, int *f, int *s, double *p,
9948
//****************************************************************************80
9952
// NEGATIVE_BINOMIAL_CDF_VALUES returns values of the negative binomial CDF.
9956
// Assume that a coin has a probability P of coming up heads on
9957
// any one trial. Suppose that we plan to flip the coin until we
9958
// achieve a total of S heads. If we let F represent the number of
9959
// tails that occur in this process, then the value of F satisfies
9960
// a negative binomial PDF:
9962
// PDF(F,S,P) = Choose ( F from F+S-1 ) * P**S * (1-P)**F
9964
// The negative binomial CDF is the probability that there are F or
9965
// fewer failures upon the attainment of the S-th success. Thus,
9967
// CDF(F,S,P) = sum ( 0 <= G <= F ) PDF(G,S,P)
9980
// Statistical Tables for Sociology, Biology and Physical Sciences,
9981
// Cambridge University Press, 1982.
9985
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
9986
// first call. On each call, the routine increments N_DATA by 1, and
9987
// returns the corresponding data; when there is no more data, the
9988
// output value of N_DATA will be 0 again.
9990
// Output, int *F, the maximum number of failures.
9992
// Output, int *S, the number of successes.
9994
// Output, double *P, the probability of a success on one trial.
9996
// Output, double *CDF, the probability of at most F failures before the
10002
double cdf_vec[N_MAX] = {
10003
0.6367, 0.3633, 0.1445,
10004
0.5000, 0.2266, 0.0625,
10005
0.3438, 0.1094, 0.0156,
10006
0.1792, 0.0410, 0.0041,
10007
0.0705, 0.0109, 0.0007,
10008
0.9862, 0.9150, 0.7472,
10009
0.8499, 0.5497, 0.2662,
10010
0.6513, 0.2639, 0.0702,
10011
1.0000, 0.0199, 0.0001 };
10012
int f_vec[N_MAX] = {
10022
double p_vec[N_MAX] = {
10031
0.01, 0.01, 0.01 };
10032
int s_vec[N_MAX] = {
10048
*n_data = *n_data + 1;
10050
if ( N_MAX < *n_data )
10060
*f = f_vec[*n_data-1];
10061
*s = s_vec[*n_data-1];
10062
*p = p_vec[*n_data-1];
10063
*cdf = cdf_vec[*n_data-1];
10069
//****************************************************************************80
10071
void normal_cdf_values ( int *n_data, double *x, double *fx )
10073
//****************************************************************************80
10077
// NORMAL_CDF_VALUES returns some values of the Normal CDF.
10089
// Milton Abramowitz and Irene Stegun,
10090
// Handbook of Mathematical Functions,
10091
// US Department of Commerce, 1964.
10095
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
10096
// first call. On each call, the routine increments N_DATA by 1, and
10097
// returns the corresponding data; when there is no more data, the
10098
// output value of N_DATA will be 0 again.
10100
// Output, double *X, the argument of the function.
10102
// Output double *FX, the value of the function.
10107
double fx_vec[N_MAX] = {
10108
0.500000000000000E+00, 0.539827837277029E+00, 0.579259709439103E+00,
10109
0.617911422188953E+00, 0.655421741610324E+00, 0.691462461274013E+00,
10110
0.725746882249927E+00, 0.758036347776927E+00, 0.788144601416604E+00,
10111
0.815939874653241E+00, 0.841344746068543E+00, 0.933192798731142E+00,
10112
0.977249868051821E+00 };
10113
double x_vec[N_MAX] = {
10114
0.00E+00, 0.10E+00, 0.20E+00,
10115
0.30E+00, 0.40E+00, 0.50E+00,
10116
0.60E+00, 0.70E+00, 0.80E+00,
10117
0.90E+00, 1.00E+00, 1.50E+00,
10125
*n_data = *n_data + 1;
10127
if ( N_MAX < *n_data )
10135
*x = x_vec[*n_data-1];
10136
*fx = fx_vec[*n_data-1];
10142
//****************************************************************************80
10144
void poisson_cdf_values ( int *n_data, double *a, int *x, double *fx )
10146
//****************************************************************************80
10150
// POISSON_CDF_VALUES returns some values of the Poisson CDF.
10154
// CDF(X)(A) is the probability of at most X successes in unit time,
10155
// given that the expected mean number of successes is A.
10167
// Milton Abramowitz and Irene Stegun,
10168
// Handbook of Mathematical Functions,
10169
// US Department of Commerce, 1964.
10171
// Daniel Zwillinger,
10172
// CRC Standard Mathematical Tables and Formulae,
10173
// 30th Edition, CRC Press, 1996, pages 653-658.
10177
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
10178
// first call. On each call, the routine increments N_DATA by 1, and
10179
// returns the corresponding data; when there is no more data, the
10180
// output value of N_DATA will be 0 again.
10182
// Output, double *A, the parameter of the function.
10184
// Output, int *X, the argument of the function.
10186
// Output, double *FX, the value of the function.
10191
double a_vec[N_MAX] = {
10192
0.02E+00, 0.10E+00, 0.10E+00, 0.50E+00,
10193
0.50E+00, 0.50E+00, 1.00E+00, 1.00E+00,
10194
1.00E+00, 1.00E+00, 2.00E+00, 2.00E+00,
10195
2.00E+00, 2.00E+00, 5.00E+00, 5.00E+00,
10196
5.00E+00, 5.00E+00, 5.00E+00, 5.00E+00,
10198
double fx_vec[N_MAX] = {
10199
0.980E+00, 0.905E+00, 0.995E+00, 0.607E+00,
10200
0.910E+00, 0.986E+00, 0.368E+00, 0.736E+00,
10201
0.920E+00, 0.981E+00, 0.135E+00, 0.406E+00,
10202
0.677E+00, 0.857E+00, 0.007E+00, 0.040E+00,
10203
0.125E+00, 0.265E+00, 0.441E+00, 0.616E+00,
10205
int x_vec[N_MAX] = {
10218
*n_data = *n_data + 1;
10220
if ( N_MAX < *n_data )
10229
*a = a_vec[*n_data-1];
10230
*x = x_vec[*n_data-1];
10231
*fx = fx_vec[*n_data-1];
10236
//****************************************************************************80
10238
double psi ( double *xx )
10240
//****************************************************************************80
10244
// PSI evaluates the psi or digamma function, d/dx ln(gamma(x)).
10248
// The main computation involves evaluation of rational Chebyshev
10249
// approximations. PSI was written at Argonne National Laboratory
10250
// for FUNPACK, and subsequently modified by A. H. Morris of NSWC.
10254
// William Cody, Strecok and Thacher,
10255
// Chebyshev Approximations for the Psi Function,
10256
// Mathematics of Computation,
10257
// Volume 27, 1973, pages 123-127.
10261
// Input, double *XX, the argument of the psi function.
10263
// Output, double PSI, the value of the psi function. PSI
10264
// is assigned the value 0 when the psi function is undefined.
10267
static double dx0 = 1.461632144968362341262659542325721325e0;
10268
static double piov4 = .785398163397448e0;
10269
static double p1[7] = {
10270
.895385022981970e-02,.477762828042627e+01,.142441585084029e+03,
10271
.118645200713425e+04,.363351846806499e+04,.413810161269013e+04,
10272
.130560269827897e+04
10274
static double p2[4] = {
10275
-.212940445131011e+01,-.701677227766759e+01,-.448616543918019e+01,
10276
-.648157123766197e+00
10278
static double q1[6] = {
10279
.448452573429826e+02,.520752771467162e+03,.221000799247830e+04,
10280
.364127349079381e+04,.190831076596300e+04,.691091682714533e-05
10282
static double q2[4] = {
10283
.322703493791143e+02,.892920700481861e+02,.546117738103215e+02,
10284
.777788548522962e+01
10288
static double psi,aug,den,sgn,upper,w,x,xmax1,xmx0,xsmall,z;
10289
static int i,m,n,nq;
10291
// MACHINE DEPENDENT CONSTANTS ...
10292
// XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
10293
// WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED
10294
// AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
10295
// ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
10296
// PSI MAY BE REPRESENTED AS ALOG(X).
10297
// XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
10298
// MAY BE REPRESENTED BY 1/X.
10300
xmax1 = ipmpar(&K1);
10301
xmax1 = fifdmin1(xmax1,1.0e0/dpmpar(&K2));
10305
if(x >= 0.5e0) goto S50;
10307
// X .LT. 0.5, USE REFLECTION FORMULA
10308
// PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
10310
if(fabs(x) > xsmall) goto S10;
10311
if(x == 0.0e0) goto S100;
10313
// 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE
10314
// FOR PI*COTAN(PI*X)
10320
// REDUCTION OF ARGUMENT FOR COTAN
10324
if(w > 0.0e0) goto S20;
10329
// MAKE AN ERROR EXIT IF X .LE. -XMAX1
10331
if(w >= xmax1) goto S100;
10334
nq = fifidint(w*4.0e0);
10335
w = 4.0e0*(w-(double)nq*.25e0);
10337
// W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X.
10338
// ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
10339
// QUADRANT AND DETERMINE SIGN
10342
if(n+n != nq) w = 1.0e0-w;
10345
if(m+m != n) sgn = -sgn;
10347
// DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X)
10352
if(m != n) goto S30;
10354
// CHECK FOR SINGULARITY
10356
if(z == 0.0e0) goto S100;
10358
// USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
10359
// SIN/COS AS A SUBSTITUTE FOR TAN
10361
aug = sgn*(cos(z)/sin(z)*4.0e0);
10364
aug = sgn*(sin(z)/cos(z)*4.0e0);
10368
if(x > 3.0e0) goto S70;
10370
// 0.5 .LE. X .LE. 3.0
10374
for ( i = 1; i <= 5; i++ )
10376
den = (den+q1[i-1])*x;
10377
upper = (upper+p1[i+1-1])*x;
10379
den = (upper+p1[6])/(den+q1[5]);
10381
psi = den*xmx0+aug;
10385
// IF X .GE. XMAX1, PSI = LN(X)
10387
if(x >= xmax1) goto S90;
10389
// 3.0 .LT. X .LT. XMAX1
10394
for ( i = 1; i <= 3; i++ )
10396
den = (den+q2[i-1])*w;
10397
upper = (upper+p2[i+1-1])*w;
10399
aug = upper/(den+q2[3])-0.5e0/x+aug;
10410
//****************************************************************************80
10412
void psi_values ( int *n_data, double *x, double *fx )
10414
//****************************************************************************80
10418
// PSI_VALUES returns some values of the Psi or Digamma function.
10422
// PSI(X) = d LN ( Gamma ( X ) ) / d X = Gamma'(X) / Gamma(X)
10424
// PSI(1) = - Euler's constant.
10426
// PSI(X+1) = PSI(X) + 1 / X.
10438
// Milton Abramowitz and Irene Stegun,
10439
// Handbook of Mathematical Functions,
10440
// US Department of Commerce, 1964.
10444
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
10445
// first call. On each call, the routine increments N_DATA by 1, and
10446
// returns the corresponding data; when there is no more data, the
10447
// output value of N_DATA will be 0 again.
10449
// Output, double *X, the argument of the function.
10451
// Output, double *FX, the value of the function.
10456
double fx_vec[N_MAX] = {
10457
-0.5772156649E+00, -0.4237549404E+00, -0.2890398966E+00,
10458
-0.1691908889E+00, -0.0613845446E+00, -0.0364899740E+00,
10459
0.1260474528E+00, 0.2085478749E+00, 0.2849914333E+00,
10460
0.3561841612E+00, 0.4227843351E+00 };
10461
double x_vec[N_MAX] = {
10462
1.0E+00, 1.1E+00, 1.2E+00,
10463
1.3E+00, 1.4E+00, 1.5E+00,
10464
1.6E+00, 1.7E+00, 1.8E+00,
10465
1.9E+00, 2.0E+00 };
10472
*n_data = *n_data + 1;
10474
if ( N_MAX < *n_data )
10482
*x = x_vec[*n_data-1];
10483
*fx = fx_vec[*n_data-1];
10488
//****************************************************************************80
10490
double rcomp ( double *a, double *x )
10492
//****************************************************************************80
10496
// RCOMP evaluates exp(-X) * X**A / Gamma(A).
10500
// Input, double *A, *X, arguments of the quantity to be computed.
10502
// Output, double RCOMP, the value of exp(-X) * X**A / Gamma(A).
10504
// Local parameters:
10506
// RT2PIN = 1/SQRT(2*PI)
10509
static double rt2pin = .398942280401433e0;
10510
static double rcomp,t,t1,u;
10512
if(*a >= 20.0e0) goto S20;
10514
if(*a >= 1.0e0) goto S10;
10515
rcomp = *a*exp(t)*(1.0e0+gam1(a));
10518
rcomp = exp(t)/ gamma_x(a);
10522
if(u == 0.0e0) return rcomp;
10523
t = pow(1.0e0/ *a,2.0);
10524
t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
10525
t1 -= (*a*rlog(&u));
10526
rcomp = rt2pin*sqrt(*a)*exp(t1);
10529
//****************************************************************************80
10531
double rexp ( double *x )
10533
//****************************************************************************80
10537
// REXP evaluates the function EXP(X) - 1.
10541
// 09 December 1999
10545
// Input, double *X, the argument of the function.
10547
// Output, double REXP, the value of EXP(X)-1.
10550
static double p1 = .914041914819518e-09;
10551
static double p2 = .238082361044469e-01;
10552
static double q1 = -.499999999085958e+00;
10553
static double q2 = .107141568980644e+00;
10554
static double q3 = -.119041179760821e-01;
10555
static double q4 = .595130811860248e-03;
10556
static double rexp,w;
10558
if(fabs(*x) > 0.15e0) goto S10;
10559
rexp = *x*(((p2**x+p1)**x+1.0e0)/((((q4**x+q3)**x+q2)**x+q1)**x+1.0e0));
10563
if(*x > 0.0e0) goto S20;
10564
rexp = w-0.5e0-0.5e0;
10567
rexp = w*(0.5e0+(0.5e0-1.0e0/w));
10570
//****************************************************************************80
10572
double rlog ( double *x )
10574
//****************************************************************************80
10578
// RLOG computes X - 1 - LN(X).
10582
// 09 December 1999
10586
// Input, double *X, the argument of the function.
10588
// Output, double RLOG, the value of the function.
10591
static double a = .566749439387324e-01;
10592
static double b = .456512608815524e-01;
10593
static double p0 = .333333333333333e+00;
10594
static double p1 = -.224696413112536e+00;
10595
static double p2 = .620886815375787e-02;
10596
static double q1 = -.127408923933623e+01;
10597
static double q2 = .354508718369557e+00;
10598
static double rlog,r,t,u,w,w1;
10600
if(*x < 0.61e0 || *x > 1.57e0) goto S40;
10601
if(*x < 0.82e0) goto S10;
10602
if(*x > 1.18e0) goto S20;
10604
// ARGUMENT REDUCTION
10606
u = *x-0.5e0-0.5e0;
10615
u = 0.75e0**x-1.e0;
10619
// SERIES EXPANSION
10623
w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
10624
rlog = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
10627
r = *x-0.5e0-0.5e0;
10631
//****************************************************************************80
10633
double rlog1 ( double *x )
10635
//****************************************************************************80
10639
// RLOG1 evaluates the function X - ln ( 1 + X ).
10643
// Input, double *X, the argument.
10645
// Output, double RLOG1, the value of X - ln ( 1 + X ).
10648
static double a = .566749439387324e-01;
10649
static double b = .456512608815524e-01;
10650
static double p0 = .333333333333333e+00;
10651
static double p1 = -.224696413112536e+00;
10652
static double p2 = .620886815375787e-02;
10653
static double q1 = -.127408923933623e+01;
10654
static double q2 = .354508718369557e+00;
10655
static double rlog1,h,r,t,w,w1;
10657
if(*x < -0.39e0 || *x > 0.57e0) goto S40;
10658
if(*x < -0.18e0) goto S10;
10659
if(*x > 0.18e0) goto S20;
10661
// ARGUMENT REDUCTION
10672
h = 0.75e0**x-0.25e0;
10676
// SERIES EXPANSION
10680
w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
10681
rlog1 = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
10684
w = *x+0.5e0+0.5e0;
10688
//****************************************************************************80
10690
void student_cdf_values ( int *n_data, int *a, double *x, double *fx )
10692
//****************************************************************************80
10696
// STUDENT_CDF_VALUES returns some values of the Student CDF.
10708
// Milton Abramowitz and Irene Stegun,
10709
// Handbook of Mathematical Functions,
10710
// US Department of Commerce, 1964.
10714
// Input/output, int *N_DATA. The user sets N_DATA to 0 before the
10715
// first call. On each call, the routine increments N_DATA by 1, and
10716
// returns the corresponding data; when there is no more data, the
10717
// output value of N_DATA will be 0 again.
10719
// Output, int *A, the parameter of the function.
10721
// Output, double *X, the argument of the function.
10723
// Output, double *FX, the value of the function.
10728
int a_vec[N_MAX] = {
10733
double fx_vec[N_MAX] = {
10734
0.60E+00, 0.60E+00, 0.60E+00, 0.60E+00,
10735
0.60E+00, 0.75E+00, 0.75E+00, 0.95E+00,
10736
0.95E+00, 0.99E+00, 0.99E+00, 0.99E+00,
10738
double x_vec[N_MAX] = {
10739
0.325E+00, 0.289E+00, 0.277E+00, 0.271E+00,
10740
0.267E+00, 0.816E+00, 0.727E+00, 2.920E+00,
10741
2.015E+00, 6.965E+00, 4.541E+00, 3.747E+00,
10749
*n_data = *n_data + 1;
10751
if ( N_MAX < *n_data )
10760
*a = a_vec[*n_data-1];
10761
*x = x_vec[*n_data-1];
10762
*fx = fx_vec[*n_data-1];
10768
//****************************************************************************80
10770
double stvaln ( double *p )
10772
//****************************************************************************80
10776
// STVALN provides starting values for the inverse of the normal distribution.
10780
// The routine returns X such that
10783
// P = Integral from -infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU.
10787
// Kennedy and Gentle,
10788
// Statistical Computing,
10789
// Marcel Dekker, NY, 1980, page 95,
10794
// Input, double *P, the probability whose normal deviate
10797
// Output, double STVALN, the normal deviate whose probability
10801
static double xden[5] = {
10802
0.993484626060e-1,0.588581570495e0,0.531103462366e0,0.103537752850e0,
10805
static double xnum[5] = {
10806
-0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1,
10810
static double stvaln,sign,y,z;
10812
if(!(*p <= 0.5e0)) goto S10;
10820
y = sqrt(-(2.0e0*log(z)));
10821
stvaln = y+ eval_pol ( xnum, &K1, &y ) / eval_pol ( xden, &K1, &y );
10822
stvaln = sign*stvaln;
10825
//**************************************************************************80
10827
#if !defined(TIMESTAMP)
10831
//**************************************************************************80
10835
// TIMESTAMP prints the current YMDHMS date as a time stamp.
10839
// May 31 2001 09:45:54 AM
10843
// 24 September 2003
10854
# define TIME_SIZE 40
10856
static char time_buffer[TIME_SIZE];
10857
const struct tm *tm;
10861
now = time ( NULL );
10862
tm = localtime ( &now );
10864
len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
10866
cout << time_buffer << "\n";