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 */
|
