0001-Add-gamma-and-lgamma.patch

From 4301f2adf0bb762f70995dacc2c57ae478a602b4 Mon Sep 17 00:00:00 2001
From: Alan Carvalho de Assis <acassis@gmail.com>
Date: Wed, 8 Mar 2017 16:55:49 -0300
Subject: [PATCH] Add gamma() and lgamma()

---
 include/cxx/cmath        |   2 +
 include/nuttx/lib/math.h |   5 ++
 libc/math/Make.defs      |   2 +-
 libc/math/lib_gamma.c    | 177 +++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 185 insertions(+), 1 deletion(-)
 create mode 100644 libc/math/lib_gamma.c

diff --git a/include/cxx/cmath b/include/cxx/cmath
index 998406c..b635505 100644
--- a/include/cxx/cmath
+++ b/include/cxx/cmath
@@ -103,6 +103,8 @@ namespace std
   using ::sqrt;
   using ::tan;
   using ::tanh;
+  using ::gamma;
+  using ::lgamma;
 #endif
 
 #ifdef CONFIG_HAVE_LONG_DOUBLE
diff --git a/include/nuttx/lib/math.h b/include/nuttx/lib/math.h
index 64db8d7..a394d54 100644
--- a/include/nuttx/lib/math.h
+++ b/include/nuttx/lib/math.h
@@ -212,6 +212,11 @@ long double expl  (long double x);
 #define expm1l(x) (expl(x) - 1.0)
 #endif
 
+#ifdef CONFIG_HAVE_DOUBLE
+double      gamma(double x);
+double      lgamma(double x);
+#endif
+
 float       logf  (float x);
 #ifdef CONFIG_HAVE_DOUBLE
 double      log   (double x);
diff --git a/libc/math/Make.defs b/libc/math/Make.defs
index 7ead148..bb97849 100644
--- a/libc/math/Make.defs
+++ b/libc/math/Make.defs
@@ -57,7 +57,7 @@ CSRCS += lib_fmodl.c lib_frexpl.c lib_ldexpl.c lib_logl.c lib_log10l.c
 CSRCS += lib_log2l.c lib_modfl.c lib_powl.c lib_rintl.c lib_roundl.c
 CSRCS += lib_sinl.c lib_sinhl.c lib_sqrtl.c lib_tanl.c lib_tanhl.c
 CSRCS += lib_asinhl.c lib_acoshl.c lib_atanhl.c lib_erfl.c lib_copysignl.c
-CSRCS += lib_truncl.c
+CSRCS += lib_truncl.c lib_gamma.c
 
 CSRCS += lib_libexpi.c lib_libsqrtapprox.c
 CSRCS += lib_libexpif.c
diff --git a/libc/math/lib_gamma.c b/libc/math/lib_gamma.c
new file mode 100644
index 0000000..d6586c2
--- /dev/null
+++ b/libc/math/lib_gamma.c
@@ -0,0 +1,177 @@
+// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.
+
+#include <math.h>
+#include <debug.h>
+#include <nuttx/lib/float.h>
+
+// Note that the functions Gamma and LogGamma are mutually dependent.
+
+double gamma(double x)
+{
+  if (x <= 0.0)
+    {
+      _err("x needs to be >= 0\n");
+      return 0;
+    }
+
+  /* Split the function domain into three intervals:
+   * (0, 0.001), [0.001, 12), and (12, infinity)
+   */
+
+  /* First interval: (0, 0.001)
+   *
+   * For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
+   * So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
+   * The relative error over this interval is less than 6e-7.
+   */
+
+  /* Euler's gamma constant */
+  const double gamma = 0.577215664901532860606512090;
+
+  if (x < 0.001)
+    return 1.0/(x*(1.0 + gamma*x));
+
+  /* Second interval: [0.001, 12) */
+    
+  if (x < 12.0)
+    {
+      /* The algorithm directly approximates gamma over (1,2) and uses
+       * reduction identities to reduce other arguments to this interval.
+       */
+		
+      double y = x;
+      int n = 0;
+      bool arg_was_less_than_one = (y < 1.0);
+
+      /* Add or subtract integers as necessary to bring y into (1,2)
+       * Will correct for this below
+       */
+
+      if (arg_was_less_than_one)
+        {
+          y += 1.0;
+        }
+      else
+        {
+          /* will use n later */
+          n = (int) (floor(y)) - 1;
+          y -= n;
+        }
+
+      /* numerator coefficients for approximation over the interval (1,2) */
+      static const double p[] =
+        {
+            -1.71618513886549492533811E+0,
+             2.47656508055759199108314E+1,
+            -3.79804256470945635097577E+2,
+             6.29331155312818442661052E+2,
+             8.66966202790413211295064E+2,
+            -3.14512729688483675254357E+4,
+            -3.61444134186911729807069E+4,
+             6.64561438202405440627855E+4
+        };
+
+      /* denominator coefficients for approximation over the interval (1,2) */
+      static const double q[] =
+        {
+            -3.08402300119738975254353E+1,
+             3.15350626979604161529144E+2,
+            -1.01515636749021914166146E+3,
+            -3.10777167157231109440444E+3,
+             2.25381184209801510330112E+4,
+             4.75584627752788110767815E+3,
+            -1.34659959864969306392456E+5,
+            -1.15132259675553483497211E+5
+        };
+
+      double num = 0.0;
+      double den = 1.0;
+      int i;
+
+      double z = y - 1;
+      for (i = 0; i < 8; i++)
+        {
+            num = (num + p[i])*z;
+            den = den*z + q[i];
+        }
+      double result = num/den + 1.0;
+
+      /* Apply correction if argument was not initially in (1,2) */
+      if (arg_was_less_than_one)
+        {
+          /* Use identity gamma(z) = gamma(z+1)/z
+           * The variable "result" now holds gamma of the original y + 1
+           * Thus we use y-1 to get back the orginal y.
+           */
+          result /= (y-1.0);
+        }
+      else
+        {
+          /* Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z) */
+          for (i = 0; i < n; i++)
+            {
+              result *= y++;
+            }
+        }
+
+      return result;
+    }
+
+    /* Third interval: [12, infinity) */
+
+    if (x > 171.624)
+      {
+        /* Correct answer too large to display. Force +infinity. */
+        double temp = DBL_MAX;
+        return temp*2.0;
+    }
+
+    return exp(lgamma(x));
+}
+
+double lgamma(double x)
+{
+  if (x <= 0.0)
+    {
+      _err("x needs to be >= 0\n");
+      return 0;
+    }
+
+  if (x < 12.0)
+    {
+      return log(fabs(gamma(x)));
+    }
+
+  /* Abramowitz and Stegun 6.1.41
+   * Asymptotic series should be good to at least 11 or 12 figures
+   * For error analysis, see Whittiker and Watson
+   * A Course in Modern Analysis (1927), page 252
+   */
+
+  static const double c[8] =
+    {
+		 1.0/12.0,
+		-1.0/360.0,
+		1.0/1260.0,
+		-1.0/1680.0,
+		1.0/1188.0,
+		-691.0/360360.0,
+		1.0/156.0,
+		-3617.0/122400.0
+    };
+
+  int i;
+  double z = 1.0/(x*x);
+  double sum = c[7];
+  for (i=6; i >= 0; i--)
+    {
+      sum *= z;
+      sum += c[i];
+    }
+  double series = sum/x;
+
+  static const double halflogtwopi = 0.91893853320467274178032973640562;
+  double loggamma = (x - 0.5)*log(x) - x + halflogtwopi + series;
+  return loggamma;
+}
+
-- 
2.7.4