summaryrefslogtreecommitdiff
path: root/math.c
diff options
context:
space:
mode:
Diffstat (limited to 'math.c')
-rw-r--r--math.c148
1 files changed, 83 insertions, 65 deletions
diff --git a/math.c b/math.c
index 4996718b10..2394fe9f58 100644
--- a/math.c
+++ b/math.c
@@ -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));
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-19 21:30:31 +0000