西松です.

[ruby-list:37740] でまつもとさん:
(erf, erfc について)
>「きりがないから」というのがその理由です。で、一応POSIXをひ
>とつの目安にしてます。が、C99に入ってるなら問題ないでしょう。
>
>|まずその行は実行されないハズで, 実行されたとしても最小
>|の桁が浮動小数点数演算のため振動した結果だと思うので,
>|[1]その行は削除.
>|[2]printf("missing/erf.c:q_gamma(): could not converge.");
>|      などにしておく.
>|もしくは,
>|[3]http://ftp.rutgers.edu/pub/FreeBSD/FreeBSD-current/src/lib/libm/common_source/erf.c
>|      のような近似曲線を使ったものにする.(これをその
>|      まま使うのはライセンス的にまずいのですよね.)
>|という策があります.
>
>どれが良いのかは、私がerfを使わないので分かりません。
>
>少なくとも[3]のものにライセンス的な問題はありません、BSDライ
>センスなら。FreeBSDの一部なのできっとそうでしょう。


結局, [2]を選択して, [ruby-list:37737] の missing/erf.c を
    fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
などと変更しました. これは実行されることのまずない保険です.

[3]にしなかったのは, いろいろ変更を試みましたが, いろいろな
OSでコンパイルできるか自信がなかったからです.

     ちなみに, 教えをあおいだP氏によれば, 現在の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
     でした. こいつは下から5行目の
     	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 <stdio.h>
+#include <math.h>
+
+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