[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [PATCHES] Numeric improvements
|
From: |
Mark H Weaver |
|
Subject: |
Re: [PATCHES] Numeric improvements |
|
Date: |
Thu, 07 Mar 2013 01:03:34 -0500 |
|
User-agent: |
Gnus/5.13 (Gnus v5.13) Emacs/24.3 (gnu/linux) |
Here are improved versions of the patches needed to enable mini-gmp
integration. I think these are ready to commit. Reviews welcome.
Mark
>From f32e8c5ffd789a6dbee48be74f5bbf32978382c3 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
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
>From 5c9c098176718fdf9bbe675750258babce9ead69 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
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
address@hidden {Scheme Procedure} ash n cnt
address@hidden {C Function} scm_ash (n, cnt)
-Return @var{n} shifted left by @var{cnt} bits, or shifted right if
address@hidden is negative. This is an ``arithmetic'' shift.
address@hidden {Scheme Procedure} ash n count
address@hidden {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * address@hidden)}.
address@hidden and @var{count} must be exact integers.
-This is effectively a multiplication by @m{2^{cnt}, address@hidden, and
-when @var{cnt} is negative it's a division, rounded towards negative
-infinity. (Note that this is not the same rounding as @code{quotient}
-does.)
-
-With @var{n} viewed as an infinite precision twos complement,
address@hidden means a left shift introducing zero bits, or a right shift
-dropping bits.
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{ash} means a left shift introducing zero bits
+when @var{count} is positive, or a right shift dropping bits
+when @var{count} is negative. This is an ``arithmetic'' shift.
@lisp
(number->string (ash #b1 3) 2) @result{} "1000"
@@ -1709,6 +1705,28 @@ dropping bits.
@end lisp
@end deffn
address@hidden {Scheme Procedure} round-ash n count
address@hidden {C Function} scm_round_ash (n, count)
+Return @math{round(@var{n} * address@hidden)}.
address@hidden and @var{count} must be exact integers.
+
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{round-ash} means a left shift introducing zero
+bits when @var{count} is positive, or a right shift rounding
+to the nearest integer (with ties going to the nearest even
+integer) when @var{count} is negative. This is a rounded
+``arithmetic'' shift.
+
address@hidden
+(number->string (round-ash #b1 3) 2) @result{} \"1000\"
+(number->string (round-ash #b1010 -1) 2) @result{} \"101\"
+(number->string (round-ash #b1010 -2) 2) @result{} \"10\"
+(number->string (round-ash #b1011 -2) 2) @result{} \"11\"
+(number->string (round-ash #b1101 -2) 2) @result{} \"11\"
+(number->string (round-ash #b1110 -2) 2) @result{} \"100\"
address@hidden lisp
address@hidden deffn
+
@deffn {Scheme Procedure} logcount n
@deffnx {C Function} scm_logcount (n)
Return the number of bits in integer @var{n}. If @var{n} is
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 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} * address@hidden)}.\n"
+ "@var{n} and @var{count} must be exact integers.\n"
"\n"
- "This is effectively a multiplication by address@hidden, and when\n"
- "@var{cnt} is negative it's a division, rounded towards negative\n"
- "infinity. (Note that this is not the same rounding as\n"
- "@code{quotient} does.)\n"
- "\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} * address@hidden)}.\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
>From 098aba10f7be8e0a7dc966b8525bef1ca50789e8 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
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
>From ea0421d6937d39c62ac3c10abaa9bad0b230b28a Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
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
>From 6e6991d9f7043aa70dc53c12afb97eeaf3cd63ba Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
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