[ruby-dev:33620] Math.lgamma and Math.gamma

ふと、log2(nCm) のグラフを書いてみようと思ったのですが、nCm
の計算に階乗が必要で、これをふつうに計算するのはループが必要
でかなり遅く、Γ関数が (正確にはオーバーフローを避ける必要上、
Γ関数の対数が) あれば高速に計算できるのに、と残念な思いをし
ました。

というわけで、
  Math.gamma(x)   => float
  Math.lgamma(x)  => [float, sign_of_gamma_x]
というのを実装してみたんですがどうでしょうか。

欲しかったのは lgamma なんですが、lgamma だけで gamma がない
のもなんなので gamma も実装してあります。

使用する関数は tgamma と lgamma_r です。

tgamma は C99 で定義されてます。

lgamma_r は規格にありません。規格にあるのは lgamma で、これ
はスレッドセーフじゃないので避けて、lgamma_r を使っています。
lgamma_r は GNU/Linux, NetBSD, Solaris などにあります。

missing/tgamma.c は「C言語による最新アルゴリズム事典」の
サポートページの http://oku.edu.mie-u.ac.jp/~okumura/algo/
から取得した algo.tar.gz 内の gamma.c を元にして関数名やコメ
ントなどを調整したものです。

missing/lgamma_r.c は tgamma.c をさらにいじって作りました。

algo.tar.gz 中の README に
| ご利用についての制限はありません。ただしバグによる損害の
| 賠償などには応じかねますのでご了承ください。
とあったのと、やはり「C言語による最新アルゴリズム事典」から
とられている erf.c に public domain と書いてあったので、
missing/tgamma.c, missing/lgamma_r.c も public domain として
あります。

Index: math.c
===================================================================
--- math.c	(revision 15384)
+++ math.c	(working copy)
@@ -487,6 +487,42 @@ math_erfc(VALUE obj, VALUE x)
 }
 
 /*
+ * call-seq:
+ *    Math.gamma(x)  => float
+ *
+ *  Calculates the gamma function of x.
+ */
+
+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
  *  <code>Float</code> for a list of constants that
@@ -541,4 +577,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);
 }
Index: include/ruby/missing.h
===================================================================
--- include/ruby/missing.h	(revision 15384)
+++ include/ruby/missing.h	(working copy)
@@ -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)
Index: configure.in
===================================================================
--- configure.in	(revision 15384)
+++ configure.in	(working copy)
@@ -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\
Index: missing/tgamma.c
===================================================================
--- missing/tgamma.c	(revision 0)
+++ missing/tgamma.c	(revision 0)
@@ -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));
+}
+
Index: missing/lgamma_r.c
===================================================================
--- missing/lgamma_r.c	(revision 0)
+++ missing/lgamma_r.c	(revision 0)
@@ -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);
+}
Index: LEGAL
===================================================================
--- LEGAL	(revision 15384)
+++ LEGAL	(working copy)
@@ -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:
 
-- 
[田中 哲][たなか あきら][Tanaka Akira]