diff options
| -rw-r--r-- | src/mpfr-mini-gmp.c | 2 | ||||
| -rw-r--r-- | src/mul_ui.c | 119 |
2 files changed, 63 insertions, 58 deletions
diff --git a/src/mpfr-mini-gmp.c b/src/mpfr-mini-gmp.c index 7ae33d5a4..a01520574 100644 --- a/src/mpfr-mini-gmp.c +++ b/src/mpfr-mini-gmp.c @@ -224,6 +224,7 @@ mpn_tdiv_qr (mp_limb_t *qp, mp_limb_t *rp, mp_size_t qxn, } #endif +#if 0 /* this function is useful for debugging, thus please keep it here */ void mpz_dump (mpz_t z) { @@ -264,5 +265,6 @@ mpz_dump (mpz_t z) } printf ("\n"); } +#endif #endif /* MPFR_USE_MINI_GMP */ diff --git a/src/mul_ui.c b/src/mul_ui.c index bf9bf3977..83ff816c1 100644 --- a/src/mul_ui.c +++ b/src/mul_ui.c @@ -27,10 +27,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., MPFR_HOT_FUNCTION_ATTR int mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode) { - mp_limb_t *yp; - mp_size_t xn; - int cnt, inexact; - MPFR_TMP_DECL (marker); + int inexact; if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { @@ -76,60 +73,66 @@ mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode return mpfr_mul_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode); #ifdef MPFR_LONG_WITHIN_LIMB - yp = MPFR_MANT (y); - xn = MPFR_LIMB_SIZE (x); - - MPFR_ASSERTD (xn < MP_SIZE_T_MAX); - MPFR_TMP_MARK(marker); - yp = MPFR_TMP_LIMBS_ALLOC (xn + 1); - - MPFR_ASSERTN (u == (mp_limb_t) u); - yp[xn] = mpn_mul_1 (yp, MPFR_MANT (x), xn, u); - - /* x * u is stored in yp[xn], ..., yp[0] */ - - /* since the case u=1 was treated above, we have u >= 2, thus - yp[xn] >= 1 since x was msb-normalized */ - MPFR_ASSERTD (yp[xn] != 0); - if (MPFR_LIKELY (MPFR_LIMB_MSB (yp[xn]) == 0)) - { - count_leading_zeros (cnt, yp[xn]); - mpn_lshift (yp, yp, xn + 1, cnt); - } - else - { - cnt = 0; - } - - /* now yp[xn], ..., yp[0] is msb-normalized too, and has at most - PREC(x) + (GMP_NUMB_BITS - cnt) non-zero bits */ - MPFR_RNDRAW (inexact, y, yp, (mpfr_prec_t) (xn + 1) * GMP_NUMB_BITS, - rnd_mode, MPFR_SIGN (x), cnt -- ); - - MPFR_TMP_FREE (marker); - - cnt = GMP_NUMB_BITS - cnt; - if (MPFR_UNLIKELY (__gmpfr_emax < MPFR_EMAX_MIN + cnt - || MPFR_GET_EXP (x) > __gmpfr_emax - cnt)) - return mpfr_overflow (y, rnd_mode, MPFR_SIGN(x)); - - MPFR_SET_EXP (y, MPFR_GET_EXP (x) + cnt); - MPFR_SET_SAME_SIGN (y, x); + { + mp_limb_t *yp; + mp_size_t xn; + int cnt; + MPFR_TMP_DECL (marker); + + yp = MPFR_MANT (y); + xn = MPFR_LIMB_SIZE (x); + + MPFR_ASSERTD (xn < MP_SIZE_T_MAX); + MPFR_TMP_MARK(marker); + yp = MPFR_TMP_LIMBS_ALLOC (xn + 1); + + MPFR_ASSERTN (u == (mp_limb_t) u); + yp[xn] = mpn_mul_1 (yp, MPFR_MANT (x), xn, u); + + /* x * u is stored in yp[xn], ..., yp[0] */ + + /* since the case u=1 was treated above, we have u >= 2, thus + yp[xn] >= 1 since x was msb-normalized */ + MPFR_ASSERTD (yp[xn] != 0); + if (MPFR_LIKELY (MPFR_LIMB_MSB (yp[xn]) == 0)) + { + count_leading_zeros (cnt, yp[xn]); + mpn_lshift (yp, yp, xn + 1, cnt); + } + else + { + cnt = 0; + } + + /* now yp[xn], ..., yp[0] is msb-normalized too, and has at most + PREC(x) + (GMP_NUMB_BITS - cnt) non-zero bits */ + MPFR_RNDRAW (inexact, y, yp, (mpfr_prec_t) (xn + 1) * GMP_NUMB_BITS, + rnd_mode, MPFR_SIGN (x), cnt -- ); + + MPFR_TMP_FREE (marker); + + cnt = GMP_NUMB_BITS - cnt; + if (MPFR_UNLIKELY (__gmpfr_emax < MPFR_EMAX_MIN + cnt + || MPFR_GET_EXP (x) > __gmpfr_emax - cnt)) + return mpfr_overflow (y, rnd_mode, MPFR_SIGN(x)); + + MPFR_SET_EXP (y, MPFR_GET_EXP (x) + cnt); + MPFR_SET_SAME_SIGN (y, x); + MPFR_RET (inexact); + } #else - { - mpfr_t uu; - MPFR_SAVE_EXPO_DECL (expo); - - mpfr_init2 (uu, sizeof (unsigned long) * CHAR_BIT); - /* Warning: u might be outside the current exponent range! */ - MPFR_SAVE_EXPO_MARK (expo); - mpfr_set_ui (uu, u, MPFR_RNDZ); - inexact = mpfr_mul (y, x, uu, rnd_mode); - mpfr_clear (uu); - MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range (y, inexact, rnd_mode); - } + { + mpfr_t uu; + MPFR_SAVE_EXPO_DECL (expo); + + mpfr_init2 (uu, sizeof (unsigned long) * CHAR_BIT); + /* Warning: u might be outside the current exponent range! */ + MPFR_SAVE_EXPO_MARK (expo); + mpfr_set_ui (uu, u, MPFR_RNDZ); + inexact = mpfr_mul (y, x, uu, rnd_mode); + mpfr_clear (uu); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (y, inexact, rnd_mode); + } #endif /* MPFR_LONG_WITHIN_LIMB */ - - MPFR_RET (inexact); } |