summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ChangeLog9
-rw-r--r--bcc32/Makefile.sub5
-rw-r--r--missing.h4
-rw-r--r--missing/erf.c91
-rw-r--r--win32/Makefile.sub5
-rw-r--r--wince/Makefile.sub5
6 files changed, 111 insertions, 8 deletions
diff --git a/ChangeLog b/ChangeLog
index 4224fbc..329cdc0 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -3,6 +3,15 @@ Thu Jun 5 18:33:46 2003 WATANABE Hirofumi <eban@ruby-lang.org>
* ext/curses/curses.c (window_s_allocate,curses_finalize):
avoid VC++ warnings.
+Thu Jun 5 16:11:50 2003 NAKAMURA Usaku <usa@ruby-lang.org>
+
+ * bcc32/Makefile.sub, win32/Makefile.sub, wince/Makefile.sub
+ (MISSING): link with missing/erf.c.
+
+ * missing.h (erf, erfc): fix prototype.
+
+ * missing/erf.c: new. [ruby-list:37753]
+
Thu Jun 5 15:09:06 2003 Yukihiro Matsumoto <matz@ruby-lang.org>
* math.c (math_erf,math_erfc): new function. [ruby-list:37753]
diff --git a/bcc32/Makefile.sub b/bcc32/Makefile.sub
index 695422b..dda52b9 100644
--- a/bcc32/Makefile.sub
+++ b/bcc32/Makefile.sub
@@ -102,7 +102,7 @@ RFLAGS = $(iconinc)
EXTLIBS =
!endif
LIBS = cw32.lib import32.lib ws2_32.lib $(EXTLIBS)
-MISSING = acosh.obj crypt.obj win32.obj
+MISSING = acosh.obj crypt.obj erf.obj win32.obj
!ifndef STACK
STACK = 0x2000000
@@ -334,7 +334,7 @@ s,@AR@,$(AR),;t t
s,@ARFLAGS@,$(ARFLAGS) ,;t t
s,@LN_S@,$(LN_S),;t t
s,@SET_MAKE@,$(SET_MAKE),;t t
-s,@LIBOBJS@, acosh.obj crypt.obj win32.obj,;t t
+s,@LIBOBJS@, acosh.obj crypt.obj erf.obj win32.obj,;t t
s,@ALLOCA@,$(ALLOCA),;t t
s,@DEFAULT_KCODE@,$(DEFAULT_KCODE),;t t
s,@EXEEXT@,.exe,;t t
@@ -518,6 +518,7 @@ acosh.obj: acosh.c win32.h
alloca.obj: alloca.c win32.h
crypt.obj: crypt.c win32.h
dup2.obj: dup2.c win32.h
+erf.obj: erf.c win32.h
finite.obj: finite.c win32.h
flock.obj: flock.c win32.h
memcmp.obj: memcmp.c win32.h
diff --git a/missing.h b/missing.h
index 5416500..0603db1 100644
--- a/missing.h
+++ b/missing.h
@@ -55,8 +55,8 @@ extern double hypot _((double, double));
#endif
#ifndef HAVE_ERF
-extern double erf _((double, double));
-extern double erfc _((double, double));
+extern double erf _((double));
+extern double erfc _((double));
#endif
#ifndef HAVE_ISINF
diff --git a/missing/erf.c b/missing/erf.c
new file mode 100644
index 0000000..e30a900
--- /dev/null
+++ b/missing/erf.c
@@ -0,0 +1,91 @@
+/* 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>
+
+#ifdef _WIN32
+# include <float.h>
+# if !defined __MINGW32__ || defined __NO_ISOCEXT
+# ifndef isnan
+# define isnan(x) _isnan(x)
+# endif
+# ifndef isinf
+# define isinf(x) (!_finite(x) && !_isnan(x))
+# endif
+# ifndef finite
+# define finite(x) _finite(x)
+# endif
+# endif
+#endif
+
+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 --git a/win32/Makefile.sub b/win32/Makefile.sub
index 8b00c72..8e81404 100644
--- a/win32/Makefile.sub
+++ b/win32/Makefile.sub
@@ -93,7 +93,7 @@ RFLAGS = -r
EXTLIBS =
!endif
LIBS = oldnames.lib user32.lib advapi32.lib wsock32.lib $(EXTLIBS)
-MISSING = acosh.obj crypt.obj win32.obj
+MISSING = acosh.obj crypt.obj erf.obj win32.obj
ARFLAGS = -machine:$(MACHINE) -out:
CC = $(CC) -nologo
@@ -331,7 +331,7 @@ s,@AR@,$(AR),;t t
s,@ARFLAGS@,$(ARFLAGS),;t t
s,@LN_S@,$(LN_S),;t t
s,@SET_MAKE@,$(SET_MAKE),;t t
-s,@LIBOBJS@, acosh.obj crypt.obj win32.obj,;t t
+s,@LIBOBJS@, acosh.obj crypt.obj erf.obj win32.obj,;t t
s,@ALLOCA@,$(ALLOCA),;t t
s,@DEFAULT_KCODE@,$(DEFAULT_KCODE),;t t
s,@EXEEXT@,.exe,;t t
@@ -517,6 +517,7 @@ acosh.obj: {$(srcdir)}missing/acosh.c
alloca.obj: {$(srcdir)}missing/alloca.c
crypt.obj: {$(srcdir)}missing/crypt.c
dup2.obj: {$(srcdir)}missing/dup2.c
+erf.obj: {$(srcdir)}missing/erf.c
finite.obj: {$(srcdir)}missing/finite.c
flock.obj: {$(srcdir)}missing/flock.c
memcmp.obj: {$(srcdir)}missing/memcmp.c
diff --git a/wince/Makefile.sub b/wince/Makefile.sub
index d4b058b..d56b217 100644
--- a/wince/Makefile.sub
+++ b/wince/Makefile.sub
@@ -97,7 +97,7 @@ RFLAGS = -r
EXTLIBS =
!endif
LIBS = coredll.lib winsock.lib $(EXTLIBS)
-MISSING = acosh.obj crypt.obj dup2.obj hypot.obj \
+MISSING = acosh.obj crypt.obj dup2.obj erf.obj hypot.obj \
isinf.obj isnan.obj strftime.obj win32.obj
WINCEOBJ= assert.obj direct.obj errno.obj io_wce.obj process_wce.obj \
signal_wce.obj stdio.obj stdlib.obj string_wce.obj \
@@ -361,7 +361,7 @@ s,@AR@,$(AR),;t t
s,@ARFLAGS@,$(ARFLAGS),;t t
s,@LN_S@,$(LN_S),;t t
s,@SET_MAKE@,$(SET_MAKE),;t t
-s,@LIBOBJS@, acosh.obj crypt.obj win32.obj isinf.obj isnan.obj,;t t
+s,@LIBOBJS@, acosh.obj crypt.obj erf.obj win32.obj isinf.obj isnan.obj,;t t
s,@ALLOCA@,$(ALLOCA),;t t
s,@DEFAULT_KCODE@,$(DEFAULT_KCODE),;t t
s,@EXEEXT@,.exe,;t t
@@ -552,6 +552,7 @@ acosh.obj: {$(srcdir)}missing/acosh.c
alloca.obj: {$(srcdir)}missing/alloca.c
crypt.obj: {$(srcdir)}missing/crypt.c
dup2.obj: {$(srcdir)}missing/dup2.c
+erf.obj: {$(srcdir)}missing/erf.c
finite.obj: {$(srcdir)}missing/finite.c
flock.obj: {$(srcdir)}missing/flock.c
isinf.obj: {$(srcdir)}missing/isinf.c