Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_06.c

00001 #include "cdflib.h"
00002 double betaln(double *a0,double *b0)
00003 /*
00004 -----------------------------------------------------------------------
00005      EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
00006 -----------------------------------------------------------------------
00007      E = 0.5*LN(2*PI)
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      .. 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 */
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

cdf_06.c
