diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/div.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/div.c | 861 |
1 files changed, 807 insertions, 54 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/div.c b/Build/source/libs/mpfr/mpfr-src/src/div.c index 120df1bbfb7..7d208621f24 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/div.c +++ b/Build/source/libs/mpfr/mpfr-src/src/div.c @@ -24,30 +24,639 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., [1] Short Division of Long Integers, David Harvey and Paul Zimmermann, Proceedings of the 20th Symposium on Computer Arithmetic (ARITH-20), July 25-27, 2011, pages 7-14. + [2] Improved Division by Invariant Integers, Niels Möller and Torbjörn Granlund, + IEEE Transactions on Computers, volume 60, number 2, pages 165-175, 2011. */ #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -#ifdef DEBUG2 -#define mpfr_mpn_print(ap,n) mpfr_mpn_print3 (ap,n,MPFR_LIMB_ZERO) +#if !defined(MPFR_GENERIC_ABI) + +#if GMP_NUMB_BITS == 64 + +#include "invert_limb.h" + +/* Given u = u1*B+u0 < v = v1*B+v0 with v normalized (high bit of v1 set), + put in q = Q1*B+Q0 an approximation of floor(u*B^2/v), with: + B = 2^GMP_NUMB_BITS and q <= floor(u*B^2/v) <= q + 21. + Note: this function requires __gmpfr_invert_limb_approx (from invert_limb.h) + which is only provided so far for 64-bit limb. + Note: __gmpfr_invert_limb_approx can be replaced by __gmpfr_invert_limb, + in that case the bound 21 reduces to 16. */ static void -mpfr_mpn_print3 (mpfr_limb_ptr ap, mp_size_t n, mp_limb_t cy) +mpfr_div2_approx (mpfr_limb_ptr Q1, mpfr_limb_ptr Q0, + mp_limb_t u1, mp_limb_t u0, + mp_limb_t v1, mp_limb_t v0) +{ + mp_limb_t inv, q1, q0, r1, r0, cy, xx, yy; + + /* first compute an approximation of q1, using a lower approximation of + B^2/(v1+1) - B */ + if (MPFR_UNLIKELY(v1 == MPFR_LIMB_MAX)) + inv = MPFR_LIMB_ZERO; + else + __gmpfr_invert_limb_approx (inv, v1 + 1); + /* now inv <= B^2/(v1+1) - B */ + umul_ppmm (q1, q0, u1, inv); + q1 += u1; + /* now q1 <= u1*B/(v1+1) < (u1*B+u0)*B/(v1*B+v0) */ + + /* compute q1*(v1*B+v0) into r1:r0:yy and subtract from u1:u0:0 */ + umul_ppmm (r1, r0, q1, v1); + umul_ppmm (xx, yy, q1, v0); + + ADD_LIMB (r0, xx, cy); + r1 += cy; + + /* we ignore yy below, but first increment r0, to ensure we get a lower + approximation of the remainder */ + r0 += yy != 0; + r1 += r0 == 0 && yy != 0; + r0 = u0 - r0; + r1 = u1 - r1 - (r0 > u0); + + /* r1:r0 should be nonnegative */ + MPFR_ASSERTD((r1 & MPFR_LIMB_HIGHBIT) == 0); + + /* the second quotient limb is approximated by (r1*B^2+r0*B) / v1, + and since (B+inv)/B approximates B/v1, this is in turn approximated + by (r1*B+r0)*(B+inv)/B = r1*B*r1*inv+r0+(r0*inv/B) */ + + q0 = r0; + q1 += r1; + /* add floor(r0*inv/B) to q0 */ + umul_ppmm (xx, yy, r0, inv); + ADD_LIMB (q0, xx, cy); + q1 += cy; + MPFR_ASSERTD (r1 <= 4); + /* TODO: use value coverage on r1 to check that the 5 cases are tested. */ + while (r1) /* the number of loops is at most 4 */ + { + /* add inv to q0 */ + ADD_LIMB (q0, inv, cy); + q1 += cy; + r1 --; + } + + *Q1 = q1; + *Q0 = q0; +} + +#endif /* GMP_NUMB_BITS == 64 */ + +/* Special code for PREC(q) = PREC(u) = PREC(v) = p < GMP_NUMB_BITS */ +static int +mpfr_div_1 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) { - mp_size_t i; - for (i = 0; i < n; i++) - printf ("+%lu*2^%lu", (unsigned long) ap[i], (unsigned long) - (GMP_NUMB_BITS * i)); - if (cy) - printf ("+2^%lu", (unsigned long) (GMP_NUMB_BITS * n)); - printf ("\n"); + mpfr_prec_t p = MPFR_GET_PREC(q); + mpfr_limb_ptr qp = MPFR_MANT(q); + mpfr_exp_t qx = MPFR_GET_EXP(u) - MPFR_GET_EXP(v); + mpfr_prec_t sh = GMP_NUMB_BITS - p; + mp_limb_t u0 = MPFR_MANT(u)[0]; + mp_limb_t v0 = MPFR_MANT(v)[0]; + mp_limb_t q0, rb, sb, mask = MPFR_LIMB_MASK(sh); + int extra; + + if ((extra = (u0 >= v0))) + u0 -= v0; + +#if GMP_NUMB_BITS == 64 /* __gmpfr_invert_limb_approx only exists for 64-bit */ + /* First try with an approximate quotient. + FIXME: for p<=62 we have sh-1<2 and will never be able to round correctly. + Even for p=61 we have sh-1=2 and we can round correctly only when the two + last bist of q0 are 01, which happens with probability 25% only. */ + { + mp_limb_t inv; + __gmpfr_invert_limb_approx (inv, v0); + umul_ppmm (rb, sb, u0, inv); + } + rb += u0; + q0 = rb >> extra; + /* rb does not exceed the true quotient floor(u0*2^GMP_NUMB_BITS/v0), + with error at most 2, which means the rational quotient q satisfies + rb <= q < rb + 3. We can round correctly except when the last sh-1 bits + of q0 are 000..000 or 111..111 or 111..110. */ + if (MPFR_LIKELY(((q0 + 2) & (mask >> 1)) > 2)) + { + rb = q0 & (MPFR_LIMB_ONE << (sh - 1)); + sb = 1; /* result cannot be exact in this case */ + } + else /* the true quotient is rb, rb+1 or rb+2 */ + { + mp_limb_t h, l; + q0 = rb; + umul_ppmm (h, l, q0, v0); + MPFR_ASSERTD(h < u0 || (h == u0 && l == MPFR_LIMB_ZERO)); + /* subtract {h,l} from {u0,0} */ + sub_ddmmss (h, l, u0, 0, h, l); + /* the remainder {h, l} should be < v0 */ + /* This while loop is executed at most two times, but does not seem + slower than two consecutive identical if-statements. */ + while (h || l >= v0) + { + q0 ++; + h -= (l < v0); + l -= v0; + } + MPFR_ASSERTD(h == 0 && l < v0); + sb = l | (q0 & extra); + q0 >>= extra; + rb = q0 & (MPFR_LIMB_ONE << (sh - 1)); + sb |= q0 & (mask >> 1); + } +#else + udiv_qrnnd (q0, sb, u0, 0, v0); + sb |= q0 & extra; + q0 >>= extra; + rb = q0 & (MPFR_LIMB_ONE << (sh - 1)); + sb |= q0 & (mask >> 1); +#endif + + qp[0] = (MPFR_LIMB_HIGHBIT | q0) & ~mask; + qx += extra; + MPFR_SIGN(q) = MPFR_MULT_SIGN (MPFR_SIGN (u), MPFR_SIGN (v)); + + /* rounding */ + if (MPFR_UNLIKELY(qx > __gmpfr_emax)) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + + /* Warning: underflow should be checked *after* rounding, thus when rounding + away and when q > 0.111...111*2^(emin-1), or when rounding to nearest and + q >= 0.111...111[1]*2^(emin-1), there is no underflow. */ + if (MPFR_UNLIKELY(qx < __gmpfr_emin)) + { + /* Note: the case 0.111...111*2^(emin-1) < q < 2^(emin-1) is not possible + here since (up to exponent) this would imply 1 - 2^(-p) < u/v < 1, + thus v - 2^(-p)*v < u < v, and since we can assume 1/2 <= v < 1, it + would imply v - 2^(-p) = v - ulp(v) < u < v, which has no solution. */ + + /* For RNDN, mpfr_underflow always rounds away, thus for |q|<=2^(emin-2) + we have to change to RNDZ. This corresponds to: + (a) either qx < emin - 1 + (b) or qx = emin - 1 and qp[0] = 1000....000 and rb = sb = 0. + Note: in case (b), it suffices to check whether sb = 0, since rb = 1 + and sb = 0 is not possible (the exact quotient would have p+1 bits, + thus u would need at least p+1 bits). */ + if (rnd_mode == MPFR_RNDN && + (qx < __gmpfr_emin - 1 || (qp[0] == MPFR_LIMB_HIGHBIT && sb == 0))) + rnd_mode = MPFR_RNDZ; + return mpfr_underflow (q, rnd_mode, MPFR_SIGN(q)); + } + + MPFR_EXP (q) = qx; /* Don't use MPFR_SET_EXP since qx might be < __gmpfr_emin + in the cases "goto rounding" above. */ + if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) + { + MPFR_ASSERTD(qx >= __gmpfr_emin); + MPFR_RET (0); + } + else if (rnd_mode == MPFR_RNDN) + { + /* It is not possible to have rb <> 0 and sb = 0 here, since it would + mean a n-bit by n-bit division gives an exact (n+1)-bit number. + And since the case rb = sb = 0 was already dealt with, we cannot + have sb = 0. Thus we cannot be in the middle of two numbers. */ + MPFR_ASSERTD(sb != 0); + if (rb == 0) + goto truncate; + else + goto add_one_ulp; + } + else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(q))) + { + truncate: + MPFR_ASSERTD(qx >= __gmpfr_emin); + MPFR_RET(-MPFR_SIGN(q)); + } + else /* round away from zero */ + { + add_one_ulp: + qp[0] += MPFR_LIMB_ONE << sh; + if (qp[0] == 0) + { + qp[0] = MPFR_LIMB_HIGHBIT; + if (MPFR_UNLIKELY(qx + 1 > __gmpfr_emax)) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + MPFR_ASSERTD(qx + 1 <= __gmpfr_emax); + MPFR_ASSERTD(qx + 1 >= __gmpfr_emin); + MPFR_SET_EXP (q, qx + 1); + } + MPFR_RET(MPFR_SIGN(q)); + } } + +/* Special code for PREC(q) = GMP_NUMB_BITS, + with PREC(u), PREC(v) <= GMP_NUMB_BITS. */ +static int +mpfr_div_1n (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) +{ + mpfr_limb_ptr qp = MPFR_MANT(q); + mpfr_exp_t qx = MPFR_GET_EXP(u) - MPFR_GET_EXP(v); + mp_limb_t u0 = MPFR_MANT(u)[0]; + mp_limb_t v0 = MPFR_MANT(v)[0]; + mp_limb_t q0, rb, sb, l; + int extra; + + MPFR_ASSERTD(MPFR_PREC(q) == GMP_NUMB_BITS); + MPFR_ASSERTD(MPFR_PREC(u) <= GMP_NUMB_BITS); + MPFR_ASSERTD(MPFR_PREC(v) <= GMP_NUMB_BITS); + + if ((extra = (u0 >= v0))) + u0 -= v0; + +#if GMP_NUMB_BITS == 64 /* __gmpfr_invert_limb_approx only exists for 64-bit */ + { + mp_limb_t inv, h; + + /* First compute an approximate quotient. */ + __gmpfr_invert_limb_approx (inv, v0); + umul_ppmm (rb, sb, u0, inv); + q0 = u0 + rb; + /* rb does not exceed the true quotient floor(u0*2^GMP_NUMB_BITS/v0), + with error at most 2, which means the rational quotient q satisfies + rb <= q < rb + 3, thus the true quotient is rb, rb+1 or rb+2 */ + umul_ppmm (h, l, q0, v0); + MPFR_ASSERTD(h < u0 || (h == u0 && l == MPFR_LIMB_ZERO)); + /* subtract {h,l} from {u0,0} */ + sub_ddmmss (h, l, u0, 0, h, l); + /* the remainder {h, l} should be < v0 */ + /* This while loop is executed at most two times, but does not seem + slower than two consecutive identical if-statements. */ + while (h || l >= v0) + { + q0 ++; + h -= (l < v0); + l -= v0; + } + MPFR_ASSERTD(h == 0 && l < v0); + } +#else + udiv_qrnnd (q0, l, u0, 0, v0); +#endif + + /* now (u0 - extra*v0) * 2^GMP_NUMB_BITS = q0*v0 + l with 0 <= l < v0 */ + + /* If extra=0, the quotient is q0, the round bit is 1 if l >= v0/2, + and sb are the remaining bits from l. + If extra=1, the quotient is MPFR_LIMB_HIGHBIT + (q0 >> 1), the round bit + is the least significant bit of q0, and sb is l. */ + + if (extra == 0) + { + qp[0] = q0; + /* If "l + l < l", then there is a carry in l + l, thus 2*l > v0. + Otherwise if there is no carry, we check whether 2*l >= v0. */ + rb = (l + l < l) || (l + l >= v0); + sb = (rb) ? l + l - v0 : l; + } + else + { + qp[0] = MPFR_LIMB_HIGHBIT | (q0 >> 1); + rb = q0 & MPFR_LIMB_ONE; + sb = l; + qx ++; + } + + MPFR_SIGN(q) = MPFR_MULT_SIGN (MPFR_SIGN (u), MPFR_SIGN (v)); + + /* rounding */ + if (MPFR_UNLIKELY(qx > __gmpfr_emax)) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + + /* Warning: underflow should be checked *after* rounding, thus when rounding + away and when q > 0.111...111*2^(emin-1), or when rounding to nearest and + q >= 0.111...111[1]*2^(emin-1), there is no underflow. */ + if (MPFR_UNLIKELY(qx < __gmpfr_emin)) + { + /* Note: the case 0.111...111*2^(emin-1) < q < 2^(emin-1) is not possible + here since (up to exponent) this would imply 1 - 2^(-p) < u/v < 1, + thus v - 2^(-p)*v < u < v, and since we can assume 1/2 <= v < 1, it + would imply v - 2^(-p) = v - ulp(v) < u < v, which has no solution. */ + + /* For RNDN, mpfr_underflow always rounds away, thus for |q|<=2^(emin-2) + we have to change to RNDZ. This corresponds to: + (a) either qx < emin - 1 + (b) or qx = emin - 1 and qp[0] = 1000....000 and rb = sb = 0. + Note: in case (b), it suffices to check whether sb = 0, since rb = 1 + and sb = 0 is not possible (the exact quotient would have p+1 bits, + thus u would need at least p+1 bits). */ + if (rnd_mode == MPFR_RNDN && + (qx < __gmpfr_emin - 1 || (qp[0] == MPFR_LIMB_HIGHBIT && sb == 0))) + rnd_mode = MPFR_RNDZ; + return mpfr_underflow (q, rnd_mode, MPFR_SIGN(q)); + } + + MPFR_EXP (q) = qx; /* Don't use MPFR_SET_EXP since qx might be < __gmpfr_emin + in the cases "goto rounding" above. */ + if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) + { + MPFR_ASSERTD(qx >= __gmpfr_emin); + MPFR_RET (0); + } + else if (rnd_mode == MPFR_RNDN) + { + /* It is not possible to have rb <> 0 and sb = 0 here, since it would + mean a n-bit by n-bit division gives an exact (n+1)-bit number. + And since the case rb = sb = 0 was already dealt with, we cannot + have sb = 0. Thus we cannot be in the middle of two numbers. */ + MPFR_ASSERTD(sb != 0); + if (rb == 0) + goto truncate; + else + goto add_one_ulp; + } + else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(q))) + { + truncate: + MPFR_ASSERTD(qx >= __gmpfr_emin); + MPFR_RET(-MPFR_SIGN(q)); + } + else /* round away from zero */ + { + add_one_ulp: + qp[0] += MPFR_LIMB_ONE; + if (qp[0] == 0) + { + qp[0] = MPFR_LIMB_HIGHBIT; + if (MPFR_UNLIKELY(qx + 1 > __gmpfr_emax)) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + MPFR_ASSERTD(qx + 1 <= __gmpfr_emax); + MPFR_ASSERTD(qx + 1 >= __gmpfr_emin); + MPFR_SET_EXP (q, qx + 1); + } + MPFR_RET(MPFR_SIGN(q)); + } +} + +/* Special code for GMP_NUMB_BITS < PREC(q) < 2*GMP_NUMB_BITS and + PREC(u) = PREC(v) = PREC(q) */ +static int +mpfr_div_2 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) +{ + mpfr_prec_t p = MPFR_GET_PREC(q); + mpfr_limb_ptr qp = MPFR_MANT(q); + mpfr_exp_t qx = MPFR_GET_EXP(u) - MPFR_GET_EXP(v); + mpfr_prec_t sh = 2*GMP_NUMB_BITS - p; + mp_limb_t h, rb, sb, mask = MPFR_LIMB_MASK(sh); + mp_limb_t v1 = MPFR_MANT(v)[1], v0 = MPFR_MANT(v)[0]; + mp_limb_t q1, q0, r3, r2, r1, r0, l, t; + int extra; + + r3 = MPFR_MANT(u)[1]; + r2 = MPFR_MANT(u)[0]; + extra = r3 > v1 || (r3 == v1 && r2 >= v0); + if (extra) + sub_ddmmss (r3, r2, r3, r2, v1, v0); + + MPFR_ASSERTD(r3 < v1 || (r3 == v1 && r2 < v0)); + +#if GMP_NUMB_BITS == 64 + mpfr_div2_approx (&q1, &q0, r3, r2, v1, v0); + /* we know q1*B+q0 is smaller or equal to the exact quotient, with + difference at most 21 */ + if (MPFR_LIKELY(((q0 + 21) & (mask >> 1)) > 21)) + sb = 1; /* result is not exact when we can round with an approximation */ + else + { + /* we know q1:q0 is a good-enough approximation, use it! */ + mp_limb_t s0, s1, s2, h, l; + + /* Since we know the difference should be at most 21*(v1:v0) after the + subtraction below, thus at most 21*2^128, it suffices to compute the + lower 3 limbs of (q1:q0) * (v1:v0). */ + umul_ppmm (s1, s0, q0, v0); + umul_ppmm (s2, l, q0, v1); + s1 += l; + s2 += (s1 < l); + umul_ppmm (h, l, q1, v0); + s1 += l; + s2 += h + (s1 < l); + s2 += q1 * v1; + /* Subtract s2:s1:s0 from r2:0:0, with result in s2:s1:s0. */ + s2 = r2 - s2; + /* now negate s1:s0 */ + s0 = -s0; + s1 = -s1 - (s0 != 0); + /* there is a borrow in s2 when s0 and s1 are not both zero */ + s2 -= (s1 != 0 || s0 != 0); + while (s2 > 0 || (s1 > v1) || (s1 == v1 && s0 >= v0)) + { + /* add 1 to q1:q0 */ + q0 ++; + q1 += (q0 == 0); + /* subtract v1:v0 to s2:s1:s0 */ + s2 -= (s1 < v1) || (s1 == v1 && s0 < v0); + sub_ddmmss (s1, s0, s1, s0, v1, v0); + } + sb = s1 | s0; + } + goto round_div2; +#endif + + /* now r3:r2 < v1:v0 */ + if (MPFR_UNLIKELY(r3 == v1)) /* can occur in some rare cases */ + { + /* This can only occur in case extra=0, since otherwise we would have + u_old >= u_new + v >= B^2/2 + B^2/2 = B^2. In this case we have + r3 = u1 and r2 = u0, thus the remainder u*B-q1*v is + v1*B^2+u0*B-(B-1)*(v1*B+v0) = (u0-v0+v1)*B+v0. + Warning: in this case q1 = B-1 can be too large, for example with + u = B^2/2 and v = B^2/2 + B - 1, then u*B-(B-1)*u = -1/2*B^2+2*B-1. */ + MPFR_ASSERTD(extra == 0); + q1 = MPFR_LIMB_MAX; + r1 = v0; + t = v0 - r2; /* t > 0 since r3:r2 < v1:v0 */ + r2 = v1 - t; + if (t > v1) /* q1 = B-1 is too large, we need q1 = B-2, which is ok + since u*B - q1*v >= v1*B^2-(B-2)*(v1*B+B-1) = + -B^2 + 2*B*v1 + 3*B - 2 >= 0 since v1>=B/2 and B>=2 */ + { + q1 --; + /* add v to r2:r1 */ + r1 += v0; + r2 += v1 + (r1 < v0); + } + } + else + { + /* divide r3:r2 by v1: requires r3 < v1 */ + udiv_qrnnd (q1, r2, r3, r2, v1); + /* u-extra*v = q1 * v1 + r2 */ + + /* now subtract q1*v0 to r2:0 */ + umul_ppmm (h, l, q1, v0); + t = r2; /* save old value of r2 */ + r1 = -l; + r2 -= h + (l != 0); + /* Note: h + (l != 0) < 2^GMP_NUMB_BITS. */ + + /* we have r2:r1 = oldr2:0 - q1*v0 mod 2^(2*GMP_NUMB_BITS) + thus (u-extra*v)*B = q1 * v + r2:r1 mod 2^(2*GMP_NUMB_BITS) */ + + /* this while loop should be run at most twice */ + while (r2 > t) /* borrow when subtracting h + (l != 0), q1 too large */ + { + q1 --; + /* add v1:v0 to r2:r1 */ + t = r2; + r1 += v0; + r2 += v1 + (r1 < v0); + /* note: since 2^(GMP_NUMB_BITS-1) <= v1 + (r1 < v0) + <= 2^GMP_NUMB_BITS, it suffices to check if r2 <= t to see + if there was a carry or not. */ + } + } + + /* now (u-extra*v)*B = q1 * v + r2:r1 with 0 <= r2:r1 < v */ + + MPFR_ASSERTD(r2 < v1 || (r2 == v1 && r1 < v0)); + + if (MPFR_UNLIKELY(r2 == v1)) + { + q0 = MPFR_LIMB_MAX; + /* r2:r1:0 - q0*(v1:v0) = v1:r1:0 - (B-1)*(v1:v0) + = r1:0 - v0:0 + v1:v0 */ + r0 = v0; + t = v0 - r1; /* t > 0 since r2:r1 < v1:v0 */ + r1 = v1 - t; + if (t > v1) + { + q0 --; + /* add v to r1:r0 */ + r0 += v0; + r1 += v1 + (r0 < v0); + } + } + else + { + /* divide r2:r1 by v1: requires r2 < v1 */ + udiv_qrnnd (q0, r1, r2, r1, v1); + + /* r2:r1 = q0*v1 + r1 */ + + /* subtract q0*v0 to r1:0 */ + umul_ppmm (h, l, q0, v0); + t = r1; + r0 = -l; + r1 -= h + (l != 0); + + /* this while loop should be run at most twice */ + while (r1 > t) /* borrow when subtracting h + (l != 0), + q0 was too large */ + { + q0 --; + /* add v1:v0 to r1:r0 */ + t = r1; + r0 += v0; + r1 += v1 + (r0 < v0); + /* note: since 2^(GMP_NUMB_BITS-1) <= v1 + (r0 < v0) + <= 2^GMP_NUMB_BITS, it suffices to check if r1 <= t to see + if there was a carry or not. */ + } + } + + MPFR_ASSERTD(r1 < v1 || (r1 == v1 && r0 < v0)); + + /* now (u-extra*v)*B^2 = (q1:q0) * v + r1:r0 */ + + sb = r1 | r0; + + /* here, q1:q0 should be an approximation of the quotient (or the exact + quotient), and sb the sticky bit */ + +#if GMP_NUMB_BITS == 64 + round_div2: #endif + if (extra) + { + qx ++; + sb |= q0 & 1; + q0 = (q1 << (GMP_NUMB_BITS - 1)) | (q0 >> 1); + q1 = MPFR_LIMB_HIGHBIT | (q1 >> 1); + } + rb = q0 & (MPFR_LIMB_ONE << (sh - 1)); + sb |= (q0 & mask) ^ rb; + qp[1] = q1; + qp[0] = q0 & ~mask; + + MPFR_SIGN(q) = MPFR_MULT_SIGN (MPFR_SIGN (u), MPFR_SIGN (v)); + + /* rounding */ + if (qx > __gmpfr_emax) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + + /* Warning: underflow should be checked *after* rounding, thus when rounding + away and when q > 0.111...111*2^(emin-1), or when rounding to nearest and + q >= 0.111...111[1]*2^(emin-1), there is no underflow. */ + if (qx < __gmpfr_emin) + { + /* Note: the case 0.111...111*2^(emin-1) < q < 2^(emin-1) is not possible + here since (up to exponent) this would imply 1 - 2^(-p) < u/v < 1, + thus v - 2^(-p)*v < u < v, and since we can assume 1/2 <= v < 1, it + would imply v - 2^(-p) = v - ulp(v) < u < v, which has no solution. */ + + /* For RNDN, mpfr_underflow always rounds away, thus for |q|<=2^(emin-2) + we have to change to RNDZ. This corresponds to: + (a) either qx < emin - 1 + (b) or qx = emin - 1 and qp[1] = 100....000, qp[0] = 0 and rb = sb = 0. + Note: in case (b), it suffices to check whether sb = 0, since rb = 1 + and sb = 0 is not possible (the exact quotient would have p+1 bits, thus + u would need at least p+1 bits). */ + if (rnd_mode == MPFR_RNDN && + (qx < __gmpfr_emin - 1 || + (qp[1] == MPFR_LIMB_HIGHBIT && qp[0] == MPFR_LIMB_ZERO && sb == 0))) + rnd_mode = MPFR_RNDZ; + return mpfr_underflow (q, rnd_mode, MPFR_SIGN(q)); + } + + MPFR_EXP (q) = qx; /* Don't use MPFR_SET_EXP since qx might be < __gmpfr_emin + in the cases "goto rounding" above. */ + if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) + { + MPFR_ASSERTD(qx >= __gmpfr_emin); + MPFR_RET (0); + } + else if (rnd_mode == MPFR_RNDN) + { + /* See the comment in mpfr_div_1. */ + MPFR_ASSERTD(sb != 0); + if (rb == 0) + goto truncate; + else + goto add_one_ulp; + } + else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(q))) + { + truncate: + MPFR_ASSERTD(qx >= __gmpfr_emin); + MPFR_RET(-MPFR_SIGN(q)); + } + else /* round away from zero */ + { + add_one_ulp: + qp[0] += MPFR_LIMB_ONE << sh; + qp[1] += (qp[0] == 0); + if (qp[1] == 0) + { + qp[1] = MPFR_LIMB_HIGHBIT; + if (MPFR_UNLIKELY(qx + 1 > __gmpfr_emax)) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + MPFR_ASSERTD(qx + 1 <= __gmpfr_emax); + MPFR_ASSERTD(qx + 1 >= __gmpfr_emin); + MPFR_SET_EXP (q, qx + 1); + } + MPFR_RET(MPFR_SIGN(q)); + } +} + +#endif /* !defined(MPFR_GENERIC_ABI) */ /* check if {ap, an} is zero */ static int mpfr_mpn_cmpzero (mpfr_limb_ptr ap, mp_size_t an) { + MPFR_ASSERTD (an >= 0); while (an > 0) if (MPFR_LIKELY(ap[--an] != MPFR_LIMB_ZERO)) return 1; @@ -66,6 +675,10 @@ mpfr_mpn_cmp_aux (mpfr_limb_ptr ap, mp_size_t an, mp_size_t k; mp_limb_t bb; + MPFR_ASSERTD (an >= 0); + MPFR_ASSERTD (bn >= 0); + MPFR_ASSERTD (extra == 0 || extra == 1); + if (an >= bn) { k = an - bn; @@ -122,6 +735,8 @@ mpfr_mpn_sub_aux (mpfr_limb_ptr ap, mpfr_limb_ptr bp, mp_size_t n, mp_limb_t bb, rp; MPFR_ASSERTD (cy <= 1); + MPFR_ASSERTD (n >= 0); + while (n--) { bb = (extra) ? ((bp[1] << (GMP_NUMB_BITS-1)) | (bp[0] >> 1)) : bp[0]; @@ -136,27 +751,130 @@ mpfr_mpn_sub_aux (mpfr_limb_ptr ap, mpfr_limb_ptr bp, mp_size_t n, return cy; } -int +/* For large precision, mpz_tdiv_q (which computes only quotient) + is faster than mpn_divrem (which computes also the remainder). + Unfortunately as of GMP 6.0.0 the corresponding mpn_div_q function + is not in the public interface, thus we call mpz_tdiv_q. + + If this function succeeds in computing the correct rounding, return 1, + and put the ternary value in inex. + + Otherwise return 0 (and inex is undefined). +*/ +static int +mpfr_div_with_mpz_tdiv_q (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, + mpfr_rnd_t rnd_mode, int *inex) +{ + mpz_t qm, um, vm; + mpfr_exp_t ue, ve; + mpfr_prec_t qp = MPFR_PREC(q), wp = qp + GMP_NUMB_BITS; + mp_size_t up, vp, k; + int ok; + + mpz_init (qm); + mpz_init (um); + mpz_init (vm); + + ue = mpfr_get_z_2exp (um, u); /* u = um * 2^ue */ + ve = mpfr_get_z_2exp (vm, v); /* v = vm * 2^ve */ + + vp = mpz_sizeinbase (vm, 2); + if (vp > wp) + { + k = vp - wp; /* truncate k bits of vm */ + mpz_tdiv_q_2exp (vm, vm, k); + ve += k; + vp -= k; + } + + /* we want about qp + GMP_NUMB_BITS bits of the quotient, thus um should + have qp + GMP_NUMB_BITS more bits than vm */ + + up = mpz_sizeinbase (um, 2); + if (up > vp + wp) + { + k = up - (vp + wp); /* truncate k bits of um */ + mpz_tdiv_q_2exp (um, um, k); + ue += k; + up -= k; + } + else if (up < vp + wp) /* we need more bits */ + { + k = (vp + wp) - up; + mpz_mul_2exp (um, um, k); + ue -= k; + up += k; + } + + /* now um has exactly wp more bits than vp */ + mpz_tdiv_q (qm, um, vm); + /* qm has either wp or wp+1 bits, and we have: + (a) um = u/2^ue*(1-tu) with tu=0 if no truncation of um, + and 0 <= tu < 2^(1-wp) otherwise; + (b) vm = v/2^ve*(1-tv) with tv=0 if no truncation of vm, + and 0 <= tv < 2^(1-wp) otherwise; + (c) um/vm - 1 < qm <= um/vm, thus qm = um/vm*(1-tq) with + 0 <= tw < 2^(1-wp) since um/vm >= 2^(wp-1) + Altogether we have: + q = u/v*2^(ve-ue)*(1-tu)/(1-tv)*(1-tq) + Thus: + u/v*2^(ve-ue)*(1-2^(2-wp)) < q < u/v*2^(ve-ue)*(1+2^(2-wp)). + If q has wp bits, the error is less than 2^(wp-1)*2^(2-wp) <= 2. + If q has wp+1 bits, the error is less than 2^wp*2^(2-wp) <= 4. + */ + + k = mpz_sizeinbase (qm, 2) - wp; /* 0 or 1 */ + /* Assume qm has wp bits (i.e. k=0) and a directed rounding: if the first + set bit after position 1 has position less than GMP_NUMB_BITS, then + subtracting 2 to qm will not change the bits beyond the GMP_NUMB_BITS + low ones, thus we get correct rounding. + For k=1, we need to start at position 2, and the first set bit has to be + in posiiton less than GMP_NUMB_BITS+1. + For rounding to nearest, the first set bit has to be in position less + than GMP_NUMB_BITS-1 for k=0 (or less than GMP_NUMB_BITS for k=1). + */ + if (mpz_scan1 (qm, k + 1) < GMP_NUMB_BITS + k - (rnd_mode == MPFR_RNDN) && + mpz_scan0 (qm, k + 1) < GMP_NUMB_BITS + k - (rnd_mode == MPFR_RNDN)) + { + MPFR_SAVE_EXPO_DECL (expo); + ok = 1; + MPFR_SAVE_EXPO_MARK (expo); + *inex = mpfr_set_z (q, qm, rnd_mode); + MPFR_SAVE_EXPO_FREE (expo); + /* if we got an underflow or overflow, the result is not valid */ + if (MPFR_IS_SINGULAR(q) || MPFR_EXP(q) == MPFR_EXT_EMIN || + MPFR_EXP(q) == MPFR_EXT_EMAX) + ok = 0; + MPFR_EXP(q) += ue - ve; + *inex = mpfr_check_range (q, *inex, rnd_mode); + } + else + ok = 0; + + mpz_clear (qm); + mpz_clear (um); + mpz_clear (vm); + + return ok; +} + +MPFR_HOT_FUNCTION_ATTR int mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) { - mp_size_t q0size = MPFR_LIMB_SIZE(q); /* number of limbs of destination */ - mp_size_t usize = MPFR_LIMB_SIZE(u); - mp_size_t vsize = MPFR_LIMB_SIZE(v); + mp_size_t q0size, usize, vsize; mp_size_t qsize; /* number of limbs wanted for the computed quotient */ mp_size_t qqsize; mp_size_t k; - mpfr_limb_ptr q0p = MPFR_MANT(q), qp; - mpfr_limb_ptr up = MPFR_MANT(u); - mpfr_limb_ptr vp = MPFR_MANT(v); + mpfr_limb_ptr q0p, qp; + mpfr_limb_ptr up, vp; mpfr_limb_ptr ap; mpfr_limb_ptr bp; mp_limb_t qh; - mp_limb_t sticky_u = MPFR_LIMB_ZERO; + mp_limb_t sticky_u, sticky_v; mp_limb_t low_u; - mp_limb_t sticky_v = MPFR_LIMB_ZERO; mp_limb_t sticky; mp_limb_t sticky3; - mp_limb_t round_bit = MPFR_LIMB_ZERO; + mp_limb_t round_bit; mpfr_exp_t qexp; int sign_quotient; int extra_bit; @@ -215,7 +933,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) { MPFR_ASSERTD (! MPFR_IS_INF (u)); MPFR_SET_INF(q); - mpfr_set_divby0 (); + MPFR_SET_DIVBY0 (); MPFR_RET(0); } } @@ -227,12 +945,72 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) } } + /* When MPFR_GENERIC_ABI is defined, we don't use special code. */ +#if !defined(MPFR_GENERIC_ABI) + if (MPFR_GET_PREC(u) == MPFR_GET_PREC(q) && + MPFR_GET_PREC(v) == MPFR_GET_PREC(q)) + { + if (MPFR_GET_PREC(q) < GMP_NUMB_BITS) + return mpfr_div_1 (q, u, v, rnd_mode); + + if (GMP_NUMB_BITS < MPFR_GET_PREC(q) && + MPFR_GET_PREC(q) < 2 * GMP_NUMB_BITS) + return mpfr_div_2 (q, u, v, rnd_mode); + + if (MPFR_GET_PREC(q) == GMP_NUMB_BITS) + return mpfr_div_1n (q, u, v, rnd_mode); + } +#endif /* !defined(MPFR_GENERIC_ABI) */ + + usize = MPFR_LIMB_SIZE(u); + vsize = MPFR_LIMB_SIZE(v); + q0size = MPFR_LIMB_SIZE(q); /* number of limbs of destination */ + q0p = MPFR_MANT(q); + up = MPFR_MANT(u); + vp = MPFR_MANT(v); + sticky_u = MPFR_LIMB_ZERO; + sticky_v = MPFR_LIMB_ZERO; + round_bit = MPFR_LIMB_ZERO; + /************************************************************************** * * * End of the part concerning special values. * * * **************************************************************************/ + /* 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 + && vp[0] <= ULONG_MAX) + { + mpfr_exp_t exp_v = MPFR_EXP(v); /* save it in case q=v */ + if (MPFR_IS_POS (v)) + 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; + /* The following test is needed, otherwise the next addition + on the exponent may overflow, e.g. when dividing the + largest finite MPFR number by the smallest positive one. */ + if (MPFR_UNLIKELY (MPFR_EXP(q) > __gmpfr_emax - GMP_NUMB_BITS)) + return mpfr_overflow (q, rnd_mode, MPFR_SIGN(q)); + MPFR_EXP(q) += GMP_NUMB_BITS; + return mpfr_check_range (q, inex, rnd_mode); + } + + /* for large precisions, try using truncated division first */ + if (q0size >= 32 && mpfr_div_with_mpz_tdiv_q (q, u, v, rnd_mode, &inex)) + return inex; + MPFR_TMP_MARK(marker); /* set sign */ @@ -252,19 +1030,14 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) l = vsize - 1; while (k != 0 && l != 0 && up[--k] == vp[--l]); /* now k=0 or l=0 or up[k] != vp[l] */ - if (up[k] > vp[l]) - extra_bit = 1; - else if (up[k] < vp[l]) - extra_bit = 0; + if (up[k] != vp[l]) + extra_bit = (up[k] > vp[l]); /* now up[k] = vp[l], thus either k=0 or l=0 */ else if (l == 0) /* no more divisor limb */ extra_bit = 1; else /* k=0: no more dividend limb */ extra_bit = mpfr_mpn_cmpzero (vp, l) == 0; } -#ifdef DEBUG - printf ("extra_bit=%d\n", extra_bit); -#endif /* set exponent */ qexp = MPFR_GET_EXP (u) - MPFR_GET_EXP (v) + extra_bit; @@ -338,6 +1111,9 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) round_bit = (qp[1] >> (sh - 1)) & 1; else round_bit = qp[0] >> (GMP_NUMB_BITS - 1); + /* TODO: add value coverage tests in tdiv to check that + we reach this part with different values of qh and + round_bit (4 cases). */ if (round_bit == 0) { inex = -1; @@ -346,10 +1122,9 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) else /* round_bit = 1 */ goto add_one_ulp; } - else if (like_rndz == 0) /* round away */ + else if (! like_rndz) /* round away */ goto add_one_ulp; - /* else round to zero: nothing to do */ - else + else /* round to zero: nothing to do */ { inex = -1; goto truncate; @@ -433,10 +1208,6 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) /* if Mulders' short division failed, we revert to division with remainder */ qh = mpn_divrem (qp, 0, ap + k, qqsize - k, bp, qsize - k); /* warning: qh may be 1 if u1 == v1, but u < v */ -#ifdef DEBUG2 - printf ("q="); mpfr_mpn_print (qp, qsize); - printf ("r="); mpfr_mpn_print (ap, qsize); -#endif k = qsize; sticky_u = sticky_u || mpfr_mpn_cmpzero (ap, k); @@ -464,20 +1235,12 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) including the round bit, and 1 <= sh2 <= GMP_NUMB_BITS is the number of bits in sticky3 */ inex = (sticky != MPFR_LIMB_ZERO) || (sticky3 != MPFR_LIMB_ZERO); -#ifdef DEBUG - printf ("sticky=%lu sticky3=%lu inex=%d\n", - (unsigned long) sticky, (unsigned long) sticky3, inex); -#endif /* to round, we distinguish two cases: (a) vsize <= qsize: we used the full divisor (b) vsize > qsize: the divisor was truncated */ -#ifdef DEBUG - printf ("vsize=%lu qsize=%lu\n", - (unsigned long) vsize, (unsigned long) qsize); -#endif if (MPFR_LIKELY(vsize <= qsize)) /* use the full divisor */ { if (MPFR_LIKELY(rnd_mode == MPFR_RNDN)) @@ -510,10 +1273,6 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) { round_bit = sticky3 & (MPFR_LIMB_ONE << (sh2 - 1)); sticky3 = sticky3 ^ round_bit; -#ifdef DEBUG - printf ("rb=%lu sb=%lu\n", - (unsigned long) round_bit, (unsigned long) sticky3); -#endif } if (sticky3 != MPFR_LIMB_ZERO && sticky3 != MPFR_LIMB_ONE) { @@ -540,7 +1299,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) if (qh) qh2 = mpn_add_n (sp + qsize, sp + qsize, vp, k); else - qh2 = (mp_limb_t) 0; + qh2 = MPFR_LIMB_ZERO; qp[0] ^= sticky3orig; /* restore truncated quotient */ /* compare qh2 + {sp, k + qsize} to {ap, qsize} + low(u) */ @@ -551,9 +1310,6 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) mpfr_mpn_cmp_aux (sp, k, up, usize - qqsize, extra_bit) : mpfr_mpn_cmpzero (sp, k); } -#ifdef DEBUG - printf ("cmp(q*v0,r+u0)=%d\n", cmp_s_r); -#endif /* now cmp_s_r > 0 if {sp, vsize} > {ap, qsize} + low(u) cmp_s_r = 0 if {sp, vsize} = {ap, qsize} + low(u) cmp_s_r < 0 if {sp, vsize} < {ap, qsize} + low(u) */ @@ -611,9 +1367,6 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) cmp_s_r = 1; /* since in fact we subtracted less than 1 */ } -#ifdef DEBUG - printf ("cmp(q*v0-(r+u0),v)=%d\n", cmp_s_r); -#endif if (cmp_s_r <= 0) /* q1-1 <= u/v < q1 */ { if (sticky3 == MPFR_LIMB_ONE) |