From 04815be9cdf3216ce925122c82d8de6db41ddded Mon Sep 17 00:00:00 2001 From: Paolo Bosetti Date: Tue, 15 May 2012 17:41:07 -0700 Subject: Re-added Math finctions not supported under MSVC, and added preprocessor directives to support these functions under MSVC too --- src/math.c | 179 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) (limited to 'src/math.c') diff --git a/src/math.c b/src/math.c index adf02206c..2d52faaca 100644 --- a/src/math.c +++ b/src/math.c @@ -10,6 +10,39 @@ #define domain_error(msg) \ mrb_raise(mrb, E_RANGE_ERROR, "Numerical argument is out of domain - " #msg); +/* math functions not provided under Microsoft Visual C++ */ +#ifdef _MSC_VER +#define asinh(x) log(x + sqrt(pow(x,2.0) + 1)) +#define acosh(x) log(x + sqrt(pow(x,2.0) - 1)) +#define atanh(x) (log(1+x) - log(1-x))/2.0 +#define cbrt(x) pow(x,1.0/3.0) + + +double erf(double x) +{ + // constants + double a1 = 0.254829592; + double a2 = -0.284496736; + double a3 = 1.421413741; + double a4 = -1.453152027; + double a5 = 1.061405429; + double p = 0.3275911; + + // Save the sign of x + int sign = 1; + if (x < 0) + sign = -1; + x = fabs(x); + + // A&S formula 7.1.26 + double t = 1.0/(1.0 + p*x); + double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x); + + return sign*y; +} + +#endif + mrb_value mrb_assoc_new(mrb_state *mrb, mrb_value car, mrb_value cdr); @@ -215,6 +248,62 @@ math_tanh(mrb_state *mrb, mrb_value obj) } +/* + INVERSE HYPERBOLIC TRIG FUNCTIONS +*/ + +/* + * call-seq: + * Math.asinh(x) -> float + * + * Computes the inverse hyperbolic sine of x. + */ +static mrb_value +math_asinh(mrb_state *mrb, mrb_value obj) +{ + mrb_float x; + + mrb_get_args(mrb, "f", &x); + + x = asinh(x); + + return mrb_float_value(x); +} + +/* + * call-seq: + * Math.acosh(x) -> float + * + * Computes the inverse hyperbolic cosine of x. + */ +static mrb_value +math_acosh(mrb_state *mrb, mrb_value obj) +{ + mrb_float x; + + mrb_get_args(mrb, "f", &x); + x = acosh(x); + + return mrb_float_value(x); +} + +/* + * call-seq: + * Math.atanh(x) -> float + * + * Computes the inverse hyperbolic tangent of x. + */ +static mrb_value +math_atanh(mrb_state *mrb, mrb_value obj) +{ + mrb_float x; + + mrb_get_args(mrb, "f", &x); + x = atanh(x); + + return mrb_float_value(x); +} + /* EXPONENTIALS AND LOGARITHMS */ @@ -332,6 +421,49 @@ math_log10(mrb_state *mrb, mrb_value obj) return mrb_float_value(x); } +/* + * call-seq: + * Math.cbrt(numeric) -> float + * + * Returns the cube root of numeric. + * + * -9.upto(9) {|x| + * p [x, Math.cbrt(x), Math.cbrt(x)**3] + * } + * #=> + * [-9, -2.0800838230519, -9.0] + * [-8, -2.0, -8.0] + * [-7, -1.91293118277239, -7.0] + * [-6, -1.81712059283214, -6.0] + * [-5, -1.7099759466767, -5.0] + * [-4, -1.5874010519682, -4.0] + * [-3, -1.44224957030741, -3.0] + * [-2, -1.25992104989487, -2.0] + * [-1, -1.0, -1.0] + * [0, 0.0, 0.0] + * [1, 1.0, 1.0] + * [2, 1.25992104989487, 2.0] + * [3, 1.44224957030741, 3.0] + * [4, 1.5874010519682, 4.0] + * [5, 1.7099759466767, 5.0] + * [6, 1.81712059283214, 6.0] + * [7, 1.91293118277239, 7.0] + * [8, 2.0, 8.0] + * [9, 2.0800838230519, 9.0] + * + */ +static mrb_value +math_cbrt(mrb_state *mrb, mrb_value obj) +{ + mrb_float x; + + mrb_get_args(mrb, "f", &x); + x = cbrt(x); + + return mrb_float_value(x); +} + + /* * call-seq: * Math.frexp(numeric) -> [ fraction, exponent ] @@ -396,6 +528,43 @@ math_hypot(mrb_state *mrb, mrb_value obj) return mrb_float_value(x); } +/* + * call-seq: + * Math.erf(x) -> float + * + * Calculates the error function of x. + */ +static mrb_value +math_erf(mrb_state *mrb, mrb_value obj) +{ + mrb_float x; + + mrb_get_args(mrb, "f", &x); + x = erf(x); + + return mrb_float_value(x); +} + + +#ifndef _MSC_VER +/* + * call-seq: + * Math.erfc(x) -> float + * + * Calculates the complementary error function of x. + */ +static mrb_value +math_erfc(mrb_state *mrb, mrb_value obj) +{ + mrb_float x; + + mrb_get_args(mrb, "f", &x); + x = erfc(x); + + return mrb_float_value(x); +} +#endif + /* ------------------------------------------------------------------------*/ void mrb_init_math(mrb_state *mrb) @@ -428,13 +597,23 @@ mrb_init_math(mrb_state *mrb) mrb_define_class_method(mrb, mrb_math, "cosh", math_cosh, 1); mrb_define_class_method(mrb, mrb_math, "tanh", math_tanh, 1); + mrb_define_class_method(mrb, mrb_math, "asinh", math_asinh, 1); + mrb_define_class_method(mrb, mrb_math, "acosh", math_acosh, 1); + mrb_define_class_method(mrb, mrb_math, "atanh", math_atanh, 1); + mrb_define_class_method(mrb, mrb_math, "exp", math_exp, 1); mrb_define_class_method(mrb, mrb_math, "log", math_log, -1); mrb_define_class_method(mrb, mrb_math, "log2", math_log2, 1); mrb_define_class_method(mrb, mrb_math, "log10", math_log10, 1); + mrb_define_class_method(mrb, mrb_math, "cbrt", math_cbrt, 1); mrb_define_class_method(mrb, mrb_math, "frexp", math_frexp, 1); mrb_define_class_method(mrb, mrb_math, "ldexp", math_ldexp, 2); mrb_define_class_method(mrb, mrb_math, "hypot", math_hypot, 2); + + mrb_define_class_method(mrb, mrb_math, "erf", math_erf, 1); +#ifndef _MSC_VER + mrb_define_class_method(mrb, mrb_math, "erfc", math_erfc, 1); +#endif } -- cgit v1.2.3 From 8d0738e7b4d6260c3b357c1fe37db22fed250337 Mon Sep 17 00:00:00 2001 From: Paolo Bosetti Date: Tue, 15 May 2012 17:58:40 -0700 Subject: Added also support for erfc under MSVC --- src/math.c | 71 +++++++++++++++++++++++++++++++++++++++++--------------------- 1 file changed, 47 insertions(+), 24 deletions(-) (limited to 'src/math.c') diff --git a/src/math.c b/src/math.c index 2d52faaca..56cde6bf5 100644 --- a/src/math.c +++ b/src/math.c @@ -17,28 +17,55 @@ #define atanh(x) (log(1+x) - log(1-x))/2.0 #define cbrt(x) pow(x,1.0/3.0) +#define REL_ERROR 1E-12 +/* Implementation of Error function */ +double +erf(double x) +{ + static const double two_sqrtpi= 1.128379167095512574; + if (fabs(x) > 2.2) { + return 1.0 - erfc(x); + } + double sum= x, term= x, xsqr= x*x; + int j= 1; + do { + term*= xsqr/j; + sum-= term/(2*j+1); + ++j; + term*= xsqr/j; + sum+= term/(2*j+1); + ++j; + } while (fabs(term)/sum > REL_ERROR); + return two_sqrtpi*sum; +} -double erf(double x) +/* Implementation of complementary Error function */ +double +erfc(double x) { - // constants - double a1 = 0.254829592; - double a2 = -0.284496736; - double a3 = 1.421413741; - double a4 = -1.453152027; - double a5 = 1.061405429; - double p = 0.3275911; - - // Save the sign of x - int sign = 1; - if (x < 0) - sign = -1; - x = fabs(x); - - // A&S formula 7.1.26 - double t = 1.0/(1.0 + p*x); - double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x); - - return sign*y; + static const double one_sqrtpi= 0.564189583547756287; + if (fabs(x) < 2.2) { + return 1.0 - erf(x); + } + if (signbit(x)) { + return 2.0 - erfc(-x); + } + double a=1, b=x; + double c=x, d=x*x+0.5; + double q1,q2; + double n= 1.0, t; + do { + t= a*n+b*x; + a= b; + b= t; + t= c*n+d*x; + c= d; + d= t; + n+= 0.5; + q1= q2; + q2= b/d; + } while (fabs(q1-q2)/q2 > REL_ERROR); + return one_sqrtpi*exp(-x*x)*q2; } #endif @@ -546,7 +573,6 @@ math_erf(mrb_state *mrb, mrb_value obj) } -#ifndef _MSC_VER /* * call-seq: * Math.erfc(x) -> float @@ -563,7 +589,6 @@ math_erfc(mrb_state *mrb, mrb_value obj) return mrb_float_value(x); } -#endif /* ------------------------------------------------------------------------*/ void @@ -613,7 +638,5 @@ mrb_init_math(mrb_state *mrb) mrb_define_class_method(mrb, mrb_math, "hypot", math_hypot, 2); mrb_define_class_method(mrb, mrb_math, "erf", math_erf, 1); -#ifndef _MSC_VER mrb_define_class_method(mrb, mrb_math, "erfc", math_erfc, 1); -#endif } -- cgit v1.2.3 From 6971b0e9ab025ddf0e5cfa16ab79ef0acc14acad Mon Sep 17 00:00:00 2001 From: Paolo Bosetti Date: Tue, 15 May 2012 18:04:40 -0700 Subject: Code formatting; Added credits for erf/erfc implementation --- src/math.c | 78 +++++++++++++++++++++++++++++++++----------------------------- 1 file changed, 41 insertions(+), 37 deletions(-) (limited to 'src/math.c') diff --git a/src/math.c b/src/math.c index 56cde6bf5..c6f412055 100644 --- a/src/math.c +++ b/src/math.c @@ -17,55 +17,59 @@ #define atanh(x) (log(1+x) - log(1-x))/2.0 #define cbrt(x) pow(x,1.0/3.0) +/* +** Implementations of error functions +** credits to http://www.digitalmars.com/archives/cplusplus/3634.html +*/ #define REL_ERROR 1E-12 /* Implementation of Error function */ double erf(double x) { static const double two_sqrtpi= 1.128379167095512574; - if (fabs(x) > 2.2) { - return 1.0 - erfc(x); - } - double sum= x, term= x, xsqr= x*x; - int j= 1; - do { - term*= xsqr/j; - sum-= term/(2*j+1); - ++j; - term*= xsqr/j; - sum+= term/(2*j+1); - ++j; - } while (fabs(term)/sum > REL_ERROR); - return two_sqrtpi*sum; + if (fabs(x) > 2.2) { + return 1.0 - erfc(x); + } + double sum= x, term= x, xsqr= x*x; + int j= 1; + do { + term*= xsqr/j; + sum-= term/(2*j+1); + ++j; + term*= xsqr/j; + sum+= term/(2*j+1); + ++j; + } while (fabs(term)/sum > REL_ERROR); + return two_sqrtpi*sum; } /* Implementation of complementary Error function */ double erfc(double x) { - static const double one_sqrtpi= 0.564189583547756287; - if (fabs(x) < 2.2) { - return 1.0 - erf(x); - } - if (signbit(x)) { - return 2.0 - erfc(-x); - } - double a=1, b=x; - double c=x, d=x*x+0.5; - double q1,q2; - double n= 1.0, t; - do { - t= a*n+b*x; - a= b; - b= t; - t= c*n+d*x; - c= d; - d= t; - n+= 0.5; - q1= q2; - q2= b/d; - } while (fabs(q1-q2)/q2 > REL_ERROR); - return one_sqrtpi*exp(-x*x)*q2; + static const double one_sqrtpi= 0.564189583547756287; + if (fabs(x) < 2.2) { + return 1.0 - erf(x); + } + if (signbit(x)) { + return 2.0 - erfc(-x); + } + double a=1, b=x; + double c=x, d=x*x+0.5; + double q1,q2; + double n= 1.0, t; + do { + t= a*n+b*x; + a= b; + b= t; + t= c*n+d*x; + c= d; + d= t; + n+= 0.5; + q1= q2; + q2= b/d; + } while (fabs(q1-q2)/q2 > REL_ERROR); + return one_sqrtpi*exp(-x*x)*q2; } #endif -- cgit v1.2.3 From 36aaaafae4437dda5664c99b3497976fecb642fa Mon Sep 17 00:00:00 2001 From: nkshigeru Date: Wed, 16 May 2012 15:14:23 +0900 Subject: fix for MSVC move declaration in the middle no signbit --- src/math.c | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) (limited to 'src/math.c') diff --git a/src/math.c b/src/math.c index c6f412055..a4451ea9d 100644 --- a/src/math.c +++ b/src/math.c @@ -17,6 +17,10 @@ #define atanh(x) (log(1+x) - log(1-x))/2.0 #define cbrt(x) pow(x,1.0/3.0) +/* Declaration of complementary Error function */ +double +erfc(double x); + /* ** Implementations of error functions ** credits to http://www.digitalmars.com/archives/cplusplus/3634.html @@ -27,11 +31,11 @@ double erf(double x) { static const double two_sqrtpi= 1.128379167095512574; + double sum= x, term= x, xsqr= x*x; + int j= 1; if (fabs(x) > 2.2) { return 1.0 - erfc(x); } - double sum= x, term= x, xsqr= x*x; - int j= 1; do { term*= xsqr/j; sum-= term/(2*j+1); @@ -48,16 +52,16 @@ double erfc(double x) { static const double one_sqrtpi= 0.564189583547756287; + double a=1, b=x; + double c=x, d=x*x+0.5; + double q1,q2; + double n= 1.0, t; if (fabs(x) < 2.2) { return 1.0 - erf(x); } - if (signbit(x)) { + if (x < 0.0) { /*signbit(x)*/ return 2.0 - erfc(-x); } - double a=1, b=x; - double c=x, d=x*x+0.5; - double q1,q2; - double n= 1.0, t; do { t= a*n+b*x; a= b; -- cgit v1.2.3 From b9b20c8be8edba01f72c378e3e626e93408d8089 Mon Sep 17 00:00:00 2001 From: Paolo Bosetti Date: Tue, 15 May 2012 23:55:44 -0700 Subject: Tested under MSVC by nkshigeru. Minor code formatting --- src/math.c | 44 +++++++++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 19 deletions(-) (limited to 'src/math.c') diff --git a/src/math.c b/src/math.c index a4451ea9d..96f803580 100644 --- a/src/math.c +++ b/src/math.c @@ -11,6 +11,7 @@ mrb_raise(mrb, E_RANGE_ERROR, "Numerical argument is out of domain - " #msg); /* math functions not provided under Microsoft Visual C++ */ +#define _MSC_VER #ifdef _MSC_VER #define asinh(x) log(x + sqrt(pow(x,2.0) + 1)) #define acosh(x) log(x + sqrt(pow(x,2.0) - 1)) @@ -30,18 +31,20 @@ erfc(double x); double erf(double x) { - static const double two_sqrtpi= 1.128379167095512574; - double sum= x, term= x, xsqr= x*x; + static const double two_sqrtpi = 1.128379167095512574; + double sum = x; + double term = x; + double xsqr = x*x; int j= 1; if (fabs(x) > 2.2) { return 1.0 - erfc(x); } do { - term*= xsqr/j; - sum-= term/(2*j+1); + term *= xsqr/j; + sum -= term/(2*j+1); ++j; - term*= xsqr/j; - sum+= term/(2*j+1); + term *= xsqr/j; + sum += term/(2*j+1); ++j; } while (fabs(term)/sum > REL_ERROR); return two_sqrtpi*sum; @@ -52,10 +55,13 @@ double erfc(double x) { static const double one_sqrtpi= 0.564189583547756287; - double a=1, b=x; - double c=x, d=x*x+0.5; - double q1,q2; - double n= 1.0, t; + double a = 1; + double b = x; + double c = x; + double d = x*x+0.5; + double q1, q2; + double n = 1.0; + double t; if (fabs(x) < 2.2) { return 1.0 - erf(x); } @@ -63,15 +69,15 @@ erfc(double x) return 2.0 - erfc(-x); } do { - t= a*n+b*x; - a= b; - b= t; - t= c*n+d*x; - c= d; - d= t; - n+= 0.5; - q1= q2; - q2= b/d; + t = a*n+b*x; + a = b; + b = t; + t = c*n+d*x; + c = d; + d = t; + n += 0.5; + q1 = q2; + q2 = b/d; } while (fabs(q1-q2)/q2 > REL_ERROR); return one_sqrtpi*exp(-x*x)*q2; } -- cgit v1.2.3