mirror/guile
mirror/guile
1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-06-27 13:30:31 +02:00
Code Activity

Improvements to log' and log10'

Browse source 
* 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 2011-02-15 10:37:03 -05:00 • committed by Ludovic Courtès
parent 14a01ec1a8
commit 05f167beb2
2 changed files with 150 additions and 34 deletions
Download patch file Download diff file Expand all files Collapse all files

98
libguile/numbers.c
View file

@ -111,6 +111,7 @@ typedef scm_t_signed_bits scm_t_inum;
static SCM flo0; static SCM flo0;
static SCM exactly_one_half; static SCM exactly_one_half;
static SCM flo_log10e;
#define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0) #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() /* 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 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 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)); atan2 (im, re));
#endif #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 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
although the value itself overflows. */ if (scm_is_eq (z, SCM_INUM0))
double re = scm_to_double (z); scm_num_overflow (s_scm_log);
double l = log (fabs (re)); #endif
if (re >= 0.0) return log_of_shifted_double (SCM_I_INUM (z), 0);
return scm_from_double (l);
else
return scm_c_make_rectangular (l, M_PI);
} }
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 else
SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log); 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)); M_LOG10E * atan2 (im, re));
#endif #endif
} }
else if (SCM_NUMBERP (z)) else if (SCM_REALP (z) || SCM_I_INUMP (z))
{
#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
if (scm_is_eq (z, SCM_INUM0))
scm_num_overflow (s_scm_log10);
#endif
{ {
/* 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 re = scm_to_double (z);
double l = log10 (fabs (re)); double l = log10 (fabs (re));
if (re >= 0.0) if (re > 0.0 || double_is_non_negative_zero (re))
return scm_from_double (l); return scm_from_double (l);
else else
return scm_c_make_rectangular (l, M_LOG10E * M_PI); 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 else
SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10); 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 ("complex");
scm_add_feature ("inexact"); scm_add_feature ("inexact");
flo0 = scm_from_double (0.0); flo0 = scm_from_double (0.0);
flo_log10e = scm_from_double (M_LOG10E);
/* determine floating point precision */ /* determine floating point precision */
for (i=2; i <= SCM_MAX_DBL_RADIX; ++i) for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)

72
test-suite/tests/numbers.test
View file

@ -4323,14 +4323,36 @@
(log)) (log))
(pass-if-exception "two args" exception:wrong-num-args (pass-if-exception "two args" exception:wrong-num-args
(log 123 456)) (log 123 456))
(pass-if-exception "(log 0)" exception:numerical-overflow
(log 0))
(pass-if (negative-infinity? (log 0))) (pass-if (test-eqv? -inf.0 (log 0.0)))
(pass-if (negative-infinity? (log 0.0))) (pass-if (test-eqv? +inf.0 (log +inf.0)))
(pass-if (eqv? 0.0 (log 1))) (pass-if (test-eqv? -inf.0+3.14159265358979i (log -0.0)))
(pass-if (eqv? 0.0 (log 1.0))) (pass-if (test-eqv? +inf.0+3.14159265358979i (log -inf.0)))
(pass-if (eqv-loosely? 1.0 (log const-e))) (pass-if (test-eqv? 0.0 (log 1 )))
(pass-if (eqv-loosely? 2.0 (log const-e^2))) (pass-if (test-eqv? 0.0 (log 1.0)))
(pass-if (eqv-loosely? -1.0 (log const-1/e))) (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)))
(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,14 +4372,36 @@
(log10)) (log10))
(pass-if-exception "two args" exception:wrong-num-args (pass-if-exception "two args" exception:wrong-num-args
(log10 123 456)) (log10 123 456))
(pass-if-exception "(log10 0)" exception:numerical-overflow
(log10 0))
(pass-if (negative-infinity? (log10 0))) (pass-if (test-eqv? -inf.0 (log10 0.0)))
(pass-if (negative-infinity? (log10 0.0))) (pass-if (test-eqv? +inf.0 (log10 +inf.0)))
(pass-if (eqv? 0.0 (log10 1))) (pass-if (test-eqv? -inf.0+1.36437635384184i (log10 -0.0)))
(pass-if (eqv? 0.0 (log10 1.0))) (pass-if (test-eqv? +inf.0+1.36437635384184i (log10 -inf.0)))
(pass-if (eqv-loosely? 1.0 (log10 10.0))) (pass-if (test-eqv? 0.0 (log10 1 )))
(pass-if (eqv-loosely? 2.0 (log10 100.0))) (pass-if (test-eqv? 0.0 (log10 1.0)))
(pass-if (eqv-loosely? -1.0 (log10 0.1))) (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? 1.0-0.68218i (log10 0-10.0i))) (pass-if (eqv-loosely? 1.0-0.68218i (log10 0-10.0i)))