* [PATCHES] Numeric improvements
@ 2013-03-05 11:04 Mark H Weaver
2013-03-06 14:05 ` Ludovic Courtès
2013-03-17 23:49 ` Mark H Weaver
0 siblings, 2 replies; 8+ messages in thread
From: Mark H Weaver @ 2013-03-05 11:04 UTC (permalink / raw)
To: guile-devel
[-- Attachment #1: Type: text/plain, Size: 569 bytes --]
Hello all,
Here are ten patches to improve numerics. Among other things, this
eliminates the known obstacles to linking with mini-gmp
<http://bugs.gnu.org/10519>, fixes several problems with our number
printer <http://bugs.gnu.org/13757>, adds 'round-ash', and speeds up
handling of exact rationals.
Still to be done: add more tests to ensure full coverage of all code
paths. It's possible that there are some bugs lurking. However, I
wanted to post it now to allow review and possible work on mini-gmp
integration.
Comments and suggestions welcome.
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: [PATCH 01/10] Improve code in scm_gcd for inum/inum case --]
[-- Type: text/x-diff, Size: 2480 bytes --]
From cd9784ed33d78e6647a752123bf7be91d65b5c96 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:34:17 -0500
Subject: [PATCH 01/10] Improve code in scm_gcd for inum/inum case
* libguile/numbers.c (scm_gcd): Improve implementation of inum/inum case
to be more clear and efficient.
---
libguile/numbers.c | 54 +++++++++++++++++++++++++++++-----------------------
1 file changed, 30 insertions(+), 24 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 66c95db..9c28a79 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -3889,52 +3889,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: [PATCH 02/10] Optimize and simplify fractions code --]
[-- Type: text/x-diff, Size: 13959 bytes --]
From f6201616f7304979a31ab41814e2b297f74a3484 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:34:50 -0500
Subject: [PATCH 02/10] 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, then a^n and b^n have no
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.
* test-suite/tests/numbers.test (expt, integer-expt): Add tests.
---
libguile/numbers.c | 242 +++++++++++++++++++++++++----------------
test-suite/tests/numbers.test | 6 +
2 files changed, 157 insertions(+), 91 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 9c28a79..b5ba1a0 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -439,96 +439,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_t_inum yy = SCM_I_INUM (denominator);
- if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy))
- return scm_divide (numerator, denominator);
- }
- else
+ SCM the_gcd = scm_gcd (numerator, denominator);
+ if (!(scm_is_eq (the_gcd, SCM_INUM1)))
{
- 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
@@ -820,8 +782,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);
@@ -889,6 +852,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
@@ -4672,6 +4712,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);
@@ -7327,8 +7384,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);
}
@@ -7876,14 +7934,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))
{
@@ -7913,8 +7971,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);
}
@@ -8877,8 +8935,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);
@@ -8979,8 +9038,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 66aa01a..20b7eb1 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4044,6 +4044,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)))
@@ -4317,6 +4320,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #4: [PATCH 03/10] Add 'round-ash', a rounding arithmetic shift operator --]
[-- Type: text/x-diff, Size: 18198 bytes --]
From 09b12cb11a64a3cecc15fd54a8e3bb20ead0845c Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:35:09 -0500
Subject: [PATCH 03/10] 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'.
---
doc/ref/api-data.texi | 44 +++++---
libguile/numbers.c | 228 ++++++++++++++++++++++++++++-------------
libguile/numbers.h | 3 +-
test-suite/tests/numbers.test | 114 +++++++++------------
4 files changed, 237 insertions(+), 152 deletions(-)
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index 9da17d8..b33270c 100644
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1686,19 +1686,16 @@ starts from 0 for the least significant bit.
@end lisp
@end deffn
-@deffn {Scheme Procedure} ash n cnt
-@deffnx {C Function} scm_ash (n, cnt)
-Return @var{n} shifted left by @var{cnt} bits, or shifted right if
-@var{cnt} is negative. This is an ``arithmetic'' shift.
-
-This is effectively a multiplication by @m{2^{cnt}, 2^@var{cnt}}, 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.)
+@deffn {Scheme Procedure} ash n count
+@deffnx {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * 2^@var{count})}, but computed
+more efficiently. @var{n} and @var{count} must be exact
+integers.
-With @var{n} viewed as an infinite precision twos complement,
-@code{ash} 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"
@@ -1709,6 +1706,29 @@ dropping bits.
@end lisp
@end deffn
+@deffn {Scheme Procedure} round-ash n count
+@deffnx {C Function} scm_round_ash (n, count)
+Return @math{round(@var{n} * 2^@var{count})}, 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.
+
+@lisp
+(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\"
+@end lisp
+@end 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 b5ba1a0..aa63fbd 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4786,19 +4786,120 @@ 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))
+ {
+ if (count >= SCM_I_FIXNUM_BIT)
+ return SCM_INUM0;
+ else
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+ scm_t_inum qq = SCM_SRS (nn, count);
+
+ 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)
+ && (mpz_odd_p (SCM_I_BIG_MPZ (q))
+ || (mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1)))
+ 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"
+ (SCM n, SCM count),
+ "Return @math{floor(@var{n} * 2^@var{count})}, but computed\n"
+ "more efficiently. @var{n} and @var{count} must be exact\n"
+ "integers.\n"
"\n"
- "This is effectively a multiplication by 2^@var{cnt}, 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"
- "\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"
@@ -4809,79 +4910,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} * 2^@var{count})}, 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 2c8b260..912f287 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -206,7 +206,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 20b7eb1..8dab29c 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?
;;;
@@ -4904,3 +4839,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #5: [PATCH 04/10] Verify that FLT_RADIX is 2. --]
[-- Type: text/x-diff, Size: 777 bytes --]
From e163f7fbf1b3af95fef5ca13f8a36166cbdb0233 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 4 Mar 2013 18:37:23 -0500
Subject: [PATCH 04/10] Verify that FLT_RADIX is 2.
* libguile/numbers.c: Trigger a compilation error if FLT_RADIX is not 2.
This has long been assumed by code in numbers.c.
---
libguile/numbers.c | 3 +++
1 file changed, 3 insertions(+)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index aa63fbd..f9c65a8 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -81,6 +81,9 @@
#define M_PI 3.14159265358979323846
#endif
+/* FIXME: We assume that FLT_RADIX is 2 */
+verify (FLT_RADIX == 2);
+
typedef scm_t_signed_bits scm_t_inum;
#define scm_from_inum(x) (scm_from_signed_integer (x))
--
1.7.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #6: [PATCH 05/10] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp --]
[-- Type: text/x-diff, Size: 4755 bytes --]
From ccd155332040dff20dff0dde63fd2a7656f0900a Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:58:55 -0500
Subject: [PATCH 05/10] Simplify and improve scm_i_big2dbl, and add
scm_i_big2dbl_2exp
* libguile/numbers.c: Verify that FLT_RADIX == 2, which has long
been assumed in various places.
(scm_i_big2dbl_2exp): New static function.
(scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp'.
Ties now go to even.
---
libguile/numbers.c | 101 +++++++++++++++++++---------------------------------
1 file changed, 36 insertions(+), 65 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index f9c65a8..11911d7 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -330,81 +330,52 @@ 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.
+static SCM round_right_shift_exact_integer (SCM n, long count);
- 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.
+/* 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. */
- 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 DBL_MANT_DIG 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)
+static double
+scm_i_big2dbl_2exp (SCM b, long *expon_p)
{
- 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 */
-
- /* 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)
- {
- size_t shift = bits - DBL_MANT_DIG;
- mpz_init2 (tmp, DBL_MANT_DIG);
- 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
+ size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
+ size_t shift = 0;
if (bits > DBL_MANT_DIG)
{
- unsigned long pos = bits - DBL_MANT_DIG - 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)))
+ shift = bits - DBL_MANT_DIG;
+ b = round_right_shift_exact_integer (b, shift);
+ if (SCM_I_INUMP (b))
{
- result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
+ int expon;
+ double signif = frexp (SCM_I_INUM (b), &expon);
+ *expon_p = expon + shift;
+ return signif;
}
}
- scm_remember_upto_here_1 (b);
- return result;
+ {
+ long expon;
+ double signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b));
+ scm_remember_upto_here_1 (b);
+ *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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #7: [PATCH 06/10] Optimize logarithms using scm_i_big2dbl_2exp --]
[-- Type: text/x-diff, Size: 2344 bytes --]
From 5fc88dd6f30a693eff5bea2faadb69c12d3a8359 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 05:02:53 -0500
Subject: [PATCH 06/10] 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 file changed, 12 insertions(+), 18 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 11911d7..e9059a3 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9576,26 +9576,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 * 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 */
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 */
@@ -9606,8 +9600,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #8: [PATCH 07/10] Reimplement 'inexact->exact' to avoid mpq functions. --]
[-- Type: text/x-diff, Size: 2165 bytes --]
From fcd504b656a47ef0b92d916c2503214ed0220002 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 4 Mar 2013 18:42:27 -0500
Subject: [PATCH 07/10] Reimplement 'inexact->exact' to avoid mpq functions.
* libguile/numbers.c (scm_inexact_to_exact): Implement conversion of a
double to an exact rational without using the mpq functions.
---
libguile/numbers.c | 41 +++++++++++++++++++++++++++--------------
1 file changed, 27 insertions(+), 14 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index e9059a3..5feed70 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9085,22 +9085,35 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
if (!SCM_LIKELY (DOUBLE_IS_FINITE (val)))
SCM_OUT_OF_RANGE (1, z);
+ else if (val == 0.0)
+ return SCM_INUM0;
else
{
- mpq_t frac;
- SCM q;
-
- mpq_init (frac);
- mpq_set_d (frac, val);
- 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...
- */
- mpq_clear (frac);
- return q;
+ int expon;
+ SCM numerator;
+
+ numerator = scm_i_dbl2big (ldexp (frexp (val, &expon),
+ DBL_MANT_DIG));
+ expon -= DBL_MANT_DIG;
+ if (expon < 0)
+ {
+ int shift = mpz_scan1 (SCM_I_BIG_MPZ (numerator), 0);
+
+ if (shift > -expon)
+ shift = -expon;
+ mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (numerator),
+ SCM_I_BIG_MPZ (numerator),
+ shift);
+ expon += shift;
+ }
+ numerator = scm_i_normbig (numerator);
+ if (expon < 0)
+ return scm_i_make_ratio_already_reduced
+ (numerator, left_shift_exact_integer (SCM_INUM1, -expon));
+ else if (expon > 0)
+ return left_shift_exact_integer (numerator, expon);
+ else
+ return numerator;
}
}
}
--
1.7.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #9: [PATCH 08/10] Improve inexact division of exact integers. --]
[-- Type: text/x-diff, Size: 13863 bytes --]
From 3f5f479b73744d939d15aa1a5736e013ded3ea1d Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 4 Mar 2013 18:46:33 -0500
Subject: [PATCH 08/10] Improve inexact division of exact integers.
* libguile/numbers.c (scm_i_divide2double): New function.
(scm_i_divide2double_lo2b): New variable.
(scm_i_fraction2double, log_of_fraction): Use 'scm_i_divide2double'.
(do_divide): Removed. Its code is now in 'scm_divide'.
(scm_divide2real): Removed. Superceded by 'scm_i_divide2double'.
(scm_divide): Inherit code from 'do_divide', but without support for
forcing a 'double' result (that functionality is now implemented by
'scm_i_divide2double'). Add FIXME comments in cases where divisions
might not be as precise as they should be.
(scm_init_numbers): Initialize 'scm_i_divide2double_lo2b'.
---
libguile/numbers.c | 257 ++++++++++++++++++++++++++++++++++++----------------
1 file changed, 177 insertions(+), 80 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 5feed70..601b2a6 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -410,9 +410,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:
@@ -468,11 +465,122 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
}
#undef FUNC_NAME
+static mpz_t scm_i_divide2double_lo2b;
+
+/* Return n/d as a double, where n and d are exact integers. */
+static double
+scm_i_divide2double (SCM n, SCM d)
+{
+ int neg;
+ mpz_t nn, dd, lo, hi, x;
+ ssize_t e;
+
+ if (SCM_I_INUMP (d))
+ {
+ if (SCM_UNLIKELY (scm_is_eq (d, SCM_INUM0)))
+ {
+ if (scm_is_true (scm_positive_p (n)))
+ return 1.0 / 0.0;
+ else if (scm_is_true (scm_negative_p (n)))
+ return -1.0 / 0.0;
+ else
+ return 0.0 / 0.0;
+ }
+ mpz_init_set_si (dd, SCM_I_INUM (d));
+ }
+ else
+ mpz_init_set (dd, SCM_I_BIG_MPZ (d));
+
+ if (SCM_I_INUMP (n))
+ mpz_init_set_si (nn, SCM_I_INUM (n));
+ else
+ mpz_init_set (nn, SCM_I_BIG_MPZ (n));
+
+ neg = (mpz_sgn (nn) < 0) ^ (mpz_sgn (dd) < 0);
+ mpz_abs (nn, nn);
+ mpz_abs (dd, dd);
+
+ /* Now we need to find the value of e such that:
+
+ For e <= 0:
+ b^{p-1} - 1/2b <= b^-e n / d < b^p - 1/2
+ (2 b^p - 1) <= 2 b b^-e n / d < (2 b^p - 1) b
+ (2 b^p - 1) d <= 2 b b^-e n < (2 b^p - 1) d b
+
+ For e >= 0:
+ b^{p-1} - 1/2b <= n / b^e d < b^p - 1/2
+ (2 b^p - 1) <= 2 b n / b^e d < (2 b^p - 1) b
+ (2 b^p - 1) d b^e <= 2 b n < (2 b^p - 1) d b b^e
+
+ p == DBL_MANT_DIG
+ b == FLT_RADIX (assumed 2 here) */
+
+ /* Initial guess for e */
+ e = mpz_sizeinbase (nn, 2) - mpz_sizeinbase (dd, 2) - (DBL_MANT_DIG-1);
+ if (e < DBL_MIN_EXP - DBL_MANT_DIG)
+ e = DBL_MIN_EXP - DBL_MANT_DIG;
+
+ /* Compute initial lo, hi, and x */
+ mpz_inits (lo, hi, x, NULL);
+ mpz_mul_2exp (x, nn, 2 + ((e < 0) ? -e : 0));
+ mpz_mul (lo, dd, scm_i_divide2double_lo2b);
+ if (e > 0)
+ mpz_mul_2exp (lo, lo, e);
+ mpz_mul_2exp (hi, lo, 1);
+
+ /* Adjust e as needed */
+ while (mpz_cmp (x, lo) < 0 && e > DBL_MIN_EXP - DBL_MANT_DIG)
+ {
+ mpz_mul_2exp (x, x, 1);
+ e--;
+ }
+ while (mpz_cmp (x, hi) >= 0)
+ {
+ mpz_mul_2exp (hi, hi, 1);
+ e++;
+ }
+
+ /* Now compute the rounded mantissa:
+ n / b^e d (if e >= 0)
+ n b^-e / d (if e <= 0) */
+ {
+ int cmp, needs_adjustment;
+ double result;
+
+ if (e < 0)
+ mpz_mul_2exp (nn, nn, -e);
+ else
+ mpz_mul_2exp (dd, dd, e);
+
+ /* mpz does not directly support rounded right
+ shifts, so we have to do it the hard way.
+ hi == quotient, lo == remainder */
+ mpz_fdiv_qr (hi, lo, nn, dd);
+ mpz_mul_2exp (lo, lo, 1);
+
+ cmp = mpz_cmp (lo, dd);
+ if (mpz_odd_p (hi))
+ needs_adjustment = (cmp >= 0);
+ else
+ needs_adjustment = (cmp > 0);
+
+ if (needs_adjustment)
+ mpz_add_ui (hi, hi, 1);
+
+ if (neg)
+ mpz_neg (hi, hi);
+
+ result = ldexp (mpz_get_d (hi), e);
+ mpz_clears (nn, dd, lo, hi, x, NULL);
+ return result;
+ }
+}
+
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
@@ -7965,8 +8073,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;
@@ -7985,18 +8093,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);
@@ -8046,11 +8146,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;
@@ -8061,11 +8157,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);
@@ -8074,6 +8166,9 @@ do_divide (SCM x, SCM y, int inexact)
scm_num_overflow (s_divide);
else
#endif
+ /* FIXME: Precision may be lost here due to:
+ (1) The cast from 'scm_t_inum' to 'double'
+ (2) Double rounding */
return scm_from_double ((double) xx / yy);
}
else if (SCM_COMPLEXP (y))
@@ -8100,7 +8195,7 @@ do_divide (SCM x, SCM y, int inexact)
else if (SCM_FRACTIONP (y))
/* a / b/c = ac / b */
return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
- SCM_FRACTION_NUMERATOR (y));
+ SCM_FRACTION_NUMERATOR (y));
else
SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
}
@@ -8144,43 +8239,24 @@ do_divide (SCM x, SCM y, int inexact)
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);
- }
+ 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
- {
- 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);
- }
+ 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);
}
else if (SCM_REALP (y))
{
@@ -8190,6 +8266,8 @@ do_divide (SCM x, SCM y, int inexact)
scm_num_overflow (s_divide);
else
#endif
+ /* FIXME: Precision may be lost here due to:
+ (1) scm_i_big2dbl (2) Double rounding */
return scm_from_double (scm_i_big2dbl (x) / yy);
}
else if (SCM_COMPLEXP (y))
@@ -8199,7 +8277,7 @@ do_divide (SCM x, SCM y, int inexact)
}
else if (SCM_FRACTIONP (y))
return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
- SCM_FRACTION_NUMERATOR (y));
+ SCM_FRACTION_NUMERATOR (y));
else
SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
}
@@ -8214,10 +8292,16 @@ do_divide (SCM x, SCM y, int inexact)
scm_num_overflow (s_divide);
else
#endif
+ /* FIXME: Precision may be lost here due to:
+ (1) The cast from 'scm_t_inum' to 'double'
+ (2) Double rounding */
return scm_from_double (rx / (double) yy);
}
else if (SCM_BIGP (y))
{
+ /* FIXME: Precision may be lost here due to:
+ (1) The conversion from bignum to double
+ (2) Double rounding */
double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
scm_remember_upto_here_1 (y);
return scm_from_double (rx / dby);
@@ -8255,12 +8339,18 @@ do_divide (SCM x, SCM y, int inexact)
else
#endif
{
+ /* FIXME: Precision may be lost here due to:
+ (1) The conversion from 'scm_t_inum' to double
+ (2) Double rounding */
double d = yy;
return scm_c_make_rectangular (rx / d, ix / d);
}
}
else if (SCM_BIGP (y))
{
+ /* FIXME: Precision may be lost here due to:
+ (1) The conversion from bignum to double
+ (2) Double rounding */
double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
scm_remember_upto_here_1 (y);
return scm_c_make_rectangular (rx / dby, ix / dby);
@@ -8294,6 +8384,9 @@ do_divide (SCM x, SCM y, int inexact)
}
else if (SCM_FRACTIONP (y))
{
+ /* FIXME: Precision may be lost here due to:
+ (1) The conversion from fraction to double
+ (2) Double rounding */
double yy = scm_i_fraction2double (y);
return scm_c_make_rectangular (rx / yy, ix / yy);
}
@@ -8311,12 +8404,12 @@ do_divide (SCM x, SCM y, int inexact)
else
#endif
return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
- scm_product (SCM_FRACTION_DENOMINATOR (x), y));
+ scm_product (SCM_FRACTION_DENOMINATOR (x), y));
}
else if (SCM_BIGP (y))
{
return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
- scm_product (SCM_FRACTION_DENOMINATOR (x), y));
+ scm_product (SCM_FRACTION_DENOMINATOR (x), y));
}
else if (SCM_REALP (y))
{
@@ -8326,33 +8419,28 @@ do_divide (SCM x, SCM y, int inexact)
scm_num_overflow (s_divide);
else
#endif
+ /* FIXME: Precision may be lost here due to:
+ (1) The conversion from fraction to double
+ (2) Double rounding */
return scm_from_double (scm_i_fraction2double (x) / yy);
}
else if (SCM_COMPLEXP (y))
{
+ /* FIXME: Precision may be lost here due to:
+ (1) The conversion from fraction to double
+ (2) Double rounding */
a = scm_i_fraction2double (x);
goto complex_div;
}
else if (SCM_FRACTIONP (y))
return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
- scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x)));
+ scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x)));
else
SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
}
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
@@ -9617,12 +9705,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);
}
@@ -9890,6 +9977,16 @@ scm_init_numbers ()
#endif
exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2));
+
+ {
+ /* Set scm_i_divide2double_lo2b to (2 b^p - 1) */
+ mpz_init_set_ui (scm_i_divide2double_lo2b, 1);
+ mpz_mul_2exp (scm_i_divide2double_lo2b,
+ scm_i_divide2double_lo2b,
+ DBL_MANT_DIG + 1); /* 2 b^p */
+ mpz_sub_ui (scm_i_divide2double_lo2b, scm_i_divide2double_lo2b, 1);
+ }
+
#include "libguile/numbers.x"
}
--
1.7.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #10: [PATCH 09/10] Support reading negative exponents larger than SCM_MAXEXP. --]
[-- Type: text/x-diff, Size: 896 bytes --]
From 8d3f8e1b1fcc1abac1967bb9cf9643b63f3dea6d Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Tue, 5 Mar 2013 05:45:37 -0500
Subject: [PATCH 09/10] Support reading negative exponents larger than
SCM_MAXEXP.
* libguile/numbers.c (mem2decimal_from_point): Accept negative exponents
larger than SCM_MAXEXP that produce subnormals.
---
libguile/numbers.c | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 601b2a6..c19b380 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5929,7 +5929,7 @@ mem2decimal_from_point (SCM result, SCM mem,
break;
}
- if (exponent > SCM_MAXEXP)
+ if (exponent > ((sign == 1) ? SCM_MAXEXP : SCM_MAXEXP + DBL_DIG + 1))
{
size_t exp_len = idx - start;
SCM exp_string = scm_i_substring_copy (mem, start, start + exp_len);
--
1.7.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #11: [PATCH 10/10] Reimplement idbl2str number printer. --]
[-- Type: text/x-diff, Size: 8819 bytes --]
From 9fbdffa44db46dfabb940949f73b142a6a033b0d Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Tue, 5 Mar 2013 05:47:56 -0500
Subject: [PATCH 10/10] Reimplement idbl2str number printer.
Fixes <http://bugs.gnu.org/13757>.
* libguile/numbers.c (idbl2str): Reimplement.
---
libguile/numbers.c | 339 +++++++++++++++++++++++++++-------------------------
1 file changed, 179 insertions(+), 160 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index c19b380..6b3cb16 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5267,185 +5267,196 @@ void init_fx_radix(double *fx_list, int radix)
/* use this array as a way to generate a single digit */
static const char number_chars[] = "0123456789abcdefghijklmnopqrstuvwxyz";
+static mpz_t dbl_minimum_normal_mantissa;
+
static size_t
-idbl2str (double f, char *a, int radix)
-{
- int efmt, dpt, d, i, wp;
- double *fx;
-#ifdef DBL_MIN_10_EXP
- double f_cpy;
- int exp_cpy;
-#endif /* DBL_MIN_10_EXP */
- size_t ch = 0;
- int exp = 0;
-
- if(radix < 2 ||
- radix > SCM_MAX_DBL_RADIX)
- {
- /* revert to existing behavior */
- radix = 10;
- }
+idbl2str (double dbl, char *a, int radix)
+{
+ int ch = 0;
- wp = scm_dblprec[radix-2];
- fx = fx_per_radix[radix-2];
+ if (radix < 2 || radix > SCM_MAX_DBL_RADIX)
+ /* revert to existing behavior */
+ radix = 10;
- if (f == 0.0)
+ if (isinf (dbl))
{
-#ifdef HAVE_COPYSIGN
- double sgn = copysign (1.0, f);
-
- if (sgn < 0.0)
- a[ch++] = '-';
-#endif
- goto zero; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
+ strcpy (a, (dbl > 0.0) ? "+inf.0" : "-inf.0");
+ return 6;
}
-
- if (isinf (f))
+ else if (dbl > 0.0)
+ ;
+ else if (dbl < 0.0)
{
- if (f < 0)
- strcpy (a, "-inf.0");
- else
- strcpy (a, "+inf.0");
- return ch+6;
+ dbl = -dbl;
+ a[ch++] = '-';
}
- else if (isnan (f))
+ else if (dbl == 0.0)
{
- strcpy (a, "+nan.0");
- return ch+6;
+ if (!double_is_non_negative_zero (dbl))
+ a[ch++] = '-';
+ strcpy (a + ch, "0.0");
+ return ch + 3;
}
-
- if (f < 0.0)
+ else if (isnan (dbl))
{
- f = -f;
- a[ch++] = '-';
+ strcpy (a, "+nan.0");
+ return 6;
}
-#ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
- make-uniform-vector, from causing infinite loops. */
- /* just do the checking...if it passes, we do the conversion for our
- radix again below */
- f_cpy = f;
- exp_cpy = exp;
+ /* Algorithm taken from "Printing Floating-Point Numbers Quickly and
+ Accurately" by Robert G. Burger and R. Kent Dybvig */
+ {
+ int e, k;
+ mpz_t f, r, s, mplus, mminus, hi, digit;
+ int f_is_even, f_is_odd;
+ int show_exp = 0;
+
+ mpz_inits (f, r, s, mplus, mminus, hi, digit, NULL);
+ mpz_set_d (f, ldexp (frexp (dbl, &e), DBL_MANT_DIG));
+ if (e < DBL_MIN_EXP)
+ {
+ mpz_tdiv_q_2exp (f, f, DBL_MIN_EXP - e);
+ e = DBL_MIN_EXP;
+ }
+ e -= DBL_MANT_DIG;
- while (f_cpy < 1.0)
- {
- f_cpy *= 10.0;
- if (exp_cpy-- < DBL_MIN_10_EXP)
- {
- a[ch++] = '#';
- a[ch++] = '.';
- a[ch++] = '#';
- return ch;
- }
- }
- while (f_cpy > 10.0)
- {
- f_cpy *= 0.10;
- if (exp_cpy++ > DBL_MAX_10_EXP)
- {
- a[ch++] = '#';
- a[ch++] = '.';
- a[ch++] = '#';
- return ch;
- }
- }
-#endif
+ f_is_even = !mpz_odd_p (f);
+ f_is_odd = !f_is_even;
- while (f < 1.0)
- {
- f *= radix;
- exp--;
- }
- while (f > radix)
- {
- f /= radix;
- exp++;
- }
+ /* Initialize r, s, mplus, and mminus according
+ to Table 1 from the paper. */
+ if (e < 0)
+ {
+ mpz_set_ui (mminus, 1);
+ if (mpz_cmp (f, dbl_minimum_normal_mantissa) != 0
+ || e == DBL_MIN_EXP - DBL_MANT_DIG)
+ {
+ mpz_set_ui (mplus, 1);
+ mpz_mul_2exp (r, f, 1);
+ mpz_mul_2exp (s, mminus, 1 - e);
+ }
+ else
+ {
+ mpz_set_ui (mplus, 2);
+ mpz_mul_2exp (r, f, 2);
+ mpz_mul_2exp (s, mminus, 2 - e);
+ }
+ }
+ else
+ {
+ mpz_set_ui (mminus, 1);
+ mpz_mul_2exp (mminus, mminus, e);
+ if (mpz_cmp (f, dbl_minimum_normal_mantissa) != 0)
+ {
+ mpz_set (mplus, mminus);
+ mpz_mul_2exp (r, f, 1 + e);
+ mpz_set_ui (s, 2);
+ }
+ else
+ {
+ mpz_mul_2exp (mplus, mminus, 1);
+ mpz_mul_2exp (r, f, 2 + e);
+ mpz_set_ui (s, 4);
+ }
+ }
- if (f + fx[wp] >= radix)
+ /* Find the smallest k such that:
+ (r + mplus) / s < radix^k (if f is even)
+ (r + mplus) / s <= radix^k (if f is odd) */
{
- f = 1.0;
- exp++;
- }
- zero:
-#ifdef ENGNOT
- /* adding 9999 makes this equivalent to abs(x) % 3 */
- dpt = (exp + 9999) % 3;
- exp -= dpt++;
- efmt = 1;
-#else
- efmt = (exp < -3) || (exp > wp + 2);
- if (!efmt)
- {
- if (exp < 0)
- {
- a[ch++] = '0';
- a[ch++] = '.';
- dpt = exp;
- while (++dpt)
- a[ch++] = '0';
- }
- else
- dpt = exp + 1;
+ /* IMPROVE-ME: Make an initial guess to speed this up */
+ mpz_add (hi, r, mplus);
+ k = 0;
+ while (mpz_cmp (hi, s) >= f_is_odd)
+ {
+ mpz_mul_ui (s, s, radix);
+ k++;
+ }
+ if (k == 0)
+ {
+ mpz_mul_ui (hi, hi, radix);
+ while (mpz_cmp (hi, s) < f_is_odd)
+ {
+ mpz_mul_ui (r, r, radix);
+ mpz_mul_ui (mplus, mplus, radix);
+ mpz_mul_ui (mminus, mminus, radix);
+ mpz_mul_ui (hi, hi, radix);
+ k--;
+ }
+ }
}
- else
- dpt = 1;
-#endif
- do
- {
- d = f;
- f -= d;
- a[ch++] = number_chars[d];
- if (f < fx[wp])
- break;
- if (f + fx[wp] >= 1.0)
- {
- a[ch - 1] = number_chars[d+1];
- break;
- }
- f *= radix;
- if (!(--dpt))
- a[ch++] = '.';
- }
- while (wp--);
+ if (k >= 8 || k <= -3)
+ {
+ /* Use scientific notation */
+ show_exp = k - 1;
+ k = 1;
+ }
+ else if (k <= 0)
+ {
+ int i;
- if (dpt > 0)
- {
-#ifndef ENGNOT
- if ((dpt > 4) && (exp > 6))
- {
- d = (a[0] == '-' ? 2 : 1);
- for (i = ch++; i > d; i--)
- a[i] = a[i - 1];
- a[d] = '.';
- efmt = 1;
- }
- else
-#endif
- {
- while (--dpt)
- a[ch++] = '0';
- a[ch++] = '.';
- }
- }
- if (a[ch - 1] == '.')
- a[ch++] = '0'; /* trailing zero */
- if (efmt && exp)
- {
- a[ch++] = 'e';
- if (exp < 0)
- {
- exp = -exp;
- a[ch++] = '-';
- }
- for (i = radix; i <= exp; i *= radix);
- for (i /= radix; i; i /= radix)
- {
- a[ch++] = number_chars[exp / i];
- exp %= i;
- }
- }
+ /* Print leading zeroes */
+ a[ch++] = '0';
+ a[ch++] = '.';
+ for (i = 0; i > k; i--)
+ a[ch++] = '0';
+ }
+
+ for (;;)
+ {
+ int end_1_p, end_2_p;
+ int d;
+
+ mpz_mul_ui (mplus, mplus, radix);
+ mpz_mul_ui (mminus, mminus, radix);
+ mpz_mul_ui (r, r, radix);
+ mpz_fdiv_qr (digit, r, r, s);
+ d = mpz_get_ui (digit);
+
+ mpz_add (hi, r, mplus);
+ end_1_p = (mpz_cmp (r, mminus) < f_is_even);
+ end_2_p = (mpz_cmp (s, hi) < f_is_even);
+ if (end_1_p || end_2_p)
+ {
+ mpz_mul_2exp (r, r, 1);
+ if (!end_2_p)
+ ;
+ else if (!end_1_p)
+ d++;
+ else if (mpz_cmp (r, s) >= !(d & 1))
+ d++;
+ a[ch++] = number_chars[d];
+ if (--k == 0)
+ a[ch++] = '.';
+ break;
+ }
+ else
+ {
+ a[ch++] = number_chars[d];
+ if (--k == 0)
+ a[ch++] = '.';
+ }
+ }
+
+ if (k > 0)
+ {
+ for (; k > 0; k--)
+ a[ch++] = '0';
+ a[ch++] = '.';
+ }
+
+ if (k == 0)
+ a[ch++] = '0';
+
+ if (show_exp)
+ {
+ a[ch++] = 'e';
+ ch += scm_iint2str (show_exp, radix, a + ch);
+ }
+
+ mpz_clears (f, r, s, mplus, mminus, hi, digit, NULL);
+ }
return ch;
}
@@ -9987,6 +9998,14 @@ scm_init_numbers ()
mpz_sub_ui (scm_i_divide2double_lo2b, scm_i_divide2double_lo2b, 1);
}
+ {
+ /* Set dbl_minimum_normal_mantissa to b^{p-1} */
+ mpz_init_set_ui (dbl_minimum_normal_mantissa, 1);
+ mpz_mul_2exp (dbl_minimum_normal_mantissa,
+ dbl_minimum_normal_mantissa,
+ DBL_MANT_DIG - 1);
+ }
+
#include "libguile/numbers.x"
}
--
1.7.10.4
^ permalink raw reply related [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-05 11:04 [PATCHES] Numeric improvements Mark H Weaver
@ 2013-03-06 14:05 ` Ludovic Courtès
2013-03-06 20:30 ` Mark H Weaver
2013-03-07 0:16 ` Mark H Weaver
2013-03-17 23:49 ` Mark H Weaver
1 sibling, 2 replies; 8+ messages in thread
From: Ludovic Courtès @ 2013-03-06 14:05 UTC (permalink / raw)
To: guile-devel
Hi,
Mark H Weaver <mhw@netris.org> skribis:
> Here are ten patches to improve numerics. Among other things, this
> eliminates the known obstacles to linking with mini-gmp
> <http://bugs.gnu.org/10519>, fixes several problems with our number
> printer <http://bugs.gnu.org/13757>, adds 'round-ash', and speeds up
> handling of exact rationals.
Ten patches!
> Still to be done: add more tests to ensure full coverage of all code
> paths. It's possible that there are some bugs lurking. However, I
> wanted to post it now to allow review and possible work on mini-gmp
> integration.
Which of these patches are needed for mini-gmp integration? It would
probably be easier to discuss things separately, by small chunks.
Overall I think I’m mostly incompetent on the numeric stuff, so I’d
mostly comment on form.
At first sight, it seems that there are few tests added, in particular
for the number printer. I think tests must be added along with the
patches that claim to fix something.
> From cd9784ed33d78e6647a752123bf7be91d65b5c96 Mon Sep 17 00:00:00 2001
> From: Mark H Weaver <mhw@netris.org>
> Date: Sun, 3 Mar 2013 04:34:17 -0500
> Subject: [PATCH 01/10] Improve code in scm_gcd for inum/inum case
>
> * libguile/numbers.c (scm_gcd): Improve implementation of inum/inum case
> to be more clear and efficient.
This one looks OK, and should be covered by the tests, according to
<http://hydra.nixos.org/build/4268423/download/2/coverage/libguile/numbers.c.gcov.html>.
Did you measure the performance difference?
> From f6201616f7304979a31ab41814e2b297f74a3484 Mon Sep 17 00:00:00 2001
> From: Mark H Weaver <mhw@netris.org>
> Date: Sun, 3 Mar 2013 04:34:50 -0500
> Subject: [PATCH 02/10] 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, then a^n and b^n have no
> 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.
Can you fold the explanations as comments? In particular those that
explain the implementation strategy, and why it’s more efficient this way.
> +/* 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)
Please use imperative-tense sentences above functions to describe what
they do, without repeating the function name; also make sure to
introduce the arguments in the description, written in capital letters
(info "(standards) Comments").
Some helper functions introduced by the patches lack a comment.
> +@deffn {Scheme Procedure} round-ash n count
> +@deffnx {C Function} scm_round_ash (n, count)
> +Return @math{round(@var{n} * 2^@var{count})}, 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.
> +
> +@lisp
> +(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\"
> +@end lisp
> +@end deffn
What was the rationale for ‘round-ash’?
It seems to address to do two things differently from ‘ash’: it’s more
efficient, and it rounds when right-shifting. The “more efficient” bit
is an implementation detail that should also benefit to ‘ash’ (as a user
I don’t want to have to choose between efficient and non-rounding.)
WDYT?
OK, these are just random comments for now.
Thanks for making Guile benefit from your expertise in this domain!
Ludo’.
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-06 14:05 ` Ludovic Courtès
@ 2013-03-06 20:30 ` Mark H Weaver
2013-03-07 19:02 ` Ludovic Courtès
2013-03-07 0:16 ` Mark H Weaver
1 sibling, 1 reply; 8+ messages in thread
From: Mark H Weaver @ 2013-03-06 20:30 UTC (permalink / raw)
To: Ludovic Courtès; +Cc: guile-devel
Hi Ludovic,
ludo@gnu.org (Ludovic Courtès) writes:
> Which of these patches are needed for mini-gmp integration? It would
> probably be easier to discuss things separately, by small chunks.
Mini-gmp integration depends directly on patches 5 and 7, and those
depend on patches 2, 3, and 4.
> Overall I think I’m mostly incompetent on the numeric stuff, so I’d
> mostly comment on form.
>
> At first sight, it seems that there are few tests added, in particular
> for the number printer. I think tests must be added along with the
> patches that claim to fix something.
Agreed. I'll work on that.
>> From cd9784ed33d78e6647a752123bf7be91d65b5c96 Mon Sep 17 00:00:00 2001
>> From: Mark H Weaver <mhw@netris.org>
>> Date: Sun, 3 Mar 2013 04:34:17 -0500
>> Subject: [PATCH 01/10] Improve code in scm_gcd for inum/inum case
>>
>> * libguile/numbers.c (scm_gcd): Improve implementation of inum/inum case
>> to be more clear and efficient.
>
> This one looks OK, and should be covered by the tests, according to
> <http://hydra.nixos.org/build/4268423/download/2/coverage/libguile/numbers.c.gcov.html>.
>
> Did you measure the performance difference?
Yes. On my x86_64 machine it improves gcd performance on fixnums by
about 4.25%. I went ahead and pushed this one.
> Can you fold the explanations as comments?
[...]
> Please use imperative-tense sentences above functions to describe what
> they do, without repeating the function name; also make sure to
> introduce the arguments in the description, written in capital letters
> (info "(standards) Comments").
>
> Some helper functions introduced by the patches lack a comment.
Will do.
>> +@deffn {Scheme Procedure} round-ash n count
>> +@deffnx {C Function} scm_round_ash (n, count)
>> +Return @math{round(@var{n} * 2^@var{count})}, 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.
>> +
>> +@lisp
>> +(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\"
>> +@end lisp
>> +@end deffn
>
> What was the rationale for ‘round-ash’?
Well, we need it internally to efficiently convert large integers to
floating-point with proper rounding (see the call to
'round_right_shift_exact_integer' in patch 5). Given that, it seemed
reasonable to export it to the user, since it's not possible to do that
job efficiently with our current primitives. However, I don't feel
strongly about it.
> It seems to address to do two things differently from ‘ash’: it’s more
> efficient, and it rounds when right-shifting. The “more efficient” bit
> is an implementation detail that should also benefit to ‘ash’ (as a user
> I don’t want to have to choose between efficient and non-rounding.)
No, 'round-ash' is not more efficient than 'ash'. It's more efficient
than the equivalent (round (* n (expt 2 count))).
Thanks for looking through the patches. I'll work on improving them.
Mark
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-06 14:05 ` Ludovic Courtès
2013-03-06 20:30 ` Mark H Weaver
@ 2013-03-07 0:16 ` Mark H Weaver
2013-03-07 6:03 ` Mark H Weaver
1 sibling, 1 reply; 8+ messages in thread
From: Mark H Weaver @ 2013-03-07 0:16 UTC (permalink / raw)
To: Ludovic Courtès; +Cc: guile-devel
[-- Attachment #1: Type: text/plain, Size: 214 bytes --]
I pushed two of the patches: the 'gcd' patch and the one that verifies
that FLT_RADIX is 2. Of the remaining 8 patches, I've attached improved
versions of the first 4. Reviews welcome.
Regards,
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: [PATCH 1/4] Optimize and simplify fractions code --]
[-- Type: text/x-diff, Size: 12913 bytes --]
From f32e8c5ffd789a6dbee48be74f5bbf32978382c3 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:34:50 -0500
Subject: [PATCH 1/4] Optimize and simplify fractions code.
* libguile/numbers.c (scm_exact_integer_quotient,
scm_i_make_ratio_already_reduced): New static functions.
(scm_i_make_ratio): Rewrite in terms of
'scm_i_make_ratio_already_reduced'.
(scm_integer_expt): Optimize fraction case.
(scm_abs, scm_magnitude, scm_difference, do_divide): Use
'scm_i_make_ratio_already_reduced'.
* test-suite/tests/numbers.test (expt, integer-expt): Add tests.
---
libguile/numbers.c | 244 ++++++++++++++++++++++++++---------------
test-suite/tests/numbers.test | 6 +
2 files changed, 159 insertions(+), 91 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 393cf64..e63c17a 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -442,96 +442,56 @@ 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);
+/* Make the ratio NUMERATOR/DENOMINATOR, where:
+ 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);
+
+/* Make the ratio NUMERATOR/DENOMINATOR */
+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
@@ -823,8 +783,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);
@@ -892,6 +853,84 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0,
}
#undef FUNC_NAME
+/* Return 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). For large integers, this can be computed more
+ efficiently than when the remainder is unknown. */
+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
@@ -4675,6 +4714,26 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
else /* return NaN for (0 ^ k) for negative k per R6RS */
return scm_nan ();
}
+ else if (SCM_FRACTIONP (n))
+ {
+ /* Optimize the fraction case by (a/b)^k ==> (a^k)/(b^k), to avoid
+ needless reduction of intermediate products to lowest terms.
+ If a and b have no common factors, then a^k and b^k have no
+ 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. */
+ 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);
@@ -7330,8 +7389,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);
}
@@ -7879,14 +7939,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))
{
@@ -7916,8 +7976,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);
}
@@ -8880,8 +8940,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);
@@ -8982,8 +9043,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 66aa01a..20b7eb1 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4044,6 +4044,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)))
@@ -4317,6 +4320,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: [PATCH 2/4] Add 'round-ash', a rounding arithmetic shift operator --]
[-- Type: text/x-diff, Size: 17644 bytes --]
From da13ef66508cd76b1c173ca355856b1e1b7667f0 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:35:09 -0500
Subject: [PATCH 2/4] 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
static functions.
(scm_round_ash): New procedure.
(scm_ash): Reimplement in terms of 'left_shift_exact_integer' and
'floor_right_shift_exact_integer'.
* libguile/numbers.h: Add prototype for scm_round_ash. Rename the
second argument of 'scm_ash' 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'.
---
doc/ref/api-data.texi | 44 +++++---
libguile/numbers.c | 226 +++++++++++++++++++++++++++--------------
libguile/numbers.h | 3 +-
test-suite/tests/numbers.test | 114 +++++++++------------
4 files changed, 235 insertions(+), 152 deletions(-)
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index 9da17d8..b33270c 100644
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1686,19 +1686,16 @@ starts from 0 for the least significant bit.
@end lisp
@end deffn
-@deffn {Scheme Procedure} ash n cnt
-@deffnx {C Function} scm_ash (n, cnt)
-Return @var{n} shifted left by @var{cnt} bits, or shifted right if
-@var{cnt} is negative. This is an ``arithmetic'' shift.
-
-This is effectively a multiplication by @m{2^{cnt}, 2^@var{cnt}}, 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.)
+@deffn {Scheme Procedure} ash n count
+@deffnx {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * 2^@var{count})}, but computed
+more efficiently. @var{n} and @var{count} must be exact
+integers.
-With @var{n} viewed as an infinite precision twos complement,
-@code{ash} 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"
@@ -1709,6 +1706,29 @@ dropping bits.
@end lisp
@end deffn
+@deffn {Scheme Procedure} round-ash n count
+@deffnx {C Function} scm_round_ash (n, count)
+Return @math{round(@var{n} * 2^@var{count})}, 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.
+
+@lisp
+(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\"
+@end lisp
+@end 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 e63c17a..b99a04b 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4791,19 +4791,119 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
}
#undef FUNC_NAME
+/* Efficiently compute (N * 2^COUNT),
+ where N is an exact integer, and COUNT > 0. */
+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");
+}
+
+/* Efficiently compute floor (N / 2^COUNT),
+ where N is an exact integer and COUNT > 0. */
+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");
+}
+
+/* Efficiently compute round (N / 2^COUNT),
+ where N is an exact integer and COUNT > 0. */
+static SCM
+round_right_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ if (count >= SCM_I_FIXNUM_BIT)
+ return SCM_INUM0;
+ else
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+ scm_t_inum qq = SCM_SRS (nn, count);
+
+ 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)
+ && (mpz_odd_p (SCM_I_BIG_MPZ (q))
+ || (mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1)))
+ 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"
+ (SCM n, SCM count),
+ "Return @math{floor(@var{n} * 2^@var{count})}.\n"
+ "@var{n} and @var{count} must be exact integers.\n"
"\n"
- "This is effectively a multiplication by 2^@var{cnt}, 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"
- "\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"
@@ -4814,79 +4914,57 @@ 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} * 2^@var{count})}.\n"
+ "@var{n} and @var{count} must be exact 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 2c8b260..912f287 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -206,7 +206,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 20b7eb1..8dab29c 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?
;;;
@@ -4904,3 +4839,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #4: [PATCH 3/4] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp --]
[-- Type: text/x-diff, Size: 7119 bytes --]
From f48e23eba996566adbe112731ceda4761070fbce Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:58:55 -0500
Subject: [PATCH 3/4] Simplify and improve scm_i_big2dbl, and add
scm_i_big2dbl_2exp
* libguile/numbers.c (scm_i_big2dbl_2exp): New static function.
(scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp',
with proper rounding.
* test-suite/tests/numbers.test ("inexact->exact"): Add tests.
---
libguile/numbers.c | 101 +++++++++++++++--------------------------
test-suite/tests/numbers.test | 45 ++++++++++++++++++
2 files changed, 81 insertions(+), 65 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index b99a04b..49b4a50 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -330,81 +330,52 @@ 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.
+static SCM round_right_shift_exact_integer (SCM n, long count);
- 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.
+/* 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. */
- 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 DBL_MANT_DIG 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)
+static double
+scm_i_big2dbl_2exp (SCM b, long *expon_p)
{
- 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 */
-
- /* 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)
- {
- size_t shift = bits - DBL_MANT_DIG;
- mpz_init2 (tmp, DBL_MANT_DIG);
- 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
+ size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
+ size_t shift = 0;
if (bits > DBL_MANT_DIG)
{
- unsigned long pos = bits - DBL_MANT_DIG - 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)))
+ shift = bits - DBL_MANT_DIG;
+ b = round_right_shift_exact_integer (b, shift);
+ if (SCM_I_INUMP (b))
{
- result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
+ int expon;
+ double signif = frexp (SCM_I_INUM (b), &expon);
+ *expon_p = expon + shift;
+ return signif;
}
}
- scm_remember_upto_here_1 (b);
- return result;
+ {
+ long expon;
+ double signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b));
+ scm_remember_upto_here_1 (b);
+ *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
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 8dab29c..e3a1099 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4220,6 +4220,51 @@
(pass-if-exception "nan" exception:out-of-range
(inexact->exact +nan.0))
+
+ (pass-if-equal "round up to odd"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111000101 ->
+ ;; 11111111111111111111111111111111111111111111111111001000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001000)
+ (inexact->exact
+ (exact->inexact
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000101))))
+
+ (pass-if-equal "round down to odd"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111001011 ->
+ ;; 11111111111111111111111111111111111111111111111111001000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001000)
+ (inexact->exact
+ (exact->inexact
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001011))))
+
+ (pass-if-equal "round tie up to even"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111011100 ->
+ ;; 11111111111111111111111111111111111111111111111111100000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b100000)
+ (inexact->exact
+ (exact->inexact
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b011100))))
+
+ (pass-if-equal "round tie down to even"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111000100 ->
+ ;; 11111111111111111111111111111111111111111111111111000000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000000)
+ (inexact->exact
+ (exact->inexact
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000100))))
+
+ (pass-if-equal "round tie up to next power of two"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111111100 ->
+ ;; 100000000000000000000000000000000000000000000000000000000
+ (expt 2 (+ dbl-mant-dig 3))
+ (inexact->exact
+ (exact->inexact
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b111100))))
(with-test-prefix "2.0**i to exact and back"
(do ((i 0 (1+ i))
--
1.7.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #5: [PATCH 4/4] Optimize logarithms using scm_i_big2dbl_2exp --]
[-- Type: text/x-diff, Size: 2354 bytes --]
From b21e3d03049e986b0ae6e0ad59aa79c6c38de842 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 05:02:53 -0500
Subject: [PATCH 4/4] 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 file changed, 12 insertions(+), 18 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 49b4a50..0b13d69 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9576,26 +9576,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 * 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 */
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 */
@@ -9606,8 +9600,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.10.4
^ permalink raw reply related [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-07 0:16 ` Mark H Weaver
@ 2013-03-07 6:03 ` Mark H Weaver
2013-03-12 20:58 ` Mark H Weaver
0 siblings, 1 reply; 8+ messages in thread
From: Mark H Weaver @ 2013-03-07 6:03 UTC (permalink / raw)
To: Ludovic Courtès; +Cc: guile-devel
[-- Attachment #1: Type: text/plain, Size: 148 bytes --]
Here are improved versions of the patches needed to enable mini-gmp
integration. I think these are ready to commit. Reviews welcome.
Mark
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: [PATCH 1/5] Optimize and simplify fractions code --]
[-- Type: text/x-diff, Size: 12913 bytes --]
From f32e8c5ffd789a6dbee48be74f5bbf32978382c3 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:34:50 -0500
Subject: [PATCH 1/5] Optimize and simplify fractions code.
* libguile/numbers.c (scm_exact_integer_quotient,
scm_i_make_ratio_already_reduced): New static functions.
(scm_i_make_ratio): Rewrite in terms of
'scm_i_make_ratio_already_reduced'.
(scm_integer_expt): Optimize fraction case.
(scm_abs, scm_magnitude, scm_difference, do_divide): Use
'scm_i_make_ratio_already_reduced'.
* test-suite/tests/numbers.test (expt, integer-expt): Add tests.
---
libguile/numbers.c | 244 ++++++++++++++++++++++++++---------------
test-suite/tests/numbers.test | 6 +
2 files changed, 159 insertions(+), 91 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 393cf64..e63c17a 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -442,96 +442,56 @@ 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);
+/* Make the ratio NUMERATOR/DENOMINATOR, where:
+ 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);
+
+/* Make the ratio NUMERATOR/DENOMINATOR */
+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
@@ -823,8 +783,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);
@@ -892,6 +853,84 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0,
}
#undef FUNC_NAME
+/* Return 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). For large integers, this can be computed more
+ efficiently than when the remainder is unknown. */
+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
@@ -4675,6 +4714,26 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
else /* return NaN for (0 ^ k) for negative k per R6RS */
return scm_nan ();
}
+ else if (SCM_FRACTIONP (n))
+ {
+ /* Optimize the fraction case by (a/b)^k ==> (a^k)/(b^k), to avoid
+ needless reduction of intermediate products to lowest terms.
+ If a and b have no common factors, then a^k and b^k have no
+ 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. */
+ 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);
@@ -7330,8 +7389,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);
}
@@ -7879,14 +7939,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))
{
@@ -7916,8 +7976,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);
}
@@ -8880,8 +8940,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);
@@ -8982,8 +9043,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 66aa01a..20b7eb1 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4044,6 +4044,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)))
@@ -4317,6 +4320,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: [PATCH 2/5] Add 'round-ash', a rounding arithmetic shift operator --]
[-- Type: text/x-diff, Size: 17578 bytes --]
From 5c9c098176718fdf9bbe675750258babce9ead69 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:35:09 -0500
Subject: [PATCH 2/5] 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
static functions.
(scm_round_ash): New procedure.
(scm_ash): Reimplement in terms of 'left_shift_exact_integer' and
'floor_right_shift_exact_integer'.
* libguile/numbers.h: Add prototype for scm_round_ash. Rename the
second argument of 'scm_ash' 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'.
---
doc/ref/api-data.texi | 42 +++++---
libguile/numbers.c | 226 +++++++++++++++++++++++++++--------------
libguile/numbers.h | 3 +-
test-suite/tests/numbers.test | 114 +++++++++------------
4 files changed, 233 insertions(+), 152 deletions(-)
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index 9da17d8..c376eb9 100644
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1686,19 +1686,15 @@ starts from 0 for the least significant bit.
@end lisp
@end deffn
-@deffn {Scheme Procedure} ash n cnt
-@deffnx {C Function} scm_ash (n, cnt)
-Return @var{n} shifted left by @var{cnt} bits, or shifted right if
-@var{cnt} is negative. This is an ``arithmetic'' shift.
+@deffn {Scheme Procedure} ash n count
+@deffnx {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * 2^@var{count})}.
+@var{n} and @var{count} must be exact integers.
-This is effectively a multiplication by @m{2^{cnt}, 2^@var{cnt}}, 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.)
-
-With @var{n} viewed as an infinite precision twos complement,
-@code{ash} 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"
@@ -1709,6 +1705,28 @@ dropping bits.
@end lisp
@end deffn
+@deffn {Scheme Procedure} round-ash n count
+@deffnx {C Function} scm_round_ash (n, count)
+Return @math{round(@var{n} * 2^@var{count})}.
+@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.
+
+@lisp
+(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\"
+@end lisp
+@end 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 e63c17a..b99a04b 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4791,19 +4791,119 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
}
#undef FUNC_NAME
+/* Efficiently compute (N * 2^COUNT),
+ where N is an exact integer, and COUNT > 0. */
+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");
+}
+
+/* Efficiently compute floor (N / 2^COUNT),
+ where N is an exact integer and COUNT > 0. */
+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");
+}
+
+/* Efficiently compute round (N / 2^COUNT),
+ where N is an exact integer and COUNT > 0. */
+static SCM
+round_right_shift_exact_integer (SCM n, long count)
+{
+ if (SCM_I_INUMP (n))
+ {
+ if (count >= SCM_I_FIXNUM_BIT)
+ return SCM_INUM0;
+ else
+ {
+ scm_t_inum nn = SCM_I_INUM (n);
+ scm_t_inum qq = SCM_SRS (nn, count);
+
+ 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)
+ && (mpz_odd_p (SCM_I_BIG_MPZ (q))
+ || (mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1)))
+ 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"
+ (SCM n, SCM count),
+ "Return @math{floor(@var{n} * 2^@var{count})}.\n"
+ "@var{n} and @var{count} must be exact integers.\n"
"\n"
- "This is effectively a multiplication by 2^@var{cnt}, 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"
- "\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"
@@ -4814,79 +4914,57 @@ 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} * 2^@var{count})}.\n"
+ "@var{n} and @var{count} must be exact 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 2c8b260..912f287 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -206,7 +206,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 20b7eb1..8dab29c 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?
;;;
@@ -4904,3 +4839,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #4: [PATCH 3/5] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp --]
[-- Type: text/x-diff, Size: 7935 bytes --]
From 098aba10f7be8e0a7dc966b8525bef1ca50789e8 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 04:58:55 -0500
Subject: [PATCH 3/5] Simplify and improve scm_i_big2dbl, and add
scm_i_big2dbl_2exp
* libguile/numbers.c (scm_i_big2dbl_2exp): New static function.
(scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp',
with proper rounding.
* test-suite/tests/numbers.test ("exact->inexact"): Add tests.
---
libguile/numbers.c | 101 +++++++++++++++--------------------------
test-suite/tests/numbers.test | 57 +++++++++++++++++------
2 files changed, 80 insertions(+), 78 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index b99a04b..49b4a50 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -330,81 +330,52 @@ 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.
+static SCM round_right_shift_exact_integer (SCM n, long count);
- 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.
+/* 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. */
- 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 DBL_MANT_DIG 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)
+static double
+scm_i_big2dbl_2exp (SCM b, long *expon_p)
{
- 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 */
-
- /* 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)
- {
- size_t shift = bits - DBL_MANT_DIG;
- mpz_init2 (tmp, DBL_MANT_DIG);
- 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
+ size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
+ size_t shift = 0;
if (bits > DBL_MANT_DIG)
{
- unsigned long pos = bits - DBL_MANT_DIG - 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)))
+ shift = bits - DBL_MANT_DIG;
+ b = round_right_shift_exact_integer (b, shift);
+ if (SCM_I_INUMP (b))
{
- result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
+ int expon;
+ double signif = frexp (SCM_I_INUM (b), &expon);
+ *expon_p = expon + shift;
+ return signif;
}
}
- scm_remember_upto_here_1 (b);
- return result;
+ {
+ long expon;
+ double signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b));
+ scm_remember_upto_here_1 (b);
+ *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
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 8dab29c..58ff7b8 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -3848,21 +3848,17 @@
;;;
(with-test-prefix "exact->inexact"
-
+
+ ;; Test "(exact->inexact n)", expect "want".
+ (define (test name n want)
+ (with-test-prefix name
+ (pass-if-equal "pos" want (exact->inexact n))
+ (pass-if-equal "neg" (- want) (exact->inexact (- n)))))
+
;; Test "(exact->inexact n)", expect "want".
;; "i" is a index, for diagnostic purposes.
(define (try-i i n want)
- (with-test-prefix (list i n want)
- (with-test-prefix "pos"
- (let ((got (exact->inexact n)))
- (pass-if "inexact?" (inexact? got))
- (pass-if (list "=" got) (= want got))))
- (set! n (- n))
- (set! want (- want))
- (with-test-prefix "neg"
- (let ((got (exact->inexact n)))
- (pass-if "inexact?" (inexact? got))
- (pass-if (list "=" got) (= want got))))))
+ (test (list i n want) n want))
(with-test-prefix "2^i, no round"
(do ((i 0 (1+ i))
@@ -3935,7 +3931,42 @@
;; convert the num and den to doubles, resulting in infs.
(pass-if "frac big/big, exceeding double"
(let ((big (ash 1 4096)))
- (= 1.0 (exact->inexact (/ (1+ big) big))))))
+ (= 1.0 (exact->inexact (/ (1+ big) big)))))
+
+ (test "round up to odd"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111000101 ->
+ ;; 11111111111111111111111111111111111111111111111111001000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000101)
+ (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b001000))
+
+ (test "round down to odd"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111001011 ->
+ ;; 11111111111111111111111111111111111111111111111111001000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b001011)
+ (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b001000))
+
+ (test "round tie up to even"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111011100 ->
+ ;; 11111111111111111111111111111111111111111111111111100000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b011100)
+ (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b100000))
+
+ (test "round tie down to even"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111000100 ->
+ ;; 11111111111111111111111111111111111111111111111111000000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b000100)
+ (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b000000))
+
+ (test "round tie up to next power of two"
+ ;; =====================================================
+ ;; 11111111111111111111111111111111111111111111111111111100 ->
+ ;; 100000000000000000000000000000000000000000000000000000000
+ (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b111100)
+ (expt 2.0 (+ dbl-mant-dig 3))))
;;;
;;; expt
--
1.7.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #5: [PATCH 4/5] Optimize logarithms using scm_i_big2dbl_2exp --]
[-- Type: text/x-diff, Size: 2354 bytes --]
From ea0421d6937d39c62ac3c10abaa9bad0b230b28a Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Sun, 3 Mar 2013 05:02:53 -0500
Subject: [PATCH 4/5] 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 file changed, 12 insertions(+), 18 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 49b4a50..0b13d69 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9576,26 +9576,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 * 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 */
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 */
@@ -9606,8 +9600,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.10.4
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #6: [PATCH 5/5] Reimplement 'inexact->exact' to avoid mpq functions --]
[-- Type: text/x-diff, Size: 5621 bytes --]
From 6e6991d9f7043aa70dc53c12afb97eeaf3cd63ba Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 4 Mar 2013 18:42:27 -0500
Subject: [PATCH 5/5] Reimplement 'inexact->exact' to avoid mpq functions.
* libguile/numbers.c (scm_inexact_to_exact): Implement conversion of a
double to an exact rational without using the mpq functions.
* test-suite/tests/numbers.test (dbl-mant-dig): Simplify initializer.
(dbl-epsilon, dbl-min-exp): New variables.
("inexact->exact"): Add tests. Fix broken "2.0**i to exact and back"
test, and change it to "2.0**i to exact", to avoid use of
'exact->inexact'.
---
libguile/numbers.c | 41 +++++++++++++--------
test-suite/tests/numbers.test | 80 +++++++++++++++++++++++++++++++++--------
2 files changed, 93 insertions(+), 28 deletions(-)
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 0b13d69..845b37a 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9085,22 +9085,35 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
if (!SCM_LIKELY (DOUBLE_IS_FINITE (val)))
SCM_OUT_OF_RANGE (1, z);
+ else if (val == 0.0)
+ return SCM_INUM0;
else
{
- mpq_t frac;
- SCM q;
-
- mpq_init (frac);
- mpq_set_d (frac, val);
- 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...
- */
- mpq_clear (frac);
- return q;
+ int expon;
+ SCM numerator;
+
+ numerator = scm_i_dbl2big (ldexp (frexp (val, &expon),
+ DBL_MANT_DIG));
+ expon -= DBL_MANT_DIG;
+ if (expon < 0)
+ {
+ int shift = mpz_scan1 (SCM_I_BIG_MPZ (numerator), 0);
+
+ if (shift > -expon)
+ shift = -expon;
+ mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (numerator),
+ SCM_I_BIG_MPZ (numerator),
+ shift);
+ expon += shift;
+ }
+ numerator = scm_i_normbig (numerator);
+ if (expon < 0)
+ return scm_i_make_ratio_already_reduced
+ (numerator, left_shift_exact_integer (SCM_INUM1, -expon));
+ else if (expon > 0)
+ return left_shift_exact_integer (numerator, expon);
+ else
+ return numerator;
}
}
}
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 58ff7b8..93caf04 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -46,15 +46,24 @@
;; the usual 53.
;;
(define dbl-mant-dig
- (let more ((i 1)
- (d 2.0))
- (if (> i 1024)
- (error "Oops, cannot determine number of bits in mantissa of inexact"))
- (let* ((sum (+ 1.0 d))
- (diff (- sum d)))
- (if (= diff 1.0)
- (more (1+ i) (* 2.0 d))
- i))))
+ (do ((prec 0 (+ prec 1))
+ (eps 1.0 (/ eps 2.0)))
+ ((begin (when (> prec 1000000)
+ (error "Unable to determine dbl-mant-dig"))
+ (= 1.0 (+ 1.0 eps)))
+ prec)))
+
+(define dbl-epsilon
+ (expt 0.5 (- dbl-mant-dig 1)))
+
+(define dbl-min-exp
+ (do ((x 1.0 (/ x 2.0))
+ (y (+ 1.0 dbl-epsilon) (/ y 2.0))
+ (e 2 (- e 1)))
+ ((begin (when (< e -100000000)
+ (error "Unable to determine dbl-min-exp"))
+ (= x y))
+ e)))
;; like ash, but working on a flonum
(define (ash-flo x n)
@@ -4241,6 +4250,13 @@
;;;
(with-test-prefix "inexact->exact"
+
+ ;; Test "(inexact->exact f)", expect "want".
+ (define (test name f want)
+ (with-test-prefix name
+ (pass-if-equal "pos" want (inexact->exact f))
+ (pass-if-equal "neg" (- want) (inexact->exact (- f)))))
+
(pass-if (documented? inexact->exact))
(pass-if-exception "+inf" exception:out-of-range
@@ -4251,13 +4267,49 @@
(pass-if-exception "nan" exception:out-of-range
(inexact->exact +nan.0))
-
- (with-test-prefix "2.0**i to exact and back"
+
+ (test "0.0" 0.0 0)
+ (test "small even integer" 72.0 72)
+ (test "small odd integer" 73.0 73)
+
+ (test "largest inexact odd integer"
+ (- (expt 2.0 dbl-mant-dig) 1)
+ (- (expt 2 dbl-mant-dig) 1))
+
+ (test "largest inexact odd integer - 1"
+ (- (expt 2.0 dbl-mant-dig) 2)
+ (- (expt 2 dbl-mant-dig) 2))
+
+ (test "largest inexact odd integer + 3"
+ (+ (expt 2.0 dbl-mant-dig) 2)
+ (+ (expt 2 dbl-mant-dig) 2))
+
+ (test "largest inexact odd integer * 2^48"
+ (* (expt 2.0 48) (- (expt 2.0 dbl-mant-dig) 1))
+ (* (expt 2 48) (- (expt 2 dbl-mant-dig) 1)))
+
+ (test "largest inexact odd integer / 2^48"
+ (* (expt 0.5 48) (- (expt 2.0 dbl-mant-dig) 1))
+ (* (expt 1/2 48) (- (expt 2 dbl-mant-dig) 1)))
+
+ (test "smallest inexact"
+ (expt 2.0 (- dbl-min-exp dbl-mant-dig))
+ (expt 2 (- dbl-min-exp dbl-mant-dig)))
+
+ (test "smallest inexact * 2"
+ (* 2.0 (expt 2.0 (- dbl-min-exp dbl-mant-dig)))
+ (* 2 (expt 2 (- dbl-min-exp dbl-mant-dig))))
+
+ (test "smallest inexact * 3"
+ (* 3.0 (expt 2.0 (- dbl-min-exp dbl-mant-dig)))
+ (* 3 (expt 2 (- dbl-min-exp dbl-mant-dig))))
+
+ (with-test-prefix "2.0**i to exact"
(do ((i 0 (1+ i))
- (n 1.0 (* 2.0 n)))
+ (n 1 (* 2 n))
+ (f 1.0 (* 2.0 f)))
((> i 100))
- (pass-if (list i n)
- (= n (inexact->exact (exact->inexact n)))))))
+ (test (list i n) f n))))
;;;
;;; integer-expt
--
1.7.10.4
^ permalink raw reply related [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-06 20:30 ` Mark H Weaver
@ 2013-03-07 19:02 ` Ludovic Courtès
0 siblings, 0 replies; 8+ messages in thread
From: Ludovic Courtès @ 2013-03-07 19:02 UTC (permalink / raw)
To: Mark H Weaver; +Cc: guile-devel
Mark H Weaver <mhw@netris.org> skribis:
>> What was the rationale for ‘round-ash’?
>
> Well, we need it internally to efficiently convert large integers to
> floating-point with proper rounding (see the call to
> 'round_right_shift_exact_integer' in patch 5). Given that, it seemed
> reasonable to export it to the user, since it's not possible to do that
> job efficiently with our current primitives. However, I don't feel
> strongly about it.
OK, the rationale makes sense to me.
>> It seems to address to do two things differently from ‘ash’: it’s more
>> efficient, and it rounds when right-shifting. The “more efficient” bit
>> is an implementation detail that should also benefit to ‘ash’ (as a user
>> I don’t want to have to choose between efficient and non-rounding.)
>
> No, 'round-ash' is not more efficient than 'ash'. It's more efficient
> than the equivalent (round (* n (expt 2 count))).
I see. I misunderstood the documentation’s reference to efficiency.
Thanks,
Ludo’.
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-07 6:03 ` Mark H Weaver
@ 2013-03-12 20:58 ` Mark H Weaver
0 siblings, 0 replies; 8+ messages in thread
From: Mark H Weaver @ 2013-03-12 20:58 UTC (permalink / raw)
To: guile-devel
I wrote:
> Here are improved versions of the patches needed to enable mini-gmp
> integration. I think these are ready to commit. Reviews welcome.
I went ahead and pushed these first 5 patches. However, I've since
realized that these 5 aren't quite enough for mini-gmp integration. One
of the 3 remaining patches is also needed. I hope to get those ready
soon. They lack only tests.
Thanks,
Mark
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [PATCHES] Numeric improvements
2013-03-05 11:04 [PATCHES] Numeric improvements Mark H Weaver
2013-03-06 14:05 ` Ludovic Courtès
@ 2013-03-17 23:49 ` Mark H Weaver
1 sibling, 0 replies; 8+ messages in thread
From: Mark H Weaver @ 2013-03-17 23:49 UTC (permalink / raw)
To: guile-devel
I've pushed improved versions of these patches to stable-2.0.
Thanks,
Mark
^ permalink raw reply [flat|nested] 8+ messages in thread
end of thread, other threads:[~2013-03-17 23:49 UTC | newest]
Thread overview: 8+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2013-03-05 11:04 [PATCHES] Numeric improvements Mark H Weaver
2013-03-06 14:05 ` Ludovic Courtès
2013-03-06 20:30 ` Mark H Weaver
2013-03-07 19:02 ` Ludovic Courtès
2013-03-07 0:16 ` Mark H Weaver
2013-03-07 6:03 ` Mark H Weaver
2013-03-12 20:58 ` Mark H Weaver
2013-03-17 23:49 ` Mark H Weaver
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).