e$B$U$H!“e(Blog2(nCm) e$B$N%0%i%U$r=q$$$F$_$h$&$H;W$C$?$N$G$9$,!“e(BnCm
e$B$N7W;;$K3,>h$,I,MW$G!”$3$l$r$U$D$&$K7W;;$9$k$N$O%k!<%W$,I,MWe(B
e$B$G$+$J$jCY$/!”X?t$,e(B
(e$B@53N$K$O%*!<%P!<%U%m!<$rHr$1$kI,MW>e!“e(B
e$BX?t$NBP?t$,e(B) e$B$”$l$P9bB.$K7W;;$G$-$k$N$K!"$H;DG0$J;W$$$r$7e(B
e$B$^$7$?!#e(B
e$B$H$$$&$o$1$G!"e(B
Math.gamma(x) => float
Math.lgamma(x) => [float, sign_of_gamma_x]
e$B$H$$$&$N$r<BAu$7$F$_$?$s$G$9$,$I$&$G$7$g$&$+!#e(B
e$BM_$7$+$C$?$N$Oe(B lgamma e$B$J$s$G$9$,!“e(Blgamma e$B$@$1$Ge(B gamma
e$B$,$J$$e(B
e$B$N$b$J$s$J$N$Ge(B gamma e$B$b<BAu$7$F$”$j$^$9!#e(B
e$B;HMQ$9$k4X?t$Oe(B tgamma e$B$He(B lgamma_r e$B$G$9!#e(B
tgamma e$B$Oe(B C99 e$B$GDj5A$5$l$F$^$9!#e(B
lgamma_r e$B$O5,3J$K$“$j$^$;$s!#5,3J$K$”$k$N$Oe(B lgamma e$B$G!"$3$le(B
e$B$O%9%l%C%I%;!<%U$8$c$J$$$N$GHr$1$F!“e(Blgamma_r
e$B$r;H$C$F$$$^$9!#e(B
lgamma_r e$B$Oe(B GNU/Linux, NetBSD, Solaris e$B$J$I$K$”$j$^$9!#e(B
missing/tgamma.c e$B$O!Ve(BCe$B8@8l$K$h$k:G?7%"%k%4%j%:%`;vE5!W$Ne(B
e$B%5%]!<%H%Z!<%8$Ne(B 『C言語による最新アルゴリズム事典』
e$B$+$i<hF@$7$?e(B algo.tar.gz e$BFb$Ne(B gamma.c
e$B$r85$K$7$F4X?tL>$d%3%ae(B
e$B%s%H$J$I$rD4@0$7$?$b$N$G$9!#e(B
missing/lgamma_r.c e$B$Oe(B tgamma.c e$B$r$5$i$K$$$8$C$F:n$j$^$7$?!#e(B
algo.tar.gz e$BCf$Ne(B README e$B$Ke(B
| e$B$4MxMQ$K$D$$$F$N@)8B$O$“$j$^$;$s!#$?$@$7%P%0$K$h$kB;32$Ne(B
| e$BGe=~$J$I$K$O1~$8$+$M$^$9$N$G$4N;>5$/$@$5$$!#e(B
e$B$H$”$C$?$N$H!“$d$O$j!Ve(BCe$B8@8l$K$h$k:G?7%”%k%4%j%:%`;vE5!W$+$ie(B
e$B$H$i$l$F$$$ke(B erf.c e$B$Ke(B public domain
e$B$H=q$$$F$"$C$?$N$G!“e(B
missing/tgamma.c, missing/lgamma_r.c e$B$be(B public domain e$B$H$7$Fe(B
e$B$”$j$^$9!#e(B
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
Mathmodule contains module functions for basic -
trigonometric and transcendental functions. See class
-
Floatfor 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 A…
+
+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: