Doxygen Source Code Documentation
cdf_06.c File Reference
#include "cdflib.h"
Go to the source code of this file.
Functions | |
double | betaln (double *a0, double *b0) |
Function Documentation
|
Definition at line 2 of file cdf_06.c. References a, algdiv(), alnrel(), bcorr(), betaln(), c, fifdmax1(), fifdmin1(), gamln(), gsumln(), i, and v. Referenced by betaln(), bpser(), brcmp1(), and brcomp().
00010 { 00011 static double e = .918938533204673e0; 00012 static double betaln,a,b,c,h,u,v,w,z; 00013 static int i,n; 00014 static double T1; 00015 /* 00016 .. 00017 .. Executable Statements .. 00018 */ 00019 a = fifdmin1(*a0,*b0); 00020 b = fifdmax1(*a0,*b0); 00021 if(a >= 8.0e0) goto S100; 00022 if(a >= 1.0e0) goto S20; 00023 /* 00024 ----------------------------------------------------------------------- 00025 PROCEDURE WHEN A .LT. 1 00026 ----------------------------------------------------------------------- 00027 */ 00028 if(b >= 8.0e0) goto S10; 00029 T1 = a+b; 00030 betaln = gamln(&a)+(gamln(&b)-gamln(&T1)); 00031 return betaln; 00032 S10: 00033 betaln = gamln(&a)+algdiv(&a,&b); 00034 return betaln; 00035 S20: 00036 /* 00037 ----------------------------------------------------------------------- 00038 PROCEDURE WHEN 1 .LE. A .LT. 8 00039 ----------------------------------------------------------------------- 00040 */ 00041 if(a > 2.0e0) goto S40; 00042 if(b > 2.0e0) goto S30; 00043 betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b); 00044 return betaln; 00045 S30: 00046 w = 0.0e0; 00047 if(b < 8.0e0) goto S60; 00048 betaln = gamln(&a)+algdiv(&a,&b); 00049 return betaln; 00050 S40: 00051 /* 00052 REDUCTION OF A WHEN B .LE. 1000 00053 */ 00054 if(b > 1000.0e0) goto S80; 00055 n = a-1.0e0; 00056 w = 1.0e0; 00057 for(i=1; i<=n; i++) { 00058 a -= 1.0e0; 00059 h = a/b; 00060 w *= (h/(1.0e0+h)); 00061 } 00062 w = log(w); 00063 if(b < 8.0e0) goto S60; 00064 betaln = w+gamln(&a)+algdiv(&a,&b); 00065 return betaln; 00066 S60: 00067 /* 00068 REDUCTION OF B WHEN B .LT. 8 00069 */ 00070 n = b-1.0e0; 00071 z = 1.0e0; 00072 for(i=1; i<=n; i++) { 00073 b -= 1.0e0; 00074 z *= (b/(a+b)); 00075 } 00076 betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b))); 00077 return betaln; 00078 S80: 00079 /* 00080 REDUCTION OF A WHEN B .GT. 1000 00081 */ 00082 n = a-1.0e0; 00083 w = 1.0e0; 00084 for(i=1; i<=n; i++) { 00085 a -= 1.0e0; 00086 w *= (a/(1.0e0+a/b)); 00087 } 00088 betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b)); 00089 return betaln; 00090 S100: 00091 /* 00092 ----------------------------------------------------------------------- 00093 PROCEDURE WHEN A .GE. 8 00094 ----------------------------------------------------------------------- 00095 */ 00096 w = bcorr(&a,&b); 00097 h = a/b; 00098 c = h/(1.0e0+h); 00099 u = -((a-0.5e0)*log(c)); 00100 v = b*alnrel(&h); 00101 if(u <= v) goto S110; 00102 betaln = -(0.5e0*log(b))+e+w-v-u; 00103 return betaln; 00104 S110: 00105 betaln = -(0.5e0*log(b))+e+w-u-v; 00106 return betaln; 00107 } /* END */ |