mirror of
https://git.savannah.gnu.org/git/guile.git
synced 2025-05-20 11:40:18 +02:00
Improvements to log' and log10'
* libguile/numbers.c (log_of_shifted_double, log_of_exact_integer, log_of_exact_integer_with_size, log_of_fraction): New internal static functions used by scm_log and scm_log10. (scm_log, scm_log10): Robustly handle large integers, large and small fractions, and fractions close to 1. Previously, computing logarithms of fractions close to 1 yielded grossly inaccurate results, and the other cases yielded infinities even though the answer could easily fit in a double. (log -0.0) now returns -inf.0+<PI>i, where previously it returned -inf.0. (log 0) now throws a numerical overflow exception, where previously it returned -inf.0. (log 0.0) still returns -inf.0. Analogous changes made to `log10'. * test-suite/tests/numbers.test (log, log10): Add tests. Signed-off-by: Ludovic Courtès <ludo@gnu.org>
This commit is contained in:
Mark H Weaver
Ludovic Courtès
parent
c05696aa94
commit
a5f6b751be
2 changed files with 150 additions and 34 deletions
|
|
@ -111,6 +111,7 @@ typedef scm_t_signed_bits scm_t_inum;
|
|||
|
||||
static SCM flo0;
|
||||
static SCM exactly_one_half;
|
||||
static SCM flo_log10e;
|
||||
|
||||
#define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
|
||||
|
||||
|
|
@ -9372,6 +9373,62 @@ scm_is_number (SCM z)
|
|||
}
|
||||
|
||||
|
||||
/* Returns log(x * 2^shift) */
|
||||
static SCM
|
||||
log_of_shifted_double (double x, long shift)
|
||||
{
|
||||
double ans = log (fabs (x)) + shift * M_LN2;
|
||||
|
||||
if (x > 0.0 || double_is_non_negative_zero (x))
|
||||
return scm_from_double (ans);
|
||||
else
|
||||
return scm_c_make_rectangular (ans, M_PI);
|
||||
}
|
||||
|
||||
/* Returns log(n), for exact integer n of integer-length size */
|
||||
static SCM
|
||||
log_of_exact_integer_with_size (SCM n, long size)
|
||||
{
|
||||
long shift = size - 2 * scm_dblprec[0];
|
||||
|
||||
if (shift > 0)
|
||||
return log_of_shifted_double
|
||||
(scm_to_double (scm_ash (n, scm_from_long(-shift))),
|
||||
shift);
|
||||
else
|
||||
return log_of_shifted_double (scm_to_double (n), 0);
|
||||
}
|
||||
|
||||
/* Returns log(n), for exact integer n of integer-length size */
|
||||
static SCM
|
||||
log_of_exact_integer (SCM n)
|
||||
{
|
||||
return log_of_exact_integer_with_size
|
||||
(n, scm_to_long (scm_integer_length (n)));
|
||||
}
|
||||
|
||||
/* Returns log(n/d), for exact non-zero integers n and d */
|
||||
static SCM
|
||||
log_of_fraction (SCM n, SCM d)
|
||||
{
|
||||
long n_size = scm_to_long (scm_integer_length (n));
|
||||
long d_size = scm_to_long (scm_integer_length (d));
|
||||
|
||||
if (abs (n_size - d_size) > 1)
|
||||
return (scm_difference (log_of_exact_integer_with_size (n, n_size),
|
||||
log_of_exact_integer_with_size (d, d_size)));
|
||||
else if (scm_is_false (scm_negative_p (n)))
|
||||
return scm_from_double
|
||||
(log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
|
||||
else
|
||||
return scm_c_make_rectangular
|
||||
(log1p (scm_to_double (scm_divide2real
|
||||
(scm_difference (scm_abs (n), d),
|
||||
d))),
|
||||
M_PI);
|
||||
}
|
||||
|
||||
|
||||
/* In the following functions we dispatch to the real-arg funcs like log()
|
||||
when we know the arg is real, instead of just handing everything to
|
||||
clog() for instance. This is in case clog() doesn't optimize for a
|
||||
|
|
@ -9394,17 +9451,21 @@ SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0,
|
|||
atan2 (im, re));
|
||||
#endif
|
||||
}
|
||||
else if (SCM_NUMBERP (z))
|
||||
else if (SCM_REALP (z))
|
||||
return log_of_shifted_double (SCM_REAL_VALUE (z), 0);
|
||||
else if (SCM_I_INUMP (z))
|
||||
{
|
||||
/* ENHANCE-ME: When z is a bignum the logarithm will fit a double
|
||||
although the value itself overflows. */
|
||||
double re = scm_to_double (z);
|
||||
double l = log (fabs (re));
|
||||
if (re >= 0.0)
|
||||
return scm_from_double (l);
|
||||
else
|
||||
return scm_c_make_rectangular (l, M_PI);
|
||||
#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
|
||||
if (scm_is_eq (z, SCM_INUM0))
|
||||
scm_num_overflow (s_scm_log);
|
||||
#endif
|
||||
return log_of_shifted_double (SCM_I_INUM (z), 0);
|
||||
}
|
||||
else if (SCM_BIGP (z))
|
||||
return log_of_exact_integer (z);
|
||||
else if (SCM_FRACTIONP (z))
|
||||
return log_of_fraction (SCM_FRACTION_NUMERATOR (z),
|
||||
SCM_FRACTION_DENOMINATOR (z));
|
||||
else
|
||||
SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log);
|
||||
}
|
||||
|
|
@ -9431,17 +9492,27 @@ SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0,
|
|||
M_LOG10E * atan2 (im, re));
|
||||
#endif
|
||||
}
|
||||
else if (SCM_NUMBERP (z))
|
||||
else if (SCM_REALP (z) || SCM_I_INUMP (z))
|
||||
{
|
||||
/* ENHANCE-ME: When z is a bignum the logarithm will fit a double
|
||||
although the value itself overflows. */
|
||||
double re = scm_to_double (z);
|
||||
double l = log10 (fabs (re));
|
||||
if (re >= 0.0)
|
||||
return scm_from_double (l);
|
||||
else
|
||||
return scm_c_make_rectangular (l, M_LOG10E * M_PI);
|
||||
#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
|
||||
if (scm_is_eq (z, SCM_INUM0))
|
||||
scm_num_overflow (s_scm_log10);
|
||||
#endif
|
||||
{
|
||||
double re = scm_to_double (z);
|
||||
double l = log10 (fabs (re));
|
||||
if (re > 0.0 || double_is_non_negative_zero (re))
|
||||
return scm_from_double (l);
|
||||
else
|
||||
return scm_c_make_rectangular (l, M_LOG10E * M_PI);
|
||||
}
|
||||
}
|
||||
else if (SCM_BIGP (z))
|
||||
return scm_product (flo_log10e, log_of_exact_integer (z));
|
||||
else if (SCM_FRACTIONP (z))
|
||||
return scm_product (flo_log10e,
|
||||
log_of_fraction (SCM_FRACTION_NUMERATOR (z),
|
||||
SCM_FRACTION_DENOMINATOR (z)));
|
||||
else
|
||||
SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10);
|
||||
}
|
||||
|
|
@ -9536,6 +9607,7 @@ scm_init_numbers ()
|
|||
scm_add_feature ("complex");
|
||||
scm_add_feature ("inexact");
|
||||
flo0 = scm_from_double (0.0);
|
||||
flo_log10e = scm_from_double (M_LOG10E);
|
||||
|
||||
/* determine floating point precision */
|
||||
for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
|
||||
|
|
|
|||
|
|
@ -4323,14 +4323,36 @@
|
|||
(log))
|
||||
(pass-if-exception "two args" exception:wrong-num-args
|
||||
(log 123 456))
|
||||
(pass-if-exception "(log 0)" exception:numerical-overflow
|
||||
(log 0))
|
||||
|
||||
(pass-if (negative-infinity? (log 0)))
|
||||
(pass-if (negative-infinity? (log 0.0)))
|
||||
(pass-if (eqv? 0.0 (log 1)))
|
||||
(pass-if (eqv? 0.0 (log 1.0)))
|
||||
(pass-if (eqv-loosely? 1.0 (log const-e)))
|
||||
(pass-if (eqv-loosely? 2.0 (log const-e^2)))
|
||||
(pass-if (eqv-loosely? -1.0 (log const-1/e)))
|
||||
(pass-if (test-eqv? -inf.0 (log 0.0)))
|
||||
(pass-if (test-eqv? +inf.0 (log +inf.0)))
|
||||
(pass-if (test-eqv? -inf.0+3.14159265358979i (log -0.0)))
|
||||
(pass-if (test-eqv? +inf.0+3.14159265358979i (log -inf.0)))
|
||||
(pass-if (test-eqv? 0.0 (log 1 )))
|
||||
(pass-if (test-eqv? 0.0 (log 1.0)))
|
||||
(pass-if (test-eqv? 1.0 (log const-e)))
|
||||
(pass-if (test-eqv? 2.0 (log const-e^2)))
|
||||
(pass-if (test-eqv? -1.0 (log const-1/e)))
|
||||
(pass-if (test-eqv? -1.0+3.14159265358979i (log (- const-1/e))))
|
||||
(pass-if (test-eqv? 2.30258509299405 (log 10)))
|
||||
(pass-if (test-eqv? 2.30258509299405+3.14159265358979i (log -10)))
|
||||
|
||||
(pass-if (test-eqv? 1.0+0.0i (log (+ const-e +0.0i))))
|
||||
(pass-if (test-eqv? 1.0-0.0i (log (+ const-e -0.0i))))
|
||||
|
||||
(pass-if (eqv-loosely? 230258.509299405 (log (expt 10 100000))))
|
||||
(pass-if (eqv-loosely? -230258.509299405 (log (expt 10 -100000))))
|
||||
(pass-if (eqv-loosely? 230257.410687116 (log (/ (expt 10 100000) 3))))
|
||||
(pass-if (eqv-loosely? 230258.509299405+3.14159265358979i
|
||||
(log (- (expt 10 100000)))))
|
||||
(pass-if (eqv-loosely? -230258.509299405+3.14159265358979i
|
||||
(log (- (expt 10 -100000)))))
|
||||
(pass-if (eqv-loosely? 230257.410687116+3.14159265358979i
|
||||
(log (- (/ (expt 10 100000) 3)))))
|
||||
(pass-if (test-eqv? 3.05493636349961e-151
|
||||
(log (/ (1+ (expt 2 500)) (expt 2 500)))))
|
||||
|
||||
(pass-if (eqv-loosely? 1.0+1.57079i (log 0+2.71828i)))
|
||||
(pass-if (eqv-loosely? 1.0-1.57079i (log 0-2.71828i)))
|
||||
|
|
@ -4350,20 +4372,42 @@
|
|||
(log10))
|
||||
(pass-if-exception "two args" exception:wrong-num-args
|
||||
(log10 123 456))
|
||||
(pass-if-exception "(log10 0)" exception:numerical-overflow
|
||||
(log10 0))
|
||||
|
||||
(pass-if (negative-infinity? (log10 0)))
|
||||
(pass-if (negative-infinity? (log10 0.0)))
|
||||
(pass-if (eqv? 0.0 (log10 1)))
|
||||
(pass-if (eqv? 0.0 (log10 1.0)))
|
||||
(pass-if (eqv-loosely? 1.0 (log10 10.0)))
|
||||
(pass-if (eqv-loosely? 2.0 (log10 100.0)))
|
||||
(pass-if (eqv-loosely? -1.0 (log10 0.1)))
|
||||
(pass-if (test-eqv? -inf.0 (log10 0.0)))
|
||||
(pass-if (test-eqv? +inf.0 (log10 +inf.0)))
|
||||
(pass-if (test-eqv? -inf.0+1.36437635384184i (log10 -0.0)))
|
||||
(pass-if (test-eqv? +inf.0+1.36437635384184i (log10 -inf.0)))
|
||||
(pass-if (test-eqv? 0.0 (log10 1 )))
|
||||
(pass-if (test-eqv? 0.0 (log10 1.0)))
|
||||
(pass-if (test-eqv? 1.0 (log10 10 )))
|
||||
(pass-if (test-eqv? 1.0 (log10 10.0)))
|
||||
(pass-if (test-eqv? 2.0 (log10 100.0)))
|
||||
(pass-if (test-eqv? -1.0 (log10 0.1)))
|
||||
(pass-if (test-eqv? -1.0+1.36437635384184i (log10 -0.1)))
|
||||
(pass-if (test-eqv? 1.0+1.36437635384184i (log10 -10 )))
|
||||
|
||||
(pass-if (test-eqv? 1.0+0.0i (log10 10.0+0.0i)))
|
||||
(pass-if (test-eqv? 1.0-0.0i (log10 10.0-0.0i)))
|
||||
|
||||
(pass-if (eqv-loosely? 100000.0 (log10 (expt 10 100000))))
|
||||
(pass-if (eqv-loosely? -100000.0 (log10 (expt 10 -100000))))
|
||||
(pass-if (eqv-loosely? 99999.5228787453 (log10 (/ (expt 10 100000) 3))))
|
||||
(pass-if (eqv-loosely? 100000.0+1.36437635384184i
|
||||
(log10 (- (expt 10 100000)))))
|
||||
(pass-if (eqv-loosely? -100000.0+1.36437635384184i
|
||||
(log10 (- (expt 10 -100000)))))
|
||||
(pass-if (eqv-loosely? 99999.5228787453+1.36437635384184i
|
||||
(log10 (- (/ (expt 10 100000) 3)))))
|
||||
(pass-if (test-eqv? 1.32674200523347e-151
|
||||
(log10 (/ (1+ (expt 2 500)) (expt 2 500)))))
|
||||
|
||||
(pass-if (eqv-loosely? 1.0+0.68218i (log10 0+10.0i)))
|
||||
(pass-if (eqv-loosely? 1.0-0.68218i (log10 0-10.0i)))
|
||||
|
||||
(pass-if (eqv-loosely? 0.0+1.36437i (log10 -1)))
|
||||
(pass-if (eqv-loosely? 1.0+1.36437i (log10 -10)))
|
||||
(pass-if (eqv-loosely? 0.0+1.36437i (log10 -1)))
|
||||
(pass-if (eqv-loosely? 1.0+1.36437i (log10 -10)))
|
||||
(pass-if (eqv-loosely? 2.0+1.36437i (log10 -100))))
|
||||
|
||||
;;;
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue