6 static void E0000(int,int*,double*,double*,unsigned long*,
7 unsigned long*,double*,double*,double*,
8 double*,double*,double*,double*);
9 static void E0001(int,int*,double*,double*,double*,double*,
10 unsigned long*,unsigned long*,double*,double*,
14 * A comment about ints and longs - whether ints or longs are used should
15 * make no difference, but where double r-values are assigned to ints the
16 * r-value is cast converted to a long, which is then assigned to the int
17 * to be compatible with the operation of fifidint.
20 -----------------------------------------------------------------------
22 COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
26 IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
27 LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
29 -----------------------------------------------------------------------
31 double algdiv(double *a,double *b)
33 static double c0 = .833333333333333e-01;
34 static double c1 = -.277777777760991e-02;
35 static double c2 = .793650666825390e-03;
36 static double c3 = -.595202931351870e-03;
37 static double c4 = .837308034031215e-03;
38 static double c5 = -.165322962780713e-02;
39 static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;
42 .. Executable Statements ..
44 if(*a <= *b) goto S10;
57 SET SN = (1 - X**N)/(1 - X)
64 s11 = 1.0e0+(x+x2*s9);
66 SET W = DEL(B) - DEL(A + B)
68 t = pow(1.0e0/ *b,2.0);
69 w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
76 v = *a*(log(*b)-1.0e0);
84 double alngam(double *x)
86 **********************************************************************
88 double alngam(double *x)
89 double precision LN of the GAMma function
95 Returns the natural logarithm of GAMMA(X).
101 X --> value at which scaled log gamma is to be returned
102 X is DOUBLE PRECISION
108 If X .le. 6.0, then use recursion to get X below 3
109 then apply rational approximation number 5236 of
110 Hart et al, Computer Approximations, John Wiley and
113 If X .gt. 6.0, then use recursion to get X to at least 12 and
114 then use formula 5423 of the same source.
116 **********************************************************************
119 #define hln2pi 0.91893853320467274178e0
120 static double coef[5] = {
121 0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
122 -0.594997310889e-3,0.8065880899e-3
124 static double scoefd[4] = {
125 0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
126 0.1000000000000000000e1
128 static double scoefn[9] = {
129 0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
130 0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
131 0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
136 static double alngam,offset,prod,xx;
138 static double T2,T4,T6;
141 .. Executable Statements ..
143 if(!(*x <= 6.0e0)) goto S70;
146 if(!(*x > 3.0e0)) goto S30;
148 if(!(xx > 3.0e0)) goto S20;
154 if(!(*x < 2.0e0)) goto S60;
156 if(!(xx < 2.0e0)) goto S50;
164 alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
166 COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
169 alngam = log(alngam);
174 IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
176 n = fifidint(12.0e0-*x);
177 if(!(n > 0)) goto S90;
179 for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));
189 T6 = 1.0e0/pow(xx,2.0);
190 alngam = devlpl(coef,&K5,&T6)/xx;
191 alngam += (offset+(xx-0.5e0)*log(xx)-xx);
196 double alnrel(double *a)
198 -----------------------------------------------------------------------
199 EVALUATION OF THE FUNCTION LN(1 + A)
200 -----------------------------------------------------------------------
203 static double p1 = -.129418923021993e+01;
204 static double p2 = .405303492862024e+00;
205 static double p3 = -.178874546012214e-01;
206 static double q1 = -.162752256355323e+01;
207 static double q2 = .747811014037616e+00;
208 static double q3 = -.845104217945565e-01;
209 static double alnrel,t,t2,w,x;
212 .. Executable Statements ..
214 if(fabs(*a) > 0.375e0) goto S10;
217 w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);
225 double apser(double *a,double *b,double *x,double *eps)
227 -----------------------------------------------------------------------
228 APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
229 A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
230 A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
231 -----------------------------------------------------------------------
234 static double g = .577215664901533e0;
235 static double apser,aj,bx,c,j,s,t,tol;
238 .. Executable Statements ..
242 if(*b**eps > 2.e-2) goto S10;
243 c = log(*x)+psi(b)+g+t;
248 tol = 5.0e0**eps*fabs(c);
256 if(fabs(aj) > tol) goto S30;
260 double basym(double *a,double *b,double *lambda,double *eps)
262 -----------------------------------------------------------------------
263 ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
264 LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED.
265 IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
266 A AND B ARE GREATER THAN OR EQUAL TO 15.
267 -----------------------------------------------------------------------
270 static double e0 = 1.12837916709551e0;
271 static double e1 = .353553390593274e0;
274 ------------------------
275 ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
276 ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
277 THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
278 ------------------------
281 ------------------------
284 static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
286 static int i,im1,imj,j,m,mm1,mmj,n,np1;
287 static double a0[21],b0[21],c[21],d[21],T1,T2;
290 .. Executable Statements ..
293 if(*a >= *b) goto S10;
295 r0 = 1.0e0/(1.0e0+h);
297 w0 = 1.0e0/sqrt(*a*(1.0e0+h));
301 r0 = 1.0e0/(1.0e0+h);
303 w0 = 1.0e0/sqrt(*b*(1.0e0+h));
307 f = *a*rlog1(&T1)+*b*rlog1(&T2);
309 if(t == 0.0e0) return basym;
313 a0[0] = 2.0e0/3.0e0*r1;
314 c[0] = -(0.5e0*a0[0]);
316 j0 = 0.5e0/e0*erfc1(&K3,&z0);
325 for(n=2; n<=num; n+=2) {
327 a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
330 a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
331 for(i=n; i<=np1; i++) {
332 r = -(0.5e0*((double)i+1.0e0));
334 for(m=2; m<=i; m++) {
337 for(j=1; j<=mm1; j++) {
339 bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
341 b0[m-1] = r*a0[m-1]+bsum/(double)m;
343 c[i-1] = b0[i-1]/((double)i+1.0e0);
346 for(j=1; j<=im1; j++) {
348 dsum += (d[imj-1]*c[j-1]);
350 d[i-1] = -(dsum+c[i-1]);
352 j0 = e1*znm1+((double)n-1.0e0)*j0;
353 j1 = e1*zn+(double)n*j1;
361 if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
364 u = exp(-bcorr(a,b));
368 double bcorr(double *a0,double *b0)
370 -----------------------------------------------------------------------
372 EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE
373 LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
374 IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8.
376 -----------------------------------------------------------------------
379 static double c0 = .833333333333333e-01;
380 static double c1 = -.277777777760991e-02;
381 static double c2 = .793650666825390e-03;
382 static double c3 = -.595202931351870e-03;
383 static double c4 = .837308034031215e-03;
384 static double c5 = -.165322962780713e-02;
385 static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;
388 .. Executable Statements ..
390 a = fifdmin1(*a0,*b0);
391 b = fifdmax1(*a0,*b0);
397 SET SN = (1 - X**N)/(1 - X)
400 s5 = 1.0e0+(x+x2*s3);
401 s7 = 1.0e0+(x+x2*s5);
402 s9 = 1.0e0+(x+x2*s7);
403 s11 = 1.0e0+(x+x2*s9);
405 SET W = DEL(B) - DEL(A + B)
407 t = pow(1.0e0/b,2.0);
408 w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
413 t = pow(1.0e0/a,2.0);
414 bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w;
417 double betaln(double *a0,double *b0)
419 -----------------------------------------------------------------------
420 EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
421 -----------------------------------------------------------------------
423 --------------------------
426 static double e = .918938533204673e0;
427 static double betaln,a,b,c,h,u,v,w,z;
432 .. Executable Statements ..
434 a = fifdmin1(*a0,*b0);
435 b = fifdmax1(*a0,*b0);
436 if(a >= 8.0e0) goto S100;
437 if(a >= 1.0e0) goto S20;
439 -----------------------------------------------------------------------
440 PROCEDURE WHEN A .LT. 1
441 -----------------------------------------------------------------------
443 if(b >= 8.0e0) goto S10;
445 betaln = gamln(&a)+(gamln(&b)-gamln(&T1));
448 betaln = gamln(&a)+algdiv(&a,&b);
452 -----------------------------------------------------------------------
453 PROCEDURE WHEN 1 .LE. A .LT. 8
454 -----------------------------------------------------------------------
456 if(a > 2.0e0) goto S40;
457 if(b > 2.0e0) goto S30;
458 betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b);
462 if(b < 8.0e0) goto S60;
463 betaln = gamln(&a)+algdiv(&a,&b);
467 REDUCTION OF A WHEN B .LE. 1000
469 if(b > 1000.0e0) goto S80;
470 n = (long)(a - 1.0e0);
472 for(i=1; i<=n; i++) {
478 if(b < 8.0e0) goto S60;
479 betaln = w+gamln(&a)+algdiv(&a,&b);
483 REDUCTION OF B WHEN B .LT. 8
485 n = (long)(b - 1.0e0);
487 for(i=1; i<=n; i++) {
491 betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
495 REDUCTION OF A WHEN B .GT. 1000
497 n = (long)(a - 1.0e0);
499 for(i=1; i<=n; i++) {
501 w *= (a/(1.0e0+a/b));
503 betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
507 -----------------------------------------------------------------------
508 PROCEDURE WHEN A .GE. 8
509 -----------------------------------------------------------------------
514 u = -((a-0.5e0)*log(c));
516 if(u <= v) goto S110;
517 betaln = -(0.5e0*log(b))+e+w-v-u;
520 betaln = -(0.5e0*log(b))+e+w-u-v;
523 double bfrac(double *a,double *b,double *x,double *y,double *lambda,
526 -----------------------------------------------------------------------
527 CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
528 IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B.
529 -----------------------------------------------------------------------
532 static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;
535 .. Executable Statements ..
537 bfrac = brcomp(a,b,x,y);
538 if(bfrac == 0.0e0) return bfrac;
541 c1 = 1.0e0+1.0e0/ *a;
552 CONTINUED FRACTION CALCULATION
558 alpha = p*(p+c0)*e*e*(w**x);
559 e = (1.0e0+t)/(c1+t+t);
560 beta = n+w/s+e*(c+n*yp1);
564 UPDATE AN, BN, ANP1, AND BNP1
566 t = alpha*an+beta*anp1;
569 t = alpha*bn+beta*bnp1;
574 if(fabs(r-r0) <= *eps*r) goto S20;
576 RESCALE AN, BN, ANP1, AND BNP1
590 void bgrat(double *a,double *b,double *x,double *y,double *w,
591 double *eps,int *ierr)
593 -----------------------------------------------------------------------
594 ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
595 THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
596 THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED.
597 IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
598 -----------------------------------------------------------------------
601 static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
603 static double c[30],d[30],T1;
606 .. Executable Statements ..
608 bm1 = *b-0.5e0-0.5e0;
610 if(*y > 0.375e0) goto S10;
618 if(*b*z == 0.0e0) goto S70;
620 COMPUTATION OF THE EXPANSION
621 SET R = EXP(-Z)*Z**B/GAMMA(B)
623 r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
624 r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
625 u = algdiv(b,a)+*b*log(nu);
627 if(u == 0.0e0) goto S70;
628 grat1(b,&z,&r,&p,&q,eps);
629 v = 0.25e0*pow(1.0e0/nu,2.0);
636 for(n=1; n<=30; n++) {
638 j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
641 cn /= (n2*(n2+1.0e0));
647 for(i=1; i<=nm1; i++) {
648 s += (coef*c[i-1]*d[n-i-1]);
652 d[n-1] = bm1*cn+s/(double)n;
655 if(sum <= 0.0e0) goto S70;
656 if(fabs(dj) <= *eps*(sum+l)) goto S60;
667 THE EXPANSION CANNOT BE COMPUTED
672 double bpser(double *a,double *b,double *x,double *eps)
674 -----------------------------------------------------------------------
675 POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
676 OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED.
677 -----------------------------------------------------------------------
680 static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
684 .. Executable Statements ..
687 if(*x == 0.0e0) return bpser;
689 -----------------------------------------------------------------------
690 COMPUTE THE FACTOR X**A/(A*BETA(A,B))
691 -----------------------------------------------------------------------
693 a0 = fifdmin1(*a,*b);
694 if(a0 < 1.0e0) goto S10;
695 z = *a*log(*x)-betaln(a,b);
699 b0 = fifdmax1(*a,*b);
700 if(b0 >= 8.0e0) goto S90;
701 if(b0 > 1.0e0) goto S40;
703 PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
706 if(bpser == 0.0e0) return bpser;
708 if(apb > 1.0e0) goto S20;
709 z = 1.0e0+gam1(&apb);
713 z = (1.0e0+gam1(&u))/apb;
715 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
716 bpser *= (c*(*b/apb));
720 PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
723 m = (long)(b0 - 1.0e0);
726 for(i=1; i<=m; i++) {
735 if(apb > 1.0e0) goto S70;
736 t = 1.0e0+gam1(&apb);
740 t = (1.0e0+gam1(&u))/apb;
742 bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
746 PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
748 u = gamln1(&a0)+algdiv(&a0,&b0);
750 bpser = a0/ *a*exp(z);
752 if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
754 -----------------------------------------------------------------------
756 -----------------------------------------------------------------------
763 c *= ((0.5e0+(0.5e0-*b/n))**x);
766 if(fabs(w) > tol) goto S110;
767 bpser *= (1.0e0+*a*sum);
770 void bratio(double *a,double *b,double *x,double *y,double *w,
771 double *w1,int *ierr)
773 -----------------------------------------------------------------------
775 EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
779 IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
780 AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES
785 IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
786 IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
787 W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
788 THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
789 ONE OF THE FOLLOWING VALUES ...
791 IERR = 1 IF A OR B IS NEGATIVE
792 IERR = 2 IF A = B = 0
793 IERR = 3 IF X .LT. 0 OR X .GT. 1
794 IERR = 4 IF Y .LT. 0 OR Y .GT. 1
795 IERR = 5 IF X + Y .NE. 1
796 IERR = 6 IF X = A = 0
797 IERR = 7 IF Y = B = 0
800 WRITTEN BY ALFRED H. MORRIS, JR.
801 NAVAL SURFACE WARFARE CENTER
804 -----------------------------------------------------------------------
808 static double a0,b0,eps,lambda,t,x0,y0,z;
809 static int ierr1,ind,n;
810 static double T2,T3,T4,T5;
813 .. Executable Statements ..
816 ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
817 FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
821 if(*a < 0.0e0 || *b < 0.0e0) goto S270;
822 if(*a == 0.0e0 && *b == 0.0e0) goto S280;
823 if(*x < 0.0e0 || *x > 1.0e0) goto S290;
824 if(*y < 0.0e0 || *y > 1.0e0) goto S300;
825 z = *x+*y-0.5e0-0.5e0;
826 if(fabs(z) > 3.0e0*eps) goto S310;
828 if(*x == 0.0e0) goto S210;
829 if(*y == 0.0e0) goto S230;
830 if(*a == 0.0e0) goto S240;
831 if(*b == 0.0e0) goto S220;
832 eps = fifdmax1(eps,1.e-15);
833 if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
839 if(fifdmin1(a0,b0) > 1.0e0) goto S40;
841 PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
843 if(*x <= 0.5e0) goto S10;
850 if(b0 < fifdmin1(eps,eps*a0)) goto S90;
851 if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
852 if(fifdmax1(a0,b0) > 1.0e0) goto S20;
853 if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
854 if(pow(x0,a0) <= 0.9e0) goto S110;
855 if(x0 >= 0.3e0) goto S120;
859 if(b0 <= 1.0e0) goto S110;
860 if(x0 >= 0.3e0) goto S120;
861 if(x0 >= 0.1e0) goto S30;
862 if(pow(x0*b0,a0) <= 0.7e0) goto S110;
864 if(b0 > 15.0e0) goto S150;
869 PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
871 if(*a > *b) goto S50;
872 lambda = *a-(*a+*b)**x;
875 lambda = (*a+*b)**y-*b;
877 if(lambda >= 0.0e0) goto S70;
883 lambda = fabs(lambda);
885 if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
886 if(b0 < 40.0e0) goto S160;
887 if(a0 > b0) goto S80;
888 if(a0 <= 100.0e0) goto S130;
889 if(lambda > 0.03e0*a0) goto S130;
892 if(b0 <= 100.0e0) goto S130;
893 if(lambda > 0.03e0*b0) goto S130;
897 EVALUATION OF THE APPROPRIATE ALGORITHM
899 *w = fpser(&a0,&b0,&x0,&eps);
900 *w1 = 0.5e0+(0.5e0-*w);
903 *w1 = apser(&a0,&b0,&x0,&eps);
904 *w = 0.5e0+(0.5e0-*w1);
907 *w = bpser(&a0,&b0,&x0,&eps);
908 *w1 = 0.5e0+(0.5e0-*w);
911 *w1 = bpser(&b0,&a0,&y0,&eps);
912 *w = 0.5e0+(0.5e0-*w1);
916 *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
917 *w1 = 0.5e0+(0.5e0-*w);
920 *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
924 bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
925 *w = 0.5e0+(0.5e0-*w1);
930 if(b0 != 0.0e0) goto S170;
934 *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
935 if(x0 > 0.7e0) goto S180;
936 *w += bpser(&a0,&b0,&x0,&eps);
937 *w1 = 0.5e0+(0.5e0-*w);
940 if(a0 > 15.0e0) goto S190;
942 *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
946 bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
947 *w1 = 0.5e0+(0.5e0-*w);
951 *w = basym(&a0,&b0,&lambda,&T5);
952 *w1 = 0.5e0+(0.5e0-*w);
956 TERMINATION OF THE PROCEDURE
958 if(*a == 0.0e0) goto S320;
964 if(*b == 0.0e0) goto S330;
977 PROCEDURE FOR A AND B .LT. 1.E-3*EPS
1007 double brcmp1(int *mu,double *a,double *b,double *x,double *y)
1009 -----------------------------------------------------------------------
1010 EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B))
1011 -----------------------------------------------------------------------
1014 static double Const = .398942280401433e0;
1015 static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
1019 CONST = 1/SQRT(2*PI)
1022 static double T1,T2,T3,T4;
1025 .. Executable Statements ..
1027 a0 = fifdmin1(*a,*b);
1028 if(a0 >= 8.0e0) goto S130;
1029 if(*x > 0.375e0) goto S10;
1035 if(*y > 0.375e0) goto S20;
1045 if(a0 < 1.0e0) goto S40;
1047 brcmp1 = esum(mu,&z);
1051 -----------------------------------------------------------------------
1052 PROCEDURE FOR A .LT. 1 OR B .LT. 1
1053 -----------------------------------------------------------------------
1055 b0 = fifdmax1(*a,*b);
1056 if(b0 >= 8.0e0) goto S120;
1057 if(b0 > 1.0e0) goto S70;
1059 ALGORITHM FOR B0 .LE. 1
1061 brcmp1 = esum(mu,&z);
1062 if(brcmp1 == 0.0e0) return brcmp1;
1064 if(apb > 1.0e0) goto S50;
1065 z = 1.0e0+gam1(&apb);
1069 z = (1.0e0+gam1(&u))/apb;
1071 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1072 brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
1076 ALGORITHM FOR 1 .LT. B0 .LT. 8
1079 n = (long)(b0 - 1.0e0);
1082 for(i=1; i<=n; i++) {
1091 if(apb > 1.0e0) goto S100;
1092 t = 1.0e0+gam1(&apb);
1096 t = (1.0e0+gam1(&u))/apb;
1098 brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
1102 ALGORITHM FOR B0 .GE. 8
1104 u = gamln1(&a0)+algdiv(&a0,&b0);
1106 brcmp1 = a0*esum(mu,&T3);
1110 -----------------------------------------------------------------------
1111 PROCEDURE FOR A .GE. 8 AND B .GE. 8
1112 -----------------------------------------------------------------------
1114 if(*a > *b) goto S140;
1117 y0 = 1.0e0/(1.0e0+h);
1118 lambda = *a-(*a+*b)**x;
1122 x0 = 1.0e0/(1.0e0+h);
1124 lambda = (*a+*b)**y-*b;
1127 if(fabs(e) > 0.6e0) goto S160;
1134 if(fabs(e) > 0.6e0) goto S180;
1142 brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
1145 double brcomp(double *a,double *b,double *x,double *y)
1147 -----------------------------------------------------------------------
1148 EVALUATION OF X**A*Y**B/BETA(A,B)
1149 -----------------------------------------------------------------------
1152 static double Const = .398942280401433e0;
1153 static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
1157 CONST = 1/SQRT(2*PI)
1160 static double T1,T2;
1163 .. Executable Statements ..
1166 if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
1167 a0 = fifdmin1(*a,*b);
1168 if(a0 >= 8.0e0) goto S130;
1169 if(*x > 0.375e0) goto S10;
1175 if(*y > 0.375e0) goto S20;
1185 if(a0 < 1.0e0) goto S40;
1191 -----------------------------------------------------------------------
1192 PROCEDURE FOR A .LT. 1 OR B .LT. 1
1193 -----------------------------------------------------------------------
1195 b0 = fifdmax1(*a,*b);
1196 if(b0 >= 8.0e0) goto S120;
1197 if(b0 > 1.0e0) goto S70;
1199 ALGORITHM FOR B0 .LE. 1
1202 if(brcomp == 0.0e0) return brcomp;
1204 if(apb > 1.0e0) goto S50;
1205 z = 1.0e0+gam1(&apb);
1209 z = (1.0e0+gam1(&u))/apb;
1211 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1212 brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
1216 ALGORITHM FOR 1 .LT. B0 .LT. 8
1219 n = (long)(b0 - 1.0e0);
1222 for(i=1; i<=n; i++) {
1231 if(apb > 1.0e0) goto S100;
1232 t = 1.0e0+gam1(&apb);
1236 t = (1.0e0+gam1(&u))/apb;
1238 brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
1242 ALGORITHM FOR B0 .GE. 8
1244 u = gamln1(&a0)+algdiv(&a0,&b0);
1245 brcomp = a0*exp(z-u);
1249 -----------------------------------------------------------------------
1250 PROCEDURE FOR A .GE. 8 AND B .GE. 8
1251 -----------------------------------------------------------------------
1253 if(*a > *b) goto S140;
1256 y0 = 1.0e0/(1.0e0+h);
1257 lambda = *a-(*a+*b)**x;
1261 x0 = 1.0e0/(1.0e0+h);
1263 lambda = (*a+*b)**y-*b;
1266 if(fabs(e) > 0.6e0) goto S160;
1273 if(fabs(e) > 0.6e0) goto S180;
1279 z = exp(-(*a*u+*b*v));
1280 brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
1283 double bup(double *a,double *b,double *x,double *y,int *n,double *eps)
1285 -----------------------------------------------------------------------
1286 EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
1287 EPS IS THE TOLERANCE USED.
1288 -----------------------------------------------------------------------
1293 static double bup,ap1,apb,d,l,r,t,w;
1294 static int i,k,kp1,mu,nm1;
1297 .. Executable Statements ..
1300 OBTAIN THE SCALING FACTOR EXP(-MU) AND
1301 EXP(MU)*(X**A*Y**B/BETA(A,B))/A
1307 if(*n == 1 || *a < 1.0e0) goto S10;
1308 if(apb < 1.1e0*ap1) goto S10;
1309 mu = (long)(fabs(exparg(&K1)));
1310 k = (long)(exparg(&K2));
1315 bup = brcmp1(&mu,a,b,x,y)/ *a;
1316 if(*n == 1 || bup == 0.0e0) return bup;
1320 LET K BE THE INDEX OF THE MAXIMUM TERM
1323 if(*b <= 1.0e0) goto S50;
1324 if(*y > 1.e-4) goto S20;
1328 r = (*b-1.0e0)**x/ *y-*a;
1329 if(r < 1.0e0) goto S50;
1332 if(r < t) k = (long)(r);
1335 ADD THE INCREASING TERMS OF THE SERIES
1337 for(i=1; i<=k; i++) {
1339 d = (apb+l)/(ap1+l)**x*d;
1342 if(k == nm1) goto S70;
1345 ADD THE REMAINING TERMS OF THE SERIES
1348 for(i=kp1; i<=nm1; i++) {
1350 d = (apb+l)/(ap1+l)**x*d;
1352 if(d <= *eps*w) goto S70;
1356 TERMINATE THE PROCEDURE
1361 void cdfbet(int *which,double *p,double *q,double *x,double *y,
1362 double *a,double *b,int *status,double *bound)
1363 /**********************************************************************
1365 void cdfbet(int *which,double *p,double *q,double *x,double *y,
1366 double *a,double *b,int *status,double *bound)
1368 Cumulative Distribution Function
1375 Calculates any one parameter of the beta distribution given
1376 values for the others.
1382 WHICH --> Integer indicating which of the next four argument
1383 values is to be calculated from the others.
1385 iwhich = 1 : Calculate P and Q from X,Y,A and B
1386 iwhich = 2 : Calculate X and Y from P,Q,A and B
1387 iwhich = 3 : Calculate A from P,Q,X,Y and B
1388 iwhich = 4 : Calculate B from P,Q,X,Y and A
1390 P <--> The integral from 0 to X of the chi-square
1392 Input range: [0, 1].
1395 Input range: [0, 1].
1398 X <--> Upper limit of integration of beta density.
1407 A <--> The first parameter of the beta density.
1408 Input range: (0, +infinity).
1409 Search range: [1D-100,1D100]
1411 B <--> The second parameter of the beta density.
1412 Input range: (0, +infinity).
1413 Search range: [1D-100,1D100]
1415 STATUS <-- 0 if calculation completed correctly
1416 -I if input parameter number I is out of range
1417 1 if answer appears to be lower than lowest
1419 2 if answer appears to be higher than greatest
1424 BOUND <-- Undefined if STATUS is 0
1426 Bound exceeded by parameter number I if STATUS
1429 Lower search bound if STATUS is 1.
1431 Upper search bound if STATUS is 2.
1437 Cumulative distribution function (P) is calculated directly by
1438 code associated with the following reference.
1440 DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant
1441 Digit Computation of the Incomplete Beta Function Ratios. ACM
1442 Trans. Math. Softw. 18 (1993), 360-373.
1444 Computation of other parameters involve a seach for a value that
1445 produces the desired value of P. The search relies on the
1446 monotinicity of P with the other parameter.
1452 The beta density is proportional to
1453 t^(A-1) * (1-t)^(B-1)
1455 **********************************************************************/
1458 #define atol 1.0e-50
1459 #define zero 1.0e-100
1463 static double K2 = 0.0e0;
1464 static double K3 = 1.0e0;
1465 static double K8 = 0.5e0;
1466 static double K9 = 5.0e0;
1467 static double fx,xhi,xlo,cum,ccum,xy,pq;
1468 static unsigned long qhi,qleft,qporq;
1469 static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
1472 .. Executable Statements ..
1477 if(!(*which < 1 || *which > 4)) goto S30;
1478 if(!(*which < 1)) goto S10;
1487 if(*which == 1) goto S70;
1491 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
1492 if(!(*p < 0.0e0)) goto S40;
1502 if(*which == 1) goto S110;
1506 if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
1507 if(!(*q < 0.0e0)) goto S80;
1517 if(*which == 2) goto S150;
1521 if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
1522 if(!(*x < 0.0e0)) goto S120;
1532 if(*which == 2) goto S190;
1536 if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
1537 if(!(*y < 0.0e0)) goto S160;
1547 if(*which == 3) goto S210;
1551 if(!(*a <= 0.0e0)) goto S200;
1557 if(*which == 4) goto S230;
1561 if(!(*b <= 0.0e0)) goto S220;
1567 if(*which == 1) goto S270;
1572 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
1573 if(!(pq < 0.0e0)) goto S240;
1583 if(*which == 2) goto S310;
1588 if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
1589 if(!(xy < 0.0e0)) goto S280;
1599 if(!(*which == 1)) qporq = *p <= *q;
1601 Select the minimum of P or Q
1608 cumbet(x,y,a,b,p,q);
1611 else if(2 == *which) {
1617 dstzr(&K2,&K3,&T4,&T5);
1618 if(!qporq) goto S340;
1620 dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
1623 if(!(*status == 1)) goto S330;
1624 cumbet(x,y,a,b,&cum,&ccum);
1626 dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
1633 dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
1636 if(!(*status == 1)) goto S360;
1637 cumbet(x,y,a,b,&cum,&ccum);
1639 dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
1644 if(!(*status == -1)) goto S400;
1645 if(!qleft) goto S380;
1656 else if(3 == *which) {
1665 dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
1667 dinvr(status,a,&fx,&qleft,&qhi);
1669 if(!(*status == 1)) goto S440;
1670 cumbet(x,y,a,b,&cum,&ccum);
1671 if(!qporq) goto S420;
1677 dinvr(status,a,&fx,&qleft,&qhi);
1680 if(!(*status == -1)) goto S470;
1681 if(!qleft) goto S450;
1692 else if(4 == *which) {
1701 dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
1703 dinvr(status,b,&fx,&qleft,&qhi);
1705 if(!(*status == 1)) goto S510;
1706 cumbet(x,y,a,b,&cum,&ccum);
1707 if(!qporq) goto S490;
1713 dinvr(status,b,&fx,&qleft,&qhi);
1716 if(!(*status == -1)) goto S540;
1717 if(!qleft) goto S520;
1735 void cdfbin(int *which,double *p,double *q,double *s,double *xn,
1736 double *pr,double *ompr,int *status,double *bound)
1737 /**********************************************************************
1739 void cdfbin(int *which,double *p,double *q,double *s,double *xn,
1740 double *pr,double *ompr,int *status,double *bound)
1742 Cumulative Distribution Function
1743 BINomial distribution
1749 Calculates any one parameter of the binomial
1750 distribution given values for the others.
1756 WHICH --> Integer indicating which of the next four argument
1757 values is to be calculated from the others.
1759 iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
1760 iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
1761 iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
1762 iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
1764 P <--> The cumulation from 0 to S of the binomial distribution.
1765 (Probablility of S or fewer successes in XN trials each
1766 with probability of success PR.)
1770 Input range: [0, 1].
1773 S <--> The number of successes observed.
1774 Input range: [0, XN]
1775 Search range: [0, XN]
1777 XN <--> The number of binomial trials.
1778 Input range: (0, +infinity).
1779 Search range: [1E-100, 1E100]
1781 PR <--> The probability of success in each binomial trial.
1790 STATUS <-- 0 if calculation completed correctly
1791 -I if input parameter number I is out of range
1792 1 if answer appears to be lower than lowest
1794 2 if answer appears to be higher than greatest
1797 4 if PR + OMPR .ne. 1
1799 BOUND <-- Undefined if STATUS is 0
1801 Bound exceeded by parameter number I if STATUS
1804 Lower search bound if STATUS is 1.
1806 Upper search bound if STATUS is 2.
1812 Formula 26.5.24 of Abramowitz and Stegun, Handbook of
1813 Mathematical Functions (1966) is used to reduce the binomial
1814 distribution to the cumulative incomplete beta distribution.
1816 Computation of other parameters involve a seach for a value that
1817 produces the desired value of P. The search relies on the
1818 monotinicity of P with the other parameter.
1821 **********************************************************************/
1823 #define atol 1.0e-50
1825 #define zero 1.0e-100
1829 static double K2 = 0.0e0;
1830 static double K3 = 0.5e0;
1831 static double K4 = 5.0e0;
1832 static double K11 = 1.0e0;
1833 static double fx,xhi,xlo,cum,ccum,pq,prompr;
1834 static unsigned long qhi,qleft,qporq;
1835 static double T5,T6,T7,T8,T9,T10,T12,T13;
1838 .. Executable Statements ..
1843 if(!(*which < 1 && *which > 4)) goto S30;
1844 if(!(*which < 1)) goto S10;
1853 if(*which == 1) goto S70;
1857 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
1858 if(!(*p < 0.0e0)) goto S40;
1868 if(*which == 1) goto S110;
1872 if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
1873 if(!(*q < 0.0e0)) goto S80;
1883 if(*which == 3) goto S130;
1887 if(!(*xn <= 0.0e0)) goto S120;
1893 if(*which == 2) goto S170;
1897 if(!(*s < 0.0e0 || (*which != 3 && *s > *xn))) goto S160;
1898 if(!(*s < 0.0e0)) goto S140;
1908 if(*which == 4) goto S210;
1912 if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S200;
1913 if(!(*pr < 0.0e0)) goto S180;
1923 if(*which == 4) goto S250;
1927 if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S240;
1928 if(!(*ompr < 0.0e0)) goto S220;
1938 if(*which == 1) goto S290;
1943 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S280;
1944 if(!(pq < 0.0e0)) goto S260;
1954 if(*which == 4) goto S330;
1959 if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S320;
1960 if(!(prompr < 0.0e0)) goto S300;
1970 if(!(*which == 1)) qporq = *p <= *q;
1972 Select the minimum of P or Q
1979 cumbin(s,xn,pr,ompr,p,q);
1982 else if(2 == *which) {
1989 dstinv(&K2,xn,&K3,&K3,&K4,&T5,&T6);
1991 dinvr(status,s,&fx,&qleft,&qhi);
1993 if(!(*status == 1)) goto S370;
1994 cumbin(s,xn,pr,ompr,&cum,&ccum);
1995 if(!qporq) goto S350;
2001 dinvr(status,s,&fx,&qleft,&qhi);
2004 if(!(*status == -1)) goto S400;
2005 if(!qleft) goto S380;
2016 else if(3 == *which) {
2025 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
2027 dinvr(status,xn,&fx,&qleft,&qhi);
2029 if(!(*status == 1)) goto S440;
2030 cumbin(s,xn,pr,ompr,&cum,&ccum);
2031 if(!qporq) goto S420;
2037 dinvr(status,xn,&fx,&qleft,&qhi);
2040 if(!(*status == -1)) goto S470;
2041 if(!qleft) goto S450;
2052 else if(4 == *which) {
2054 Calculating PR and OMPR
2058 dstzr(&K2,&K11,&T12,&T13);
2059 if(!qporq) goto S500;
2061 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
2064 if(!(*status == 1)) goto S490;
2065 cumbin(s,xn,pr,ompr,&cum,&ccum);
2067 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
2074 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
2077 if(!(*status == 1)) goto S520;
2078 cumbin(s,xn,pr,ompr,&cum,&ccum);
2080 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
2085 if(!(*status == -1)) goto S560;
2086 if(!qleft) goto S540;
2104 void cdfchi(int *which,double *p,double *q,double *x,double *df,
2105 int *status,double *bound)
2106 /**********************************************************************
2108 void cdfchi(int *which,double *p,double *q,double *x,double *df,
2109 int *status,double *bound)
2111 Cumulative Distribution Function
2112 CHI-Square distribution
2118 Calculates any one parameter of the chi-square
2119 distribution given values for the others.
2125 WHICH --> Integer indicating which of the next three argument
2126 values is to be calculated from the others.
2128 iwhich = 1 : Calculate P and Q from X and DF
2129 iwhich = 2 : Calculate X from P,Q and DF
2130 iwhich = 3 : Calculate DF from P,Q and X
2132 P <--> The integral from 0 to X of the chi-square
2134 Input range: [0, 1].
2137 Input range: (0, 1].
2140 X <--> Upper limit of integration of the non-central
2141 chi-square distribution.
2142 Input range: [0, +infinity).
2143 Search range: [0,1E100]
2145 DF <--> Degrees of freedom of the
2146 chi-square distribution.
2147 Input range: (0, +infinity).
2148 Search range: [ 1E-100, 1E100]
2150 STATUS <-- 0 if calculation completed correctly
2151 -I if input parameter number I is out of range
2152 1 if answer appears to be lower than lowest
2154 2 if answer appears to be higher than greatest
2157 10 indicates error returned from cumgam. See
2158 references in cdfgam
2160 BOUND <-- Undefined if STATUS is 0
2162 Bound exceeded by parameter number I if STATUS
2165 Lower search bound if STATUS is 1.
2167 Upper search bound if STATUS is 2.
2173 Formula 26.4.19 of Abramowitz and Stegun, Handbook of
2174 Mathematical Functions (1966) is used to reduce the chisqure
2175 distribution to the incomplete distribution.
2177 Computation of other parameters involve a seach for a value that
2178 produces the desired value of P. The search relies on the
2179 monotinicity of P with the other parameter.
2181 **********************************************************************/
2184 #define atol 1.0e-50
2185 #define zero 1.0e-100
2188 static double K2 = 0.0e0;
2189 static double K4 = 0.5e0;
2190 static double K5 = 5.0e0;
2191 static double fx,cum,ccum,pq,porq;
2192 static unsigned long qhi,qleft,qporq;
2193 static double T3,T6,T7,T8,T9,T10,T11;
2196 .. Executable Statements ..
2201 if(!(*which < 1 || *which > 3)) goto S30;
2202 if(!(*which < 1)) goto S10;
2211 if(*which == 1) goto S70;
2215 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
2216 if(!(*p < 0.0e0)) goto S40;
2226 if(*which == 1) goto S110;
2230 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
2231 if(!(*q <= 0.0e0)) goto S80;
2241 if(*which == 2) goto S130;
2245 if(!(*x < 0.0e0)) goto S120;
2251 if(*which == 3) goto S150;
2255 if(!(*df <= 0.0e0)) goto S140;
2261 if(*which == 1) goto S190;
2266 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
2267 if(!(pq < 0.0e0)) goto S160;
2277 if(*which == 1) goto S220;
2279 Select the minimum of P or Q
2282 if(!qporq) goto S200;
2303 else if(2 == *which) {
2311 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
2313 dinvr(status,x,&fx,&qleft,&qhi);
2315 if(!(*status == 1)) goto S270;
2316 cumchi(x,df,&cum,&ccum);
2317 if(!qporq) goto S240;
2323 if(!(fx+porq > 1.5e0)) goto S260;
2327 dinvr(status,x,&fx,&qleft,&qhi);
2330 if(!(*status == -1)) goto S300;
2331 if(!qleft) goto S280;
2342 else if(3 == *which) {
2351 dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
2353 dinvr(status,df,&fx,&qleft,&qhi);
2355 if(!(*status == 1)) goto S350;
2356 cumchi(x,df,&cum,&ccum);
2357 if(!qporq) goto S320;
2363 if(!(fx+porq > 1.5e0)) goto S340;
2367 dinvr(status,df,&fx,&qleft,&qhi);
2370 if(!(*status == -1)) goto S380;
2371 if(!qleft) goto S360;
2388 void cdfchn(int *which,double *p,double *q,double *x,double *df,
2389 double *pnonc,int *status,double *bound)
2390 /**********************************************************************
2392 void cdfchn(int *which,double *p,double *q,double *x,double *df,
2393 double *pnonc,int *status,double *bound)
2395 Cumulative Distribution Function
2396 Non-central Chi-Square
2402 Calculates any one parameter of the non-central chi-square
2403 distribution given values for the others.
2409 WHICH --> Integer indicating which of the next three argument
2410 values is to be calculated from the others.
2412 iwhich = 1 : Calculate P and Q from X and DF
2413 iwhich = 2 : Calculate X from P,DF and PNONC
2414 iwhich = 3 : Calculate DF from P,X and PNONC
2415 iwhich = 3 : Calculate PNONC from P,X and DF
2417 P <--> The integral from 0 to X of the non-central chi-square
2419 Input range: [0, 1-1E-16).
2422 Q is not used by this subroutine and is only included
2423 for similarity with other cdf* routines.
2425 X <--> Upper limit of integration of the non-central
2426 chi-square distribution.
2427 Input range: [0, +infinity).
2428 Search range: [0,1E100]
2430 DF <--> Degrees of freedom of the non-central
2431 chi-square distribution.
2432 Input range: (0, +infinity).
2433 Search range: [ 1E-100, 1E100]
2435 PNONC <--> Non-centrality parameter of the non-central
2436 chi-square distribution.
2437 Input range: [0, +infinity).
2438 Search range: [0,1E4]
2440 STATUS <-- 0 if calculation completed correctly
2441 -I if input parameter number I is out of range
2442 1 if answer appears to be lower than lowest
2444 2 if answer appears to be higher than greatest
2447 BOUND <-- Undefined if STATUS is 0
2449 Bound exceeded by parameter number I if STATUS
2452 Lower search bound if STATUS is 1.
2454 Upper search bound if STATUS is 2.
2460 Formula 26.4.25 of Abramowitz and Stegun, Handbook of
2461 Mathematical Functions (1966) is used to compute the cumulative
2462 distribution function.
2464 Computation of other parameters involve a seach for a value that
2465 produces the desired value of P. The search relies on the
2466 monotinicity of P with the other parameter.
2471 The computation time required for this routine is proportional
2472 to the noncentrality parameter (PNONC). Very large values of
2473 this parameter can consume immense computer resources. This is
2474 why the search range is bounded by 10,000.
2476 **********************************************************************/
2480 #define atol 1.0e-50
2481 #define zero 1.0e-100
2482 #define one ( 1.0e0 - 1.0e-16 )
2484 static double K1 = 0.0e0;
2485 static double K3 = 0.5e0;
2486 static double K4 = 5.0e0;
2487 static double fx,cum,ccum;
2488 static unsigned long qhi,qleft;
2489 static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;
2492 .. Executable Statements ..
2497 if(!(*which < 1 || *which > 4)) goto S30;
2498 if(!(*which < 1)) goto S10;
2507 if(*which == 1) goto S70;
2511 if(!(*p < 0.0e0 || *p > one)) goto S60;
2512 if(!(*p < 0.0e0)) goto S40;
2522 if(*which == 2) goto S90;
2526 if(!(*x < 0.0e0)) goto S80;
2532 if(*which == 3) goto S110;
2536 if(!(*df <= 0.0e0)) goto S100;
2542 if(*which == 4) goto S130;
2546 if(!(*pnonc < 0.0e0)) goto S120;
2559 cumchn(x,df,pnonc,p,q);
2562 else if(2 == *which) {
2570 dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
2572 dinvr(status,x,&fx,&qleft,&qhi);
2574 if(!(*status == 1)) goto S150;
2575 cumchn(x,df,pnonc,&cum,&ccum);
2577 dinvr(status,x,&fx,&qleft,&qhi);
2580 if(!(*status == -1)) goto S180;
2581 if(!qleft) goto S160;
2592 else if(3 == *which) {
2601 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
2603 dinvr(status,df,&fx,&qleft,&qhi);
2605 if(!(*status == 1)) goto S200;
2606 cumchn(x,df,pnonc,&cum,&ccum);
2608 dinvr(status,df,&fx,&qleft,&qhi);
2611 if(!(*status == -1)) goto S230;
2612 if(!qleft) goto S210;
2623 else if(4 == *which) {
2631 dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);
2633 dinvr(status,pnonc,&fx,&qleft,&qhi);
2635 if(!(*status == 1)) goto S250;
2636 cumchn(x,df,pnonc,&cum,&ccum);
2638 dinvr(status,pnonc,&fx,&qleft,&qhi);
2641 if(!(*status == -1)) goto S280;
2642 if(!qleft) goto S260;
2661 void cdff(int *which,double *p,double *q,double *f,double *dfn,
2662 double *dfd,int *status,double *bound)
2663 /**********************************************************************
2665 void cdff(int *which,double *p,double *q,double *f,double *dfn,
2666 double *dfd,int *status,double *bound)
2668 Cumulative Distribution Function
2675 Calculates any one parameter of the F distribution
2676 given values for the others.
2682 WHICH --> Integer indicating which of the next four argument
2683 values is to be calculated from the others.
2685 iwhich = 1 : Calculate P and Q from F,DFN and DFD
2686 iwhich = 2 : Calculate F from P,Q,DFN and DFD
2687 iwhich = 3 : Calculate DFN from P,Q,F and DFD
2688 iwhich = 4 : Calculate DFD from P,Q,F and DFN
2690 P <--> The integral from 0 to F of the f-density.
2694 Input range: (0, 1].
2697 F <--> Upper limit of integration of the f-density.
2698 Input range: [0, +infinity).
2699 Search range: [0,1E100]
2701 DFN < --> Degrees of freedom of the numerator sum of squares.
2702 Input range: (0, +infinity).
2703 Search range: [ 1E-100, 1E100]
2705 DFD < --> Degrees of freedom of the denominator sum of squares.
2706 Input range: (0, +infinity).
2707 Search range: [ 1E-100, 1E100]
2709 STATUS <-- 0 if calculation completed correctly
2710 -I if input parameter number I is out of range
2711 1 if answer appears to be lower than lowest
2713 2 if answer appears to be higher than greatest
2717 BOUND <-- Undefined if STATUS is 0
2719 Bound exceeded by parameter number I if STATUS
2722 Lower search bound if STATUS is 1.
2724 Upper search bound if STATUS is 2.
2730 Formula 26.6.2 of Abramowitz and Stegun, Handbook of
2731 Mathematical Functions (1966) is used to reduce the computation
2732 of the cumulative distribution function for the F variate to
2733 that of an incomplete beta.
2735 Computation of other parameters involve a seach for a value that
2736 produces the desired value of P. The search relies on the
2737 monotinicity of P with the other parameter.
2741 The value of the cumulative F distribution is not necessarily
2742 monotone in either degrees of freedom. There thus may be two
2743 values that provide a given CDF value. This routine assumes
2744 monotonicity and will find an arbitrary one of the two values.
2746 **********************************************************************/
2749 #define atol 1.0e-50
2750 #define zero 1.0e-100
2753 static double K2 = 0.0e0;
2754 static double K4 = 0.5e0;
2755 static double K5 = 5.0e0;
2756 static double pq,fx,cum,ccum;
2757 static unsigned long qhi,qleft,qporq;
2758 static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
2761 .. Executable Statements ..
2766 if(!(*which < 1 || *which > 4)) goto S30;
2767 if(!(*which < 1)) goto S10;
2776 if(*which == 1) goto S70;
2780 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
2781 if(!(*p < 0.0e0)) goto S40;
2791 if(*which == 1) goto S110;
2795 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
2796 if(!(*q <= 0.0e0)) goto S80;
2806 if(*which == 2) goto S130;
2810 if(!(*f < 0.0e0)) goto S120;
2816 if(*which == 3) goto S150;
2820 if(!(*dfn <= 0.0e0)) goto S140;
2826 if(*which == 4) goto S170;
2830 if(!(*dfd <= 0.0e0)) goto S160;
2836 if(*which == 1) goto S210;
2841 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
2842 if(!(pq < 0.0e0)) goto S180;
2852 if(!(*which == 1)) qporq = *p <= *q;
2854 Select the minimum of P or Q
2861 cumf(f,dfn,dfd,p,q);
2864 else if(2 == *which) {
2872 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
2874 dinvr(status,f,&fx,&qleft,&qhi);
2876 if(!(*status == 1)) goto S250;
2877 cumf(f,dfn,dfd,&cum,&ccum);
2878 if(!qporq) goto S230;
2884 dinvr(status,f,&fx,&qleft,&qhi);
2887 if(!(*status == -1)) goto S280;
2888 if(!qleft) goto S260;
2899 else if(3 == *which) {
2908 dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
2910 dinvr(status,dfn,&fx,&qleft,&qhi);
2912 if(!(*status == 1)) goto S320;
2913 cumf(f,dfn,dfd,&cum,&ccum);
2914 if(!qporq) goto S300;
2920 dinvr(status,dfn,&fx,&qleft,&qhi);
2923 if(!(*status == -1)) goto S350;
2924 if(!qleft) goto S330;
2935 else if(4 == *which) {
2944 dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15);
2946 dinvr(status,dfd,&fx,&qleft,&qhi);
2948 if(!(*status == 1)) goto S390;
2949 cumf(f,dfn,dfd,&cum,&ccum);
2950 if(!qporq) goto S370;
2956 dinvr(status,dfd,&fx,&qleft,&qhi);
2959 if(!(*status == -1)) goto S420;
2960 if(!qleft) goto S400;
2977 void cdffnc(int *which,double *p,double *q,double *f,double *dfn,
2978 double *dfd,double *phonc,int *status,double *bound)
2979 /**********************************************************************
2981 void cdffnc(int *which,double *p,double *q,double *f,double *dfn,
2982 double *dfd,double *phonc,int *status,double *bound)
2984 Cumulative Distribution Function
2985 Non-central F distribution
2991 Calculates any one parameter of the Non-central F
2992 distribution given values for the others.
2998 WHICH --> Integer indicating which of the next five argument
2999 values is to be calculated from the others.
3001 iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
3002 iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
3003 iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
3004 iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
3005 iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
3007 P <--> The integral from 0 to F of the non-central f-density.
3008 Input range: [0,1-1E-16).
3011 Q is not used by this subroutine and is only included
3012 for similarity with other cdf* routines.
3014 F <--> Upper limit of integration of the non-central f-density.
3015 Input range: [0, +infinity).
3016 Search range: [0,1E100]
3018 DFN < --> Degrees of freedom of the numerator sum of squares.
3019 Input range: (0, +infinity).
3020 Search range: [ 1E-100, 1E100]
3022 DFD < --> Degrees of freedom of the denominator sum of squares.
3023 Must be in range: (0, +infinity).
3024 Input range: (0, +infinity).
3025 Search range: [ 1E-100, 1E100]
3027 PNONC <-> The non-centrality parameter
3028 Input range: [0,infinity)
3029 Search range: [0,1E4]
3031 STATUS <-- 0 if calculation completed correctly
3032 -I if input parameter number I is out of range
3033 1 if answer appears to be lower than lowest
3035 2 if answer appears to be higher than greatest
3039 BOUND <-- Undefined if STATUS is 0
3041 Bound exceeded by parameter number I if STATUS
3044 Lower search bound if STATUS is 1.
3046 Upper search bound if STATUS is 2.
3052 Formula 26.6.20 of Abramowitz and Stegun, Handbook of
3053 Mathematical Functions (1966) is used to compute the cumulative
3054 distribution function.
3056 Computation of other parameters involve a seach for a value that
3057 produces the desired value of P. The search relies on the
3058 monotinicity of P with the other parameter.
3062 The computation time required for this routine is proportional
3063 to the noncentrality parameter (PNONC). Very large values of
3064 this parameter can consume immense computer resources. This is
3065 why the search range is bounded by 10,000.
3069 The value of the cumulative noncentral F distribution is not
3070 necessarily monotone in either degrees of freedom. There thus
3071 may be two values that provide a given CDF value. This routine
3072 assumes monotonicity and will find an arbitrary one of the two
3075 **********************************************************************/
3079 #define atol 1.0e-50
3080 #define zero 1.0e-100
3081 #define one ( 1.0e0 - 1.0e-16 )
3083 static double K1 = 0.0e0;
3084 static double K3 = 0.5e0;
3085 static double K4 = 5.0e0;
3086 static double fx,cum,ccum;
3087 static unsigned long qhi,qleft;
3088 static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17;
3091 .. Executable Statements ..
3096 if(!(*which < 1 || *which > 5)) goto S30;
3097 if(!(*which < 1)) goto S10;
3106 if(*which == 1) goto S70;
3110 if(!(*p < 0.0e0 || *p > one)) goto S60;
3111 if(!(*p < 0.0e0)) goto S40;
3121 if(*which == 2) goto S90;
3125 if(!(*f < 0.0e0)) goto S80;
3131 if(*which == 3) goto S110;
3135 if(!(*dfn <= 0.0e0)) goto S100;
3141 if(*which == 4) goto S130;
3145 if(!(*dfd <= 0.0e0)) goto S120;
3151 if(*which == 5) goto S150;
3155 if(!(*phonc < 0.0e0)) goto S140;
3168 cumfnc(f,dfn,dfd,phonc,p,q);
3171 else if(2 == *which) {
3179 dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
3181 dinvr(status,f,&fx,&qleft,&qhi);
3183 if(!(*status == 1)) goto S170;
3184 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3186 dinvr(status,f,&fx,&qleft,&qhi);
3189 if(!(*status == -1)) goto S200;
3190 if(!qleft) goto S180;
3201 else if(3 == *which) {
3210 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
3212 dinvr(status,dfn,&fx,&qleft,&qhi);
3214 if(!(*status == 1)) goto S220;
3215 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3217 dinvr(status,dfn,&fx,&qleft,&qhi);
3220 if(!(*status == -1)) goto S250;
3221 if(!qleft) goto S230;
3232 else if(4 == *which) {
3241 dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
3243 dinvr(status,dfd,&fx,&qleft,&qhi);
3245 if(!(*status == 1)) goto S270;
3246 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3248 dinvr(status,dfd,&fx,&qleft,&qhi);
3251 if(!(*status == -1)) goto S300;
3252 if(!qleft) goto S280;
3263 else if(5 == *which) {
3271 dstinv(&K1,&T15,&K3,&K3,&K4,&T16,&T17);
3273 dinvr(status,phonc,&fx,&qleft,&qhi);
3275 if(!(*status == 1)) goto S320;
3276 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3278 dinvr(status,phonc,&fx,&qleft,&qhi);
3281 if(!(*status == -1)) goto S350;
3282 if(!qleft) goto S330;
3301 void cdfgam(int *which,double *p,double *q,double *x,double *shape,
3302 double *scale,int *status,double *bound)
3303 /**********************************************************************
3305 void cdfgam(int *which,double *p,double *q,double *x,double *shape,
3306 double *scale,int *status,double *bound)
3308 Cumulative Distribution Function
3315 Calculates any one parameter of the gamma
3316 distribution given values for the others.
3322 WHICH --> Integer indicating which of the next four argument
3323 values is to be calculated from the others.
3325 iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE
3326 iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE
3327 iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE
3328 iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE
3330 P <--> The integral from 0 to X of the gamma density.
3334 Input range: (0, 1].
3337 X <--> The upper limit of integration of the gamma density.
3338 Input range: [0, +infinity).
3339 Search range: [0,1E100]
3341 SHAPE <--> The shape parameter of the gamma density.
3342 Input range: (0, +infinity).
3343 Search range: [1E-100,1E100]
3345 SCALE <--> The scale parameter of the gamma density.
3346 Input range: (0, +infinity).
3347 Search range: (1E-100,1E100]
3349 STATUS <-- 0 if calculation completed correctly
3350 -I if input parameter number I is out of range
3351 1 if answer appears to be lower than lowest
3353 2 if answer appears to be higher than greatest
3356 10 if the gamma or inverse gamma routine cannot
3357 compute the answer. Usually happens only for
3358 X and SHAPE very large (gt 1E10 or more)
3360 BOUND <-- Undefined if STATUS is 0
3362 Bound exceeded by parameter number I if STATUS
3365 Lower search bound if STATUS is 1.
3367 Upper search bound if STATUS is 2.
3373 Cumulative distribution function (P) is calculated directly by
3374 the code associated with:
3376 DiDinato, A. R. and Morris, A. H. Computation of the incomplete
3377 gamma function ratios and their inverse. ACM Trans. Math.
3378 Softw. 12 (1986), 377-393.
3380 Computation of other parameters involve a seach for a value that
3381 produces the desired value of P. The search relies on the
3382 monotinicity of P with the other parameter.
3389 The gamma density is proportional to
3390 T**(SHAPE - 1) * EXP(- SCALE * T)
3392 **********************************************************************/
3395 #define atol 1.0e-50
3396 #define zero 1.0e-100
3399 static double K5 = 0.5e0;
3400 static double K6 = 5.0e0;
3401 static double xx,fx,xscale,cum,ccum,pq,porq;
3403 static unsigned long qhi,qleft,qporq;
3404 static double T2,T3,T4,T7,T8,T9;
3407 .. Executable Statements ..
3412 if(!(*which < 1 || *which > 4)) goto S30;
3413 if(!(*which < 1)) goto S10;
3422 if(*which == 1) goto S70;
3426 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
3427 if(!(*p < 0.0e0)) goto S40;
3437 if(*which == 1) goto S110;
3441 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
3442 if(!(*q <= 0.0e0)) goto S80;
3452 if(*which == 2) goto S130;
3456 if(!(*x < 0.0e0)) goto S120;
3462 if(*which == 3) goto S150;
3466 if(!(*shape <= 0.0e0)) goto S140;
3472 if(*which == 4) goto S170;
3476 if(!(*scale <= 0.0e0)) goto S160;
3482 if(*which == 1) goto S210;
3487 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
3488 if(!(pq < 0.0e0)) goto S180;
3498 if(*which == 1) goto S240;
3500 Select the minimum of P or Q
3503 if(!qporq) goto S220;
3519 cumgam(&xscale,shape,p,q);
3520 if(porq > 1.5e0) *status = 10;
3522 else if(2 == *which) {
3527 gaminv(shape,&xx,&T2,p,q,&ierr);
3537 else if(3 == *which) {
3547 dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
3549 dinvr(status,shape,&fx,&qleft,&qhi);
3551 if(!(*status == 1)) goto S290;
3552 cumgam(&xscale,shape,&cum,&ccum);
3553 if(!qporq) goto S260;
3559 if(!((qporq && cum > 1.5e0) || (!qporq && ccum > 1.5e0))) goto S280;
3563 dinvr(status,shape,&fx,&qleft,&qhi);
3566 if(!(*status == -1)) goto S320;
3567 if(!qleft) goto S300;
3578 else if(4 == *which) {
3583 gaminv(shape,&xx,&T9,p,q,&ierr);
3599 void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
3600 double *pr,double *ompr,int *status,double *bound)
3601 /**********************************************************************
3603 void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
3604 double *pr,double *ompr,int *status,double *bound)
3606 Cumulative Distribution Function
3607 Negative BiNomial distribution
3613 Calculates any one parameter of the negative binomial
3614 distribution given values for the others.
3616 The cumulative negative binomial distribution returns the
3617 probability that there will be F or fewer failures before the
3618 XNth success in binomial trials each of which has probability of
3621 The individual term of the negative binomial is the probability of
3622 S failures before XN successes and is
3623 Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S
3629 WHICH --> Integer indicating which of the next four argument
3630 values is to be calculated from the others.
3632 iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
3633 iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
3634 iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
3635 iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
3637 P <--> The cumulation from 0 to S of the negative
3638 binomial distribution.
3642 Input range: (0, 1].
3645 S <--> The upper limit of cumulation of the binomial distribution.
3646 There are F or fewer failures before the XNth success.
3647 Input range: [0, +infinity).
3648 Search range: [0, 1E100]
3650 XN <--> The number of successes.
3651 Input range: [0, +infinity).
3652 Search range: [0, 1E100]
3654 PR <--> The probability of success in each binomial trial.
3656 Search range: [0,1].
3663 STATUS <-- 0 if calculation completed correctly
3664 -I if input parameter number I is out of range
3665 1 if answer appears to be lower than lowest
3667 2 if answer appears to be higher than greatest
3670 4 if PR + OMPR .ne. 1
3672 BOUND <-- Undefined if STATUS is 0
3674 Bound exceeded by parameter number I if STATUS
3677 Lower search bound if STATUS is 1.
3679 Upper search bound if STATUS is 2.
3685 Formula 26.5.26 of Abramowitz and Stegun, Handbook of
3686 Mathematical Functions (1966) is used to reduce calculation of
3687 the cumulative distribution function to that of an incomplete
3690 Computation of other parameters involve a seach for a value that
3691 produces the desired value of P. The search relies on the
3692 monotinicity of P with the other parameter.
3694 **********************************************************************/
3697 #define atol 1.0e-50
3701 static double K2 = 0.0e0;
3702 static double K4 = 0.5e0;
3703 static double K5 = 5.0e0;
3704 static double K11 = 1.0e0;
3705 static double fx,xhi,xlo,pq,prompr,cum,ccum;
3706 static unsigned long qhi,qleft,qporq;
3707 static double T3,T6,T7,T8,T9,T10,T12,T13;
3710 .. Executable Statements ..
3715 if(!(*which < 1 || *which > 4)) goto S30;
3716 if(!(*which < 1)) goto S10;
3725 if(*which == 1) goto S70;
3729 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
3730 if(!(*p < 0.0e0)) goto S40;
3740 if(*which == 1) goto S110;
3744 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
3745 if(!(*q <= 0.0e0)) goto S80;
3755 if(*which == 2) goto S130;
3759 if(!(*s < 0.0e0)) goto S120;
3765 if(*which == 3) goto S150;
3769 if(!(*xn < 0.0e0)) goto S140;
3775 if(*which == 4) goto S190;
3779 if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
3780 if(!(*pr < 0.0e0)) goto S160;
3790 if(*which == 4) goto S230;
3794 if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
3795 if(!(*ompr < 0.0e0)) goto S200;
3805 if(*which == 1) goto S270;
3810 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
3811 if(!(pq < 0.0e0)) goto S240;
3821 if(*which == 4) goto S310;
3826 if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
3827 if(!(prompr < 0.0e0)) goto S280;
3837 if(!(*which == 1)) qporq = *p <= *q;
3839 Select the minimum of P or Q
3846 cumnbn(s,xn,pr,ompr,p,q);
3849 else if(2 == *which) {
3857 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
3859 dinvr(status,s,&fx,&qleft,&qhi);
3861 if(!(*status == 1)) goto S350;
3862 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3863 if(!qporq) goto S330;
3869 dinvr(status,s,&fx,&qleft,&qhi);
3872 if(!(*status == -1)) goto S380;
3873 if(!qleft) goto S360;
3884 else if(3 == *which) {
3892 dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
3894 dinvr(status,xn,&fx,&qleft,&qhi);
3896 if(!(*status == 1)) goto S420;
3897 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3898 if(!qporq) goto S400;
3904 dinvr(status,xn,&fx,&qleft,&qhi);
3907 if(!(*status == -1)) goto S450;
3908 if(!qleft) goto S430;
3919 else if(4 == *which) {
3921 Calculating PR and OMPR
3925 dstzr(&K2,&K11,&T12,&T13);
3926 if(!qporq) goto S480;
3928 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
3931 if(!(*status == 1)) goto S470;
3932 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3934 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
3941 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
3944 if(!(*status == 1)) goto S500;
3945 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3947 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
3952 if(!(*status == -1)) goto S540;
3953 if(!qleft) goto S520;
3970 void cdfnor(int *which,double *p,double *q,double *x,double *mean,
3971 double *sd,int *status,double *bound)
3972 /**********************************************************************
3974 void cdfnor(int *which,double *p,double *q,double *x,double *mean,
3975 double *sd,int *status,double *bound)
3977 Cumulative Distribution Function
3984 Calculates any one parameter of the normal
3985 distribution given values for the others.
3991 WHICH --> Integer indicating which of the next parameter
3992 values is to be calculated using values of the others.
3994 iwhich = 1 : Calculate P and Q from X,MEAN and SD
3995 iwhich = 2 : Calculate X from P,Q,MEAN and SD
3996 iwhich = 3 : Calculate MEAN from P,Q,X and SD
3997 iwhich = 4 : Calculate SD from P,Q,X and MEAN
3999 P <--> The integral from -infinity to X of the normal density.
4003 Input range: (0, 1].
4006 X < --> Upper limit of integration of the normal-density.
4007 Input range: ( -infinity, +infinity)
4009 MEAN <--> The mean of the normal density.
4010 Input range: (-infinity, +infinity)
4012 SD <--> Standard Deviation of the normal density.
4013 Input range: (0, +infinity).
4015 STATUS <-- 0 if calculation completed correctly
4016 -I if input parameter number I is out of range
4017 1 if answer appears to be lower than lowest
4019 2 if answer appears to be higher than greatest
4023 BOUND <-- Undefined if STATUS is 0
4025 Bound exceeded by parameter number I if STATUS
4028 Lower search bound if STATUS is 1.
4030 Upper search bound if STATUS is 2.
4038 A slightly modified version of ANORM from
4040 Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
4041 Package of Special Function Routines and Test Drivers"
4042 acm Transactions on Mathematical Software. 19, 22-32.
4044 is used to calulate the cumulative standard normal distribution.
4046 The rational functions from pages 90-95 of Kennedy and Gentle,
4047 Statistical Computing, Marcel Dekker, NY, 1980 are used as
4048 starting values to Newton's Iterations which compute the inverse
4049 standard normal. Therefore no searches are necessary for any
4052 For X < -15, the asymptotic expansion for the normal is used as
4053 the starting value in finding the inverse standard normal.
4054 This is formula 26.2.12 of Abramowitz and Stegun.
4060 The normal density is proportional to
4061 exp( - 0.5 * (( X - MEAN)/SD)**2)
4063 **********************************************************************/
4069 .. Executable Statements ..
4075 if(!(*which < 1 || *which > 4)) goto S30;
4076 if(!(*which < 1)) goto S10;
4085 if(*which == 1) goto S70;
4089 if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
4090 if(!(*p <= 0.0e0)) goto S40;
4100 if(*which == 1) goto S110;
4104 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4105 if(!(*q <= 0.0e0)) goto S80;
4115 if(*which == 1) goto S150;
4120 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S140;
4121 if(!(pq < 0.0e0)) goto S120;
4131 if(*which == 4) goto S170;
4135 if(!(*sd <= 0.0e0)) goto S160;
4148 z = (*x-*mean)/ *sd;
4151 else if(2 == *which) {
4158 else if(3 == *which) {
4165 else if(4 == *which) {
4174 void cdfpoi(int *which,double *p,double *q,double *s,double *xlam,
4175 int *status,double *bound)
4176 /**********************************************************************
4178 void cdfpoi(int *which,double *p,double *q,double *s,double *xlam,
4179 int *status,double *bound)
4181 Cumulative Distribution Function
4182 POIsson distribution
4188 Calculates any one parameter of the Poisson
4189 distribution given values for the others.
4195 WHICH --> Integer indicating which argument
4196 value is to be calculated from the others.
4198 iwhich = 1 : Calculate P and Q from S and XLAM
4199 iwhich = 2 : Calculate A from P,Q and XLAM
4200 iwhich = 3 : Calculate XLAM from P,Q and S
4202 P <--> The cumulation from 0 to S of the poisson density.
4206 Input range: (0, 1].
4209 S <--> Upper limit of cumulation of the Poisson.
4210 Input range: [0, +infinity).
4211 Search range: [0,1E100]
4213 XLAM <--> Mean of the Poisson distribution.
4214 Input range: [0, +infinity).
4215 Search range: [0,1E100]
4217 STATUS <-- 0 if calculation completed correctly
4218 -I if input parameter number I is out of range
4219 1 if answer appears to be lower than lowest
4221 2 if answer appears to be higher than greatest
4225 BOUND <-- Undefined if STATUS is 0
4227 Bound exceeded by parameter number I if STATUS
4230 Lower search bound if STATUS is 1.
4232 Upper search bound if STATUS is 2.
4238 Formula 26.4.21 of Abramowitz and Stegun, Handbook of
4239 Mathematical Functions (1966) is used to reduce the computation
4240 of the cumulative distribution function to that of computing a
4241 chi-square, hence an incomplete gamma function.
4243 Cumulative distribution function (P) is calculated directly.
4244 Computation of other parameters involve a seach for a value that
4245 produces the desired value of P. The search relies on the
4246 monotinicity of P with the other parameter.
4248 **********************************************************************/
4251 #define atol 1.0e-50
4254 static double K2 = 0.0e0;
4255 static double K4 = 0.5e0;
4256 static double K5 = 5.0e0;
4257 static double fx,cum,ccum,pq;
4258 static unsigned long qhi,qleft,qporq;
4259 static double T3,T6,T7,T8,T9,T10;
4262 .. Executable Statements ..
4267 if(!(*which < 1 || *which > 3)) goto S30;
4268 if(!(*which < 1)) goto S10;
4277 if(*which == 1) goto S70;
4281 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
4282 if(!(*p < 0.0e0)) goto S40;
4292 if(*which == 1) goto S110;
4296 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4297 if(!(*q <= 0.0e0)) goto S80;
4307 if(*which == 2) goto S130;
4311 if(!(*s < 0.0e0)) goto S120;
4317 if(*which == 3) goto S150;
4321 if(!(*xlam < 0.0e0)) goto S140;
4327 if(*which == 1) goto S190;
4332 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
4333 if(!(pq < 0.0e0)) goto S160;
4343 if(!(*which == 1)) qporq = *p <= *q;
4345 Select the minimum of P or Q
4355 else if(2 == *which) {
4363 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
4365 dinvr(status,s,&fx,&qleft,&qhi);
4367 if(!(*status == 1)) goto S230;
4368 cumpoi(s,xlam,&cum,&ccum);
4369 if(!qporq) goto S210;
4375 dinvr(status,s,&fx,&qleft,&qhi);
4378 if(!(*status == -1)) goto S260;
4379 if(!qleft) goto S240;
4390 else if(3 == *which) {
4398 dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
4400 dinvr(status,xlam,&fx,&qleft,&qhi);
4402 if(!(*status == 1)) goto S300;
4403 cumpoi(s,xlam,&cum,&ccum);
4404 if(!qporq) goto S280;
4410 dinvr(status,xlam,&fx,&qleft,&qhi);
4413 if(!(*status == -1)) goto S330;
4414 if(!qleft) goto S310;
4430 void cdft(int *which,double *p,double *q,double *t,double *df,
4431 int *status,double *bound)
4432 /**********************************************************************
4434 void cdft(int *which,double *p,double *q,double *t,double *df,
4435 int *status,double *bound)
4437 Cumulative Distribution Function
4444 Calculates any one parameter of the t distribution given
4445 values for the others.
4451 WHICH --> Integer indicating which argument
4452 values is to be calculated from the others.
4454 iwhich = 1 : Calculate P and Q from T and DF
4455 iwhich = 2 : Calculate T from P,Q and DF
4456 iwhich = 3 : Calculate DF from P,Q and T
4458 P <--> The integral from -infinity to t of the t-density.
4462 Input range: (0, 1].
4465 T <--> Upper limit of integration of the t-density.
4466 Input range: ( -infinity, +infinity).
4467 Search range: [ -1E100, 1E100 ]
4469 DF <--> Degrees of freedom of the t-distribution.
4470 Input range: (0 , +infinity).
4471 Search range: [1e-100, 1E10]
4473 STATUS <-- 0 if calculation completed correctly
4474 -I if input parameter number I is out of range
4475 1 if answer appears to be lower than lowest
4477 2 if answer appears to be higher than greatest
4481 BOUND <-- Undefined if STATUS is 0
4483 Bound exceeded by parameter number I if STATUS
4486 Lower search bound if STATUS is 1.
4488 Upper search bound if STATUS is 2.
4494 Formula 26.5.27 of Abramowitz and Stegun, Handbook of
4495 Mathematical Functions (1966) is used to reduce the computation
4496 of the cumulative distribution function to that of an incomplete
4499 Computation of other parameters involve a seach for a value that
4500 produces the desired value of P. The search relies on the
4501 monotinicity of P with the other parameter.
4503 **********************************************************************/
4506 #define atol 1.0e-50
4507 #define zero 1.0e-100
4509 #define rtinf 1.0e100
4510 #define maxdf 1.0e10
4512 static double K4 = 0.5e0;
4513 static double K5 = 5.0e0;
4514 static double fx,cum,ccum,pq;
4515 static unsigned long qhi,qleft,qporq;
4516 static double T2,T3,T6,T7,T8,T9,T10,T11;
4519 .. Executable Statements ..
4524 if(!(*which < 1 || *which > 3)) goto S30;
4525 if(!(*which < 1)) goto S10;
4534 if(*which == 1) goto S70;
4538 if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
4539 if(!(*p <= 0.0e0)) goto S40;
4549 if(*which == 1) goto S110;
4553 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4554 if(!(*q <= 0.0e0)) goto S80;
4564 if(*which == 3) goto S130;
4568 if(!(*df <= 0.0e0)) goto S120;
4574 if(*which == 1) goto S170;
4579 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S160;
4580 if(!(pq < 0.0e0)) goto S140;
4590 if(!(*which == 1)) qporq = *p <= *q;
4592 Select the minimum of P or Q
4602 else if(2 == *which) {
4605 .. Get initial approximation for T
4612 dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
4614 dinvr(status,t,&fx,&qleft,&qhi);
4616 if(!(*status == 1)) goto S210;
4617 cumt(t,df,&cum,&ccum);
4618 if(!qporq) goto S190;
4624 dinvr(status,t,&fx,&qleft,&qhi);
4627 if(!(*status == -1)) goto S240;
4628 if(!qleft) goto S220;
4639 else if(3 == *which) {
4648 dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
4650 dinvr(status,df,&fx,&qleft,&qhi);
4652 if(!(*status == 1)) goto S280;
4653 cumt(t,df,&cum,&ccum);
4654 if(!qporq) goto S260;
4660 dinvr(status,df,&fx,&qleft,&qhi);
4663 if(!(*status == -1)) goto S310;
4664 if(!qleft) goto S290;
4683 void cdftnc(int *which,double *p,double *q,double *t,double *df,
4684 double *pnonc,int *status,double *bound)
4685 /**********************************************************************
4687 void cdftnc(int *which,double *p,double *q,double *t,double *df,
4688 double *pnonc,int *status,double *bound)
4690 Cumulative Distribution Function
4691 Non-Central T distribution
4695 Calculates any one parameter of the noncentral t distribution give
4696 values for the others.
4700 WHICH --> Integer indicating which argument
4701 values is to be calculated from the others.
4703 iwhich = 1 : Calculate P and Q from T,DF,PNONC
4704 iwhich = 2 : Calculate T from P,Q,DF,PNONC
4705 iwhich = 3 : Calculate DF from P,Q,T
4706 iwhich = 4 : Calculate PNONC from P,Q,DF,T
4708 P <--> The integral from -infinity to t of the noncentral t-den
4712 Input range: (0, 1].
4715 T <--> Upper limit of integration of the noncentral t-density.
4716 Input range: ( -infinity, +infinity).
4717 Search range: [ -1E100, 1E100 ]
4719 DF <--> Degrees of freedom of the noncentral t-distribution.
4720 Input range: (0 , +infinity).
4721 Search range: [1e-100, 1E10]
4723 PNONC <--> Noncentrality parameter of the noncentral t-distributio
4724 Input range: [-infinity , +infinity).
4725 Search range: [-1e4, 1E4]
4727 STATUS <-- 0 if calculation completed correctly
4728 -I if input parameter number I is out of range
4729 1 if answer appears to be lower than lowest
4731 2 if answer appears to be higher than greatest
4735 BOUND <-- Undefined if STATUS is 0
4737 Bound exceeded by parameter number I if STATUS
4740 Lower search bound if STATUS is 1.
4742 Upper search bound if STATUS is 2.
4746 Upper tail of the cumulative noncentral t is calculated usin
4747 formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuou
4748 Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
4750 Computation of other parameters involve a seach for a value that
4751 produces the desired value of P. The search relies on the
4752 monotinicity of P with the other parameter.
4754 **********************************************************************/
4758 #define atol 1.0e-50
4759 #define zero 1.0e-100
4760 #define one ( 1.0e0 - 1.0e-16 )
4762 static double K3 = 0.5e0;
4763 static double K4 = 5.0e0;
4764 static double ccum,cum,fx;
4765 static unsigned long qhi,qleft;
4766 static double T1,T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14;
4769 .. Executable Statements ..
4771 if(!(*which < 1 || *which > 4)) goto S30;
4772 if(!(*which < 1)) goto S10;
4781 if(*which == 1) goto S70;
4782 if(!(*p < 0.0e0 || *p > one)) goto S60;
4783 if(!(*p < 0.0e0)) goto S40;
4793 if(*which == 3) goto S90;
4794 if(!(*df <= 0.0e0)) goto S80;
4800 if(*which == 4) goto S100;
4803 cumtnc(t,df,pnonc,p,q);
4806 else if(2 == *which) {
4812 dstinv(&T1,&T2,&K3,&K3,&K4,&T5,&T6);
4814 dinvr(status,t,&fx,&qleft,&qhi);
4816 if(!(*status == 1)) goto S120;
4817 cumtnc(t,df,pnonc,&cum,&ccum);
4819 dinvr(status,t,&fx,&qleft,&qhi);
4822 if(!(*status == -1)) goto S150;
4823 if(!qleft) goto S130;
4834 else if(3 == *which) {
4840 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
4842 dinvr(status,df,&fx,&qleft,&qhi);
4844 if(!(*status == 1)) goto S170;
4845 cumtnc(t,df,pnonc,&cum,&ccum);
4847 dinvr(status,df,&fx,&qleft,&qhi);
4850 if(!(*status == -1)) goto S200;
4851 if(!qleft) goto S180;
4862 else if(4 == *which) {
4868 dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
4870 dinvr(status,pnonc,&fx,&qleft,&qhi);
4872 if(!(*status == 1)) goto S220;
4873 cumtnc(t,df,pnonc,&cum,&ccum);
4875 dinvr(status,pnonc,&fx,&qleft,&qhi);
4878 if(!(*status == -1)) goto S250;
4879 if(!qleft) goto S230;
4898 void cumbet(double *x,double *y,double *a,double *b,double *cum,
4901 **********************************************************************
4903 void cumbet(double *x,double *y,double *a,double *b,double *cum,
4906 Double precision cUMulative incomplete BETa distribution
4912 Calculates the cdf to X of the incomplete beta distribution
4913 with parameters a and b. This is the integral from 0 to x
4914 of (1/B(a,b))*f(t)) where f(t) = t**(a-1) * (1-t)**(b-1)
4920 X --> Upper limit of integration.
4921 X is DOUBLE PRECISION
4924 Y is DOUBLE PRECISION
4926 A --> First parameter of the beta distribution.
4927 A is DOUBLE PRECISION
4929 B --> Second parameter of the beta distribution.
4930 B is DOUBLE PRECISION
4932 CUM <-- Cumulative incomplete beta distribution.
4933 CUM is DOUBLE PRECISION
4935 CCUM <-- Compliment of Cumulative incomplete beta distribution.
4936 CCUM is DOUBLE PRECISION
4942 Calls the routine BRATIO.
4946 Didonato, Armido R. and Morris, Alfred H. Jr. (1992) Algorithim
4947 708 Significant Digit Computation of the Incomplete Beta Function
4948 Ratios. ACM ToMS, Vol.18, No. 3, Sept. 1992, 360-373.
4950 **********************************************************************
4956 .. Executable Statements ..
4958 if(!(*x <= 0.0e0)) goto S10;
4963 if(!(*y <= 0.0e0)) goto S20;
4968 bratio(a,b,x,y,cum,ccum,&ierr);
4974 void cumbin(double *s,double *xn,double *pr,double *ompr,
4975 double *cum,double *ccum)
4977 **********************************************************************
4979 void cumbin(double *s,double *xn,double *pr,double *ompr,
4980 double *cum,double *ccum)
4982 CUmulative BINomial distribution
4988 Returns the probability of 0 to S successes in XN binomial
4989 trials, each of which has a probability of success, PBIN.
4995 S --> The upper limit of cumulation of the binomial distribution.
4996 S is DOUBLE PRECISION
4998 XN --> The number of binomial trials.
4999 XN is DOUBLE PRECISIO
5001 PBIN --> The probability of success in each binomial trial.
5002 PBIN is DOUBLE PRECIS
5005 OMPR is DOUBLE PRECIS
5007 CUM <-- Cumulative binomial distribution.
5008 CUM is DOUBLE PRECISI
5010 CCUM <-- Compliment of Cumulative binomial distribution.
5011 CCUM is DOUBLE PRECIS
5017 Formula 26.5.24 of Abramowitz and Stegun, Handbook of
5018 Mathematical Functions (1966) is used to reduce the binomial
5019 distribution to the cumulative beta distribution.
5021 **********************************************************************
5024 static double T1,T2;
5027 .. Executable Statements ..
5029 if(!(*s < *xn)) goto S10;
5032 cumbet(pr,ompr,&T1,&T2,ccum,cum);
5040 void cumchi(double *x,double *df,double *cum,double *ccum)
5042 **********************************************************************
5044 void cumchi(double *x,double *df,double *cum,double *ccum)
5045 CUMulative of the CHi-square distribution
5051 Calculates the cumulative chi-square distribution.
5057 X --> Upper limit of integration of the
5058 chi-square distribution.
5059 X is DOUBLE PRECISION
5061 DF --> Degrees of freedom of the
5062 chi-square distribution.
5063 DF is DOUBLE PRECISION
5065 CUM <-- Cumulative chi-square distribution.
5066 CUM is DOUBLE PRECISIO
5068 CCUM <-- Compliment of Cumulative chi-square distribution.
5069 CCUM is DOUBLE PRECISI
5075 Calls incomplete gamma function (CUMGAM)
5077 **********************************************************************
5083 .. Executable Statements ..
5087 cumgam(&xx,&a,cum,ccum);
5090 void cumchn(double *x,double *df,double *pnonc,double *cum,
5092 /**********************************************************************
5094 void cumchn(double *x,double *df,double *pnonc,double *cum,
5097 CUMulative of the Non-central CHi-square distribution
5101 Calculates the cumulative non-central chi-square
5102 distribution, i.e., the probability that a random variable
5103 which follows the non-central chi-square distribution, with
5104 non-centrality parameter PNONC and continuous degrees of
5105 freedom DF, is less than or equal to X.
5109 X --> Upper limit of integration of the non-central
5110 chi-square distribution.
5112 DF --> Degrees of freedom of the non-central
5113 chi-square distribution.
5115 PNONC --> Non-centrality parameter of the non-central
5116 chi-square distribution.
5118 CUM <-- Cumulative non-central chi-square distribution.
5120 CCUM <-- Compliment of Cumulative non-central chi-square distribut
5125 Uses formula 26.4.25 of Abramowitz and Stegun, Handbook of
5126 Mathematical Functions, US NBS (1966) to calculate the
5127 non-central chi-square.
5131 EPS --- Convergence criterion. The sum stops when a
5132 term is less than EPS*SUM.
5134 CCUM <-- Compliment of Cumulative non-central
5135 chi-square distribution.
5137 **********************************************************************/
5139 #define dg(i) (*df + 2.0e0 * (double)(i))
5140 #define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps * sum)
5141 static double eps = 1.0e-5;
5142 static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
5143 sumadj,term,wt,xnonc;
5145 static double T1,T2,T3;
5148 .. Executable Statements ..
5150 if(!(*x <= 0.0e0)) goto S10;
5155 if(!(*pnonc <= 1.0e-10 )) goto S20;
5157 When non-centrality parameter is (essentially) zero,
5158 use cumulative chi-square distribution
5160 cumchi(x,df,cum,ccum);
5163 xnonc = *pnonc / 2.0e0;
5165 ***********************************************************************
5166 The following code calcualtes the weight, chi-square, and
5167 adjustment term for the central term in the infinite series.
5168 The central term is the one in which the poisson weight is
5169 greatest. The adjustment term is the amount that must
5170 be subtracted from the chi-square to move up two degrees
5172 ***********************************************************************
5174 icent = fifidint(xnonc);
5175 if(icent == 0) icent = 1;
5178 Calculate central weight term
5180 T1 = (double)(icent + 1);
5181 lfact = alngam(&T1);
5182 lcntwt = -xnonc + (double)icent * log(xnonc) - lfact;
5183 centwt = exp(lcntwt);
5185 Calculate central chi-square
5188 cumchi(x,&T2,&pcent,ccum);
5190 Calculate central adjustment term
5192 dfd2 = dg(icent) / 2.0e0;
5194 lfact = alngam(&T3);
5195 lcntaj = dfd2 * log(chid2) - chid2 - lfact;
5196 centaj = exp(lcntaj);
5197 sum = centwt * pcent;
5199 ***********************************************************************
5200 Sum backwards from the central term towards zero.
5201 Quit whenever either
5202 (1) the zero term is reached, or
5203 (2) the term gets small relative to the sum
5204 ***********************************************************************
5212 if(qsmall(term) || i == 0) goto S50;
5214 dfd2 = dg(i) / 2.0e0;
5216 Adjust chi-square for two fewer degrees of freedom.
5217 The adjusted value ends up in PTERM.
5219 adj = adj * dfd2 / chid2;
5221 pterm = pcent + sumadj;
5223 Adjust poisson weight for J decreased by one
5225 wt *= ((double)i / xnonc);
5232 ***********************************************************************
5233 Now sum forward from the central term towards infinity.
5235 (1) the term gets small relative to the sum, or
5236 ***********************************************************************
5238 sumadj = adj = centaj;
5243 if(qsmall(term)) goto S80;
5246 Update weights for next higher J
5248 wt *= (xnonc / (double)(i + 1));
5250 Calculate PTERM and add term to sum
5252 pterm = pcent - sumadj;
5256 Update adjustment term for DF for next iteration
5259 dfd2 = dg(i) / 2.0e0;
5260 adj = adj * chid2 / dfd2;
5265 *ccum = 0.5e0 + (0.5e0 - *cum);
5270 void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
5272 **********************************************************************
5274 void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
5275 CUMulative F distribution
5281 Computes the integral from 0 to F of the f-density with DFN
5282 and DFD degrees of freedom.
5288 F --> Upper limit of integration of the f-density.
5289 F is DOUBLE PRECISION
5291 DFN --> Degrees of freedom of the numerator sum of squares.
5292 DFN is DOUBLE PRECISI
5294 DFD --> Degrees of freedom of the denominator sum of squares.
5295 DFD is DOUBLE PRECISI
5297 CUM <-- Cumulative f distribution.
5298 CUM is DOUBLE PRECISI
5300 CCUM <-- Compliment of Cumulative f distribution.
5301 CCUM is DOUBLE PRECIS
5307 Formula 26.5.28 of Abramowitz and Stegun is used to reduce
5308 the cumulative F to a cumulative beta distribution.
5314 If F is less than or equal to 0, 0 is returned.
5316 **********************************************************************
5321 static double dsum,prod,xx,yy;
5323 static double T1,T2;
5326 .. Executable Statements ..
5328 if(!(*f <= 0.0e0)) goto S10;
5335 XX is such that the incomplete beta with parameters
5336 DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM
5338 Calculate the smaller of XX and YY accurately
5349 bratio(&T1,&T2,&xx,&yy,ccum,cum,&ierr);
5354 void cumfnc(double *f,double *dfn,double *dfd,double *pnonc,
5355 double *cum,double *ccum)
5357 **********************************************************************
5359 F -NON- -C-ENTRAL F DISTRIBUTION
5366 COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
5367 DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC
5373 X --> UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION
5375 DFN --> DEGREES OF FREEDOM OF NUMERATOR
5377 DFD --> DEGREES OF FREEDOM OF DENOMINATOR
5379 PNONC --> NONCENTRALITY PARAMETER.
5381 CUM <-- CUMULATIVE NONCENTRAL F DISTRIBUTION
5383 CCUM <-- COMPLIMENT OF CUMMULATIVE
5389 USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
5390 SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
5391 (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
5392 THE CONVERGENCE CRITERION IS MET.
5394 FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
5401 HANDBOOD OF MATHEMATICAL FUNCTIONS
5402 EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
5403 NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
5405 P 947, EQUATIONS 26.6.17, 26.6.18
5411 THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
5412 TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20). EPS IS
5413 SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED.
5415 **********************************************************************
5418 #define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
5421 static double eps = 1.0e-4;
5422 static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
5424 static int i,icent,ierr;
5425 static double T1,T2,T3,T4,T5,T6;
5428 .. Executable Statements ..
5430 if(!(*f <= 0.0e0)) goto S10;
5435 if(!(*pnonc < 1.0e-10)) goto S20;
5437 Handle case in which the non-centrality parameter is
5440 cumf(f,dfn,dfd,cum,ccum);
5443 xnonc = *pnonc/2.0e0;
5445 Calculate the central term of the poisson weighting factor.
5447 icent = (long)(xnonc);
5448 if(icent == 0) icent = 1;
5450 Compute central weight term
5452 T1 = (double)(icent+1);
5453 centwt = exp(-xnonc+(double)icent*log(xnonc)-alngam(&T1));
5455 Compute central incomplete beta term
5456 Assure that minimum of arg to beta and 1 - arg is computed
5467 T2 = *dfn*half+(double)icent;
5469 bratio(&T2,&T3,&xx,&yy,&betdn,&dummy,&ierr);
5470 adn = *dfn/2.0e0+(double)icent;
5476 Now sum terms backward from icent until convergence or all done
5482 dnterm = exp(alngam(&T4)-alngam(&T5)-alngam(&b)+adn*log(xx)+b*log(yy));
5484 if(qsmall(xmult*betdn) || i <= 0) goto S40;
5485 xmult *= ((double)i/xnonc);
5488 dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
5490 sum += (xmult*betdn);
5495 Now sum forwards until convergence
5498 if(aup-1.0+b == 0) upterm = exp(-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+
5502 upterm = exp(alngam(&T6)-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+b*
5507 if(qsmall(xmult*betup)) goto S70;
5509 xmult *= (xnonc/(double)i);
5512 upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
5514 sum += (xmult*betup);
5518 *ccum = 0.5e0+(0.5e0-*cum);
5524 void cumgam(double *x,double *a,double *cum,double *ccum)
5526 **********************************************************************
5528 void cumgam(double *x,double *a,double *cum,double *ccum)
5529 Double precision cUMulative incomplete GAMma distribution
5535 Computes the cumulative of the incomplete gamma
5536 distribution, i.e., the integral from 0 to X of
5537 (1/GAM(A))*EXP(-T)*T**(A-1) DT
5538 where GAM(A) is the complete gamma function of A, i.e.,
5539 GAM(A) = integral from 0 to infinity of
5546 X --> The upper limit of integration of the incomplete gamma.
5547 X is DOUBLE PRECISION
5549 A --> The shape parameter of the incomplete gamma.
5550 A is DOUBLE PRECISION
5552 CUM <-- Cumulative incomplete gamma distribution.
5553 CUM is DOUBLE PRECISION
5555 CCUM <-- Compliment of Cumulative incomplete gamma distribution.
5556 CCUM is DOUBLE PRECISIO
5562 Calls the routine GRATIO.
5564 **********************************************************************
5570 .. Executable Statements ..
5572 if(!(*x <= 0.0e0)) goto S10;
5577 gratio(a,x,cum,ccum,&K1);
5583 void cumnbn(double *s,double *xn,double *pr,double *ompr,
5584 double *cum,double *ccum)
5586 **********************************************************************
5588 void cumnbn(double *s,double *xn,double *pr,double *ompr,
5589 double *cum,double *ccum)
5591 CUmulative Negative BINomial distribution
5597 Returns the probability that it there will be S or fewer failures
5598 before there are XN successes, with each binomial trial having
5599 a probability of success PR.
5601 Prob(# failures = S | XN successes, PR) =
5603 ( ) * PR^XN * (1-PR)^S
5610 S --> The number of failures
5611 S is DOUBLE PRECISION
5613 XN --> The number of successes
5614 XN is DOUBLE PRECISIO
5616 PR --> The probability of success in each binomial trial.
5617 PR is DOUBLE PRECISIO
5620 OMPR is DOUBLE PRECIS
5622 CUM <-- Cumulative negative binomial distribution.
5623 CUM is DOUBLE PRECISI
5625 CCUM <-- Compliment of Cumulative negative binomial distribution.
5626 CCUM is DOUBLE PRECIS
5632 Formula 26.5.26 of Abramowitz and Stegun, Handbook of
5633 Mathematical Functions (1966) is used to reduce the negative
5634 binomial distribution to the cumulative beta distribution.
5636 **********************************************************************
5642 .. Executable Statements ..
5645 cumbet(pr,ompr,xn,&T1,cum,ccum);
5648 void cumnor(double *arg,double *result,double *ccum)
5650 **********************************************************************
5652 void cumnor(double *arg,double *result,double *ccum)
5658 Computes the cumulative of the normal distribution, i.e.,
5659 the integral from -infinity to x of
5660 (1/sqrt(2*pi)) exp(-u*u/2) du
5662 X --> Upper limit of integration.
5663 X is DOUBLE PRECISION
5665 RESULT <-- Cumulative normal distribution.
5666 RESULT is DOUBLE PRECISION
5668 CCUM <-- Compliment of Cumulative normal distribution.
5669 CCUM is DOUBLE PRECISION
5671 Renaming of function ANORM from:
5673 Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
5674 Package of Special Function Routines and Test Drivers"
5675 acm Transactions on Mathematical Software. 19, 22-32.
5677 with slight modifications to return ccum and to deal with
5680 **********************************************************************
5682 ------------------------------------------------------------------
5684 This function evaluates the normal distribution function:
5688 P(x) = ----------- | e dt
5692 The main computation evaluates near-minimax approximations
5693 derived from those in "Rational Chebyshev approximations for
5694 the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
5695 This transportable program uses rational functions that
5696 theoretically approximate the normal distribution function to
5697 at least 18 significant decimal digits. The accuracy achieved
5698 depends on the arithmetic system, the compiler, the intrinsic
5699 functions, and proper selection of the machine-dependent
5702 *******************************************************************
5703 *******************************************************************
5705 Explanation of machine-dependent constants.
5707 MIN = smallest machine representable number.
5709 EPS = argument below which anorm(x) may be represented by
5710 0.5 and above which x*x will not underflow.
5711 A conservative value is the largest machine number X
5712 such that 1.0 + X = 1.0 to machine precision.
5713 *******************************************************************
5714 *******************************************************************
5718 The program returns ANORM = 0 for ARG .LE. XLOW.
5721 Intrinsic functions required are:
5727 Mathematics and Computer Science Division
5728 Argonne National Laboratory
5731 Latest modification: March 15, 1992
5733 ------------------------------------------------------------------
5736 static double a[5] = {
5737 2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
5738 1.8154981253343561249e04,6.5682337918207449113e-2
5740 static double b[4] = {
5741 4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
5742 4.5507789335026729956e04
5744 static double c[9] = {
5745 3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
5746 5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
5747 1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
5749 static double d[8] = {
5750 2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
5751 6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
5752 3.8912003286093271411e04,1.9685429676859990727e04
5754 static double half = 0.5e0;
5755 static double p[6] = {
5756 2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
5757 1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
5759 static double one = 1.0e0;
5760 static double q[5] = {
5761 1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
5762 3.78239633202758244e-3,7.29751555083966205e-5
5764 static double sixten = 1.60e0;
5765 static double sqrpi = 3.9894228040143267794e-1;
5766 static double thrsh = 0.66291e0;
5767 static double root32 = 5.656854248e0;
5768 static double zero = 0.0e0;
5772 static double del,eps,temp,x,xden,xnum,y,xsq,min;
5774 ------------------------------------------------------------------
5775 Machine dependent constants
5776 ------------------------------------------------------------------
5778 eps = spmpar(&K1)*0.5e0;
5784 ------------------------------------------------------------------
5785 Evaluate anorm for |X| <= 0.66291
5786 ------------------------------------------------------------------
5789 if(y > eps) xsq = x*x;
5792 for(i=0; i<3; i++) {
5793 xnum = (xnum+a[i])*xsq;
5794 xden = (xden+b[i])*xsq;
5796 *result = x*(xnum+a[3])/(xden+b[3]);
5798 *result = half+temp;
5802 ------------------------------------------------------------------
5803 Evaluate anorm for 0.66291 <= |X| <= sqrt(32)
5804 ------------------------------------------------------------------
5806 else if(y <= root32) {
5809 for(i=0; i<7; i++) {
5810 xnum = (xnum+c[i])*y;
5811 xden = (xden+d[i])*y;
5813 *result = (xnum+c[7])/(xden+d[7]);
5814 xsq = fifdint(y*sixten)/sixten;
5815 del = (y-xsq)*(y+xsq);
5816 *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
5817 *ccum = one-*result;
5825 ------------------------------------------------------------------
5826 Evaluate anorm for |X| > sqrt(32)
5827 ------------------------------------------------------------------
5834 for(i=0; i<4; i++) {
5835 xnum = (xnum+p[i])*xsq;
5836 xden = (xden+q[i])*xsq;
5838 *result = xsq*(xnum+p[4])/(xden+q[4]);
5839 *result = (sqrpi-*result)/y;
5840 xsq = fifdint(x*sixten)/sixten;
5841 del = (x-xsq)*(x+xsq);
5842 *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
5843 *ccum = one-*result;
5850 if(*result < min) *result = 0.0e0;
5852 ------------------------------------------------------------------
5853 Fix up for negative argument, erf, etc.
5854 ------------------------------------------------------------------
5855 ----------Last card of ANORM ----------
5857 if(*ccum < min) *ccum = 0.0e0;
5859 void cumpoi(double *s,double *xlam,double *cum,double *ccum)
5861 **********************************************************************
5863 void cumpoi(double *s,double *xlam,double *cum,double *ccum)
5864 CUMulative POIsson distribution
5870 Returns the probability of S or fewer events in a Poisson
5871 distribution with mean XLAM.
5877 S --> Upper limit of cumulation of the Poisson.
5878 S is DOUBLE PRECISION
5880 XLAM --> Mean of the Poisson distribution.
5881 XLAM is DOUBLE PRECIS
5883 CUM <-- Cumulative poisson distribution.
5884 CUM is DOUBLE PRECISION
5886 CCUM <-- Compliment of Cumulative poisson distribution.
5887 CCUM is DOUBLE PRECIS
5893 Uses formula 26.4.21 of Abramowitz and Stegun, Handbook of
5894 Mathematical Functions to reduce the cumulative Poisson to
5895 the cumulative chi-square distribution.
5897 **********************************************************************
5900 static double chi,df;
5903 .. Executable Statements ..
5905 df = 2.0e0*(*s+1.0e0);
5907 cumchi(&chi,&df,ccum,cum);
5910 void cumt(double *t,double *df,double *cum,double *ccum)
5912 **********************************************************************
5914 void cumt(double *t,double *df,double *cum,double *ccum)
5915 CUMulative T-distribution
5921 Computes the integral from -infinity to T of the t-density.
5927 T --> Upper limit of integration of the t-density.
5928 T is DOUBLE PRECISION
5930 DF --> Degrees of freedom of the t-distribution.
5931 DF is DOUBLE PRECISIO
5933 CUM <-- Cumulative t-distribution.
5934 CCUM is DOUBLE PRECIS
5936 CCUM <-- Compliment of Cumulative t-distribution.
5937 CCUM is DOUBLE PRECIS
5943 Formula 26.5.27 of Abramowitz and Stegun, Handbook of
5944 Mathematical Functions is used to reduce the t-distribution
5945 to an incomplete beta.
5947 **********************************************************************
5950 static double K2 = 0.5e0;
5951 static double xx,a,oma,tt,yy,dfptt,T1;
5954 .. Executable Statements ..
5961 cumbet(&xx,&yy,&T1,&K2,&a,&oma);
5962 if(!(*t <= 0.0e0)) goto S10;
5972 void cumtnc(double *t,double *df,double *pnonc,double *cum,
5974 /**********************************************************************
5976 void cumtnc(double *t,double *df,double *pnonc,double *cum,
5979 CUMulative Non-Central T-distribution
5985 Computes the integral from -infinity to T of the non-central
5992 T --> Upper limit of integration of the non-central t-density.
5994 DF --> Degrees of freedom of the non-central t-distribution.
5996 PNONC --> Non-centrality parameter of the non-central t distibutio
5998 CUM <-- Cumulative t-distribution.
6000 CCUM <-- Compliment of Cumulative t-distribution.
6005 Upper tail of the cumulative noncentral t using
6006 formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuous
6007 Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
6009 This implementation starts the calculation at i = lambda,
6010 which is near the largest Di. It then sums forward and backward.
6011 **********************************************************************/
6019 #define tiny 1.0e-10
6020 static double alghdf,b,bb,bbcent,bcent,cent,d,dcent,dpnonc,dum1,dum2,e,ecent,
6021 halfdf,lambda,lnomx,lnx,omx,pnonc2,s,scent,ss,sscent,t2,term,tt,twoi,x,xi,
6024 static unsigned long qrevs;
6025 static double T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
6028 .. Executable Statements ..
6031 Case pnonc essentially zero
6033 if(fabs(*pnonc) <= tiny) {
6034 cumt(t,df,cum,ccum);
6046 pnonc2 = dpnonc * dpnonc;
6048 if(fabs(tt) <= tiny) {
6050 cumnor(&T1,cum,ccum);
6053 lambda = half * pnonc2;
6054 x = *df / (*df + t2);
6058 halfdf = half * *df;
6059 alghdf = gamln(&halfdf);
6061 ******************** Case i = lambda
6063 cent = fifidint(lambda);
6064 if(cent < one) cent = one;
6066 Compute d=T(2i) in log space and offset by exp(-lambda)
6069 xlnd = cent * log(lambda) - gamln(&T2) - lambda;
6072 Compute e=t(2i+1) in log space offset by exp(-lambda)
6075 xlne = (cent + half) * log(lambda) - gamln(&T3) - lambda;
6077 if(dpnonc < zero) ecent = -ecent;
6079 Compute bcent=B(2*cent)
6082 bratio(&halfdf,&T4,&x,&omx,&bcent,&dum1,&ierr);
6084 compute bbcent=B(2*cent+1)
6087 bratio(&halfdf,&T5,&x,&omx,&bbcent,&dum2,&ierr);
6089 Case bcent and bbcent are essentially zero
6090 Thus t is effectively infinite
6092 if(bcent + bbcent < tiny) {
6104 Case bcent and bbcent are essentially one
6105 Thus t is effectively zero
6107 if(dum1 + dum2 < tiny) {
6109 cumnor(&T6,cum,ccum);
6113 First term in ccum is D*B + E*BB
6115 *ccum = dcent * bcent + ecent * bbcent;
6117 compute s(cent) = B(2*(cent+1)) - B(2*cent))
6119 T7 = halfdf + cent + half;
6121 scent = gamln(&T7) - gamln(&T8) - alghdf + halfdf * lnx + (cent + half) *
6125 compute ss(cent) = B(2*cent+3) - B(2*cent+1)
6127 T9 = halfdf + cent + one;
6129 sscent = gamln(&T9) - gamln(&T10) - alghdf + halfdf * lnx + (cent + one) *
6131 sscent = exp(sscent);
6133 ******************** Sum Forward
6146 d = lambda / xi * d;
6147 e = lambda / (xi + half) * e;
6148 term = d * b + e * bb;
6150 s = s * omx * (*df + twoi - one) / (twoi + one);
6151 ss = ss * omx * (*df + twoi) / (twoi + two);
6154 if(fabs(term) > conv * *ccum) goto S10;
6156 ******************** Sum Backward
6164 s = scent * (one + twoi) / ((*df + twoi - one) * omx);
6165 ss = sscent * (two + twoi) / ((*df + twoi) * omx);
6170 e *= ((xi + half) / lambda);
6171 term = d * b + e * bb;
6174 if(xi < half) goto S30;
6176 s = s * (one + twoi) / ((*df + twoi - one) * omx);
6177 ss = ss * (two + twoi) / ((*df + twoi) * omx);
6178 if(fabs(term) > conv * *ccum) goto S20;
6181 *cum = half * *ccum;
6185 *ccum = half * *ccum;
6189 Due to roundoff error the answer may not lie between zero and one
6192 *cum = fifdmax1(fifdmin1(*cum,one),zero);
6193 *ccum = fifdmax1(fifdmin1(*ccum,one),zero);
6203 double devlpl(double a[],int *n,double *x)
6205 **********************************************************************
6207 double devlpl(double a[],int *n,double *x)
6208 Double precision EVALuate a PoLynomial at X
6215 A(1) + A(2)*X + ... + A(N)*X**(N-1)
6221 A --> Array of coefficients of the polynomial.
6222 A is DOUBLE PRECISION(N)
6224 N --> Length of A, also degree of polynomial - 1.
6227 X --> Point at which the polynomial is to be evaluated.
6228 X is DOUBLE PRECISION
6230 **********************************************************************
6233 static double devlpl,term;
6237 .. Executable Statements ..
6240 for(i= *n-1-1; i>=0; i--) term = a[i]+term**x;
6244 double dinvnr(double *p,double *q)
6246 **********************************************************************
6248 double dinvnr(double *p,double *q)
6249 Double precision NoRmal distribution INVerse
6255 Returns X such that CUMNOR(X) = P, i.e., the integral from -
6256 infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
6262 P --> The probability whose normal deviate is sought.
6263 P is DOUBLE PRECISION
6266 P is DOUBLE PRECISION
6272 The rational function on page 95 of Kennedy and Gentle,
6273 Statistical Computing, Marcel Dekker, NY , 1980 is used as a start
6274 value for the Newton method of finding roots.
6280 If P or Q .lt. machine EPS returns +/- DINVNR(EPS)
6282 **********************************************************************
6287 #define r2pi 0.3989422804014326e0
6288 #define nhalf -0.5e0
6289 #define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
6290 static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
6292 static unsigned long qporq;
6295 .. Executable Statements ..
6298 FIND MINIMUM OF P AND Q
6301 if(!qporq) goto S10;
6310 strtx = stvaln(&pp);
6315 for(i=1; i<=maxit; i++) {
6316 cumnor(&xcur,&cum,&ccum);
6317 dx = (cum-pp)/dennor(xcur);
6319 if(fabs(dx/xcur) < eps) goto S40;
6323 IF WE GET HERE, NEWTON HAS FAILED
6325 if(!qporq) dinvnr = -dinvnr;
6329 IF WE GET HERE, NEWTON HAS SUCCEDED
6332 if(!qporq) dinvnr = -dinvnr;
6341 static void E0000(int IENTRY,int *status,double *x,double *fx,
6342 unsigned long *qleft,unsigned long *qhi,double *zabsst,
6343 double *zabsto,double *zbig,double *zrelst,
6344 double *zrelto,double *zsmall,double *zstpmu)
6346 #define qxmon(zx,zy,zz) (int)((zx) <= (zy) && (zy) <= (zz))
6347 static double absstp,abstol,big,fbig,fsmall,relstp,reltol,small,step,stpmul,xhi,
6348 xlb,xlo,xsave,xub,yy;
6350 static unsigned long qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup;
6351 switch(IENTRY){case 0: goto DINVR; case 1: goto DSTINV;}
6353 if(*status > 0) goto S310;
6354 qcond = !qxmon(small,*x,big);
6355 if(qcond) ftnstop((char *) " SMALL, X, BIG not monotone in INVR");
6358 See that SMALL and BIG bound the zero and set QINCR
6376 qincr = fbig > fsmall;
6377 if(!qincr) goto S50;
6378 if(fsmall <= 0.0e0) goto S30;
6383 if(fbig >= 0.0e0) goto S40;
6390 if(fsmall >= 0.0e0) goto S60;
6396 if(fbig <= 0.0e0) goto S70;
6404 step = fifdmax1(absstp,relstp*fabs(*x));
6413 if(!(yy == 0.0e0)) goto S100;
6418 qup = (qincr && yy < 0.0e0) || (!qincr && yy > 0.0e0);
6420 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6421 HANDLE CASE IN WHICH WE MUST STEP HIGHER
6422 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6426 xub = fifdmin1(xlb+step,big);
6429 if(qcond) goto S150;
6442 qbdd = (qincr && yy >= 0.0e0) || (!qincr && yy <= 0.0e0);
6444 qcond = qbdd || qlim;
6445 if(qcond) goto S140;
6448 xub = fifdmin1(xlb+step,big);
6452 if(!(qlim && !qbdd)) goto S160;
6462 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6463 HANDLE CASE IN WHICH WE MUST STEP LOWER
6464 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6467 xlb = fifdmax1(xub-step,small);
6470 if(qcond) goto S220;
6483 qbdd = (qincr && yy <= 0.0e0) || (!qincr && yy >= 0.0e0);
6484 qlim = xlb <= small;
6485 qcond = qbdd || qlim;
6486 if(qcond) goto S210;
6489 xlb = fifdmax1(xub-step,small);
6493 if(!(qlim && !qbdd)) goto S230;
6501 dstzr(&xlb,&xub,&abstol,&reltol);
6503 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6504 IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
6505 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6510 if(!(*status == 1)) goto S290;
6512 dzror(status,x,fx,&xlo,&xhi,&qdum1,&qdum2);
6513 if(!(*status == 1)) goto S280;
6537 TO GET-FUNCTION-VALUE
6542 switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S90;case
6543 4: goto S130;case 5: goto S200;case 6: goto S270;default: break;}
6546 void dinvr(int *status,double *x,double *fx,
6547 unsigned long *qleft,unsigned long *qhi)
6549 **********************************************************************
6551 void dinvr(int *status,double *x,double *fx,
6552 unsigned long *qleft,unsigned long *qhi)
6555 bounds the zero of the function and invokes zror
6556 Reverse Communication
6562 Bounds the function and invokes ZROR to perform the zero
6563 finding. STINVR must have been called before this routine
6564 in order to set its parameters.
6570 STATUS <--> At the beginning of a zero finding problem, STATUS
6571 should be set to 0 and INVR invoked. (The value
6572 of parameters other than X will be ignored on this cal
6574 When INVR needs the function evaluated, it will set
6575 STATUS to 1 and return. The value of the function
6576 should be set in FX and INVR again called without
6577 changing any of its other parameters.
6579 When INVR has finished without error, it will return
6580 with STATUS 0. In that case X is approximately a root
6583 If INVR cannot bound the function, it returns status
6584 -1 and sets QLEFT and QHI.
6587 X <-- The value of X at which F(X) is to be evaluated.
6590 FX --> The value of F(X) calculated when INVR returns with
6594 QLEFT <-- Defined only if QMFINV returns .FALSE. In that
6595 case it is .TRUE. If the stepping search terminated
6596 unsucessfully at SMALL. If it is .FALSE. the search
6597 terminated unsucessfully at BIG.
6600 QHI <-- Defined only if QMFINV returns .FALSE. In that
6601 case it is .TRUE. if F(X) .GT. Y at the termination
6602 of the search and .FALSE. if F(X) .LT. Y at the
6603 termination of the search.
6606 **********************************************************************
6609 E0000(0,status,x,fx,qleft,qhi,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
6611 void dstinv(double *zsmall,double *zbig,double *zabsst,
6612 double *zrelst,double *zstpmu,double *zabsto,
6615 **********************************************************************
6616 void dstinv(double *zsmall,double *zbig,double *zabsst,
6617 double *zrelst,double *zstpmu,double *zabsto,
6620 Double Precision - SeT INverse finder - Reverse Communication
6622 Concise Description - Given a monotone function F finds X
6623 such that F(X) = Y. Uses Reverse communication -- see invr.
6624 This routine sets quantities needed by INVR.
6625 More Precise Description of INVR -
6626 F must be a monotone function, the results of QMFINV are
6627 otherwise undefined. QINCR must be .TRUE. if F is non-
6628 decreasing and .FALSE. if F is non-increasing.
6629 QMFINV will return .TRUE. if and only if F(SMALL) and
6630 F(BIG) bracket Y, i. e.,
6631 QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
6632 QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
6633 if QMFINV returns .TRUE., then the X returned satisfies
6634 the following condition. let
6635 TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
6636 then if QINCR is .TRUE.,
6637 F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
6638 and if QINCR is .FALSE.
6639 F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
6641 SMALL --> The left endpoint of the interval to be
6642 searched for a solution.
6643 SMALL is DOUBLE PRECISION
6644 BIG --> The right endpoint of the interval to be
6645 searched for a solution.
6646 BIG is DOUBLE PRECISION
6647 ABSSTP, RELSTP --> The initial step size in the search
6648 is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
6649 ABSSTP is DOUBLE PRECISION
6650 RELSTP is DOUBLE PRECISION
6651 STPMUL --> When a step doesn't bound the zero, the step
6652 size is multiplied by STPMUL and another step
6653 taken. A popular value is 2.0
6654 DOUBLE PRECISION STPMUL
6655 ABSTOL, RELTOL --> Two numbers that determine the accuracy
6656 of the solution. See function for a precise definition.
6657 ABSTOL is DOUBLE PRECISION
6658 RELTOL is DOUBLE PRECISION
6660 Compares F(X) with Y for the input value of X then uses QINCR
6661 to determine whether to step left or right to bound the
6662 desired x. the initial step size is
6663 MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
6664 Iteratively steps right or left until it bounds X.
6665 At each step which doesn't bound X, the step size is doubled.
6666 The routine is careful never to step beyond SMALL or BIG. If
6667 it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
6668 after setting QLEFT and QHI.
6669 If X is successfully bounded then Algorithm R of the paper
6670 'Two Efficient Algorithms with Guaranteed Convergence for
6671 Finding a Zero of a Function' by J. C. P. Bus and
6672 T. J. Dekker in ACM Transactions on Mathematical
6673 Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
6674 to find the zero of the function F(X)-Y. This is routine
6676 **********************************************************************
6679 E0000(1,NULL,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
6682 double dt1(double *p,double *q,double *df)
6684 **********************************************************************
6686 double dt1(double *p,double *q,double *df)
6687 Double precision Initalize Approximation to
6688 INVerse of the cumulative T distribution
6694 Returns the inverse of the T distribution function, i.e.,
6695 the integral from 0 to INVT of the T density is P. This is an
6696 initial approximation
6702 P --> The p-value whose inverse from the T distribution is
6704 P is DOUBLE PRECISION
6707 Q is DOUBLE PRECISION
6709 DF --> Degrees of freedom of the T distribution.
6710 DF is DOUBLE PRECISION
6712 **********************************************************************
6715 static double coef[4][5] = {
6716 {1.0e0,1.0e0,0.0e0,0.0e0,0.0e0},
6717 {3.0e0,16.0e0,5.0e0,0.0e0,0.0e0},
6718 {-15.0e0,17.0e0,19.0e0,3.0e0,0.0e0},
6719 {-945.0e0,-1920.0e0,1482.0e0,776.0e0,79.0e0}
6721 static double denom[4] = {
6722 4.0e0,96.0e0,384.0e0,92160.0e0
6724 static int ideg[4] = {
6727 static double dt1,denpow,sum,term,x,xp,xx;
6731 .. Executable Statements ..
6733 x = fabs(dinvnr(p,q));
6737 for(i=0; i<4; i++) {
6738 term = devlpl(&coef[i][0],&ideg[i],&xx)*x;
6740 sum += (term/(denpow*denom[i]));
6742 if(!(*p >= 0.5e0)) goto S20;
6752 static void E0001(int IENTRY,int *status,double *x,double *fx,
6753 double *xlo,double *xhi,unsigned long *qleft,
6754 unsigned long *qhi,double *zabstl,double *zreltl,
6755 double *zxhi,double *zxlo)
6757 #define ftol(zx) (0.5e0*fifdmax1(abstol,reltol*fabs((zx))))
6758 static double a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q,reltol,tol,w,xxhi,xxlo;
6759 static int ext,i99999;
6760 static unsigned long first,qrzero;
6761 switch(IENTRY){case 0: goto DZROR; case 1: goto DSTZR;}
6763 if(*status > 0) goto S280;
6783 Check that F(ZXLO) < 0 < F(ZXHI) or
6784 F(ZXLO) > 0 > F(ZXHI)
6786 if(!(fb < 0.0e0)) goto S40;
6787 if(!(*fx < 0.0e0)) goto S30;
6794 if(!(fb > 0.0e0)) goto S60;
6795 if(!(*fx > 0.0e0)) goto S50;
6809 if(!(fabs(fc) < fabs(fb))) goto S100;
6810 if(!(c != a)) goto S90;
6825 if(!(fabs(mb) > tol)) goto S240;
6826 if(!(ext > 3)) goto S110;
6830 tol = fifdsign(tol,mb);
6832 if(!first) goto S120;
6837 fdb = (fd-fb)/(d-b);
6838 fda = (fd-fa)/(d-a);
6842 if(!(p < 0.0e0)) goto S140;
6846 if(ext == 3) p *= 2.0e0;
6847 if(!(p*1.0e0 == 0.0e0 || p <= q*tol)) goto S150;
6851 if(!(p < mb*q)) goto S160;
6873 if(!(fc*fb >= 0.0e0)) goto S210;
6876 if(!(w == mb)) goto S220;
6885 qrzero = (fc >= 0.0e0 && fb <= 0.0e0) || (fc < 0.0e0 && fb >= 0.0e0);
6886 if(!qrzero) goto S250;
6901 TO GET-FUNCTION-VALUE
6906 switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S200;
6910 void dzror(int *status,double *x,double *fx,double *xlo,
6911 double *xhi,unsigned long *qleft,unsigned long *qhi)
6913 **********************************************************************
6915 void dzror(int *status,double *x,double *fx,double *xlo,
6916 double *xhi,unsigned long *qleft,unsigned long *qhi)
6918 Double precision ZeRo of a function -- Reverse Communication
6924 Performs the zero finding. STZROR must have been called before
6925 this routine in order to set its parameters.
6931 STATUS <--> At the beginning of a zero finding problem, STATUS
6932 should be set to 0 and ZROR invoked. (The value
6933 of other parameters will be ignored on this call.)
6935 When ZROR needs the function evaluated, it will set
6936 STATUS to 1 and return. The value of the function
6937 should be set in FX and ZROR again called without
6938 changing any of its other parameters.
6940 When ZROR has finished without error, it will return
6941 with STATUS 0. In that case (XLO,XHI) bound the answe
6943 If ZROR finds an error (which implies that F(XLO)-Y an
6944 F(XHI)-Y have the same sign, it returns STATUS -1. In
6945 this case, XLO and XHI are undefined.
6948 X <-- The value of X at which F(X) is to be evaluated.
6951 FX --> The value of F(X) calculated when ZROR returns with
6955 XLO <-- When ZROR returns with STATUS = 0, XLO bounds the
6956 inverval in X containing the solution below.
6957 DOUBLE PRECISION XLO
6959 XHI <-- When ZROR returns with STATUS = 0, XHI bounds the
6960 inverval in X containing the solution above.
6961 DOUBLE PRECISION XHI
6963 QLEFT <-- .TRUE. if the stepping search terminated unsucessfully
6964 at XLO. If it is .FALSE. the search terminated
6965 unsucessfully at XHI.
6968 QHI <-- .TRUE. if F(X) .GT. Y at the termination of the
6969 search and .FALSE. if F(X) .LT. Y at the
6970 termination of the search.
6973 **********************************************************************
6976 E0001(0,status,x,fx,xlo,xhi,qleft,qhi,NULL,NULL,NULL,NULL);
6978 void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl)
6980 **********************************************************************
6981 void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl)
6982 Double precision SeT ZeRo finder - Reverse communication version
6984 Sets quantities needed by ZROR. The function of ZROR
6985 and the quantities set is given here.
6986 Concise Description - Given a function F
6987 find XLO such that F(XLO) = 0.
6988 More Precise Description -
6989 Input condition. F is a double precision function of a single
6990 double precision argument and XLO and XHI are such that
6991 F(XLO)*F(XHI) .LE. 0.0
6992 If the input condition is met, QRZERO returns .TRUE.
6993 and output values of XLO and XHI satisfy the following
6994 F(XLO)*F(XHI) .LE. 0.
6995 ABS(F(XLO) .LE. ABS(F(XHI)
6996 ABS(XLO-XHI) .LE. TOL(X)
6998 TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
6999 If this algorithm does not find XLO and XHI satisfying
7000 these conditions then QRZERO returns .FALSE. This
7001 implies that the input condition was not met.
7003 XLO --> The left endpoint of the interval to be
7004 searched for a solution.
7005 XLO is DOUBLE PRECISION
7006 XHI --> The right endpoint of the interval to be
7008 XHI is DOUBLE PRECISION
7009 ABSTOL, RELTOL --> Two numbers that determine the accuracy
7010 of the solution. See function for a
7012 ABSTOL is DOUBLE PRECISION
7013 RELTOL is DOUBLE PRECISION
7015 Algorithm R of the paper 'Two Efficient Algorithms with
7016 Guaranteed Convergence for Finding a Zero of a Function'
7017 by J. C. P. Bus and T. J. Dekker in ACM Transactions on
7018 Mathematical Software, Volume 1, no. 4 page 330
7019 (Dec. '75) is employed to find the zero of F(X)-Y.
7020 **********************************************************************
7023 E0001(1,NULL,NULL,NULL,NULL,NULL,NULL,NULL,zabstl,zreltl,zxhi,zxlo);
7025 double erf1(double *x)
7027 -----------------------------------------------------------------------
7028 EVALUATION OF THE REAL ERROR FUNCTION
7029 -----------------------------------------------------------------------
7032 static double c = .564189583547756e0;
7033 static double a[5] = {
7034 .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
7035 .479137145607681e-01,.128379167095513e+00
7037 static double b[3] = {
7038 .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
7040 static double p[8] = {
7041 -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
7042 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
7043 4.51918953711873e+02,3.00459261020162e+02
7045 static double q[8] = {
7046 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
7047 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7048 7.90950925327898e+02,3.00459260956983e+02
7050 static double r[5] = {
7051 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
7052 4.65807828718470e+00,2.82094791773523e-01
7054 static double s[4] = {
7055 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
7056 1.80124575948747e+01
7058 static double erf1,ax,bot,t,top,x2;
7061 .. Executable Statements ..
7064 if(ax > 0.5e0) goto S10;
7066 top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
7067 bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
7068 erf1 = *x*(top/bot);
7071 if(ax > 4.0e0) goto S20;
7072 top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7074 bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7076 erf1 = 0.5e0+(0.5e0-exp(-(*x**x))*top/bot);
7077 if(*x < 0.0e0) erf1 = -erf1;
7080 if(ax >= 5.8e0) goto S30;
7083 top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
7084 bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
7085 erf1 = (c-top/(x2*bot))/ax;
7086 erf1 = 0.5e0+(0.5e0-exp(-x2)*erf1);
7087 if(*x < 0.0e0) erf1 = -erf1;
7090 erf1 = fifdsign(1.0e0,*x);
7093 double erfc1(int *ind,double *x)
7095 -----------------------------------------------------------------------
7096 EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
7098 ERFC1(IND,X) = ERFC(X) IF IND = 0
7099 ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE
7100 -----------------------------------------------------------------------
7103 static double c = .564189583547756e0;
7104 static double a[5] = {
7105 .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
7106 .479137145607681e-01,.128379167095513e+00
7108 static double b[3] = {
7109 .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
7111 static double p[8] = {
7112 -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
7113 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
7114 4.51918953711873e+02,3.00459261020162e+02
7116 static double q[8] = {
7117 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
7118 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7119 7.90950925327898e+02,3.00459260956983e+02
7121 static double r[5] = {
7122 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
7123 4.65807828718470e+00,2.82094791773523e-01
7125 static double s[4] = {
7126 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
7127 1.80124575948747e+01
7130 static double erfc1,ax,bot,e,t,top,w;
7133 .. Executable Statements ..
7139 if(ax > 0.5e0) goto S10;
7141 top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
7142 bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
7143 erfc1 = 0.5e0+(0.5e0-*x*(top/bot));
7144 if(*ind != 0) erfc1 = exp(t)*erfc1;
7148 0.5 .LT. ABS(X) .LE. 4
7150 if(ax > 4.0e0) goto S20;
7151 top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7153 bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7161 if(*x <= -5.6e0) goto S60;
7162 if(*ind != 0) goto S30;
7163 if(*x > 100.0e0) goto S70;
7164 if(*x**x > -exparg(&K1)) goto S70;
7166 t = pow(1.0e0/ *x,2.0);
7167 top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
7168 bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
7169 erfc1 = (c-t*top/bot)/ax;
7174 if(*ind == 0) goto S50;
7175 if(*x < 0.0e0) erfc1 = 2.0e0*exp(*x**x)-erfc1;
7181 erfc1 = (0.5e0+(0.5e0-e))*exp(-t)*erfc1;
7182 if(*x < 0.0e0) erfc1 = 2.0e0-erfc1;
7186 LIMIT VALUE FOR LARGE NEGATIVE X
7189 if(*ind != 0) erfc1 = 2.0e0*exp(*x**x);
7193 LIMIT VALUE FOR LARGE POSITIVE X
7199 double esum(int *mu,double *x)
7201 -----------------------------------------------------------------------
7202 EVALUATION OF EXP(MU + X)
7203 -----------------------------------------------------------------------
7206 static double esum,w;
7209 .. Executable Statements ..
7211 if(*x > 0.0e0) goto S10;
7212 if(*mu < 0) goto S20;
7214 if(w > 0.0e0) goto S20;
7218 if(*mu > 0) goto S20;
7220 if(w < 0.0e0) goto S20;
7225 esum = exp(w)*exp(*x);
7228 double exparg(int *l)
7230 --------------------------------------------------------------------
7231 IF L = 0 THEN EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
7232 EXP(W) CAN BE COMPUTED.
7234 IF L IS NONZERO THEN EXPARG(L) = THE LARGEST NEGATIVE W FOR
7235 WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.
7237 NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.
7238 --------------------------------------------------------------------
7244 static double exparg,lnb;
7248 .. Executable Statements ..
7251 if(b != 2) goto S10;
7252 lnb = .69314718055995e0;
7255 if(b != 8) goto S20;
7256 lnb = 2.0794415416798e0;
7259 if(b != 16) goto S30;
7260 lnb = 2.7725887222398e0;
7263 lnb = log((double)b);
7265 if(*l == 0) goto S50;
7267 exparg = 0.99999e0*((double)m*lnb);
7271 exparg = 0.99999e0*((double)m*lnb);
7274 double fpser(double *a,double *b,double *x,double *eps)
7276 -----------------------------------------------------------------------
7278 EVALUATION OF I (A,B)
7281 FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
7283 -----------------------------------------------------------------------
7289 static double fpser,an,c,s,t,tol;
7292 .. Executable Statements ..
7295 if(*a <= 1.e-3**eps) goto S10;
7298 if(t < exparg(&K1)) return fpser;
7302 NOTE THAT 1/B(A,B) = B
7304 fpser = *b/ *a*fpser;
7314 if(fabs(c) > tol) goto S20;
7315 fpser *= (1.0e0+*a*s);
7318 double gam1(double *a)
7320 ------------------------------------------------------------------
7321 COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5
7322 ------------------------------------------------------------------
7325 static double s1 = .273076135303957e+00;
7326 static double s2 = .559398236957378e-01;
7327 static double p[7] = {
7328 .577215664901533e+00,-.409078193005776e+00,-.230975380857675e+00,
7329 .597275330452234e-01,.766968181649490e-02,-.514889771323592e-02,
7330 .589597428611429e-03
7332 static double q[5] = {
7333 .100000000000000e+01,.427569613095214e+00,.158451672430138e+00,
7334 .261132021441447e-01,.423244297896961e-02
7336 static double r[9] = {
7337 -.422784335098468e+00,-.771330383816272e+00,-.244757765222226e+00,
7338 .118378989872749e+00,.930357293360349e-03,-.118290993445146e-01,
7339 .223047661158249e-02,.266505979058923e-03,-.132674909766242e-03
7341 static double gam1,bot,d,t,top,w,T1;
7344 .. Executable Statements ..
7348 if(d > 0.0e0) t = d-0.5e0;
7350 if(T1 < 0) goto S40;
7351 else if(T1 == 0) goto S10;
7357 top = (((((p[6]*t+p[5])*t+p[4])*t+p[3])*t+p[2])*t+p[1])*t+p[0];
7358 bot = (((q[4]*t+q[3])*t+q[2])*t+q[1])*t+1.0e0;
7360 if(d > 0.0e0) goto S30;
7364 gam1 = t/ *a*(w-0.5e0-0.5e0);
7367 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+
7369 bot = (s2*t+s1)*t+1.0e0;
7371 if(d > 0.0e0) goto S50;
7372 gam1 = *a*(w+0.5e0+0.5e0);
7378 void gaminv(double *a,double *x,double *x0,double *p,double *q,
7381 ----------------------------------------------------------------------
7382 INVERSE INCOMPLETE GAMMA RATIO FUNCTION
7384 GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1.
7385 THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER
7386 ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X
7387 TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE
7388 PARTICULAR COMPUTER ARITHMETIC BEING USED.
7392 X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0,
7393 AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT
7394 NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN
7395 A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE
7396 IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X.
7398 X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER
7399 DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET
7402 IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
7403 WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING
7406 IERR = 0 THE SOLUTION WAS OBTAINED. ITERATION WAS
7408 IERR.GT.0 THE SOLUTION WAS OBTAINED. IERR ITERATIONS
7410 IERR = -2 (INPUT ERROR) A .LE. 0
7411 IERR = -3 NO SOLUTION WAS OBTAINED. THE RATIO Q/A
7413 IERR = -4 (INPUT ERROR) P + Q .NE. 1
7414 IERR = -6 20 ITERATIONS WERE PERFORMED. THE MOST
7415 RECENT VALUE OBTAINED FOR X IS GIVEN.
7416 THIS CANNOT OCCUR IF X0 .LE. 0.
7417 IERR = -7 ITERATION FAILED. NO VALUE IS GIVEN FOR X.
7418 THIS MAY OCCUR WHEN X IS APPROXIMATELY 0.
7419 IERR = -8 A VALUE FOR X HAS BEEN OBTAINED, BUT THE
7420 ROUTINE IS NOT CERTAIN OF ITS ACCURACY.
7421 ITERATION CANNOT BE PERFORMED IN THIS
7422 CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY
7423 WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS
7424 POSITIVE THEN THIS CAN OCCUR WHEN A IS
7425 EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY
7426 LARGE (SAY A .GE. 1.E20).
7427 ----------------------------------------------------------------------
7428 WRITTEN BY ALFRED H. MORRIS, JR.
7429 NAVAL SURFACE WEAPONS CENTER
7434 static double a0 = 3.31125922108741e0;
7435 static double a1 = 11.6616720288968e0;
7436 static double a2 = 4.28342155967104e0;
7437 static double a3 = .213623493715853e0;
7438 static double b1 = 6.61053765625462e0;
7439 static double b2 = 6.40691597760039e0;
7440 static double b3 = 1.27364489782223e0;
7441 static double b4 = .036117081018842e0;
7442 static double c = .577215664901533e0;
7443 static double ln10 = 2.302585e0;
7444 static double tol = 1.e-5;
7445 static double amin[2] = {
7448 static double bmin[2] = {
7451 static double dmin[2] = {
7454 static double emin[2] = {
7457 static double eps0[2] = {
7464 static double am1,amax,ap1,ap2,ap3,apn,b,c1,c2,c3,c4,c5,d,e,e2,eps,g,h,pn,qg,qn,
7465 r,rta,s,s2,sum,t,u,w,xmax,xmin,xn,y,z;
7467 static double T4,T5,T6,T7,T9;
7470 .. Executable Statements ..
7473 ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
7474 E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
7475 XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
7476 LARGEST POSITIVE NUMBER.
7482 if(*a <= 0.0e0) goto S300;
7484 if(fabs(t) > e) goto S320;
7486 if(*p == 0.0e0) return;
7487 if(*q == 0.0e0) goto S270;
7488 if(*a == 1.0e0) goto S280;
7490 amax = 0.4e-10/(e*e);
7492 if(e > 1.e-10) iop = 2;
7495 if(*x0 > 0.0e0) goto S160;
7497 SELECTION OF THE INITIAL APPROXIMATION XN OF X
7500 if(*a > 1.0e0) goto S80;
7504 if(qg == 0.0e0) goto S360;
7506 if(qg > 0.6e0**a) goto S40;
7507 if(*a >= 0.30e0 || b < 0.35e0) goto S10;
7513 if(b >= 0.45e0) goto S40;
7514 if(b == 0.0e0) goto S360;
7516 s = 0.5e0+(0.5e0-*a);
7519 if(b < 0.15e0) goto S20;
7520 xn = y-s*log(t)-log(1.0e0+s/(t+1.0e0));
7523 if(b <= 0.01e0) goto S30;
7524 u = ((t+2.0e0*(3.0e0-*a))*t+(2.0e0-*a)*(3.0e0-*a))/((t+(5.0e0-*a))*t+2.0e0);
7525 xn = y-s*log(t)-log(u);
7529 c2 = -(s*(1.0e0+c1));
7530 c3 = s*((0.5e0*c1+(2.0e0-*a))*c1+(2.5e0-1.5e0**a));
7531 c4 = -(s*(((c1/3.0e0+(2.5e0-1.5e0**a))*c1+((*a-6.0e0)**a+7.0e0))*c1+(
7532 (11.0e0**a-46.0)**a+47.0e0)/6.0e0));
7533 c5 = -(s*((((-(c1/4.0e0)+(11.0e0**a-17.0e0)/6.0e0)*c1+((-(3.0e0**a)+13.0e0)*
7534 *a-13.0e0))*c1+0.5e0*(((2.0e0**a-25.0e0)**a+72.0e0)**a-61.0e0))*c1+((
7535 (25.0e0**a-195.0e0)**a+477.0e0)**a-379.0e0)/12.0e0));
7536 xn = (((c5/y+c4)/y+c3)/y+c2)/y+c1+y;
7537 if(*a > 1.0e0) goto S220;
7538 if(b > bmin[iop-1]) goto S220;
7542 if(b**q > 1.e-8) goto S50;
7543 xn = exp(-(*q/ *a+c));
7546 if(*p <= 0.9e0) goto S60;
7548 xn = exp((alnrel(&T5)+gamln1(a))/ *a);
7551 xn = exp(log(*p*g)/ *a);
7553 if(xn == 0.0e0) goto S310;
7554 t = 0.5e0+(0.5e0-xn/(*a+1.0e0));
7559 SELECTION OF THE INITIAL APPROXIMATION XN OF X
7562 if(*q <= 0.5e0) goto S90;
7568 t = sqrt(-(2.0e0*w));
7569 s = t-(((a3*t+a2)*t+a1)*t+a0)/((((b4*t+b3)*t+b2)*t+b1)*t+1.0e0);
7570 if(*q > 0.5e0) s = -s;
7573 xn = *a+s*rta+(s2-1.0e0)/3.0e0+s*(s2-7.0e0)/(36.0e0*rta)-((3.0e0*s2+7.0e0)*
7574 s2-16.0e0)/(810.0e0**a)+s*((9.0e0*s2+256.0e0)*s2-433.0e0)/(38880.0e0**a*
7576 xn = fifdmax1(xn,0.0e0);
7577 if(*a < amin[iop-1]) goto S110;
7579 d = 0.5e0+(0.5e0-*x/ *a);
7580 if(fabs(d) <= dmin[iop-1]) return;
7582 if(*p <= 0.5e0) goto S130;
7583 if(xn < 3.0e0**a) goto S220;
7585 d = fifdmax1(2.0e0,*a*(*a-1.0e0));
7586 if(y < ln10*d) goto S120;
7592 T6 = -(t/(xn+1.0e0));
7593 xn = y+t*log(xn)-alnrel(&T6);
7594 T7 = -(t/(xn+1.0e0));
7595 xn = y+t*log(xn)-alnrel(&T7);
7599 if(xn > 0.70e0*ap1) goto S170;
7601 if(xn > 0.15e0*ap1) goto S140;
7604 *x = exp((w+*x)/ *a);
7605 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
7606 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
7607 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2*(1.0e0+*x/ap3))))/ *a);
7609 if(xn > 1.e-2*ap1) goto S140;
7610 if(xn <= emin[iop-1]*ap1) return;
7620 if(t > 1.e-4) goto S150;
7622 xn = exp((xn+t)/ *a);
7623 xn *= (1.0e0-(*a*log(xn)-xn-t)/(*a-xn));
7627 SCHRODER ITERATION USING P
7629 if(*p > 0.5e0) goto S220;
7631 if(*p <= 1.e10*xmin) goto S350;
7632 am1 = *a-0.5e0-0.5e0;
7634 if(*a <= amax) goto S190;
7635 d = 0.5e0+(0.5e0-xn/ *a);
7636 if(fabs(d) <= e2) goto S350;
7638 if(*ierr >= 20) goto S330;
7640 gratio(a,&xn,&pn,&qn,&K8);
7641 if(pn == 0.0e0 || qn == 0.0e0) goto S350;
7643 if(r == 0.0e0) goto S350;
7646 if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S200;
7648 if(*x <= 0.0e0) goto S340;
7654 if(*x <= 0.0e0) goto S340;
7655 if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
7659 if(d > tol) goto S180;
7660 if(d <= eps) return;
7661 if(fabs(*p-pn) <= tol**p) return;
7665 SCHRODER ITERATION USING Q
7667 if(*q <= 1.e10*xmin) goto S350;
7668 am1 = *a-0.5e0-0.5e0;
7670 if(*a <= amax) goto S240;
7671 d = 0.5e0+(0.5e0-xn/ *a);
7672 if(fabs(d) <= e2) goto S350;
7674 if(*ierr >= 20) goto S330;
7676 gratio(a,&xn,&pn,&qn,&K8);
7677 if(pn == 0.0e0 || qn == 0.0e0) goto S350;
7679 if(r == 0.0e0) goto S350;
7682 if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S250;
7684 if(*x <= 0.0e0) goto S340;
7690 if(*x <= 0.0e0) goto S340;
7691 if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
7695 if(d > tol) goto S230;
7696 if(d <= eps) return;
7697 if(fabs(*q-qn) <= tol**q) return;
7706 if(*q < 0.9e0) goto S290;
7740 double gamln(double *a)
7742 -----------------------------------------------------------------------
7743 EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
7744 -----------------------------------------------------------------------
7745 WRITTEN BY ALFRED H. MORRIS
7746 NAVAL SURFACE WARFARE CENTER
7748 --------------------------
7749 D = 0.5*(LN(2*PI) - 1)
7750 --------------------------
7753 static double c0 = .833333333333333e-01;
7754 static double c1 = -.277777777760991e-02;
7755 static double c2 = .793650666825390e-03;
7756 static double c3 = -.595202931351870e-03;
7757 static double c4 = .837308034031215e-03;
7758 static double c5 = -.165322962780713e-02;
7759 static double d = .418938533204673e0;
7760 static double gamln,t,w;
7765 .. Executable Statements ..
7767 if(*a > 0.8e0) goto S10;
7768 gamln = gamln1(a)-log(*a);
7771 if(*a > 2.25e0) goto S20;
7776 if(*a >= 10.0e0) goto S40;
7777 n = (long)(*a - 1.25e0);
7780 for(i=1; i<=n; i++) {
7785 gamln = gamln1(&T1)+log(w);
7788 t = pow(1.0e0/ *a,2.0);
7789 w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a;
7790 gamln = d+w+(*a-0.5e0)*(log(*a)-1.0e0);
7793 double gamln1(double *a)
7795 -----------------------------------------------------------------------
7796 EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
7797 -----------------------------------------------------------------------
7800 static double p0 = .577215664901533e+00;
7801 static double p1 = .844203922187225e+00;
7802 static double p2 = -.168860593646662e+00;
7803 static double p3 = -.780427615533591e+00;
7804 static double p4 = -.402055799310489e+00;
7805 static double p5 = -.673562214325671e-01;
7806 static double p6 = -.271935708322958e-02;
7807 static double q1 = .288743195473681e+01;
7808 static double q2 = .312755088914843e+01;
7809 static double q3 = .156875193295039e+01;
7810 static double q4 = .361951990101499e+00;
7811 static double q5 = .325038868253937e-01;
7812 static double q6 = .667465618796164e-03;
7813 static double r0 = .422784335098467e+00;
7814 static double r1 = .848044614534529e+00;
7815 static double r2 = .565221050691933e+00;
7816 static double r3 = .156513060486551e+00;
7817 static double r4 = .170502484022650e-01;
7818 static double r5 = .497958207639485e-03;
7819 static double s1 = .124313399877507e+01;
7820 static double s2 = .548042109832463e+00;
7821 static double s3 = .101552187439830e+00;
7822 static double s4 = .713309612391000e-02;
7823 static double s5 = .116165475989616e-03;
7824 static double gamln1,w,x;
7827 .. Executable Statements ..
7829 if(*a >= 0.6e0) goto S10;
7830 w = ((((((p6**a+p5)**a+p4)**a+p3)**a+p2)**a+p1)**a+p0)/((((((q6**a+q5)**a+
7831 q4)**a+q3)**a+q2)**a+q1)**a+1.0e0);
7836 w = (((((r5*x+r4)*x+r3)*x+r2)*x+r1)*x+r0)/(((((s5*x+s4)*x+s3)*x+s2)*x+s1)*x
7841 double Xgamm(double *a)
7843 -----------------------------------------------------------------------
7845 EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS
7849 GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
7852 -----------------------------------------------------------------------
7853 WRITTEN BY ALFRED H. MORRIS, JR.
7854 NAVAL SURFACE WEAPONS CENTER
7856 -----------------------------------------------------------------------
7859 static double d = .41893853320467274178e0;
7860 static double pi = 3.1415926535898e0;
7861 static double r1 = .820756370353826e-03;
7862 static double r2 = -.595156336428591e-03;
7863 static double r3 = .793650663183693e-03;
7864 static double r4 = -.277777777770481e-02;
7865 static double r5 = .833333333333333e-01;
7866 static double p[7] = {
7867 .539637273585445e-03,.261939260042690e-02,.204493667594920e-01,
7868 .730981088720487e-01,.279648642639792e+00,.553413866010467e+00,1.0e0
7870 static double q[7] = {
7871 -.832979206704073e-03,.470059485860584e-02,.225211131035340e-01,
7872 -.170458969313360e+00,-.567902761974940e-01,.113062953091122e+01,1.0e0
7876 static double Xgamm,bot,g,lnx,s,t,top,w,x,z;
7877 static int i,j,m,n,T1;
7880 .. Executable Statements ..
7884 if(fabs(*a) >= 15.0e0) goto S110;
7886 -----------------------------------------------------------------------
7887 EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15
7888 -----------------------------------------------------------------------
7893 LET T BE THE PRODUCT OF A-J WHEN A .GE. 2
7896 if(T1 < 0) goto S40;
7897 else if(T1 == 0) goto S30;
7900 for(j=1; j<=m; j++) {
7909 LET T BE THE PRODUCT OF A+J WHEN A .LT. 1
7912 if(*a > 0.0e0) goto S70;
7914 if(m == 0) goto S60;
7915 for(j=1; j<=m; j++) {
7922 if(t == 0.0e0) return Xgamm;
7925 THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
7926 CODE MAY BE OMITTED IF DESIRED.
7928 if(fabs(t) >= 1.e-30) goto S80;
7929 if(fabs(t)*spmpar(&K2) <= 1.0001e0) return Xgamm;
7934 COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1
7938 for(i=1; i<7; i++) {
7946 if(*a < 1.0e0) goto S100;
7954 -----------------------------------------------------------------------
7955 EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15
7956 -----------------------------------------------------------------------
7958 if(fabs(*a) >= 1.e3) return Xgamm;
7959 if(*a > 0.0e0) goto S120;
7963 if(t > 0.9e0) t = 1.0e0-t;
7965 if(fifmod(n,2) == 0) s = -s;
7966 if(s == 0.0e0) return Xgamm;
7969 COMPUTE THE MODIFIED ASYMPTOTIC SUM
7972 g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x;
7974 ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X)
7975 BUT LESS ACCURACY WILL NORMALLY BE OBTAINED.
7982 g = d+g+(z-0.5e0)*(lnx-1.e0);
7985 if(w > 0.99999e0*exparg(&K3)) return Xgamm;
7986 Xgamm = exp(w)*(1.0e0+t);
7987 if(*a < 0.0e0) Xgamm = 1.0e0/(Xgamm*s)/x;
7990 void grat1(double *a,double *x,double *r,double *p,double *q,
7994 static double a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,t,tol,w,z,T1,T3;
7997 .. Executable Statements ..
8000 -----------------------------------------------------------------------
8001 EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
8003 IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED.
8004 THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A).
8005 -----------------------------------------------------------------------
8007 if(*a**x == 0.0e0) goto S120;
8008 if(*a == 0.5e0) goto S100;
8009 if(*x < 1.1e0) goto S10;
8013 TAYLOR SERIES FOR P(A,X)/X**A
8017 sum = *x/(*a+3.0e0);
8018 tol = 0.1e0**eps/(*a+1.0e0);
8024 if(fabs(t) > tol) goto S20;
8025 j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
8029 if(*x < 0.25e0) goto S30;
8030 if(*a < *x/2.59e0) goto S50;
8033 if(z > -.13394e0) goto S50;
8036 *p = w*g*(0.5e0+(0.5e0-j));
8037 *q = 0.5e0+(0.5e0-*p);
8041 w = 0.5e0+(0.5e0+l);
8043 if(*q < 0.0e0) goto S90;
8044 *p = 0.5e0+(0.5e0-*q);
8048 CONTINUED FRACTION EXPANSION
8050 a2nm1 = a2n = 1.0e0;
8052 b2n = *x+(1.0e0-*a);
8055 a2nm1 = *x*a2n+c*a2nm1;
8056 b2nm1 = *x*b2n+c*b2nm1;
8060 a2n = a2nm1+cma*a2n;
8061 b2n = b2nm1+cma*b2n;
8063 if(fabs(an0-am0) >= *eps*an0) goto S70;
8065 *p = 0.5e0+(0.5e0-*q);
8079 if(*x >= 0.25e0) goto S110;
8082 *q = 0.5e0+(0.5e0-*p);
8086 *q = erfc1(&K2,&T3);
8087 *p = 0.5e0+(0.5e0-*q);
8090 if(*x <= *a) goto S80;
8093 void gratio(double *a,double *x,double *ans,double *qans,int *ind)
8095 ----------------------------------------------------------------------
8096 EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
8101 IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
8104 ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
8105 P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
8106 IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
8107 POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
8108 IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
8109 6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
8110 IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.
8113 ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
8114 WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
8115 P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
8116 X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.
8117 ----------------------------------------------------------------------
8118 WRITTEN BY ALFRED H. MORRIS, JR.
8119 NAVAL SURFACE WEAPONS CENTER
8121 --------------------
8124 static double alog10 = 2.30258509299405e0;
8125 static double d10 = -.185185185185185e-02;
8126 static double d20 = .413359788359788e-02;
8127 static double d30 = .649434156378601e-03;
8128 static double d40 = -.861888290916712e-03;
8129 static double d50 = -.336798553366358e-03;
8130 static double d60 = .531307936463992e-03;
8131 static double d70 = .344367606892378e-03;
8132 static double rt2pin = .398942280401433e0;
8133 static double rtpi = 1.77245385090552e0;
8134 static double third = .333333333333333e0;
8135 static double acc0[3] = {
8138 static double big[3] = {
8139 20.0e0,14.0e0,10.0e0
8141 static double d0[13] = {
8142 .833333333333333e-01,-.148148148148148e-01,.115740740740741e-02,
8143 .352733686067019e-03,-.178755144032922e-03,.391926317852244e-04,
8144 -.218544851067999e-05,-.185406221071516e-05,.829671134095309e-06,
8145 -.176659527368261e-06,.670785354340150e-08,.102618097842403e-07,
8146 -.438203601845335e-08
8148 static double d1[12] = {
8149 -.347222222222222e-02,.264550264550265e-02,-.990226337448560e-03,
8150 .205761316872428e-03,-.401877572016461e-06,-.180985503344900e-04,
8151 .764916091608111e-05,-.161209008945634e-05,.464712780280743e-08,
8152 .137863344691572e-06,-.575254560351770e-07,.119516285997781e-07
8154 static double d2[10] = {
8155 -.268132716049383e-02,.771604938271605e-03,.200938786008230e-05,
8156 -.107366532263652e-03,.529234488291201e-04,-.127606351886187e-04,
8157 .342357873409614e-07,.137219573090629e-05,-.629899213838006e-06,
8158 .142806142060642e-06
8160 static double d3[8] = {
8161 .229472093621399e-03,-.469189494395256e-03,.267720632062839e-03,
8162 -.756180167188398e-04,-.239650511386730e-06,.110826541153473e-04,
8163 -.567495282699160e-05,.142309007324359e-05
8165 static double d4[6] = {
8166 .784039221720067e-03,-.299072480303190e-03,-.146384525788434e-05,
8167 .664149821546512e-04,-.396836504717943e-04,.113757269706784e-04
8169 static double d5[4] = {
8170 -.697281375836586e-04,.277275324495939e-03,-.199325705161888e-03,
8171 .679778047793721e-04
8173 static double d6[2] = {
8174 -.592166437353694e-03,.270878209671804e-03
8176 static double e00[3] = {
8179 static double x00[3] = {
8184 static double a2n,a2nm1,acc,am0,amn,an,an0,apn,b2n,b2nm1,c,c0,c1,c2,c3,c4,c5,c6,
8185 cma,e,e0,g,h,j,l,r,rta,rtx,s,sum,t,t1,tol,twoa,u,w,x0,y,z;
8186 static int i,iop,m,max,n;
8187 static double wk[20],T3;
8189 static double T6,T7;
8192 .. Executable Statements ..
8195 --------------------
8196 ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
8197 FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
8200 if(*a < 0.0e0 || *x < 0.0e0) goto S430;
8201 if(*a == 0.0e0 && *x == 0.0e0) goto S430;
8202 if(*a**x == 0.0e0) goto S420;
8204 if(iop != 1 && iop != 2) iop = 3;
8205 acc = fifdmax1(acc0[iop-1],e);
8209 SELECT THE APPROPRIATE ALGORITHM
8211 if(*a >= 1.0e0) goto S10;
8212 if(*a == 0.5e0) goto S390;
8213 if(*x < 1.1e0) goto S160;
8216 if(u == 0.0e0) goto S380;
8217 r = u*(1.0e0+gam1(a));
8220 if(*a >= big[iop-1]) goto S30;
8221 if(*a > *x || *x >= x0) goto S20;
8224 if(twoa != (double)m) goto S20;
8226 if(*a == (double)i) goto S210;
8230 r = exp(t1)/Xgamm(a);
8234 if(l == 0.0e0) goto S370;
8235 s = 0.5e0+(0.5e0-l);
8237 if(z >= 700.0e0/ *a) goto S410;
8240 if(fabs(s) <= e0/rta) goto S330;
8241 if(fabs(s) <= 0.4e0) goto S270;
8242 t = pow(1.0e0/ *a,2.0);
8243 t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
8245 r = rt2pin*rta*exp(t1);
8247 if(r == 0.0e0) goto S420;
8248 if(*x <= fifdmax1(*a,alog10)) goto S50;
8249 if(*x < x0) goto S250;
8253 TAYLOR SERIES FOR P/R
8258 for(n=2; n<=20; n++) {
8261 if(t <= 1.e-3) goto S70;
8272 if(t > tol) goto S80;
8274 for(m=1; m<=max; m++) {
8278 *ans = r/ *a*(1.0e0+sum);
8279 *qans = 0.5e0+(0.5e0-*ans);
8283 ASYMPTOTIC EXPANSION
8288 for(n=2; n<=20; n++) {
8291 if(fabs(t) <= 1.e-3) goto S120;
8298 if(fabs(t) <= acc) goto S140;
8305 for(m=1; m<=max; m++) {
8309 *qans = r/ *x*(1.0e0+sum);
8310 *ans = 0.5e0+(0.5e0-*qans);
8314 TAYLOR SERIES FOR P(A,X)/X**A
8318 sum = *x/(*a+3.0e0);
8319 tol = 3.0e0*acc/(*a+1.0e0);
8325 if(fabs(t) > tol) goto S170;
8326 j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
8330 if(*x < 0.25e0) goto S180;
8331 if(*a < *x/2.59e0) goto S200;
8334 if(z > -.13394e0) goto S200;
8337 *ans = w*g*(0.5e0+(0.5e0-j));
8338 *qans = 0.5e0+(0.5e0-*ans);
8342 w = 0.5e0+(0.5e0+l);
8343 *qans = (w*j-l)*g-h;
8344 if(*qans < 0.0e0) goto S380;
8345 *ans = 0.5e0+(0.5e0-*qans);
8349 FINITE SUMS FOR Q WHEN A .GE. 1
8350 AND 2*A IS AN INTEGER
8359 sum = erfc1(&K2,&rtx);
8360 t = exp(-*x)/(rtpi*rtx);
8364 if(n == i) goto S240;
8372 *ans = 0.5e0+(0.5e0-*qans);
8376 CONTINUED FRACTION EXPANSION
8378 tol = fifdmax1(5.0e0*e,acc);
8379 a2nm1 = a2n = 1.0e0;
8381 b2n = *x+(1.0e0-*a);
8384 a2nm1 = *x*a2n+c*a2nm1;
8385 b2nm1 = *x*b2n+c*b2nm1;
8389 a2n = a2nm1+cma*a2n;
8390 b2n = b2nm1+cma*b2n;
8392 if(fabs(an0-am0) >= tol*an0) goto S260;
8394 *ans = 0.5e0+(0.5e0-*qans);
8398 GENERAL TEMME EXPANSION
8400 if(fabs(s) <= 2.0e0*e && *a*e*e > 3.28e-3) goto S430;
8403 w = 0.5e0*erfc1(&K1,&T3);
8406 if(l < 1.0e0) z = -z;
8408 if(T4 < 0) goto S280;
8409 else if(T4 == 0) goto S290;
8412 if(fabs(s) <= 1.e-3) goto S340;
8413 c0 = ((((((((((((d0[12]*z+d0[11])*z+d0[10])*z+d0[9])*z+d0[8])*z+d0[7])*z+d0[
8414 6])*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
8415 c1 = (((((((((((d1[11]*z+d1[10])*z+d1[9])*z+d1[8])*z+d1[7])*z+d1[6])*z+d1[5]
8416 )*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8417 c2 = (((((((((d2[9]*z+d2[8])*z+d2[7])*z+d2[6])*z+d2[5])*z+d2[4])*z+d2[3])*z+
8418 d2[2])*z+d2[1])*z+d2[0])*z+d20;
8419 c3 = (((((((d3[7]*z+d3[6])*z+d3[5])*z+d3[4])*z+d3[3])*z+d3[2])*z+d3[1])*z+
8421 c4 = (((((d4[5]*z+d4[4])*z+d4[3])*z+d4[2])*z+d4[1])*z+d4[0])*z+d40;
8422 c5 = (((d5[3]*z+d5[2])*z+d5[1])*z+d5[0])*z+d50;
8423 c6 = (d6[1]*z+d6[0])*z+d60;
8424 t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
8427 c0 = (((((d0[5]*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
8428 c1 = (((d1[3]*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8433 t = ((d0[2]*z+d0[1])*z+d0[0])*z-third;
8435 if(l < 1.0e0) goto S320;
8436 *qans = c*(w+rt2pin*t/rta);
8437 *ans = 0.5e0+(0.5e0-*qans);
8440 *ans = c*(w-rt2pin*t/rta);
8441 *qans = 0.5e0+(0.5e0-*ans);
8445 TEMME EXPANSION FOR L = 1
8447 if(*a*e*e > 3.28e-3) goto S430;
8448 c = 0.5e0+(0.5e0-y);
8449 w = (0.5e0-sqrt(y)*(0.5e0+(0.5e0-y/3.0e0))/rtpi)/c;
8452 if(l < 1.0e0) z = -z;
8454 if(T5 < 0) goto S340;
8455 else if(T5 == 0) goto S350;
8458 c0 = ((((((d0[6]*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-
8460 c1 = (((((d1[5]*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8461 c2 = ((((d2[4]*z+d2[3])*z+d2[2])*z+d2[1])*z+d2[0])*z+d20;
8462 c3 = (((d3[3]*z+d3[2])*z+d3[1])*z+d3[0])*z+d30;
8463 c4 = (d4[1]*z+d4[0])*z+d40;
8464 c5 = (d5[1]*z+d5[0])*z+d50;
8466 t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
8469 c0 = (d0[1]*z+d0[0])*z-third;
8471 t = (d20*u+c1)*u+c0;
8488 if(*x >= 0.25e0) goto S400;
8491 *qans = 0.5e0+(0.5e0-*ans);
8495 *qans = erfc1(&K2,&T7);
8496 *ans = 0.5e0+(0.5e0-*qans);
8499 if(fabs(s) <= 2.0e0*e) goto S430;
8501 if(*x <= *a) goto S370;
8510 double gsumln(double *a,double *b)
8512 -----------------------------------------------------------------------
8513 EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
8514 FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2
8515 -----------------------------------------------------------------------
8518 static double gsumln,x,T1,T2;
8521 .. Executable Statements ..
8524 if(x > 0.25e0) goto S10;
8526 gsumln = gamln1(&T1);
8529 if(x > 1.25e0) goto S20;
8530 gsumln = gamln1(&x)+alnrel(&x);
8534 gsumln = gamln1(&T2)+log(x*(1.0e0+x));
8537 double psi(double *xx)
8539 ---------------------------------------------------------------------
8541 EVALUATION OF THE DIGAMMA FUNCTION
8545 PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
8548 THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
8549 APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
8550 CODY, STRECOK AND THACHER.
8552 ---------------------------------------------------------------------
8553 PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
8554 PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
8556 ---------------------------------------------------------------------
8559 static double dx0 = 1.461632144968362341262659542325721325e0;
8560 static double piov4 = .785398163397448e0;
8561 static double p1[7] = {
8562 .895385022981970e-02,.477762828042627e+01,.142441585084029e+03,
8563 .118645200713425e+04,.363351846806499e+04,.413810161269013e+04,
8564 .130560269827897e+04
8566 static double p2[4] = {
8567 -.212940445131011e+01,-.701677227766759e+01,-.448616543918019e+01,
8568 -.648157123766197e+00
8570 static double q1[6] = {
8571 .448452573429826e+02,.520752771467162e+03,.221000799247830e+04,
8572 .364127349079381e+04,.190831076596300e+04,.691091682714533e-05
8574 static double q2[4] = {
8575 .322703493791143e+02,.892920700481861e+02,.546117738103215e+02,
8576 .777788548522962e+01
8580 static double psi,aug,den,sgn,upper,w,x,xmax1,xmx0,xsmall,z;
8581 static int i,m,n,nq;
8584 .. Executable Statements ..
8587 ---------------------------------------------------------------------
8588 MACHINE DEPENDENT CONSTANTS ...
8589 XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
8590 WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED
8591 AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
8592 ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
8593 PSI MAY BE REPRESENTED AS ALOG(X).
8594 XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
8595 MAY BE REPRESENTED BY 1/X.
8596 ---------------------------------------------------------------------
8598 xmax1 = ipmpar(&K1);
8599 xmax1 = fifdmin1(xmax1,1.0e0/spmpar(&K2));
8603 if(x >= 0.5e0) goto S50;
8605 ---------------------------------------------------------------------
8606 X .LT. 0.5, USE REFLECTION FORMULA
8607 PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
8608 ---------------------------------------------------------------------
8610 if(fabs(x) > xsmall) goto S10;
8611 if(x == 0.0e0) goto S100;
8613 ---------------------------------------------------------------------
8614 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE
8616 ---------------------------------------------------------------------
8622 ---------------------------------------------------------------------
8623 REDUCTION OF ARGUMENT FOR COTAN
8624 ---------------------------------------------------------------------
8628 if(w > 0.0e0) goto S20;
8633 ---------------------------------------------------------------------
8634 MAKE AN ERROR EXIT IF X .LE. -XMAX1
8635 ---------------------------------------------------------------------
8637 if(w >= xmax1) goto S100;
8640 nq = fifidint(w*4.0e0);
8641 w = 4.0e0*(w-(double)nq*.25e0);
8643 ---------------------------------------------------------------------
8644 W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X.
8645 ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
8646 QUADRANT AND DETERMINE SIGN
8647 ---------------------------------------------------------------------
8650 if(n+n != nq) w = 1.0e0-w;
8653 if(m+m != n) sgn = -sgn;
8655 ---------------------------------------------------------------------
8656 DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X)
8657 ---------------------------------------------------------------------
8662 if(m != n) goto S30;
8664 ---------------------------------------------------------------------
8665 CHECK FOR SINGULARITY
8666 ---------------------------------------------------------------------
8668 if(z == 0.0e0) goto S100;
8670 ---------------------------------------------------------------------
8671 USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
8672 SIN/COS AS A SUBSTITUTE FOR TAN
8673 ---------------------------------------------------------------------
8675 aug = sgn*(cos(z)/sin(z)*4.0e0);
8678 aug = sgn*(sin(z)/cos(z)*4.0e0);
8682 if(x > 3.0e0) goto S70;
8684 ---------------------------------------------------------------------
8686 ---------------------------------------------------------------------
8690 for(i=1; i<=5; i++) {
8691 den = (den+q1[i-1])*x;
8692 upper = (upper+p1[i+1-1])*x;
8694 den = (upper+p1[6])/(den+q1[5]);
8700 ---------------------------------------------------------------------
8701 IF X .GE. XMAX1, PSI = LN(X)
8702 ---------------------------------------------------------------------
8704 if(x >= xmax1) goto S90;
8706 ---------------------------------------------------------------------
8707 3.0 .LT. X .LT. XMAX1
8708 ---------------------------------------------------------------------
8713 for(i=1; i<=3; i++) {
8714 den = (den+q2[i-1])*w;
8715 upper = (upper+p2[i+1-1])*w;
8717 aug = upper/(den+q2[3])-0.5e0/x+aug;
8723 ---------------------------------------------------------------------
8725 ---------------------------------------------------------------------
8730 double rcomp(double *a,double *x)
8733 EVALUATION OF EXP(-X)*X**A/GAMMA(A)
8735 RT2PIN = 1/SQRT(2*PI)
8739 static double rt2pin = .398942280401433e0;
8740 static double rcomp,t,t1,u;
8743 .. Executable Statements ..
8746 if(*a >= 20.0e0) goto S20;
8748 if(*a >= 1.0e0) goto S10;
8749 rcomp = *a*exp(t)*(1.0e0+gam1(a));
8752 rcomp = exp(t)/Xgamm(a);
8756 if(u == 0.0e0) return rcomp;
8757 t = pow(1.0e0/ *a,2.0);
8758 t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
8759 t1 -= (*a*rlog(&u));
8760 rcomp = rt2pin*sqrt(*a)*exp(t1);
8763 double rexp(double *x)
8765 -----------------------------------------------------------------------
8766 EVALUATION OF THE FUNCTION EXP(X) - 1
8767 -----------------------------------------------------------------------
8770 static double p1 = .914041914819518e-09;
8771 static double p2 = .238082361044469e-01;
8772 static double q1 = -.499999999085958e+00;
8773 static double q2 = .107141568980644e+00;
8774 static double q3 = -.119041179760821e-01;
8775 static double q4 = .595130811860248e-03;
8776 static double rexp,w;
8779 .. Executable Statements ..
8781 if(fabs(*x) > 0.15e0) goto S10;
8782 rexp = *x*(((p2**x+p1)**x+1.0e0)/((((q4**x+q3)**x+q2)**x+q1)**x+1.0e0));
8786 if(*x > 0.0e0) goto S20;
8787 rexp = w-0.5e0-0.5e0;
8790 rexp = w*(0.5e0+(0.5e0-1.0e0/w));
8793 double rlog(double *x)
8796 COMPUTATION OF X - 1 - LN(X)
8800 static double a = .566749439387324e-01;
8801 static double b = .456512608815524e-01;
8802 static double p0 = .333333333333333e+00;
8803 static double p1 = -.224696413112536e+00;
8804 static double p2 = .620886815375787e-02;
8805 static double q1 = -.127408923933623e+01;
8806 static double q2 = .354508718369557e+00;
8807 static double rlog,r,t,u,w,w1;
8810 .. Executable Statements ..
8812 if(*x < 0.61e0 || *x > 1.57e0) goto S40;
8813 if(*x < 0.82e0) goto S10;
8814 if(*x > 1.18e0) goto S20;
8835 w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
8836 rlog = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
8843 double rlog1(double *x)
8845 -----------------------------------------------------------------------
8846 EVALUATION OF THE FUNCTION X - LN(1 + X)
8847 -----------------------------------------------------------------------
8850 static double a = .566749439387324e-01;
8851 static double b = .456512608815524e-01;
8852 static double p0 = .333333333333333e+00;
8853 static double p1 = -.224696413112536e+00;
8854 static double p2 = .620886815375787e-02;
8855 static double q1 = -.127408923933623e+01;
8856 static double q2 = .354508718369557e+00;
8857 static double rlog1,h,r,t,w,w1;
8860 .. Executable Statements ..
8862 if(*x < -0.39e0 || *x > 0.57e0) goto S40;
8863 if(*x < -0.18e0) goto S10;
8864 if(*x > 0.18e0) goto S20;
8877 h = 0.75e0**x-0.25e0;
8885 w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
8886 rlog1 = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
8893 double spmpar(int *i)
8895 -----------------------------------------------------------------------
8897 SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
8898 THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
8899 I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
8900 SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
8901 ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
8903 SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
8905 SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
8907 SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
8909 -----------------------------------------------------------------------
8911 ALFRED H. MORRIS, JR.
8912 NAVAL SURFACE WARFARE CENTER
8914 -----------------------------------------------------------------------
8915 -----------------------------------------------------------------------
8916 MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE
8917 CONSTANTS FOR THE COMPUTER BEING USED. THIS MODIFICATION WAS
8918 MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION
8919 -----------------------------------------------------------------------
8926 static double spmpar,b,binv,bm1,one,w,z;
8927 static int emax,emin,ibeta,m;
8930 .. Executable Statements ..
8932 if(*i > 1) goto S10;
8935 spmpar = pow(b,(double)(1-m));
8938 if(*i > 2) goto S20;
8943 w = pow(b,(double)(emin+2));
8944 spmpar = w*binv*binv*binv;
8947 ibeta = ipmpar(&K1);
8953 z = pow(b,(double)(m-1));
8954 w = ((z-one)*b+bm1)/(b*z);
8955 z = pow(b,(double)(emax-2));
8959 double stvaln(double *p)
8961 **********************************************************************
8963 double stvaln(double *p)
8964 STarting VALue for Neton-Raphon
8965 calculation of Normal distribution Inverse
8971 Returns X such that CUMNOR(X) = P, i.e., the integral from -
8972 infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
8978 P --> The probability whose normal deviate is sought.
8979 P is DOUBLE PRECISION
8985 The rational function on page 95 of Kennedy and Gentle,
8986 Statistical Computing, Marcel Dekker, NY , 1980.
8988 **********************************************************************
8991 static double xden[5] = {
8992 0.993484626060e-1,0.588581570495e0,0.531103462366e0,0.103537752850e0,
8995 static double xnum[5] = {
8996 -0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1,
9000 static double stvaln,sign,y,z;
9003 .. Executable Statements ..
9005 if(!(*p <= 0.5e0)) goto S10;
9013 y = sqrt(-(2.0e0*log(z)));
9014 stvaln = y+devlpl(xnum,&K1,&y)/devlpl(xden,&K1,&y);
9015 stvaln = sign*stvaln;
9018 /************************************************************************
9020 Truncates a double precision number to an integer and returns the
9022 ************************************************************************/
9023 double fifdint(double a)
9024 /* a - number to be truncated */
9028 return (double)(temp);
9030 /************************************************************************
9032 returns the maximum of two numbers a and b
9033 ************************************************************************/
9034 double fifdmax1(double a,double b)
9035 /* a - first number */
9036 /* b - second number */
9038 if (a < b) return b;
9041 /************************************************************************
9043 returns the minimum of two numbers a and b
9044 ************************************************************************/
9045 double fifdmin1(double a,double b)
9046 /* a - first number */
9047 /* b - second number */
9049 if (a < b) return a;
9052 /************************************************************************
9054 transfers the sign of the variable "sign" to the variable "mag"
9055 ************************************************************************/
9056 double fifdsign(double mag,double sign)
9057 /* mag - magnitude */
9058 /* sign - sign to be transfered */
9060 if (mag < 0) mag = -mag;
9061 if (sign < 0) mag = -mag;
9065 /************************************************************************
9067 Truncates a double precision number to a long integer
9068 ************************************************************************/
9069 long fifidint(double a)
9070 /* a - number to be truncated */
9074 /************************************************************************
9076 returns the modulo of a and b
9077 ************************************************************************/
9078 long fifmod(long a,long b)
9080 /* b - denominator */
9084 /************************************************************************
9086 Prints msg to standard error and then exits
9087 ************************************************************************/
9088 void ftnstop(char* msg)
9089 /* msg - error message */
9091 if (msg != NULL) fprintf(stderr,"%s\n",msg);
9092 exit(EXIT_FAILURE); /* EXIT_FAILURE from stdlib.h, or use an int */