[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Guile-commits] GNU Guile branch, stable-2.0, updated. v2.0.7-193-g24475
|
From: |
Mark H Weaver |
|
Subject: |
[Guile-commits] GNU Guile branch, stable-2.0, updated. v2.0.7-193-g24475b8 |
|
Date: |
Tue, 12 Mar 2013 20:50:46 +0000 |
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU Guile".
http://git.savannah.gnu.org/cgit/guile.git/commit/?id=24475b860b02880b1cfdf4e03f9659a8af09eb72
The branch, stable-2.0 has been updated
via 24475b860b02880b1cfdf4e03f9659a8af09eb72 (commit)
via 7f34acd8a48198c7fec2daf8d2f4161eaa9963ec (commit)
via 1eb6a33a30ea27f97fc401a25a3014e10e3c6f98 (commit)
via e08a12b5356c20ed0418bcaee136eb3632c5616f (commit)
via a285b18ca820e089e2e5d02f8ed07a1e341dffc3 (commit)
from d2df3950a905f7acab70633717beddfd90455b68 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 24475b860b02880b1cfdf4e03f9659a8af09eb72
Author: Mark H Weaver <address@hidden>
Date: Mon Mar 4 18:42:27 2013 -0500
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'.
commit 7f34acd8a48198c7fec2daf8d2f4161eaa9963ec
Author: Mark H Weaver <address@hidden>
Date: Sun Mar 3 05:02:53 2013 -0500
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'.
commit 1eb6a33a30ea27f97fc401a25a3014e10e3c6f98
Author: Mark H Weaver <address@hidden>
Date: Sun Mar 3 04:58:55 2013 -0500
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.
commit e08a12b5356c20ed0418bcaee136eb3632c5616f
Author: Mark H Weaver <address@hidden>
Date: Sun Mar 3 04:35:09 2013 -0500
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'.
commit a285b18ca820e089e2e5d02f8ed07a1e341dffc3
Author: Mark H Weaver <address@hidden>
Date: Sun Mar 3 04:34:50 2013 -0500
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.
-----------------------------------------------------------------------
Summary of changes:
doc/ref/api-data.texi | 42 ++-
libguile/numbers.c | 634 ++++++++++++++++++++++++-----------------
libguile/numbers.h | 3 +-
test-suite/tests/numbers.test | 257 +++++++++++------
4 files changed, 573 insertions(+), 363 deletions(-)
diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index fb12d2c..81c6d5b 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 bb1ecf5..f0f7236 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
@@ -442,96 +413,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_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
@@ -823,8 +754,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 +824,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 +4685,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);
@@ -4732,19 +4762,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"
@@ -4755,79 +4885,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
@@ -7354,8 +7462,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);
}
@@ -7903,14 +8012,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))
{
@@ -7940,8 +8049,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);
}
@@ -8904,8 +9013,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);
@@ -8999,21 +9109,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 (scm_i_mpz2num (mpq_numref (frac)),
- scm_i_mpz2num (mpq_denref (frac)));
+ int expon;
+ SCM numerator;
- /* When scm_i_make_ratio throws, we leak the memory allocated
- for frac...
- */
- mpq_clear (frac);
- return q;
+ 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;
}
}
}
@@ -9489,26 +9613,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 */
@@ -9519,8 +9637,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))));
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 be378b7..550dc50 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)
@@ -201,71 +210,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?
;;;
@@ -3923,21 +3867,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))
@@ -4010,7 +3950,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
@@ -4054,6 +4029,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)))
@@ -4282,6 +4260,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
@@ -4292,13 +4277,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
@@ -4327,6 +4348,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))))
@@ -4908,3 +4932,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))
hooks/post-receive
--
GNU Guile
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Guile-commits] GNU Guile branch, stable-2.0, updated. v2.0.7-193-g24475b8,
Mark H Weaver <=