diff options
Diffstat (limited to 'math.c')
| -rw-r--r-- | math.c | 156 |
1 files changed, 122 insertions, 34 deletions
@@ -15,7 +15,6 @@ # define _USE_MATH_DEFINES 1 #endif -#include <errno.h> #include <float.h> #include <math.h> @@ -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) { @@ -543,16 +553,51 @@ get_double_rshift(VALUE x, size_t *pnumbits) } 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 @@ -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); /* @@ -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 @@ -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>. @@ -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); |