Doxygen Source Code Documentation
cdf_46.c File Reference
#include "cdflib.h"Go to the source code of this file.
Functions | |
| double | dlnbet (double *a0, double *b0) |
Function Documentation
|
||||||||||||
|
Definition at line 2 of file cdf_46.c. References a, algdiv(), alnrel(), bcorr(), c, dlnbet(), fifdmax1(), fifdmin1(), gamln(), gsumln(), i, and v. Referenced by dlnbet().
00029 : 00030 DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant 00031 Digit Computation of the Incomplete Beta Function Ratios. ACM 00032 Trans. Math. Softw. 18 (1993), 360-373. 00033 00034 ********************************************************************** 00035 ----------------------------------------------------------------------- 00036 EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION 00037 ----------------------------------------------------------------------- 00038 E = 0.5*LN(2*PI) 00039 -------------------------- 00040 */ 00041 { 00042 static double e = .918938533204673e0; 00043 static double dlnbet,a,b,c,h,u,v,w,z; 00044 static int i,n; 00045 static double T1; 00046 /* 00047 .. 00048 .. Executable Statements .. 00049 */ 00050 a = fifdmin1(*a0,*b0); 00051 b = fifdmax1(*a0,*b0); 00052 if(a >= 8.0e0) goto S100; 00053 if(a >= 1.0e0) goto S20; 00054 /* 00055 ----------------------------------------------------------------------- 00056 PROCEDURE WHEN A .LT. 1 00057 ----------------------------------------------------------------------- 00058 */ 00059 if(b >= 8.0e0) goto S10; 00060 T1 = a+b; 00061 dlnbet = gamln(&a)+(gamln(&b)-gamln(&T1)); 00062 return dlnbet; 00063 S10: 00064 dlnbet = gamln(&a)+algdiv(&a,&b); 00065 return dlnbet; 00066 S20: 00067 /* 00068 ----------------------------------------------------------------------- 00069 PROCEDURE WHEN 1 .LE. A .LT. 8 00070 ----------------------------------------------------------------------- 00071 */ 00072 if(a > 2.0e0) goto S40; 00073 if(b > 2.0e0) goto S30; 00074 dlnbet = gamln(&a)+gamln(&b)-gsumln(&a,&b); 00075 return dlnbet; 00076 S30: 00077 w = 0.0e0; 00078 if(b < 8.0e0) goto S60; 00079 dlnbet = gamln(&a)+algdiv(&a,&b); 00080 return dlnbet; 00081 S40: 00082 /* 00083 REDUCTION OF A WHEN B .LE. 1000 00084 */ 00085 if(b > 1000.0e0) goto S80; 00086 n = a-1.0e0; 00087 w = 1.0e0; 00088 for(i=1; i<=n; i++) { 00089 a -= 1.0e0; 00090 h = a/b; 00091 w *= (h/(1.0e0+h)); 00092 } 00093 w = log(w); 00094 if(b < 8.0e0) goto S60; 00095 dlnbet = w+gamln(&a)+algdiv(&a,&b); 00096 return dlnbet; 00097 S60: 00098 /* 00099 REDUCTION OF B WHEN B .LT. 8 00100 */ 00101 n = b-1.0e0; 00102 z = 1.0e0; 00103 for(i=1; i<=n; i++) { 00104 b -= 1.0e0; 00105 z *= (b/(a+b)); 00106 } 00107 dlnbet = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b))); 00108 return dlnbet; 00109 S80: 00110 /* 00111 REDUCTION OF A WHEN B .GT. 1000 00112 */ 00113 n = a-1.0e0; 00114 w = 1.0e0; 00115 for(i=1; i<=n; i++) { 00116 a -= 1.0e0; 00117 w *= (a/(1.0e0+a/b)); 00118 } 00119 dlnbet = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b)); 00120 return dlnbet; 00121 S100: 00122 /* 00123 ----------------------------------------------------------------------- 00124 PROCEDURE WHEN A .GE. 8 00125 ----------------------------------------------------------------------- 00126 */ 00127 w = bcorr(&a,&b); 00128 h = a/b; 00129 c = h/(1.0e0+h); 00130 u = -((a-0.5e0)*log(c)); 00131 v = b*alnrel(&h); 00132 if(u <= v) goto S110; 00133 dlnbet = -(0.5e0*log(b))+e+w-v-u; 00134 return dlnbet; 00135 S110: 00136 dlnbet = -(0.5e0*log(b))+e+w-u-v; 00137 return dlnbet; } /* END */ |
