diff --git a/libguile/integers.c b/libguile/integers.c index a20fdf7db..820f19ddf 100644 --- a/libguile/integers.c +++ b/libguile/integers.c @@ -1,4 +1,4 @@ -/* Copyright 1995-2016,2018-2021 +/* Copyright 1995-2016,2018-2022 Free Software Foundation, Inc. This file is part of Guile. @@ -2102,3 +2102,105 @@ scm_integer_modulo_expt_nnn (SCM n, SCM k, SCM m) return take_mpz (n_tmp); } + +/* Efficiently compute (N * 2^COUNT), where N is an exact integer, and + COUNT > 0. */ +SCM +scm_integer_lsh_iu (scm_t_inum n, unsigned long count) +{ + ASSERT (count > 0); + /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will almost[*] 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)". + + [*] There's one exception: + (-1) << SCM_I_FIXNUM_BIT-1 == SCM_MOST_NEGATIVE_FIXNUM */ + + if (n == 0) + return SCM_I_MAKINUM (n); + else if (count < SCM_I_FIXNUM_BIT-1 && + ((scm_t_bits) (SCM_SRS (n, (SCM_I_FIXNUM_BIT-1 - count)) + 1) + <= 1)) + return SCM_I_MAKINUM (n < 0 ? -(-n << count) : (n << count)); + else + { + mpz_t result; + mpz_init_set_si (result, n); + mpz_mul_2exp (result, result, count); + return take_mpz (result); + } +} + +SCM +scm_integer_lsh_zu (SCM n, unsigned long count) +{ + ASSERT (count > 0); + mpz_t result, zn; + mpz_init (result); + alias_bignum_to_mpz (scm_bignum (n), zn); + mpz_mul_2exp (result, zn, count); + scm_remember_upto_here_1 (n); + return take_mpz (result); +} + +/* Efficiently compute floor (N / 2^COUNT), where N is an exact integer + and COUNT > 0. */ +SCM +scm_integer_floor_rsh_iu (scm_t_inum n, unsigned long count) +{ + ASSERT (count > 0); + if (count >= SCM_I_FIXNUM_BIT) + return (n >= 0 ? SCM_INUM0 : SCM_I_MAKINUM (-1)); + else + return SCM_I_MAKINUM (SCM_SRS (n, count)); +} + +SCM +scm_integer_floor_rsh_zu (SCM n, unsigned long count) +{ + ASSERT (count > 0); + mpz_t result, zn; + mpz_init (result); + alias_bignum_to_mpz (scm_bignum (n), zn); + mpz_fdiv_q_2exp (result, zn, count); + scm_remember_upto_here_1 (n); + return take_mpz (result); +} + +/* Efficiently compute round (N / 2^COUNT), where N is an exact integer + and COUNT > 0. */ +SCM +scm_integer_round_rsh_iu (scm_t_inum n, unsigned long count) +{ + ASSERT (count > 0); + if (count >= SCM_I_FIXNUM_BIT) + return SCM_INUM0; + else + { + scm_t_inum q = SCM_SRS (n, count); + + if (0 == (n & (1L << (count-1)))) + return SCM_I_MAKINUM (q); /* round down */ + else if (n & ((1L << (count-1)) - 1)) + return SCM_I_MAKINUM (q + 1); /* round up */ + else + return SCM_I_MAKINUM ((~1L) & (q + 1)); /* round to even */ + } +} + +SCM +scm_integer_round_rsh_zu (SCM n, unsigned long count) +{ + ASSERT (count > 0); + mpz_t q, zn; + mpz_init (q); + alias_bignum_to_mpz (scm_bignum (n), zn); + mpz_fdiv_q_2exp (q, zn, count); + if (mpz_tstbit (zn, count-1) + && (mpz_odd_p (q) || mpz_scan1 (zn, 0) < count-1)) + mpz_add_ui (q, q, 1); + scm_remember_upto_here_1 (n); + return take_mpz (q); +} diff --git a/libguile/integers.h b/libguile/integers.h index 74624f143..dea4c2235 100644 --- a/libguile/integers.h +++ b/libguile/integers.h @@ -156,6 +156,13 @@ SCM_INTERNAL SCM scm_integer_lognot_z (SCM n); SCM_INTERNAL SCM scm_integer_modulo_expt_nnn (SCM n, SCM k, SCM m); +SCM_INTERNAL SCM scm_integer_lsh_iu (scm_t_inum n, unsigned long count); +SCM_INTERNAL SCM scm_integer_lsh_zu (SCM n, unsigned long count); +SCM_INTERNAL SCM scm_integer_floor_rsh_iu (scm_t_inum n, unsigned long count); +SCM_INTERNAL SCM scm_integer_floor_rsh_zu (SCM n, unsigned long count); +SCM_INTERNAL SCM scm_integer_round_rsh_iu (scm_t_inum n, unsigned long count); +SCM_INTERNAL SCM scm_integer_round_rsh_zu (SCM n, unsigned long count); + #endif /* SCM_INTEGERS_H */ diff --git a/libguile/numbers.c b/libguile/numbers.c index ab7a659c8..46f7b21d2 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -387,7 +387,7 @@ scm_i_dbl2num (double u) return scm_i_dbl2big (u); } -static SCM round_right_shift_exact_integer (SCM n, long count); +static SCM round_rsh (SCM n, SCM count); /* scm_i_big2dbl_2exp() is like frexp for bignums: it converts the bignum b into a normalized significand and exponent such that @@ -406,7 +406,7 @@ scm_i_big2dbl_2exp (SCM b, long *expon_p) if (bits > DBL_MANT_DIG) { shift = bits - DBL_MANT_DIG; - b = round_right_shift_exact_integer (b, shift); + b = round_rsh (b, scm_from_size_t (shift)); if (SCM_I_INUMP (b)) { int expon; @@ -3229,118 +3229,52 @@ scm_integer_expt (SCM n, SCM k) return scm_call_2 (scm_variable_ref (integer_expt_var), n, k); } -/* Efficiently compute (N * 2^COUNT), - where N is an exact integer, and COUNT > 0. */ static SCM -left_shift_exact_integer (SCM n, long count) +lsh (SCM n, SCM count, const char *fn) { + if (scm_is_eq (n, SCM_INUM0)) + return n; + if (!scm_is_unsigned_integer (count, 0, ULONG_MAX)) + scm_num_overflow (fn); + + unsigned long ucount = scm_to_ulong (count); + if (ucount == 0) + return n; + if (ucount / (sizeof (int) * 8) >= (unsigned long) INT_MAX) + scm_num_overflow (fn); if (SCM_I_INUMP (n)) - { - scm_t_inum nn = SCM_I_INUM (n); - - /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will almost[*] 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)". - - [*] There's one exception: - (-1) << SCM_I_FIXNUM_BIT-1 == SCM_MOST_NEGATIVE_FIXNUM */ - - 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 < 0 ? -(-nn << count) : (nn << count)); - else - { - SCM result = scm_i_inum2big (nn); - mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), - count); - return scm_i_normbig (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 - assert (0); + return scm_integer_lsh_iu (SCM_I_INUM (n), ucount); + return scm_integer_lsh_zu (n, ucount); } -/* 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) +floor_rsh (SCM n, SCM count) { - if (SCM_I_INUMP (n)) - { - scm_t_inum nn = SCM_I_INUM (n); + if (!scm_is_unsigned_integer (count, 0, ULONG_MAX)) + return scm_is_false (scm_negative_p (n)) ? SCM_INUM0 : SCM_I_MAKINUM (-1); - 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 - assert (0); + unsigned long ucount = scm_to_ulong (count); + if (ucount == 0) + return n; + if (SCM_I_INUMP (n)) + return scm_integer_floor_rsh_iu (SCM_I_INUM (n), ucount); + return scm_integer_floor_rsh_zu (n, ucount); } -/* 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) +round_rsh (SCM n, SCM count) { + if (!scm_is_unsigned_integer (count, 0, ULONG_MAX)) + return SCM_INUM0; + + unsigned long ucount = scm_to_ulong (count); + if (ucount == 0) + return n; 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 - assert (0); + return scm_integer_round_rsh_iu (SCM_I_INUM (n), ucount); + return scm_integer_round_rsh_zu (n, ucount); } -/* 'scm_ash' and 'scm_round_ash' assume that fixnums fit within a long, - and moreover that they can be negated without overflow. */ -verify (SCM_MOST_NEGATIVE_FIXNUM >= LONG_MIN + 1 - && SCM_MOST_POSITIVE_FIXNUM <= LONG_MAX); - SCM_DEFINE (scm_ash, "ash", 2, 0, 0, (SCM n, SCM count), "Return @math{floor(@var{n} * 2^@var{count})}.\n" @@ -3360,30 +3294,15 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_ash { - if (!SCM_I_INUMP (n) && !SCM_BIGP (n)) + if (!scm_is_exact_integer (n)) SCM_WRONG_TYPE_ARG (SCM_ARG1, n); + if (!scm_is_exact_integer (count)) + SCM_WRONG_TYPE_ARG (SCM_ARG2, count); - if (scm_is_false (scm_positive_p (scm_sum (scm_integer_length (n), - count)))) - /* Huge right shift that eliminates all but the sign bit */ - return scm_is_false (scm_negative_p (n)) - ? SCM_INUM0 : SCM_I_MAKINUM (-1); - else if (scm_is_true (scm_zero_p (n))) - return SCM_INUM0; - else if (scm_is_signed_integer (count, INT32_MIN + 1, INT32_MAX)) - { - /* We exclude MIN to ensure that 'bits_to_shift' can be - negated without overflowing, if INT32_MIN happens to be LONG_MIN */ - long bits_to_shift = scm_to_long (count); - if (bits_to_shift > 0) - 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 - return n; - } - else - scm_num_overflow ("ash"); + if (scm_is_false (scm_negative_p (count))) + return lsh (n, count, "ash"); + + return floor_rsh (n, scm_difference (count, SCM_UNDEFINED)); } #undef FUNC_NAME @@ -3409,29 +3328,15 @@ SCM_DEFINE (scm_round_ash, "round-ash", 2, 0, 0, "@end lisp") #define FUNC_NAME s_scm_round_ash { - if (!SCM_I_INUMP (n) && !SCM_BIGP (n)) + if (!scm_is_exact_integer (n)) SCM_WRONG_TYPE_ARG (SCM_ARG1, n); + if (!scm_is_exact_integer (count)) + SCM_WRONG_TYPE_ARG (SCM_ARG2, count); - if (scm_is_true (scm_negative_p (scm_sum (scm_integer_length (n), - count))) - || scm_is_true (scm_zero_p (n))) - /* If N is zero, or the right shift count exceeds the integer - length, the result is zero. */ - return SCM_INUM0; - else if (scm_is_signed_integer (count, INT32_MIN + 1, INT32_MAX)) - { - /* We exclude MIN to ensure that 'bits_to_shift' can be - negated without overflowing, if INT32_MIN happens to be LONG_MIN */ - long bits_to_shift = scm_to_long (count); - 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 - return n; - } - else - scm_num_overflow ("round-ash"); + if (scm_is_false (scm_negative_p (count))) + return lsh (n, count, "round-ash"); + + return round_rsh (n, scm_difference (count, SCM_UNDEFINED)); } #undef FUNC_NAME @@ -7696,9 +7601,9 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0, numerator = scm_i_normbig (numerator); if (expon < 0) return scm_i_make_ratio_already_reduced - (numerator, left_shift_exact_integer (SCM_INUM1, -expon)); + (numerator, scm_integer_lsh_iu (1, -expon)); else if (expon > 0) - return left_shift_exact_integer (numerator, expon); + return lsh (numerator, scm_from_int (expon), FUNC_NAME); else return numerator; } @@ -8621,9 +8526,9 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0, shift = (scm_to_long (scm_integer_length (n)) - scm_to_long (scm_integer_length (d))) / 2; if (shift > 0) - d = left_shift_exact_integer (d, 2 * shift); + d = lsh (d, scm_from_long (2 * shift), FUNC_NAME); else - n = left_shift_exact_integer (n, -2 * shift); + n = lsh (n, scm_from_long (-2 * shift), FUNC_NAME); xx = scm_i_divide2double (n, d); }