Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_06.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 double betaln(double *a0,double *b0)
00003
00004
00005
00006
00007
00008
00009
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
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
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
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
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
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
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
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 }