[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[PATCH] Improvements to exact rationals et al
|
From: |
Mark H Weaver |
|
Subject: |
[PATCH] Improvements to exact rationals et al |
|
Date: |
Fri, 07 Oct 2011 01:25:02 -0400 |
Hello all,
After a long hiatus, and newly inspired by both peval and the videos
from GHM 2011 (especially Andy's excellent talk on the importance of
extensibility) I'd like to try to find more time to hack on Guile.
However, before embarking on more interesting work, I'd like to submit
some numerics patches that have been lingering in my personal repository
for far too long.
This patch set makes some cleanups and improvements to the numerics
code, especially to the handling of exact rationals. Here are some
examples of the visible improvements:
Before:
(exact->inexact (/ 1000001 (expt 10 310))) ==> 0.0
(exact->inexact (/ (expt 10 310) 1000001)) ==> +inf.0
After:
(exact->inexact (/ 1000001 (expt 10 310))) ==> 1.000001e-304
(exact->inexact (/ (expt 10 310) 1000001)) ==> 9.99999000001e303
Before: (on the Yeeloong)
scheme@(guile-user)> ,time (define r (expt 11/13 100000))
;; 3.506000s real time, 3.145000s run time. 0.000000s spent in GC.
scheme@(guile-user)> ,time (define r (expt 11/13 100000))
;; 3.250000s real time, 3.140000s run time. 0.000000s spent in GC.
scheme@(guile-user)> ,time (define r (expt 11/13 100000))
;; 3.242000s real time, 3.141000s run time. 0.000000s spent in GC.
After:
scheme@(guile-user)> ,time (define r (expt 11/13 100000))
;; 0.305000s real time, 0.273000s run time. 0.000000s spent in GC.
scheme@(guile-user)> ,time (define r (expt 11/13 100000))
;; 0.202000s real time, 0.187000s run time. 0.000000s spent in GC.
scheme@(guile-user)> ,time (define r (expt 11/13 100000))
;; 0.208000s real time, 0.190000s run time. 0.000000s spent in GC.
Also, this patch set adds the user-visible procedure `round-ash' which I
proposed several months ago, because its functionality was needed
internally to implement proper rounding in exact->inexact for bignums.
While it doesn't strictly need to be made public, it will certainly be
needed to implement many kinds of inexact GOOPS numbers efficiently such
as fixed-point and arbitrary-precision floats, and I suspect it has
wider utility as well.
I'd like to apply these patches to stable-2.0. Comments and suggestions
are welcome.
Mark
>From 80f12d13984dc01ff74620e46c46dda9e3dea226 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Tue, 1 Mar 2011 20:12:51 -0500
Subject: [PATCH 1/7] Improve code in scm_gcd for inum/inum case
* libguile/numbers.c (scm_gcd): Improve implementation of inum/inum case
in scm_gcd to be more clear and efficient.
---
libguile/numbers.c | 54 ++++++++++++++++++++++++++++-----------------------
1 files changed, 30 insertions(+), 24 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index b01af9f..aa768b4 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -3853,52 +3853,58 @@ SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1,
SCM
scm_gcd (SCM x, SCM y)
{
- if (SCM_UNBNDP (y))
+ if (SCM_UNLIKELY (SCM_UNBNDP (y)))
return SCM_UNBNDP (x) ? SCM_INUM0 : scm_abs (x);
- if (SCM_I_INUMP (x))
+ if (SCM_LIKELY (SCM_I_INUMP (x)))
{
- if (SCM_I_INUMP (y))
+ if (SCM_LIKELY (SCM_I_INUMP (y)))
{
scm_t_inum xx = SCM_I_INUM (x);
scm_t_inum yy = SCM_I_INUM (y);
scm_t_inum u = xx < 0 ? -xx : xx;
scm_t_inum v = yy < 0 ? -yy : yy;
scm_t_inum result;
- if (xx == 0)
+ if (SCM_UNLIKELY (xx == 0))
result = v;
- else if (yy == 0)
+ else if (SCM_UNLIKELY (yy == 0))
result = u;
else
{
- scm_t_inum k = 1;
- scm_t_inum t;
+ int k = 0;
/* Determine a common factor 2^k */
- while (!(1 & (u | v)))
+ while (((u | v) & 1) == 0)
{
- k <<= 1;
+ k++;
u >>= 1;
v >>= 1;
}
/* Now, any factor 2^n can be eliminated */
- if (u & 1)
- t = -v;
+ if ((u & 1) == 0)
+ while ((u & 1) == 0)
+ u >>= 1;
else
+ while ((v & 1) == 0)
+ v >>= 1;
+ /* Both u and v are now odd. Subtract the smaller one
+ from the larger one to produce an even number, remove
+ more factors of two, and repeat. */
+ while (u != v)
{
- t = u;
- b3:
- t = SCM_SRS (t, 1);
+ if (u > v)
+ {
+ u -= v;
+ while ((u & 1) == 0)
+ u >>= 1;
+ }
+ else
+ {
+ v -= u;
+ while ((v & 1) == 0)
+ v >>= 1;
+ }
}
- if (!(1 & t))
- goto b3;
- if (t > 0)
- u = t;
- else
- v = -t;
- t = u - v;
- if (t != 0)
- goto b3;
- result = u * k;
+ result = u << k;
}
return (SCM_POSFIXABLE (result)
? SCM_I_MAKINUM (result)
--
1.7.2.5
>From 807de5ec6dab40fc699f7e3353de36ae1ee36604 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Tue, 1 Mar 2011 16:03:59 -0500
Subject: [PATCH 2/7] Optimize and simplify fractions code
* libguile/numbers.c (scm_exact_integer_quotient): New internal static
function that computes the quotient of two exact integers when the
remainder is known in advance to be zero. For large integers this can
be implemented more efficiently than when the remainder is unknown.
(scm_i_make_ratio_already_reduced): New internal static function that
creates a ratio when the numerator and denominator are already known
to be reduced to lowest terms (i.e. when their gcd is 1). This can be
used in several places to avoid unnecessary gcd computations.
(scm_i_make_ratio): Rewrite in terms of
scm_i_make_ratio_already_reduced. Don't bother checking to see if the
denominator divides the numerator evenly. This is wasted effort in
the common case. Instead, compute the gcd, reduce to lowest terms
(using scm_exact_integer_quotient), and let
scm_i_make_ratio_already_reduced do the integer check (by checking for
a unit denominator).
(scm_integer_expt): Optimize fraction case by (a/b)^n ==> (a^n)/(b^n),
to avoid needless effort to reduce intermediate products to lowest
terms. If a and b have no common factors, it follows that a^n and b^n
don't have common factors. Use scm_i_make_ratio_already_reduced to
construct the final result, so that no gcd computations are needed to
exponentiate a fraction.
(scm_abs, scm_magnitude, scm_difference): Use
scm_i_make_ratio_already_reduced to avoid wasteful gcd computations
when negating fractions.
(do_divide): Use scm_i_make_ratio_already_reduced to avoid wasteful
gcd computations when taking the reciprocal of a fraction or integer.
When dividing exact integers, don't bother checking to see if the
denominator divides the numerator evenly. Instead, let
scm_i_make_ratio do the job.
* test-suite/tests/numbers.test (expt, integer-expt): Add tests.
---
libguile/numbers.c | 287 ++++++++++++++++++++++-------------------
test-suite/tests/numbers.test | 6 +
2 files changed, 161 insertions(+), 132 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index aa768b4..80b46df 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -413,96 +413,58 @@ scm_i_mpz2num (mpz_t b)
/* this is needed when we want scm_divide to make a float, not a ratio, even
if passed two ints */
static SCM scm_divide2real (SCM x, SCM y);
+/* scm_i_make_ratio_already_reduced makes a ratio more efficiently,
+ when the following conditions are known in advance:
+
+ 1. numerator and denominator are exact integers
+ 2. numerator and denominator are reduced to lowest terms: gcd(n,d) == 1
+*/
static SCM
-scm_i_make_ratio (SCM numerator, SCM denominator)
-#define FUNC_NAME "make-ratio"
+scm_i_make_ratio_already_reduced (SCM numerator, SCM denominator)
{
- /* First make sure the arguments are proper.
- */
- if (SCM_I_INUMP (denominator))
+ /* Flip signs so that the denominator is positive. */
+ if (scm_is_false (scm_positive_p (denominator)))
{
- if (scm_is_eq (denominator, SCM_INUM0))
+ if (SCM_UNLIKELY (scm_is_eq (denominator, SCM_INUM0)))
scm_num_overflow ("make-ratio");
- if (scm_is_eq (denominator, SCM_INUM1))
- return numerator;
- }
- else
- {
- if (!(SCM_BIGP(denominator)))
- SCM_WRONG_TYPE_ARG (2, denominator);
- }
- if (!SCM_I_INUMP (numerator) && !SCM_BIGP (numerator))
- SCM_WRONG_TYPE_ARG (1, numerator);
-
- /* Then flip signs so that the denominator is positive.
- */
- if (scm_is_true (scm_negative_p (denominator)))
- {
- numerator = scm_difference (numerator, SCM_UNDEFINED);
- denominator = scm_difference (denominator, SCM_UNDEFINED);
- }
-
- /* Now consider for each of the four fixnum/bignum combinations
- whether the rational number is really an integer.
- */
- if (SCM_I_INUMP (numerator))
- {
- scm_t_inum x = SCM_I_INUM (numerator);
- if (scm_is_eq (numerator, SCM_INUM0))
- return SCM_INUM0;
- if (SCM_I_INUMP (denominator))
+ else
{
- scm_t_inum y;
- y = SCM_I_INUM (denominator);
- if (x == y)
- return SCM_INUM1;
- if ((x % y) == 0)
- return SCM_I_MAKINUM (x / y);
+ numerator = scm_difference (numerator, SCM_UNDEFINED);
+ denominator = scm_difference (denominator, SCM_UNDEFINED);
}
- else
- {
- /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
- of that value for the denominator, as a bignum. Apart from
- that case, abs(bignum) > abs(inum) so inum/bignum is not an
- integer. */
- if (x == SCM_MOST_NEGATIVE_FIXNUM
- && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator),
- - SCM_MOST_NEGATIVE_FIXNUM) == 0)
- return SCM_I_MAKINUM(-1);
- }
}
- else if (SCM_BIGP (numerator))
+
+ /* Check for the integer case */
+ if (scm_is_eq (denominator, SCM_INUM1))
+ return numerator;
+
+ return scm_double_cell (scm_tc16_fraction,
+ SCM_UNPACK (numerator),
+ SCM_UNPACK (denominator), 0);
+}
+
+static SCM scm_exact_integer_quotient (SCM x, SCM y);
+
+static SCM
+scm_i_make_ratio (SCM numerator, SCM denominator)
+#define FUNC_NAME "make-ratio"
+{
+ /* Make sure the arguments are proper */
+ if (!SCM_LIKELY (SCM_I_INUMP (numerator) || SCM_BIGP (numerator)))
+ SCM_WRONG_TYPE_ARG (1, numerator);
+ else if (!SCM_LIKELY (SCM_I_INUMP (denominator) || SCM_BIGP (denominator)))
+ SCM_WRONG_TYPE_ARG (2, denominator);
+ else
{
- if (SCM_I_INUMP (denominator))
+ SCM the_gcd = scm_gcd (numerator, denominator);
+ if (!(scm_is_eq (the_gcd, SCM_INUM1)))
{
- scm_t_inum yy = SCM_I_INUM (denominator);
- if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy))
- return scm_divide (numerator, denominator);
- }
- else
- {
- if (scm_is_eq (numerator, denominator))
- return SCM_INUM1;
- if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator),
- SCM_I_BIG_MPZ (denominator)))
- return scm_divide(numerator, denominator);
+ /* Reduce to lowest terms */
+ numerator = scm_exact_integer_quotient (numerator, the_gcd);
+ denominator = scm_exact_integer_quotient (denominator, the_gcd);
}
+ return scm_i_make_ratio_already_reduced (numerator, denominator);
}
-
- /* No, it's a proper fraction.
- */
- {
- SCM divisor = scm_gcd (numerator, denominator);
- if (!(scm_is_eq (divisor, SCM_INUM1)))
- {
- numerator = scm_divide (numerator, divisor);
- denominator = scm_divide (denominator, divisor);
- }
-
- return scm_double_cell (scm_tc16_fraction,
- SCM_UNPACK (numerator),
- SCM_UNPACK (denominator), 0);
- }
}
#undef FUNC_NAME
@@ -784,8 +746,9 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
{
if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x))))
return x;
- return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x),
SCM_UNDEFINED),
- SCM_FRACTION_DENOMINATOR (x));
+ return scm_i_make_ratio_already_reduced
+ (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
+ SCM_FRACTION_DENOMINATOR (x));
}
else
SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs);
@@ -853,6 +816,83 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0,
}
#undef FUNC_NAME
+/* scm_exact_integer_quotient returns the exact integer q such that
+ n = q*d, for exact integers n and d, where d is known in advance to
+ divide n evenly (with zero remainder). */
+static SCM
+scm_exact_integer_quotient (SCM n, SCM d)
+#define FUNC_NAME "exact-integer-quotient"
+{
+ if (SCM_LIKELY (SCM_I_INUMP (n)))
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+ if (SCM_LIKELY (SCM_I_INUMP (d)))
+ {
+ scm_t_inum dd = SCM_I_INUM (d);
+ if (SCM_UNLIKELY (dd == 0))
+ scm_num_overflow ("exact-integer-quotient");
+ else
+ {
+ scm_t_inum qq = nn / dd;
+ if (SCM_LIKELY (SCM_FIXABLE (qq)))
+ return SCM_I_MAKINUM (qq);
+ else
+ return scm_i_inum2big (qq);
+ }
+ }
+ else if (SCM_LIKELY (SCM_BIGP (d)))
+ {
+ /* n is an inum and d is a bignum. Given that d is known to
+ divide n evenly, there are only two possibilities: n is 0,
+ or else n is fixnum-min and d is abs(fixnum-min). */
+ if (nn == 0)
+ return SCM_INUM0;
+ else
+ return SCM_I_MAKINUM (-1);
+ }
+ else
+ SCM_WRONG_TYPE_ARG (2, d);
+ }
+ else if (SCM_LIKELY (SCM_BIGP (n)))
+ {
+ if (SCM_LIKELY (SCM_I_INUMP (d)))
+ {
+ scm_t_inum dd = SCM_I_INUM (d);
+ if (SCM_UNLIKELY (dd == 0))
+ scm_num_overflow ("exact-integer-quotient");
+ else if (SCM_UNLIKELY (dd == 1))
+ return n;
+ else
+ {
+ SCM q = scm_i_mkbig ();
+ if (dd > 0)
+ mpz_divexact_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), dd);
+ else
+ {
+ mpz_divexact_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), -dd);
+ mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+ }
+ scm_remember_upto_here_1 (n);
+ return scm_i_normbig (q);
+ }
+ }
+ else if (SCM_LIKELY (SCM_BIGP (d)))
+ {
+ SCM q = scm_i_mkbig ();
+ mpz_divexact (SCM_I_BIG_MPZ (q),
+ SCM_I_BIG_MPZ (n),
+ SCM_I_BIG_MPZ (d));
+ scm_remember_upto_here_2 (n, d);
+ return scm_i_normbig (q);
+ }
+ else
+ SCM_WRONG_TYPE_ARG (2, d);
+ }
+ else
+ SCM_WRONG_TYPE_ARG (1, n);
+}
+#undef FUNC_NAME
+
/* two_valued_wta_dispatch_2 is a version of SCM_WTA_DISPATCH_2 for
two-valued functions. It is called from primitive generics that take
two arguments and return two values, when the core procedure is
@@ -4636,6 +4676,23 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
else /* return NaN for (0 ^ k) for negative k per R6RS */
return scm_nan ();
}
+ /* Optimize the fraction case, to avoid needless reduction of
+ intermediate products to lowest terms. There will never be any
+ common factor to cancel out, so it would be a waste. */
+ else if (SCM_FRACTIONP (n))
+ {
+ if (scm_is_true (scm_positive_p (k)))
+ return scm_i_make_ratio_already_reduced
+ (scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k),
+ scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k));
+ else
+ {
+ k = scm_difference (k, SCM_UNDEFINED);
+ return scm_i_make_ratio_already_reduced
+ (scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k),
+ scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k));
+ }
+ }
if (SCM_I_INUMP (k))
i2 = SCM_I_INUM (k);
@@ -7282,8 +7339,9 @@ scm_difference (SCM x, SCM y)
return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x),
-SCM_COMPLEX_IMAG (x));
else if (SCM_FRACTIONP (x))
- return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x),
SCM_UNDEFINED),
- SCM_FRACTION_DENOMINATOR (x));
+ return scm_i_make_ratio_already_reduced
+ (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
+ SCM_FRACTION_DENOMINATOR (x));
else
SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference);
}
@@ -7825,14 +7883,14 @@ do_divide (SCM x, SCM y, int inexact)
{
if (inexact)
return scm_from_double (1.0 / (double) xx);
- else return scm_i_make_ratio (SCM_INUM1, x);
+ else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
}
}
else if (SCM_BIGP (x))
{
if (inexact)
return scm_from_double (1.0 / scm_i_big2dbl (x));
- else return scm_i_make_ratio (SCM_INUM1, x);
+ else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
}
else if (SCM_REALP (x))
{
@@ -7862,8 +7920,8 @@ do_divide (SCM x, SCM y, int inexact)
}
}
else if (SCM_FRACTIONP (x))
- return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x),
- SCM_FRACTION_NUMERATOR (x));
+ return scm_i_make_ratio_already_reduced (SCM_FRACTION_DENOMINATOR (x),
+ SCM_FRACTION_NUMERATOR (x));
else
SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide);
}
@@ -7960,32 +8018,9 @@ do_divide (SCM x, SCM y, int inexact)
return x;
else
{
- /* FIXME: HMM, what are the relative performance issues here?
- We need to test. Is it faster on average to test
- divisible_p, then perform whichever operation, or is it
- faster to perform the integer div opportunistically and
- switch to real if there's a remainder? For now we take the
- middle ground: test, then if divisible, use the faster div
- func. */
-
- scm_t_inum abs_yy = yy < 0 ? -yy : yy;
- int divisible_p = mpz_divisible_ui_p (SCM_I_BIG_MPZ (x), abs_yy);
-
- if (divisible_p)
- {
- SCM result = scm_i_mkbig ();
- mpz_divexact_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x),
abs_yy);
- scm_remember_upto_here_1 (x);
- if (yy < 0)
- mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
- return scm_i_normbig (result);
- }
- else
- {
- if (inexact)
- return scm_from_double (scm_i_big2dbl (x) / (double) yy);
- else return scm_i_make_ratio (x, y);
- }
+ if (inexact)
+ return scm_from_double (scm_i_big2dbl (x) / (double) yy);
+ else return scm_i_make_ratio (x, y);
}
}
else if (SCM_BIGP (y))
@@ -8003,21 +8038,7 @@ do_divide (SCM x, SCM y, int inexact)
return scm_from_double (mpq_get_d (q));
}
else
- {
- int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
- if (divisible_p)
- {
- SCM result = scm_i_mkbig ();
- mpz_divexact (SCM_I_BIG_MPZ (result),
- SCM_I_BIG_MPZ (x),
- SCM_I_BIG_MPZ (y));
- scm_remember_upto_here_2 (x, y);
- return scm_i_normbig (result);
- }
- else
- return scm_i_make_ratio (x, y);
- }
+ return scm_i_make_ratio (x, y);
}
else if (SCM_REALP (y))
{
@@ -8826,8 +8847,9 @@ SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0,
0,
{
if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
return z;
- return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z),
SCM_UNDEFINED),
- SCM_FRACTION_DENOMINATOR (z));
+ return scm_i_make_ratio_already_reduced
+ (scm_difference (SCM_FRACTION_NUMERATOR (z), SCM_UNDEFINED),
+ SCM_FRACTION_DENOMINATOR (z));
}
else
SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude);
@@ -8927,8 +8949,9 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact,
"inexact->exact", 1, 0, 0,
mpq_init (frac);
mpq_set_d (frac, val);
- q = scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
- scm_i_mpz2num (mpq_denref (frac)));
+ q = scm_i_make_ratio_already_reduced
+ (scm_i_mpz2num (mpq_numref (frac)),
+ scm_i_mpz2num (mpq_denref (frac)));
/* When scm_i_make_ratio throws, we leak the memory allocated
for frac...
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index b1f3d8b..7cb4694 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4034,6 +4034,9 @@
(pass-if (eqv? -0.125 (expt -2 -3.0)))
(pass-if (eqv? -0.125 (expt -2.0 -3.0)))
(pass-if (eqv? 0.25 (expt 2.0 -2.0)))
+ (pass-if (eqv? 32/243 (expt 2/3 5)))
+ (pass-if (eqv? 243/32 (expt 2/3 -5)))
+ (pass-if (eqv? 32 (expt 1/2 -5)))
(pass-if (test-eqv? (* -1.0+0.0i 12398 12398) (expt +12398i 2.0)))
(pass-if (eqv-loosely? +i (expt -1 0.5)))
(pass-if (eqv-loosely? +i (expt -1 1/2)))
@@ -4304,6 +4307,9 @@
(pass-if (eqv? -1/8 (integer-expt -2 -3)))
(pass-if (eqv? -0.125 (integer-expt -2.0 -3)))
(pass-if (eqv? 0.25 (integer-expt 2.0 -2)))
+ (pass-if (eqv? 32/243 (integer-expt 2/3 5)))
+ (pass-if (eqv? 243/32 (integer-expt 2/3 -5)))
+ (pass-if (eqv? 32 (integer-expt 1/2 -5)))
(pass-if (test-eqv? (* -1.0+0.0i 12398 12398) (integer-expt +12398.0i 2))))
--
1.7.2.5
>From c1e2e76872dad0fb53c4f389b387019b136e36e1 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Wed, 6 Apr 2011 17:29:06 -0400
Subject: [PATCH 3/7] Add `round-ash', a rounding arithmetic shift operator
* libguile/numbers.c (left_shift_exact_integer,
floor_right_shift_exact_integer, round_right_shift_exact_integer):
New internal static functions to efficiently compute n * 2^count,
floor(n / 2^count), and round(n / 2^count), respectively, for
exact integer N and positive exact integer COUNT. Used to combine
the implementations of `ash' and `round-ash' with minimal code
duplication.
(scm_ash): Reimplement in terms of left_shift_exact_integer and
floor_right_shift_exact_integer. Note that this function efficiently
computes floor(n * 2^count).
(scm_round_ash): New procedure to efficiently compute
round(n * 2^count), for exact integers N and COUNT. Implemented
in terms of left_shift_exact_integer and
round_right_shift_exact_integer.
* libguile/numbers.h: Add prototype for scm_round_ash. Change the
name of scm_ash's second formal parameter from `cnt' to `count'.
* test-suite/tests/numbers.test (round-ash, ash): Add new unified
testing framework for `ash' and `round-ash'. Previously, the tests
for `ash' were not very comprehensive; for example, they did not
include a single test where the number to be shifted was a bignum.
* doc/ref/api-data.texi (Bitwise Operations): Add documentation for
`round-ash'. Improve documentation for `ash'.
* NEWS: Add NEWS entry.
---
NEWS | 8 ++
doc/ref/api-data.texi | 44 ++++++--
libguile/numbers.c | 230 ++++++++++++++++++++++++++++-------------
libguile/numbers.h | 3 +-
test-suite/tests/numbers.test | 114 +++++++++------------
5 files changed, 247 insertions(+), 152 deletions(-)
diff --git a/NEWS b/NEWS
index f5760cf..68168b0 100644
--- a/NEWS
+++ b/NEWS
@@ -4,7 +4,15 @@ See the end for copying conditions.
Please send Guile bug reports to address@hidden
+Changes in 2.0.3 (since 2.0.2):
+* Notable changes
+
+** New procedure: `round-ash', a rounding arithmetic shift operator
+
+See "Bitwise Operations" in the manual, for more.
+
+
Changes in 2.0.2 (since 2.0.1):
* Notable changes
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index 9825bef..3a14d29 100644
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1671,19 +1671,16 @@ starts from 0 for the least significant bit.
@end lisp
@end deffn
address@hidden {Scheme Procedure} ash n cnt
address@hidden {C Function} scm_ash (n, cnt)
-Return @var{n} shifted left by @var{cnt} bits, or shifted right if
address@hidden is negative. This is an ``arithmetic'' shift.
-
-This is effectively a multiplication by @m{2^{cnt}, address@hidden, and
-when @var{cnt} is negative it's a division, rounded towards negative
-infinity. (Note that this is not the same rounding as @code{quotient}
-does.)
address@hidden {Scheme Procedure} ash n count
address@hidden {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * address@hidden)}, but computed
+more efficiently. @var{n} and @var{count} must be exact
+integers.
-With @var{n} viewed as an infinite precision twos complement,
address@hidden means a left shift introducing zero bits, or a right shift
-dropping bits.
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{ash} means a left shift introducing zero bits
+when @var{count} is positive, or a right shift dropping bits
+when @var{count} is negative. This is an ``arithmetic'' shift.
@lisp
(number->string (ash #b1 3) 2) @result{} "1000"
@@ -1694,6 +1691,29 @@ dropping bits.
@end lisp
@end deffn
address@hidden {Scheme Procedure} round-ash n count
address@hidden {C Function} scm_round_ash (n, count)
+Return @math{round(@var{n} * address@hidden)}, but computed
+more efficiently. @var{n} and @var{count} must be exact
+integers.
+
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{round-ash} means a left shift introducing zero
+bits when @var{count} is positive, or a right shift rounding
+to the nearest integer (with ties going to the nearest even
+integer) when @var{count} is negative. This is a rounded
+``arithmetic'' shift.
+
address@hidden
+(number->string (round-ash #b1 3) 2) @result{} \"1000\"
+(number->string (round-ash #b1010 -1) 2) @result{} \"101\"
+(number->string (round-ash #b1010 -2) 2) @result{} \"10\"
+(number->string (round-ash #b1011 -2) 2) @result{} \"11\"
+(number->string (round-ash #b1101 -2) 2) @result{} \"11\"
+(number->string (round-ash #b1110 -2) 2) @result{} \"100\"
address@hidden lisp
address@hidden deffn
+
@deffn {Scheme Procedure} logcount n
@deffnx {C Function} scm_logcount (n)
Return the number of bits in integer @var{n}. If @var{n} is
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 80b46df..e2c74d2 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4750,19 +4750,122 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
}
#undef FUNC_NAME
+/* n must be an exact integer, and count > 0.
+ Returns n * 2^count. */
+static SCM
+left_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+
+ /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will always
+ overflow a non-zero fixnum. For smaller shifts we check the
+ bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
+ all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
+ Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
+ count)". */
+
+ if (nn == 0)
+ return n;
+ else if (count < SCM_I_FIXNUM_BIT-1 &&
+ ((scm_t_bits) (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - count)) + 1)
+ <= 1))
+ {
+ return SCM_I_MAKINUM (nn << count);
+ }
+ else
+ {
+ SCM result = scm_i_inum2big (nn);
+ mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
+ count);
+ return result;
+ }
+ }
+ else if (SCM_BIGP (n))
+ {
+ SCM result = scm_i_mkbig ();
+ mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
+ count);
+ scm_remember_upto_here_1 (n);
+ return result;
+ }
+ else
+ scm_syserror ("left_shift_exact_integer");
+}
+
+/* n must be an exact integer, and count > 0.
+ Returns floor(n / 2^count). */
+static SCM
+floor_right_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+
+ if (count >= SCM_I_FIXNUM_BIT)
+ return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM (-1));
+ else
+ return SCM_I_MAKINUM (SCM_SRS (nn, count));
+ }
+ else if (SCM_BIGP (n))
+ {
+ SCM result = scm_i_mkbig ();
+ mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
+ count);
+ scm_remember_upto_here_1 (n);
+ return scm_i_normbig (result);
+ }
+ else
+ scm_syserror ("floor_right_shift_exact_integer");
+}
+
+/* n must be an exact integer, and count > 0.
+ Returns round(n / 2^count). */
+static SCM
+round_right_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+ scm_t_inum qq = SCM_SRS (nn, count);
+
+ if (count >= SCM_I_FIXNUM_BIT)
+ return SCM_INUM0;
+ else if (0 == (nn & (1L << (count-1))))
+ return SCM_I_MAKINUM (qq); /* round down */
+ else if (nn & ((1L << (count-1)) - 1))
+ return SCM_I_MAKINUM (qq + 1); /* round up */
+ else
+ return SCM_I_MAKINUM ((~1L) & (qq + 1)); /* round to even */
+ }
+ else if (SCM_BIGP (n))
+ {
+ SCM q = scm_i_mkbig ();
+ mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), count);
+ if (mpz_tstbit (SCM_I_BIG_MPZ (n), count-1))
+ {
+ if ((mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1) ||
+ (mpz_odd_p (SCM_I_BIG_MPZ (q))))
+ mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
+ }
+ scm_remember_upto_here_1 (n);
+ return scm_i_normbig (q);
+ }
+ else
+ scm_syserror ("round_right_shift_exact_integer");
+}
+
SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
- (SCM n, SCM cnt),
- "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
- "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
- "\n"
- "This is effectively a multiplication by address@hidden, and when\n"
- "@var{cnt} is negative it's a division, rounded towards negative\n"
- "infinity. (Note that this is not the same rounding as\n"
- "@code{quotient} does.)\n"
+ (SCM n, SCM count),
+ "Return @math{floor(@var{n} * address@hidden)}, but computed\n"
+ "more efficiently. @var{n} and @var{count} must be exact\n"
+ "integers.\n"
"\n"
- "With @var{n} viewed as an infinite precision twos complement,\n"
- "@code{ash} means a left shift introducing zero bits, or a right\n"
- "shift dropping bits.\n"
+ "With @var{n} viewed as an infinite-precision twos-complement\n"
+ "integer, @code{ash} means a left shift introducing zero bits\n"
+ "when @var{count} is positive, or a right shift dropping bits\n"
+ "when @var{count} is negative. This is an ``arithmetic'' shift.\n"
"\n"
"@lisp\n"
"(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
@@ -4773,79 +4876,58 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
"@end lisp")
#define FUNC_NAME s_scm_ash
{
- long bits_to_shift;
- bits_to_shift = scm_to_long (cnt);
-
- if (SCM_I_INUMP (n))
+ if (SCM_I_INUMP (n) || SCM_BIGP (n))
{
- scm_t_inum nn = SCM_I_INUM (n);
+ long bits_to_shift = scm_to_long (count);
if (bits_to_shift > 0)
- {
- /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
- overflow a non-zero fixnum. For smaller shifts we check the
- bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
- all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
- Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
- bits_to_shift)". */
-
- if (nn == 0)
- return n;
-
- if (bits_to_shift < SCM_I_FIXNUM_BIT-1
- && ((scm_t_bits)
- (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1)
- <= 1))
- {
- return SCM_I_MAKINUM (nn << bits_to_shift);
- }
- else
- {
- SCM result = scm_i_inum2big (nn);
- mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
- bits_to_shift);
- return result;
- }
- }
+ return left_shift_exact_integer (n, bits_to_shift);
+ else if (SCM_LIKELY (bits_to_shift < 0))
+ return floor_right_shift_exact_integer (n, -bits_to_shift);
else
- {
- bits_to_shift = -bits_to_shift;
- if (bits_to_shift >= SCM_LONG_BIT)
- return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM(-1));
- else
- return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift));
- }
-
+ return n;
}
- else if (SCM_BIGP (n))
- {
- SCM result;
+ else
+ SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
+}
+#undef FUNC_NAME
- if (bits_to_shift == 0)
- return n;
+SCM_DEFINE (scm_round_ash, "round-ash", 2, 0, 0,
+ (SCM n, SCM count),
+ "Return @math{round(@var{n} * address@hidden)}, but computed\n"
+ "more efficiently. @var{n} and @var{count} must be exact\n"
+ "integers.\n"
+ "\n"
+ "With @var{n} viewed as an infinite-precision twos-complement\n"
+ "integer, @code{round-ash} means a left shift introducing zero\n"
+ "bits when @var{count} is positive, or a right shift rounding\n"
+ "to the nearest integer (with ties going to the nearest even\n"
+ "integer) when @var{count} is negative. This is a rounded\n"
+ "``arithmetic'' shift.\n"
+ "\n"
+ "@lisp\n"
+ "(number->string (round-ash #b1 3) 2) @result{} \"1000\"\n"
+ "(number->string (round-ash #b1010 -1) 2) @result{} \"101\"\n"
+ "(number->string (round-ash #b1010 -2) 2) @result{} \"10\"\n"
+ "(number->string (round-ash #b1011 -2) 2) @result{} \"11\"\n"
+ "(number->string (round-ash #b1101 -2) 2) @result{} \"11\"\n"
+ "(number->string (round-ash #b1110 -2) 2) @result{} \"100\"\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_round_ash
+{
+ if (SCM_I_INUMP (n) || SCM_BIGP (n))
+ {
+ long bits_to_shift = scm_to_long (count);
- result = scm_i_mkbig ();
- if (bits_to_shift >= 0)
- {
- mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
- bits_to_shift);
- return result;
- }
+ if (bits_to_shift > 0)
+ return left_shift_exact_integer (n, bits_to_shift);
+ else if (SCM_LIKELY (bits_to_shift < 0))
+ return round_right_shift_exact_integer (n, -bits_to_shift);
else
- {
- /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
- we have to allocate a bignum even if the result is going to be a
- fixnum. */
- mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
- -bits_to_shift);
- return scm_i_normbig (result);
- }
-
+ return n;
}
else
- {
- SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
- }
+ SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
}
#undef FUNC_NAME
diff --git a/libguile/numbers.h b/libguile/numbers.h
index d985830..a0cc8db 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -204,7 +204,8 @@ SCM_API SCM scm_logbit_p (SCM n1, SCM n2);
SCM_API SCM scm_lognot (SCM n);
SCM_API SCM scm_modulo_expt (SCM n, SCM k, SCM m);
SCM_API SCM scm_integer_expt (SCM z1, SCM z2);
-SCM_API SCM scm_ash (SCM n, SCM cnt);
+SCM_API SCM scm_ash (SCM n, SCM count);
+SCM_API SCM scm_round_ash (SCM n, SCM count);
SCM_API SCM scm_bit_extract (SCM n, SCM start, SCM end);
SCM_API SCM scm_logcount (SCM n);
SCM_API SCM scm_integer_length (SCM n);
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 7cb4694..c4979fd 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -201,71 +201,6 @@
(eqv? -2305843009213693953 (1- -2305843009213693952))))
;;;
-;;; ash
-;;;
-
-(with-test-prefix "ash"
-
- (pass-if "documented?"
- (documented? ash))
-
- (pass-if (eqv? 0 (ash 0 0)))
- (pass-if (eqv? 0 (ash 0 1)))
- (pass-if (eqv? 0 (ash 0 1000)))
- (pass-if (eqv? 0 (ash 0 -1)))
- (pass-if (eqv? 0 (ash 0 -1000)))
-
- (pass-if (eqv? 1 (ash 1 0)))
- (pass-if (eqv? 2 (ash 1 1)))
- (pass-if (eqv? 340282366920938463463374607431768211456 (ash 1 128)))
- (pass-if (eqv? 0 (ash 1 -1)))
- (pass-if (eqv? 0 (ash 1 -1000)))
-
- (pass-if (eqv? -1 (ash -1 0)))
- (pass-if (eqv? -2 (ash -1 1)))
- (pass-if (eqv? -340282366920938463463374607431768211456 (ash -1 128)))
- (pass-if (eqv? -1 (ash -1 -1)))
- (pass-if (eqv? -1 (ash -1 -1000)))
-
- (pass-if (eqv? -3 (ash -3 0)))
- (pass-if (eqv? -6 (ash -3 1)))
- (pass-if (eqv? -1020847100762815390390123822295304634368 (ash -3 128)))
- (pass-if (eqv? -2 (ash -3 -1)))
- (pass-if (eqv? -1 (ash -3 -1000)))
-
- (pass-if (eqv? -6 (ash -23 -2)))
-
- (pass-if (eqv? most-positive-fixnum (ash most-positive-fixnum 0)))
- (pass-if (eqv? (* 2 most-positive-fixnum) (ash most-positive-fixnum 1)))
- (pass-if (eqv? (* 4 most-positive-fixnum) (ash most-positive-fixnum 2)))
- (pass-if
- (eqv? (* most-positive-fixnum 340282366920938463463374607431768211456)
- (ash most-positive-fixnum 128)))
- (pass-if (eqv? (quotient most-positive-fixnum 2)
- (ash most-positive-fixnum -1)))
- (pass-if (eqv? 0 (ash most-positive-fixnum -1000)))
-
- (let ((mpf4 (quotient most-positive-fixnum 4)))
- (pass-if (eqv? (* 2 mpf4) (ash mpf4 1)))
- (pass-if (eqv? (* 4 mpf4) (ash mpf4 2)))
- (pass-if (eqv? (* 8 mpf4) (ash mpf4 3))))
-
- (pass-if (eqv? most-negative-fixnum (ash most-negative-fixnum 0)))
- (pass-if (eqv? (* 2 most-negative-fixnum) (ash most-negative-fixnum 1)))
- (pass-if (eqv? (* 4 most-negative-fixnum) (ash most-negative-fixnum 2)))
- (pass-if
- (eqv? (* most-negative-fixnum 340282366920938463463374607431768211456)
- (ash most-negative-fixnum 128)))
- (pass-if (eqv? (quotient-floor most-negative-fixnum 2)
- (ash most-negative-fixnum -1)))
- (pass-if (eqv? -1 (ash most-negative-fixnum -1000)))
-
- (let ((mnf4 (quotient-floor most-negative-fixnum 4)))
- (pass-if (eqv? (* 2 mnf4) (ash mnf4 1)))
- (pass-if (eqv? (* 4 mnf4) (ash mnf4 2)))
- (pass-if (eqv? (* 8 mnf4) (ash mnf4 3)))))
-
-;;;
;;; exact?
;;;
@@ -4891,3 +4826,52 @@
round-quotient
round-remainder
valid-round-answer?)))
+
+;;;
+;;; ash
+;;; round-ash
+;;;
+
+(let ()
+ (define (test-ash-variant name ash-variant round-variant)
+ (with-test-prefix name
+ (define (test n count)
+ (pass-if (list n count)
+ (eqv? (ash-variant n count)
+ (round-variant (* n (expt 2 count))))))
+
+ (pass-if "documented?"
+ (documented? ash-variant))
+
+ (for-each (lambda (n)
+ (for-each (lambda (count) (test n count))
+ '(-1000 -3 -2 -1 0 1 2 3 1000)))
+ (list 0 1 3 23 -1 -3 -23
+ fixnum-max
+ (1+ fixnum-max)
+ (1- fixnum-max)
+ (* fixnum-max 4)
+ (quotient fixnum-max 4)
+ fixnum-min
+ (1+ fixnum-min)
+ (1- fixnum-min)
+ (* fixnum-min 4)
+ (quotient fixnum-min 4)))
+
+ (do ((count -2 (1- count))
+ (vals '(1 3 5 7 9 11)
+ (map (lambda (n) (* 2 n)) vals)))
+ ((> (car vals) (* 2 fixnum-max)) 'done)
+ (for-each (lambda (n)
+ (test n count)
+ (test (- n) count))
+ vals))
+
+ ;; Test rounding
+ (for-each (lambda (base)
+ (for-each (lambda (offset) (test (+ base offset) -3))
+ '(#b11001 #b11100 #b11101 #b10001 #b10100
#b10101)))
+ (list 0 64 -64 (* 64 fixnum-max) (* 64 fixnum-min)))))
+
+ (test-ash-variant 'ash ash floor)
+ (test-ash-variant 'round-ash round-ash round))
--
1.7.2.5
>From 42ecb421749b7b5ab64943c402c92587509e9f04 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Wed, 23 Feb 2011 00:04:12 -0500
Subject: [PATCH 4/7] Introduce double_precision and use it
* libguile/numbers.c (double_precision): New internal static global
which contains the number of bits of precision in the C double type.
(scm_i_dbl2num, scm_num_eq_p): Use double_precision instead of
DBL_MANT_DIG. Previously, the code improperly assumed that
DBL_MANT_DIG is a number of bits, when in fact it is a number of
base-`FLT_RADIX' digits.
(log_of_shifted_double): Use double_precision instead of
scm_dblprec[0]. (Note that scm_dblprec[0] is used by the idbl2str,
and is either 60 or the number of bits of precision of doubles,
whichever is less. This limitation cannot be lifted without
increasing the size of the fx_per_radix array, which is already over
16 kilobytes.)
---
libguile/numbers.c | 64 ++++++++++++++++++++++++++++++++++++---------------
1 files changed, 45 insertions(+), 19 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index e2c74d2..1274306 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -168,10 +168,40 @@ scm_from_complex_double (complex double z)
#endif /* GUILE_I */
+/* double_precision is the smallest positive integer k
+ such that 1.0 == 1.0 + 2.0^(-k).
-static mpz_t z_negative_one;
+ For floating-point numbers with significands represented
+ in base 2, it is the number of bits in the significand
+ (including the most-significant 1 bit which is sometimes
+ left implicit and not actually stored).
+
+ When the C "double" type corresponds to the IEEE 754
+ binary64 format, double_precision should be 53.
+*/
+static long double_precision;
+
+static void
+init_double_precision (void)
+{
+ double sum;
+ long i;
+
+ for (i = 1; i < 1000000; i++)
+ {
+ sum = 1.0 + ldexp (1.0, -i);
+ if (sum == 1.0)
+ {
+ double_precision = i;
+ return;
+ }
+ }
+ scm_syserror ("init_double_precision");
+}
+static mpz_t z_negative_one;
+
/* Clear the `mpz_t' embedded in bignum PTR. */
static void
finalize_bignum (GC_PTR ptr, GC_PTR data)
@@ -304,10 +334,10 @@ scm_i_dbl2num (double u)
/* scm_i_big2dbl() rounds to the closest representable double, in accordance
with R5RS exact->inexact.
- The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
- (ie. truncate towards zero), then adjust to get the closest double by
- examining the next lower bit and adding 1 (to the absolute value) if
- necessary.
+ The approach is to use mpz_get_d to pick out the high
+ double_precision bits (ie. truncate towards zero), then adjust to
+ get the closest double by examining the next lower bit and adding 1
+ (to the absolute value) if necessary.
Bignums exactly half way between representable doubles are rounded to the
next higher absolute value (ie. away from zero). This seems like an
@@ -321,7 +351,7 @@ scm_i_dbl2num (double u)
In GMP before 4.2, mpz_get_d rounding was unspecified. It ended up
following the hardware rounding mode, but applied to the absolute
value of the mpz_t operand. This is not what we want so we put the
- high DBL_MANT_DIG bits into a temporary. Starting with GMP 4.2
+ high double_precision bits into a temporary. Starting with GMP 4.2
(released in March 2006) mpz_get_d now always truncates towards zero.
ENHANCE-ME: The temporary init+clear to force the rounding in GMP
@@ -339,16 +369,11 @@ scm_i_big2dbl (SCM b)
#if 1
{
/* For GMP earlier than 4.2, force truncation towards zero */
-
- /* FIXME: DBL_MANT_DIG is the number of base-`FLT_RADIX' digits,
- _not_ the number of bits, so this code will break badly on a
- system with non-binary doubles. */
-
mpz_t tmp;
- if (bits > DBL_MANT_DIG)
+ if (bits > double_precision)
{
- size_t shift = bits - DBL_MANT_DIG;
- mpz_init2 (tmp, DBL_MANT_DIG);
+ size_t shift = bits - double_precision;
+ mpz_init2 (tmp, double_precision);
mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift);
result = ldexp (mpz_get_d (tmp), shift);
mpz_clear (tmp);
@@ -363,9 +388,9 @@ scm_i_big2dbl (SCM b)
result = mpz_get_d (SCM_I_BIG_MPZ (b));
#endif
- if (bits > DBL_MANT_DIG)
+ if (bits > double_precision)
{
- unsigned long pos = bits - DBL_MANT_DIG - 1;
+ unsigned long pos = bits - double_precision - 1;
/* test bit number "pos" in absolute value */
if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS)
& ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS)))
@@ -6332,7 +6357,7 @@ scm_num_eq_p (SCM x, SCM y)
double yy = SCM_REAL_VALUE (y);
return scm_from_bool ((double) xx == yy
- && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
+ && (double_precision >= SCM_I_FIXNUM_BIT-1
|| xx == (scm_t_signed_bits) yy));
}
else if (SCM_COMPLEXP (y))
@@ -6386,7 +6411,7 @@ scm_num_eq_p (SCM x, SCM y)
/* see comments with inum/real above */
scm_t_signed_bits yy = SCM_I_INUM (y);
return scm_from_bool (xx == (double) yy
- && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
+ && (double_precision >= SCM_I_FIXNUM_BIT-1
|| (scm_t_signed_bits) xx == yy));
}
else if (SCM_BIGP (y))
@@ -9519,7 +9544,7 @@ log_of_shifted_double (double x, long shift)
static SCM
log_of_exact_integer_with_size (SCM n, long size)
{
- long shift = size - 2 * scm_dblprec[0];
+ long shift = size - 2 * double_precision;
if (shift > 0)
return log_of_shifted_double
@@ -9806,6 +9831,7 @@ scm_init_numbers ()
flo_log10e = scm_from_double (M_LOG10E);
/* determine floating point precision */
+ init_double_precision ();
for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
{
init_dblprec(&scm_dblprec[i-2],i);
--
1.7.2.5
>From 54f194f285114ad046fb67fd3027053f9544080f Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Fri, 4 Mar 2011 15:28:51 -0500
Subject: [PATCH 5/7] Simplify and improve scm_i_big2dbl, and add
scm_i_big2dbl_2exp
* libguile/numbers.c (scm_i_big2dbl_2exp): New internal function which
converts a bignum into a normalized significand and an exponent,
similar to frexp, or mpz_get_d_2exp but with correct rounding. This
can be used on bignums that are too large or small to be represented
by a double.
(scm_i_big2dbl): Simplify and improve scm_i_big2dbl, which is now
implemented in terms of scm_i_big2dbl_2exp. Ties now go to even.
Previously they were rounded away from zero.
---
libguile/numbers.c | 94 ++++++++++++++++-----------------------------------
1 files changed, 30 insertions(+), 64 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 1274306..1773715 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -331,76 +331,42 @@ scm_i_dbl2num (double u)
return scm_i_dbl2big (u);
}
-/* scm_i_big2dbl() rounds to the closest representable double, in accordance
- with R5RS exact->inexact.
-
- The approach is to use mpz_get_d to pick out the high
- double_precision bits (ie. truncate towards zero), then adjust to
- get the closest double by examining the next lower bit and adding 1
- (to the absolute value) if necessary.
-
- Bignums exactly half way between representable doubles are rounded to the
- next higher absolute value (ie. away from zero). This seems like an
- adequate interpretation of R5RS "numerically closest", and it's easier
- and faster than a full "nearest-even" style.
-
- The bit test must be done on the absolute value of the mpz_t, which means
- we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
- negatives as twos complement.
-
- In GMP before 4.2, mpz_get_d rounding was unspecified. It ended up
- following the hardware rounding mode, but applied to the absolute
- value of the mpz_t operand. This is not what we want so we put the
- high double_precision bits into a temporary. Starting with GMP 4.2
- (released in March 2006) mpz_get_d now always truncates towards zero.
-
- ENHANCE-ME: The temporary init+clear to force the rounding in GMP
- before 4.2 is a slowdown. It'd be faster to pick out the relevant
- high bits with mpz_getlimbn. */
-
-double
-scm_i_big2dbl (SCM b)
-{
- double result;
- size_t bits;
-
- bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
-
-#if 1
- {
- /* For GMP earlier than 4.2, force truncation towards zero */
- mpz_t tmp;
- if (bits > double_precision)
- {
- size_t shift = bits - double_precision;
- mpz_init2 (tmp, double_precision);
- mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift);
- result = ldexp (mpz_get_d (tmp), shift);
- mpz_clear (tmp);
- }
- else
- {
- result = mpz_get_d (SCM_I_BIG_MPZ (b));
- }
- }
-#else
- /* GMP 4.2 or later */
- result = mpz_get_d (SCM_I_BIG_MPZ (b));
-#endif
+static SCM round_right_shift_exact_integer (SCM n, long count);
+
+/* scm_i_big2dbl_2exp() is like frexp for bignums: it converts the
+ bignum b into a normalized significand and exponent such that
+ b = significand * 2^exponent and 1/2 <= abs(significand) < 1.
+ The return value is the significand rounded to the closest
+ representable double, and the exponent is placed into *expon_p.
+ If b is zero, then the returned exponent and significand are both
+ zero. */
+static double
+scm_i_big2dbl_2exp (SCM b, long *expon_p)
+{
+ size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
+ size_t shift = 0;
+ double signif;
+ long expon;
if (bits > double_precision)
{
- unsigned long pos = bits - double_precision - 1;
- /* test bit number "pos" in absolute value */
- if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS)
- & ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS)))
- {
- result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
- }
+ shift = bits - double_precision;
+ b = round_right_shift_exact_integer (b, shift);
}
-
+ signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b));
scm_remember_upto_here_1 (b);
- return result;
+ *expon_p = expon + shift;
+ return signif;
+}
+
+/* scm_i_big2dbl() rounds to the closest representable double,
+ in accordance with R5RS exact->inexact. */
+double
+scm_i_big2dbl (SCM b)
+{
+ long expon;
+ double signif = scm_i_big2dbl_2exp (b, &expon);
+ return ldexp (signif, expon);
}
SCM
--
1.7.2.5
>From b2f2efbfa10dcb00d59cb4567ba5b6a217b1ed81 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Wed, 23 Feb 2011 02:09:34 -0500
Subject: [PATCH 6/7] Optimize logarithms using scm_i_big2dbl_2exp
* libguile/numbers.c (log_of_exact_integer_with_size): Removed.
(log_of_exact_integer): Handle bignums too large to fit in a double
using scm_i_big2dbl_2exp instead of scm_integer_length and scm_ash.
(log_of_fraction): Use log_of_exact_integer instead of
log_of_exact_integer_with_size.
---
libguile/numbers.c | 30 ++++++++++++------------------
1 files changed, 12 insertions(+), 18 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 1773715..c9be26d 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9506,26 +9506,20 @@ log_of_shifted_double (double x, long shift)
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 * double_precision;
-
- 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 */
static SCM
log_of_exact_integer (SCM n)
{
- return log_of_exact_integer_with_size
- (n, scm_to_long (scm_integer_length (n)));
+ if (SCM_I_INUMP (n))
+ return log_of_shifted_double (SCM_I_INUM (n), 0);
+ else if (SCM_BIGP (n))
+ {
+ long expon;
+ double signif = scm_i_big2dbl_2exp (n, &expon);
+ return log_of_shifted_double (signif, expon);
+ }
+ else
+ scm_wrong_type_arg ("log_of_exact_integer", SCM_ARG1, n);
}
/* Returns log(n/d), for exact non-zero integers n and d */
@@ -9536,8 +9530,8 @@ log_of_fraction (SCM n, SCM d)
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)));
+ return (scm_difference (log_of_exact_integer (n),
+ log_of_exact_integer (d)));
else if (scm_is_false (scm_negative_p (n)))
return scm_from_double
(log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
--
1.7.2.5
>From d754e3f1de92fa2d33edc610f37ff24b58cddbcd Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Fri, 4 Mar 2011 17:40:56 -0500
Subject: [PATCH 7/7] Improve inexact division of exact integers et al
* libguile/numbers.c (scm_to_double_2exp): New internal procedure which
converts an arbitrary real number to a normalized C double significand
and a base-2 exponent, similar to frexp. This is the basis for making
several other inexact numerical operations robust when working with
large exact integers and rationals.
(scm_divide2real): Removed; superceded by scm_i_divide2double.
(scm_i_divide2double): Returns the ratio of two arbitrary real number
objects as a C double.
(scm_i_fraction2double, log_of_fraction): Use scm_i_divide2double
instead of scm_divide2real.
(scm_divide, do_divide): The core implementation is now in scm_divide.
do_divide has been removed. Previously, do_divide was the basis for
both scm_divide and scm_divide2real. Remove the `inexact' flag
argument, which (when non-zero) requested an inexact answer. The
inexact implementation was flawed in several respects, and is now
done by scm_i_divide2double.
(scm_i_inum_length): New internal procedure that handles the inum case
of scm_integer_length and returns the result as a C int. Used by the
inum case of scm_to_double_2exp.
(scm_integer_length): Moved inum case to scm_i_inum_length.
---
libguile/numbers.c | 170 ++++++++++++++++++++++++++++-----------------------
1 files changed, 93 insertions(+), 77 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index c9be26d..b35f66c 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -401,9 +401,6 @@ scm_i_mpz2num (mpz_t b)
}
}
-/* this is needed when we want scm_divide to make a float, not a ratio, even
if passed two ints */
-static SCM scm_divide2real (SCM x, SCM y);
-
/* scm_i_make_ratio_already_reduced makes a ratio more efficiently,
when the following conditions are known in advance:
@@ -459,11 +456,74 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
}
#undef FUNC_NAME
+static int scm_i_inum_length (scm_t_inum n);
+
+/* scm_to_double_2exp() an inexact frexp for SCM real numbers: it
+ converts x into a normalized significand and exponent such that
+ x = significand * 2^exponent and 1/2 <= abs(significand) < 1.
+ The return value is the significand rounded to the closest
+ representable double, and the exponent is placed into *expon_p.
+ If x is zero, then the returned exponent and significand are both
+ zero. */
+static double
+scm_to_double_2exp (SCM x, long *expon_p)
+{
+ if (SCM_I_INUMP (x))
+ {
+ int expon = scm_i_inum_length (abs (SCM_I_INUM (x)));
+ *expon_p = expon;
+ return ldexp (SCM_I_INUM (x), -expon);
+ }
+ else if (SCM_BIGP (x))
+ return scm_i_big2dbl_2exp (x, expon_p);
+ else if (SCM_REALP (x))
+ {
+ int expon;
+ double signif = frexp (SCM_REAL_VALUE (x), &expon);
+ *expon_p = expon;
+ return signif;
+ }
+ else if (SCM_FRACTIONP (x))
+ {
+ long n_expon, d_expon;
+ double n_signif = scm_to_double_2exp (SCM_FRACTION_NUMERATOR (x),
+ &n_expon);
+ double d_signif = scm_to_double_2exp (SCM_FRACTION_DENOMINATOR (x),
+ &d_expon);
+ long expon = n_expon - d_expon;
+ double signif = n_signif / d_signif;
+ if (signif >= 1.0 || signif <= -1.0)
+ {
+ signif /= 2.0;
+ expon++;
+ }
+ *expon_p = expon;
+ return signif;
+ }
+ else
+ scm_wrong_type_arg ("scm_to_double_2exp", SCM_ARG1, x);
+}
+
+/* Return n/d as a double, where n and d are exact integers. */
+static double
+scm_i_divide2double (SCM n, SCM d)
+{
+ if (SCM_LIKELY (SCM_I_INUMP (n) && SCM_I_INUMP (d)))
+ return ((double) SCM_I_INUM (n)) / ((double) SCM_I_INUM (d));
+ else
+ {
+ long n_expon, d_expon;
+ double n_signif = scm_to_double_2exp (n, &n_expon);
+ double d_signif = scm_to_double_2exp (d, &d_expon);
+ return ldexp (n_signif / d_signif, n_expon - d_expon);
+ }
+}
+
double
scm_i_fraction2double (SCM z)
{
- return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z),
- SCM_FRACTION_DENOMINATOR (z)));
+ return scm_i_divide2double (SCM_FRACTION_NUMERATOR (z),
+ SCM_FRACTION_DENOMINATOR (z));
}
static int
@@ -5049,6 +5109,22 @@ static const char scm_ilentab[] = {
0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
};
+/* Return the number of bits necessary to represent n */
+static int scm_i_inum_length (scm_t_inum n)
+{
+ unsigned long c = 0;
+ unsigned int l = 4;
+
+ if (n < 0)
+ n = -1 - n;
+ while (n)
+ {
+ c += 4;
+ l = scm_ilentab [15 & n];
+ n >>= 4;
+ }
+ return c - 4 + l;
+}
SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0,
(SCM n),
@@ -5065,20 +5141,7 @@ SCM_DEFINE (scm_integer_length, "integer-length", 1, 0,
0,
#define FUNC_NAME s_scm_integer_length
{
if (SCM_I_INUMP (n))
- {
- unsigned long c = 0;
- unsigned int l = 4;
- scm_t_inum nn = SCM_I_INUM (n);
- if (nn < 0)
- nn = -1 - nn;
- while (nn)
- {
- c += 4;
- l = scm_ilentab [15 & nn];
- nn >>= 4;
- }
- return SCM_I_MAKINUM (c - 4 + l);
- }
+ return SCM_I_MAKINUM (scm_i_inum_length (SCM_I_INUM (n)));
else if (SCM_BIGP (n))
{
/* mpz_sizeinbase looks at the absolute value of negatives, whereas we
@@ -7933,8 +7996,8 @@ SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1,
#define s_divide s_scm_i_divide
#define g_divide g_scm_i_divide
-static SCM
-do_divide (SCM x, SCM y, int inexact)
+SCM
+scm_divide (SCM x, SCM y)
#define FUNC_NAME s_divide
{
double a;
@@ -7953,18 +8016,10 @@ do_divide (SCM x, SCM y, int inexact)
scm_num_overflow (s_divide);
#endif
else
- {
- if (inexact)
- return scm_from_double (1.0 / (double) xx);
- else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
- }
+ return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
}
else if (SCM_BIGP (x))
- {
- if (inexact)
- return scm_from_double (1.0 / scm_i_big2dbl (x));
- else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
- }
+ return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
else if (SCM_REALP (x))
{
double xx = SCM_REAL_VALUE (x);
@@ -8014,11 +8069,7 @@ do_divide (SCM x, SCM y, int inexact)
#endif
}
else if (xx % yy != 0)
- {
- if (inexact)
- return scm_from_double ((double) xx / (double) yy);
- else return scm_i_make_ratio (x, y);
- }
+ return scm_i_make_ratio (x, y);
else
{
scm_t_inum z = xx / yy;
@@ -8029,11 +8080,7 @@ do_divide (SCM x, SCM y, int inexact)
}
}
else if (SCM_BIGP (y))
- {
- if (inexact)
- return scm_from_double ((double) xx / scm_i_big2dbl (y));
- else return scm_i_make_ratio (x, y);
- }
+ return scm_i_make_ratio (x, y);
else if (SCM_REALP (y))
{
double yy = SCM_REAL_VALUE (y);
@@ -8090,29 +8137,10 @@ do_divide (SCM x, SCM y, int inexact)
else if (yy == 1)
return x;
else
- {
- if (inexact)
- return scm_from_double (scm_i_big2dbl (x) / (double) yy);
- else return scm_i_make_ratio (x, y);
- }
- }
- else if (SCM_BIGP (y))
- {
- /* big_x / big_y */
- if (inexact)
- {
- /* It's easily possible for the ratio x/y to fit a double
- but one or both x and y be too big to fit a double,
- hence the use of mpq_get_d rather than converting and
- dividing. */
- mpq_t q;
- *mpq_numref(q) = *SCM_I_BIG_MPZ (x);
- *mpq_denref(q) = *SCM_I_BIG_MPZ (y);
- return scm_from_double (mpq_get_d (q));
- }
- else
return scm_i_make_ratio (x, y);
}
+ else if (SCM_BIGP (y))
+ return scm_i_make_ratio (x, y);
else if (SCM_REALP (y))
{
double yy = SCM_REAL_VALUE (y);
@@ -8273,17 +8301,6 @@ do_divide (SCM x, SCM y, int inexact)
else
SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide);
}
-
-SCM
-scm_divide (SCM x, SCM y)
-{
- return do_divide (x, y, 0);
-}
-
-static SCM scm_divide2real (SCM x, SCM y)
-{
- return do_divide (x, y, 1);
-}
#undef FUNC_NAME
@@ -9534,12 +9551,11 @@ log_of_fraction (SCM n, SCM d)
log_of_exact_integer (d)));
else if (scm_is_false (scm_negative_p (n)))
return scm_from_double
- (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
+ (log1p (scm_i_divide2double (scm_difference (n, d), d)));
else
return scm_c_make_rectangular
- (log1p (scm_to_double (scm_divide2real
- (scm_difference (scm_abs (n), d),
- d))),
+ (log1p (scm_i_divide2double (scm_difference (scm_abs (n), d),
+ d)),
M_PI);
}
--
1.7.2.5
- [PATCH] Improvements to exact rationals et al,
Mark H Weaver <=