From 3f64e2e88b7879e079a2dc65b9856158c11a2a63 Mon Sep 17 00:00:00 2001 From: Marius Gerbershagen Date: Mon, 2 Jan 2023 16:16:32 +0100 Subject: [PATCH] numbers/log.d: prevent unnecessary loss of precision For (log x y) where one of the two arguments is a double or long float and the other a rational number, defining (log x y) as (/ (log y) (log x)) is imprecise since intermediate single float results will be used for the rational argument. Prevent this by computing the logarithm of the rational to the same precision as that of the other argument. Fixes #688. --- src/c/numbers/log.d | 127 +++++++++++++++++++++++++++++++ src/tests/normal-tests/mixed.lsp | 19 +++++ 2 files changed, 146 insertions(+) diff --git a/src/c/numbers/log.d b/src/c/numbers/log.d index d2a87fb42..811b4f4a0 100644 --- a/src/c/numbers/log.d +++ b/src/c/numbers/log.d @@ -189,9 +189,136 @@ MATH_DEF_DISPATCH1(log1, @[log], @[number], ecl_log1_complex, ecl_log1_csfloat, ecl_log1_cdfloat, ecl_log1_clfloat); +static cl_object ecl_log1_double_precision(cl_object x); + +static cl_object +ecl_log1_bignum_double_precision(cl_object x) +{ + cl_fixnum l = ecl_integer_length(x) - 1; + cl_object r = ecl_make_ratio(x, ecl_ash(ecl_make_fixnum(1), l)); + double f = ecl_to_double(r); + if (f < 0) { +#ifdef ECL_COMPLEX_FLOAT + return ecl_make_cdfloat(clog(f) + l * log(2.0)); +#else + return ecl_make_complex(ecl_make_double_float(log(-f) + l * log(2.0)), ecl_make_double_float(ECL_PI_D)); +#endif + } else { + return ecl_make_double_float(log(f) + l * log(2.0)); + } +} + +static cl_object +ecl_log1_simple_double_precision(cl_object x) +{ + double f = ecl_to_double(x); + if (f < 0) { +#ifdef ECL_COMPLEX_FLOAT + return ecl_make_cdfloat(clog(f)); +#else + return ecl_make_complex(ecl_make_double_float(log(-f)), ecl_make_double_float(ECL_PI_D)); +#endif + } + return ecl_make_double_float(log(f)); +} + +static cl_object +ecl_log1_ratio_double_precision(cl_object x) +{ + cl_object num = x->ratio.num; + cl_object den = x->ratio.den; + cl_index lnum = ecl_integer_length(num); + cl_index lden = ecl_integer_length(den); + if ((lnum > lden) ? (lnum - lden >= DBL_MAX_EXP) : (lden - lnum >= -DBL_MIN_EXP)) { + cl_object numlog = ecl_log1_double_precision(num); + cl_object denlog = ecl_log1_double_precision(den); + return ecl_minus(numlog, denlog); + } else { + return ecl_log1_simple_double_precision(x); + } +} + +MATH_DEF_DISPATCH1(log1_double_precision, @[log], @[number], + ecl_log1_simple_double_precision, ecl_log1_bignum_double_precision, ecl_log1_ratio_double_precision, + ecl_log1_single_float, ecl_log1_double_float, ecl_log1_long_float, + ecl_log1_complex, + ecl_log1_csfloat, ecl_log1_cdfloat, ecl_log1_clfloat); + +static cl_object ecl_log1_long_precision(cl_object x); + +static cl_object +ecl_log1_bignum_long_precision(cl_object x) +{ + cl_fixnum l = ecl_integer_length(x) - 1; + cl_object r = ecl_make_ratio(x, ecl_ash(ecl_make_fixnum(1), l)); + long double f = ecl_to_long_double(r); + if (f < 0) { +#ifdef ECL_COMPLEX_FLOAT + return ecl_make_clfloat(clogl(f) + l * logl(2.0)); +#else + return ecl_make_complex(ecl_make_long_float(logl(-f) + l * logl(2.0)), ecl_make_long_float(ECL_PI_L)); +#endif + } else { + return ecl_make_long_float(logl(f) + l * logl(2.0)); + } +} + +static cl_object +ecl_log1_simple_long_precision(cl_object x) +{ + long double f = ecl_to_long_double(x); + if (f < 0) { +#ifdef ECL_COMPLEX_FLOAT + return ecl_make_clfloat(clogl(f)); +#else + return ecl_make_complex(ecl_make_long_float(logl(-f)), ecl_make_long_float(ECL_PI_L)); +#endif + } + return ecl_make_long_float(logl(f)); +} + +static cl_object +ecl_log1_ratio_long_precision(cl_object x) +{ + cl_object num = x->ratio.num; + cl_object den = x->ratio.den; + cl_index lnum = ecl_integer_length(num); + cl_index lden = ecl_integer_length(den); + if ((lnum > lden) ? (lnum - lden >= LDBL_MAX_EXP) : (lden - lnum >= -LDBL_MIN_EXP)) { + cl_object numlog = ecl_log1_long_precision(num); + cl_object denlog = ecl_log1_long_precision(den); + return ecl_minus(numlog, denlog); + } else { + return ecl_log1_simple_long_precision(x); + } +} + +MATH_DEF_DISPATCH1(log1_long_precision, @[log], @[number], + ecl_log1_simple_long_precision, ecl_log1_bignum_long_precision, ecl_log1_ratio_long_precision, + ecl_log1_single_float, ecl_log1_double_float, ecl_log1_long_float, + ecl_log1_complex, + ecl_log1_csfloat, ecl_log1_cdfloat, ecl_log1_clfloat); + cl_object ecl_log2(cl_object x, cl_object y) { + cl_type tx = ecl_t_of(x); + cl_type ty = ecl_t_of(y); + /* Prevent loss of precision from intermediate single float results */ + if (tx == t_doublefloat || ty == t_doublefloat +#ifdef ECL_COMPLEX_FLOAT + || tx == t_cdfloat || ty == t_cdfloat +#endif + ) { + return ecl_divide(ecl_log1_double_precision(y), ecl_log1_double_precision(x)); + } + if (tx == t_longfloat || ty == t_longfloat +#ifdef ECL_COMPLEX_FLOAT + || tx == t_clfloat || ty == t_clfloat +#endif + ) { + return ecl_divide(ecl_log1_long_precision(y), ecl_log1_long_precision(x)); + } return ecl_divide(ecl_log1(y), ecl_log1(x)); } diff --git a/src/tests/normal-tests/mixed.lsp b/src/tests/normal-tests/mixed.lsp index 037e9892a..d8549b716 100644 --- a/src/tests/normal-tests/mixed.lsp +++ b/src/tests/normal-tests/mixed.lsp @@ -441,3 +441,22 @@ (finishes (log (- x))) (finishes (log (/ 1 x))) (finishes (log (- (/ 1 x)))))) + +;;; Created: 2023-01-02 +;;; Contains: tests checking that (log x y) does not unnecessarily +;;; lose precision through intermediate single-float calculations when +;;; the final result is a double or long float. +(test mix.0024.log-loss-of-precision + (is (eql (log 2 2d0) 1d0)) + (is (eql (realpart (log -2 2d0)) 1d0)) + (is (eql (log (ash 1 1024) 2d0) 1024d0)) + (is (eql (realpart (log (- (ash 1 1024)) 2d0)) 1024d0)) + (is (eql (log 1/2 2d0) -1d0)) + (is (eql (realpart (log -1/2 2d0)) -1d0)) + + (is (eql (log 2 2l0) 1l0)) + (is (eql (realpart (log -2 2l0)) 1l0)) + (is (eql (log (ash 1 1024) 2l0) 1024l0)) + (is (eql (realpart (log (- (ash 1 1024)) 2l0)) 1024l0)) + (is (eql (log 1/2 2l0) -1l0)) + (is (eql (realpart (log -1/2 2l0)) -1l0)))