diff options
Diffstat (limited to 'mrbgems/mruby-math/src/math.c')
| -rw-r--r-- | mrbgems/mruby-math/src/math.c | 127 |
1 files changed, 96 insertions, 31 deletions
diff --git a/mrbgems/mruby-math/src/math.c b/mrbgems/mruby-math/src/math.c index 871cba301..79fdb7c6f 100644 --- a/mrbgems/mruby-math/src/math.c +++ b/mrbgems/mruby-math/src/math.c @@ -4,29 +4,94 @@ ** See Copyright Notice in mruby.h */ -#include "mruby.h" -#include "mruby/array.h" +#include <mruby.h> -#include <math.h> +#ifdef MRB_NO_FLOAT +# error Math conflicts with 'MRB_NO_FLOAT' configuration +#endif + +#include <mruby/array.h> +#include <mruby/presym.h> static void domain_error(mrb_state *mrb, const char *func) { - struct RClass *math = mrb_module_get(mrb, "Math"); - struct RClass *domainerror = mrb_class_get_under(mrb, math, "DomainError"); - mrb_value str = mrb_str_new_cstr(mrb, func); - mrb_raisef(mrb, domainerror, "Numerical argument is out of domain - %S", str); + struct RClass *math = mrb_module_get_id(mrb, MRB_SYM(Math)); + struct RClass *domainerror = mrb_class_get_under_id(mrb, math, MRB_SYM(DomainError)); + mrb_raisef(mrb, domainerror, "Numerical argument is out of domain - %s", func); } /* math functions not provided by Microsoft Visual C++ 2012 or older */ -#if defined _MSC_VER && _MSC_VER < 1800 +#if defined _MSC_VER && _MSC_VER <= 1700 -#define MATH_TOLERANCE 1E-12 +double +asinh(double x) +{ + double xa, ya, y; -#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) + /* Basic formula loses precision for x < 0, but asinh is an odd function */ + xa = fabs(x); + if (xa > 3.16227E+18) { + /* Prevent x*x from overflowing; basic formula reduces to log(2*x) */ + ya = log(xa) + 0.69314718055994530942; + } + else { + /* Basic formula for asinh */ + ya = log(xa + sqrt(xa*xa + 1.0)); + } + + y = _copysign(ya, x); + return y; +} + +double +acosh(double x) +{ + double y; + + if (x > 3.16227E+18) { + /* Prevent x*x from overflowing; basic formula reduces to log(2*x) */ + y = log(x) + 0.69314718055994530942; + } + else { + /* Basic formula for acosh */ + y = log(x + sqrt(x*x - 1.0)); + } + + return y; +} + +double +atanh(double x) +{ + double y; + + if (fabs(x) < 1E-2) { + /* The sums 1+x and 1-x lose precision for small x. Use the polynomial + instead. */ + double x2 = x * x; + y = x*(1.0 + x2*(1.0/3.0 + x2*(1.0/5.0 + x2*(1.0/7.0)))); + } + else { + /* Basic formula for atanh */ + y = 0.5 * (log(1.0+x) - log(1.0-x)); + } + + return y; +} + +double +cbrt(double x) +{ + double xa, ya, y; + + /* pow(x, y) is undefined for x < 0 and y not an integer, but cbrt is an + odd function */ + xa = fabs(x); + ya = pow(xa, 1.0/3.0); + y = _copysign(ya, x); + return y; +} /* Declaration of complementary Error function */ double @@ -56,7 +121,8 @@ erf(double x) term *= xsqr/j; sum += term/(2*j+1); ++j; - } while (fabs(term/sum) > MATH_TOLERANCE); + if (sum == 0) break; + } while (fabs(term/sum) > DBL_EPSILON); return two_sqrtpi*sum; } @@ -89,12 +155,16 @@ erfc(double x) n += 0.5; q1 = q2; q2 = b/d; - } while (fabs(q1-q2)/q2 > MATH_TOLERANCE); + } while (fabs(q1-q2)/q2 > DBL_EPSILON); return one_sqrtpi*exp(-x*x)*q2; } #endif +#if defined __FreeBSD__ && !defined __FreeBSD_version +#include <osreldate.h> /* for __FreeBSD_version */ +#endif + #if (defined _MSC_VER && _MSC_VER < 1800) || defined __ANDROID__ || (defined __FreeBSD__ && __FreeBSD_version < 803000) double @@ -170,7 +240,8 @@ math_tan(mrb_state *mrb, mrb_value obj) * call-seq: * Math.asin(x) -> float * - * Computes the arc sine of <i>x</i>. Returns -{PI/2} .. {PI/2}. + * Computes the arc sine of <i>x</i>. + * @return computed value between `-(PI/2)` and `(PI/2)`. */ static mrb_value math_asin(mrb_state *mrb, mrb_value obj) @@ -210,7 +281,7 @@ math_acos(mrb_state *mrb, mrb_value obj) * call-seq: * Math.atan(x) -> float * - * Computes the arc tangent of <i>x</i>. Returns -{PI/2} .. {PI/2}. + * Computes the arc tangent of <i>x</i>. Returns `-(PI/2) .. (PI/2)`. */ static mrb_value math_atan(mrb_state *mrb, mrb_value obj) @@ -419,7 +490,7 @@ static mrb_value math_log(mrb_state *mrb, mrb_value obj) { mrb_float x, base; - int argc; + mrb_int argc; argc = mrb_get_args(mrb, "f|f", &x, &base); if (x < 0.0) { @@ -556,7 +627,7 @@ math_cbrt(mrb_state *mrb, mrb_value obj) * Math.frexp(numeric) -> [ fraction, exponent ] * * Returns a two-element array containing the normalized fraction (a - * <code>Float</code>) and exponent (a <code>Fixnum</code>) of + * <code>Float</code>) and exponent (a <code>Integer</code>) of * <i>numeric</i>. * * fraction, exponent = Math.frexp(1234) #=> [0.6025390625, 11] @@ -590,7 +661,7 @@ math_ldexp(mrb_state *mrb, mrb_value obj) mrb_int i; mrb_get_args(mrb, "fi", &x, &i); - x = ldexp(x, i); + x = ldexp(x, (int)i); return mrb_float_value(mrb, x); } @@ -657,24 +728,18 @@ mrb_mruby_math_gem_init(mrb_state* mrb) struct RClass *mrb_math; mrb_math = mrb_define_module(mrb, "Math"); - mrb_define_class_under(mrb, mrb_math, "DomainError", mrb->eStandardError_class); + mrb_define_class_under_id(mrb, mrb_math, MRB_SYM(DomainError), mrb->eStandardError_class); #ifdef M_PI - mrb_define_const(mrb, mrb_math, "PI", mrb_float_value(mrb, M_PI)); + mrb_define_const_id(mrb, mrb_math, MRB_SYM(PI), mrb_float_value(mrb, M_PI)); #else - mrb_define_const(mrb, mrb_math, "PI", mrb_float_value(mrb, atan(1.0)*4.0)); + mrb_define_const_id(mrb, mrb_math, MRB_SYM(PI), mrb_float_value(mrb, atan(1.0)*4.0)); #endif #ifdef M_E - mrb_define_const(mrb, mrb_math, "E", mrb_float_value(mrb, M_E)); -#else - mrb_define_const(mrb, mrb_math, "E", mrb_float_value(mrb, exp(1.0))); -#endif - -#ifdef MRB_USE_FLOAT - mrb_define_const(mrb, mrb_math, "TOLERANCE", mrb_float_value(mrb, 1e-5)); + mrb_define_const_id(mrb, mrb_math, MRB_SYM(E), mrb_float_value(mrb, M_E)); #else - mrb_define_const(mrb, mrb_math, "TOLERANCE", mrb_float_value(mrb, 1e-12)); + mrb_define_const_id(mrb, mrb_math, MRB_SYM(E), mrb_float_value(mrb, exp(1.0))); #endif mrb_define_module_function(mrb, mrb_math, "sin", math_sin, MRB_ARGS_REQ(1)); |