Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_62.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 double Xgamm(double *a)
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 {
00020 static double d = .41893853320467274178e0;
00021 static double pi = 3.1415926535898e0;
00022 static double r1 = .820756370353826e-03;
00023 static double r2 = -.595156336428591e-03;
00024 static double r3 = .793650663183693e-03;
00025 static double r4 = -.277777777770481e-02;
00026 static double r5 = .833333333333333e-01;
00027 static double p[7] = {
00028 .539637273585445e-03,.261939260042690e-02,.204493667594920e-01,
00029 .730981088720487e-01,.279648642639792e+00,.553413866010467e+00,1.0e0
00030 };
00031 static double q[7] = {
00032 -.832979206704073e-03,.470059485860584e-02,.225211131035340e-01,
00033 -.170458969313360e+00,-.567902761974940e-01,.113062953091122e+01,1.0e0
00034 };
00035 static int K2 = 3;
00036 static int K3 = 0;
00037 static double Xgamm,bot,g,lnx,s,t,top,w,x,z;
00038 static int i,j,m,n,T1;
00039
00040
00041
00042
00043 Xgamm = 0.0e0;
00044 x = *a;
00045 if(fabs(*a) >= 15.0e0) goto S110;
00046
00047
00048
00049
00050
00051 t = 1.0e0;
00052 m = fifidint(*a)-1;
00053
00054
00055
00056 T1 = m;
00057 if(T1 < 0) goto S40;
00058 else if(T1 == 0) goto S30;
00059 else goto S10;
00060 S10:
00061 for(j=1; j<=m; j++) {
00062 x -= 1.0e0;
00063 t = x*t;
00064 }
00065 S30:
00066 x -= 1.0e0;
00067 goto S80;
00068 S40:
00069
00070
00071
00072 t = *a;
00073 if(*a > 0.0e0) goto S70;
00074 m = -m-1;
00075 if(m == 0) goto S60;
00076 for(j=1; j<=m; j++) {
00077 x += 1.0e0;
00078 t = x*t;
00079 }
00080 S60:
00081 x += (0.5e0+0.5e0);
00082 t = x*t;
00083 if(t == 0.0e0) return Xgamm;
00084 S70:
00085
00086
00087
00088
00089 if(fabs(t) >= 1.e-30) goto S80;
00090 if(fabs(t)*spmpar(&K2) <= 1.0001e0) return Xgamm;
00091 Xgamm = 1.0e0/t;
00092 return Xgamm;
00093 S80:
00094
00095
00096
00097 top = p[0];
00098 bot = q[0];
00099 for(i=1; i<7; i++) {
00100 top = p[i]+x*top;
00101 bot = q[i]+x*bot;
00102 }
00103 Xgamm = top/bot;
00104
00105
00106
00107 if(*a < 1.0e0) goto S100;
00108 Xgamm *= t;
00109 return Xgamm;
00110 S100:
00111 Xgamm /= t;
00112 return Xgamm;
00113 S110:
00114
00115
00116
00117
00118
00119 if(fabs(*a) >= 1.e3) return Xgamm;
00120 if(*a > 0.0e0) goto S120;
00121 x = -*a;
00122 n = x;
00123 t = x-(double)n;
00124 if(t > 0.9e0) t = 1.0e0-t;
00125 s = sin(pi*t)/pi;
00126 if(fifmod(n,2) == 0) s = -s;
00127 if(s == 0.0e0) return Xgamm;
00128 S120:
00129
00130
00131
00132 t = 1.0e0/(x*x);
00133 g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x;
00134
00135
00136
00137
00138 lnx = log(x);
00139
00140
00141
00142 z = x;
00143 g = d+g+(z-0.5e0)*(lnx-1.e0);
00144 w = g;
00145 t = g-w;
00146 if(w > 0.99999e0*exparg(&K3)) return Xgamm;
00147 Xgamm = exp(w)*(1.0e0+t);
00148 if(*a < 0.0e0) Xgamm = 1.0e0/(Xgamm*s)/x;
00149 return Xgamm;
00150 }