From: NISHIMATSU Takeshi Date: 2003-06-05T03:42:50+09:00 Subject: [ruby-list:37753] Re: missing/erf.c 西松です. [ruby-list:37740] でまつもとさん: （erf, erfc について） >「きりがないから」というのがその理由です。で、一応POSIXをひ >とつの目安にしてます。が、C99に入ってるなら問題ないでしょう。 > >|まずその行は実行されないハズで, 実行されたとしても最小 >|の桁が浮動小数点数演算のため振動した結果だと思うので, >|［１］その行は削除. >|［２］printf("missing/erf.c:q_gamma(): could not converge."); >| などにしておく. >|もしくは, >|［３］http://ftp.rutgers.edu/pub/FreeBSD/FreeBSD-current/src/lib/libm/common_source/erf.c >| のような近似曲線を使ったものにする.（これをその >| まま使うのはライセンス的にまずいのですよね.） >|という策があります. > >どれが良いのかは、私がerfを使わないので分かりません。 > >少なくとも[3]のものにライセンス的な問題はありません、BSDライ >センスなら。FreeBSDの一部なのできっとそうでしょう。 結局, ［２］を選択して, [ruby-list:37737] の missing/erf.c を fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__); などと変更しました. これは実行されることのまずない保険です. ［３］にしなかったのは, いろいろ変更を試みましたが, いろいろな ＯＳでコンパイルできるか自信がなかったからです. ちなみに, 教えをあおいだＰ氏によれば, 現在のFreeBSDの erf.c（やその他のlibmの数学関数のソース）は http://www.jp.freebsd.org/cgi/cvsweb.cgi/src/lib/msun/src/s_erf.c であり, ぼくが参照した erf.c は「屋根裏」に置いてある http://www.jp.freebsd.org/cgi/cvsweb.cgi/src/lib/libm/common_source/Attic/exp.c でした. こいつは下から５行目の r = __exp__D(s, y)/x; を r = __exp__D(s, y)/ax; と直さないと x<-1.25 で erfc(x)>2 になるというバグがあります. 2003-06-03の1.8.0のsnapshotのconfigure.in, math.c, missing.h, LEGAL, MANIFESTへのパッチと missing/erf.c (パッチに含めました) を最後に添付します. 採用していただけるととてもうれしいです. -- love && peace && free_software NISHIMATSU Takeshi t-nissie@imr.edu OR t-nissie@imr.tohoku.ac.jp 西松 毅 ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ diff -u -r -P ../ruby-1.8.0-2003-06-03-snapshot.original/LEGAL ./LEGAL --- ../ruby-1.8.0-2003-06-03-snapshot.original/LEGAL Fri Jul 26 15:12:38 2002 +++ ./LEGAL Wed Jun 4 19:21:55 2003 @@ -134,6 +134,7 @@ x68/*: missing/alloca.c: missing/dup2.c: +missing/erf.c: missing/finite.c: missing/hypot.c: missing/isinf.c: diff -u -r -P ../ruby-1.8.0-2003-06-03-snapshot.original/MANIFEST ./MANIFEST --- ../ruby-1.8.0-2003-06-03-snapshot.original/MANIFEST Wed May 21 17:48:02 2003 +++ ./MANIFEST Wed Jun 4 19:21:20 2003 @@ -262,6 +262,7 @@ missing/alloca.c missing/crypt.c missing/dup2.c +missing/erf.c missing/file.h missing/fileblocks.c missing/finite.c diff -u -r -P ../ruby-1.8.0-2003-06-03-snapshot.original/configure.in ./configure.in --- ../ruby-1.8.0-2003-06-03-snapshot.original/configure.in Sun Jun 1 21:58:56 2003 +++ ./configure.in Wed Jun 4 19:19:49 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 -u -r -P ../ruby-1.8.0-2003-06-03-snapshot.original/math.c ./math.c --- ../ruby-1.8.0-2003-06-03-snapshot.original/math.c Thu Jan 16 16:34:01 2003 +++ ./math.c Wed Jun 4 19:19:49 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 -u -r -P ../ruby-1.8.0-2003-06-03-snapshot.original/missing/erf.c ./missing/erf.c --- ../ruby-1.8.0-2003-06-03-snapshot.original/missing/erf.c Thu Jan 1 09:00:00 1970 +++ ./missing/erf.c Wed Jun 4 19:27:13 2003 @@ -0,0 +1,76 @@ +/* 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 +#include + +static double q_gamma(double, double, double); + +/* Incomplete gamma function + 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; + } + fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__); + return result; +} + +/* Incomplete gamma function + 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; + } + fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__); + return result; +} + +#define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */ + +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_OVER_2); + else return - p_gamma(0.5, x * x, LOG_PI_OVER_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_OVER_2); + else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2); +} diff -u -r -P ../ruby-1.8.0-2003-06-03-snapshot.original/missing.h ./missing.h --- ../ruby-1.8.0-2003-06-03-snapshot.original/missing.h Fri Mar 21 01:43:24 2003 +++ ./missing.h Wed Jun 4 19:19:49 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