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 | 1446 |
1 files changed, 0 insertions, 1446 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/div.c b/Build/source/libs/mpfr/mpfr-src/src/div.c deleted file mode 100644 index a04f2539d91..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/div.c +++ /dev/null @@ -1,1446 +0,0 @@ -/* mpfr_div -- divide two floating-point numbers - -Copyright 1999, 2001-2020 Free Software Foundation, Inc. -Contributed by the AriC and Caramba projects, INRIA. - -This file is part of the GNU MPFR Library. - -The GNU MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Lesser General Public License as published by -the Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -The GNU MPFR Library is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public -License for more details. - -You should have received a copy of the GNU Lesser General Public License -along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see -https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., -51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ - -/* References: - [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" - -#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_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) -{ - 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; - MPFR_ASSERTD(qp[0] != 0); - /* It is not possible to have an overflow in the addition above. - Proof: if p is the precision of the inputs, it would mean we have two - integers n and d with 2^(p-1) <= n, d < 2^p, such that the binary - expansion of n/d starts with p '1', and has at least one '1' later. - We distinguish two cases: - (1) if n/d < 1, it would mean 1-2^(-p) < n/d < 1 - (2) if n/d >= 1, it would mean 2-2^(1-p) < n/d < 1 - In case (1), multiplying by d we get 1-d/2^p < n < d, - which has no integer solution since d/2^p < 1. - In case (2), multiplying by d we get 2d-2d/2^p < n < 2d: - (2a) if d=2^(p-1), we get 2^p-1 < n < 2^p which has no solution; - if d>=2^(p-1)+1, then 2d-2d/2^p >= 2^p+2-2 = 2^p, thus there is - solution n < 2^p either. */ - 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; - /* there can be no overflow in the addition above, - see the analysis of mpfr_div_1 */ - MPFR_ASSERTD(qp[0] != 0); - 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; - /* there can be no overflow in the addition above, - see the analysis of mpfr_div_1 */ - MPFR_ASSERTD(qp[1] != 0); - 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; - return 0; -} - -/* compare {ap, an} and {bp, bn} >> extra, - aligned by the more significant limbs. - Takes into account bp[0] for extra=1. -*/ -static int -mpfr_mpn_cmp_aux (mpfr_limb_ptr ap, mp_size_t an, - mpfr_limb_ptr bp, mp_size_t bn, int extra) -{ - int cmp = 0; - 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; - while (cmp == 0 && bn > 0) - { - bn --; - bb = (extra) ? ((bp[bn+1] << (GMP_NUMB_BITS - 1)) | (bp[bn] >> 1)) - : bp[bn]; - cmp = (ap[k + bn] > bb) ? 1 : ((ap[k + bn] < bb) ? -1 : 0); - } - bb = (extra) ? bp[0] << (GMP_NUMB_BITS - 1) : MPFR_LIMB_ZERO; - while (cmp == 0 && k > 0) - { - k--; - cmp = (ap[k] > bb) ? 1 : ((ap[k] < bb) ? -1 : 0); - bb = MPFR_LIMB_ZERO; /* ensure we consider only once bp[0] & 1 */ - } - if (cmp == 0 && bb != MPFR_LIMB_ZERO) - cmp = -1; - } - else /* an < bn */ - { - k = bn - an; - while (cmp == 0 && an > 0) - { - an --; - bb = (extra) ? ((bp[k+an+1] << (GMP_NUMB_BITS - 1)) | (bp[k+an] >> 1)) - : bp[k+an]; - if (ap[an] > bb) - cmp = 1; - else if (ap[an] < bb) - cmp = -1; - } - while (cmp == 0 && k > 0) - { - k--; - bb = (extra) ? ((bp[k+1] << (GMP_NUMB_BITS - 1)) | (bp[k] >> 1)) - : bp[k]; - cmp = (bb != MPFR_LIMB_ZERO) ? -1 : 0; - } - if (cmp == 0 && extra && (bp[0] & MPFR_LIMB_ONE)) - cmp = -1; - } - return cmp; -} - -/* {ap, n} <- {ap, n} - {bp, n} >> extra - cy, with cy = 0 or 1. - Return borrow out. -*/ -static mp_limb_t -mpfr_mpn_sub_aux (mpfr_limb_ptr ap, mpfr_limb_ptr bp, mp_size_t n, - mp_limb_t cy, int extra) -{ - mp_limb_t bb, rp; - - MPFR_ASSERTD (cy <= 1); - MPFR_ASSERTD (n >= 0); - - while (n--) - { - bb = (extra) ? (MPFR_LIMB_LSHIFT(bp[1],GMP_NUMB_BITS-1) | (bp[0] >> 1)) : bp[0]; - rp = ap[0] - bb - cy; - cy = (ap[0] < bb) || (cy && rp == MPFR_LIMB_MAX) ? - MPFR_LIMB_ONE : MPFR_LIMB_ZERO; - ap[0] = rp; - ap ++; - bp ++; - } - MPFR_ASSERTD (cy <= 1); - return cy; -} - -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, 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, qp; - mpfr_limb_ptr up, vp; - mpfr_limb_ptr ap; - mpfr_limb_ptr bp; - mp_limb_t qh; - mp_limb_t sticky_u, sticky_v; - mp_limb_t low_u; - mp_limb_t sticky; - mp_limb_t sticky3; - mp_limb_t round_bit; - mpfr_exp_t qexp; - int sign_quotient; - int extra_bit; - int sh, sh2; - int inex; - int like_rndz; - MPFR_TMP_DECL(marker); - - MPFR_LOG_FUNC ( - ("u[%Pu]=%.*Rg v[%Pu]=%.*Rg rnd=%d", - mpfr_get_prec(u), mpfr_log_prec, u, - mpfr_get_prec (v),mpfr_log_prec, v, rnd_mode), - ("q[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(q), mpfr_log_prec, q, inex)); - - /************************************************************************** - * * - * This part of the code deals with special cases * - * * - **************************************************************************/ - - if (MPFR_UNLIKELY(MPFR_ARE_SINGULAR(u,v))) - { - if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) - { - MPFR_SET_NAN(q); - MPFR_RET_NAN; - } - sign_quotient = MPFR_MULT_SIGN( MPFR_SIGN(u) , MPFR_SIGN(v) ); - MPFR_SET_SIGN(q, sign_quotient); - if (MPFR_IS_INF(u)) - { - if (MPFR_IS_INF(v)) - { - MPFR_SET_NAN(q); - MPFR_RET_NAN; - } - else - { - MPFR_SET_INF(q); - MPFR_RET(0); - } - } - else if (MPFR_IS_INF(v)) - { - MPFR_SET_ZERO (q); - MPFR_RET (0); - } - else if (MPFR_IS_ZERO (v)) - { - if (MPFR_IS_ZERO (u)) - { - MPFR_SET_NAN(q); - MPFR_RET_NAN; - } - else - { - MPFR_ASSERTD (! MPFR_IS_INF (u)); - MPFR_SET_INF(q); - MPFR_SET_DIVBY0 (); - MPFR_RET(0); - } - } - else - { - MPFR_ASSERTD (MPFR_IS_ZERO (u)); - MPFR_SET_ZERO (q); - MPFR_RET (0); - } - } - - /* 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 and MPFR_LONG_WITHIN_LIMB is defined, - 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. */ -#ifdef MPFR_LONG_WITHIN_LIMB - 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); - } -#endif - - MPFR_TMP_MARK(marker); - - /* set sign */ - sign_quotient = MPFR_MULT_SIGN( MPFR_SIGN(u) , MPFR_SIGN(v) ); - MPFR_SET_SIGN(q, sign_quotient); - - /* determine if an extra bit comes from the division, i.e. if the - significand of u (as a fraction in [1/2, 1[) is larger than that - of v */ - if (MPFR_LIKELY(up[usize - 1] != vp[vsize - 1])) - extra_bit = (up[usize - 1] > vp[vsize - 1]) ? 1 : 0; - else /* most significant limbs are equal, must look at further limbs */ - { - mp_size_t l; - - k = usize - 1; - 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 = (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; - } - - /* set exponent */ - qexp = MPFR_GET_EXP (u) - MPFR_GET_EXP (v) + extra_bit; - - /* sh is the number of zero bits in the low limb of the quotient */ - MPFR_UNSIGNED_MINUS_MODULO(sh, MPFR_PREC(q)); - - like_rndz = rnd_mode == MPFR_RNDZ || - rnd_mode == (sign_quotient < 0 ? MPFR_RNDU : MPFR_RNDD); - - /************************************************************************** - * * - * We first try Mulders' short division (for large operands) * - * * - **************************************************************************/ - - if (MPFR_UNLIKELY(q0size >= MPFR_DIV_THRESHOLD && - vsize >= MPFR_DIV_THRESHOLD)) - { - mp_size_t n = q0size + 1; /* we will perform a short (2n)/n division */ - mpfr_limb_ptr ap, bp, qp; - mpfr_prec_t p; - - /* since Mulders' short division clobbers the dividend, we have to - copy it */ - ap = MPFR_TMP_LIMBS_ALLOC (n + n); - if (usize >= n + n) /* truncate the dividend */ - MPN_COPY(ap, up + usize - (n + n), n + n); - else /* zero-pad the dividend */ - { - MPN_COPY(ap + (n + n) - usize, up, usize); - MPN_ZERO(ap, (n + n) - usize); - } - - if (vsize >= n) /* truncate the divisor */ - bp = vp + vsize - n; - else /* zero-pad the divisor */ - { - bp = MPFR_TMP_LIMBS_ALLOC (n); - MPN_COPY(bp + n - vsize, vp, vsize); - MPN_ZERO(bp, n - vsize); - } - - qp = MPFR_TMP_LIMBS_ALLOC (n); - /* since n = q0size + 1, we have n >= 2 here */ - qh = mpfr_divhigh_n (qp, ap, bp, n); - MPFR_ASSERTD (qh == 0 || qh == 1); - /* in all cases, the error is at most (2n+2) ulps on qh*B^n+{qp,n}, - cf algorithms.tex */ - - p = n * GMP_NUMB_BITS - MPFR_INT_CEIL_LOG2 (2 * n + 2); - /* If rnd=RNDN, we need to be able to round with a directed rounding - and one more bit. */ - if (qh == 1) - { - mpn_rshift (qp, qp, n, 1); - qp[n - 1] |= MPFR_LIMB_HIGHBIT; - } - if (MPFR_LIKELY (mpfr_round_p (qp, n, p, - MPFR_PREC(q) + (rnd_mode == MPFR_RNDN)))) - { - /* we can round correctly whatever the rounding mode */ - MPN_COPY (q0p, qp + 1, q0size); - q0p[0] &= ~MPFR_LIMB_MASK(sh); /* put to zero low sh bits */ - - if (rnd_mode == MPFR_RNDN) /* round to nearest */ - { - /* we know we can round, thus we are never in the even rule case: - if the round bit is 0, we truncate - if the round bit is 1, we add 1 */ - if (sh > 0) - 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; - goto truncate; - } - else /* round_bit = 1 */ - goto add_one_ulp; - } - else if (! like_rndz) /* round away */ - goto add_one_ulp; - else /* round to zero: nothing to do */ - { - inex = -1; - goto truncate; - } - } - } - - /************************************************************************** - * * - * Mulders' short division failed: we revert to integer division * - * * - **************************************************************************/ - - if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDN && sh == 0)) - { /* we compute the quotient with one more limb, in order to get - the round bit in the quotient, and the remainder only contains - sticky bits */ - qsize = q0size + 1; - /* need to allocate memory for the quotient */ - qp = MPFR_TMP_LIMBS_ALLOC (qsize); - } - else - { - qsize = q0size; - qp = q0p; /* directly put the quotient in the destination */ - } - qqsize = qsize + qsize; - - /* prepare the dividend */ - ap = MPFR_TMP_LIMBS_ALLOC (qqsize); - if (MPFR_LIKELY(qqsize > usize)) /* use the full dividend */ - { - k = qqsize - usize; /* k > 0 */ - MPN_ZERO(ap, k); - if (extra_bit) - ap[k - 1] = mpn_rshift (ap + k, up, usize, 1); - else - MPN_COPY(ap + k, up, usize); - } - else /* truncate the dividend */ - { - k = usize - qqsize; - if (extra_bit) - sticky_u = mpn_rshift (ap, up + k, qqsize, 1); - else - MPN_COPY(ap, up + k, qqsize); - sticky_u = sticky_u || mpfr_mpn_cmpzero (up, k); - } - low_u = sticky_u; - - /* now sticky_u is non-zero iff the truncated part of u is non-zero */ - - /* prepare the divisor */ - if (MPFR_LIKELY(vsize >= qsize)) - { - k = vsize - qsize; - if (qp != vp) - bp = vp + k; /* avoid copying the divisor */ - else /* need to copy, since mpn_divrem doesn't allow overlap - between quotient and divisor, necessarily k = 0 - since quotient and divisor are the same mpfr variable */ - { - bp = MPFR_TMP_LIMBS_ALLOC (qsize); - MPN_COPY(bp, vp, vsize); - } - sticky_v = sticky_v || mpfr_mpn_cmpzero (vp, k); - k = 0; - } - else /* vsize < qsize: small divisor case */ - { - bp = vp; - k = qsize - vsize; - } - - /************************************************************************** - * * - * Here we perform the real division of {ap+k,qqsize-k} by {bp,qsize-k} * - * * - **************************************************************************/ - - /* In the general case (usize > 2*qsize and vsize > qsize), we have: - ______________________________________ - | | | u1 has 2*qsize limbs - | u1 | u0 | u0 has usize-2*qsize limbs - |__________________________|___________| - - ____________________ - | | | v1 has qsize limbs - | v1 | v0 | v0 has vsize-qsize limbs - |___________|________| - - We divide u1 by v1, with quotient in qh + {qp, qsize} and - remainder (denoted r below) stored in place of the low qsize limbs of u1. - */ - - /* if Mulders' short division failed, we revert to division with remainder */ - qh = mpn_divrem (qp, 0, ap + k, qqsize - k, bp, qsize - k); - /* let u1 be the upper part of u, and v1 the upper part of v (with sticky_u - and sticky_v representing the lower parts), then the quotient of u1 by v1 - is now in {qp, qsize}, with possible carry in qh, and the remainder in - {ap + k, qsize - k} */ - /* warning: qh may be 1 if u1 == v1, but u < v */ - - k = qsize; - sticky_u = sticky_u || mpfr_mpn_cmpzero (ap, k); - - sticky = sticky_u | sticky_v; - - /* now sticky is non-zero iff one of the following holds: - (a) the truncated part of u is non-zero - (b) the truncated part of v is non-zero - (c) the remainder from division is non-zero */ - - if (MPFR_LIKELY(qsize == q0size)) - { - sticky3 = qp[0] & MPFR_LIMB_MASK(sh); /* does nothing when sh=0 */ - sh2 = sh; - } - else /* qsize = q0size + 1: only happens when rnd_mode=MPFR_RNDN and sh=0 */ - { - MPN_COPY (q0p, qp + 1, q0size); - sticky3 = qp[0]; - sh2 = GMP_NUMB_BITS; - } - qp[0] ^= sticky3; - /* sticky3 contains the truncated bits from the quotient, - 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); - - /* to round, we distinguish two cases: - (a) vsize <= qsize: we used the full divisor - (b) vsize > qsize: the divisor was truncated - */ - - if (MPFR_LIKELY(vsize <= qsize)) /* use the full divisor */ - { - if (MPFR_LIKELY(rnd_mode == MPFR_RNDN)) - { - round_bit = sticky3 & (MPFR_LIMB_ONE << (sh2 - 1)); - sticky = (sticky3 ^ round_bit) | sticky_u; - } - else if (like_rndz || inex == 0) - sticky = (inex == 0) ? MPFR_LIMB_ZERO : MPFR_LIMB_ONE; - else /* round away from zero */ - sticky = MPFR_LIMB_ONE; - goto case_1; - } - else /* vsize > qsize: need to truncate the divisor */ - { - if (inex == 0) - goto truncate; - else - { - /* We know the estimated quotient is an upper bound of the exact - quotient (with rounding toward zero), with a difference of at - most 2 in qp[0]. - Thus we can round except when sticky3 is 000...000 or 000...001 - for directed rounding, and 100...000 or 100...001 for rounding - to nearest. (For rounding to nearest, we cannot determine the - inexact flag for 000...000 or 000...001.) - */ - mp_limb_t sticky3orig = sticky3; - if (rnd_mode == MPFR_RNDN) - { - round_bit = sticky3 & (MPFR_LIMB_ONE << (sh2 - 1)); - sticky3 = sticky3 ^ round_bit; - } - if (sticky3 != MPFR_LIMB_ZERO && sticky3 != MPFR_LIMB_ONE) - { - sticky = sticky3; - goto case_1; - } - else /* hard case: we have to compare q1 * v0 and r + u0, - where q1 * v0 has qsize + (vsize-qsize) = vsize limbs, and - r + u0 has qsize + (usize-2*qsize) = usize-qsize limbs */ - { - mp_size_t l; - mpfr_limb_ptr sp; - int cmp_s_r; - mp_limb_t qh2; - - sp = MPFR_TMP_LIMBS_ALLOC (vsize); - k = vsize - qsize; - /* sp <- {qp, qsize} * {vp, vsize-qsize} */ - qp[0] ^= sticky3orig; /* restore original quotient */ - if (qsize >= k) - mpn_mul (sp, qp, qsize, vp, k); - else - mpn_mul (sp, vp, k, qp, qsize); - if (qh) - qh2 = mpn_add_n (sp + qsize, sp + qsize, vp, k); - else - qh2 = MPFR_LIMB_ZERO; - qp[0] ^= sticky3orig; /* restore truncated quotient */ - - /* compare qh2 + {sp, k + qsize} to {ap, qsize} + u0 */ - cmp_s_r = (qh2 != 0) ? 1 : mpn_cmp (sp + k, ap, qsize); - if (cmp_s_r == 0) /* compare {sp, k} and u0 */ - { - cmp_s_r = (usize >= qqsize) ? - mpfr_mpn_cmp_aux (sp, k, up, usize - qqsize, extra_bit) : - mpfr_mpn_cmpzero (sp, k); - } - /* now cmp_s_r > 0 if {sp, vsize} > {ap, qsize} + u0 - cmp_s_r = 0 if {sp, vsize} = {ap, qsize} + u0 - cmp_s_r < 0 if {sp, vsize} < {ap, qsize} + u0 */ - if (cmp_s_r <= 0) /* quotient is in [q1, q1+1) */ - { - sticky = (cmp_s_r == 0) ? sticky3 : MPFR_LIMB_ONE; - goto case_1; - } - else /* cmp_s_r > 0, quotient is < q1: to determine if it is - in [q1-2,q1-1] or in [q1-1,q1], we need to subtract - the low part u0 of the dividend from q*v0 */ - { - mp_limb_t cy = MPFR_LIMB_ZERO; - - /* subtract u0 >> extra_bit if non-zero */ - if (qh2 != 0) /* whatever the value of {up, m + k}, it - will be smaller than qh2 + {sp, k} */ - cmp_s_r = 1; - else - { - if (low_u != MPFR_LIMB_ZERO) - { - mp_size_t m; - l = usize - qqsize; /* number of limbs in u0 */ - m = (l > k) ? l - k : 0; - cy = (extra_bit) ? - (up[m] & MPFR_LIMB_ONE) : MPFR_LIMB_ZERO; - if (l >= k) /* u0 has at least as many limbs than s: - first look if {up, m} is not zero, - and compare {sp, k} and {up + m, k} */ - { - cy = cy || mpfr_mpn_cmpzero (up, m); - low_u = cy; - cy = mpfr_mpn_sub_aux (sp, up + m, k, - cy, extra_bit); - } - else /* l < k: s has more limbs than u0 */ - { - low_u = MPFR_LIMB_ZERO; - if (cy != MPFR_LIMB_ZERO) - cy = mpn_sub_1 (sp + k - l - 1, sp + k - l - 1, - 1, MPFR_LIMB_HIGHBIT); - cy = mpfr_mpn_sub_aux (sp + k - l, up, l, - cy, extra_bit); - } - } - MPFR_ASSERTD (cy <= 1); - cy = mpn_sub_1 (sp + k, sp + k, qsize, cy); - /* subtract r */ - cy += mpn_sub_n (sp + k, sp + k, ap, qsize); - MPFR_ASSERTD (cy <= 1); - /* now compare {sp, ssize} to v */ - cmp_s_r = mpn_cmp (sp, vp, vsize); - if (cmp_s_r == 0 && low_u != MPFR_LIMB_ZERO) - cmp_s_r = 1; /* since in fact we subtracted - less than 1 */ - } - if (cmp_s_r <= 0) /* q1-1 <= u/v < q1 */ - { - if (sticky3 == MPFR_LIMB_ONE) - { /* q1-1 is either representable (directed rounding), - or the middle of two numbers (nearest) */ - sticky = (cmp_s_r) ? MPFR_LIMB_ONE : MPFR_LIMB_ZERO; - goto case_1; - } - /* now necessarily sticky3=0 */ - else if (round_bit == MPFR_LIMB_ZERO) - { /* round_bit=0, sticky3=0: q1-1 is exact only - when sh=0 */ - inex = (cmp_s_r || sh) ? -1 : 0; - if (rnd_mode == MPFR_RNDN || - (! like_rndz && inex != 0)) - { - inex = 1; - goto truncate_check_qh; - } - else /* round down */ - goto sub_one_ulp; - } - else /* sticky3=0, round_bit=1 ==> rounding to nearest */ - { - inex = cmp_s_r; - goto truncate; - } - } - else /* q1-2 < u/v < q1-1 */ - { - /* if rnd=MPFR_RNDN, the result is q1 when - q1-2 >= q1-2^(sh-1), i.e. sh >= 2, - otherwise (sh=1) it is q1-2 */ - if (rnd_mode == MPFR_RNDN) /* sh > 0 */ - { - /* Case sh=1: sb=0 always, and q1-rb is exactly - representable, like q1-rb-2. - rb action - 0 subtract two ulps, inex=-1 - 1 truncate, inex=1 - - Case sh>1: one ulp is 2^(sh-1) >= 2 - rb sb action - 0 0 truncate, inex=1 - 0 1 truncate, inex=1 - 1 x truncate, inex=-1 - */ - if (sh == 1) - { - if (round_bit == MPFR_LIMB_ZERO) - { - inex = -1; - sh = 0; - goto sub_two_ulp; - } - else - { - inex = 1; - goto truncate_check_qh; - } - } - else /* sh > 1 */ - { - inex = (round_bit == MPFR_LIMB_ZERO) ? 1 : -1; - goto truncate_check_qh; - } - } - else if (like_rndz) - { - /* the result is down(q1-2), i.e. subtract one - ulp if sh > 0, and two ulps if sh=0 */ - inex = -1; - if (sh > 0) - goto sub_one_ulp; - else - goto sub_two_ulp; - } - /* if round away from zero, the result is up(q1-1), - which is q1 unless sh = 0, where it is q1-1 */ - else - { - inex = 1; - if (sh > 0) - goto truncate_check_qh; - else /* sh = 0 */ - goto sub_one_ulp; - } - } - } - } - } - } - - case_1: /* quotient is in [q1, q1+1), - round_bit is the round_bit (0 for directed rounding), - sticky the sticky bit */ - if (like_rndz || (round_bit == MPFR_LIMB_ZERO && sticky == MPFR_LIMB_ZERO)) - { - inex = round_bit == MPFR_LIMB_ZERO && sticky == MPFR_LIMB_ZERO ? 0 : -1; - goto truncate; - } - else if (rnd_mode == MPFR_RNDN) /* sticky <> 0 or round <> 0 */ - { - if (round_bit == MPFR_LIMB_ZERO) /* necessarily sticky <> 0 */ - { - inex = -1; - goto truncate; - } - /* round_bit = 1 */ - else if (sticky != MPFR_LIMB_ZERO) - goto add_one_ulp; /* inex=1 */ - else /* round_bit=1, sticky=0 */ - goto even_rule; - } - else /* round away from zero, sticky <> 0 */ - goto add_one_ulp; /* with inex=1 */ - - sub_two_ulp: - /* we cannot subtract MPFR_LIMB_MPFR_LIMB_ONE << (sh+1) since this is - undefined for sh = GMP_NUMB_BITS */ - qh -= mpn_sub_1 (q0p, q0p, q0size, MPFR_LIMB_ONE << sh); - /* go through */ - - sub_one_ulp: - qh -= mpn_sub_1 (q0p, q0p, q0size, MPFR_LIMB_ONE << sh); - /* go through truncate_check_qh */ - - truncate_check_qh: - if (qh) - { - if (MPFR_LIKELY (qexp < MPFR_EXP_MAX)) - qexp ++; - /* else qexp is now incorrect, but one will still get an overflow */ - q0p[q0size - 1] = MPFR_LIMB_HIGHBIT; - } - goto truncate; - - even_rule: /* has to set inex */ - inex = (q0p[0] & (MPFR_LIMB_ONE << sh)) ? 1 : -1; - if (inex < 0) - goto truncate; - /* else go through add_one_ulp */ - - add_one_ulp: - inex = 1; /* always here */ - if (mpn_add_1 (q0p, q0p, q0size, MPFR_LIMB_ONE << sh)) - { - if (MPFR_LIKELY (qexp < MPFR_EXP_MAX)) - qexp ++; - /* else qexp is now incorrect, but one will still get an overflow */ - q0p[q0size - 1] = MPFR_LIMB_HIGHBIT; - } - - truncate: /* inex already set */ - - MPFR_TMP_FREE(marker); - - /* check for underflow/overflow */ - if (MPFR_UNLIKELY(qexp > __gmpfr_emax)) - return mpfr_overflow (q, rnd_mode, sign_quotient); - else if (MPFR_UNLIKELY(qexp < __gmpfr_emin)) - { - if (rnd_mode == MPFR_RNDN && ((qexp < __gmpfr_emin - 1) || - (inex >= 0 && mpfr_powerof2_raw (q)))) - rnd_mode = MPFR_RNDZ; - return mpfr_underflow (q, rnd_mode, sign_quotient); - } - MPFR_SET_EXP(q, qexp); - - inex *= sign_quotient; - MPFR_RET (inex); -} |