7 static void E0000(int,int*,double*,double*,unsigned long*,
8 unsigned long*,double*,double*,double*,
9 double*,double*,double*,double*);
10 static void E0001(int,int*,double*,double*,double*,double*,
11 unsigned long*,unsigned long*,double*,double*,
15 * A comment about ints and longs - whether ints or longs are used should
16 * make no difference, but where double r-values are assigned to ints the
17 * r-value is cast converted to a long, which is then assigned to the int
18 * to be compatible with the operation of fifidint.
21 -----------------------------------------------------------------------
23 COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
27 IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
28 LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
30 -----------------------------------------------------------------------
32 double algdiv(double *a,double *b)
34 static double c0 = .833333333333333e-01;
35 static double c1 = -.277777777760991e-02;
36 static double c2 = .793650666825390e-03;
37 static double c3 = -.595202931351870e-03;
38 static double c4 = .837308034031215e-03;
39 static double c5 = -.165322962780713e-02;
40 static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;
43 .. Executable Statements ..
45 if(*a <= *b) goto S10;
58 SET SN = (1 - X**N)/(1 - X)
65 s11 = 1.0e0+(x+x2*s9);
67 SET W = DEL(B) - DEL(A + B)
69 t = pow(1.0e0/ *b,2.0);
70 w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
77 v = *a*(log(*b)-1.0e0);
85 double alngam(double *x)
87 **********************************************************************
89 double alngam(double *x)
90 double precision LN of the GAMma function
96 Returns the natural logarithm of GAMMA(X).
102 X --> value at which scaled log gamma is to be returned
103 X is DOUBLE PRECISION
109 If X .le. 6.0, then use recursion to get X below 3
110 then apply rational approximation number 5236 of
111 Hart et al, Computer Approximations, John Wiley and
114 If X .gt. 6.0, then use recursion to get X to at least 12 and
115 then use formula 5423 of the same source.
117 **********************************************************************
120 #define hln2pi 0.91893853320467274178e0
121 static double coef[5] = {
122 0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
123 -0.594997310889e-3,0.8065880899e-3
125 static double scoefd[4] = {
126 0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
127 0.1000000000000000000e1
129 static double scoefn[9] = {
130 0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
131 0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
132 0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
137 static double alngam,offset,prod,xx;
139 static double T2,T4,T6;
142 .. Executable Statements ..
144 if(!(*x <= 6.0e0)) goto S70;
147 if(!(*x > 3.0e0)) goto S30;
149 if(!(xx > 3.0e0)) goto S20;
155 if(!(*x < 2.0e0)) goto S60;
157 if(!(xx < 2.0e0)) goto S50;
165 alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
167 COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
170 alngam = log(alngam);
175 IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
177 n = fifidint(12.0e0-*x);
178 if(!(n > 0)) goto S90;
180 for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));
190 T6 = 1.0e0/pow(xx,2.0);
191 alngam = devlpl(coef,&K5,&T6)/xx;
192 alngam += (offset+(xx-0.5e0)*log(xx)-xx);
197 double alnrel(double *a)
199 -----------------------------------------------------------------------
200 EVALUATION OF THE FUNCTION LN(1 + A)
201 -----------------------------------------------------------------------
204 static double p1 = -.129418923021993e+01;
205 static double p2 = .405303492862024e+00;
206 static double p3 = -.178874546012214e-01;
207 static double q1 = -.162752256355323e+01;
208 static double q2 = .747811014037616e+00;
209 static double q3 = -.845104217945565e-01;
210 static double alnrel,t,t2,w,x;
213 .. Executable Statements ..
215 if(fabs(*a) > 0.375e0) goto S10;
218 w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);
226 double apser(double *a,double *b,double *x,double *eps)
228 -----------------------------------------------------------------------
229 APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
230 A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
231 A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
232 -----------------------------------------------------------------------
235 static double g = .577215664901533e0;
236 static double apser,aj,bx,c,j,s,t,tol;
239 .. Executable Statements ..
243 if(*b**eps > 2.e-2) goto S10;
244 c = log(*x)+psi(b)+g+t;
249 tol = 5.0e0**eps*fabs(c);
257 if(fabs(aj) > tol) goto S30;
261 double basym(double *a,double *b,double *lambda,double *eps)
263 -----------------------------------------------------------------------
264 ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
265 LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED.
266 IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
267 A AND B ARE GREATER THAN OR EQUAL TO 15.
268 -----------------------------------------------------------------------
271 static double e0 = 1.12837916709551e0;
272 static double e1 = .353553390593274e0;
275 ------------------------
276 ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
277 ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
278 THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
279 ------------------------
282 ------------------------
285 static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
287 static int i,im1,imj,j,m,mm1,mmj,n,np1;
288 static double a0[21],b0[21],c[21],d[21],T1,T2;
291 .. Executable Statements ..
294 if(*a >= *b) goto S10;
296 r0 = 1.0e0/(1.0e0+h);
298 w0 = 1.0e0/sqrt(*a*(1.0e0+h));
302 r0 = 1.0e0/(1.0e0+h);
304 w0 = 1.0e0/sqrt(*b*(1.0e0+h));
308 f = *a*rlog1(&T1)+*b*rlog1(&T2);
310 if(t == 0.0e0) return basym;
314 a0[0] = 2.0e0/3.0e0*r1;
315 c[0] = -(0.5e0*a0[0]);
317 j0 = 0.5e0/e0*erfc1(&K3,&z0);
326 for(n=2; n<=num; n+=2) {
328 a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
331 a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
332 for(i=n; i<=np1; i++) {
333 r = -(0.5e0*((double)i+1.0e0));
335 for(m=2; m<=i; m++) {
338 for(j=1; j<=mm1; j++) {
340 bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
342 b0[m-1] = r*a0[m-1]+bsum/(double)m;
344 c[i-1] = b0[i-1]/((double)i+1.0e0);
347 for(j=1; j<=im1; j++) {
349 dsum += (d[imj-1]*c[j-1]);
351 d[i-1] = -(dsum+c[i-1]);
353 j0 = e1*znm1+((double)n-1.0e0)*j0;
354 j1 = e1*zn+(double)n*j1;
362 if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
365 u = exp(-bcorr(a,b));
369 double bcorr(double *a0,double *b0)
371 -----------------------------------------------------------------------
373 EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE
374 LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
375 IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8.
377 -----------------------------------------------------------------------
380 static double c0 = .833333333333333e-01;
381 static double c1 = -.277777777760991e-02;
382 static double c2 = .793650666825390e-03;
383 static double c3 = -.595202931351870e-03;
384 static double c4 = .837308034031215e-03;
385 static double c5 = -.165322962780713e-02;
386 static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;
389 .. Executable Statements ..
391 a = fifdmin1(*a0,*b0);
392 b = fifdmax1(*a0,*b0);
398 SET SN = (1 - X**N)/(1 - X)
401 s5 = 1.0e0+(x+x2*s3);
402 s7 = 1.0e0+(x+x2*s5);
403 s9 = 1.0e0+(x+x2*s7);
404 s11 = 1.0e0+(x+x2*s9);
406 SET W = DEL(B) - DEL(A + B)
408 t = pow(1.0e0/b,2.0);
409 w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
414 t = pow(1.0e0/a,2.0);
415 bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w;
418 double betaln(double *a0,double *b0)
420 -----------------------------------------------------------------------
421 EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
422 -----------------------------------------------------------------------
424 --------------------------
427 static double e = .918938533204673e0;
428 static double betaln,a,b,c,h,u,v,w,z;
433 .. Executable Statements ..
435 a = fifdmin1(*a0,*b0);
436 b = fifdmax1(*a0,*b0);
437 if(a >= 8.0e0) goto S100;
438 if(a >= 1.0e0) goto S20;
440 -----------------------------------------------------------------------
441 PROCEDURE WHEN A .LT. 1
442 -----------------------------------------------------------------------
444 if(b >= 8.0e0) goto S10;
446 betaln = gamln(&a)+(gamln(&b)-gamln(&T1));
449 betaln = gamln(&a)+algdiv(&a,&b);
453 -----------------------------------------------------------------------
454 PROCEDURE WHEN 1 .LE. A .LT. 8
455 -----------------------------------------------------------------------
457 if(a > 2.0e0) goto S40;
458 if(b > 2.0e0) goto S30;
459 betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b);
463 if(b < 8.0e0) goto S60;
464 betaln = gamln(&a)+algdiv(&a,&b);
468 REDUCTION OF A WHEN B .LE. 1000
470 if(b > 1000.0e0) goto S80;
471 n = (long)(a - 1.0e0);
473 for(i=1; i<=n; i++) {
479 if(b < 8.0e0) goto S60;
480 betaln = w+gamln(&a)+algdiv(&a,&b);
484 REDUCTION OF B WHEN B .LT. 8
486 n = (long)(b - 1.0e0);
488 for(i=1; i<=n; i++) {
492 betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
496 REDUCTION OF A WHEN B .GT. 1000
498 n = (long)(a - 1.0e0);
500 for(i=1; i<=n; i++) {
502 w *= (a/(1.0e0+a/b));
504 betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
508 -----------------------------------------------------------------------
509 PROCEDURE WHEN A .GE. 8
510 -----------------------------------------------------------------------
515 u = -((a-0.5e0)*log(c));
517 if(u <= v) goto S110;
518 betaln = -(0.5e0*log(b))+e+w-v-u;
521 betaln = -(0.5e0*log(b))+e+w-u-v;
524 double bfrac(double *a,double *b,double *x,double *y,double *lambda,
527 -----------------------------------------------------------------------
528 CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
529 IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B.
530 -----------------------------------------------------------------------
533 static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;
536 .. Executable Statements ..
538 bfrac = brcomp(a,b,x,y);
539 if(bfrac == 0.0e0) return bfrac;
542 c1 = 1.0e0+1.0e0/ *a;
553 CONTINUED FRACTION CALCULATION
559 alpha = p*(p+c0)*e*e*(w**x);
560 e = (1.0e0+t)/(c1+t+t);
561 beta = n+w/s+e*(c+n*yp1);
565 UPDATE AN, BN, ANP1, AND BNP1
567 t = alpha*an+beta*anp1;
570 t = alpha*bn+beta*bnp1;
575 if(fabs(r-r0) <= *eps*r) goto S20;
577 RESCALE AN, BN, ANP1, AND BNP1
591 void bgrat(double *a,double *b,double *x,double *y,double *w,
592 double *eps,int *ierr)
594 -----------------------------------------------------------------------
595 ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
596 THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
597 THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED.
598 IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
599 -----------------------------------------------------------------------
602 static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
604 static double c[30],d[30],T1;
607 .. Executable Statements ..
609 bm1 = *b-0.5e0-0.5e0;
611 if(*y > 0.375e0) goto S10;
619 if(*b*z == 0.0e0) goto S70;
621 COMPUTATION OF THE EXPANSION
622 SET R = EXP(-Z)*Z**B/GAMMA(B)
624 r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
625 r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
626 u = algdiv(b,a)+*b*log(nu);
628 if(u == 0.0e0) goto S70;
629 grat1(b,&z,&r,&p,&q,eps);
630 v = 0.25e0*pow(1.0e0/nu,2.0);
637 for(n=1; n<=30; n++) {
639 j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
642 cn /= (n2*(n2+1.0e0));
648 for(i=1; i<=nm1; i++) {
649 s += (coef*c[i-1]*d[n-i-1]);
653 d[n-1] = bm1*cn+s/(double)n;
656 if(sum <= 0.0e0) goto S70;
657 if(fabs(dj) <= *eps*(sum+l)) goto S60;
668 THE EXPANSION CANNOT BE COMPUTED
673 double bpser(double *a,double *b,double *x,double *eps)
675 -----------------------------------------------------------------------
676 POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
677 OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED.
678 -----------------------------------------------------------------------
681 static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
685 .. Executable Statements ..
688 if(*x == 0.0e0) return bpser;
690 -----------------------------------------------------------------------
691 COMPUTE THE FACTOR X**A/(A*BETA(A,B))
692 -----------------------------------------------------------------------
694 a0 = fifdmin1(*a,*b);
695 if(a0 < 1.0e0) goto S10;
696 z = *a*log(*x)-betaln(a,b);
700 b0 = fifdmax1(*a,*b);
701 if(b0 >= 8.0e0) goto S90;
702 if(b0 > 1.0e0) goto S40;
704 PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
707 if(bpser == 0.0e0) return bpser;
709 if(apb > 1.0e0) goto S20;
710 z = 1.0e0+gam1(&apb);
714 z = (1.0e0+gam1(&u))/apb;
716 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
717 bpser *= (c*(*b/apb));
721 PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
724 m = (long)(b0 - 1.0e0);
727 for(i=1; i<=m; i++) {
736 if(apb > 1.0e0) goto S70;
737 t = 1.0e0+gam1(&apb);
741 t = (1.0e0+gam1(&u))/apb;
743 bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
747 PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
749 u = gamln1(&a0)+algdiv(&a0,&b0);
751 bpser = a0/ *a*exp(z);
753 if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
755 -----------------------------------------------------------------------
757 -----------------------------------------------------------------------
764 c *= ((0.5e0+(0.5e0-*b/n))**x);
767 if(fabs(w) > tol) goto S110;
768 bpser *= (1.0e0+*a*sum);
771 void bratio(double *a,double *b,double *x,double *y,double *w,
772 double *w1,int *ierr)
774 -----------------------------------------------------------------------
776 EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
780 IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
781 AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES
786 IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
787 IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
788 W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
789 THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
790 ONE OF THE FOLLOWING VALUES ...
792 IERR = 1 IF A OR B IS NEGATIVE
793 IERR = 2 IF A = B = 0
794 IERR = 3 IF X .LT. 0 OR X .GT. 1
795 IERR = 4 IF Y .LT. 0 OR Y .GT. 1
796 IERR = 5 IF X + Y .NE. 1
797 IERR = 6 IF X = A = 0
798 IERR = 7 IF Y = B = 0
801 WRITTEN BY ALFRED H. MORRIS, JR.
802 NAVAL SURFACE WARFARE CENTER
805 -----------------------------------------------------------------------
809 static double a0,b0,eps,lambda,t,x0,y0,z;
810 static int ierr1,ind,n;
811 static double T2,T3,T4,T5;
814 .. Executable Statements ..
817 ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
818 FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
822 if(*a < 0.0e0 || *b < 0.0e0) goto S270;
823 if(*a == 0.0e0 && *b == 0.0e0) goto S280;
824 if(*x < 0.0e0 || *x > 1.0e0) goto S290;
825 if(*y < 0.0e0 || *y > 1.0e0) goto S300;
826 z = *x+*y-0.5e0-0.5e0;
827 if(fabs(z) > 3.0e0*eps) goto S310;
829 if(*x == 0.0e0) goto S210;
830 if(*y == 0.0e0) goto S230;
831 if(*a == 0.0e0) goto S240;
832 if(*b == 0.0e0) goto S220;
833 eps = fifdmax1(eps,1.e-15);
834 if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
840 if(fifdmin1(a0,b0) > 1.0e0) goto S40;
842 PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
844 if(*x <= 0.5e0) goto S10;
851 if(b0 < fifdmin1(eps,eps*a0)) goto S90;
852 if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
853 if(fifdmax1(a0,b0) > 1.0e0) goto S20;
854 if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
855 if(pow(x0,a0) <= 0.9e0) goto S110;
856 if(x0 >= 0.3e0) goto S120;
860 if(b0 <= 1.0e0) goto S110;
861 if(x0 >= 0.3e0) goto S120;
862 if(x0 >= 0.1e0) goto S30;
863 if(pow(x0*b0,a0) <= 0.7e0) goto S110;
865 if(b0 > 15.0e0) goto S150;
870 PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
872 if(*a > *b) goto S50;
873 lambda = *a-(*a+*b)**x;
876 lambda = (*a+*b)**y-*b;
878 if(lambda >= 0.0e0) goto S70;
884 lambda = fabs(lambda);
886 if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
887 if(b0 < 40.0e0) goto S160;
888 if(a0 > b0) goto S80;
889 if(a0 <= 100.0e0) goto S130;
890 if(lambda > 0.03e0*a0) goto S130;
893 if(b0 <= 100.0e0) goto S130;
894 if(lambda > 0.03e0*b0) goto S130;
898 EVALUATION OF THE APPROPRIATE ALGORITHM
900 *w = fpser(&a0,&b0,&x0,&eps);
901 *w1 = 0.5e0+(0.5e0-*w);
904 *w1 = apser(&a0,&b0,&x0,&eps);
905 *w = 0.5e0+(0.5e0-*w1);
908 *w = bpser(&a0,&b0,&x0,&eps);
909 *w1 = 0.5e0+(0.5e0-*w);
912 *w1 = bpser(&b0,&a0,&y0,&eps);
913 *w = 0.5e0+(0.5e0-*w1);
917 *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
918 *w1 = 0.5e0+(0.5e0-*w);
921 *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
925 bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
926 *w = 0.5e0+(0.5e0-*w1);
931 if(b0 != 0.0e0) goto S170;
935 *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
936 if(x0 > 0.7e0) goto S180;
937 *w += bpser(&a0,&b0,&x0,&eps);
938 *w1 = 0.5e0+(0.5e0-*w);
941 if(a0 > 15.0e0) goto S190;
943 *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
947 bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
948 *w1 = 0.5e0+(0.5e0-*w);
952 *w = basym(&a0,&b0,&lambda,&T5);
953 *w1 = 0.5e0+(0.5e0-*w);
957 TERMINATION OF THE PROCEDURE
959 if(*a == 0.0e0) goto S320;
965 if(*b == 0.0e0) goto S330;
978 PROCEDURE FOR A AND B .LT. 1.E-3*EPS
1008 double brcmp1(int *mu,double *a,double *b,double *x,double *y)
1010 -----------------------------------------------------------------------
1011 EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B))
1012 -----------------------------------------------------------------------
1015 static double Const = .398942280401433e0;
1016 static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
1020 CONST = 1/SQRT(2*PI)
1023 static double T1,T2,T3,T4;
1026 .. Executable Statements ..
1028 a0 = fifdmin1(*a,*b);
1029 if(a0 >= 8.0e0) goto S130;
1030 if(*x > 0.375e0) goto S10;
1036 if(*y > 0.375e0) goto S20;
1046 if(a0 < 1.0e0) goto S40;
1048 brcmp1 = esum(mu,&z);
1052 -----------------------------------------------------------------------
1053 PROCEDURE FOR A .LT. 1 OR B .LT. 1
1054 -----------------------------------------------------------------------
1056 b0 = fifdmax1(*a,*b);
1057 if(b0 >= 8.0e0) goto S120;
1058 if(b0 > 1.0e0) goto S70;
1060 ALGORITHM FOR B0 .LE. 1
1062 brcmp1 = esum(mu,&z);
1063 if(brcmp1 == 0.0e0) return brcmp1;
1065 if(apb > 1.0e0) goto S50;
1066 z = 1.0e0+gam1(&apb);
1070 z = (1.0e0+gam1(&u))/apb;
1072 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1073 brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
1077 ALGORITHM FOR 1 .LT. B0 .LT. 8
1080 n = (long)(b0 - 1.0e0);
1083 for(i=1; i<=n; i++) {
1092 if(apb > 1.0e0) goto S100;
1093 t = 1.0e0+gam1(&apb);
1097 t = (1.0e0+gam1(&u))/apb;
1099 brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
1103 ALGORITHM FOR B0 .GE. 8
1105 u = gamln1(&a0)+algdiv(&a0,&b0);
1107 brcmp1 = a0*esum(mu,&T3);
1111 -----------------------------------------------------------------------
1112 PROCEDURE FOR A .GE. 8 AND B .GE. 8
1113 -----------------------------------------------------------------------
1115 if(*a > *b) goto S140;
1118 y0 = 1.0e0/(1.0e0+h);
1119 lambda = *a-(*a+*b)**x;
1123 x0 = 1.0e0/(1.0e0+h);
1125 lambda = (*a+*b)**y-*b;
1128 if(fabs(e) > 0.6e0) goto S160;
1135 if(fabs(e) > 0.6e0) goto S180;
1143 brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
1146 double brcomp(double *a,double *b,double *x,double *y)
1148 -----------------------------------------------------------------------
1149 EVALUATION OF X**A*Y**B/BETA(A,B)
1150 -----------------------------------------------------------------------
1153 static double Const = .398942280401433e0;
1154 static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
1158 CONST = 1/SQRT(2*PI)
1161 static double T1,T2;
1164 .. Executable Statements ..
1167 if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
1168 a0 = fifdmin1(*a,*b);
1169 if(a0 >= 8.0e0) goto S130;
1170 if(*x > 0.375e0) goto S10;
1176 if(*y > 0.375e0) goto S20;
1186 if(a0 < 1.0e0) goto S40;
1192 -----------------------------------------------------------------------
1193 PROCEDURE FOR A .LT. 1 OR B .LT. 1
1194 -----------------------------------------------------------------------
1196 b0 = fifdmax1(*a,*b);
1197 if(b0 >= 8.0e0) goto S120;
1198 if(b0 > 1.0e0) goto S70;
1200 ALGORITHM FOR B0 .LE. 1
1203 if(brcomp == 0.0e0) return brcomp;
1205 if(apb > 1.0e0) goto S50;
1206 z = 1.0e0+gam1(&apb);
1210 z = (1.0e0+gam1(&u))/apb;
1212 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
1213 brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
1217 ALGORITHM FOR 1 .LT. B0 .LT. 8
1220 n = (long)(b0 - 1.0e0);
1223 for(i=1; i<=n; i++) {
1232 if(apb > 1.0e0) goto S100;
1233 t = 1.0e0+gam1(&apb);
1237 t = (1.0e0+gam1(&u))/apb;
1239 brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
1243 ALGORITHM FOR B0 .GE. 8
1245 u = gamln1(&a0)+algdiv(&a0,&b0);
1246 brcomp = a0*exp(z-u);
1250 -----------------------------------------------------------------------
1251 PROCEDURE FOR A .GE. 8 AND B .GE. 8
1252 -----------------------------------------------------------------------
1254 if(*a > *b) goto S140;
1257 y0 = 1.0e0/(1.0e0+h);
1258 lambda = *a-(*a+*b)**x;
1262 x0 = 1.0e0/(1.0e0+h);
1264 lambda = (*a+*b)**y-*b;
1267 if(fabs(e) > 0.6e0) goto S160;
1274 if(fabs(e) > 0.6e0) goto S180;
1280 z = exp(-(*a*u+*b*v));
1281 brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
1284 double bup(double *a,double *b,double *x,double *y,int *n,double *eps)
1286 -----------------------------------------------------------------------
1287 EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
1288 EPS IS THE TOLERANCE USED.
1289 -----------------------------------------------------------------------
1294 static double bup,ap1,apb,d,l,r,t,w;
1295 static int i,k,kp1,mu,nm1;
1298 .. Executable Statements ..
1301 OBTAIN THE SCALING FACTOR EXP(-MU) AND
1302 EXP(MU)*(X**A*Y**B/BETA(A,B))/A
1308 if(*n == 1 || *a < 1.0e0) goto S10;
1309 if(apb < 1.1e0*ap1) goto S10;
1310 mu = (long)(fabs(exparg(&K1)));
1311 k = (long)(exparg(&K2));
1316 bup = brcmp1(&mu,a,b,x,y)/ *a;
1317 if(*n == 1 || bup == 0.0e0) return bup;
1321 LET K BE THE INDEX OF THE MAXIMUM TERM
1324 if(*b <= 1.0e0) goto S50;
1325 if(*y > 1.e-4) goto S20;
1329 r = (*b-1.0e0)**x/ *y-*a;
1330 if(r < 1.0e0) goto S50;
1333 if(r < t) k = (long)(r);
1336 ADD THE INCREASING TERMS OF THE SERIES
1338 for(i=1; i<=k; i++) {
1340 d = (apb+l)/(ap1+l)**x*d;
1343 if(k == nm1) goto S70;
1346 ADD THE REMAINING TERMS OF THE SERIES
1349 for(i=kp1; i<=nm1; i++) {
1351 d = (apb+l)/(ap1+l)**x*d;
1353 if(d <= *eps*w) goto S70;
1357 TERMINATE THE PROCEDURE
1362 void cdfbet(int *which,double *p,double *q,double *x,double *y,
1363 double *a,double *b,int *status,double *bound)
1364 /**********************************************************************
1366 void cdfbet(int *which,double *p,double *q,double *x,double *y,
1367 double *a,double *b,int *status,double *bound)
1369 Cumulative Distribution Function
1376 Calculates any one parameter of the beta distribution given
1377 values for the others.
1383 WHICH --> Integer indicating which of the next four argument
1384 values is to be calculated from the others.
1386 iwhich = 1 : Calculate P and Q from X,Y,A and B
1387 iwhich = 2 : Calculate X and Y from P,Q,A and B
1388 iwhich = 3 : Calculate A from P,Q,X,Y and B
1389 iwhich = 4 : Calculate B from P,Q,X,Y and A
1391 P <--> The integral from 0 to X of the chi-square
1393 Input range: [0, 1].
1396 Input range: [0, 1].
1399 X <--> Upper limit of integration of beta density.
1408 A <--> The first parameter of the beta density.
1409 Input range: (0, +infinity).
1410 Search range: [1D-100,1D100]
1412 B <--> The second parameter of the beta density.
1413 Input range: (0, +infinity).
1414 Search range: [1D-100,1D100]
1416 STATUS <-- 0 if calculation completed correctly
1417 -I if input parameter number I is out of range
1418 1 if answer appears to be lower than lowest
1420 2 if answer appears to be higher than greatest
1425 BOUND <-- Undefined if STATUS is 0
1427 Bound exceeded by parameter number I if STATUS
1430 Lower search bound if STATUS is 1.
1432 Upper search bound if STATUS is 2.
1438 Cumulative distribution function (P) is calculated directly by
1439 code associated with the following reference.
1441 DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant
1442 Digit Computation of the Incomplete Beta Function Ratios. ACM
1443 Trans. Math. Softw. 18 (1993), 360-373.
1445 Computation of other parameters involve a seach for a value that
1446 produces the desired value of P. The search relies on the
1447 monotinicity of P with the other parameter.
1453 The beta density is proportional to
1454 t^(A-1) * (1-t)^(B-1)
1456 **********************************************************************/
1459 #define atol 1.0e-50
1460 #define zero 1.0e-100
1464 static double K2 = 0.0e0;
1465 static double K3 = 1.0e0;
1466 static double K8 = 0.5e0;
1467 static double K9 = 5.0e0;
1468 static double fx,xhi,xlo,cum,ccum,xy,pq;
1469 static unsigned long qhi,qleft,qporq;
1470 static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
1473 .. Executable Statements ..
1478 if(!(*which < 1 || *which > 4)) goto S30;
1479 if(!(*which < 1)) goto S10;
1488 if(*which == 1) goto S70;
1492 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
1493 if(!(*p < 0.0e0)) goto S40;
1503 if(*which == 1) goto S110;
1507 if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
1508 if(!(*q < 0.0e0)) goto S80;
1518 if(*which == 2) goto S150;
1522 if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
1523 if(!(*x < 0.0e0)) goto S120;
1533 if(*which == 2) goto S190;
1537 if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
1538 if(!(*y < 0.0e0)) goto S160;
1548 if(*which == 3) goto S210;
1552 if(!(*a <= 0.0e0)) goto S200;
1558 if(*which == 4) goto S230;
1562 if(!(*b <= 0.0e0)) goto S220;
1568 if(*which == 1) goto S270;
1573 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
1574 if(!(pq < 0.0e0)) goto S240;
1584 if(*which == 2) goto S310;
1589 if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
1590 if(!(xy < 0.0e0)) goto S280;
1600 if(!(*which == 1)) qporq = *p <= *q;
1602 Select the minimum of P or Q
1609 cumbet(x,y,a,b,p,q);
1612 else if(2 == *which) {
1618 dstzr(&K2,&K3,&T4,&T5);
1619 if(!qporq) goto S340;
1621 dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
1624 if(!(*status == 1)) goto S330;
1625 cumbet(x,y,a,b,&cum,&ccum);
1627 dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
1634 dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
1637 if(!(*status == 1)) goto S360;
1638 cumbet(x,y,a,b,&cum,&ccum);
1640 dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
1645 if(!(*status == -1)) goto S400;
1646 if(!qleft) goto S380;
1657 else if(3 == *which) {
1666 dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
1668 dinvr(status,a,&fx,&qleft,&qhi);
1670 if(!(*status == 1)) goto S440;
1671 cumbet(x,y,a,b,&cum,&ccum);
1672 if(!qporq) goto S420;
1678 dinvr(status,a,&fx,&qleft,&qhi);
1681 if(!(*status == -1)) goto S470;
1682 if(!qleft) goto S450;
1693 else if(4 == *which) {
1702 dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
1704 dinvr(status,b,&fx,&qleft,&qhi);
1706 if(!(*status == 1)) goto S510;
1707 cumbet(x,y,a,b,&cum,&ccum);
1708 if(!qporq) goto S490;
1714 dinvr(status,b,&fx,&qleft,&qhi);
1717 if(!(*status == -1)) goto S540;
1718 if(!qleft) goto S520;
1736 void cdfbin(int *which,double *p,double *q,double *s,double *xn,
1737 double *pr,double *ompr,int *status,double *bound)
1738 /**********************************************************************
1740 void cdfbin(int *which,double *p,double *q,double *s,double *xn,
1741 double *pr,double *ompr,int *status,double *bound)
1743 Cumulative Distribution Function
1744 BINomial distribution
1750 Calculates any one parameter of the binomial
1751 distribution given values for the others.
1757 WHICH --> Integer indicating which of the next four argument
1758 values is to be calculated from the others.
1760 iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
1761 iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
1762 iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
1763 iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
1765 P <--> The cumulation from 0 to S of the binomial distribution.
1766 (Probablility of S or fewer successes in XN trials each
1767 with probability of success PR.)
1771 Input range: [0, 1].
1774 S <--> The number of successes observed.
1775 Input range: [0, XN]
1776 Search range: [0, XN]
1778 XN <--> The number of binomial trials.
1779 Input range: (0, +infinity).
1780 Search range: [1E-100, 1E100]
1782 PR <--> The probability of success in each binomial trial.
1791 STATUS <-- 0 if calculation completed correctly
1792 -I if input parameter number I is out of range
1793 1 if answer appears to be lower than lowest
1795 2 if answer appears to be higher than greatest
1798 4 if PR + OMPR .ne. 1
1800 BOUND <-- Undefined if STATUS is 0
1802 Bound exceeded by parameter number I if STATUS
1805 Lower search bound if STATUS is 1.
1807 Upper search bound if STATUS is 2.
1813 Formula 26.5.24 of Abramowitz and Stegun, Handbook of
1814 Mathematical Functions (1966) is used to reduce the binomial
1815 distribution to the cumulative incomplete beta distribution.
1817 Computation of other parameters involve a seach for a value that
1818 produces the desired value of P. The search relies on the
1819 monotinicity of P with the other parameter.
1822 **********************************************************************/
1824 #define atol 1.0e-50
1826 #define zero 1.0e-100
1830 static double K2 = 0.0e0;
1831 static double K3 = 0.5e0;
1832 static double K4 = 5.0e0;
1833 static double K11 = 1.0e0;
1834 static double fx,xhi,xlo,cum,ccum,pq,prompr;
1835 static unsigned long qhi,qleft,qporq;
1836 static double T5,T6,T7,T8,T9,T10,T12,T13;
1839 .. Executable Statements ..
1844 if(!(*which < 1 && *which > 4)) goto S30;
1845 if(!(*which < 1)) goto S10;
1854 if(*which == 1) goto S70;
1858 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
1859 if(!(*p < 0.0e0)) goto S40;
1869 if(*which == 1) goto S110;
1873 if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
1874 if(!(*q < 0.0e0)) goto S80;
1884 if(*which == 3) goto S130;
1888 if(!(*xn <= 0.0e0)) goto S120;
1894 if(*which == 2) goto S170;
1898 if(!(*s < 0.0e0 || (*which != 3 && *s > *xn))) goto S160;
1899 if(!(*s < 0.0e0)) goto S140;
1909 if(*which == 4) goto S210;
1913 if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S200;
1914 if(!(*pr < 0.0e0)) goto S180;
1924 if(*which == 4) goto S250;
1928 if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S240;
1929 if(!(*ompr < 0.0e0)) goto S220;
1939 if(*which == 1) goto S290;
1944 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S280;
1945 if(!(pq < 0.0e0)) goto S260;
1955 if(*which == 4) goto S330;
1960 if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S320;
1961 if(!(prompr < 0.0e0)) goto S300;
1971 if(!(*which == 1)) qporq = *p <= *q;
1973 Select the minimum of P or Q
1980 cumbin(s,xn,pr,ompr,p,q);
1983 else if(2 == *which) {
1990 dstinv(&K2,xn,&K3,&K3,&K4,&T5,&T6);
1992 dinvr(status,s,&fx,&qleft,&qhi);
1994 if(!(*status == 1)) goto S370;
1995 cumbin(s,xn,pr,ompr,&cum,&ccum);
1996 if(!qporq) goto S350;
2002 dinvr(status,s,&fx,&qleft,&qhi);
2005 if(!(*status == -1)) goto S400;
2006 if(!qleft) goto S380;
2017 else if(3 == *which) {
2026 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
2028 dinvr(status,xn,&fx,&qleft,&qhi);
2030 if(!(*status == 1)) goto S440;
2031 cumbin(s,xn,pr,ompr,&cum,&ccum);
2032 if(!qporq) goto S420;
2038 dinvr(status,xn,&fx,&qleft,&qhi);
2041 if(!(*status == -1)) goto S470;
2042 if(!qleft) goto S450;
2053 else if(4 == *which) {
2055 Calculating PR and OMPR
2059 dstzr(&K2,&K11,&T12,&T13);
2060 if(!qporq) goto S500;
2062 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
2065 if(!(*status == 1)) goto S490;
2066 cumbin(s,xn,pr,ompr,&cum,&ccum);
2068 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
2075 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
2078 if(!(*status == 1)) goto S520;
2079 cumbin(s,xn,pr,ompr,&cum,&ccum);
2081 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
2086 if(!(*status == -1)) goto S560;
2087 if(!qleft) goto S540;
2105 void cdfchi(int *which,double *p,double *q,double *x,double *df,
2106 int *status,double *bound)
2107 /**********************************************************************
2109 void cdfchi(int *which,double *p,double *q,double *x,double *df,
2110 int *status,double *bound)
2112 Cumulative Distribution Function
2113 CHI-Square distribution
2119 Calculates any one parameter of the chi-square
2120 distribution given values for the others.
2126 WHICH --> Integer indicating which of the next three argument
2127 values is to be calculated from the others.
2129 iwhich = 1 : Calculate P and Q from X and DF
2130 iwhich = 2 : Calculate X from P,Q and DF
2131 iwhich = 3 : Calculate DF from P,Q and X
2133 P <--> The integral from 0 to X of the chi-square
2135 Input range: [0, 1].
2138 Input range: (0, 1].
2141 X <--> Upper limit of integration of the non-central
2142 chi-square distribution.
2143 Input range: [0, +infinity).
2144 Search range: [0,1E100]
2146 DF <--> Degrees of freedom of the
2147 chi-square distribution.
2148 Input range: (0, +infinity).
2149 Search range: [ 1E-100, 1E100]
2151 STATUS <-- 0 if calculation completed correctly
2152 -I if input parameter number I is out of range
2153 1 if answer appears to be lower than lowest
2155 2 if answer appears to be higher than greatest
2158 10 indicates error returned from cumgam. See
2159 references in cdfgam
2161 BOUND <-- Undefined if STATUS is 0
2163 Bound exceeded by parameter number I if STATUS
2166 Lower search bound if STATUS is 1.
2168 Upper search bound if STATUS is 2.
2174 Formula 26.4.19 of Abramowitz and Stegun, Handbook of
2175 Mathematical Functions (1966) is used to reduce the chisqure
2176 distribution to the incomplete distribution.
2178 Computation of other parameters involve a seach for a value that
2179 produces the desired value of P. The search relies on the
2180 monotinicity of P with the other parameter.
2182 **********************************************************************/
2185 #define atol 1.0e-50
2186 #define zero 1.0e-100
2189 static double K2 = 0.0e0;
2190 static double K4 = 0.5e0;
2191 static double K5 = 5.0e0;
2192 static double fx,cum,ccum,pq,porq;
2193 static unsigned long qhi,qleft,qporq;
2194 static double T3,T6,T7,T8,T9,T10,T11;
2197 .. Executable Statements ..
2202 if(!(*which < 1 || *which > 3)) goto S30;
2203 if(!(*which < 1)) goto S10;
2212 if(*which == 1) goto S70;
2216 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
2217 if(!(*p < 0.0e0)) goto S40;
2227 if(*which == 1) goto S110;
2231 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
2232 if(!(*q <= 0.0e0)) goto S80;
2242 if(*which == 2) goto S130;
2246 if(!(*x < 0.0e0)) goto S120;
2252 if(*which == 3) goto S150;
2256 if(!(*df <= 0.0e0)) goto S140;
2262 if(*which == 1) goto S190;
2267 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
2268 if(!(pq < 0.0e0)) goto S160;
2278 if(*which == 1) goto S220;
2280 Select the minimum of P or Q
2283 if(!qporq) goto S200;
2304 else if(2 == *which) {
2312 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
2314 dinvr(status,x,&fx,&qleft,&qhi);
2316 if(!(*status == 1)) goto S270;
2317 cumchi(x,df,&cum,&ccum);
2318 if(!qporq) goto S240;
2324 if(!(fx+porq > 1.5e0)) goto S260;
2328 dinvr(status,x,&fx,&qleft,&qhi);
2331 if(!(*status == -1)) goto S300;
2332 if(!qleft) goto S280;
2343 else if(3 == *which) {
2352 dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
2354 dinvr(status,df,&fx,&qleft,&qhi);
2356 if(!(*status == 1)) goto S350;
2357 cumchi(x,df,&cum,&ccum);
2358 if(!qporq) goto S320;
2364 if(!(fx+porq > 1.5e0)) goto S340;
2368 dinvr(status,df,&fx,&qleft,&qhi);
2371 if(!(*status == -1)) goto S380;
2372 if(!qleft) goto S360;
2389 void cdfchn(int *which,double *p,double *q,double *x,double *df,
2390 double *pnonc,int *status,double *bound)
2391 /**********************************************************************
2393 void cdfchn(int *which,double *p,double *q,double *x,double *df,
2394 double *pnonc,int *status,double *bound)
2396 Cumulative Distribution Function
2397 Non-central Chi-Square
2403 Calculates any one parameter of the non-central chi-square
2404 distribution given values for the others.
2410 WHICH --> Integer indicating which of the next three argument
2411 values is to be calculated from the others.
2413 iwhich = 1 : Calculate P and Q from X and DF
2414 iwhich = 2 : Calculate X from P,DF and PNONC
2415 iwhich = 3 : Calculate DF from P,X and PNONC
2416 iwhich = 3 : Calculate PNONC from P,X and DF
2418 P <--> The integral from 0 to X of the non-central chi-square
2420 Input range: [0, 1-1E-16).
2423 Q is not used by this subroutine and is only included
2424 for similarity with other cdf* routines.
2426 X <--> Upper limit of integration of the non-central
2427 chi-square distribution.
2428 Input range: [0, +infinity).
2429 Search range: [0,1E100]
2431 DF <--> Degrees of freedom of the non-central
2432 chi-square distribution.
2433 Input range: (0, +infinity).
2434 Search range: [ 1E-100, 1E100]
2436 PNONC <--> Non-centrality parameter of the non-central
2437 chi-square distribution.
2438 Input range: [0, +infinity).
2439 Search range: [0,1E4]
2441 STATUS <-- 0 if calculation completed correctly
2442 -I if input parameter number I is out of range
2443 1 if answer appears to be lower than lowest
2445 2 if answer appears to be higher than greatest
2448 BOUND <-- Undefined if STATUS is 0
2450 Bound exceeded by parameter number I if STATUS
2453 Lower search bound if STATUS is 1.
2455 Upper search bound if STATUS is 2.
2461 Formula 26.4.25 of Abramowitz and Stegun, Handbook of
2462 Mathematical Functions (1966) is used to compute the cumulative
2463 distribution function.
2465 Computation of other parameters involve a seach for a value that
2466 produces the desired value of P. The search relies on the
2467 monotinicity of P with the other parameter.
2472 The computation time required for this routine is proportional
2473 to the noncentrality parameter (PNONC). Very large values of
2474 this parameter can consume immense computer resources. This is
2475 why the search range is bounded by 10,000.
2477 **********************************************************************/
2481 #define atol 1.0e-50
2482 #define zero 1.0e-100
2483 #define one ( 1.0e0 - 1.0e-16 )
2485 static double K1 = 0.0e0;
2486 static double K3 = 0.5e0;
2487 static double K4 = 5.0e0;
2488 static double fx,cum,ccum;
2489 static unsigned long qhi,qleft;
2490 static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;
2493 .. Executable Statements ..
2498 if(!(*which < 1 || *which > 4)) goto S30;
2499 if(!(*which < 1)) goto S10;
2508 if(*which == 1) goto S70;
2512 if(!(*p < 0.0e0 || *p > one)) goto S60;
2513 if(!(*p < 0.0e0)) goto S40;
2523 if(*which == 2) goto S90;
2527 if(!(*x < 0.0e0)) goto S80;
2533 if(*which == 3) goto S110;
2537 if(!(*df <= 0.0e0)) goto S100;
2543 if(*which == 4) goto S130;
2547 if(!(*pnonc < 0.0e0)) goto S120;
2560 cumchn(x,df,pnonc,p,q);
2563 else if(2 == *which) {
2571 dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
2573 dinvr(status,x,&fx,&qleft,&qhi);
2575 if(!(*status == 1)) goto S150;
2576 cumchn(x,df,pnonc,&cum,&ccum);
2578 dinvr(status,x,&fx,&qleft,&qhi);
2581 if(!(*status == -1)) goto S180;
2582 if(!qleft) goto S160;
2593 else if(3 == *which) {
2602 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
2604 dinvr(status,df,&fx,&qleft,&qhi);
2606 if(!(*status == 1)) goto S200;
2607 cumchn(x,df,pnonc,&cum,&ccum);
2609 dinvr(status,df,&fx,&qleft,&qhi);
2612 if(!(*status == -1)) goto S230;
2613 if(!qleft) goto S210;
2624 else if(4 == *which) {
2632 dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);
2634 dinvr(status,pnonc,&fx,&qleft,&qhi);
2636 if(!(*status == 1)) goto S250;
2637 cumchn(x,df,pnonc,&cum,&ccum);
2639 dinvr(status,pnonc,&fx,&qleft,&qhi);
2642 if(!(*status == -1)) goto S280;
2643 if(!qleft) goto S260;
2662 void cdff(int *which,double *p,double *q,double *f,double *dfn,
2663 double *dfd,int *status,double *bound)
2664 /**********************************************************************
2666 void cdff(int *which,double *p,double *q,double *f,double *dfn,
2667 double *dfd,int *status,double *bound)
2669 Cumulative Distribution Function
2676 Calculates any one parameter of the F distribution
2677 given values for the others.
2683 WHICH --> Integer indicating which of the next four argument
2684 values is to be calculated from the others.
2686 iwhich = 1 : Calculate P and Q from F,DFN and DFD
2687 iwhich = 2 : Calculate F from P,Q,DFN and DFD
2688 iwhich = 3 : Calculate DFN from P,Q,F and DFD
2689 iwhich = 4 : Calculate DFD from P,Q,F and DFN
2691 P <--> The integral from 0 to F of the f-density.
2695 Input range: (0, 1].
2698 F <--> Upper limit of integration of the f-density.
2699 Input range: [0, +infinity).
2700 Search range: [0,1E100]
2702 DFN < --> Degrees of freedom of the numerator sum of squares.
2703 Input range: (0, +infinity).
2704 Search range: [ 1E-100, 1E100]
2706 DFD < --> Degrees of freedom of the denominator sum of squares.
2707 Input range: (0, +infinity).
2708 Search range: [ 1E-100, 1E100]
2710 STATUS <-- 0 if calculation completed correctly
2711 -I if input parameter number I is out of range
2712 1 if answer appears to be lower than lowest
2714 2 if answer appears to be higher than greatest
2718 BOUND <-- Undefined if STATUS is 0
2720 Bound exceeded by parameter number I if STATUS
2723 Lower search bound if STATUS is 1.
2725 Upper search bound if STATUS is 2.
2731 Formula 26.6.2 of Abramowitz and Stegun, Handbook of
2732 Mathematical Functions (1966) is used to reduce the computation
2733 of the cumulative distribution function for the F variate to
2734 that of an incomplete beta.
2736 Computation of other parameters involve a seach for a value that
2737 produces the desired value of P. The search relies on the
2738 monotinicity of P with the other parameter.
2742 The value of the cumulative F distribution is not necessarily
2743 monotone in either degrees of freedom. There thus may be two
2744 values that provide a given CDF value. This routine assumes
2745 monotonicity and will find an arbitrary one of the two values.
2747 **********************************************************************/
2750 #define atol 1.0e-50
2751 #define zero 1.0e-100
2754 static double K2 = 0.0e0;
2755 static double K4 = 0.5e0;
2756 static double K5 = 5.0e0;
2757 static double pq,fx,cum,ccum;
2758 static unsigned long qhi,qleft,qporq;
2759 static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
2762 .. Executable Statements ..
2767 if(!(*which < 1 || *which > 4)) goto S30;
2768 if(!(*which < 1)) goto S10;
2777 if(*which == 1) goto S70;
2781 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
2782 if(!(*p < 0.0e0)) goto S40;
2792 if(*which == 1) goto S110;
2796 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
2797 if(!(*q <= 0.0e0)) goto S80;
2807 if(*which == 2) goto S130;
2811 if(!(*f < 0.0e0)) goto S120;
2817 if(*which == 3) goto S150;
2821 if(!(*dfn <= 0.0e0)) goto S140;
2827 if(*which == 4) goto S170;
2831 if(!(*dfd <= 0.0e0)) goto S160;
2837 if(*which == 1) goto S210;
2842 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
2843 if(!(pq < 0.0e0)) goto S180;
2853 if(!(*which == 1)) qporq = *p <= *q;
2855 Select the minimum of P or Q
2862 cumf(f,dfn,dfd,p,q);
2865 else if(2 == *which) {
2873 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
2875 dinvr(status,f,&fx,&qleft,&qhi);
2877 if(!(*status == 1)) goto S250;
2878 cumf(f,dfn,dfd,&cum,&ccum);
2879 if(!qporq) goto S230;
2885 dinvr(status,f,&fx,&qleft,&qhi);
2888 if(!(*status == -1)) goto S280;
2889 if(!qleft) goto S260;
2900 else if(3 == *which) {
2909 dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
2911 dinvr(status,dfn,&fx,&qleft,&qhi);
2913 if(!(*status == 1)) goto S320;
2914 cumf(f,dfn,dfd,&cum,&ccum);
2915 if(!qporq) goto S300;
2921 dinvr(status,dfn,&fx,&qleft,&qhi);
2924 if(!(*status == -1)) goto S350;
2925 if(!qleft) goto S330;
2936 else if(4 == *which) {
2945 dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15);
2947 dinvr(status,dfd,&fx,&qleft,&qhi);
2949 if(!(*status == 1)) goto S390;
2950 cumf(f,dfn,dfd,&cum,&ccum);
2951 if(!qporq) goto S370;
2957 dinvr(status,dfd,&fx,&qleft,&qhi);
2960 if(!(*status == -1)) goto S420;
2961 if(!qleft) goto S400;
2978 void cdffnc(int *which,double *p,double *q,double *f,double *dfn,
2979 double *dfd,double *phonc,int *status,double *bound)
2980 /**********************************************************************
2982 void cdffnc(int *which,double *p,double *q,double *f,double *dfn,
2983 double *dfd,double *phonc,int *status,double *bound)
2985 Cumulative Distribution Function
2986 Non-central F distribution
2992 Calculates any one parameter of the Non-central F
2993 distribution given values for the others.
2999 WHICH --> Integer indicating which of the next five argument
3000 values is to be calculated from the others.
3002 iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
3003 iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
3004 iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
3005 iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
3006 iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
3008 P <--> The integral from 0 to F of the non-central f-density.
3009 Input range: [0,1-1E-16).
3012 Q is not used by this subroutine and is only included
3013 for similarity with other cdf* routines.
3015 F <--> Upper limit of integration of the non-central f-density.
3016 Input range: [0, +infinity).
3017 Search range: [0,1E100]
3019 DFN < --> Degrees of freedom of the numerator sum of squares.
3020 Input range: (0, +infinity).
3021 Search range: [ 1E-100, 1E100]
3023 DFD < --> Degrees of freedom of the denominator sum of squares.
3024 Must be in range: (0, +infinity).
3025 Input range: (0, +infinity).
3026 Search range: [ 1E-100, 1E100]
3028 PNONC <-> The non-centrality parameter
3029 Input range: [0,infinity)
3030 Search range: [0,1E4]
3032 STATUS <-- 0 if calculation completed correctly
3033 -I if input parameter number I is out of range
3034 1 if answer appears to be lower than lowest
3036 2 if answer appears to be higher than greatest
3040 BOUND <-- Undefined if STATUS is 0
3042 Bound exceeded by parameter number I if STATUS
3045 Lower search bound if STATUS is 1.
3047 Upper search bound if STATUS is 2.
3053 Formula 26.6.20 of Abramowitz and Stegun, Handbook of
3054 Mathematical Functions (1966) is used to compute the cumulative
3055 distribution function.
3057 Computation of other parameters involve a seach for a value that
3058 produces the desired value of P. The search relies on the
3059 monotinicity of P with the other parameter.
3063 The computation time required for this routine is proportional
3064 to the noncentrality parameter (PNONC). Very large values of
3065 this parameter can consume immense computer resources. This is
3066 why the search range is bounded by 10,000.
3070 The value of the cumulative noncentral F distribution is not
3071 necessarily monotone in either degrees of freedom. There thus
3072 may be two values that provide a given CDF value. This routine
3073 assumes monotonicity and will find an arbitrary one of the two
3076 **********************************************************************/
3080 #define atol 1.0e-50
3081 #define zero 1.0e-100
3082 #define one ( 1.0e0 - 1.0e-16 )
3084 static double K1 = 0.0e0;
3085 static double K3 = 0.5e0;
3086 static double K4 = 5.0e0;
3087 static double fx,cum,ccum;
3088 static unsigned long qhi,qleft;
3089 static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17;
3092 .. Executable Statements ..
3097 if(!(*which < 1 || *which > 5)) goto S30;
3098 if(!(*which < 1)) goto S10;
3107 if(*which == 1) goto S70;
3111 if(!(*p < 0.0e0 || *p > one)) goto S60;
3112 if(!(*p < 0.0e0)) goto S40;
3122 if(*which == 2) goto S90;
3126 if(!(*f < 0.0e0)) goto S80;
3132 if(*which == 3) goto S110;
3136 if(!(*dfn <= 0.0e0)) goto S100;
3142 if(*which == 4) goto S130;
3146 if(!(*dfd <= 0.0e0)) goto S120;
3152 if(*which == 5) goto S150;
3156 if(!(*phonc < 0.0e0)) goto S140;
3169 cumfnc(f,dfn,dfd,phonc,p,q);
3172 else if(2 == *which) {
3180 dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
3182 dinvr(status,f,&fx,&qleft,&qhi);
3184 if(!(*status == 1)) goto S170;
3185 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3187 dinvr(status,f,&fx,&qleft,&qhi);
3190 if(!(*status == -1)) goto S200;
3191 if(!qleft) goto S180;
3202 else if(3 == *which) {
3211 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
3213 dinvr(status,dfn,&fx,&qleft,&qhi);
3215 if(!(*status == 1)) goto S220;
3216 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3218 dinvr(status,dfn,&fx,&qleft,&qhi);
3221 if(!(*status == -1)) goto S250;
3222 if(!qleft) goto S230;
3233 else if(4 == *which) {
3242 dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
3244 dinvr(status,dfd,&fx,&qleft,&qhi);
3246 if(!(*status == 1)) goto S270;
3247 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3249 dinvr(status,dfd,&fx,&qleft,&qhi);
3252 if(!(*status == -1)) goto S300;
3253 if(!qleft) goto S280;
3264 else if(5 == *which) {
3272 dstinv(&K1,&T15,&K3,&K3,&K4,&T16,&T17);
3274 dinvr(status,phonc,&fx,&qleft,&qhi);
3276 if(!(*status == 1)) goto S320;
3277 cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
3279 dinvr(status,phonc,&fx,&qleft,&qhi);
3282 if(!(*status == -1)) goto S350;
3283 if(!qleft) goto S330;
3302 void cdfgam(int *which,double *p,double *q,double *x,double *shape,
3303 double *scale,int *status,double *bound)
3304 /**********************************************************************
3306 void cdfgam(int *which,double *p,double *q,double *x,double *shape,
3307 double *scale,int *status,double *bound)
3309 Cumulative Distribution Function
3316 Calculates any one parameter of the gamma
3317 distribution given values for the others.
3323 WHICH --> Integer indicating which of the next four argument
3324 values is to be calculated from the others.
3326 iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE
3327 iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE
3328 iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE
3329 iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE
3331 P <--> The integral from 0 to X of the gamma density.
3335 Input range: (0, 1].
3338 X <--> The upper limit of integration of the gamma density.
3339 Input range: [0, +infinity).
3340 Search range: [0,1E100]
3342 SHAPE <--> The shape parameter of the gamma density.
3343 Input range: (0, +infinity).
3344 Search range: [1E-100,1E100]
3346 SCALE <--> The scale parameter of the gamma density.
3347 Input range: (0, +infinity).
3348 Search range: (1E-100,1E100]
3350 STATUS <-- 0 if calculation completed correctly
3351 -I if input parameter number I is out of range
3352 1 if answer appears to be lower than lowest
3354 2 if answer appears to be higher than greatest
3357 10 if the gamma or inverse gamma routine cannot
3358 compute the answer. Usually happens only for
3359 X and SHAPE very large (gt 1E10 or more)
3361 BOUND <-- Undefined if STATUS is 0
3363 Bound exceeded by parameter number I if STATUS
3366 Lower search bound if STATUS is 1.
3368 Upper search bound if STATUS is 2.
3374 Cumulative distribution function (P) is calculated directly by
3375 the code associated with:
3377 DiDinato, A. R. and Morris, A. H. Computation of the incomplete
3378 gamma function ratios and their inverse. ACM Trans. Math.
3379 Softw. 12 (1986), 377-393.
3381 Computation of other parameters involve a seach for a value that
3382 produces the desired value of P. The search relies on the
3383 monotinicity of P with the other parameter.
3390 The gamma density is proportional to
3391 T**(SHAPE - 1) * EXP(- SCALE * T)
3393 **********************************************************************/
3396 #define atol 1.0e-50
3397 #define zero 1.0e-100
3400 static double K5 = 0.5e0;
3401 static double K6 = 5.0e0;
3402 static double xx,fx,xscale,cum,ccum,pq,porq;
3404 static unsigned long qhi,qleft,qporq;
3405 static double T2,T3,T4,T7,T8,T9;
3408 .. Executable Statements ..
3413 if(!(*which < 1 || *which > 4)) goto S30;
3414 if(!(*which < 1)) goto S10;
3423 if(*which == 1) goto S70;
3427 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
3428 if(!(*p < 0.0e0)) goto S40;
3438 if(*which == 1) goto S110;
3442 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
3443 if(!(*q <= 0.0e0)) goto S80;
3453 if(*which == 2) goto S130;
3457 if(!(*x < 0.0e0)) goto S120;
3463 if(*which == 3) goto S150;
3467 if(!(*shape <= 0.0e0)) goto S140;
3473 if(*which == 4) goto S170;
3477 if(!(*scale <= 0.0e0)) goto S160;
3483 if(*which == 1) goto S210;
3488 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
3489 if(!(pq < 0.0e0)) goto S180;
3499 if(*which == 1) goto S240;
3501 Select the minimum of P or Q
3504 if(!qporq) goto S220;
3520 cumgam(&xscale,shape,p,q);
3521 if(porq > 1.5e0) *status = 10;
3523 else if(2 == *which) {
3528 gaminv(shape,&xx,&T2,p,q,&ierr);
3538 else if(3 == *which) {
3548 dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
3550 dinvr(status,shape,&fx,&qleft,&qhi);
3552 if(!(*status == 1)) goto S290;
3553 cumgam(&xscale,shape,&cum,&ccum);
3554 if(!qporq) goto S260;
3560 if(!((qporq && cum > 1.5e0) || (!qporq && ccum > 1.5e0))) goto S280;
3564 dinvr(status,shape,&fx,&qleft,&qhi);
3567 if(!(*status == -1)) goto S320;
3568 if(!qleft) goto S300;
3579 else if(4 == *which) {
3584 gaminv(shape,&xx,&T9,p,q,&ierr);
3600 void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
3601 double *pr,double *ompr,int *status,double *bound)
3602 /**********************************************************************
3604 void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
3605 double *pr,double *ompr,int *status,double *bound)
3607 Cumulative Distribution Function
3608 Negative BiNomial distribution
3614 Calculates any one parameter of the negative binomial
3615 distribution given values for the others.
3617 The cumulative negative binomial distribution returns the
3618 probability that there will be F or fewer failures before the
3619 XNth success in binomial trials each of which has probability of
3622 The individual term of the negative binomial is the probability of
3623 S failures before XN successes and is
3624 Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S
3630 WHICH --> Integer indicating which of the next four argument
3631 values is to be calculated from the others.
3633 iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
3634 iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
3635 iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
3636 iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
3638 P <--> The cumulation from 0 to S of the negative
3639 binomial distribution.
3643 Input range: (0, 1].
3646 S <--> The upper limit of cumulation of the binomial distribution.
3647 There are F or fewer failures before the XNth success.
3648 Input range: [0, +infinity).
3649 Search range: [0, 1E100]
3651 XN <--> The number of successes.
3652 Input range: [0, +infinity).
3653 Search range: [0, 1E100]
3655 PR <--> The probability of success in each binomial trial.
3657 Search range: [0,1].
3664 STATUS <-- 0 if calculation completed correctly
3665 -I if input parameter number I is out of range
3666 1 if answer appears to be lower than lowest
3668 2 if answer appears to be higher than greatest
3671 4 if PR + OMPR .ne. 1
3673 BOUND <-- Undefined if STATUS is 0
3675 Bound exceeded by parameter number I if STATUS
3678 Lower search bound if STATUS is 1.
3680 Upper search bound if STATUS is 2.
3686 Formula 26.5.26 of Abramowitz and Stegun, Handbook of
3687 Mathematical Functions (1966) is used to reduce calculation of
3688 the cumulative distribution function to that of an incomplete
3691 Computation of other parameters involve a seach for a value that
3692 produces the desired value of P. The search relies on the
3693 monotinicity of P with the other parameter.
3695 **********************************************************************/
3698 #define atol 1.0e-50
3702 static double K2 = 0.0e0;
3703 static double K4 = 0.5e0;
3704 static double K5 = 5.0e0;
3705 static double K11 = 1.0e0;
3706 static double fx,xhi,xlo,pq,prompr,cum,ccum;
3707 static unsigned long qhi,qleft,qporq;
3708 static double T3,T6,T7,T8,T9,T10,T12,T13;
3711 .. Executable Statements ..
3716 if(!(*which < 1 || *which > 4)) goto S30;
3717 if(!(*which < 1)) goto S10;
3726 if(*which == 1) goto S70;
3730 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
3731 if(!(*p < 0.0e0)) goto S40;
3741 if(*which == 1) goto S110;
3745 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
3746 if(!(*q <= 0.0e0)) goto S80;
3756 if(*which == 2) goto S130;
3760 if(!(*s < 0.0e0)) goto S120;
3766 if(*which == 3) goto S150;
3770 if(!(*xn < 0.0e0)) goto S140;
3776 if(*which == 4) goto S190;
3780 if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
3781 if(!(*pr < 0.0e0)) goto S160;
3791 if(*which == 4) goto S230;
3795 if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
3796 if(!(*ompr < 0.0e0)) goto S200;
3806 if(*which == 1) goto S270;
3811 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
3812 if(!(pq < 0.0e0)) goto S240;
3822 if(*which == 4) goto S310;
3827 if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
3828 if(!(prompr < 0.0e0)) goto S280;
3838 if(!(*which == 1)) qporq = *p <= *q;
3840 Select the minimum of P or Q
3847 cumnbn(s,xn,pr,ompr,p,q);
3850 else if(2 == *which) {
3858 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
3860 dinvr(status,s,&fx,&qleft,&qhi);
3862 if(!(*status == 1)) goto S350;
3863 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3864 if(!qporq) goto S330;
3870 dinvr(status,s,&fx,&qleft,&qhi);
3873 if(!(*status == -1)) goto S380;
3874 if(!qleft) goto S360;
3885 else if(3 == *which) {
3893 dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
3895 dinvr(status,xn,&fx,&qleft,&qhi);
3897 if(!(*status == 1)) goto S420;
3898 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3899 if(!qporq) goto S400;
3905 dinvr(status,xn,&fx,&qleft,&qhi);
3908 if(!(*status == -1)) goto S450;
3909 if(!qleft) goto S430;
3920 else if(4 == *which) {
3922 Calculating PR and OMPR
3926 dstzr(&K2,&K11,&T12,&T13);
3927 if(!qporq) goto S480;
3929 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
3932 if(!(*status == 1)) goto S470;
3933 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3935 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
3942 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
3945 if(!(*status == 1)) goto S500;
3946 cumnbn(s,xn,pr,ompr,&cum,&ccum);
3948 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
3953 if(!(*status == -1)) goto S540;
3954 if(!qleft) goto S520;
3971 void cdfnor(int *which,double *p,double *q,double *x,double *mean,
3972 double *sd,int *status,double *bound)
3973 /**********************************************************************
3975 void cdfnor(int *which,double *p,double *q,double *x,double *mean,
3976 double *sd,int *status,double *bound)
3978 Cumulative Distribution Function
3985 Calculates any one parameter of the normal
3986 distribution given values for the others.
3992 WHICH --> Integer indicating which of the next parameter
3993 values is to be calculated using values of the others.
3995 iwhich = 1 : Calculate P and Q from X,MEAN and SD
3996 iwhich = 2 : Calculate X from P,Q,MEAN and SD
3997 iwhich = 3 : Calculate MEAN from P,Q,X and SD
3998 iwhich = 4 : Calculate SD from P,Q,X and MEAN
4000 P <--> The integral from -infinity to X of the normal density.
4004 Input range: (0, 1].
4007 X < --> Upper limit of integration of the normal-density.
4008 Input range: ( -infinity, +infinity)
4010 MEAN <--> The mean of the normal density.
4011 Input range: (-infinity, +infinity)
4013 SD <--> Standard Deviation of the normal density.
4014 Input range: (0, +infinity).
4016 STATUS <-- 0 if calculation completed correctly
4017 -I if input parameter number I is out of range
4018 1 if answer appears to be lower than lowest
4020 2 if answer appears to be higher than greatest
4024 BOUND <-- Undefined if STATUS is 0
4026 Bound exceeded by parameter number I if STATUS
4029 Lower search bound if STATUS is 1.
4031 Upper search bound if STATUS is 2.
4039 A slightly modified version of ANORM from
4041 Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
4042 Package of Special Function Routines and Test Drivers"
4043 acm Transactions on Mathematical Software. 19, 22-32.
4045 is used to calulate the cumulative standard normal distribution.
4047 The rational functions from pages 90-95 of Kennedy and Gentle,
4048 Statistical Computing, Marcel Dekker, NY, 1980 are used as
4049 starting values to Newton's Iterations which compute the inverse
4050 standard normal. Therefore no searches are necessary for any
4053 For X < -15, the asymptotic expansion for the normal is used as
4054 the starting value in finding the inverse standard normal.
4055 This is formula 26.2.12 of Abramowitz and Stegun.
4061 The normal density is proportional to
4062 exp( - 0.5 * (( X - MEAN)/SD)**2)
4064 **********************************************************************/
4070 .. Executable Statements ..
4076 if(!(*which < 1 || *which > 4)) goto S30;
4077 if(!(*which < 1)) goto S10;
4086 if(*which == 1) goto S70;
4090 if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
4091 if(!(*p <= 0.0e0)) goto S40;
4101 if(*which == 1) goto S110;
4105 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4106 if(!(*q <= 0.0e0)) goto S80;
4116 if(*which == 1) goto S150;
4121 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S140;
4122 if(!(pq < 0.0e0)) goto S120;
4132 if(*which == 4) goto S170;
4136 if(!(*sd <= 0.0e0)) goto S160;
4149 z = (*x-*mean)/ *sd;
4152 else if(2 == *which) {
4159 else if(3 == *which) {
4166 else if(4 == *which) {
4175 void cdfpoi(int *which,double *p,double *q,double *s,double *xlam,
4176 int *status,double *bound)
4177 /**********************************************************************
4179 void cdfpoi(int *which,double *p,double *q,double *s,double *xlam,
4180 int *status,double *bound)
4182 Cumulative Distribution Function
4183 POIsson distribution
4189 Calculates any one parameter of the Poisson
4190 distribution given values for the others.
4196 WHICH --> Integer indicating which argument
4197 value is to be calculated from the others.
4199 iwhich = 1 : Calculate P and Q from S and XLAM
4200 iwhich = 2 : Calculate A from P,Q and XLAM
4201 iwhich = 3 : Calculate XLAM from P,Q and S
4203 P <--> The cumulation from 0 to S of the poisson density.
4207 Input range: (0, 1].
4210 S <--> Upper limit of cumulation of the Poisson.
4211 Input range: [0, +infinity).
4212 Search range: [0,1E100]
4214 XLAM <--> Mean of the Poisson distribution.
4215 Input range: [0, +infinity).
4216 Search range: [0,1E100]
4218 STATUS <-- 0 if calculation completed correctly
4219 -I if input parameter number I is out of range
4220 1 if answer appears to be lower than lowest
4222 2 if answer appears to be higher than greatest
4226 BOUND <-- Undefined if STATUS is 0
4228 Bound exceeded by parameter number I if STATUS
4231 Lower search bound if STATUS is 1.
4233 Upper search bound if STATUS is 2.
4239 Formula 26.4.21 of Abramowitz and Stegun, Handbook of
4240 Mathematical Functions (1966) is used to reduce the computation
4241 of the cumulative distribution function to that of computing a
4242 chi-square, hence an incomplete gamma function.
4244 Cumulative distribution function (P) is calculated directly.
4245 Computation of other parameters involve a seach for a value that
4246 produces the desired value of P. The search relies on the
4247 monotinicity of P with the other parameter.
4249 **********************************************************************/
4252 #define atol 1.0e-50
4255 static double K2 = 0.0e0;
4256 static double K4 = 0.5e0;
4257 static double K5 = 5.0e0;
4258 static double fx,cum,ccum,pq;
4259 static unsigned long qhi,qleft,qporq;
4260 static double T3,T6,T7,T8,T9,T10;
4263 .. Executable Statements ..
4268 if(!(*which < 1 || *which > 3)) goto S30;
4269 if(!(*which < 1)) goto S10;
4278 if(*which == 1) goto S70;
4282 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
4283 if(!(*p < 0.0e0)) goto S40;
4293 if(*which == 1) goto S110;
4297 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4298 if(!(*q <= 0.0e0)) goto S80;
4308 if(*which == 2) goto S130;
4312 if(!(*s < 0.0e0)) goto S120;
4318 if(*which == 3) goto S150;
4322 if(!(*xlam < 0.0e0)) goto S140;
4328 if(*which == 1) goto S190;
4333 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
4334 if(!(pq < 0.0e0)) goto S160;
4344 if(!(*which == 1)) qporq = *p <= *q;
4346 Select the minimum of P or Q
4356 else if(2 == *which) {
4364 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
4366 dinvr(status,s,&fx,&qleft,&qhi);
4368 if(!(*status == 1)) goto S230;
4369 cumpoi(s,xlam,&cum,&ccum);
4370 if(!qporq) goto S210;
4376 dinvr(status,s,&fx,&qleft,&qhi);
4379 if(!(*status == -1)) goto S260;
4380 if(!qleft) goto S240;
4391 else if(3 == *which) {
4399 dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
4401 dinvr(status,xlam,&fx,&qleft,&qhi);
4403 if(!(*status == 1)) goto S300;
4404 cumpoi(s,xlam,&cum,&ccum);
4405 if(!qporq) goto S280;
4411 dinvr(status,xlam,&fx,&qleft,&qhi);
4414 if(!(*status == -1)) goto S330;
4415 if(!qleft) goto S310;
4431 void cdft(int *which,double *p,double *q,double *t,double *df,
4432 int *status,double *bound)
4433 /**********************************************************************
4435 void cdft(int *which,double *p,double *q,double *t,double *df,
4436 int *status,double *bound)
4438 Cumulative Distribution Function
4445 Calculates any one parameter of the t distribution given
4446 values for the others.
4452 WHICH --> Integer indicating which argument
4453 values is to be calculated from the others.
4455 iwhich = 1 : Calculate P and Q from T and DF
4456 iwhich = 2 : Calculate T from P,Q and DF
4457 iwhich = 3 : Calculate DF from P,Q and T
4459 P <--> The integral from -infinity to t of the t-density.
4463 Input range: (0, 1].
4466 T <--> Upper limit of integration of the t-density.
4467 Input range: ( -infinity, +infinity).
4468 Search range: [ -1E100, 1E100 ]
4470 DF <--> Degrees of freedom of the t-distribution.
4471 Input range: (0 , +infinity).
4472 Search range: [1e-100, 1E10]
4474 STATUS <-- 0 if calculation completed correctly
4475 -I if input parameter number I is out of range
4476 1 if answer appears to be lower than lowest
4478 2 if answer appears to be higher than greatest
4482 BOUND <-- Undefined if STATUS is 0
4484 Bound exceeded by parameter number I if STATUS
4487 Lower search bound if STATUS is 1.
4489 Upper search bound if STATUS is 2.
4495 Formula 26.5.27 of Abramowitz and Stegun, Handbook of
4496 Mathematical Functions (1966) is used to reduce the computation
4497 of the cumulative distribution function to that of an incomplete
4500 Computation of other parameters involve a seach for a value that
4501 produces the desired value of P. The search relies on the
4502 monotinicity of P with the other parameter.
4504 **********************************************************************/
4507 #define atol 1.0e-50
4508 #define zero 1.0e-100
4510 #define rtinf 1.0e100
4511 #define maxdf 1.0e10
4513 static double K4 = 0.5e0;
4514 static double K5 = 5.0e0;
4515 static double fx,cum,ccum,pq;
4516 static unsigned long qhi,qleft,qporq;
4517 static double T2,T3,T6,T7,T8,T9,T10,T11;
4520 .. Executable Statements ..
4525 if(!(*which < 1 || *which > 3)) goto S30;
4526 if(!(*which < 1)) goto S10;
4535 if(*which == 1) goto S70;
4539 if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
4540 if(!(*p <= 0.0e0)) goto S40;
4550 if(*which == 1) goto S110;
4554 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
4555 if(!(*q <= 0.0e0)) goto S80;
4565 if(*which == 3) goto S130;
4569 if(!(*df <= 0.0e0)) goto S120;
4575 if(*which == 1) goto S170;
4580 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S160;
4581 if(!(pq < 0.0e0)) goto S140;
4591 if(!(*which == 1)) qporq = *p <= *q;
4593 Select the minimum of P or Q
4603 else if(2 == *which) {
4606 .. Get initial approximation for T
4613 dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
4615 dinvr(status,t,&fx,&qleft,&qhi);
4617 if(!(*status == 1)) goto S210;
4618 cumt(t,df,&cum,&ccum);
4619 if(!qporq) goto S190;
4625 dinvr(status,t,&fx,&qleft,&qhi);
4628 if(!(*status == -1)) goto S240;
4629 if(!qleft) goto S220;
4640 else if(3 == *which) {
4649 dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
4651 dinvr(status,df,&fx,&qleft,&qhi);
4653 if(!(*status == 1)) goto S280;
4654 cumt(t,df,&cum,&ccum);
4655 if(!qporq) goto S260;
4661 dinvr(status,df,&fx,&qleft,&qhi);
4664 if(!(*status == -1)) goto S310;
4665 if(!qleft) goto S290;
4684 void cdftnc(int *which,double *p,double *q,double *t,double *df,
4685 double *pnonc,int *status,double *bound)
4686 /**********************************************************************
4688 void cdftnc(int *which,double *p,double *q,double *t,double *df,
4689 double *pnonc,int *status,double *bound)
4691 Cumulative Distribution Function
4692 Non-Central T distribution
4696 Calculates any one parameter of the noncentral t distribution give
4697 values for the others.
4701 WHICH --> Integer indicating which argument
4702 values is to be calculated from the others.
4704 iwhich = 1 : Calculate P and Q from T,DF,PNONC
4705 iwhich = 2 : Calculate T from P,Q,DF,PNONC
4706 iwhich = 3 : Calculate DF from P,Q,T
4707 iwhich = 4 : Calculate PNONC from P,Q,DF,T
4709 P <--> The integral from -infinity to t of the noncentral t-den
4713 Input range: (0, 1].
4716 T <--> Upper limit of integration of the noncentral t-density.
4717 Input range: ( -infinity, +infinity).
4718 Search range: [ -1E100, 1E100 ]
4720 DF <--> Degrees of freedom of the noncentral t-distribution.
4721 Input range: (0 , +infinity).
4722 Search range: [1e-100, 1E10]
4724 PNONC <--> Noncentrality parameter of the noncentral t-distributio
4725 Input range: [-infinity , +infinity).
4726 Search range: [-1e4, 1E4]
4728 STATUS <-- 0 if calculation completed correctly
4729 -I if input parameter number I is out of range
4730 1 if answer appears to be lower than lowest
4732 2 if answer appears to be higher than greatest
4736 BOUND <-- Undefined if STATUS is 0
4738 Bound exceeded by parameter number I if STATUS
4741 Lower search bound if STATUS is 1.
4743 Upper search bound if STATUS is 2.
4747 Upper tail of the cumulative noncentral t is calculated usin
4748 formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuou
4749 Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
4751 Computation of other parameters involve a seach for a value that
4752 produces the desired value of P. The search relies on the
4753 monotinicity of P with the other parameter.
4755 **********************************************************************/
4759 #define atol 1.0e-50
4760 #define zero 1.0e-100
4761 #define one ( 1.0e0 - 1.0e-16 )
4763 static double K3 = 0.5e0;
4764 static double K4 = 5.0e0;
4765 static double ccum,cum,fx;
4766 static unsigned long qhi,qleft;
4767 static double T1,T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14;
4770 .. Executable Statements ..
4772 if(!(*which < 1 || *which > 4)) goto S30;
4773 if(!(*which < 1)) goto S10;
4782 if(*which == 1) goto S70;
4783 if(!(*p < 0.0e0 || *p > one)) goto S60;
4784 if(!(*p < 0.0e0)) goto S40;
4794 if(*which == 3) goto S90;
4795 if(!(*df <= 0.0e0)) goto S80;
4801 if(*which == 4) goto S100;
4804 cumtnc(t,df,pnonc,p,q);
4807 else if(2 == *which) {
4813 dstinv(&T1,&T2,&K3,&K3,&K4,&T5,&T6);
4815 dinvr(status,t,&fx,&qleft,&qhi);
4817 if(!(*status == 1)) goto S120;
4818 cumtnc(t,df,pnonc,&cum,&ccum);
4820 dinvr(status,t,&fx,&qleft,&qhi);
4823 if(!(*status == -1)) goto S150;
4824 if(!qleft) goto S130;
4835 else if(3 == *which) {
4841 dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
4843 dinvr(status,df,&fx,&qleft,&qhi);
4845 if(!(*status == 1)) goto S170;
4846 cumtnc(t,df,pnonc,&cum,&ccum);
4848 dinvr(status,df,&fx,&qleft,&qhi);
4851 if(!(*status == -1)) goto S200;
4852 if(!qleft) goto S180;
4863 else if(4 == *which) {
4869 dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
4871 dinvr(status,pnonc,&fx,&qleft,&qhi);
4873 if(!(*status == 1)) goto S220;
4874 cumtnc(t,df,pnonc,&cum,&ccum);
4876 dinvr(status,pnonc,&fx,&qleft,&qhi);
4879 if(!(*status == -1)) goto S250;
4880 if(!qleft) goto S230;
4899 void cumbet(double *x,double *y,double *a,double *b,double *cum,
4902 **********************************************************************
4904 void cumbet(double *x,double *y,double *a,double *b,double *cum,
4907 Double precision cUMulative incomplete BETa distribution
4913 Calculates the cdf to X of the incomplete beta distribution
4914 with parameters a and b. This is the integral from 0 to x
4915 of (1/B(a,b))*f(t)) where f(t) = t**(a-1) * (1-t)**(b-1)
4921 X --> Upper limit of integration.
4922 X is DOUBLE PRECISION
4925 Y is DOUBLE PRECISION
4927 A --> First parameter of the beta distribution.
4928 A is DOUBLE PRECISION
4930 B --> Second parameter of the beta distribution.
4931 B is DOUBLE PRECISION
4933 CUM <-- Cumulative incomplete beta distribution.
4934 CUM is DOUBLE PRECISION
4936 CCUM <-- Compliment of Cumulative incomplete beta distribution.
4937 CCUM is DOUBLE PRECISION
4943 Calls the routine BRATIO.
4947 Didonato, Armido R. and Morris, Alfred H. Jr. (1992) Algorithim
4948 708 Significant Digit Computation of the Incomplete Beta Function
4949 Ratios. ACM ToMS, Vol.18, No. 3, Sept. 1992, 360-373.
4951 **********************************************************************
4957 .. Executable Statements ..
4959 if(!(*x <= 0.0e0)) goto S10;
4964 if(!(*y <= 0.0e0)) goto S20;
4969 bratio(a,b,x,y,cum,ccum,&ierr);
4975 void cumbin(double *s,double *xn,double *pr,double *ompr,
4976 double *cum,double *ccum)
4978 **********************************************************************
4980 void cumbin(double *s,double *xn,double *pr,double *ompr,
4981 double *cum,double *ccum)
4983 CUmulative BINomial distribution
4989 Returns the probability of 0 to S successes in XN binomial
4990 trials, each of which has a probability of success, PBIN.
4996 S --> The upper limit of cumulation of the binomial distribution.
4997 S is DOUBLE PRECISION
4999 XN --> The number of binomial trials.
5000 XN is DOUBLE PRECISIO
5002 PBIN --> The probability of success in each binomial trial.
5003 PBIN is DOUBLE PRECIS
5006 OMPR is DOUBLE PRECIS
5008 CUM <-- Cumulative binomial distribution.
5009 CUM is DOUBLE PRECISI
5011 CCUM <-- Compliment of Cumulative binomial distribution.
5012 CCUM is DOUBLE PRECIS
5018 Formula 26.5.24 of Abramowitz and Stegun, Handbook of
5019 Mathematical Functions (1966) is used to reduce the binomial
5020 distribution to the cumulative beta distribution.
5022 **********************************************************************
5025 static double T1,T2;
5028 .. Executable Statements ..
5030 if(!(*s < *xn)) goto S10;
5033 cumbet(pr,ompr,&T1,&T2,ccum,cum);
5041 void cumchi(double *x,double *df,double *cum,double *ccum)
5043 **********************************************************************
5045 void cumchi(double *x,double *df,double *cum,double *ccum)
5046 CUMulative of the CHi-square distribution
5052 Calculates the cumulative chi-square distribution.
5058 X --> Upper limit of integration of the
5059 chi-square distribution.
5060 X is DOUBLE PRECISION
5062 DF --> Degrees of freedom of the
5063 chi-square distribution.
5064 DF is DOUBLE PRECISION
5066 CUM <-- Cumulative chi-square distribution.
5067 CUM is DOUBLE PRECISIO
5069 CCUM <-- Compliment of Cumulative chi-square distribution.
5070 CCUM is DOUBLE PRECISI
5076 Calls incomplete gamma function (CUMGAM)
5078 **********************************************************************
5084 .. Executable Statements ..
5088 cumgam(&xx,&a,cum,ccum);
5091 void cumchn(double *x,double *df,double *pnonc,double *cum,
5093 /**********************************************************************
5095 void cumchn(double *x,double *df,double *pnonc,double *cum,
5098 CUMulative of the Non-central CHi-square distribution
5102 Calculates the cumulative non-central chi-square
5103 distribution, i.e., the probability that a random variable
5104 which follows the non-central chi-square distribution, with
5105 non-centrality parameter PNONC and continuous degrees of
5106 freedom DF, is less than or equal to X.
5110 X --> Upper limit of integration of the non-central
5111 chi-square distribution.
5113 DF --> Degrees of freedom of the non-central
5114 chi-square distribution.
5116 PNONC --> Non-centrality parameter of the non-central
5117 chi-square distribution.
5119 CUM <-- Cumulative non-central chi-square distribution.
5121 CCUM <-- Compliment of Cumulative non-central chi-square distribut
5126 Uses formula 26.4.25 of Abramowitz and Stegun, Handbook of
5127 Mathematical Functions, US NBS (1966) to calculate the
5128 non-central chi-square.
5132 EPS --- Convergence criterion. The sum stops when a
5133 term is less than EPS*SUM.
5135 CCUM <-- Compliment of Cumulative non-central
5136 chi-square distribution.
5138 **********************************************************************/
5140 #define dg(i) (*df + 2.0e0 * (double)(i))
5141 #define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps * sum)
5142 static double eps = 1.0e-5;
5143 static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
5144 sumadj,term,wt,xnonc;
5146 static double T1,T2,T3;
5149 .. Executable Statements ..
5151 if(!(*x <= 0.0e0)) goto S10;
5156 if(!(*pnonc <= 1.0e-10 )) goto S20;
5158 When non-centrality parameter is (essentially) zero,
5159 use cumulative chi-square distribution
5161 cumchi(x,df,cum,ccum);
5164 xnonc = *pnonc / 2.0e0;
5166 ***********************************************************************
5167 The following code calcualtes the weight, chi-square, and
5168 adjustment term for the central term in the infinite series.
5169 The central term is the one in which the poisson weight is
5170 greatest. The adjustment term is the amount that must
5171 be subtracted from the chi-square to move up two degrees
5173 ***********************************************************************
5175 icent = fifidint(xnonc);
5176 if(icent == 0) icent = 1;
5179 Calculate central weight term
5181 T1 = (double)(icent + 1);
5182 lfact = alngam(&T1);
5183 lcntwt = -xnonc + (double)icent * log(xnonc) - lfact;
5184 centwt = exp(lcntwt);
5186 Calculate central chi-square
5189 cumchi(x,&T2,&pcent,ccum);
5191 Calculate central adjustment term
5193 dfd2 = dg(icent) / 2.0e0;
5195 lfact = alngam(&T3);
5196 lcntaj = dfd2 * log(chid2) - chid2 - lfact;
5197 centaj = exp(lcntaj);
5198 sum = centwt * pcent;
5200 ***********************************************************************
5201 Sum backwards from the central term towards zero.
5202 Quit whenever either
5203 (1) the zero term is reached, or
5204 (2) the term gets small relative to the sum
5205 ***********************************************************************
5213 if(qsmall(term) || i == 0) goto S50;
5215 dfd2 = dg(i) / 2.0e0;
5217 Adjust chi-square for two fewer degrees of freedom.
5218 The adjusted value ends up in PTERM.
5220 adj = adj * dfd2 / chid2;
5222 pterm = pcent + sumadj;
5224 Adjust poisson weight for J decreased by one
5226 wt *= ((double)i / xnonc);
5233 ***********************************************************************
5234 Now sum forward from the central term towards infinity.
5236 (1) the term gets small relative to the sum, or
5237 ***********************************************************************
5239 sumadj = adj = centaj;
5244 if(qsmall(term)) goto S80;
5247 Update weights for next higher J
5249 wt *= (xnonc / (double)(i + 1));
5251 Calculate PTERM and add term to sum
5253 pterm = pcent - sumadj;
5257 Update adjustment term for DF for next iteration
5260 dfd2 = dg(i) / 2.0e0;
5261 adj = adj * chid2 / dfd2;
5266 *ccum = 0.5e0 + (0.5e0 - *cum);
5271 void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
5273 **********************************************************************
5275 void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
5276 CUMulative F distribution
5282 Computes the integral from 0 to F of the f-density with DFN
5283 and DFD degrees of freedom.
5289 F --> Upper limit of integration of the f-density.
5290 F is DOUBLE PRECISION
5292 DFN --> Degrees of freedom of the numerator sum of squares.
5293 DFN is DOUBLE PRECISI
5295 DFD --> Degrees of freedom of the denominator sum of squares.
5296 DFD is DOUBLE PRECISI
5298 CUM <-- Cumulative f distribution.
5299 CUM is DOUBLE PRECISI
5301 CCUM <-- Compliment of Cumulative f distribution.
5302 CCUM is DOUBLE PRECIS
5308 Formula 26.5.28 of Abramowitz and Stegun is used to reduce
5309 the cumulative F to a cumulative beta distribution.
5315 If F is less than or equal to 0, 0 is returned.
5317 **********************************************************************
5322 static double dsum,prod,xx,yy;
5324 static double T1,T2;
5327 .. Executable Statements ..
5329 if(!(*f <= 0.0e0)) goto S10;
5336 XX is such that the incomplete beta with parameters
5337 DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM
5339 Calculate the smaller of XX and YY accurately
5350 bratio(&T1,&T2,&xx,&yy,ccum,cum,&ierr);
5355 void cumfnc(double *f,double *dfn,double *dfd,double *pnonc,
5356 double *cum,double *ccum)
5358 **********************************************************************
5360 F -NON- -C-ENTRAL F DISTRIBUTION
5367 COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
5368 DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC
5374 X --> UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION
5376 DFN --> DEGREES OF FREEDOM OF NUMERATOR
5378 DFD --> DEGREES OF FREEDOM OF DENOMINATOR
5380 PNONC --> NONCENTRALITY PARAMETER.
5382 CUM <-- CUMULATIVE NONCENTRAL F DISTRIBUTION
5384 CCUM <-- COMPLIMENT OF CUMMULATIVE
5390 USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
5391 SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
5392 (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
5393 THE CONVERGENCE CRITERION IS MET.
5395 FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
5402 HANDBOOD OF MATHEMATICAL FUNCTIONS
5403 EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
5404 NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
5406 P 947, EQUATIONS 26.6.17, 26.6.18
5412 THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
5413 TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20). EPS IS
5414 SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED.
5416 **********************************************************************
5419 #define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
5422 static double eps = 1.0e-4;
5423 static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
5425 static int i,icent,ierr;
5426 static double T1,T2,T3,T4,T5,T6;
5429 .. Executable Statements ..
5431 if(!(*f <= 0.0e0)) goto S10;
5436 if(!(*pnonc < 1.0e-10)) goto S20;
5438 Handle case in which the non-centrality parameter is
5441 cumf(f,dfn,dfd,cum,ccum);
5444 xnonc = *pnonc/2.0e0;
5446 Calculate the central term of the poisson weighting factor.
5448 icent = (long)(xnonc);
5449 if(icent == 0) icent = 1;
5451 Compute central weight term
5453 T1 = (double)(icent+1);
5454 centwt = exp(-xnonc+(double)icent*log(xnonc)-alngam(&T1));
5456 Compute central incomplete beta term
5457 Assure that minimum of arg to beta and 1 - arg is computed
5468 T2 = *dfn*half+(double)icent;
5470 bratio(&T2,&T3,&xx,&yy,&betdn,&dummy,&ierr);
5471 adn = *dfn/2.0e0+(double)icent;
5477 Now sum terms backward from icent until convergence or all done
5483 dnterm = exp(alngam(&T4)-alngam(&T5)-alngam(&b)+adn*log(xx)+b*log(yy));
5485 if(qsmall(xmult*betdn) || i <= 0) goto S40;
5486 xmult *= ((double)i/xnonc);
5489 dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
5491 sum += (xmult*betdn);
5496 Now sum forwards until convergence
5499 if(aup-1.0+b == 0) upterm = exp(-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+
5503 upterm = exp(alngam(&T6)-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+b*
5508 if(qsmall(xmult*betup)) goto S70;
5510 xmult *= (xnonc/(double)i);
5513 upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
5515 sum += (xmult*betup);
5519 *ccum = 0.5e0+(0.5e0-*cum);
5525 void cumgam(double *x,double *a,double *cum,double *ccum)
5527 **********************************************************************
5529 void cumgam(double *x,double *a,double *cum,double *ccum)
5530 Double precision cUMulative incomplete GAMma distribution
5536 Computes the cumulative of the incomplete gamma
5537 distribution, i.e., the integral from 0 to X of
5538 (1/GAM(A))*EXP(-T)*T**(A-1) DT
5539 where GAM(A) is the complete gamma function of A, i.e.,
5540 GAM(A) = integral from 0 to infinity of
5547 X --> The upper limit of integration of the incomplete gamma.
5548 X is DOUBLE PRECISION
5550 A --> The shape parameter of the incomplete gamma.
5551 A is DOUBLE PRECISION
5553 CUM <-- Cumulative incomplete gamma distribution.
5554 CUM is DOUBLE PRECISION
5556 CCUM <-- Compliment of Cumulative incomplete gamma distribution.
5557 CCUM is DOUBLE PRECISIO
5563 Calls the routine GRATIO.
5565 **********************************************************************
5571 .. Executable Statements ..
5573 if(!(*x <= 0.0e0)) goto S10;
5578 gratio(a,x,cum,ccum,&K1);
5584 void cumnbn(double *s,double *xn,double *pr,double *ompr,
5585 double *cum,double *ccum)
5587 **********************************************************************
5589 void cumnbn(double *s,double *xn,double *pr,double *ompr,
5590 double *cum,double *ccum)
5592 CUmulative Negative BINomial distribution
5598 Returns the probability that it there will be S or fewer failures
5599 before there are XN successes, with each binomial trial having
5600 a probability of success PR.
5602 Prob(# failures = S | XN successes, PR) =
5604 ( ) * PR^XN * (1-PR)^S
5611 S --> The number of failures
5612 S is DOUBLE PRECISION
5614 XN --> The number of successes
5615 XN is DOUBLE PRECISIO
5617 PR --> The probability of success in each binomial trial.
5618 PR is DOUBLE PRECISIO
5621 OMPR is DOUBLE PRECIS
5623 CUM <-- Cumulative negative binomial distribution.
5624 CUM is DOUBLE PRECISI
5626 CCUM <-- Compliment of Cumulative negative binomial distribution.
5627 CCUM is DOUBLE PRECIS
5633 Formula 26.5.26 of Abramowitz and Stegun, Handbook of
5634 Mathematical Functions (1966) is used to reduce the negative
5635 binomial distribution to the cumulative beta distribution.
5637 **********************************************************************
5643 .. Executable Statements ..
5646 cumbet(pr,ompr,xn,&T1,cum,ccum);
5649 void cumnor(double *arg,double *result,double *ccum)
5651 **********************************************************************
5653 void cumnor(double *arg,double *result,double *ccum)
5659 Computes the cumulative of the normal distribution, i.e.,
5660 the integral from -infinity to x of
5661 (1/sqrt(2*pi)) exp(-u*u/2) du
5663 X --> Upper limit of integration.
5664 X is DOUBLE PRECISION
5666 RESULT <-- Cumulative normal distribution.
5667 RESULT is DOUBLE PRECISION
5669 CCUM <-- Compliment of Cumulative normal distribution.
5670 CCUM is DOUBLE PRECISION
5672 Renaming of function ANORM from:
5674 Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
5675 Package of Special Function Routines and Test Drivers"
5676 acm Transactions on Mathematical Software. 19, 22-32.
5678 with slight modifications to return ccum and to deal with
5681 **********************************************************************
5683 ------------------------------------------------------------------
5685 This function evaluates the normal distribution function:
5689 P(x) = ----------- | e dt
5693 The main computation evaluates near-minimax approximations
5694 derived from those in "Rational Chebyshev approximations for
5695 the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
5696 This transportable program uses rational functions that
5697 theoretically approximate the normal distribution function to
5698 at least 18 significant decimal digits. The accuracy achieved
5699 depends on the arithmetic system, the compiler, the intrinsic
5700 functions, and proper selection of the machine-dependent
5703 *******************************************************************
5704 *******************************************************************
5706 Explanation of machine-dependent constants.
5708 MIN = smallest machine representable number.
5710 EPS = argument below which anorm(x) may be represented by
5711 0.5 and above which x*x will not underflow.
5712 A conservative value is the largest machine number X
5713 such that 1.0 + X = 1.0 to machine precision.
5714 *******************************************************************
5715 *******************************************************************
5719 The program returns ANORM = 0 for ARG .LE. XLOW.
5722 Intrinsic functions required are:
5728 Mathematics and Computer Science Division
5729 Argonne National Laboratory
5732 Latest modification: March 15, 1992
5734 ------------------------------------------------------------------
5737 static double a[5] = {
5738 2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
5739 1.8154981253343561249e04,6.5682337918207449113e-2
5741 static double b[4] = {
5742 4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
5743 4.5507789335026729956e04
5745 static double c[9] = {
5746 3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
5747 5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
5748 1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
5750 static double d[8] = {
5751 2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
5752 6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
5753 3.8912003286093271411e04,1.9685429676859990727e04
5755 static double half = 0.5e0;
5756 static double p[6] = {
5757 2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
5758 1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
5760 static double one = 1.0e0;
5761 static double q[5] = {
5762 1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
5763 3.78239633202758244e-3,7.29751555083966205e-5
5765 static double sixten = 1.60e0;
5766 static double sqrpi = 3.9894228040143267794e-1;
5767 static double thrsh = 0.66291e0;
5768 static double root32 = 5.656854248e0;
5769 static double zero = 0.0e0;
5773 static double del,eps,temp,x,xden,xnum,y,xsq,min;
5775 ------------------------------------------------------------------
5776 Machine dependent constants
5777 ------------------------------------------------------------------
5779 eps = spmpar(&K1)*0.5e0;
5785 ------------------------------------------------------------------
5786 Evaluate anorm for |X| <= 0.66291
5787 ------------------------------------------------------------------
5790 if(y > eps) xsq = x*x;
5793 for(i=0; i<3; i++) {
5794 xnum = (xnum+a[i])*xsq;
5795 xden = (xden+b[i])*xsq;
5797 *result = x*(xnum+a[3])/(xden+b[3]);
5799 *result = half+temp;
5803 ------------------------------------------------------------------
5804 Evaluate anorm for 0.66291 <= |X| <= sqrt(32)
5805 ------------------------------------------------------------------
5807 else if(y <= root32) {
5810 for(i=0; i<7; i++) {
5811 xnum = (xnum+c[i])*y;
5812 xden = (xden+d[i])*y;
5814 *result = (xnum+c[7])/(xden+d[7]);
5815 xsq = fifdint(y*sixten)/sixten;
5816 del = (y-xsq)*(y+xsq);
5817 *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
5818 *ccum = one-*result;
5826 ------------------------------------------------------------------
5827 Evaluate anorm for |X| > sqrt(32)
5828 ------------------------------------------------------------------
5835 for(i=0; i<4; i++) {
5836 xnum = (xnum+p[i])*xsq;
5837 xden = (xden+q[i])*xsq;
5839 *result = xsq*(xnum+p[4])/(xden+q[4]);
5840 *result = (sqrpi-*result)/y;
5841 xsq = fifdint(x*sixten)/sixten;
5842 del = (x-xsq)*(x+xsq);
5843 *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
5844 *ccum = one-*result;
5851 if(*result < min) *result = 0.0e0;
5853 ------------------------------------------------------------------
5854 Fix up for negative argument, erf, etc.
5855 ------------------------------------------------------------------
5856 ----------Last card of ANORM ----------
5858 if(*ccum < min) *ccum = 0.0e0;
5860 void cumpoi(double *s,double *xlam,double *cum,double *ccum)
5862 **********************************************************************
5864 void cumpoi(double *s,double *xlam,double *cum,double *ccum)
5865 CUMulative POIsson distribution
5871 Returns the probability of S or fewer events in a Poisson
5872 distribution with mean XLAM.
5878 S --> Upper limit of cumulation of the Poisson.
5879 S is DOUBLE PRECISION
5881 XLAM --> Mean of the Poisson distribution.
5882 XLAM is DOUBLE PRECIS
5884 CUM <-- Cumulative poisson distribution.
5885 CUM is DOUBLE PRECISION
5887 CCUM <-- Compliment of Cumulative poisson distribution.
5888 CCUM is DOUBLE PRECIS
5894 Uses formula 26.4.21 of Abramowitz and Stegun, Handbook of
5895 Mathematical Functions to reduce the cumulative Poisson to
5896 the cumulative chi-square distribution.
5898 **********************************************************************
5901 static double chi,df;
5904 .. Executable Statements ..
5906 df = 2.0e0*(*s+1.0e0);
5908 cumchi(&chi,&df,ccum,cum);
5911 void cumt(double *t,double *df,double *cum,double *ccum)
5913 **********************************************************************
5915 void cumt(double *t,double *df,double *cum,double *ccum)
5916 CUMulative T-distribution
5922 Computes the integral from -infinity to T of the t-density.
5928 T --> Upper limit of integration of the t-density.
5929 T is DOUBLE PRECISION
5931 DF --> Degrees of freedom of the t-distribution.
5932 DF is DOUBLE PRECISIO
5934 CUM <-- Cumulative t-distribution.
5935 CCUM is DOUBLE PRECIS
5937 CCUM <-- Compliment of Cumulative t-distribution.
5938 CCUM is DOUBLE PRECIS
5944 Formula 26.5.27 of Abramowitz and Stegun, Handbook of
5945 Mathematical Functions is used to reduce the t-distribution
5946 to an incomplete beta.
5948 **********************************************************************
5951 static double K2 = 0.5e0;
5952 static double xx,a,oma,tt,yy,dfptt,T1;
5955 .. Executable Statements ..
5962 cumbet(&xx,&yy,&T1,&K2,&a,&oma);
5963 if(!(*t <= 0.0e0)) goto S10;
5973 void cumtnc(double *t,double *df,double *pnonc,double *cum,
5975 /**********************************************************************
5977 void cumtnc(double *t,double *df,double *pnonc,double *cum,
5980 CUMulative Non-Central T-distribution
5986 Computes the integral from -infinity to T of the non-central
5993 T --> Upper limit of integration of the non-central t-density.
5995 DF --> Degrees of freedom of the non-central t-distribution.
5997 PNONC --> Non-centrality parameter of the non-central t distibutio
5999 CUM <-- Cumulative t-distribution.
6001 CCUM <-- Compliment of Cumulative t-distribution.
6006 Upper tail of the cumulative noncentral t using
6007 formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuous
6008 Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
6010 This implementation starts the calculation at i = lambda,
6011 which is near the largest Di. It then sums forward and backward.
6012 **********************************************************************/
6020 #define tiny 1.0e-10
6021 static double alghdf,b,bb,bbcent,bcent,cent,d,dcent,dpnonc,dum1,dum2,e,ecent,
6022 halfdf,lambda,lnomx,lnx,omx,pnonc2,s,scent,ss,sscent,t2,term,tt,twoi,x,xi,
6025 static unsigned long qrevs;
6026 static double T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
6029 .. Executable Statements ..
6032 Case pnonc essentially zero
6034 if(fabs(*pnonc) <= tiny) {
6035 cumt(t,df,cum,ccum);
6047 pnonc2 = dpnonc * dpnonc;
6049 if(fabs(tt) <= tiny) {
6051 cumnor(&T1,cum,ccum);
6054 lambda = half * pnonc2;
6055 x = *df / (*df + t2);
6059 halfdf = half * *df;
6060 alghdf = gamln(&halfdf);
6062 ******************** Case i = lambda
6064 cent = fifidint(lambda);
6065 if(cent < one) cent = one;
6067 Compute d=T(2i) in log space and offset by exp(-lambda)
6070 xlnd = cent * log(lambda) - gamln(&T2) - lambda;
6073 Compute e=t(2i+1) in log space offset by exp(-lambda)
6076 xlne = (cent + half) * log(lambda) - gamln(&T3) - lambda;
6078 if(dpnonc < zero) ecent = -ecent;
6080 Compute bcent=B(2*cent)
6083 bratio(&halfdf,&T4,&x,&omx,&bcent,&dum1,&ierr);
6085 compute bbcent=B(2*cent+1)
6088 bratio(&halfdf,&T5,&x,&omx,&bbcent,&dum2,&ierr);
6090 Case bcent and bbcent are essentially zero
6091 Thus t is effectively infinite
6093 if(bcent + bbcent < tiny) {
6105 Case bcent and bbcent are essentially one
6106 Thus t is effectively zero
6108 if(dum1 + dum2 < tiny) {
6110 cumnor(&T6,cum,ccum);
6114 First term in ccum is D*B + E*BB
6116 *ccum = dcent * bcent + ecent * bbcent;
6118 compute s(cent) = B(2*(cent+1)) - B(2*cent))
6120 T7 = halfdf + cent + half;
6122 scent = gamln(&T7) - gamln(&T8) - alghdf + halfdf * lnx + (cent + half) *
6126 compute ss(cent) = B(2*cent+3) - B(2*cent+1)
6128 T9 = halfdf + cent + one;
6130 sscent = gamln(&T9) - gamln(&T10) - alghdf + halfdf * lnx + (cent + one) *
6132 sscent = exp(sscent);
6134 ******************** Sum Forward
6147 d = lambda / xi * d;
6148 e = lambda / (xi + half) * e;
6149 term = d * b + e * bb;
6151 s = s * omx * (*df + twoi - one) / (twoi + one);
6152 ss = ss * omx * (*df + twoi) / (twoi + two);
6155 if(fabs(term) > conv * *ccum) goto S10;
6157 ******************** Sum Backward
6165 s = scent * (one + twoi) / ((*df + twoi - one) * omx);
6166 ss = sscent * (two + twoi) / ((*df + twoi) * omx);
6171 e *= ((xi + half) / lambda);
6172 term = d * b + e * bb;
6175 if(xi < half) goto S30;
6177 s = s * (one + twoi) / ((*df + twoi - one) * omx);
6178 ss = ss * (two + twoi) / ((*df + twoi) * omx);
6179 if(fabs(term) > conv * *ccum) goto S20;
6182 *cum = half * *ccum;
6186 *ccum = half * *ccum;
6190 Due to roundoff error the answer may not lie between zero and one
6193 *cum = fifdmax1(fifdmin1(*cum,one),zero);
6194 *ccum = fifdmax1(fifdmin1(*ccum,one),zero);
6204 double devlpl(double a[],int *n,double *x)
6206 **********************************************************************
6208 double devlpl(double a[],int *n,double *x)
6209 Double precision EVALuate a PoLynomial at X
6216 A(1) + A(2)*X + ... + A(N)*X**(N-1)
6222 A --> Array of coefficients of the polynomial.
6223 A is DOUBLE PRECISION(N)
6225 N --> Length of A, also degree of polynomial - 1.
6228 X --> Point at which the polynomial is to be evaluated.
6229 X is DOUBLE PRECISION
6231 **********************************************************************
6234 static double devlpl,term;
6238 .. Executable Statements ..
6241 for(i= *n-1-1; i>=0; i--) term = a[i]+term**x;
6245 double dinvnr(double *p,double *q)
6247 **********************************************************************
6249 double dinvnr(double *p,double *q)
6250 Double precision NoRmal distribution INVerse
6256 Returns X such that CUMNOR(X) = P, i.e., the integral from -
6257 infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
6263 P --> The probability whose normal deviate is sought.
6264 P is DOUBLE PRECISION
6267 P is DOUBLE PRECISION
6273 The rational function on page 95 of Kennedy and Gentle,
6274 Statistical Computing, Marcel Dekker, NY , 1980 is used as a start
6275 value for the Newton method of finding roots.
6281 If P or Q .lt. machine EPS returns +/- DINVNR(EPS)
6283 **********************************************************************
6288 #define r2pi 0.3989422804014326e0
6289 #define nhalf -0.5e0
6290 #define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
6291 static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
6293 static unsigned long qporq;
6296 .. Executable Statements ..
6299 FIND MINIMUM OF P AND Q
6302 if(!qporq) goto S10;
6311 strtx = stvaln(&pp);
6316 for(i=1; i<=maxit; i++) {
6317 cumnor(&xcur,&cum,&ccum);
6318 dx = (cum-pp)/dennor(xcur);
6320 if(fabs(dx/xcur) < eps) goto S40;
6324 IF WE GET HERE, NEWTON HAS FAILED
6326 if(!qporq) dinvnr = -dinvnr;
6330 IF WE GET HERE, NEWTON HAS SUCCEDED
6333 if(!qporq) dinvnr = -dinvnr;
6342 static void E0000(int IENTRY,int *status,double *x,double *fx,
6343 unsigned long *qleft,unsigned long *qhi,double *zabsst,
6344 double *zabsto,double *zbig,double *zrelst,
6345 double *zrelto,double *zsmall,double *zstpmu)
6347 #define qxmon(zx,zy,zz) (int)((zx) <= (zy) && (zy) <= (zz))
6348 static double absstp,abstol,big,fbig,fsmall,relstp,reltol,small,step,stpmul,xhi,
6349 xlb,xlo,xsave,xub,yy;
6351 static unsigned long qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup;
6352 switch(IENTRY){case 0: goto DINVR; case 1: goto DSTINV;}
6354 if(*status > 0) goto S310;
6355 qcond = !qxmon(small,*x,big);
6364 See that SMALL and BIG bound the zero and set QINCR
6382 qincr = fbig > fsmall;
6383 if(!qincr) goto S50;
6384 if(fsmall <= 0.0e0) goto S30;
6389 if(fbig >= 0.0e0) goto S40;
6396 if(fsmall >= 0.0e0) goto S60;
6402 if(fbig <= 0.0e0) goto S70;
6410 step = fifdmax1(absstp,relstp*fabs(*x));
6419 if(!(yy == 0.0e0)) goto S100;
6424 qup = (qincr && yy < 0.0e0) || (!qincr && yy > 0.0e0);
6426 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6427 HANDLE CASE IN WHICH WE MUST STEP HIGHER
6428 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6432 xub = fifdmin1(xlb+step,big);
6435 if(qcond) goto S150;
6448 qbdd = (qincr && yy >= 0.0e0) || (!qincr && yy <= 0.0e0);
6450 qcond = qbdd || qlim;
6451 if(qcond) goto S140;
6454 xub = fifdmin1(xlb+step,big);
6458 if(!(qlim && !qbdd)) goto S160;
6468 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6469 HANDLE CASE IN WHICH WE MUST STEP LOWER
6470 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6473 xlb = fifdmax1(xub-step,small);
6476 if(qcond) goto S220;
6489 qbdd = (qincr && yy <= 0.0e0) || (!qincr && yy >= 0.0e0);
6490 qlim = xlb <= small;
6491 qcond = qbdd || qlim;
6492 if(qcond) goto S210;
6495 xlb = fifdmax1(xub-step,small);
6499 if(!(qlim && !qbdd)) goto S230;
6507 dstzr(&xlb,&xub,&abstol,&reltol);
6509 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6510 IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
6511 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6516 if(!(*status == 1)) goto S290;
6518 dzror(status,x,fx,&xlo,&xhi,&qdum1,&qdum2);
6519 if(!(*status == 1)) goto S280;
6543 TO GET-FUNCTION-VALUE
6548 switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S90;case
6549 4: goto S130;case 5: goto S200;case 6: goto S270;default: break;}
6552 void dinvr(int *status,double *x,double *fx,
6553 unsigned long *qleft,unsigned long *qhi)
6555 **********************************************************************
6557 void dinvr(int *status,double *x,double *fx,
6558 unsigned long *qleft,unsigned long *qhi)
6561 bounds the zero of the function and invokes zror
6562 Reverse Communication
6568 Bounds the function and invokes ZROR to perform the zero
6569 finding. STINVR must have been called before this routine
6570 in order to set its parameters.
6576 STATUS <--> At the beginning of a zero finding problem, STATUS
6577 should be set to 0 and INVR invoked. (The value
6578 of parameters other than X will be ignored on this cal
6580 When INVR needs the function evaluated, it will set
6581 STATUS to 1 and return. The value of the function
6582 should be set in FX and INVR again called without
6583 changing any of its other parameters.
6585 When INVR has finished without error, it will return
6586 with STATUS 0. In that case X is approximately a root
6589 If INVR cannot bound the function, it returns status
6590 -1 and sets QLEFT and QHI.
6593 X <-- The value of X at which F(X) is to be evaluated.
6596 FX --> The value of F(X) calculated when INVR returns with
6600 QLEFT <-- Defined only if QMFINV returns .FALSE. In that
6601 case it is .TRUE. If the stepping search terminated
6602 unsucessfully at SMALL. If it is .FALSE. the search
6603 terminated unsucessfully at BIG.
6606 QHI <-- Defined only if QMFINV returns .FALSE. In that
6607 case it is .TRUE. if F(X) .GT. Y at the termination
6608 of the search and .FALSE. if F(X) .LT. Y at the
6609 termination of the search.
6612 **********************************************************************
6615 E0000(0,status,x,fx,qleft,qhi,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
6617 void dstinv(double *zsmall,double *zbig,double *zabsst,
6618 double *zrelst,double *zstpmu,double *zabsto,
6621 **********************************************************************
6622 void dstinv(double *zsmall,double *zbig,double *zabsst,
6623 double *zrelst,double *zstpmu,double *zabsto,
6626 Double Precision - SeT INverse finder - Reverse Communication
6628 Concise Description - Given a monotone function F finds X
6629 such that F(X) = Y. Uses Reverse communication -- see invr.
6630 This routine sets quantities needed by INVR.
6631 More Precise Description of INVR -
6632 F must be a monotone function, the results of QMFINV are
6633 otherwise undefined. QINCR must be .TRUE. if F is non-
6634 decreasing and .FALSE. if F is non-increasing.
6635 QMFINV will return .TRUE. if and only if F(SMALL) and
6636 F(BIG) bracket Y, i. e.,
6637 QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
6638 QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
6639 if QMFINV returns .TRUE., then the X returned satisfies
6640 the following condition. let
6641 TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
6642 then if QINCR is .TRUE.,
6643 F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
6644 and if QINCR is .FALSE.
6645 F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
6647 SMALL --> The left endpoint of the interval to be
6648 searched for a solution.
6649 SMALL is DOUBLE PRECISION
6650 BIG --> The right endpoint of the interval to be
6651 searched for a solution.
6652 BIG is DOUBLE PRECISION
6653 ABSSTP, RELSTP --> The initial step size in the search
6654 is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
6655 ABSSTP is DOUBLE PRECISION
6656 RELSTP is DOUBLE PRECISION
6657 STPMUL --> When a step doesn't bound the zero, the step
6658 size is multiplied by STPMUL and another step
6659 taken. A popular value is 2.0
6660 DOUBLE PRECISION STPMUL
6661 ABSTOL, RELTOL --> Two numbers that determine the accuracy
6662 of the solution. See function for a precise definition.
6663 ABSTOL is DOUBLE PRECISION
6664 RELTOL is DOUBLE PRECISION
6666 Compares F(X) with Y for the input value of X then uses QINCR
6667 to determine whether to step left or right to bound the
6668 desired x. the initial step size is
6669 MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
6670 Iteratively steps right or left until it bounds X.
6671 At each step which doesn't bound X, the step size is doubled.
6672 The routine is careful never to step beyond SMALL or BIG. If
6673 it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
6674 after setting QLEFT and QHI.
6675 If X is successfully bounded then Algorithm R of the paper
6676 'Two Efficient Algorithms with Guaranteed Convergence for
6677 Finding a Zero of a Function' by J. C. P. Bus and
6678 T. J. Dekker in ACM Transactions on Mathematical
6679 Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
6680 to find the zero of the function F(X)-Y. This is routine
6682 **********************************************************************
6686 E0000(1,&status,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
6688 assert(status == 0 );
6690 double dt1(double *p,double *q,double *df)
6692 **********************************************************************
6694 double dt1(double *p,double *q,double *df)
6695 Double precision Initalize Approximation to
6696 INVerse of the cumulative T distribution
6702 Returns the inverse of the T distribution function, i.e.,
6703 the integral from 0 to INVT of the T density is P. This is an
6704 initial approximation
6710 P --> The p-value whose inverse from the T distribution is
6712 P is DOUBLE PRECISION
6715 Q is DOUBLE PRECISION
6717 DF --> Degrees of freedom of the T distribution.
6718 DF is DOUBLE PRECISION
6720 **********************************************************************
6723 static double coef[4][5] = {
6724 {1.0e0,1.0e0,0.0e0,0.0e0,0.0e0},
6725 {3.0e0,16.0e0,5.0e0,0.0e0,0.0e0},
6726 {-15.0e0,17.0e0,19.0e0,3.0e0,0.0e0},
6727 {-945.0e0,-1920.0e0,1482.0e0,776.0e0,79.0e0}
6729 static double denom[4] = {
6730 4.0e0,96.0e0,384.0e0,92160.0e0
6732 static int ideg[4] = {
6735 static double dt1,denpow,sum,term,x,xp,xx;
6739 .. Executable Statements ..
6741 x = fabs(dinvnr(p,q));
6745 for(i=0; i<4; i++) {
6746 term = devlpl(&coef[i][0],&ideg[i],&xx)*x;
6748 sum += (term/(denpow*denom[i]));
6750 if(!(*p >= 0.5e0)) goto S20;
6760 static void E0001(int IENTRY,int *status,double *x,double *fx,
6761 double *xlo,double *xhi,unsigned long *qleft,
6762 unsigned long *qhi,double *zabstl,double *zreltl,
6763 double *zxhi,double *zxlo)
6765 #define ftol(zx) (0.5e0*fifdmax1(abstol,reltol*fabs((zx))))
6766 static double a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q,reltol,tol,w,xxhi,xxlo;
6767 static int ext,i99999;
6768 static unsigned long first,qrzero;
6769 switch(IENTRY){case 0: goto DZROR; case 1: goto DSTZR;}
6771 if(*status > 0) goto S280;
6791 Check that F(ZXLO) < 0 < F(ZXHI) or
6792 F(ZXLO) > 0 > F(ZXHI)
6794 if(!(fb < 0.0e0)) goto S40;
6795 if(!(*fx < 0.0e0)) goto S30;
6802 if(!(fb > 0.0e0)) goto S60;
6803 if(!(*fx > 0.0e0)) goto S50;
6817 if(!(fabs(fc) < fabs(fb))) goto S100;
6818 if(!(c != a)) goto S90;
6833 if(!(fabs(mb) > tol)) goto S240;
6834 if(!(ext > 3)) goto S110;
6838 tol = fifdsign(tol,mb);
6840 if(!first) goto S120;
6845 fdb = (fd-fb)/(d-b);
6846 fda = (fd-fa)/(d-a);
6850 if(!(p < 0.0e0)) goto S140;
6854 if(ext == 3) p *= 2.0e0;
6855 if(!(p*1.0e0 == 0.0e0 || p <= q*tol)) goto S150;
6859 if(!(p < mb*q)) goto S160;
6881 if(!(fc*fb >= 0.0e0)) goto S210;
6884 if(!(w == mb)) goto S220;
6893 qrzero = (fc >= 0.0e0 && fb <= 0.0e0) || (fc < 0.0e0 && fb >= 0.0e0);
6894 if(!qrzero) goto S250;
6909 TO GET-FUNCTION-VALUE
6914 switch((int)i99999){case 1: goto S10;case 2: goto S20;case 3: goto S200;
6918 void dzror(int *status,double *x,double *fx,double *xlo,
6919 double *xhi,unsigned long *qleft,unsigned long *qhi)
6921 **********************************************************************
6923 void dzror(int *status,double *x,double *fx,double *xlo,
6924 double *xhi,unsigned long *qleft,unsigned long *qhi)
6926 Double precision ZeRo of a function -- Reverse Communication
6932 Performs the zero finding. STZROR must have been called before
6933 this routine in order to set its parameters.
6939 STATUS <--> At the beginning of a zero finding problem, STATUS
6940 should be set to 0 and ZROR invoked. (The value
6941 of other parameters will be ignored on this call.)
6943 When ZROR needs the function evaluated, it will set
6944 STATUS to 1 and return. The value of the function
6945 should be set in FX and ZROR again called without
6946 changing any of its other parameters.
6948 When ZROR has finished without error, it will return
6949 with STATUS 0. In that case (XLO,XHI) bound the answe
6951 If ZROR finds an error (which implies that F(XLO)-Y an
6952 F(XHI)-Y have the same sign, it returns STATUS -1. In
6953 this case, XLO and XHI are undefined.
6956 X <-- The value of X at which F(X) is to be evaluated.
6959 FX --> The value of F(X) calculated when ZROR returns with
6963 XLO <-- When ZROR returns with STATUS = 0, XLO bounds the
6964 inverval in X containing the solution below.
6965 DOUBLE PRECISION XLO
6967 XHI <-- When ZROR returns with STATUS = 0, XHI bounds the
6968 inverval in X containing the solution above.
6969 DOUBLE PRECISION XHI
6971 QLEFT <-- .TRUE. if the stepping search terminated unsucessfully
6972 at XLO. If it is .FALSE. the search terminated
6973 unsucessfully at XHI.
6976 QHI <-- .TRUE. if F(X) .GT. Y at the termination of the
6977 search and .FALSE. if F(X) .LT. Y at the
6978 termination of the search.
6981 **********************************************************************
6984 E0001(0,status,x,fx,xlo,xhi,qleft,qhi,NULL,NULL,NULL,NULL);
6986 void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl)
6988 **********************************************************************
6989 void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl)
6990 Double precision SeT ZeRo finder - Reverse communication version
6992 Sets quantities needed by ZROR. The function of ZROR
6993 and the quantities set is given here.
6994 Concise Description - Given a function F
6995 find XLO such that F(XLO) = 0.
6996 More Precise Description -
6997 Input condition. F is a double precision function of a single
6998 double precision argument and XLO and XHI are such that
6999 F(XLO)*F(XHI) .LE. 0.0
7000 If the input condition is met, QRZERO returns .TRUE.
7001 and output values of XLO and XHI satisfy the following
7002 F(XLO)*F(XHI) .LE. 0.
7003 ABS(F(XLO) .LE. ABS(F(XHI)
7004 ABS(XLO-XHI) .LE. TOL(X)
7006 TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
7007 If this algorithm does not find XLO and XHI satisfying
7008 these conditions then QRZERO returns .FALSE. This
7009 implies that the input condition was not met.
7011 XLO --> The left endpoint of the interval to be
7012 searched for a solution.
7013 XLO is DOUBLE PRECISION
7014 XHI --> The right endpoint of the interval to be
7016 XHI is DOUBLE PRECISION
7017 ABSTOL, RELTOL --> Two numbers that determine the accuracy
7018 of the solution. See function for a
7020 ABSTOL is DOUBLE PRECISION
7021 RELTOL is DOUBLE PRECISION
7023 Algorithm R of the paper 'Two Efficient Algorithms with
7024 Guaranteed Convergence for Finding a Zero of a Function'
7025 by J. C. P. Bus and T. J. Dekker in ACM Transactions on
7026 Mathematical Software, Volume 1, no. 4 page 330
7027 (Dec. '75) is employed to find the zero of F(X)-Y.
7028 **********************************************************************
7031 E0001(1,NULL,NULL,NULL,NULL,NULL,NULL,NULL,zabstl,zreltl,zxhi,zxlo);
7033 double erf1(double *x)
7035 -----------------------------------------------------------------------
7036 EVALUATION OF THE REAL ERROR FUNCTION
7037 -----------------------------------------------------------------------
7040 static double c = .564189583547756e0;
7041 static double a[5] = {
7042 .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
7043 .479137145607681e-01,.128379167095513e+00
7045 static double b[3] = {
7046 .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
7048 static double p[8] = {
7049 -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
7050 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
7051 4.51918953711873e+02,3.00459261020162e+02
7053 static double q[8] = {
7054 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
7055 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7056 7.90950925327898e+02,3.00459260956983e+02
7058 static double r[5] = {
7059 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
7060 4.65807828718470e+00,2.82094791773523e-01
7062 static double s[4] = {
7063 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
7064 1.80124575948747e+01
7066 static double erf1,ax,bot,t,top,x2;
7069 .. Executable Statements ..
7072 if(ax > 0.5e0) goto S10;
7074 top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
7075 bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
7076 erf1 = *x*(top/bot);
7079 if(ax > 4.0e0) goto S20;
7080 top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7082 bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7084 erf1 = 0.5e0+(0.5e0-exp(-(*x**x))*top/bot);
7085 if(*x < 0.0e0) erf1 = -erf1;
7088 if(ax >= 5.8e0) goto S30;
7091 top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
7092 bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
7093 erf1 = (c-top/(x2*bot))/ax;
7094 erf1 = 0.5e0+(0.5e0-exp(-x2)*erf1);
7095 if(*x < 0.0e0) erf1 = -erf1;
7098 erf1 = fifdsign(1.0e0,*x);
7101 double erfc1(int *ind,double *x)
7103 -----------------------------------------------------------------------
7104 EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
7106 ERFC1(IND,X) = ERFC(X) IF IND = 0
7107 ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE
7108 -----------------------------------------------------------------------
7111 static double c = .564189583547756e0;
7112 static double a[5] = {
7113 .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
7114 .479137145607681e-01,.128379167095513e+00
7116 static double b[3] = {
7117 .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
7119 static double p[8] = {
7120 -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
7121 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
7122 4.51918953711873e+02,3.00459261020162e+02
7124 static double q[8] = {
7125 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
7126 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7127 7.90950925327898e+02,3.00459260956983e+02
7129 static double r[5] = {
7130 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
7131 4.65807828718470e+00,2.82094791773523e-01
7133 static double s[4] = {
7134 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
7135 1.80124575948747e+01
7138 static double erfc1,ax,bot,e,t,top,w;
7141 .. Executable Statements ..
7147 if(ax > 0.5e0) goto S10;
7149 top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
7150 bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
7151 erfc1 = 0.5e0+(0.5e0-*x*(top/bot));
7152 if(*ind != 0) erfc1 = exp(t)*erfc1;
7156 0.5 .LT. ABS(X) .LE. 4
7158 if(ax > 4.0e0) goto S20;
7159 top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7161 bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7169 if(*x <= -5.6e0) goto S60;
7170 if(*ind != 0) goto S30;
7171 if(*x > 100.0e0) goto S70;
7172 if(*x**x > -exparg(&K1)) goto S70;
7174 t = pow(1.0e0/ *x,2.0);
7175 top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
7176 bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
7177 erfc1 = (c-t*top/bot)/ax;
7182 if(*ind == 0) goto S50;
7183 if(*x < 0.0e0) erfc1 = 2.0e0*exp(*x**x)-erfc1;
7189 erfc1 = (0.5e0+(0.5e0-e))*exp(-t)*erfc1;
7190 if(*x < 0.0e0) erfc1 = 2.0e0-erfc1;
7194 LIMIT VALUE FOR LARGE NEGATIVE X
7197 if(*ind != 0) erfc1 = 2.0e0*exp(*x**x);
7201 LIMIT VALUE FOR LARGE POSITIVE X
7207 double esum(int *mu,double *x)
7209 -----------------------------------------------------------------------
7210 EVALUATION OF EXP(MU + X)
7211 -----------------------------------------------------------------------
7214 static double esum,w;
7217 .. Executable Statements ..
7219 if(*x > 0.0e0) goto S10;
7220 if(*mu < 0) goto S20;
7222 if(w > 0.0e0) goto S20;
7226 if(*mu > 0) goto S20;
7228 if(w < 0.0e0) goto S20;
7233 esum = exp(w)*exp(*x);
7236 double exparg(int *l)
7238 --------------------------------------------------------------------
7239 IF L = 0 THEN EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
7240 EXP(W) CAN BE COMPUTED.
7242 IF L IS NONZERO THEN EXPARG(L) = THE LARGEST NEGATIVE W FOR
7243 WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.
7245 NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.
7246 --------------------------------------------------------------------
7252 static double exparg,lnb;
7256 .. Executable Statements ..
7259 if(b != 2) goto S10;
7260 lnb = .69314718055995e0;
7263 if(b != 8) goto S20;
7264 lnb = 2.0794415416798e0;
7267 if(b != 16) goto S30;
7268 lnb = 2.7725887222398e0;
7271 lnb = log((double)b);
7273 if(*l == 0) goto S50;
7275 exparg = 0.99999e0*((double)m*lnb);
7279 exparg = 0.99999e0*((double)m*lnb);
7282 double fpser(double *a,double *b,double *x,double *eps)
7284 -----------------------------------------------------------------------
7286 EVALUATION OF I (A,B)
7289 FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
7291 -----------------------------------------------------------------------
7297 static double fpser,an,c,s,t,tol;
7300 .. Executable Statements ..
7303 if(*a <= 1.e-3**eps) goto S10;
7306 if(t < exparg(&K1)) return fpser;
7310 NOTE THAT 1/B(A,B) = B
7312 fpser = *b/ *a*fpser;
7322 if(fabs(c) > tol) goto S20;
7323 fpser *= (1.0e0+*a*s);
7326 double gam1(double *a)
7328 ------------------------------------------------------------------
7329 COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5
7330 ------------------------------------------------------------------
7333 static double s1 = .273076135303957e+00;
7334 static double s2 = .559398236957378e-01;
7335 static double p[7] = {
7336 .577215664901533e+00,-.409078193005776e+00,-.230975380857675e+00,
7337 .597275330452234e-01,.766968181649490e-02,-.514889771323592e-02,
7338 .589597428611429e-03
7340 static double q[5] = {
7341 .100000000000000e+01,.427569613095214e+00,.158451672430138e+00,
7342 .261132021441447e-01,.423244297896961e-02
7344 static double r[9] = {
7345 -.422784335098468e+00,-.771330383816272e+00,-.244757765222226e+00,
7346 .118378989872749e+00,.930357293360349e-03,-.118290993445146e-01,
7347 .223047661158249e-02,.266505979058923e-03,-.132674909766242e-03
7349 static double gam1,bot,d,t,top,w,T1;
7352 .. Executable Statements ..
7356 if(d > 0.0e0) t = d-0.5e0;
7358 if(T1 < 0) goto S40;
7359 else if(T1 == 0) goto S10;
7365 top = (((((p[6]*t+p[5])*t+p[4])*t+p[3])*t+p[2])*t+p[1])*t+p[0];
7366 bot = (((q[4]*t+q[3])*t+q[2])*t+q[1])*t+1.0e0;
7368 if(d > 0.0e0) goto S30;
7372 gam1 = t/ *a*(w-0.5e0-0.5e0);
7375 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+
7377 bot = (s2*t+s1)*t+1.0e0;
7379 if(d > 0.0e0) goto S50;
7380 gam1 = *a*(w+0.5e0+0.5e0);
7386 void gaminv(double *a,double *x,double *x0,double *p,double *q,
7389 ----------------------------------------------------------------------
7390 INVERSE INCOMPLETE GAMMA RATIO FUNCTION
7392 GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1.
7393 THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER
7394 ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X
7395 TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE
7396 PARTICULAR COMPUTER ARITHMETIC BEING USED.
7400 X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0,
7401 AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT
7402 NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN
7403 A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE
7404 IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X.
7406 X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER
7407 DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET
7410 IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
7411 WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING
7414 IERR = 0 THE SOLUTION WAS OBTAINED. ITERATION WAS
7416 IERR.GT.0 THE SOLUTION WAS OBTAINED. IERR ITERATIONS
7418 IERR = -2 (INPUT ERROR) A .LE. 0
7419 IERR = -3 NO SOLUTION WAS OBTAINED. THE RATIO Q/A
7421 IERR = -4 (INPUT ERROR) P + Q .NE. 1
7422 IERR = -6 20 ITERATIONS WERE PERFORMED. THE MOST
7423 RECENT VALUE OBTAINED FOR X IS GIVEN.
7424 THIS CANNOT OCCUR IF X0 .LE. 0.
7425 IERR = -7 ITERATION FAILED. NO VALUE IS GIVEN FOR X.
7426 THIS MAY OCCUR WHEN X IS APPROXIMATELY 0.
7427 IERR = -8 A VALUE FOR X HAS BEEN OBTAINED, BUT THE
7428 ROUTINE IS NOT CERTAIN OF ITS ACCURACY.
7429 ITERATION CANNOT BE PERFORMED IN THIS
7430 CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY
7431 WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS
7432 POSITIVE THEN THIS CAN OCCUR WHEN A IS
7433 EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY
7434 LARGE (SAY A .GE. 1.E20).
7435 ----------------------------------------------------------------------
7436 WRITTEN BY ALFRED H. MORRIS, JR.
7437 NAVAL SURFACE WEAPONS CENTER
7442 static double a0 = 3.31125922108741e0;
7443 static double a1 = 11.6616720288968e0;
7444 static double a2 = 4.28342155967104e0;
7445 static double a3 = .213623493715853e0;
7446 static double b1 = 6.61053765625462e0;
7447 static double b2 = 6.40691597760039e0;
7448 static double b3 = 1.27364489782223e0;
7449 static double b4 = .036117081018842e0;
7450 static double c = .577215664901533e0;
7451 static double ln10 = 2.302585e0;
7452 static double tol = 1.e-5;
7453 static double amin[2] = {
7456 static double bmin[2] = {
7459 static double dmin[2] = {
7462 static double emin[2] = {
7465 static double eps0[2] = {
7472 static double am1,amax,ap1,ap2,ap3,apn,b,c1,c2,c3,c4,c5,d,e,e2,eps,g,h,pn,qg,qn,
7473 r,rta,s,s2,sum,t,u,w,xmax,xmin,xn,y,z;
7475 static double T4,T5,T6,T7,T9;
7478 .. Executable Statements ..
7481 ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
7482 E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
7483 XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
7484 LARGEST POSITIVE NUMBER.
7490 if(*a <= 0.0e0) goto S300;
7492 if(fabs(t) > e) goto S320;
7494 if(*p == 0.0e0) return;
7495 if(*q == 0.0e0) goto S270;
7496 if(*a == 1.0e0) goto S280;
7498 amax = 0.4e-10/(e*e);
7500 if(e > 1.e-10) iop = 2;
7503 if(*x0 > 0.0e0) goto S160;
7505 SELECTION OF THE INITIAL APPROXIMATION XN OF X
7508 if(*a > 1.0e0) goto S80;
7512 if(qg == 0.0e0) goto S360;
7514 if(qg > 0.6e0**a) goto S40;
7515 if(*a >= 0.30e0 || b < 0.35e0) goto S10;
7521 if(b >= 0.45e0) goto S40;
7522 if(b == 0.0e0) goto S360;
7524 s = 0.5e0+(0.5e0-*a);
7527 if(b < 0.15e0) goto S20;
7528 xn = y-s*log(t)-log(1.0e0+s/(t+1.0e0));
7531 if(b <= 0.01e0) goto S30;
7532 u = ((t+2.0e0*(3.0e0-*a))*t+(2.0e0-*a)*(3.0e0-*a))/((t+(5.0e0-*a))*t+2.0e0);
7533 xn = y-s*log(t)-log(u);
7537 c2 = -(s*(1.0e0+c1));
7538 c3 = s*((0.5e0*c1+(2.0e0-*a))*c1+(2.5e0-1.5e0**a));
7539 c4 = -(s*(((c1/3.0e0+(2.5e0-1.5e0**a))*c1+((*a-6.0e0)**a+7.0e0))*c1+(
7540 (11.0e0**a-46.0)**a+47.0e0)/6.0e0));
7541 c5 = -(s*((((-(c1/4.0e0)+(11.0e0**a-17.0e0)/6.0e0)*c1+((-(3.0e0**a)+13.0e0)*
7542 *a-13.0e0))*c1+0.5e0*(((2.0e0**a-25.0e0)**a+72.0e0)**a-61.0e0))*c1+((
7543 (25.0e0**a-195.0e0)**a+477.0e0)**a-379.0e0)/12.0e0));
7544 xn = (((c5/y+c4)/y+c3)/y+c2)/y+c1+y;
7545 if(*a > 1.0e0) goto S220;
7546 if(b > bmin[iop-1]) goto S220;
7550 if(b**q > 1.e-8) goto S50;
7551 xn = exp(-(*q/ *a+c));
7554 if(*p <= 0.9e0) goto S60;
7556 xn = exp((alnrel(&T5)+gamln1(a))/ *a);
7559 xn = exp(log(*p*g)/ *a);
7561 if(xn == 0.0e0) goto S310;
7562 t = 0.5e0+(0.5e0-xn/(*a+1.0e0));
7567 SELECTION OF THE INITIAL APPROXIMATION XN OF X
7570 if(*q <= 0.5e0) goto S90;
7576 t = sqrt(-(2.0e0*w));
7577 s = t-(((a3*t+a2)*t+a1)*t+a0)/((((b4*t+b3)*t+b2)*t+b1)*t+1.0e0);
7578 if(*q > 0.5e0) s = -s;
7581 xn = *a+s*rta+(s2-1.0e0)/3.0e0+s*(s2-7.0e0)/(36.0e0*rta)-((3.0e0*s2+7.0e0)*
7582 s2-16.0e0)/(810.0e0**a)+s*((9.0e0*s2+256.0e0)*s2-433.0e0)/(38880.0e0**a*
7584 xn = fifdmax1(xn,0.0e0);
7585 if(*a < amin[iop-1]) goto S110;
7587 d = 0.5e0+(0.5e0-*x/ *a);
7588 if(fabs(d) <= dmin[iop-1]) return;
7590 if(*p <= 0.5e0) goto S130;
7591 if(xn < 3.0e0**a) goto S220;
7593 d = fifdmax1(2.0e0,*a*(*a-1.0e0));
7594 if(y < ln10*d) goto S120;
7600 T6 = -(t/(xn+1.0e0));
7601 xn = y+t*log(xn)-alnrel(&T6);
7602 T7 = -(t/(xn+1.0e0));
7603 xn = y+t*log(xn)-alnrel(&T7);
7607 if(xn > 0.70e0*ap1) goto S170;
7609 if(xn > 0.15e0*ap1) goto S140;
7612 *x = exp((w+*x)/ *a);
7613 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
7614 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2)))/ *a);
7615 *x = exp((w+*x-log(1.0e0+*x/ap1*(1.0e0+*x/ap2*(1.0e0+*x/ap3))))/ *a);
7617 if(xn > 1.e-2*ap1) goto S140;
7618 if(xn <= emin[iop-1]*ap1) return;
7628 if(t > 1.e-4) goto S150;
7630 xn = exp((xn+t)/ *a);
7631 xn *= (1.0e0-(*a*log(xn)-xn-t)/(*a-xn));
7635 SCHRODER ITERATION USING P
7637 if(*p > 0.5e0) goto S220;
7639 if(*p <= 1.e10*xmin) goto S350;
7640 am1 = *a-0.5e0-0.5e0;
7642 if(*a <= amax) goto S190;
7643 d = 0.5e0+(0.5e0-xn/ *a);
7644 if(fabs(d) <= e2) goto S350;
7646 if(*ierr >= 20) goto S330;
7648 gratio(a,&xn,&pn,&qn,&K8);
7649 if(pn == 0.0e0 || qn == 0.0e0) goto S350;
7651 if(r == 0.0e0) goto S350;
7654 if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S200;
7656 if(*x <= 0.0e0) goto S340;
7662 if(*x <= 0.0e0) goto S340;
7663 if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
7667 if(d > tol) goto S180;
7668 if(d <= eps) return;
7669 if(fabs(*p-pn) <= tol**p) return;
7673 SCHRODER ITERATION USING Q
7675 if(*q <= 1.e10*xmin) goto S350;
7676 am1 = *a-0.5e0-0.5e0;
7678 if(*a <= amax) goto S240;
7679 d = 0.5e0+(0.5e0-xn/ *a);
7680 if(fabs(d) <= e2) goto S350;
7682 if(*ierr >= 20) goto S330;
7684 gratio(a,&xn,&pn,&qn,&K8);
7685 if(pn == 0.0e0 || qn == 0.0e0) goto S350;
7687 if(r == 0.0e0) goto S350;
7690 if(fabs(t) <= 0.1e0 && fabs(w*t) <= 0.1e0) goto S250;
7692 if(*x <= 0.0e0) goto S340;
7698 if(*x <= 0.0e0) goto S340;
7699 if(fabs(w) >= 1.0e0 && fabs(w)*t*t <= eps) return;
7703 if(d > tol) goto S230;
7704 if(d <= eps) return;
7705 if(fabs(*q-qn) <= tol**q) return;
7714 if(*q < 0.9e0) goto S290;
7748 double gamln(double *a)
7750 -----------------------------------------------------------------------
7751 EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
7752 -----------------------------------------------------------------------
7753 WRITTEN BY ALFRED H. MORRIS
7754 NAVAL SURFACE WARFARE CENTER
7756 --------------------------
7757 D = 0.5*(LN(2*PI) - 1)
7758 --------------------------
7761 static double c0 = .833333333333333e-01;
7762 static double c1 = -.277777777760991e-02;
7763 static double c2 = .793650666825390e-03;
7764 static double c3 = -.595202931351870e-03;
7765 static double c4 = .837308034031215e-03;
7766 static double c5 = -.165322962780713e-02;
7767 static double d = .418938533204673e0;
7768 static double gamln,t,w;
7773 .. Executable Statements ..
7775 if(*a > 0.8e0) goto S10;
7776 gamln = gamln1(a)-log(*a);
7779 if(*a > 2.25e0) goto S20;
7784 if(*a >= 10.0e0) goto S40;
7785 n = (long)(*a - 1.25e0);
7788 for(i=1; i<=n; i++) {
7793 gamln = gamln1(&T1)+log(w);
7796 t = pow(1.0e0/ *a,2.0);
7797 w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a;
7798 gamln = d+w+(*a-0.5e0)*(log(*a)-1.0e0);
7801 double gamln1(double *a)
7803 -----------------------------------------------------------------------
7804 EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
7805 -----------------------------------------------------------------------
7808 static double p0 = .577215664901533e+00;
7809 static double p1 = .844203922187225e+00;
7810 static double p2 = -.168860593646662e+00;
7811 static double p3 = -.780427615533591e+00;
7812 static double p4 = -.402055799310489e+00;
7813 static double p5 = -.673562214325671e-01;
7814 static double p6 = -.271935708322958e-02;
7815 static double q1 = .288743195473681e+01;
7816 static double q2 = .312755088914843e+01;
7817 static double q3 = .156875193295039e+01;
7818 static double q4 = .361951990101499e+00;
7819 static double q5 = .325038868253937e-01;
7820 static double q6 = .667465618796164e-03;
7821 static double r0 = .422784335098467e+00;
7822 static double r1 = .848044614534529e+00;
7823 static double r2 = .565221050691933e+00;
7824 static double r3 = .156513060486551e+00;
7825 static double r4 = .170502484022650e-01;
7826 static double r5 = .497958207639485e-03;
7827 static double s1 = .124313399877507e+01;
7828 static double s2 = .548042109832463e+00;
7829 static double s3 = .101552187439830e+00;
7830 static double s4 = .713309612391000e-02;
7831 static double s5 = .116165475989616e-03;
7832 static double gamln1,w,x;
7835 .. Executable Statements ..
7837 if(*a >= 0.6e0) goto S10;
7838 w = ((((((p6**a+p5)**a+p4)**a+p3)**a+p2)**a+p1)**a+p0)/((((((q6**a+q5)**a+
7839 q4)**a+q3)**a+q2)**a+q1)**a+1.0e0);
7844 w = (((((r5*x+r4)*x+r3)*x+r2)*x+r1)*x+r0)/(((((s5*x+s4)*x+s3)*x+s2)*x+s1)*x
7849 double Xgamm(double *a)
7851 -----------------------------------------------------------------------
7853 EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS
7857 GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
7860 -----------------------------------------------------------------------
7861 WRITTEN BY ALFRED H. MORRIS, JR.
7862 NAVAL SURFACE WEAPONS CENTER
7864 -----------------------------------------------------------------------
7867 static double d = .41893853320467274178e0;
7868 static double pi = 3.1415926535898e0;
7869 static double r1 = .820756370353826e-03;
7870 static double r2 = -.595156336428591e-03;
7871 static double r3 = .793650663183693e-03;
7872 static double r4 = -.277777777770481e-02;
7873 static double r5 = .833333333333333e-01;
7874 static double p[7] = {
7875 .539637273585445e-03,.261939260042690e-02,.204493667594920e-01,
7876 .730981088720487e-01,.279648642639792e+00,.553413866010467e+00,1.0e0
7878 static double q[7] = {
7879 -.832979206704073e-03,.470059485860584e-02,.225211131035340e-01,
7880 -.170458969313360e+00,-.567902761974940e-01,.113062953091122e+01,1.0e0
7884 static double Xgamm,bot,g,lnx,s,t,top,w,x,z;
7885 static int i,j,m,n,T1;
7888 .. Executable Statements ..
7892 if(fabs(*a) >= 15.0e0) goto S110;
7894 -----------------------------------------------------------------------
7895 EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15
7896 -----------------------------------------------------------------------
7901 LET T BE THE PRODUCT OF A-J WHEN A .GE. 2
7904 if(T1 < 0) goto S40;
7905 else if(T1 == 0) goto S30;
7908 for(j=1; j<=m; j++) {
7917 LET T BE THE PRODUCT OF A+J WHEN A .LT. 1
7920 if(*a > 0.0e0) goto S70;
7922 if(m == 0) goto S60;
7923 for(j=1; j<=m; j++) {
7930 if(t == 0.0e0) return Xgamm;
7933 THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
7934 CODE MAY BE OMITTED IF DESIRED.
7936 if(fabs(t) >= 1.e-30) goto S80;
7937 if(fabs(t)*spmpar(&K2) <= 1.0001e0) return Xgamm;
7942 COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1
7946 for(i=1; i<7; i++) {
7954 if(*a < 1.0e0) goto S100;
7962 -----------------------------------------------------------------------
7963 EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15
7964 -----------------------------------------------------------------------
7966 if(fabs(*a) >= 1.e3) return Xgamm;
7967 if(*a > 0.0e0) goto S120;
7971 if(t > 0.9e0) t = 1.0e0-t;
7973 if(fifmod(n,2) == 0) s = -s;
7974 if(s == 0.0e0) return Xgamm;
7977 COMPUTE THE MODIFIED ASYMPTOTIC SUM
7980 g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x;
7982 ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X)
7983 BUT LESS ACCURACY WILL NORMALLY BE OBTAINED.
7990 g = d+g+(z-0.5e0)*(lnx-1.e0);
7993 if(w > 0.99999e0*exparg(&K3)) return Xgamm;
7994 Xgamm = exp(w)*(1.0e0+t);
7995 if(*a < 0.0e0) Xgamm = 1.0e0/(Xgamm*s)/x;
7998 void grat1(double *a,double *x,double *r,double *p,double *q,
8002 static double a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,t,tol,w,z,T1,T3;
8005 .. Executable Statements ..
8008 -----------------------------------------------------------------------
8009 EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
8011 IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED.
8012 THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A).
8013 -----------------------------------------------------------------------
8015 if(*a**x == 0.0e0) goto S120;
8016 if(*a == 0.5e0) goto S100;
8017 if(*x < 1.1e0) goto S10;
8021 TAYLOR SERIES FOR P(A,X)/X**A
8025 sum = *x/(*a+3.0e0);
8026 tol = 0.1e0**eps/(*a+1.0e0);
8032 if(fabs(t) > tol) goto S20;
8033 j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
8037 if(*x < 0.25e0) goto S30;
8038 if(*a < *x/2.59e0) goto S50;
8041 if(z > -.13394e0) goto S50;
8044 *p = w*g*(0.5e0+(0.5e0-j));
8045 *q = 0.5e0+(0.5e0-*p);
8049 w = 0.5e0+(0.5e0+l);
8051 if(*q < 0.0e0) goto S90;
8052 *p = 0.5e0+(0.5e0-*q);
8056 CONTINUED FRACTION EXPANSION
8058 a2nm1 = a2n = 1.0e0;
8060 b2n = *x+(1.0e0-*a);
8063 a2nm1 = *x*a2n+c*a2nm1;
8064 b2nm1 = *x*b2n+c*b2nm1;
8068 a2n = a2nm1+cma*a2n;
8069 b2n = b2nm1+cma*b2n;
8071 if(fabs(an0-am0) >= *eps*an0) goto S70;
8073 *p = 0.5e0+(0.5e0-*q);
8087 if(*x >= 0.25e0) goto S110;
8090 *q = 0.5e0+(0.5e0-*p);
8094 *q = erfc1(&K2,&T3);
8095 *p = 0.5e0+(0.5e0-*q);
8098 if(*x <= *a) goto S80;
8101 void gratio(double *a,double *x,double *ans,double *qans,int *ind)
8103 ----------------------------------------------------------------------
8104 EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
8109 IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
8112 ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
8113 P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
8114 IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
8115 POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
8116 IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
8117 6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
8118 IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.
8121 ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
8122 WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
8123 P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
8124 X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.
8125 ----------------------------------------------------------------------
8126 WRITTEN BY ALFRED H. MORRIS, JR.
8127 NAVAL SURFACE WEAPONS CENTER
8129 --------------------
8132 static double alog10 = 2.30258509299405e0;
8133 static double d10 = -.185185185185185e-02;
8134 static double d20 = .413359788359788e-02;
8135 static double d30 = .649434156378601e-03;
8136 static double d40 = -.861888290916712e-03;
8137 static double d50 = -.336798553366358e-03;
8138 static double d60 = .531307936463992e-03;
8139 static double d70 = .344367606892378e-03;
8140 static double rt2pin = .398942280401433e0;
8141 static double rtpi = 1.77245385090552e0;
8142 static double third = .333333333333333e0;
8143 static double acc0[3] = {
8146 static double big[3] = {
8147 20.0e0,14.0e0,10.0e0
8149 static double d0[13] = {
8150 .833333333333333e-01,-.148148148148148e-01,.115740740740741e-02,
8151 .352733686067019e-03,-.178755144032922e-03,.391926317852244e-04,
8152 -.218544851067999e-05,-.185406221071516e-05,.829671134095309e-06,
8153 -.176659527368261e-06,.670785354340150e-08,.102618097842403e-07,
8154 -.438203601845335e-08
8156 static double d1[12] = {
8157 -.347222222222222e-02,.264550264550265e-02,-.990226337448560e-03,
8158 .205761316872428e-03,-.401877572016461e-06,-.180985503344900e-04,
8159 .764916091608111e-05,-.161209008945634e-05,.464712780280743e-08,
8160 .137863344691572e-06,-.575254560351770e-07,.119516285997781e-07
8162 static double d2[10] = {
8163 -.268132716049383e-02,.771604938271605e-03,.200938786008230e-05,
8164 -.107366532263652e-03,.529234488291201e-04,-.127606351886187e-04,
8165 .342357873409614e-07,.137219573090629e-05,-.629899213838006e-06,
8166 .142806142060642e-06
8168 static double d3[8] = {
8169 .229472093621399e-03,-.469189494395256e-03,.267720632062839e-03,
8170 -.756180167188398e-04,-.239650511386730e-06,.110826541153473e-04,
8171 -.567495282699160e-05,.142309007324359e-05
8173 static double d4[6] = {
8174 .784039221720067e-03,-.299072480303190e-03,-.146384525788434e-05,
8175 .664149821546512e-04,-.396836504717943e-04,.113757269706784e-04
8177 static double d5[4] = {
8178 -.697281375836586e-04,.277275324495939e-03,-.199325705161888e-03,
8179 .679778047793721e-04
8181 static double d6[2] = {
8182 -.592166437353694e-03,.270878209671804e-03
8184 static double e00[3] = {
8187 static double x00[3] = {
8192 static double a2n,a2nm1,acc,am0,amn,an,an0,apn,b2n,b2nm1,c,c0,c1,c2,c3,c4,c5,c6,
8193 cma,e,e0,g,h,j,l,r,rta,rtx,s,sum,t,t1,tol,twoa,u,w,x0,y,z;
8194 static int i,iop,m,max,n;
8195 static double wk[20],T3;
8197 static double T6,T7;
8200 .. Executable Statements ..
8203 --------------------
8204 ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
8205 FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
8208 if(*a < 0.0e0 || *x < 0.0e0) goto S430;
8209 if(*a == 0.0e0 && *x == 0.0e0) goto S430;
8210 if(*a**x == 0.0e0) goto S420;
8212 if(iop != 1 && iop != 2) iop = 3;
8213 acc = fifdmax1(acc0[iop-1],e);
8217 SELECT THE APPROPRIATE ALGORITHM
8219 if(*a >= 1.0e0) goto S10;
8220 if(*a == 0.5e0) goto S390;
8221 if(*x < 1.1e0) goto S160;
8224 if(u == 0.0e0) goto S380;
8225 r = u*(1.0e0+gam1(a));
8228 if(*a >= big[iop-1]) goto S30;
8229 if(*a > *x || *x >= x0) goto S20;
8232 if(twoa != (double)m) goto S20;
8234 if(*a == (double)i) goto S210;
8238 r = exp(t1)/Xgamm(a);
8242 if(l == 0.0e0) goto S370;
8243 s = 0.5e0+(0.5e0-l);
8245 if(z >= 700.0e0/ *a) goto S410;
8248 if(fabs(s) <= e0/rta) goto S330;
8249 if(fabs(s) <= 0.4e0) goto S270;
8250 t = pow(1.0e0/ *a,2.0);
8251 t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
8253 r = rt2pin*rta*exp(t1);
8255 if(r == 0.0e0) goto S420;
8256 if(*x <= fifdmax1(*a,alog10)) goto S50;
8257 if(*x < x0) goto S250;
8261 TAYLOR SERIES FOR P/R
8266 for(n=2; n<=20; n++) {
8269 if(t <= 1.e-3) goto S70;
8280 if(t > tol) goto S80;
8282 for(m=1; m<=max; m++) {
8286 *ans = r/ *a*(1.0e0+sum);
8287 *qans = 0.5e0+(0.5e0-*ans);
8291 ASYMPTOTIC EXPANSION
8296 for(n=2; n<=20; n++) {
8299 if(fabs(t) <= 1.e-3) goto S120;
8306 if(fabs(t) <= acc) goto S140;
8313 for(m=1; m<=max; m++) {
8317 *qans = r/ *x*(1.0e0+sum);
8318 *ans = 0.5e0+(0.5e0-*qans);
8322 TAYLOR SERIES FOR P(A,X)/X**A
8326 sum = *x/(*a+3.0e0);
8327 tol = 3.0e0*acc/(*a+1.0e0);
8333 if(fabs(t) > tol) goto S170;
8334 j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
8338 if(*x < 0.25e0) goto S180;
8339 if(*a < *x/2.59e0) goto S200;
8342 if(z > -.13394e0) goto S200;
8345 *ans = w*g*(0.5e0+(0.5e0-j));
8346 *qans = 0.5e0+(0.5e0-*ans);
8350 w = 0.5e0+(0.5e0+l);
8351 *qans = (w*j-l)*g-h;
8352 if(*qans < 0.0e0) goto S380;
8353 *ans = 0.5e0+(0.5e0-*qans);
8357 FINITE SUMS FOR Q WHEN A .GE. 1
8358 AND 2*A IS AN INTEGER
8367 sum = erfc1(&K2,&rtx);
8368 t = exp(-*x)/(rtpi*rtx);
8372 if(n == i) goto S240;
8380 *ans = 0.5e0+(0.5e0-*qans);
8384 CONTINUED FRACTION EXPANSION
8386 tol = fifdmax1(5.0e0*e,acc);
8387 a2nm1 = a2n = 1.0e0;
8389 b2n = *x+(1.0e0-*a);
8392 a2nm1 = *x*a2n+c*a2nm1;
8393 b2nm1 = *x*b2n+c*b2nm1;
8397 a2n = a2nm1+cma*a2n;
8398 b2n = b2nm1+cma*b2n;
8400 if(fabs(an0-am0) >= tol*an0) goto S260;
8402 *ans = 0.5e0+(0.5e0-*qans);
8406 GENERAL TEMME EXPANSION
8408 if(fabs(s) <= 2.0e0*e && *a*e*e > 3.28e-3) goto S430;
8411 w = 0.5e0*erfc1(&K1,&T3);
8414 if(l < 1.0e0) z = -z;
8416 if(T4 < 0) goto S280;
8417 else if(T4 == 0) goto S290;
8420 if(fabs(s) <= 1.e-3) goto S340;
8421 c0 = ((((((((((((d0[12]*z+d0[11])*z+d0[10])*z+d0[9])*z+d0[8])*z+d0[7])*z+d0[
8422 6])*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
8423 c1 = (((((((((((d1[11]*z+d1[10])*z+d1[9])*z+d1[8])*z+d1[7])*z+d1[6])*z+d1[5]
8424 )*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8425 c2 = (((((((((d2[9]*z+d2[8])*z+d2[7])*z+d2[6])*z+d2[5])*z+d2[4])*z+d2[3])*z+
8426 d2[2])*z+d2[1])*z+d2[0])*z+d20;
8427 c3 = (((((((d3[7]*z+d3[6])*z+d3[5])*z+d3[4])*z+d3[3])*z+d3[2])*z+d3[1])*z+
8429 c4 = (((((d4[5]*z+d4[4])*z+d4[3])*z+d4[2])*z+d4[1])*z+d4[0])*z+d40;
8430 c5 = (((d5[3]*z+d5[2])*z+d5[1])*z+d5[0])*z+d50;
8431 c6 = (d6[1]*z+d6[0])*z+d60;
8432 t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
8435 c0 = (((((d0[5]*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-third;
8436 c1 = (((d1[3]*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8441 t = ((d0[2]*z+d0[1])*z+d0[0])*z-third;
8443 if(l < 1.0e0) goto S320;
8444 *qans = c*(w+rt2pin*t/rta);
8445 *ans = 0.5e0+(0.5e0-*qans);
8448 *ans = c*(w-rt2pin*t/rta);
8449 *qans = 0.5e0+(0.5e0-*ans);
8453 TEMME EXPANSION FOR L = 1
8455 if(*a*e*e > 3.28e-3) goto S430;
8456 c = 0.5e0+(0.5e0-y);
8457 w = (0.5e0-sqrt(y)*(0.5e0+(0.5e0-y/3.0e0))/rtpi)/c;
8460 if(l < 1.0e0) z = -z;
8462 if(T5 < 0) goto S340;
8463 else if(T5 == 0) goto S350;
8466 c0 = ((((((d0[6]*z+d0[5])*z+d0[4])*z+d0[3])*z+d0[2])*z+d0[1])*z+d0[0])*z-
8468 c1 = (((((d1[5]*z+d1[4])*z+d1[3])*z+d1[2])*z+d1[1])*z+d1[0])*z+d10;
8469 c2 = ((((d2[4]*z+d2[3])*z+d2[2])*z+d2[1])*z+d2[0])*z+d20;
8470 c3 = (((d3[3]*z+d3[2])*z+d3[1])*z+d3[0])*z+d30;
8471 c4 = (d4[1]*z+d4[0])*z+d40;
8472 c5 = (d5[1]*z+d5[0])*z+d50;
8474 t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u+c0;
8477 c0 = (d0[1]*z+d0[0])*z-third;
8479 t = (d20*u+c1)*u+c0;
8496 if(*x >= 0.25e0) goto S400;
8499 *qans = 0.5e0+(0.5e0-*ans);
8503 *qans = erfc1(&K2,&T7);
8504 *ans = 0.5e0+(0.5e0-*qans);
8507 if(fabs(s) <= 2.0e0*e) goto S430;
8509 if(*x <= *a) goto S370;
8518 double gsumln(double *a,double *b)
8520 -----------------------------------------------------------------------
8521 EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
8522 FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2
8523 -----------------------------------------------------------------------
8526 static double gsumln,x,T1,T2;
8529 .. Executable Statements ..
8532 if(x > 0.25e0) goto S10;
8534 gsumln = gamln1(&T1);
8537 if(x > 1.25e0) goto S20;
8538 gsumln = gamln1(&x)+alnrel(&x);
8542 gsumln = gamln1(&T2)+log(x*(1.0e0+x));
8545 double psi(double *xx)
8547 ---------------------------------------------------------------------
8549 EVALUATION OF THE DIGAMMA FUNCTION
8553 PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
8556 THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
8557 APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
8558 CODY, STRECOK AND THACHER.
8560 ---------------------------------------------------------------------
8561 PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
8562 PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
8564 ---------------------------------------------------------------------
8567 static double dx0 = 1.461632144968362341262659542325721325e0;
8568 static double piov4 = .785398163397448e0;
8569 static double p1[7] = {
8570 .895385022981970e-02,.477762828042627e+01,.142441585084029e+03,
8571 .118645200713425e+04,.363351846806499e+04,.413810161269013e+04,
8572 .130560269827897e+04
8574 static double p2[4] = {
8575 -.212940445131011e+01,-.701677227766759e+01,-.448616543918019e+01,
8576 -.648157123766197e+00
8578 static double q1[6] = {
8579 .448452573429826e+02,.520752771467162e+03,.221000799247830e+04,
8580 .364127349079381e+04,.190831076596300e+04,.691091682714533e-05
8582 static double q2[4] = {
8583 .322703493791143e+02,.892920700481861e+02,.546117738103215e+02,
8584 .777788548522962e+01
8588 static double psi,aug,den,sgn,upper,w,x,xmax1,xmx0,xsmall,z;
8589 static int i,m,n,nq;
8592 .. Executable Statements ..
8595 ---------------------------------------------------------------------
8596 MACHINE DEPENDENT CONSTANTS ...
8597 XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
8598 WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED
8599 AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
8600 ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
8601 PSI MAY BE REPRESENTED AS ALOG(X).
8602 XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
8603 MAY BE REPRESENTED BY 1/X.
8604 ---------------------------------------------------------------------
8606 xmax1 = ipmpar(&K1);
8607 xmax1 = fifdmin1(xmax1,1.0e0/spmpar(&K2));
8611 if(x >= 0.5e0) goto S50;
8613 ---------------------------------------------------------------------
8614 X .LT. 0.5, USE REFLECTION FORMULA
8615 PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
8616 ---------------------------------------------------------------------
8618 if(fabs(x) > xsmall) goto S10;
8619 if(x == 0.0e0) goto S100;
8621 ---------------------------------------------------------------------
8622 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE
8624 ---------------------------------------------------------------------
8630 ---------------------------------------------------------------------
8631 REDUCTION OF ARGUMENT FOR COTAN
8632 ---------------------------------------------------------------------
8636 if(w > 0.0e0) goto S20;
8641 ---------------------------------------------------------------------
8642 MAKE AN ERROR EXIT IF X .LE. -XMAX1
8643 ---------------------------------------------------------------------
8645 if(w >= xmax1) goto S100;
8648 nq = fifidint(w*4.0e0);
8649 w = 4.0e0*(w-(double)nq*.25e0);
8651 ---------------------------------------------------------------------
8652 W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X.
8653 ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
8654 QUADRANT AND DETERMINE SIGN
8655 ---------------------------------------------------------------------
8658 if(n+n != nq) w = 1.0e0-w;
8661 if(m+m != n) sgn = -sgn;
8663 ---------------------------------------------------------------------
8664 DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X)
8665 ---------------------------------------------------------------------
8670 if(m != n) goto S30;
8672 ---------------------------------------------------------------------
8673 CHECK FOR SINGULARITY
8674 ---------------------------------------------------------------------
8676 if(z == 0.0e0) goto S100;
8678 ---------------------------------------------------------------------
8679 USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
8680 SIN/COS AS A SUBSTITUTE FOR TAN
8681 ---------------------------------------------------------------------
8683 aug = sgn*(cos(z)/sin(z)*4.0e0);
8686 aug = sgn*(sin(z)/cos(z)*4.0e0);
8690 if(x > 3.0e0) goto S70;
8692 ---------------------------------------------------------------------
8694 ---------------------------------------------------------------------
8698 for(i=1; i<=5; i++) {
8699 den = (den+q1[i-1])*x;
8700 upper = (upper+p1[i+1-1])*x;
8702 den = (upper+p1[6])/(den+q1[5]);
8708 ---------------------------------------------------------------------
8709 IF X .GE. XMAX1, PSI = LN(X)
8710 ---------------------------------------------------------------------
8712 if(x >= xmax1) goto S90;
8714 ---------------------------------------------------------------------
8715 3.0 .LT. X .LT. XMAX1
8716 ---------------------------------------------------------------------
8721 for(i=1; i<=3; i++) {
8722 den = (den+q2[i-1])*w;
8723 upper = (upper+p2[i+1-1])*w;
8725 aug = upper/(den+q2[3])-0.5e0/x+aug;
8731 ---------------------------------------------------------------------
8733 ---------------------------------------------------------------------
8738 double rcomp(double *a,double *x)
8741 EVALUATION OF EXP(-X)*X**A/GAMMA(A)
8743 RT2PIN = 1/SQRT(2*PI)
8747 static double rt2pin = .398942280401433e0;
8748 static double rcomp,t,t1,u;
8751 .. Executable Statements ..
8754 if(*a >= 20.0e0) goto S20;
8756 if(*a >= 1.0e0) goto S10;
8757 rcomp = *a*exp(t)*(1.0e0+gam1(a));
8760 rcomp = exp(t)/Xgamm(a);
8764 if(u == 0.0e0) return rcomp;
8765 t = pow(1.0e0/ *a,2.0);
8766 t1 = (((0.75e0*t-1.0e0)*t+3.5e0)*t-105.0e0)/(*a*1260.0e0);
8767 t1 -= (*a*rlog(&u));
8768 rcomp = rt2pin*sqrt(*a)*exp(t1);
8771 double rexp(double *x)
8773 -----------------------------------------------------------------------
8774 EVALUATION OF THE FUNCTION EXP(X) - 1
8775 -----------------------------------------------------------------------
8778 static double p1 = .914041914819518e-09;
8779 static double p2 = .238082361044469e-01;
8780 static double q1 = -.499999999085958e+00;
8781 static double q2 = .107141568980644e+00;
8782 static double q3 = -.119041179760821e-01;
8783 static double q4 = .595130811860248e-03;
8784 static double rexp,w;
8787 .. Executable Statements ..
8789 if(fabs(*x) > 0.15e0) goto S10;
8790 rexp = *x*(((p2**x+p1)**x+1.0e0)/((((q4**x+q3)**x+q2)**x+q1)**x+1.0e0));
8794 if(*x > 0.0e0) goto S20;
8795 rexp = w-0.5e0-0.5e0;
8798 rexp = w*(0.5e0+(0.5e0-1.0e0/w));
8801 double rlog(double *x)
8804 COMPUTATION OF X - 1 - LN(X)
8808 static double a = .566749439387324e-01;
8809 static double b = .456512608815524e-01;
8810 static double p0 = .333333333333333e+00;
8811 static double p1 = -.224696413112536e+00;
8812 static double p2 = .620886815375787e-02;
8813 static double q1 = -.127408923933623e+01;
8814 static double q2 = .354508718369557e+00;
8815 static double rlog,r,t,u,w,w1;
8818 .. Executable Statements ..
8820 if(*x < 0.61e0 || *x > 1.57e0) goto S40;
8821 if(*x < 0.82e0) goto S10;
8822 if(*x > 1.18e0) goto S20;
8843 w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
8844 rlog = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
8851 double rlog1(double *x)
8853 -----------------------------------------------------------------------
8854 EVALUATION OF THE FUNCTION X - LN(1 + X)
8855 -----------------------------------------------------------------------
8858 static double a = .566749439387324e-01;
8859 static double b = .456512608815524e-01;
8860 static double p0 = .333333333333333e+00;
8861 static double p1 = -.224696413112536e+00;
8862 static double p2 = .620886815375787e-02;
8863 static double q1 = -.127408923933623e+01;
8864 static double q2 = .354508718369557e+00;
8865 static double rlog1,h,r,t,w,w1;
8868 .. Executable Statements ..
8870 if(*x < -0.39e0 || *x > 0.57e0) goto S40;
8871 if(*x < -0.18e0) goto S10;
8872 if(*x > 0.18e0) goto S20;
8885 h = 0.75e0**x-0.25e0;
8893 w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0);
8894 rlog1 = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1;
8901 double spmpar(int *i)
8903 -----------------------------------------------------------------------
8905 SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
8906 THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
8907 I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
8908 SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
8909 ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
8911 SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
8913 SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
8915 SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
8917 -----------------------------------------------------------------------
8919 ALFRED H. MORRIS, JR.
8920 NAVAL SURFACE WARFARE CENTER
8922 -----------------------------------------------------------------------
8923 -----------------------------------------------------------------------
8924 MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE
8925 CONSTANTS FOR THE COMPUTER BEING USED. THIS MODIFICATION WAS
8926 MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION
8927 -----------------------------------------------------------------------
8934 static double spmpar,b,binv,bm1,one,w,z;
8935 static int emax,emin,ibeta,m;
8938 .. Executable Statements ..
8940 if(*i > 1) goto S10;
8943 spmpar = pow(b,(double)(1-m));
8946 if(*i > 2) goto S20;
8951 w = pow(b,(double)(emin+2));
8952 spmpar = w*binv*binv*binv;
8955 ibeta = ipmpar(&K1);
8961 z = pow(b,(double)(m-1));
8962 w = ((z-one)*b+bm1)/(b*z);
8963 z = pow(b,(double)(emax-2));
8967 double stvaln(double *p)
8969 **********************************************************************
8971 double stvaln(double *p)
8972 STarting VALue for Neton-Raphon
8973 calculation of Normal distribution Inverse
8979 Returns X such that CUMNOR(X) = P, i.e., the integral from -
8980 infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
8986 P --> The probability whose normal deviate is sought.
8987 P is DOUBLE PRECISION
8993 The rational function on page 95 of Kennedy and Gentle,
8994 Statistical Computing, Marcel Dekker, NY , 1980.
8996 **********************************************************************
8999 static double xden[5] = {
9000 0.993484626060e-1,0.588581570495e0,0.531103462366e0,0.103537752850e0,
9003 static double xnum[5] = {
9004 -0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1,
9008 static double stvaln,sign,y,z;
9011 .. Executable Statements ..
9013 if(!(*p <= 0.5e0)) goto S10;
9021 y = sqrt(-(2.0e0*log(z)));
9022 stvaln = y+devlpl(xnum,&K1,&y)/devlpl(xden,&K1,&y);
9023 stvaln = sign*stvaln;
9026 /************************************************************************
9028 Truncates a double precision number to an integer and returns the
9030 ************************************************************************/
9031 double fifdint(double a)
9032 /* a - number to be truncated */
9036 return (double)(temp);
9038 /************************************************************************
9040 returns the maximum of two numbers a and b
9041 ************************************************************************/
9042 double fifdmax1(double a,double b)
9043 /* a - first number */
9044 /* b - second number */
9046 if (a < b) return b;
9049 /************************************************************************
9051 returns the minimum of two numbers a and b
9052 ************************************************************************/
9053 double fifdmin1(double a,double b)
9054 /* a - first number */
9055 /* b - second number */
9057 if (a < b) return a;
9060 /************************************************************************
9062 transfers the sign of the variable "sign" to the variable "mag"
9063 ************************************************************************/
9064 double fifdsign(double mag,double sign)
9065 /* mag - magnitude */
9066 /* sign - sign to be transfered */
9068 if (mag < 0) mag = -mag;
9069 if (sign < 0) mag = -mag;
9073 /************************************************************************
9075 Truncates a double precision number to a long integer
9076 ************************************************************************/
9077 long fifidint(double a)
9078 /* a - number to be truncated */
9082 /************************************************************************
9084 returns the modulo of a and b
9085 ************************************************************************/
9086 long fifmod(long a,long b)
9088 /* b - denominator */