Prefer integer as base of intermediate logarithms

As long as "floating point numbers" cannot accurately represent an
irrational number, the result of the natural logarithm cannot be
accurate.  Logarithms with an integer base may have the possibility to
represent more accurately.
This commit is contained in:
Nobuyoshi Nakada 2023-07-16 01:24:44 +09:00
parent be98bfc4ee
commit da39936ce1
No known key found for this signature in database
GPG Key ID: 3582D74E1FEE4465
Notes: Nobuyoshi Nakada 2023-07-16 13:41:21 +09:00

62
math.c
View File

@ -474,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));
/*
@ -509,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)
{
@ -541,16 +526,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;
}
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
#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) {
if (b > 0.0) return DBL2NUM(HUGE_VAL);
if (b < 0.0) return DBL2NUM(-HUGE_VAL);
return DBL2NUM(nan(""));
}
else if (b == 0.0) {
return DBL2NUM(-0.0);
}
d = log_intermediate(d) / log_intermediate(b);
numbits -= numbits_2;
}
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