-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* 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] git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@3910 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
- Loading branch information
Showing
6 changed files
with
111 additions
and
8 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -3,6 +3,15 @@ Thu Jun 5 18:33:46 2003 WATANABE Hirofumi <[email protected]> | |
* ext/curses/curses.c (window_s_allocate,curses_finalize): | ||
avoid VC++ warnings. | ||
|
||
Thu Jun 5 16:11:50 2003 NAKAMURA Usaku <[email protected]> | ||
|
||
* 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 <[email protected]> | ||
|
||
* math.c (math_erf,math_erfc): new function. [ruby-list:37753] | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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); | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters