00001 #include "cdflib.h"
00002 double Xgamm(double *a)
00003 /*
00004 -----------------------------------------------------------------------
00005  
00006          EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS
00007  
00008                            -----------
00009  
00010      GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
00011      BE COMPUTED.
00012  
00013 -----------------------------------------------------------------------
00014      WRITTEN BY ALFRED H. MORRIS, JR.
00015           NAVAL SURFACE WEAPONS CENTER
00016           DAHLGREN, VIRGINIA
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      .. Executable Statements ..
00042 */
00043     Xgamm = 0.0e0;
00044     x = *a;
00045     if(fabs(*a) >= 15.0e0) goto S110;
00046 /*
00047 -----------------------------------------------------------------------
00048             EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15
00049 -----------------------------------------------------------------------
00050 */
00051     t = 1.0e0;
00052     m = fifidint(*a)-1;
00053 /*
00054      LET T BE THE PRODUCT OF A-J WHEN A .GE. 2
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      LET T BE THE PRODUCT OF A+J WHEN A .LT. 1
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      THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
00087      CODE MAY BE OMITTED IF DESIRED.
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      COMPUTE GAMMA(1 + X) FOR  0 .LE. X .LT. 1
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      TERMINATION
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             EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15
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      COMPUTE THE MODIFIED ASYMPTOTIC SUM
00131 */
00132     t = 1.0e0/(x*x);
00133     g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x;
00134 /*
00135      ONE MAY REPLACE THE NEXT STATEMENT WITH  LNX = ALOG(X)
00136      BUT LESS ACCURACY WILL NORMALLY BE OBTAINED.
00137 */
00138     lnx = log(x);
00139 /*
00140      FINAL ASSEMBLY
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 } /* END */