From 4e4e800572be60f6c64390eeb13faa41cd7ba1ea Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Sat, 15 Oct 2011 12:03:35 +0000
Subject: [PATCH] core: implement the gamma function for reals using Spouge's
 formula.

---
 src/real.cpp | 27 +++++++++++++++++++++++++++
 src/real.h   |  1 +
 2 files changed, 28 insertions(+)

diff --git a/src/real.cpp b/src/real.cpp
index 1888b37c..4aef4154 100644
--- a/src/real.cpp
+++ b/src/real.cpp
@@ -607,6 +607,33 @@ real pow(real const &x, real const &y)
     }
 }
 
+real gamma(real const &x)
+{
+    /* We use Spouge's formula. FIXME: precision is far from acceptable,
+     * especially with large values. We need to compute this with higher
+     * precision values in order to attain the desired accuracy. It might
+     * also be useful to sort the ck values by decreasing absolute value
+     * and do the addition in this order. */
+    int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS);
+
+    real ret = sqrt(real::R_PI << 1);
+    real fact_k_1 = real::R_1;
+
+    for (int k = 1; k < a; k++)
+    {
+        real a_k = (real)(a - k);
+        real ck = pow(a_k, (real)((float)k - 0.5)) * exp(a_k)
+                / (fact_k_1 * (x + (real)(k - 1)));
+        ret += ck;
+        fact_k_1 *= (real)-k;
+    }
+
+    ret *= pow(x + (real)(a - 1), x - (real::R_1 >> 1));
+    ret *= exp(-x - (real)(a - 1));
+
+    return ret;
+}
+
 real fabs(real const &x)
 {
     real ret = x;
diff --git a/src/real.h b/src/real.h
index 9340660d..a6de7e93 100644
--- a/src/real.h
+++ b/src/real.h
@@ -91,6 +91,7 @@ public:
     friend real sqrt(real const &x);
     friend real cbrt(real const &x);
     friend real pow(real const &x, real const &y);
+    friend real gamma(real const &x);
 
     /* Rounding, absolute value, remainder etc. */
     friend real ceil(real const &x);