[ruby-list:37737] missing/erf.c
From:
NISHIMATSU Takeshi <t-nissie@...>
Date:
2003-05-30 12:06:11 UTC
List:
ruby-list #37737
西松と申します.
数学関数について <ruby-talk:37737> でまつもとさんは
> These functions and constants are defined under difference standards.
> POSIX, SVID, SVIF, etc. I'm not sure how far we should go, and how we
> treat undefined functions (yet).
と書かれていますが, missing/foo.c があればパッチを
受け付けていただけるかなと思い, とりあえず誤差関数の
erf, erfc のために missing/erf.c を書きました.
というか, ほとんど奥村晴彦著「Ｃ言語による最新
アルゴリズム辞典」からのコピー＆ペーストです.
ありがたいことに「プログラムは自由にお使いいただいて
かまわない.」と本にあるのでライセンスの問題はないです.
懸案の Bessel関数 (j0,j1,jn,y0,y1,yn) や lgamma, expm1,
log1p も missing/*.c があれば1.8.0に取り込んでもらえ
ますでしょうか. それとももう間に合わないでしょうか.
missing/*.c の精度や新しい数学関数の出す例外について
も議論が必要ですが.
今回の missing/erf.c の中には
rb_raise(rb_eRuntimeError, "missing/erf.c:q_gamma(): could not converge.");
があります. よいでしょうか.
2003-05-29のsnapshotのconfigure.in, math.c, missing.h
へのパッチと missing/erf.c を以下に添付します.
＝＝＝パッチ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
diff -r -u ../ruby-1.8.0-2003-05-29-snapshot.original/configure.in ./configure.in
--- ../ruby-1.8.0-2003-05-29-snapshot.original/configure.in Tue May 13 20:34:48 2003
+++ ./configure.in Fri May 30 00:54:03 2003
@@ -371,7 +371,7 @@
AC_CHECK_FUNCS(ftello)
AC_REPLACE_FUNCS(dup2 memmove mkdir strcasecmp strncasecmp strerror strftime\
strchr strstr strtoul crypt flock vsnprintf\
- isinf isnan finite hypot acosh)
+ isinf isnan finite hypot acosh erf)
AC_CHECK_FUNCS(fmod killpg wait4 waitpid syscall chroot fsync\
truncate chsize times utimes fcntl lockf lstat symlink readlink\
setitimer setruid seteuid setreuid setresuid setproctitle\
diff -r -u ../ruby-1.8.0-2003-05-29-snapshot.original/math.c ./math.c
--- ../ruby-1.8.0-2003-05-29-snapshot.original/math.c Thu Jan 16 16:34:01 2003
+++ ./math.c Fri May 30 00:13:48 2003
@@ -271,6 +271,22 @@
return rb_float_new(hypot(RFLOAT(x)->value, RFLOAT(y)->value));
}
+static VALUE
+math_erf(obj, x)
+ VALUE obj, x;
+{
+ Need_Float(x);
+ return rb_float_new(erf(RFLOAT(x)->value));
+}
+
+static VALUE
+math_erfc(obj, x)
+ VALUE obj, x;
+{
+ Need_Float(x);
+ return rb_float_new(erfc(RFLOAT(x)->value));
+}
+
void
Init_Math()
{
@@ -314,4 +330,7 @@
rb_define_module_function(rb_mMath, "ldexp", math_ldexp, 2);
rb_define_module_function(rb_mMath, "hypot", math_hypot, 2);
+
+ rb_define_module_function(rb_mMath, "erf", math_erf, 1);
+ rb_define_module_function(rb_mMath, "erfc", math_erfc, 1);
}
diff -r -u ../ruby-1.8.0-2003-05-29-snapshot.original/missing.h ./missing.h
--- ../ruby-1.8.0-2003-05-29-snapshot.original/missing.h Fri Mar 21 01:43:24 2003
+++ ./missing.h Thu May 29 23:59:34 2003
@@ -54,6 +54,11 @@
extern double hypot _((double, double));
#endif
+#ifndef HAVE_ERF
+extern double erf _((double, double));
+extern double erfc _((double, double));
+#endif
+
#ifndef HAVE_ISINF
extern int isinf _((double));
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝missing/erf.c＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
/* erf.c
reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
(New Algorithm handbook in C language) (Gijyutsu hyouron
sha, Tokyo, 1991) p.227 [in Japanese] */
#include "ruby.h"
#include <math.h>
static double q_gamma(double, double, double);
/* 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */
static double p_gamma(a, x, loggamma_a)
double a, x, loggamma_a;
{
int k;
double result, term, previous;
if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
if (x == 0) return 0;
result = term = exp(a * log(x) - x - loggamma_a) / a;
for (k = 1; k < 1000; k++) {
term *= x / (a + k);
previous = result; result += term;
if (result == previous) return result;
}
rb_raise(rb_eRuntimeError, "missing/erf.c:p_gamma(): could not converge.");
return result;
}
/* 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */
static double q_gamma(a, x, loggamma_a)
double a, x, loggamma_a;
{
int k;
double result, w, temp, previous;
double la = 1, lb = 1 + x - a; /* Laguerre polynomial */
if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
w = exp(a * log(x) - x - loggamma_a);
result = w / lb;
for (k = 2; k < 1000; k++) {
temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
la = lb; lb = temp;
w *= (k - 1 - a) / k;
temp = w / (la * lb);
previous = result; result += temp;
if (result == previous) return result;
}
rb_raise(rb_eRuntimeError, "missing/erf.c:q_gamma(): could not converge.");
return result;
}
#define LOG_PI 1.14472988584940017 /* log_e(PI) */
double erf(x)
double x;
{
if (!finite(x)) {
if (isnan(x)) return x; /* erf(NaN) = NaN */
return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */
}
if (x >= 0) return p_gamma(0.5, x * x, LOG_PI / 2);
else return - p_gamma(0.5, x * x, LOG_PI / 2);
}
double erfc(x)
double x;
{
if (!finite(x)) {
if (isnan(x)) return x; /* erfc(NaN) = NaN */
return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */
}
if (x >= 0) return q_gamma(0.5, x * x, LOG_PI / 2);
else return 1 + p_gamma(0.5, x * x, LOG_PI / 2);
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
--
love && peace && free_software
NISHIMATSU Takeshi t-nissie@imr.edu OR t-nissie@imr.tohoku.ac.jp
西松 毅