summaryrefslogtreecommitdiff
path: root/math.c
diff options
context:
space:
mode:
Diffstat (limited to 'math.c')
-rw-r--r--math.c222
1 files changed, 155 insertions, 67 deletions
diff --git a/math.c b/math.c
index 4996718b10..852620da20 100644
--- a/math.c
+++ b/math.c
@@ -15,7 +15,6 @@
# define _USE_MATH_DEFINES 1
#endif
-#include <errno.h>
#include <float.h>
#include <math.h>
@@ -65,23 +64,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 +169,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 +195,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 +218,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)
}
/*
@@ -460,6 +457,34 @@ math_exp(VALUE unused_obj, VALUE x)
return DBL2NUM(exp(Get_Double(x)));
}
+/*
+ * call-seq:
+ * Math.expm1(x) -> float
+ *
+ * Returns "exp(x) - 1", +e+ raised to the +x+ power, minus 1.
+ *
+ *
+ * - Domain: <tt>[-INFINITY, INFINITY]</tt>.
+ * - Range: <tt>[-1.0, INFINITY]</tt>.
+ *
+ * Examples:
+ *
+ * expm1(-INFINITY) # => 0.0
+ * expm1(-1.0) # => -0.6321205588285577 # 1.0/E - 1
+ * expm1(0.0) # => 0.0
+ * expm1(0.5) # => 0.6487212707001282 # sqrt(E) - 1
+ * expm1(1.0) # => 1.718281828459045 # E - 1
+ * expm1(2.0) # => 6.38905609893065 # E**2 - 1
+ * expm1(INFINITY) # => Infinity
+ *
+ */
+
+static VALUE
+math_expm1(VALUE unused_obj, VALUE x)
+{
+ return DBL2NUM(expm1(Get_Double(x)));
+}
+
#if defined __CYGWIN__
# include <cygwin/version.h>
# if CYGWIN_VERSION_DLL_MAJOR < 1005
@@ -476,7 +501,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 +535,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 +546,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;
+}
- return log(d) + numbits * M_LN2; /* log(d * 2 ** numbits) */
+#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;
+
+ 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
@@ -629,6 +674,47 @@ math_log10(VALUE unused_obj, VALUE x)
return DBL2NUM(log10(d) + numbits * log10(2)); /* log10(d * 2 ** numbits) */
}
+/*
+ * call-seq:
+ * Math.log1p(x) -> float
+ *
+ * Returns "log(x + 1)", the base E {logarithm}[https://en.wikipedia.org/wiki/Logarithm] of (+x+ + 1).
+ *
+ * - Domain: <tt>[-1, INFINITY]</tt>.
+ * - Range: <tt>[-INFINITY, INFINITY]</tt>.
+ *
+ * Examples:
+ *
+ * log1p(-1.0) # => -Infinity
+ * log1p(0.0) # => 0.0
+ * log1p(E - 1) # => 1.0
+ * log1p(INFINITY) # => Infinity
+ *
+ */
+
+static VALUE
+math_log1p(VALUE unused_obj, VALUE x)
+{
+ size_t numbits;
+ double d = get_double_rshift(x, &numbits);
+
+ if (numbits != 0) {
+ x = rb_big_plus(x, INT2FIX(1));
+ d = math_log_split(x, &numbits);
+ /* check for pole error */
+ if (d == 0.0) return DBL2NUM(-HUGE_VAL);
+ d = log(d);
+ d += numbits * M_LN2;
+ return DBL2NUM(d);
+ }
+
+ domain_check_min(d, -1.0, "log1p");
+ /* check for pole error */
+ if (d == -1.0) return DBL2NUM(-HUGE_VAL);
+
+ return DBL2NUM(log1p(d)); /* log10(d * 2 ** numbits) */
+}
+
static VALUE rb_math_sqrt(VALUE x);
/*
@@ -681,13 +767,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 +799,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 +813,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 +1031,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));
}
@@ -968,7 +1054,7 @@ math_gamma(VALUE unused_obj, VALUE x)
*
* [Math.log(Math.gamma(x).abs), Math.gamma(x) < 0 ? -1 : 1]
*
- * See {logarithmic gamma function}[https://en.wikipedia.org/wiki/Gamma_function#The_log-gamma_function].
+ * See {log gamma function}[https://en.wikipedia.org/wiki/Gamma_function#Log-gamma_function].
*
* - Domain: <tt>(-INFINITY, INFINITY]</tt>.
* - Range of first element: <tt>(-INFINITY, INFINITY]</tt>.
@@ -1007,12 +1093,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));
@@ -1103,9 +1189,11 @@ InitVM_Math(void)
rb_define_module_function(rb_mMath, "atanh", math_atanh, 1);
rb_define_module_function(rb_mMath, "exp", math_exp, 1);
+ rb_define_module_function(rb_mMath, "expm1", math_expm1, 1);
rb_define_module_function(rb_mMath, "log", math_log, -1);
rb_define_module_function(rb_mMath, "log2", math_log2, 1);
rb_define_module_function(rb_mMath, "log10", math_log10, 1);
+ rb_define_module_function(rb_mMath, "log1p", math_log1p, 1);
rb_define_module_function(rb_mMath, "sqrt", math_sqrt, 1);
rb_define_module_function(rb_mMath, "cbrt", math_cbrt, 1);
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-21 19:33:46 +0000