* math.c (math_gamma): new method Math.gamma.

(math_lgamma): new method Math.lgamma.

* include/ruby/missing.h (tgamma): declared unless HAVE_TGAMMA.
  (lgamma_r): declared unless HAVE_LGAMMA_R.

* configure.in (tgamma): check for replacement funtions.
  (lgamma_r): ditto.

* missing/tgamma.c: new file.  based on gamma.c from
  "C-gengo niyoru saishin algorithm jiten" (New Algorithm handbook
  in C language) (Gijyutsu hyouron sha, Tokyo, 1991)
  by Haruhiko Okumura.

* missing/lgamma_r.c: ditto.

* LEGAL (missing/tgamma.c): describe as public domain.
  (missing/lgamma_r.c): ditto.



git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@15388 b2dd03c8-39d4-4d8f-98ff-823fe69b080e

ChangeLog

@@ -1,3 +1,24 @@
+Thu Feb  7 10:39:21 2008  Tanaka Akira  <akr@fsij.org>
+
+	* math.c (math_gamma): new method Math.gamma.
+	  (math_lgamma): new method Math.lgamma.
+
+	* include/ruby/missing.h (tgamma): declared unless HAVE_TGAMMA.
+	  (lgamma_r): declared unless HAVE_LGAMMA_R.
+
+	* configure.in (tgamma): check for replacement funtions.
+	  (lgamma_r): ditto.
+	
+	* missing/tgamma.c: new file.  based on gamma.c from
+	  "C-gengo niyoru saishin algorithm jiten" (New Algorithm handbook
+	  in C language) (Gijyutsu hyouron sha, Tokyo, 1991)
+	  by Haruhiko Okumura.
+
+	* missing/lgamma_r.c: ditto.
+
+	* LEGAL (missing/tgamma.c): describe as public domain.
+	  (missing/lgamma_r.c): ditto.
+
 Thu Feb  7 09:05:57 2008  Yukihiro Matsumoto  <matz@ruby-lang.org>
 
 	* ext/nkf/nkf-utf8/nkf.c (nkf_enc_from_index): BINARY does not

LEGAL

@@ -158,6 +158,8 @@ ext/digest/sha1/sha1.[ch]:
   These files are all under public domain.
 
 missing/erf.c:
+missing/tgamma.c:
+missing/lgamma_r.c:
 missing/crypt.c:
 missing/vsnprintf.c:
 

configure.in

@@ -649,7 +649,8 @@ esac
 AC_FUNC_MEMCMP
 AC_REPLACE_FUNCS(dup2 memmove strerror strftime\
 		 strchr strstr crypt flock vsnprintf\
-		 isnan finite isinf hypot acosh erf strlcpy strlcat)
+		 isnan finite isinf hypot acosh erf tgamma lgamma_r \
+                 strlcpy strlcat)
 AC_CHECK_FUNCS(fmod killpg wait4 waitpid fork spawnv syscall chroot fsync getcwd eaccess\
 	      truncate chsize times utimes utimensat fcntl lockf lstat\
               link symlink readlink\

include/ruby/missing.h

@@ -79,6 +79,14 @@ extern double erf(double);
 extern double erfc(double);
 #endif
 
+#ifndef HAVE_TGAMMA
+extern double tgamma(double);
+#endif
+
+#ifndef HAVE_LGAMMA_R
+extern double lgamma_r(double, int *);
+#endif
+
 #ifndef isinf
 # ifndef HAVE_ISINF
 #  if defined(HAVE_FINITE) && defined(HAVE_ISNAN)

math.c

@@ -486,6 +486,76 @@ math_erfc(VALUE obj, VALUE x)
     return DOUBLE2NUM(erfc(RFLOAT_VALUE(x)));
 }
 
