1 files changed, 98 insertions, 0 deletions
diff --git a/mdk-stage1/dietlibc/libm/gamma.c b/mdk-stage1/dietlibc/libm/gamma.c
new file mode 100644
index 000000000..d5f3e4275
--- /dev/null
+++ b/mdk-stage1/dietlibc/libm/gamma.c
@@ -0,0 +1,98 @@
+#include "dietlibm.h"
+
+/*--------------------------------------------------------------------------*
+
+Name            gamma, lgamma - gamma function
+
+Usage           double gamma (double x);
+                double lgamma(double x);
+                extern int signgam;
+
+Prototype in    math.h
+
+Description     gamma returns the logarithm of the absolute value of the
+                gamma function. So it is possible ā(x) for very large x.
+                The sign is stored in signgam, a extern variable
+                overwritten during every call to gamma(). lgamma() is
+                a synonym for gamma().
+                You can calculate ā(x) by the following sequence:
+
+                double gammafunction(double x)
+                  { double y=exp(gamma(x));
+
+                    return signgam ? -y : +y;
+                  }
+
+Return value    gamma returns a value in range (-0.1208, +oo). For a input
+                value of zero, it returns +oo and errno is set to:
+
+                        ERANGE  Result out of range
+
+*---------------------------------------------------------------------------*/
+
+#include <stdlib.h>
+#include <math.h>
+
+#define B0      +            1.0l/   6/ 1/ 2
+#define B1      -            1.0l/  30/ 3/ 4
+#define B2      +            1.0l/  42/ 5/ 6
+#define B3      -            1.0l/  30/ 7/ 8
+#define B4      +            5.0l/  66/ 9/10
+#define B5      -          691.0l/2730/11/12
+#define B6      +            7.0l/   6/13/14
+#define B7      -         3617.0l/ 510/15/16
+#define B8      +        43867.0l/ 798/17/18
+#define B9      -       174611.0l/ 330/19/20
+#define B10     +       854513.0l/ 138/21/22
+#define B11     -    236364091.0l/2730/23/24
+#define B12     +      8553103.0l/   6/25/26
+
+static const double  coeff[] = { B0, B1, B2, B3, B4, B5, B6, B7, B8, B9, B10 };
+int                  signgam;
+
+#define EXPL(x) (((short *)&x)[4] & 0x7FFF)
+
+static double  logfact ( long double x )
+{
+    long double   z = 2. * M_PI * x;
+    register int  e = EXPL (x);
+   
+    static unsigned char list [] = { 6, 4, 3, 3, 2, 2 };
+
+    return (log(x) - 1) * x + 0.5*log(z) + __poly (1./(x*x), e<0x4003 ? 10 : (e>0x4008 ? 1 : list [e-0x4003] ), coeff) / x;
+}
+
+
+double  lgamma ( double x )
+{
+    register int  k = floor (x);
+    long double   w;
+    long double   y;
+    long double   z;
+   
+    signgam = 0;
+
+    if ( k >= 7 )
+        return logfact (x-1);
+       
+    if ( k == x )
+        switch (k) {
+        case 1 :
+        case 2 : return 0.000000000000000000000000000l;
+        case 3 : return 0.693147180559945309432805516l;
+        case 4 : return 1.791759469228055000858148560l;
+        case 5 : return 3.178053830347945619723759592l;
+        case 6 : return 4.787491742782045994244981560l;
+        default: return 1./0.; /* ignore the gcc warning, this is intentional */
+        }
+       
+    z = logfact (y = x - k + 7.0 - 1);
+    w = 1;
+    for ( k = 7 - k; k--; )
+        w *= y, y -= 1.;
+       
+    signgam = k >= 0  ?  0  :  k & 1;
+    return z - log (w);
+}
+
+double gamma ( double val )  __attribute__ ((weak,alias("lgamma")));