diff --git a/ChangeLog b/ChangeLog index 4224fbcc0c..329cdc0b13 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,15 @@ Thu Jun 5 18:33:46 2003 WATANABE Hirofumi * ext/curses/curses.c (window_s_allocate,curses_finalize): avoid VC++ warnings. +Thu Jun 5 16:11:50 2003 NAKAMURA Usaku + + * 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 * math.c (math_erf,math_erfc): new function. [ruby-list:37753] diff --git a/bcc32/Makefile.sub b/bcc32/Makefile.sub index 695422b382..dda52b9d7b 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 5416500cb0..0603db1a40 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 0000000000..e30a90051e --- /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 +#include + +#ifdef _WIN32 +# include +# 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 8b00c72fde..8e81404190 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 d4b058bb82..d56b217533 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