summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/div.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/div.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/div.c215
1 files changed, 60 insertions, 155 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/div.c b/Build/source/libs/mpfr/mpfr-src/src/div.c
index 012a5d7ef28..a04f2539d91 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/div.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/div.c
@@ -1,6 +1,6 @@
/* mpfr_div -- divide two floating-point numbers
-Copyright 1999, 2001-2019 Free Software Foundation, Inc.
+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.
@@ -233,15 +233,20 @@ mpfr_div_1 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
{
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_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));
}
}
@@ -376,15 +381,9 @@ mpfr_div_1n (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
{
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);
- }
+ /* 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));
}
}
@@ -636,16 +635,10 @@ mpfr_div_2 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
{
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);
- }
+ 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));
}
}
@@ -739,9 +732,9 @@ mpfr_mpn_sub_aux (mpfr_limb_ptr ap, mpfr_limb_ptr bp, mp_size_t n,
while (n--)
{
- bb = (extra) ? ((bp[1] << (GMP_NUMB_BITS-1)) | (bp[0] >> 1)) : bp[0];
+ 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_ZERO) ?
+ cy = (ap[0] < bb) || (cy && rp == MPFR_LIMB_MAX) ?
MPFR_LIMB_ONE : MPFR_LIMB_ZERO;
ap[0] = rp;
ap ++;
@@ -751,113 +744,6 @@ mpfr_mpn_sub_aux (mpfr_limb_ptr ap, mpfr_limb_ptr bp, mp_size_t n,
return cy;
}
-/* 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)
{
@@ -978,12 +864,14 @@ 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.
+ /* 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)
@@ -1006,10 +894,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
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;
+#endif
MPFR_TMP_MARK(marker);
@@ -1082,6 +967,7 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
}
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},
@@ -1205,8 +1091,27 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
* *
**************************************************************************/
+ /* 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;
@@ -1279,9 +1184,9 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
sticky = sticky3;
goto case_1;
}
- else /* hard case: we have to compare q1 * v0 and r + low(u),
+ else /* hard case: we have to compare q1 * v0 and r + u0,
where q1 * v0 has qsize + (vsize-qsize) = vsize limbs, and
- r + low(u) has qsize + (usize-2*qsize) = usize-qsize limbs */
+ r + u0 has qsize + (usize-2*qsize) = usize-qsize limbs */
{
mp_size_t l;
mpfr_limb_ptr sp;
@@ -1302,17 +1207,17 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
qh2 = MPFR_LIMB_ZERO;
qp[0] ^= sticky3orig; /* restore truncated quotient */
- /* compare qh2 + {sp, k + qsize} to {ap, qsize} + low(u) */
+ /* 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 low(u) */
+ 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} + low(u)
- cmp_s_r = 0 if {sp, vsize} = {ap, qsize} + low(u)
- cmp_s_r < 0 if {sp, vsize} < {ap, qsize} + low(u) */
+ /* 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;
@@ -1320,11 +1225,11 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
}
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 u0 from q*v0 */
+ the low part u0 of the dividend from q*v0 */
{
mp_limb_t cy = MPFR_LIMB_ZERO;
- /* subtract low(u)>>extra_bit if non-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;
@@ -1333,11 +1238,11 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
if (low_u != MPFR_LIMB_ZERO)
{
mp_size_t m;
- l = usize - qqsize; /* number of low limbs in u */
+ 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 more limbs than s:
+ 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} */
{