/* *************************************************************************** FILE DCDFLIB.C (janvier 2004) Ce sous-programme est constitué à partir d'extraits de la librairie mise en ligne par le M.D. Anderson Cancer Center (University of Texas) : fttp://odin.mdacc.tmc.edu/pub/source/dcdflib.c Déchargée en juillet 2003 Seules les lignes nécessaires au programme de segmentation ont été conservées et adaptées à ce programme Gaetan Peaquin - gaetan.peaquin@free.fr Cyril Labbe - cyril.labbe@imag.fr Dominique Labbe - dominique.labbe@iep.upmf-grenoble.fr *************************************************************************** */ #include #include #include #include "dcdflib.h" #include "ipmpar.h" double algdiv(double *a,double *b) { static double c0 = .833333333333333e-01; static double c1 = -.277777777760991e-02; static double c2 = .793650666825390e-03; static double c3 = -.595202931351870e-03; static double c4 = .837308034031215e-03; static double c5 = -.165322962780713e-02; static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1; /* .. .. Executable Statements .. */ if(*a <= *b) goto S10; h = *b/ *a; c = 1.0e0/(1.0e0+h); x = h/(1.0e0+h); d = *a+(*b-0.5e0); goto S20; S10: h = *a/ *b; c = h/(1.0e0+h); x = 1.0e0/(1.0e0+h); d = *b+(*a-0.5e0); S20: /* SET SN = (1 - X**N)/(1 - X) */ x2 = x*x; s3 = 1.0e0+(x+x2); s5 = 1.0e0+(x+x2*s3); s7 = 1.0e0+(x+x2*s5); s9 = 1.0e0+(x+x2*s7); s11 = 1.0e0+(x+x2*s9); /* SET W = DEL(B) - DEL(A + B) */ t = pow(1.0e0/ *b,2.0); w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0; w *= (c/ *b); /* COMBINE THE RESULTS */ T1 = *a/ *b; u = d*alnrel(&T1); v = *a*(log(*b)-1.0e0); if(u <= v) goto S30; algdiv = w-v-u; return algdiv; S30: algdiv = w-u-v; return algdiv; } double alnrel(double *a) /* ----------------------------------------------------------------------- EVALUATION OF THE FUNCTION LN(1 + A) ----------------------------------------------------------------------- */ { static double p1 = -.129418923021993e+01; static double p2 = .405303492862024e+00; static double p3 = -.178874546012214e-01; static double q1 = -.162752256355323e+01; static double q2 = .747811014037616e+00; static double q3 = -.845104217945565e-01; static double alnrel,t,t2,w,x; /* .. .. Executable Statements .. */ if(fabs(*a) > 0.375e0) goto S10; t = *a/(*a+2.0e0); t2 = t*t; w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0); alnrel = 2.0e0*t*w; return alnrel; S10: x = 1.e0+*a; alnrel = log(x); return alnrel; } double apser(double *a,double *b,double *x,double *eps) /* ----------------------------------------------------------------------- APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED. ----------------------------------------------------------------------- */ { static double g = .577215664901533e0; static double apser,aj,bx,c,j,s,t,tol; /* .. .. Executable Statements .. */ bx = *b**x; t = *x-bx; if(*b**eps > 2.e-2) goto S10; c = log(*x)+psi(b)+g+t; goto S20; S10: c = log(bx)+g+t; S20: tol = 5.0e0**eps*fabs(c); j = 1.0e0; s = 0.0e0; S30: j += 1.0e0; t *= (*x-bx/j); aj = t/j; s += aj; if(fabs(aj) > tol) goto S30; apser = -(*a*(c+s)); return apser; } double basym(double *a,double *b,double *lambda,double *eps) /* ----------------------------------------------------------------------- ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B. LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED. IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT A AND B ARE GREATER THAN OR EQUAL TO 15. ----------------------------------------------------------------------- */ { static double e0 = 1.12837916709551e0; static double e1 = .353553390593274e0; static int num = 20; /* ------------------------ ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN. THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1. ------------------------ E0 = 2/SQRT(PI) E1 = 2**(-3/2) ------------------------ */ static int K3 = 1; static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,z2,zn,znm1; static int i,im1,imj,j,m,mm1,mmj,n,np1; static double a0[21],b0[21],c[21],d[21],T1,T2; /* .. .. Executable Statements .. */ basym = 0.0e0; if(*a >= *b) goto S10; h = *a/ *b; r0 = 1.0e0/(1.0e0+h); r1 = (*b-*a)/ *b; w0 = 1.0e0/sqrt(*a*(1.0e0+h)); goto S20; S10: h = *b/ *a; r0 = 1.0e0/(1.0e0+h); r1 = (*b-*a)/ *a; w0 = 1.0e0/sqrt(*b*(1.0e0+h)); S20: T1 = -(*lambda/ *a); T2 = *lambda/ *b; f = *a*rlog1(&T1)+*b*rlog1(&T2); t = exp(-f); if(t == 0.0e0) return basym; z0 = sqrt(f); z = 0.5e0*(z0/e1); z2 = f+f; a0[0] = 2.0e0/3.0e0*r1; c[0] = -(0.5e0*a0[0]); d[0] = -c[0]; j0 = 0.5e0/e0*erfc1(&K3,&z0); j1 = e1; sum = j0+d[0]*w0*j1; s = 1.0e0; h2 = h*h; hn = 1.0e0; w = w0; znm1 = z; zn = z2; for(n=2; n<=num; n+=2) { hn = h2*hn; a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0); np1 = n+1; s += hn; a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0); for(i=n; i<=np1; i++) { r = -(0.5e0*((double)i+1.0e0)); b0[0] = r*a0[0]; for(m=2; m<=i; m++) { bsum = 0.0e0; mm1 = m-1; for(j=1; j<=mm1; j++) { mmj = m-j; bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]); } b0[m-1] = r*a0[m-1]+bsum/(double)m; } c[i-1] = b0[i-1]/((double)i+1.0e0); dsum = 0.0e0; im1 = i-1; for(j=1; j<=im1; j++) { imj = i-j; dsum += (d[imj-1]*c[j-1]); } d[i-1] = -(dsum+c[i-1]); } j0 = e1*znm1+((double)n-1.0e0)*j0; j1 = e1*zn+(double)n*j1; znm1 = z2*znm1; zn = z2*zn; w = w0*w; t0 = d[n-1]*w*j0; w = w0*w; t1 = d[np1-1]*w*j1; sum += (t0+t1); if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80; } S80: u = exp(-bcorr(a,b)); basym = e0*t*u*sum; return basym; } double bcorr(double *a0,double *b0) /* ----------------------------------------------------------------------- EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A). IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8. ----------------------------------------------------------------------- */ { static double c0 = .833333333333333e-01; static double c1 = -.277777777760991e-02; static double c2 = .793650666825390e-03; static double c3 = -.595202931351870e-03; static double c4 = .837308034031215e-03; static double c5 = -.165322962780713e-02; static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2; /* .. .. Executable Statements .. */ a = fifdmin1(*a0,*b0); b = fifdmax1(*a0,*b0); h = a/b; c = h/(1.0e0+h); x = 1.0e0/(1.0e0+h); x2 = x*x; /* SET SN = (1 - X**N)/(1 - X) */ s3 = 1.0e0+(x+x2); s5 = 1.0e0+(x+x2*s3); s7 = 1.0e0+(x+x2*s5); s9 = 1.0e0+(x+x2*s7); s11 = 1.0e0+(x+x2*s9); /* SET W = DEL(B) - DEL(A + B) */ t = pow(1.0e0/b,2.0); w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0; w *= (c/b); /* COMPUTE DEL(A) + W */ t = pow(1.0e0/a,2.0); bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w; return bcorr; } double betaln(double *a0,double *b0) /* ----------------------------------------------------------------------- EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION ----------------------------------------------------------------------- E = 0.5*LN(2*PI) -------------------------- */ { static double e = .918938533204673e0; static double betaln,a,b,c,h,u,v,w,z; static int i,n; static double T1; /* .. .. Executable Statements .. */ a = fifdmin1(*a0,*b0); b = fifdmax1(*a0,*b0); if(a >= 8.0e0) goto S100; if(a >= 1.0e0) goto S20; /* ----------------------------------------------------------------------- PROCEDURE WHEN A .LT. 1 ----------------------------------------------------------------------- */ if(b >= 8.0e0) goto S10; T1 = a+b; betaln = gamln(&a)+(gamln(&b)-gamln(&T1)); return betaln; S10: betaln = gamln(&a)+algdiv(&a,&b); return betaln; S20: /* ----------------------------------------------------------------------- PROCEDURE WHEN 1 .LE. A .LT. 8 ----------------------------------------------------------------------- */ if(a > 2.0e0) goto S40; if(b > 2.0e0) goto S30; betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b); return betaln; S30: w = 0.0e0; if(b < 8.0e0) goto S60; betaln = gamln(&a)+algdiv(&a,&b); return betaln; S40: /* REDUCTION OF A WHEN B .LE. 1000 */ if(b > 1000.0e0) goto S80; n = (long)(a - 1.0e0); w = 1.0e0; for(i=1; i<=n; i++) { a -= 1.0e0; h = a/b; w *= (h/(1.0e0+h)); } w = log(w); if(b < 8.0e0) goto S60; betaln = w+gamln(&a)+algdiv(&a,&b); return betaln; S60: /* REDUCTION OF B WHEN B .LT. 8 */ n = (long)(b - 1.0e0); z = 1.0e0; for(i=1; i<=n; i++) { b -= 1.0e0; z *= (b/(a+b)); } betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b))); return betaln; S80: /* REDUCTION OF A WHEN B .GT. 1000 */ n = (long)(a - 1.0e0); w = 1.0e0; for(i=1; i<=n; i++) { a -= 1.0e0; w *= (a/(1.0e0+a/b)); } betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b)); return betaln; S100: /* ----------------------------------------------------------------------- PROCEDURE WHEN A .GE. 8 ----------------------------------------------------------------------- */ w = bcorr(&a,&b); h = a/b; c = h/(1.0e0+h); u = -((a-0.5e0)*log(c)); v = b*alnrel(&h); if(u <= v) goto S110; betaln = -(0.5e0*log(b))+e+w-v-u; return betaln; S110: betaln = -(0.5e0*log(b))+e+w-u-v; return betaln; } double bpser(double *a,double *b,double *x,double *eps) /* ----------------------------------------------------------------------- POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1 OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED. ----------------------------------------------------------------------- */ { static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z; static int i,m; /* .. .. Executable Statements .. */ bpser = 0.0e0; if(*x == 0.0e0) return bpser; /* ----------------------------------------------------------------------- COMPUTE THE FACTOR X**A/(A*BETA(A,B)) ----------------------------------------------------------------------- */ a0 = fifdmin1(*a,*b); if(a0 < 1.0e0) goto S10; z = *a*log(*x)-betaln(a,b); bpser = exp(z)/ *a; goto S100; S10: b0 = fifdmax1(*a,*b); if(b0 >= 8.0e0) goto S90; if(b0 > 1.0e0) goto S40; /* PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1 */ bpser = pow(*x,*a); if(bpser == 0.0e0) return bpser; apb = *a+*b; if(apb > 1.0e0) goto S20; z = 1.0e0+gam1(&apb); goto S30; S20: u = *a+*b-1.e0; z = (1.0e0+gam1(&u))/apb; S30: c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z; bpser *= (c*(*b/apb)); goto S100; S40: /* PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8 */ u = gamln1(&a0); m = (long)(b0 - 1.0e0); if(m < 1) goto S60; c = 1.0e0; for(i=1; i<=m; i++) { b0 -= 1.0e0; c *= (b0/(a0+b0)); } u = log(c)+u; S60: z = *a*log(*x)-u; b0 -= 1.0e0; apb = a0+b0; if(apb > 1.0e0) goto S70; t = 1.0e0+gam1(&apb); goto S80; S70: u = a0+b0-1.e0; t = (1.0e0+gam1(&u))/apb; S80: bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t; goto S100; S90: /* PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8 */ u = gamln1(&a0)+algdiv(&a0,&b0); z = *a*log(*x)-u; bpser = a0/ *a*exp(z); S100: if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser; /* ----------------------------------------------------------------------- COMPUTE THE SERIES ----------------------------------------------------------------------- */ sum = n = 0.0e0; c = 1.0e0; tol = *eps/ *a; S110: //printf("%i\n",__LINE__); n += 1.0e0; c *= ((0.5e0+(0.5e0-*b/n))**x); w = c/(*a+n); sum += w; //printf("fabs(w) = %f\n",w); if(fabs(w) > tol) goto S110; bpser *= (1.0e0+*a*sum); return bpser; } double brcomp(double *a,double *b,double *x,double *y) /* ----------------------------------------------------------------------- EVALUATION OF X**A*Y**B/BETA(A,B) ----------------------------------------------------------------------- */ { static double Const = .398942280401433e0; static double brcomp,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z; static int i,n; /* ----------------- CONST = 1/SQRT(2*PI) ----------------- */ static double T1,T2; /* .. .. Executable Statements .. */ brcomp = 0.0e0; if(*x == 0.0e0 || *y == 0.0e0) return brcomp; a0 = fifdmin1(*a,*b); if(a0 >= 8.0e0) goto S130; if(*x > 0.375e0) goto S10; lnx = log(*x); T1 = -*x; lny = alnrel(&T1); goto S30; S10: if(*y > 0.375e0) goto S20; T2 = -*y; lnx = alnrel(&T2); lny = log(*y); goto S30; S20: lnx = log(*x); lny = log(*y); S30: z = *a*lnx+*b*lny; if(a0 < 1.0e0) goto S40; z -= betaln(a,b); brcomp = exp(z); return brcomp; S40: /* ----------------------------------------------------------------------- PROCEDURE FOR A .LT. 1 OR B .LT. 1 ----------------------------------------------------------------------- */ b0 = fifdmax1(*a,*b); if(b0 >= 8.0e0) goto S120; if(b0 > 1.0e0) goto S70; /* ALGORITHM FOR B0 .LE. 1 */ brcomp = exp(z); if(brcomp == 0.0e0) return brcomp; apb = *a+*b; if(apb > 1.0e0) goto S50; z = 1.0e0+gam1(&apb); goto S60; S50: u = *a+*b-1.e0; z = (1.0e0+gam1(&u))/apb; S60: c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z; brcomp = brcomp*(a0*c)/(1.0e0+a0/b0); return brcomp; S70: /* ALGORITHM FOR 1 .LT. B0 .LT. 8 */ u = gamln1(&a0); n = (long)(b0 - 1.0e0); if(n < 1) goto S90; c = 1.0e0; for(i=1; i<=n; i++) { b0 -= 1.0e0; c *= (b0/(a0+b0)); } u = log(c)+u; S90: z -= u; b0 -= 1.0e0; apb = a0+b0; if(apb > 1.0e0) goto S100; t = 1.0e0+gam1(&apb); goto S110; S100: u = a0+b0-1.e0; t = (1.0e0+gam1(&u))/apb; S110: brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t; return brcomp; S120: /* ALGORITHM FOR B0 .GE. 8 */ u = gamln1(&a0)+algdiv(&a0,&b0); brcomp = a0*exp(z-u); return brcomp; S130: /* ----------------------------------------------------------------------- PROCEDURE FOR A .GE. 8 AND B .GE. 8 ----------------------------------------------------------------------- */ if(*a > *b) goto S140; h = *a/ *b; x0 = h/(1.0e0+h); y0 = 1.0e0/(1.0e0+h); lambda = *a-(*a+*b)**x; goto S150; S140: h = *b/ *a; x0 = 1.0e0/(1.0e0+h); y0 = h/(1.0e0+h); lambda = (*a+*b)**y-*b; S150: e = -(lambda/ *a); if(fabs(e) > 0.6e0) goto S160; u = rlog1(&e); goto S170; S160: u = e-log(*x/x0); S170: e = lambda/ *b; if(fabs(e) > 0.6e0) goto S180; v = rlog1(&e); goto S190; S180: v = e-log(*y/y0); S190: z = exp(-(*a*u+*b*v)); brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b)); return brcomp; } double brcmp1(int *mu,double *a,double *b,double *x,double *y) /* ----------------------------------------------------------------------- EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B)) ----------------------------------------------------------------------- */ { static double Const = .398942280401433e0; static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z; static int i,n; /* ----------------- CONST = 1/SQRT(2*PI) ----------------- */ static double T1,T2,T3,T4; /* .. .. Executable Statements .. */ a0 = fifdmin1(*a,*b); if(a0 >= 8.0e0) goto S130; if(*x > 0.375e0) goto S10; lnx = log(*x); T1 = -*x; lny = alnrel(&T1); goto S30; S10: if(*y > 0.375e0) goto S20; T2 = -*y; lnx = alnrel(&T2); lny = log(*y); goto S30; S20: lnx = log(*x); lny = log(*y); S30: z = *a*lnx+*b*lny; if(a0 < 1.0e0) goto S40; z -= betaln(a,b); brcmp1 = esum(mu,&z); return brcmp1; S40: /* ----------------------------------------------------------------------- PROCEDURE FOR A .LT. 1 OR B .LT. 1 ----------------------------------------------------------------------- */ b0 = fifdmax1(*a,*b); if(b0 >= 8.0e0) goto S120; if(b0 > 1.0e0) goto S70; /* ALGORITHM FOR B0 .LE. 1 */ brcmp1 = esum(mu,&z); if(brcmp1 == 0.0e0) return brcmp1; apb = *a+*b; if(apb > 1.0e0) goto S50; z = 1.0e0+gam1(&apb); goto S60; S50: u = *a+*b-1.e0; z = (1.0e0+gam1(&u))/apb; S60: c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z; brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0); return brcmp1; S70: /* ALGORITHM FOR 1 .LT. B0 .LT. 8 */ u = gamln1(&a0); n = (long)(b0 - 1.0e0); if(n < 1) goto S90; c = 1.0e0; for(i=1; i<=n; i++) { b0 -= 1.0e0; c *= (b0/(a0+b0)); } u = log(c)+u; S90: z -= u; b0 -= 1.0e0; apb = a0+b0; if(apb > 1.0e0) goto S100; t = 1.0e0+gam1(&apb); goto S110; S100: u = a0+b0-1.e0; t = (1.0e0+gam1(&u))/apb; S110: brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t; return brcmp1; S120: /* ALGORITHM FOR B0 .GE. 8 */ u = gamln1(&a0)+algdiv(&a0,&b0); T3 = z-u; brcmp1 = a0*esum(mu,&T3); return brcmp1; S130: /* ----------------------------------------------------------------------- PROCEDURE FOR A .GE. 8 AND B .GE. 8 ----------------------------------------------------------------------- */ if(*a > *b) goto S140; h = *a/ *b; x0 = h/(1.0e0+h); y0 = 1.0e0/(1.0e0+h); lambda = *a-(*a+*b)**x; goto S150; S140: h = *b/ *a; x0 = 1.0e0/(1.0e0+h); y0 = h/(1.0e0+h); lambda = (*a+*b)**y-*b; S150: e = -(lambda/ *a); if(fabs(e) > 0.6e0) goto S160; u = rlog1(&e); goto S170; S160: u = e-log(*x/x0); S170: e = lambda/ *b; if(fabs(e) > 0.6e0) goto S180; v = rlog1(&e); goto S190; S180: v = e-log(*y/y0); S190: T4 = -(*a*u+*b*v); z = esum(mu,&T4); brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b)); return brcmp1; } double bfrac(double *a,double *b,double *x,double *y,double *lambda,double *eps) /* ----------------------------------------------------------------------- CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1. IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B. ----------------------------------------------------------------------- */ { static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1; /* .. .. Executable Statements .. */ bfrac = brcomp(a,b,x,y); if(bfrac == 0.0e0) return bfrac; c = 1.0e0+*lambda; c0 = *b/ *a; c1 = 1.0e0+1.0e0/ *a; yp1 = *y+1.0e0; n = 0.0e0; p = 1.0e0; s = *a+1.0e0; an = 0.0e0; bn = anp1 = 1.0e0; bnp1 = c/c1; r = c1/c; S10: /* CONTINUED FRACTION CALCULATION */ n += 1.0e0; t = n/ *a; w = n*(*b-n)**x; e = *a/s; alpha = p*(p+c0)*e*e*(w**x); e = (1.0e0+t)/(c1+t+t); beta = n+w/s+e*(c+n*yp1); p = 1.0e0+t; s += 2.0e0; /* UPDATE AN, BN, ANP1, AND BNP1 */ t = alpha*an+beta*anp1; an = anp1; anp1 = t; t = alpha*bn+beta*bnp1; bn = bnp1; bnp1 = t; r0 = r; r = anp1/bnp1; if(fabs(r-r0) <= *eps*r) goto S20; /* RESCALE AN, BN, ANP1, AND BNP1 */ an /= bnp1; bn /= bnp1; anp1 = r; bnp1 = 1.0e0; goto S10; S20: /* TERMINATION */ bfrac *= r; return bfrac; } double bup(double *a,double *b,double *x,double *y,int *n,double *eps) /* ----------------------------------------------------------------------- EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER. EPS IS THE TOLERANCE USED. ----------------------------------------------------------------------- */ { static int K1 = 1; static int K2 = 0; static double bup,ap1,apb,d,l,r,t,w; static int i,k,kp1,mu,nm1; /* .. .. Executable Statements .. */ /* OBTAIN THE SCALING FACTOR EXP(-MU) AND EXP(MU)*(X**A*Y**B/BETA(A,B))/A */ apb = *a+*b; ap1 = *a+1.0e0; mu = 0; d = 1.0e0; if(*n == 1 || *a < 1.0e0) goto S10; if(apb < 1.1e0*ap1) goto S10; mu = (long)(fabs(exparg(&K1))); k = (long)(exparg(&K2)); if(k < mu) mu = k; t = mu; d = exp(-t); S10: bup = brcmp1(&mu,a,b,x,y)/ *a; if(*n == 1 || bup == 0.0e0) return bup; nm1 = *n-1; w = d; /* LET K BE THE INDEX OF THE MAXIMUM TERM */ k = 0; if(*b <= 1.0e0) goto S50; if(*y > 1.e-4) goto S20; k = nm1; goto S30; S20: r = (*b-1.0e0)**x/ *y-*a; if(r < 1.0e0) goto S50; t = nm1; k = (long)(t); if(r < t) k = (long)(r); S30: /* ADD THE INCREASING TERMS OF THE SERIES */ for(i=1; i<=k; i++) { l = i-1; d = (apb+l)/(ap1+l)**x*d; w += d; } if(k == nm1) goto S70; S50: /* ADD THE REMAINING TERMS OF THE SERIES */ kp1 = k+1; for(i=kp1; i<=nm1; i++) { l = i-1; d = (apb+l)/(ap1+l)**x*d; w += d; if(d <= *eps*w) goto S70; } S70: /* TERMINATE THE PROCEDURE */ bup *= w; return bup; } void bgrat(double *a,double *b,double *x,double *y,double *w,double *eps,int *ierr) /* ----------------------------------------------------------------------- ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B. THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED. IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. ----------------------------------------------------------------------- */ { static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z; static int i,n,nm1; static double c[30],d[30],T1; /* .. .. Executable Statements .. */ bm1 = *b-0.5e0-0.5e0; nu = *a+0.5e0*bm1; if(*y > 0.375e0) goto S10; T1 = -*y; lnx = alnrel(&T1); goto S20; S10: lnx = log(*x); S20: z = -(nu*lnx); if(*b*z == 0.0e0) goto S70; /* COMPUTATION OF THE EXPANSION SET R = EXP(-Z)*Z**B/GAMMA(B) */ r = *b*(1.0e0+gam1(b))*exp(*b*log(z)); r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx)); u = algdiv(b,a)+*b*log(nu); u = r*exp(-u); if(u == 0.0e0) goto S70; grat1(b,&z,&r,&p,&q,eps); v = 0.25e0*pow(1.0e0/nu,2.0); t2 = 0.25e0*lnx*lnx; l = *w/u; j = q/r; sum = j; t = cn = 1.0e0; n2 = 0.0e0; for(n=1; n<=30; n++) { bp2n = *b+n2; j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v; n2 += 2.0e0; t *= t2; cn /= (n2*(n2+1.0e0)); c[n-1] = cn; s = 0.0e0; if(n == 1) goto S40; nm1 = n-1; coef = *b-(double)n; for(i=1; i<=nm1; i++) { s += (coef*c[i-1]*d[n-i-1]); coef += *b; } S40: d[n-1] = bm1*cn+s/(double)n; dj = d[n-1]*j; sum += dj; if(sum <= 0.0e0) goto S70; if(fabs(dj) <= *eps*(sum+l)) goto S60; } S60: /* ADD THE RESULTS TO W */ *ierr = 0; *w += (u*sum); return; S70: /* THE EXPANSION CANNOT BE COMPUTED */ *ierr = 1; return; } void bratio(double *a,double *b,double *x,double *y,double *w,double *w1,int *ierr) /* ----------------------------------------------------------------------- EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B) -------------------- IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1 AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES W = IX(A,B) W1 = 1 - IX(A,B) IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED, THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO ONE OF THE FOLLOWING VALUES ... IERR = 1 IF A OR B IS NEGATIVE IERR = 2 IF A = B = 0 IERR = 3 IF X .LT. 0 OR X .GT. 1 IERR = 4 IF Y .LT. 0 OR Y .GT. 1 IERR = 5 IF X + Y .NE. 1 IERR = 6 IF X = A = 0 IERR = 7 IF Y = B = 0 -------------------- WRITTEN BY ALFRED H. MORRIS, JR. NAVAL SURFACE WARFARE CENTER DAHLGREN, VIRGINIA REVISED ... NOV 1991 ----------------------------------------------------------------------- */ { static int K1 = 1; static double a0,b0,eps,lambda,t,x0,y0,z; static int ierr1,ind,n; static double T2,T3,T4,T5; /* .. .. Executable Statements .. */ /* ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0 */ eps = spmpar(&K1); *w = *w1 = 0.0e0; if(*a < 0.0e0 || *b < 0.0e0) goto S270; if(*a == 0.0e0 && *b == 0.0e0) goto S280; if(*x < 0.0e0 || *x > 1.0e0) goto S290; if(*y < 0.0e0 || *y > 1.0e0) goto S300; z = *x+*y-0.5e0-0.5e0; if(fabs(z) > 3.0e0*eps) goto S310; *ierr = 0; if(*x == 0.0e0) goto S210; if(*y == 0.0e0) goto S230; if(*a == 0.0e0) goto S240; if(*b == 0.0e0) goto S220; eps = fifdmax1(eps,1.e-15); if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260; ind = 0; a0 = *a; b0 = *b; x0 = *x; y0 = *y; if(fifdmin1(a0,b0) > 1.0e0) goto S40; /* PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1 */ if(*x <= 0.5e0) goto S10; ind = 1; a0 = *b; b0 = *a; x0 = *y; y0 = *x; S10: if(b0 < fifdmin1(eps,eps*a0)) goto S90; if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100; if(fifdmax1(a0,b0) > 1.0e0) goto S20; if(a0 >= fifdmin1(0.2e0,b0)) goto S110; if(pow(x0,a0) <= 0.9e0) goto S110; if(x0 >= 0.3e0) goto S120; n = 20; goto S140; S20: if(b0 <= 1.0e0) goto S110; if(x0 >= 0.3e0) goto S120; if(x0 >= 0.1e0) goto S30; if(pow(x0*b0,a0) <= 0.7e0) goto S110; S30: if(b0 > 15.0e0) goto S150; n = 20; goto S140; S40: /* PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1 */ if(*a > *b) goto S50; lambda = *a-(*a+*b)**x; goto S60; S50: lambda = (*a+*b)**y-*b; S60: if(lambda >= 0.0e0) goto S70; ind = 1; a0 = *b; b0 = *a; x0 = *y; y0 = *x; lambda = fabs(lambda); S70: if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110; if(b0 < 40.0e0) goto S160; if(a0 > b0) goto S80; if(a0 <= 100.0e0) goto S130; if(lambda > 0.03e0*a0) goto S130; goto S200; S80: if(b0 <= 100.0e0) goto S130; if(lambda > 0.03e0*b0) goto S130; goto S200; S90: /* EVALUATION OF THE APPROPRIATE ALGORITHM */ *w = fpser(&a0,&b0,&x0,&eps); *w1 = 0.5e0+(0.5e0-*w); goto S250; S100: *w1 = apser(&a0,&b0,&x0,&eps); *w = 0.5e0+(0.5e0-*w1); goto S250; S110: *w = bpser(&a0,&b0,&x0,&eps); *w1 = 0.5e0+(0.5e0-*w); goto S250; S120: *w1 = bpser(&b0,&a0,&y0,&eps); *w = 0.5e0+(0.5e0-*w1); goto S250; S130: T2 = 15.0e0*eps; *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2); *w1 = 0.5e0+(0.5e0-*w); goto S250; S140: *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps); b0 += (double)n; S150: T3 = 15.0e0*eps; bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1); *w = 0.5e0+(0.5e0-*w1); goto S250; S160: n = (long)(b0); b0 -= (double)n; if(b0 != 0.0e0) goto S170; n -= 1; b0 = 1.0e0; S170: *w = bup(&b0,&a0,&y0,&x0,&n,&eps); if(x0 > 0.7e0) goto S180; *w += bpser(&a0,&b0,&x0,&eps); *w1 = 0.5e0+(0.5e0-*w); goto S250; S180: if(a0 > 15.0e0) goto S190; n = 20; *w += bup(&a0,&b0,&x0,&y0,&n,&eps); a0 += (double)n; S190: T4 = 15.0e0*eps; bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1); *w1 = 0.5e0+(0.5e0-*w); goto S250; S200: T5 = 100.0e0*eps; *w = basym(&a0,&b0,&lambda,&T5); *w1 = 0.5e0+(0.5e0-*w); goto S250; S210: /* TERMINATION OF THE PROCEDURE */ if(*a == 0.0e0) goto S320; S220: *w = 0.0e0; *w1 = 1.0e0; return; S230: if(*b == 0.0e0) goto S330; S240: *w = 1.0e0; *w1 = 0.0e0; return; S250: if(ind == 0) return; t = *w; *w = *w1; *w1 = t; return; S260: /* PROCEDURE FOR A AND B .LT. 1.E-3*EPS */ *w = *b/(*a+*b); *w1 = *a/(*a+*b); return; S270: /* ERROR RETURN */ *ierr = 1; return; S280: *ierr = 2; return; S290: *ierr = 3; return; S300: *ierr = 4; return; S310: *ierr = 5; return; S320: *ierr = 6; return; S330: *ierr = 7; return; } double exparg(int *l) /* -------------------------------------------------------------------- IF L = 0 THEN EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH EXP(W) CAN BE COMPUTED. IF L IS NONZERO THEN EXPARG(L) = THE LARGEST NEGATIVE W FOR WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO. NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED. -------------------------------------------------------------------- */ { static int K1 = 4; static int K2 = 9; static int K3 = 10; static double exparg,lnb; static int b,m; /* .. .. Executable Statements .. */ b = ipmpar(&K1); if(b != 2) goto S10; lnb = .69314718055995e0; goto S40; S10: if(b != 8) goto S20; lnb = 2.0794415416798e0; goto S40; S20: if(b != 16) goto S30; lnb = 2.7725887222398e0; goto S40; S30: lnb = log((double)b); S40: if(*l == 0) goto S50; m = ipmpar(&K2)-1; exparg = 0.99999e0*((double)m*lnb); return exparg; S50: m = ipmpar(&K3); exparg = 0.99999e0*((double)m*lnb); return exparg; } void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum) /* ********************************************************************** void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum) CUMulative F distribution Function Computes the integral from 0 to F of the f-density with DFN and DFD degrees of freedom. Arguments F --> Upper limit of integration of the f-density. F is DOUBLE PRECISION DFN --> Degrees of freedom of the numerator sum of squares. DFN is DOUBLE PRECISI DFD --> Degrees of freedom of the denominator sum of squares. DFD is DOUBLE PRECISI CUM <-- Cumulative f distribution. CUM is DOUBLE PRECISI CCUM <-- Compliment of Cumulative f distribution. CCUM is DOUBLE PRECIS Method Formula 26.5.28 of Abramowitz and Stegun is used to reduce the cumulative F to a cumulative beta distribution. Note If F is less than or equal to 0, 0 is returned. ********************************************************************** */ { #define half 0.5e0 #define done 1.0e0 static double dsum,prod,xx,yy; static int ierr; static double T1,T2; /* .. .. Executable Statements .. */ if(!(*f <= 0.0e0)) goto S10; *cum = 0.0e0; *ccum = 1.0e0; return; S10: prod = *dfn**f; /* XX is such that the incomplete beta with parameters DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM YY is 1 - XX Calculate the smaller of XX and YY accurately */ dsum = *dfd+prod; xx = *dfd/dsum; if(xx > half) { yy = prod/dsum; xx = done-yy; } else yy = done-xx; T1 = *dfd*half; T2 = *dfn*half; bratio(&T1,&T2,&xx,&yy,ccum,cum,&ierr); return; #undef half #undef done } void cdff(int *which,double *p,double *q,double *f,double *dfn,double *dfd,int *status,double *bound) /********************************************************************** void cdff(int *which,double *p,double *q,double *f,double *dfn, double *dfd,int *status,double *bound) Cumulative Distribution Function F distribution Function Calculates any one parameter of the F distribution given values for the others. Arguments WHICH --> Integer indicating which of the next four argument values is to be calculated from the others. Legal range: 1..4 iwhich = 1 : Calculate P and Q from F,DFN and DFD iwhich = 2 : Calculate F from P,Q,DFN and DFD iwhich = 3 : Calculate DFN from P,Q,F and DFD iwhich = 4 : Calculate DFD from P,Q,F and DFN P <--> The integral from 0 to F of the f-density. Input range: [0,1]. Q <--> 1-P. Input range: (0, 1]. P + Q = 1.0. F <--> Upper limit of integration of the f-density. Input range: [0, +infinity). Search range: [0,1E100] DFN < --> Degrees of freedom of the numerator sum of squares. Input range: (0, +infinity). Search range: [ 1E-100, 1E100] DFD < --> Degrees of freedom of the denominator sum of squares. Input range: (0, +infinity). Search range: [ 1E-100, 1E100] STATUS <-- 0 if calculation completed correctly -I if input parameter number I is out of range 1 if answer appears to be lower than lowest search bound 2 if answer appears to be higher than greatest search bound 3 if P + Q .ne. 1 BOUND <-- Undefined if STATUS is 0 Bound exceeded by parameter number I if STATUS is negative. Lower search bound if STATUS is 1. Upper search bound if STATUS is 2. Method Formula 26.6.2 of Abramowitz and Stegun, Handbook of Mathematical Functions (1966) is used to reduce the computation of the cumulative distribution function for the F variate to that of an incomplete beta. Computation of other parameters involve a seach for a value that produces the desired value of P. The search relies on the monotinicity of P with the other parameter. WARNING The value of the cumulative F distribution is not necessarily monotone in either degrees of freedom. There thus may be two values that provide a given CDF value. This routine assumes monotonicity and will find an arbitrary one of the two values. **********************************************************************/ { #define tol 1.0e-8 #define atol 1.0e-50 #define zero 1.0e-100 #define inf 1.0e100 static int K1 = 1; static double K2 = 0.0e0; static double K4 = 0.5e0; static double K5 = 5.0e0; static double pq,fx,cum,ccum; static unsigned long qhi,qleft,qporq; static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15; /* .. .. Executable Statements .. */ /* Check arguments */ if(!(*which < 1 || *which > 4)) goto S30; if(!(*which < 1)) goto S10; *bound = 1.0e0; goto S20; S10: *bound = 4.0e0; S20: *status = -1; return; S30: if(*which == 1) goto S70; /* P */ if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60; if(!(*p < 0.0e0)) goto S40; *bound = 0.0e0; goto S50; S40: *bound = 1.0e0; S50: *status = -2; return; S70: S60: if(*which == 1) goto S110; /* Q */ if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100; if(!(*q <= 0.0e0)) goto S80; *bound = 0.0e0; goto S90; S80: *bound = 1.0e0; S90: *status = -3; return; S110: S100: if(*which == 2) goto S130; /* F */ if(!(*f < 0.0e0)) goto S120; *bound = 0.0e0; *status = -4; return; S130: S120: if(*which == 3) goto S150; /* DFN */ if(!(*dfn <= 0.0e0)) goto S140; *bound = 0.0e0; *status = -5; return; S150: S140: if(*which == 4) goto S170; /* DFD */ if(!(*dfd <= 0.0e0)) goto S160; *bound = 0.0e0; *status = -6; return; S170: S160: if(*which == 1) goto S210; /* P + Q */ pq = *p+*q; if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200; if(!(pq < 0.0e0)) goto S180; *bound = 0.0e0; goto S190; S180: *bound = 1.0e0; S190: *status = 3; return; S210: S200: if(!(*which == 1)) qporq = *p <= *q; /* Select the minimum of P or Q Calculate ANSWERS */ if(1 == *which) { /* Calculating P */ cumf(f,dfn,dfd,p,q); *status = 0; } else if(2 == *which) { /* Calculating F */ *f = 5.0e0; T3 = inf; T6 = atol; T7 = tol; dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7); *status = 0; dinvr(status,f,&fx,&qleft,&qhi); S220: if(!(*status == 1)) goto S250; cumf(f,dfn,dfd,&cum,&ccum); if(!qporq) goto S230; fx = cum-*p; goto S240; S230: fx = ccum-*q; S240: dinvr(status,f,&fx,&qleft,&qhi); goto S220; S250: if(!(*status == -1)) goto S280; if(!qleft) goto S260; *status = 1; *bound = 0.0e0; goto S270; S260: *status = 2; *bound = inf; S280: S270: ; } else if(3 == *which) { /* Calculating DFN */ *dfn = 5.0e0; T8 = zero; T9 = inf; T10 = atol; T11 = tol; dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11); *status = 0; dinvr(status,dfn,&fx,&qleft,&qhi); S290: if(!(*status == 1)) goto S320; cumf(f,dfn,dfd,&cum,&ccum); if(!qporq) goto S300; fx = cum-*p; goto S310; S300: fx = ccum-*q; S310: dinvr(status,dfn,&fx,&qleft,&qhi); goto S290; S320: if(!(*status == -1)) goto S350; if(!qleft) goto S330; *status = 1; *bound = zero; goto S340; S330: *status = 2; *bound = inf; S350: S340: ; } else if(4 == *which) { /* Calculating DFD */ *dfd = 5.0e0; T12 = zero; T13 = inf; T14 = atol; T15 = tol; dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15); *status = 0; dinvr(status,dfd,&fx,&qleft,&qhi); S360: if(!(*status == 1)) goto S390; cumf(f,dfn,dfd,&cum,&ccum); if(!qporq) goto S370; fx = cum-*p; goto S380; S370: fx = ccum-*q; S380: dinvr(status,dfd,&fx,&qleft,&qhi); goto S360; S390: if(!(*status == -1)) goto S420; if(!qleft) goto S400; *status = 1; *bound = zero; goto S410; S400: *status = 2; *bound = inf; S410: ; } S420: return; #undef tol #undef atol #undef zero #undef inf } /* DEFINE DZROR */ static void E0001(int IENTRY,int *status,double *x,double *fx,double *xlo,double *xhi,unsigned long *qleft,unsigned long *qhi,double *zabstl,double *zreltl,double *zxhi,double *zxlo) { #define ftol(zx) (0.5e0*fifdmax1(abstol,reltol*fabs((zx)))) static double a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q,reltol,tol,w,xxhi,xxlo; static int ext,i99999; static unsigned long first,qrzero; switch(IENTRY){case 0: goto DZROR; case 1: goto DSTZR;} DZROR: if(*status > 0) goto S280; *xlo = xxlo; *xhi = xxhi; b = *x = *xlo; /* GET-FUNCTION-VALUE */ i99999 = 1; goto S270; S10: fb = *fx; *xlo = *xhi; a = *x = *xlo; /* GET-FUNCTION-VALUE */ i99999 = 2; goto S270; S20: /* Check that F(ZXLO) < 0 < F(ZXHI) or F(ZXLO) > 0 > F(ZXHI) */ if(!(fb < 0.0e0)) goto S40; if(!(*fx < 0.0e0)) goto S30; *status = -1; *qleft = *fx < fb; *qhi = 0; return; S40: S30: if(!(fb > 0.0e0)) goto S60; if(!(*fx > 0.0e0)) goto S50; *status = -1; *qleft = *fx > fb; *qhi = 1; return; S60: S50: fa = *fx; first = 1; S70: c = a; fc = fa; ext = 0; S80: if(!(fabs(fc) < fabs(fb))) goto S100; if(!(c != a)) goto S90; d = a; fd = fa; S90: a = b; fa = fb; *xlo = c; b = *xlo; fb = fc; c = a; fc = fa; S100: tol = ftol(*xlo); m = (c+b)*.5e0; mb = m-b; if(!(fabs(mb) > tol)) goto S240; if(!(ext > 3)) goto S110; w = mb; goto S190; S110: tol = fifdsign(tol,mb); p = (b-a)*fb; if(!first) goto S120; q = fa-fb; first = 0; goto S130; S120: fdb = (fd-fb)/(d-b); fda = (fd-fa)/(d-a); p = fda*p; q = fdb*fa-fda*fb; S130: if(!(p < 0.0e0)) goto S140; p = -p; q = -q; S140: if(ext == 3) p *= 2.0e0; if(!(p*1.0e0 == 0.0e0 || p <= q*tol)) goto S150; w = tol; goto S180; S150: if(!(p < mb*q)) goto S160; w = p/q; goto S170; S160: w = mb; S190: S180: S170: d = a; fd = fa; a = b; fa = fb; b += w; *xlo = b; *x = *xlo; /* GET-FUNCTION-VALUE */ i99999 = 3; goto S270; S200: fb = *fx; if(!(fc*fb >= 0.0e0)) goto S210; goto S70; S210: if(!(w == mb)) goto S220; ext = 0; goto S230; S220: ext += 1; S230: goto S80; S240: *xhi = c; qrzero = (fc >= 0.0e0 && fb <= 0.0e0) || (fc < 0.0e0 && fb >= 0.0e0); if(!qrzero) goto S250; *status = 0; goto S260; S250: *status = -1; S260: return; DSTZR: xxlo = *zxlo; xxhi = *zxhi; abstol = *zabstl; reltol = *zreltl; return; S270: /* TO GET-FUNCTION-VALUE */ *status = 1; return; S280: switch((int)i99999) { case 1: goto S10;case 2: goto S20;case 3: goto S200; default: break; } #undef ftol } /* DEFINE DINVR */ static void E0000(int IENTRY,int *status,double *x,double *fx,unsigned long *qleft,unsigned long *qhi,double *zabsst,double *zabsto,double *zbig,double *zrelst,double *zrelto,double *zsmall,double *zstpmu) { #define qxmon(zx,zy,zz) (int)((zx) <= (zy) && (zy) <= (zz)) static double absstp,abstol,big,fbig,fsmall,relstp,reltol,petit,step,stpmul,xhi,xlb,xlo,xsave,xub,yy; static int i99999; static unsigned long qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup; switch(IENTRY){case 0: goto DINVR; case 1: goto DSTINV;} DINVR: if(*status > 0) goto S310; qcond = !qxmon(petit,*x,big); if(qcond) ftnstop(" SMALL, X, BIG not monotone in INVR"); xsave = *x; /* See that SMALL and BIG bound the zero and set QINCR */ *x = petit; /* GET-FUNCTION-VALUE */ i99999 = 1; goto S300; S10: fsmall = *fx; *x = big; /* GET-FUNCTION-VALUE */ i99999 = 2; goto S300; S20: fbig = *fx; qincr = fbig > fsmall; if(!qincr) goto S50; if(fsmall <= 0.0e0) goto S30; *status = -1; *qleft = *qhi = 1; return; S30: if(fbig >= 0.0e0) goto S40; *status = -1; *qleft = *qhi = 0; return; S40: goto S80; S50: if(fsmall >= 0.0e0) goto S60; *status = -1; *qleft = 1; *qhi = 0; return; S60: if(fbig <= 0.0e0) goto S70; *status = -1; *qleft = 0; *qhi = 1; return; S80: S70: *x = xsave; step = fifdmax1(absstp,relstp*fabs(*x)); /* YY = F(X) - Y GET-FUNCTION-VALUE */ i99999 = 3; goto S300; S90: yy = *fx; if(!(yy == 0.0e0)) goto S100; *status = 0; qok = 1; return; S100: qup = (qincr && yy < 0.0e0) || (!qincr && yy > 0.0e0); //## modif /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ HANDLE CASE IN WHICH WE MUST STEP HIGHER ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ if(!qup) goto S170; xlb = xsave; xub = fifdmin1(xlb+step,big); goto S120; S110: if(qcond) goto S150; S120: /* YY = F(XUB) - Y */ *x = xub; /* GET-FUNCTION-VALUE */ i99999 = 4; goto S300; S130: yy = *fx; qbdd = (qincr && yy >= 0.0e0) || (!qincr && yy <= 0.0e0) ; qlim = xub >= big; qcond = qbdd || qlim; if(qcond) goto S140; step = stpmul*step; xlb = xub; xub = fifdmin1(xlb+step,big); S140: goto S110; S150: if(!(qlim && !qbdd)) goto S160; *status = -1; *qleft = 0; *qhi = !qincr; *x = big; return; S160: goto S240; S170: /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ HANDLE CASE IN WHICH WE MUST STEP LOWER ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ xub = xsave; xlb = fifdmax1(xub-step,petit); goto S190; S180: if(qcond) goto S220; S190: /* YY = F(XLB) - Y */ *x = xlb; /* GET-FUNCTION-VALUE */ i99999 = 5; goto S300; S200: yy = *fx; qbdd = (qincr && yy <= 0.0e0) || (!qincr && yy >= 0.0e0); qlim = xlb <= petit; qcond = qbdd || qlim; if(qcond) goto S210; step = stpmul*step; xub = xlb; xlb = fifdmax1(xub-step,petit); S210: goto S180; S220: if(!(qlim && !qbdd)) goto S230; *status = -1; *qleft = 1; *qhi = qincr; *x = petit; return; S240: S230: dstzr(&xlb,&xub,&abstol,&reltol); /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ *status = 0; goto S260; S250: if(!(*status == 1)) goto S290; S260: dzror(status,x,fx,&xlo,&xhi,&qdum1,&qdum2); if(!(*status == 1)) goto S280; /* GET-FUNCTION-VALUE */ i99999 = 6; goto S300; S280: S270: goto S250; S290: *x = xlo; *status = 0; return; DSTINV: petit = *zsmall; big = *zbig; absstp = *zabsst; relstp = *zrelst; stpmul = *zstpmu; abstol = *zabsto; reltol = *zrelto; return; S300: /* TO GET-FUNCTION-VALUE */ *status = 1; return; S310: switch((int)i99999) { case 1: goto S10;case 2: goto S20;case 3: goto S90;case 4: goto S130;case 5: goto S200;case 6: goto S270;default: break; } #undef qxmon } double gam1(double *a) /* ------------------------------------------------------------------ COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5 ------------------------------------------------------------------ */ { static double s1 = .273076135303957e+00; static double s2 = .559398236957378e-01; static double p[7] = { .577215664901533e+00,-.409078193005776e+00,-.230975380857675e+00, .597275330452234e-01,.766968181649490e-02,-.514889771323592e-02, .589597428611429e-03 }; static double q[5] = { .100000000000000e+01,.427569613095214e+00,.158451672430138e+00, .261132021441447e-01,.423244297896961e-02 }; static double r[9] = { -.422784335098468e+00,-.771330383816272e+00,-.244757765222226e+00, .118378989872749e+00,.930357293360349e-03,-.118290993445146e-01, .223047661158249e-02,.266505979058923e-03,-.132674909766242e-03 }; static double gam1,bot,d,t,top,w,T1; /* .. .. Executable Statements .. */ t = *a; d = *a-0.5e0; if(d > 0.0e0) t = d-0.5e0; T1 = t; if(T1 < 0) goto S40; else if(T1 == 0) goto S10; else goto S20; S10: gam1 = 0.0e0; return gam1; S20: top = (((((p[6]*t+p[5])*t+p[4])*t+p[3])*t+p[2])*t+p[1])*t+p[0]; bot = (((q[4]*t+q[3])*t+q[2])*t+q[1])*t+1.0e0; w = top/bot; if(d > 0.0e0) goto S30; gam1 = *a*w; return gam1; S30: gam1 = t/ *a*(w-0.5e0-0.5e0); return gam1; S40: top = (((((((r[8]*t+r[7])*t+r[6])*t+r[5])*t+r[4])*t+r[3])*t+r[2])*t+r[1])*t+ r[0]; bot = (s2*t+s1)*t+1.0e0; w = top/bot; if(d > 0.0e0) goto S50; gam1 = *a*(w+0.5e0+0.5e0); return gam1; S50: gam1 = t*w/ *a; return gam1; } void dinvr(int *status,double *x,double *fx,unsigned long *qleft,unsigned long *qhi) /* ********************************************************************** void dinvr(int *status,double *x,double *fx, unsigned long *qleft,unsigned long *qhi) Double precision bounds the zero of the function and invokes zror Reverse Communication Function Bounds the function and invokes ZROR to perform the zero finding. STINVR must have been called before this routine in order to set its parameters. Arguments STATUS <--> At the beginning of a zero finding problem, STATUS should be set to 0 and INVR invoked. (The value of parameters other than X will be ignored on this cal When INVR needs the function evaluated, it will set STATUS to 1 and return. The value of the function should be set in FX and INVR again called without changing any of its other parameters. When INVR has finished without error, it will return with STATUS 0. In that case X is approximately a root of F(X). If INVR cannot bound the function, it returns status -1 and sets QLEFT and QHI. INTEGER STATUS X <-- The value of X at which F(X) is to be evaluated. DOUBLE PRECISION X FX --> The value of F(X) calculated when INVR returns with STATUS = 1. DOUBLE PRECISION FX QLEFT <-- Defined only if QMFINV returns .FALSE. In that case it is .TRUE. If the stepping search terminated unsucessfully at SMALL. If it is .FALSE. the search terminated unsucessfully at BIG. QLEFT is LOGICAL QHI <-- Defined only if QMFINV returns .FALSE. In that case it is .TRUE. if F(X) .GT. Y at the termination of the search and .FALSE. if F(X) .LT. Y at the termination of the search. QHI is LOGICAL ********************************************************************** */ { E0000(0,status,x,fx,qleft,qhi,NULL,NULL,NULL,NULL,NULL,NULL,NULL); } void dstinv(double *zsmall,double *zbig,double *zabsst,double *zrelst,double *zstpmu,double *zabsto,double *zrelto) /* ********************************************************************** void dstinv(double *zsmall,double *zbig,double *zabsst, double *zrelst,double *zstpmu,double *zabsto, double *zrelto) Double Precision - SeT INverse finder - Reverse Communication Function Concise Description - Given a monotone function F finds X such that F(X) = Y. Uses Reverse communication -- see invr. This routine sets quantities needed by INVR. More Precise Description of INVR - F must be a monotone function, the results of QMFINV are otherwise undefined. QINCR must be .TRUE. if F is non- decreasing and .FALSE. if F is non-increasing. QMFINV will return .TRUE. if and only if F(SMALL) and F(BIG) bracket Y, i. e., QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL) if QMFINV returns .TRUE., then the X returned satisfies the following condition. let TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) then if QINCR is .TRUE., F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X)) and if QINCR is .FALSE. F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X)) Arguments SMALL --> The left endpoint of the interval to be searched for a solution. SMALL is DOUBLE PRECISION BIG --> The right endpoint of the interval to be searched for a solution. BIG is DOUBLE PRECISION ABSSTP, RELSTP --> The initial step size in the search is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm. ABSSTP is DOUBLE PRECISION RELSTP is DOUBLE PRECISION STPMUL --> When a step doesn't bound the zero, the step size is multiplied by STPMUL and another step taken. A popular value is 2.0 DOUBLE PRECISION STPMUL ABSTOL, RELTOL --> Two numbers that determine the accuracy of the solution. See function for a precise definition. ABSTOL is DOUBLE PRECISION RELTOL is DOUBLE PRECISION Method Compares F(X) with Y for the input value of X then uses QINCR to determine whether to step left or right to bound the desired x. the initial step size is MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X. Iteratively steps right or left until it bounds X. At each step which doesn't bound X, the step size is doubled. The routine is careful never to step beyond SMALL or BIG. If it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE. after setting QLEFT and QHI. If X is successfully bounded then Algorithm R of the paper 'Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function' by J. C. P. Bus and T. J. Dekker in ACM Transactions on Mathematical Software, Volume 1, No. 4 page 330 (DEC. '75) is employed to find the zero of the function F(X)-Y. This is routine QRZERO. ********************************************************************** */ { E0000(1,NULL,NULL,NULL,NULL,NULL,zabsst,zabsto,zbig,zrelst,zrelto,zsmall, zstpmu); } void dzror(int *status,double *x,double *fx,double *xlo,double *xhi,unsigned long *qleft,unsigned long *qhi) /* ********************************************************************** void dzror(int *status,double *x,double *fx,double *xlo, double *xhi,unsigned long *qleft,unsigned long *qhi) Double precision ZeRo of a function -- Reverse Communication Function Performs the zero finding. STZROR must have been called before this routine in order to set its parameters. Arguments STATUS <--> At the beginning of a zero finding problem, STATUS should be set to 0 and ZROR invoked. (The value of other parameters will be ignored on this call.) When ZROR needs the function evaluated, it will set STATUS to 1 and return. The value of the function should be set in FX and ZROR again called without changing any of its other parameters. When ZROR has finished without error, it will return with STATUS 0. In that case (XLO,XHI) bound the answe If ZROR finds an error (which implies that F(XLO)-Y an F(XHI)-Y have the same sign, it returns STATUS -1. In this case, XLO and XHI are undefined. INTEGER STATUS X <-- The value of X at which F(X) is to be evaluated. DOUBLE PRECISION X FX --> The value of F(X) calculated when ZROR returns with STATUS = 1. DOUBLE PRECISION FX XLO <-- When ZROR returns with STATUS = 0, XLO bounds the inverval in X containing the solution below. DOUBLE PRECISION XLO XHI <-- When ZROR returns with STATUS = 0, XHI bounds the inverval in X containing the solution above. DOUBLE PRECISION XHI QLEFT <-- .TRUE. if the stepping search terminated unsucessfully at XLO. If it is .FALSE. the search terminated unsucessfully at XHI. QLEFT is LOGICAL QHI <-- .TRUE. if F(X) .GT. Y at the termination of the search and .FALSE. if F(X) .LT. Y at the termination of the search. QHI is LOGICAL ********************************************************************** */ { E0001(0,status,x,fx,xlo,xhi,qleft,qhi,NULL,NULL,NULL,NULL); } void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl) /* ********************************************************************** void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl) Double precision SeT ZeRo finder - Reverse communication version Function Sets quantities needed by ZROR. The function of ZROR and the quantities set is given here. Concise Description - Given a function F find XLO such that F(XLO) = 0. More Precise Description - Input condition. F is a double precision function of a single double precision argument and XLO and XHI are such that F(XLO)*F(XHI) .LE. 0.0 If the input condition is met, QRZERO returns .TRUE. and output values of XLO and XHI satisfy the following F(XLO)*F(XHI) .LE. 0. ABS(F(XLO) .LE. ABS(F(XHI) ABS(XLO-XHI) .LE. TOL(X) where TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) If this algorithm does not find XLO and XHI satisfying these conditions then QRZERO returns .FALSE. This implies that the input condition was not met. Arguments XLO --> The left endpoint of the interval to be searched for a solution. XLO is DOUBLE PRECISION XHI --> The right endpoint of the interval to be for a solution. XHI is DOUBLE PRECISION ABSTOL, RELTOL --> Two numbers that determine the accuracy of the solution. See function for a precise definition. ABSTOL is DOUBLE PRECISION RELTOL is DOUBLE PRECISION Method Algorithm R of the paper 'Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function' by J. C. P. Bus and T. J. Dekker in ACM Transactions on Mathematical Software, Volume 1, no. 4 page 330 (Dec. '75) is employed to find the zero of F(X)-Y. ********************************************************************** */ { E0001(1,NULL,NULL,NULL,NULL,NULL,NULL,NULL,zabstl,zreltl,zxhi,zxlo); } double erf1(double *x) /* ----------------------------------------------------------------------- EVALUATION OF THE REAL ERROR FUNCTION ----------------------------------------------------------------------- */ { static double c = .564189583547756e0; static double a[5] = { .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01, .479137145607681e-01,.128379167095513e+00 }; static double b[3] = { .301048631703895e-02,.538971687740286e-01,.375795757275549e+00 }; static double p[8] = { -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00, 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02, 4.51918953711873e+02,3.00459261020162e+02 }; static double q[8] = { 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01, 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02, 7.90950925327898e+02,3.00459260956983e+02 }; static double r[5] = { 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01, 4.65807828718470e+00,2.82094791773523e-01 }; static double s[4] = { 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01, 1.80124575948747e+01 }; static double erf1,ax,bot,t,top,x2; /* .. .. Executable Statements .. */ ax = fabs(*x); if(ax > 0.5e0) goto S10; t = *x**x; top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0; bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0; erf1 = *x*(top/bot); return erf1; S10: if(ax > 4.0e0) goto S20; top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[ 7]; bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[ 7]; erf1 = 0.5e0+(0.5e0-exp(-(*x**x))*top/bot); if(*x < 0.0e0) erf1 = -erf1; return erf1; S20: if(ax >= 5.8e0) goto S30; x2 = *x**x; t = 1.0e0/x2; top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4]; bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0; erf1 = (c-top/(x2*bot))/ax; erf1 = 0.5e0+(0.5e0-exp(-x2)*erf1); if(*x < 0.0e0) erf1 = -erf1; return erf1; S30: erf1 = fifdsign(1.0e0,*x); return erf1; } double erfc1(int *ind,double *x) /* ----------------------------------------------------------------------- EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION ERFC1(IND,X) = ERFC(X) IF IND = 0 ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE ----------------------------------------------------------------------- */ { static double c = .564189583547756e0; static double a[5] = { .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01, .479137145607681e-01,.128379167095513e+00 }; static double b[3] = { .301048631703895e-02,.538971687740286e-01,.375795757275549e+00 }; static double p[8] = { -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00, 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02, 4.51918953711873e+02,3.00459261020162e+02 }; static double q[8] = { 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01, 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02, 7.90950925327898e+02,3.00459260956983e+02 }; static double r[5] = { 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01, 4.65807828718470e+00,2.82094791773523e-01 }; static double s[4] = { 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01, 1.80124575948747e+01 }; static int K1 = 1; static double erfc1,ax,bot,e,t,top,w; /* .. .. Executable Statements .. */ /* ABS(X) .LE. 0.5 */ ax = fabs(*x); if(ax > 0.5e0) goto S10; t = *x**x; top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0; bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0; erfc1 = 0.5e0+(0.5e0-*x*(top/bot)); if(*ind != 0) erfc1 = exp(t)*erfc1; return erfc1; S10: /* 0.5 .LT. ABS(X) .LE. 4 */ if(ax > 4.0e0) goto S20; top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[ 7]; bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[ 7]; erfc1 = top/bot; goto S40; S20: /* ABS(X) .GT. 4 */ if(*x <= -5.6e0) goto S60; if(*ind != 0) goto S30; if(*x > 100.0e0) goto S70; if(*x**x > -exparg(&K1)) goto S70; S30: t = pow(1.0e0/ *x,2.0); top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4]; bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0; erfc1 = (c-t*top/bot)/ax; S40: /* FINAL ASSEMBLY */ if(*ind == 0) goto S50; if(*x < 0.0e0) erfc1 = 2.0e0*exp(*x**x)-erfc1; return erfc1; S50: w = *x**x; t = w; e = w-t; erfc1 = (0.5e0+(0.5e0-e))*exp(-t)*erfc1; if(*x < 0.0e0) erfc1 = 2.0e0-erfc1; return erfc1; S60: /* LIMIT VALUE FOR LARGE NEGATIVE X */ erfc1 = 2.0e0; if(*ind != 0) erfc1 = 2.0e0*exp(*x**x); return erfc1; S70: /* LIMIT VALUE FOR LARGE POSITIVE X WHEN IND = 0 */ erfc1 = 0.0e0; return erfc1; } double esum(int *mu,double *x) /* ----------------------------------------------------------------------- EVALUATION OF EXP(MU + X) ----------------------------------------------------------------------- */ { static double esum,w; /* .. .. Executable Statements .. */ if(*x > 0.0e0) goto S10; if(*mu < 0) goto S20; w = (double)*mu+*x; if(w > 0.0e0) goto S20; esum = exp(w); return esum; S10: if(*mu > 0) goto S20; w = (double)*mu+*x; if(w < 0.0e0) goto S20; esum = exp(w); return esum; S20: w = *mu; esum = exp(w)*exp(*x); return esum; } double fpser(double *a,double *b,double *x,double *eps) /* ----------------------------------------------------------------------- EVALUATION OF I (A,B) X FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5. ----------------------------------------------------------------------- SET FPSER = X**A */ { static int K1 = 1; static double fpser,an,c,s,t,tol; /* .. .. Executable Statements .. */ fpser = 1.0e0; if(*a <= 1.e-3**eps) goto S10; fpser = 0.0e0; t = *a*log(*x); if(t < exparg(&K1)) return fpser; fpser = exp(t); S10: /* NOTE THAT 1/B(A,B) = B */ fpser = *b/ *a*fpser; tol = *eps/ *a; an = *a+1.0e0; t = *x; s = t/an; S20: an += 1.0e0; t = *x*t; c = t/an; s += c; /* printf("%i\n",__LINE__); printf("tol = %f\n",tol); printf("c = %f\n",c); printf("fabs(c) = %f\n",fabs(c)); */ if(fabs(c) > tol) goto S20; fpser *= (1.0e0+*a*s); return fpser; } double gamln(double *a) /* ----------------------------------------------------------------------- EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A ----------------------------------------------------------------------- WRITTEN BY ALFRED H. MORRIS NAVAL SURFACE WARFARE CENTER DAHLGREN, VIRGINIA -------------------------- D = 0.5*(LN(2*PI) - 1) -------------------------- */ { static double c0 = .833333333333333e-01; static double c1 = -.277777777760991e-02; static double c2 = .793650666825390e-03; static double c3 = -.595202931351870e-03; static double c4 = .837308034031215e-03; static double c5 = -.165322962780713e-02; static double d = .418938533204673e0; static double gamln,t,w; static int i,n; static double T1; /* .. .. Executable Statements .. */ if(*a > 0.8e0) goto S10; gamln = gamln1(a)-log(*a); return gamln; S10: if(*a > 2.25e0) goto S20; t = *a-0.5e0-0.5e0; gamln = gamln1(&t); return gamln; S20: if(*a >= 10.0e0) goto S40; n = (long)(*a - 1.25e0); t = *a; w = 1.0e0; for(i=1; i<=n; i++) { t -= 1.0e0; w = t*w; } T1 = t-1.0e0; gamln = gamln1(&T1)+log(w); return gamln; S40: t = pow(1.0e0/ *a,2.0); w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a; gamln = d+w+(*a-0.5e0)*(log(*a)-1.0e0); return gamln; } double gamln1(double *a) /* ----------------------------------------------------------------------- EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25 ----------------------------------------------------------------------- */ { static double p0 = .577215664901533e+00; static double p1 = .844203922187225e+00; static double p2 = -.168860593646662e+00; static double p3 = -.780427615533591e+00; static double p4 = -.402055799310489e+00; static double p5 = -.673562214325671e-01; static double p6 = -.271935708322958e-02; static double q1 = .288743195473681e+01; static double q2 = .312755088914843e+01; static double q3 = .156875193295039e+01; static double q4 = .361951990101499e+00; static double q5 = .325038868253937e-01; static double q6 = .667465618796164e-03; static double r0 = .422784335098467e+00; static double r1 = .848044614534529e+00; static double r2 = .565221050691933e+00; static double r3 = .156513060486551e+00; static double r4 = .170502484022650e-01; static double r5 = .497958207639485e-03; static double s1 = .124313399877507e+01; static double s2 = .548042109832463e+00; static double s3 = .101552187439830e+00; static double s4 = .713309612391000e-02; static double s5 = .116165475989616e-03; static double gamln1,w,x; /* .. .. Executable Statements .. */ if(*a >= 0.6e0) goto S10; w = ((((((p6**a+p5)**a+p4)**a+p3)**a+p2)**a+p1)**a+p0)/((((((q6**a+q5)**a+ q4)**a+q3)**a+q2)**a+q1)**a+1.0e0); gamln1 = -(*a*w); return gamln1; S10: x = *a-0.5e0-0.5e0; w = (((((r5*x+r4)*x+r3)*x+r2)*x+r1)*x+r0)/(((((s5*x+s4)*x+s3)*x+s2)*x+s1)*x +1.0e0); gamln1 = x*w; return gamln1; } void grat1(double *a,double *x,double *r,double *p,double *q, double *eps) { static int K2 = 0; static double a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,t,tol,w,z,T1,T3; /* .. .. Executable Statements .. */ /* ----------------------------------------------------------------------- EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS P(A,X) AND Q(A,X) IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED. THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A). ----------------------------------------------------------------------- */ if(*a**x == 0.0e0) goto S120; if(*a == 0.5e0) goto S100; if(*x < 1.1e0) goto S10; goto S60; S10: /* TAYLOR SERIES FOR P(A,X)/X**A */ an = 3.0e0; c = *x; sum = *x/(*a+3.0e0); tol = 0.1e0**eps/(*a+1.0e0); S20: an += 1.0e0; c = -(c*(*x/an)); t = c/(*a+an); sum += t; if(fabs(t) > tol) goto S20; j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0)); z = *a*log(*x); h = gam1(a); g = 1.0e0+h; if(*x < 0.25e0) goto S30; if(*a < *x/2.59e0) goto S50; goto S40; S30: if(z > -.13394e0) goto S50; S40: w = exp(z); *p = w*g*(0.5e0+(0.5e0-j)); *q = 0.5e0+(0.5e0-*p); return; S50: l = rexp(&z); w = 0.5e0+(0.5e0+l); *q = (w*j-l)*g-h; if(*q < 0.0e0) goto S90; *p = 0.5e0+(0.5e0-*q); return; S60: /* CONTINUED FRACTION EXPANSION */ a2nm1 = a2n = 1.0e0; b2nm1 = *x; b2n = *x+(1.0e0-*a); c = 1.0e0; S70: a2nm1 = *x*a2n+c*a2nm1; b2nm1 = *x*b2n+c*b2nm1; am0 = a2nm1/b2nm1; c += 1.0e0; cma = c-*a; a2n = a2nm1+cma*a2n; b2n = b2nm1+cma*b2n; an0 = a2n/b2n; if(fabs(an0-am0) >= *eps*an0) goto S70; *q = *r*an0; *p = 0.5e0+(0.5e0-*q); return; S80: /* SPECIAL CASES */ *p = 0.0e0; *q = 1.0e0; return; S90: *p = 1.0e0; *q = 0.0e0; return; S100: if(*x >= 0.25e0) goto S110; T1 = sqrt(*x); *p = erf1(&T1); *q = 0.5e0+(0.5e0-*p); return; S110: T3 = sqrt(*x); *q = erfc1(&K2,&T3); *p = 0.5e0+(0.5e0-*q); return; S120: if(*x <= *a) goto S80; goto S90; } double gsumln(double *a,double *b) /* ----------------------------------------------------------------------- EVALUATION OF THE FUNCTION LN(GAMMA(A + B)) FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2 ----------------------------------------------------------------------- */ { static double gsumln,x,T1,T2; /* .. .. Executable Statements .. */ x = *a+*b-2.e0; if(x > 0.25e0) goto S10; T1 = 1.0e0+x; gsumln = gamln1(&T1); return gsumln; S10: if(x > 1.25e0) goto S20; gsumln = gamln1(&x)+alnrel(&x); return gsumln; S20: T2 = x-1.0e0; gsumln = gamln1(&T2)+log(x*(1.0e0+x)); return gsumln; } double rexp(double *x) /* ----------------------------------------------------------------------- EVALUATION OF THE FUNCTION EXP(X) - 1 ----------------------------------------------------------------------- */ { static double p1 = .914041914819518e-09; static double p2 = .238082361044469e-01; static double q1 = -.499999999085958e+00; static double q2 = .107141568980644e+00; static double q3 = -.119041179760821e-01; static double q4 = .595130811860248e-03; static double rexp,w; /* .. .. Executable Statements .. */ if(fabs(*x) > 0.15e0) goto S10; rexp = *x*(((p2**x+p1)**x+1.0e0)/((((q4**x+q3)**x+q2)**x+q1)**x+1.0e0)); return rexp; S10: w = exp(*x); if(*x > 0.0e0) goto S20; rexp = w-0.5e0-0.5e0; return rexp; S20: rexp = w*(0.5e0+(0.5e0-1.0e0/w)); return rexp; } double spmpar(int *i) /* ----------------------------------------------------------------------- SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION, SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE. ----------------------------------------------------------------------- WRITTEN BY ALFRED H. MORRIS, JR. NAVAL SURFACE WARFARE CENTER DAHLGREN VIRGINIA ----------------------------------------------------------------------- ----------------------------------------------------------------------- MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE CONSTANTS FOR THE COMPUTER BEING USED. THIS MODIFICATION WAS MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION ----------------------------------------------------------------------- */ { static int K1 = 4; static int K2 = 8; static int K3 = 9; static int K4 = 10; static double spmpar,b,binv,bm1,one,w,z; static int emax,emin,ibeta,m; /* .. .. Executable Statements .. */ if(*i > 1) goto S10; b = ipmpar(&K1); m = ipmpar(&K2); spmpar = pow(b,(double)(1-m)); return spmpar; S10: if(*i > 2) goto S20; b = ipmpar(&K1); emin = ipmpar(&K3); one = 1.0; binv = one/b; w = pow(b,(double)(emin+2)); spmpar = w*binv*binv*binv; return spmpar; S20: ibeta = ipmpar(&K1); m = ipmpar(&K2); emax = ipmpar(&K4); b = ibeta; bm1 = ibeta-1; one = 1.0; z = pow(b,(double)(m-1)); w = ((z-one)*b+bm1)/(b*z); z = pow(b,(double)(emax-2)); spmpar = w*z*b*b; return spmpar; } double psi(double *xx) /* --------------------------------------------------------------------- EVALUATION OF THE DIGAMMA FUNCTION ----------- PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT BE COMPUTED. THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY CODY, STRECOK AND THACHER. --------------------------------------------------------------------- PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY A.H. MORRIS (NSWC). --------------------------------------------------------------------- */ { static double dx0 = 1.461632144968362341262659542325721325e0; static double piov4 = .785398163397448e0; static double p1[7] = { .895385022981970e-02,.477762828042627e+01,.142441585084029e+03, .118645200713425e+04,.363351846806499e+04,.413810161269013e+04, .130560269827897e+04 }; static double p2[4] = { -.212940445131011e+01,-.701677227766759e+01,-.448616543918019e+01, -.648157123766197e+00 }; static double q1[6] = { .448452573429826e+02,.520752771467162e+03,.221000799247830e+04, .364127349079381e+04,.190831076596300e+04,.691091682714533e-05 }; static double q2[4] = { .322703493791143e+02,.892920700481861e+02,.546117738103215e+02, .777788548522962e+01 }; static int K1 = 3; static int K2 = 1; static double psi,aug,den,sgn,upper,w,x,xmax1,xmx0,xsmall,z; static int i,m,n,nq; /* .. .. Executable Statements .. */ /* --------------------------------------------------------------------- MACHINE DEPENDENT CONSTANTS ... XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH PSI MAY BE REPRESENTED AS ALOG(X). XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X) MAY BE REPRESENTED BY 1/X. --------------------------------------------------------------------- */ xmax1 = ipmpar(&K1); xmax1 = fifdmin1(xmax1,1.0e0/spmpar(&K2)); xsmall = 1.e-9; x = *xx; aug = 0.0e0; if(x >= 0.5e0) goto S50; /* --------------------------------------------------------------------- X .LT. 0.5, USE REFLECTION FORMULA PSI(1-X) = PSI(X) + PI * COTAN(PI*X) --------------------------------------------------------------------- */ if(fabs(x) > xsmall) goto S10; if(x == 0.0e0) goto S100; /* --------------------------------------------------------------------- 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE FOR PI*COTAN(PI*X) --------------------------------------------------------------------- */ aug = -(1.0e0/x); goto S40; S10: /* --------------------------------------------------------------------- REDUCTION OF ARGUMENT FOR COTAN --------------------------------------------------------------------- */ w = -x; sgn = piov4; if(w > 0.0e0) goto S20; w = -w; sgn = -sgn; S20: /* --------------------------------------------------------------------- MAKE AN ERROR EXIT IF X .LE. -XMAX1 --------------------------------------------------------------------- */ if(w >= xmax1) goto S100; nq = fifidint(w); w -= (double)nq; nq = fifidint(w*4.0e0); w = 4.0e0*(w-(double)nq*.25e0); /* --------------------------------------------------------------------- W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X. ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST QUADRANT AND DETERMINE SIGN --------------------------------------------------------------------- */ n = nq/2; if(n+n != nq) w = 1.0e0-w; z = piov4*w; m = n/2; if(m+m != n) sgn = -sgn; /* --------------------------------------------------------------------- DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X) --------------------------------------------------------------------- */ n = (nq+1)/2; m = n/2; m += m; if(m != n) goto S30; /* --------------------------------------------------------------------- CHECK FOR SINGULARITY --------------------------------------------------------------------- */ if(z == 0.0e0) goto S100; /* --------------------------------------------------------------------- USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND SIN/COS AS A SUBSTITUTE FOR TAN --------------------------------------------------------------------- */ aug = sgn*(cos(z)/sin(z)*4.0e0); goto S40; S30: aug = sgn*(sin(z)/cos(z)*4.0e0); S40: x = 1.0e0-x; S50: if(x > 3.0e0) goto S70; /* --------------------------------------------------------------------- 0.5 .LE. X .LE. 3.0 --------------------------------------------------------------------- */ den = x; upper = p1[0]*x; for(i=1; i<=5; i++) { den = (den+q1[i-1])*x; upper = (upper+p1[i+1-1])*x; } den = (upper+p1[6])/(den+q1[5]); xmx0 = x-dx0; psi = den*xmx0+aug; return psi; S70: /* --------------------------------------------------------------------- IF X .GE. XMAX1, PSI = LN(X) --------------------------------------------------------------------- */ if(x >= xmax1) goto S90; /* --------------------------------------------------------------------- 3.0 .LT. X .LT. XMAX1 --------------------------------------------------------------------- */ w = 1.0e0/(x*x); den = w; upper = p2[0]*w; for(i=1; i<=3; i++) { den = (den+q2[i-1])*w; upper = (upper+p2[i+1-1])*w; } aug = upper/(den+q2[3])-0.5e0/x+aug; S90: psi = aug+log(x); return psi; S100: /* --------------------------------------------------------------------- ERROR RETURN --------------------------------------------------------------------- */ psi = 0.0e0; return psi; } double rlog1(double *x) /* ----------------------------------------------------------------------- EVALUATION OF THE FUNCTION X - LN(1 + X) ----------------------------------------------------------------------- */ { static double a = .566749439387324e-01; static double b = .456512608815524e-01; static double p0 = .333333333333333e+00; static double p1 = -.224696413112536e+00; static double p2 = .620886815375787e-02; static double q1 = -.127408923933623e+01; static double q2 = .354508718369557e+00; static double rlog1,h,r,t,w,w1; /* .. .. Executable Statements .. */ if(*x < -0.39e0 || *x > 0.57e0) goto S40; if(*x < -0.18e0) goto S10; if(*x > 0.18e0) goto S20; /* ARGUMENT REDUCTION */ h = *x; w1 = 0.0e0; goto S30; S10: h = *x+0.3e0; h /= 0.7e0; w1 = a-h*0.3e0; goto S30; S20: h = 0.75e0**x-0.25e0; w1 = b+h/3.0e0; S30: /* SERIES EXPANSION */ r = h/(h+2.0e0); t = r*r; w = ((p2*t+p1)*t+p0)/((q2*t+q1)*t+1.0e0); rlog1 = 2.0e0*t*(1.0e0/(1.0e0-r)-r*w)+w1; return rlog1; S40: w = *x+0.5e0+0.5e0; rlog1 = *x-log(w); return rlog1; } /************************************************************************ FIFDMAX1: returns the maximum of two numbers a and b ************************************************************************/ double fifdmax1(double a,double b) /* a - first number */ /* b - second number */ { if (a < b) return b; else return a; } /************************************************************************ FIFDMIN1: returns the minimum of two numbers a and b ************************************************************************/ double fifdmin1(double a,double b) /* a - first number */ /* b - second number */ { if (a < b) return a; else return b; } /************************************************************************ FIFIDINT: Truncates a double precision number to a long integer ************************************************************************/ long fifidint(double a) /* a - number to be truncated */ { return (long)(a); } /************************************************************************ FIFDSIGN: transfers the sign of the variable "sign" to the variable "mag" ************************************************************************/ double fifdsign(double mag,double sign) /* mag - magnitude */ /* sign - sign to be transfered */ { if (mag < 0) mag = -mag; if (sign < 0) mag = -mag; return mag; } /************************************************************************ FTNSTOP: Prints msg to standard error and then exits ************************************************************************/ void ftnstop(char* msg) /* msg - error message */ { if (msg != NULL) fprintf(stderr,"%s\n",msg); exit(EXIT_FAILURE); /* EXIT_FAILURE from stdlib.h, or use an int */ }