diff options
Diffstat (limited to 'math.c')
-rw-r--r-- | math.c | 148 |
1 files changed, 83 insertions, 65 deletions
@@ -65,23 +65,23 @@ math_atan2(VALUE unused_obj, VALUE y, VALUE x) dx = Get_Double(x); dy = Get_Double(y); if (dx == 0.0 && dy == 0.0) { - if (!signbit(dx)) - return DBL2NUM(dy); + if (!signbit(dx)) + return DBL2NUM(dy); if (!signbit(dy)) - return DBL2NUM(M_PI); - return DBL2NUM(-M_PI); + return DBL2NUM(M_PI); + return DBL2NUM(-M_PI); } #ifndef ATAN2_INF_C99 if (isinf(dx) && isinf(dy)) { - /* optimization for FLONUM */ - if (dx < 0.0) { - const double dz = (3.0 * M_PI / 4.0); - return (dy < 0.0) ? DBL2NUM(-dz) : DBL2NUM(dz); - } - else { - const double dz = (M_PI / 4.0); - return (dy < 0.0) ? DBL2NUM(-dz) : DBL2NUM(dz); - } + /* optimization for FLONUM */ + if (dx < 0.0) { + const double dz = (3.0 * M_PI / 4.0); + return (dy < 0.0) ? DBL2NUM(-dz) : DBL2NUM(dz); + } + else { + const double dz = (M_PI / 4.0); + return (dy < 0.0) ? DBL2NUM(-dz) : DBL2NUM(dz); + } } #endif return DBL2NUM(atan2(dy, dx)); @@ -170,6 +170,12 @@ math_tan(VALUE unused_obj, VALUE x) return DBL2NUM(tan(Get_Double(x))); } +#define math_arc(num, func) \ + double d; \ + d = Get_Double((num)); \ + domain_check_range(d, -1.0, 1.0, #func); \ + return DBL2NUM(func(d)); + /* * call-seq: * Math.acos(x) -> float @@ -190,11 +196,7 @@ math_tan(VALUE unused_obj, VALUE x) static VALUE math_acos(VALUE unused_obj, VALUE x) { - double d; - - d = Get_Double(x); - domain_check_range(d, -1.0, 1.0, "acos"); - return DBL2NUM(acos(d)); + math_arc(x, acos) } /* @@ -217,11 +219,7 @@ math_acos(VALUE unused_obj, VALUE x) static VALUE math_asin(VALUE unused_obj, VALUE x) { - double d; - - d = Get_Double(x); - domain_check_range(d, -1.0, 1.0, "asin"); - return DBL2NUM(asin(d)); + math_arc(x, asin) } /* @@ -476,7 +474,6 @@ math_exp(VALUE unused_obj, VALUE x) # define M_LN10 2.30258509299404568401799145468436421 #endif -static double math_log1(VALUE x); FUNC_MINIMIZED(static VALUE math_log(int, const VALUE *, VALUE)); /* @@ -511,20 +508,6 @@ math_log(int argc, const VALUE *argv, VALUE unused_obj) return rb_math_log(argc, argv); } -VALUE -rb_math_log(int argc, const VALUE *argv) -{ - VALUE x, base; - double d; - - rb_scan_args(argc, argv, "11", &x, &base); - d = math_log1(x); - if (argc == 2) { - d /= math_log1(base); - } - return DBL2NUM(d); -} - static double get_double_rshift(VALUE x, size_t *pnumbits) { @@ -536,23 +519,58 @@ get_double_rshift(VALUE x, size_t *pnumbits) x = rb_big_rshift(x, SIZET2NUM(numbits)); } else { - numbits = 0; + numbits = 0; } *pnumbits = numbits; return Get_Double(x); } static double -math_log1(VALUE x) +math_log_split(VALUE x, size_t *numbits) { - size_t numbits; - double d = get_double_rshift(x, &numbits); + double d = get_double_rshift(x, numbits); domain_check_min(d, 0.0, "log"); - /* check for pole error */ - if (d == 0.0) return -HUGE_VAL; + return d; +} + +#if defined(log2) || defined(HAVE_LOG2) +# define log_intermediate log2 +#else +# define log_intermediate log10 +double log2(double x); +#endif + +VALUE +rb_math_log(int argc, const VALUE *argv) +{ + VALUE x, base; + double d; + size_t numbits; - return log(d) + numbits * M_LN2; /* log(d * 2 ** numbits) */ + argc = rb_scan_args(argc, argv, "11", &x, &base); + d = math_log_split(x, &numbits); + if (argc == 2) { + size_t numbits_2; + double b = math_log_split(base, &numbits_2); + /* check for pole error */ + if (d == 0.0) { + // Already DomainError if b < 0.0 + return b ? DBL2NUM(-HUGE_VAL) : DBL2NUM(NAN); + } + else if (b == 0.0) { + return DBL2NUM(-0.0); + } + d = log_intermediate(d) / log_intermediate(b); + d += (numbits - numbits_2) / log2(b); + } + else { + /* check for pole error */ + if (d == 0.0) return DBL2NUM(-HUGE_VAL); + d = log(d); + d += numbits * M_LN2; + } + return DBL2NUM(d); } #ifndef log2 @@ -681,13 +699,13 @@ rb_math_sqrt(VALUE x) double d; if (RB_TYPE_P(x, T_COMPLEX)) { - VALUE neg = f_signbit(RCOMPLEX(x)->imag); - double re = Get_Double(RCOMPLEX(x)->real), im; - d = Get_Double(rb_complex_abs(x)); - im = sqrt((d - re) / 2.0); - re = sqrt((d + re) / 2.0); - if (neg) im = -im; - return rb_complex_new(DBL2NUM(re), DBL2NUM(im)); + VALUE neg = f_signbit(RCOMPLEX(x)->imag); + double re = Get_Double(RCOMPLEX(x)->real), im; + d = Get_Double(rb_complex_abs(x)); + im = sqrt((d - re) / 2.0); + re = sqrt((d + re) / 2.0); + if (neg) im = -im; + return rb_complex_new(DBL2NUM(re), DBL2NUM(im)); } d = Get_Double(x); domain_check_min(d, 0.0, "sqrt"); @@ -713,7 +731,7 @@ rb_math_sqrt(VALUE x) * cbrt(1.0) # => 1.0 * cbrt(0.0) # => 0.0 * cbrt(1.0) # => 1.0 - cbrt(2.0) # => 1.2599210498948732 + * cbrt(2.0) # => 1.2599210498948732 * cbrt(8.0) # => 2.0 * cbrt(27.0) # => 3.0 * cbrt(INFINITY) # => Infinity @@ -727,7 +745,7 @@ math_cbrt(VALUE unused_obj, VALUE x) double r = cbrt(f); #if defined __GLIBC__ if (isfinite(r) && !(f == 0.0 && r == 0.0)) { - r = (2.0 * r + (f / r / r)) / 3.0; + r = (2.0 * r + (f / r / r)) / 3.0; } #endif return DBL2NUM(r); @@ -945,17 +963,17 @@ math_gamma(VALUE unused_obj, VALUE x) d = Get_Double(x); /* check for domain error */ if (isinf(d)) { - if (signbit(d)) domain_error("gamma"); - return DBL2NUM(HUGE_VAL); + if (signbit(d)) domain_error("gamma"); + return DBL2NUM(HUGE_VAL); } if (d == 0.0) { - return signbit(d) ? DBL2NUM(-HUGE_VAL) : DBL2NUM(HUGE_VAL); + return signbit(d) ? DBL2NUM(-HUGE_VAL) : DBL2NUM(HUGE_VAL); } if (d == floor(d)) { - domain_check_min(d, 0.0, "gamma"); - if (1.0 <= d && d <= (double)NFACT_TABLE) { - return DBL2NUM(fact_table[(int)d - 1]); - } + domain_check_min(d, 0.0, "gamma"); + if (1.0 <= d && d <= (double)NFACT_TABLE) { + return DBL2NUM(fact_table[(int)d - 1]); + } } return DBL2NUM(tgamma(d)); } @@ -1007,12 +1025,12 @@ math_lgamma(VALUE unused_obj, VALUE x) d = Get_Double(x); /* check for domain error */ if (isinf(d)) { - if (signbit(d)) domain_error("lgamma"); - return rb_assoc_new(DBL2NUM(HUGE_VAL), INT2FIX(1)); + if (signbit(d)) domain_error("lgamma"); + return rb_assoc_new(DBL2NUM(HUGE_VAL), INT2FIX(1)); } if (d == 0.0) { - VALUE vsign = signbit(d) ? INT2FIX(-1) : INT2FIX(+1); - return rb_assoc_new(DBL2NUM(HUGE_VAL), vsign); + VALUE vsign = signbit(d) ? INT2FIX(-1) : INT2FIX(+1); + return rb_assoc_new(DBL2NUM(HUGE_VAL), vsign); } v = DBL2NUM(lgamma_r(d, &sign)); return rb_assoc_new(v, INT2FIX(sign)); |