diff options
| -rw-r--r-- | NEWS | 5 | ||||
| -rw-r--r-- | src/div.c | 22 |
2 files changed, 25 insertions, 2 deletions
@@ -47,8 +47,9 @@ Changes from versions 3.1.* to version 3.2.0: - Improved MPFR manual. - New MPFRbench program (see the tools/bench directory). - Speedup in the mpfr_const_euler function (contributed by Fredrik Johansson), - and in the computation of Bernoulli numbers (used in mpfr_gamma, mpfr_li2, - mpfr_digamma, mpfr_lngamma and mpfr_lgamma). + in the computation of Bernoulli numbers (used in mpfr_gamma, mpfr_li2, + mpfr_digamma, mpfr_lngamma and mpfr_lgamma), and in mpfr_div when the + divisor has one limb. - Bug fixes. In particular, a speed improvement when the --enable-assert or --enable-assert=full configure option is used with GCC. @@ -240,6 +240,28 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) * * **************************************************************************/ + /* when the divisor has one limb, we can use mpfr_div_ui, which should be + faster, assuming there is no intermediate overflow or underflow. + The divisor interpreted as an integer satisfies + 2^(GMP_NUMB_BITS-1) <= vm < 2^GMP_NUMB_BITS, thus the quotient + satisfies 2^(EXP(u)-1-GMP_NUMB_BITS) < u/vm < 2^(EXP(u)-GMP_NUMB_BITS+1) + and its exponent is either EXP(u)-GMP_NUMB_BITS or one more. */ + if (vsize <= 1 && __gmpfr_emin <= MPFR_EXP(u) - GMP_NUMB_BITS + && MPFR_EXP(u) - GMP_NUMB_BITS + 1 <= __gmpfr_emax) + { + mpfr_exp_t exp_v = MPFR_EXP(v); /* save it in case q=v */ + if (MPFR_SIGN(v) > 0) + inex = mpfr_div_ui (q, u, vp[0], rnd_mode); + else + { + inex = -mpfr_div_ui (q, u, vp[0], MPFR_INVERT_RND(rnd_mode)); + MPFR_CHANGE_SIGN(q); + } + /* q did not under/overflow */ + MPFR_EXP(q) -= exp_v - GMP_NUMB_BITS; + return mpfr_check_range (q, inex, rnd_mode); + } + MPFR_TMP_MARK(marker); /* set sign */ |