Changeset 1026Tweet

Timestamp:

Oct 15, 2011, 2:03:35 PM (11 years ago)

Author:

sam

Message:

core: implement the gamma function for reals using Spouge's formula.

Location:

trunk/src

Files:

• real.cpp (modified) (1 diff)
• real.h (modified) (1 diff)

trunk/src/real.cpp

@@ -606 +606 @@
         return real::R_0;
     }
+}
+
+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;
 }
 

trunk/src/real.h

@@ -92 +92 @@
     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. */