Doxygen Source Code Documentation

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

00001 #include "cdflib.h"
00002 double alngam(double *x)
00003 /*
00004 **********************************************************************
00005  
00006      double alngam(double *x)
00007                  double precision LN of the GAMma function
00008  
00009  
00010                               Function
00011  
00012  
00013      Returns the natural logarithm of GAMMA(X).
00014  
00015  
00016                               Arguments
00017  
00018  
00019      X --> value at which scaled log gamma is to be returned
00020                     X is DOUBLE PRECISION
00021  
00022  
00023                               Method
00024  
00025  
00026      If X .le. 6.0, then use recursion to get X below 3
00027      then apply rational approximation number 5236 of
00028      Hart et al, Computer Approximations, John Wiley and
00029      Sons, NY, 1968.
00030  
00031      If X .gt. 6.0, then use recursion to get X to at least 12 and
00032      then use formula 5423 of the same source.
00033  
00034 **********************************************************************
00035 */
00036 {
00037 #define hln2pi 0.91893853320467274178e0
00038 static double coef[5] = {
00039     0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
00040     -0.594997310889e-3,0.8065880899e-3
00041 };
00042 static double scoefd[4] = {
00043     0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
00044     0.1000000000000000000e1
00045 };
00046 static double scoefn[9] = {
00047     0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
00048     0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
00049     0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
00050 };
00051 static int K1 = 9;
00052 static int K3 = 4;
00053 static int K5 = 5;
00054 static double alngam,offset,prod,xx;
00055 static int i,n;
00056 static double T2,T4,T6;
00057 /*
00058      ..
00059      .. Executable Statements ..
00060 */
00061     if(!(*x <= 6.0e0)) goto S70;
00062     prod = 1.0e0;
00063     xx = *x;
00064     if(!(*x > 3.0e0)) goto S30;
00065 S10:
00066     if(!(xx > 3.0e0)) goto S20;
00067     xx -= 1.0e0;
00068     prod *= xx;
00069     goto S10;
00070 S30:
00071 S20:
00072     if(!(*x < 2.0e0)) goto S60;
00073 S40:
00074     if(!(xx < 2.0e0)) goto S50;
00075     prod /= xx;
00076     xx += 1.0e0;
00077     goto S40;
00078 S60:
00079 S50:
00080     T2 = xx-2.0e0;
00081     T4 = xx-2.0e0;
00082     alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
00083 /*
00084      COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
00085 */
00086     alngam *= prod;
00087     alngam = log(alngam);
00088     goto S110;
00089 S70:
00090     offset = hln2pi;
00091 /*
00092      IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
00093 */
00094     n = fifidint(12.0e0-*x);
00095     if(!(n > 0)) goto S90;
00096     prod = 1.0e0;
00097     for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));
00098     offset -= log(prod);
00099     xx = *x+(double)n;
00100     goto S100;
00101 S90:
00102     xx = *x;
00103 S100:
00104 /*
00105      COMPUTE POWER SERIES
00106 */
00107     T6 = 1.0e0/pow(xx,2.0);
00108     alngam = devlpl(coef,&K5,&T6)/xx;
00109     alngam += (offset+(xx-0.5e0)*log(xx)-xx);
00110 S110:
00111     return alngam;
00112 #undef hln2pi
00113 } /* END */
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

cdf_01.c
