1 files changed, 52 insertions, 0 deletions
diff --git a/src/c/elementaryFunctions/special_functions/gamma/dgammas.c b/src/c/elementaryFunctions/special_functions/gamma/dgammas.c
new file mode 100644
index 00000000..c466f434
--- /dev/null
+++ b/src/c/elementaryFunctions/special_functions/gamma/dgammas.c
@@ -0,0 +1,52 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <gsl/gsl_sf_gamma.h>
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+ 
+/* very simple approximation */
+double st_gamma(double x)
+{
+  return sqrt(2.0*M_PI/x)*pow(x/M_E, x);
+}
+ 
+#define A 12
+double sp_gamma(double z)
+{
+  const int a = A;
+  static double c_space[A];
+  static double *c = NULL;
+  int k;
+  double accm;
+ 
+  if ( c == NULL ) {
+    double k1_factrl = 1.0; /* (k - 1)!*(-1)^k with 0!==1*/
+    c = c_space;
+    c[0] = sqrt(2.0*M_PI);
+    for(k=1; k < a; k++) {
+      c[k] = exp(a-k) * pow(a-k, k-0.5) / k1_factrl;
+	  k1_factrl *= -k;
+    }
+  }
+  accm = c[0];
+  for(k=1; k < a; k++) {
+    accm += c[k] / ( z + k );
+  }
+  accm *= exp(-(z+a)) * pow(z+a, z+0.5); /* Gamma(z+1) */
+  return accm/z;
+}
+ 
+int main()
+{
+  double x;
+ 
+ 
+  printf("%15s%15s%15s%15s\n", "Stirling", "Spouge", "GSL", "libm");
+  for(x=1.0; x <= 10.0; x+=1.0) {
+    printf("%15.8lf%15.8lf%15.8lf%15.8lf\n", st_gamma(x/3.0), sp_gamma(x/3.0), 
+	   gsl_sf_gamma(x/3.0), tgamma(x/3.0));
+  }
+  return 0;
+}