Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_01.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 double alngam(double *x)
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
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
00085
00086 alngam *= prod;
00087 alngam = log(alngam);
00088 goto S110;
00089 S70:
00090 offset = hln2pi;
00091
00092
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
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 }