+/*
+ * call-seq:
+ *    Math.gamma(x)  => float
+ *
+ *  Calculates the gamma function of x.
+ *
+ *  Note that gamma(n) is same as fact(n-1) for integer n >= 0.
+ *  However gamma(n) returns float and possibly has error in calculation.
+ *
+ *   def fact(n) (1..n).inject(1) {|r,i| r*i } end
+ *   0.upto(25) {|i| p [i, Math.gamma(i+1), fact(i)] }
+ *   =>
+ *   [0, 1.0, 1]
+ *   [1, 1.0, 1]
+ *   [2, 2.0, 2]
+ *   [3, 6.0, 6]
+ *   [4, 24.0, 24]
+ *   [5, 120.0, 120]
+ *   [6, 720.0, 720]
+ *   [7, 5040.0, 5040]
+ *   [8, 40320.0, 40320]
+ *   [9, 362880.0, 362880]
+ *   [10, 3628800.0, 3628800]
+ *   [11, 39916800.0, 39916800]
+ *   [12, 479001599.999999, 479001600]
+ *   [13, 6227020800.00001, 6227020800]
+ *   [14, 87178291199.9998, 87178291200]
+ *   [15, 1307674368000.0, 1307674368000]
+ *   [16, 20922789888000.0, 20922789888000]
+ *   [17, 3.55687428096001e+14, 355687428096000]
+ *   [18, 6.40237370572799e+15, 6402373705728000]
+ *   [19, 1.21645100408832e+17, 121645100408832000]
+ *   [20, 2.43290200817664e+18, 2432902008176640000]
+ *   [21, 5.10909421717094e+19, 51090942171709440000]
+ *   [22, 1.12400072777761e+21, 1124000727777607680000]
+ *   [23, 2.58520167388851e+22, 25852016738884976640000]
+ *   [24, 6.20448401733239e+23, 620448401733239439360000]
+ *   [25, 1.5511210043331e+25, 15511210043330985984000000]
+ *
+ */
+
+static VALUE
+math_gamma(VALUE obj, VALUE x)
+{
+    Need_Float(x);
+    return DOUBLE2NUM(tgamma(RFLOAT_VALUE(x)));
+}
+
+/*
+ * call-seq:
+ *    Math.lgamma(x)  => [float, -1 or 1]
+ *
+ *  Calculates the logarithmic gamma of x and
+ *  the sign of gamma of x.
+ *
+ *  Math.lgamma(x) is same as
+ *   [Math.log(Math.gamma(x)), Math.gamma(x) < 0 ? -1 : 1]
+ *  but avoid overflow by Math.gamma(x) for large x.
+ */
+
+static VALUE
+math_lgamma(VALUE obj, VALUE x)
+{
+    int sign;
+    VALUE v;
+    Need_Float(x);
+    v = DOUBLE2NUM(lgamma_r(RFLOAT_VALUE(x), &sign));
+    return rb_assoc_new(v, INT2FIX(sign));
+}
+
 /*
  *  The <code>Math</code> module contains module functions for basic
  *  trigonometric and transcendental functions. See class
@@ -541,4 +611,7 @@ Init_Math(void)
 
     rb_define_module_function(rb_mMath, "erf",  math_erf,  1);
     rb_define_module_function(rb_mMath, "erfc", math_erfc, 1);
+
+    rb_define_module_function(rb_mMath, "gamma", math_gamma, 1);
+    rb_define_module_function(rb_mMath, "lgamma", math_lgamma, 1);
 }

missing/lgamma_r.c

@@ -0,0 +1,64 @@
+/* lgamma_r.c  - public domain implementation of error function lgamma_r(3m)
+
+lgamma_r() is based on gamma().  modified by Tanaka Akira.
+
+reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
+            (New Algorithm handbook in C language) (Gijyutsu hyouron
+            sha, Tokyo, 1991) [in Japanese]
+            http://oku.edu.mie-u.ac.jp/~okumura/algo/
+*/
+
+/***********************************************************
+    gamma.c -- Gamma function
+***********************************************************/
+#include <math.h>
+#define PI      3.14159265358979324  /* $\pi$ */
+#define LOG_2PI 1.83787706640934548  /* $\log 2\pi$ */
+#define LOG_PI  1.14472988584940017  /* $\log_e \pi$ */
+#define N       8
+
+#define B0  1                 /* Bernoulli numbers */
+#define B1  (-1.0 / 2.0)
+#define B2  ( 1.0 / 6.0)
+#define B4  (-1.0 / 30.0)
+#define B6  ( 1.0 / 42.0)
+#define B8  (-1.0 / 30.0)
+#define B10 ( 5.0 / 66.0)
+#define B12 (-691.0 / 2730.0)
+#define B14 ( 7.0 / 6.0)
+#define B16 (-3617.0 / 510.0)
+
+static double
+loggamma(double x)  /* the natural logarithm of the Gamma function. */
+{
+    double v, w;
+
+    v = 1;
+    while (x < N) {  v *= x;  x++;  }
+    w = 1 / (x * x);
+    return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w
+                + (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w
+                + (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w
+                + (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x
+                + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x);
+}
+
+/* the natural logarithm of the absolute value of the Gamma function */
+double
+lgamma_r(double x, int *signp)
+{
+    if (x < 0) {
+        double i, f, s;
+        f = modf(-x, &i);
+        if (f == 0.0) {
+            *signp = 1;
+            return 1.0/0.0;
+        }
+        *signp = (fmod(i, 2.0) != 0.0) ? 1 : -1;
+        s = sin(PI * x);
+        if (s < 0) s = -s;
+        return LOG_PI - log(s) - loggamma(1 - x);
+    }
+    *signp = 1;
+    return loggamma(x);
+}

missing/tgamma.c

@@ -0,0 +1,49 @@
+/* tgamma.c  - public domain implementation of error function tgamma(3m)
+
+reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
+            (New Algorithm handbook in C language) (Gijyutsu hyouron
+            sha, Tokyo, 1991) [in Japanese]
+            http://oku.edu.mie-u.ac.jp/~okumura/algo/
+*/
+
+/***********************************************************
+    gamma.c -- Gamma function
+***********************************************************/
+#include <math.h>
+#define PI      3.14159265358979324  /* $\pi$ */
+#define LOG_2PI 1.83787706640934548  /* $\log 2\pi$ */
+#define N       8
+
+#define B0  1                 /* Bernoulli numbers */
+#define B1  (-1.0 / 2.0)
+#define B2  ( 1.0 / 6.0)
+#define B4  (-1.0 / 30.0)
+#define B6  ( 1.0 / 42.0)
+#define B8  (-1.0 / 30.0)
+#define B10 ( 5.0 / 66.0)
+#define B12 (-691.0 / 2730.0)
+#define B14 ( 7.0 / 6.0)
+#define B16 (-3617.0 / 510.0)
+
+static double
+loggamma(double x)  /* the natural logarithm of the Gamma function. */
+{
+    double v, w;
+
+    v = 1;
+    while (x < N) {  v *= x;  x++;  }
+    w = 1 / (x * x);
+    return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w
+                + (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w
+                + (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w
+                + (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x
+                + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x);
+}
+
+double tgamma(double x)  /* Gamma function */
+{
+    if (x < 0)
+        return PI / (sin(PI * x) * exp(loggamma(1 - x)));
+    return exp(loggamma(x));
+}
+