summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/cmp2.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/cmp2.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/cmp2.c254
1 files changed, 170 insertions, 84 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/cmp2.c b/Build/source/libs/mpfr/mpfr-src/src/cmp2.c
index 2526d43e66d..6519a5254cc 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/cmp2.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/cmp2.c
@@ -1,6 +1,6 @@
/* mpfr_cmp2 -- exponent shift when subtracting two numbers.
-Copyright 1999-2004, 2006-2019 Free Software Foundation, Inc.
+Copyright 1999-2004, 2006-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -29,18 +29,24 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
Assumes neither of b or c is NaN, +/- infinity, or +/- 0.
- In other terms, if |b| != |c|, mpfr_cmp2 (b, c) returns
- EXP(max(|b|,|c|)) - EXP(|b| - |c|).
+ In other terms, if |b| != |c|, mpfr_cmp2 (b, c) stores
+ EXP(max(|b|,|c|)) - EXP(|b| - |c|) in *cancel.
+
+ One necessarily has 0 <= cancel <= max(PREC(b),PREC(c)), so that this
+ value is representable in a mpfr_prec_t. Note that in the code, the
+ maximum intermediate value is cancel + 1, but since MPFR_PREC_MAX is
+ not the maximum value of mpfr_prec_t, there is no integer overflow.
*/
int
mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
{
- mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
+ mp_limb_t *bp, *cp, bb, cc, lastc, dif;
+ int high_dif; /* manipulated like a boolean */
mp_size_t bn, cn;
mpfr_exp_t sdiff_exp;
mpfr_uexp_t diff_exp;
- mpfr_prec_t res = 0;
+ mpfr_prec_t res = 0; /* will be the number of canceled bits (*cancel) */
int sign;
/* b=c should not happen, since cmp2 is called only from agm (with
@@ -55,24 +61,32 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
sdiff_exp = MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)) ?
mpfr_ubf_diff_exp (b, c) : MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
- /* The returned result is restricted to [MPFR_EXP_MIN,MPFR_EXP_MAX],
+ /* The returned result is saturated to [MPFR_EXP_MIN,MPFR_EXP_MAX],
which is the range of the mpfr_exp_t type. But under the condition
- below, the value of cancel will not be affected. */
+ below, since |MPFR_EXP_MIN| >= MPFR_EXP_MAX, the value of cancel
+ will not be affected: by symmetry (as done in the code), assume
+ |b| >= |c|; if EXP(b) - EXP(c) >= MPFR_EXP_MAX, then |c| < ulp(b),
+ so that the value of cancel is 0, or 1 if |b| is a power of 2,
+ whatever the exact value of EXP(b) - EXP(c). */
MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX);
if (sdiff_exp >= 0)
{
- sign = 1;
+ sign = 1; /* assumes |b| > |c|; will be changed if not. */
diff_exp = sdiff_exp;
bp = MPFR_MANT(b);
cp = MPFR_MANT(c);
+ /* index of the most significant limb of b and c */
bn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS;
- cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; /* # of limbs of c minus 1 */
+ cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS;
+ /* If diff_exp != 0, i.e. diff_exp > 0, then |b| > |c|. Otherwise... */
if (diff_exp == 0)
{
+ /* Skip the identical most significant limbs, adding GMP_NUMB_BITS
+ to the number of canceled bits at each iteration. */
while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
{
bn--;
@@ -85,11 +99,15 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
if (MPFR_LIKELY (cn < 0)) /* |b| = |c| */
return 0;
- /* b has been read entirely, but not c. Replace b by c for the
- symmetric case below (only the sign differs if not 0). */
+ /* b has been read entirely, but not c. Thus |b| <= |c|.
+ Swap (bp,bn) and (cp,cn), and take the opposite sign
+ for the symmetric case below (simulating a swap).
+ Note: cp will not be used, thus is not assigned; and
+ "cn = -1;" is necessary to enter the following "if"
+ (probably less confusing than a "goto"). */
bp = cp;
bn = cn;
- cn = -1; /* to enter the following "if" */
+ cn = -1;
sign = -1;
}
@@ -100,6 +118,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
MPFR_ASSERTD (bn >= 0);
+ /* Skip null limbs of b (= non-represented null limbs of c),
+ adding GMP_NUMB_BITS to the number of canceled bits at
+ each iteration. */
while (bp[bn] == 0)
{
if (--bn < 0) /* |b| = |c| */
@@ -107,7 +128,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
res += GMP_NUMB_BITS;
}
- count_leading_zeros (z, bp[bn]); /* bp[bn] <> 0 */
+ count_leading_zeros (z, bp[bn]); /* bp[bn] != 0 */
*cancel = res + z;
return sign;
}
@@ -115,6 +136,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
MPFR_ASSERTD (bn >= 0);
MPFR_ASSERTD (cn >= 0);
MPFR_ASSERTD (bp[bn] != cp[cn]);
+
+ /* |b| != |c|. If |b| < |c|: swap (bp,bn) and (cp,cn),
+ and take the opposite sign. */
if (bp[bn] < cp[cn])
{
mp_limb_t *tp;
@@ -128,6 +152,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
} /* MPFR_EXP(b) >= MPFR_EXP(c) */
else /* MPFR_EXP(b) < MPFR_EXP(c) */
{
+ /* We necessarily have |b| < |c|. Simulate a swap by reading the
+ parameters so that |(bp,bn)| > |(cp,cn)|. */
+
sign = -1;
diff_exp = - (mpfr_uexp_t) sdiff_exp;
@@ -138,119 +165,178 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
cn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS;
}
- /* now we have removed the identical upper limbs of b and c
- (can happen only when diff_exp = 0), and after the possible
- swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0,
- diff_exp = EXP(b) - EXP(c).
- */
+ /* Now we have removed the identical upper limbs of b and c
+ (when diff_exp = 0), and after the possible swap, we have |b| > |c|,
+ where b is represented by (bp,bn) and c is represented by (cp,cn).
+ The value diff_exp = EXP(b) - EXP(c) can be regarded as the number
+ of leading zeros of c, when aligned with b. */
+
+ /* When a limb of c is read from memory, the part that is not taken
+ into account for the operation with a limb bp[bn] of b will be put
+ in lastc, shifted to the leftmost part (for alignment with b):
+ [-------- bp[bn] --------][------- bp[bn-1] -------]
+ [-- old_lastc --][-------- cp[cn] --------]
+ [-- new_lastc --]
+ Note: if diff_exp == 0, then lastc will always remain 0. */
+ lastc = 0;
+
+ /* Compute the next limb difference, which cannot be 0 (dif >= 1). */
if (MPFR_LIKELY (diff_exp < GMP_NUMB_BITS))
{
cc = cp[cn] >> diff_exp;
/* warning: a shift by GMP_NUMB_BITS is not allowed by ISO C */
- if (diff_exp)
+ if (diff_exp != 0)
lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp);
cn--;
}
else
- diff_exp -= GMP_NUMB_BITS; /* cc = 0 */
-
- dif = bp[bn--] - cc; /* necessarily dif >= 1 */
- MPFR_ASSERTD(dif >= 1);
+ {
+ cc = 0;
+ diff_exp -= GMP_NUMB_BITS; /* remove GMP_NUMB_BITS leading zeros */
+ }
- /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */
+ MPFR_ASSERTD (bp[bn] >= cc); /* no borrow out in subtraction below */
+ dif = bp[bn--] - cc;
+ MPFR_ASSERTD (dif >= 1);
+ high_dif = 0;
+
+ /* The current difference, here and later, is expressed under the form
+ [high_dif][dif], where high_dif is 0 or 1, and dif is a limb.
+ Here, since we have computed a difference of limbs (with b >= c),
+ high_dif = 0. */
+
+ /* One needs to accumulate canceled bits for the remaining case where
+ b and c are close to each other due to a long borrow propagation:
+ b = [common part]1000...000[low(b)]
+ c = [common part]0111...111[low(c)]
+ After eliminating the common part above, we have computed a difference
+ of the most significant parts, which has been stored in [high_dif][dif]
+ with high_dif = 0. We will loop as long as the currently computed
+ difference [high_dif][dif] = 1 (it is >= 1 by construction). The
+ computation of the difference will be:
+ 1bbb...bbb
+ - ccc...ccc
+ where the leading 1 before bbb...bbb corresponds to [high_dif][dif]
+ at the beginning of the loop. We will exit the loop also when c has
+ entirely been taken into account as cancellation is no longer possible
+ in this case (it is no longer possible to cancel the leading 1).
+ Note: We can enter the loop only with diff_exp = 0 (with a non-empty
+ common part, partly or entirely removed) or with diff_exp = 1 (with
+ an empty common part). Indeed, if diff_exp > 1, then no limbs have
+ been skipped, so that bp[bn] had its MSB equal to 1 and the most two
+ significant bits of cc are 0, which implies that dif > 1. */
while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
- && (high_dif == 0) && (dif == 1)))
- { /* dif=1 implies diff_exp = 0 or 1 */
- MPFR_ASSERTD (diff_exp <= 1);
- bb = (bn >= 0) ? bp[bn] : 0;
- cc = lastc;
- if (cn >= 0)
+ && high_dif == 0 && dif == 1))
+ {
+ /* Since we consider the next limb, we assume a cancellation of
+ GMP_NUMB_BITS (the new exponent of the difference now being the
+ one of the MSB of the next limb). But if the leading 1 remains
+ 1 in the difference (i.e. high_dif = 1 at the end of the loop),
+ then we will need to decrease res. */
+ res += GMP_NUMB_BITS;
+ MPFR_ASSERTD (diff_exp <= 1); /* see comment before the loop */
+ bb = bn >= 0 ? bp[bn--] : 0; /* next limb of b or non-represented 0 */
+ if (MPFR_UNLIKELY (cn < 0))
{
- if (diff_exp == 0)
- {
- cc += cp[cn];
- }
- else
- {
- MPFR_ASSERTD (diff_exp == 1);
- cc += cp[cn] >> 1;
- lastc = cp[cn] << (GMP_NUMB_BITS - 1);
- }
+ cc = lastc;
+ lastc = 0;
+ }
+ else if (diff_exp == 0)
+ {
+ cc = cp[cn--];
}
else
- lastc = 0;
- high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
- bn--;
- cn--;
- res += GMP_NUMB_BITS;
+ {
+ MPFR_ASSERTD (diff_exp == 1);
+ MPFR_ASSERTD (lastc == 0 || lastc == MPFR_LIMB_HIGHBIT);
+ cc = lastc + (cp[cn] >> 1);
+ lastc = cp[cn--] << (GMP_NUMB_BITS - 1);
+ }
+ dif = bb - cc;
+ high_dif = bb >= cc;
}
- /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
+ /* Now, c has entirely been taken into account or [high_dif][dif] > 1.
+ In any case, [high_dif][dif] >= 1 by construction.
+ First, we determine the currently number of canceled bits,
+ corresponding to the exponent of the current difference.
+ The trailing bits of c, if any, can still decrease the exponent of
+ the difference when [high_dif][dif] is a power of two, but since
+ [high_dif][dif] > 1 in this case, by not more than 1. */
- if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */
+ if (high_dif != 0) /* high_dif == 1 */
{
- res--;
+ res--; /* see comment at the beginning of the above loop */
MPFR_ASSERTD (res >= 0);
- if (dif != 0)
- {
- *cancel = res;
- return sign;
- }
+ /* Terminate if [high_dif][dif] is not a power of two. */
+ if (MPFR_LIKELY (dif != 0))
+ goto end;
}
else /* high_dif == 0 */
{
int z;
- count_leading_zeros (z, dif); /* dif > 1 here */
+ MPFR_ASSERTD (dif >= 1); /* [high_dif][dif] >= 1 */
+ count_leading_zeros (z, dif);
res += z;
- if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (GMP_NUMB_BITS - z - 1))))
- { /* dif is not a power of two */
- *cancel = res;
- return sign;
- }
+ /* Terminate if [high_dif][dif] is not a power of two. */
+ if (MPFR_LIKELY (NOT_POW2 (dif)))
+ goto end;
}
- /* now result is res + (low(b) < low(c)) */
- while (bn >= 0 && (cn >= 0 || lastc != 0))
+ /* Now, the result will be res + (low(b) < low(c)). */
+
+ /* If c has entirely been taken into account, it can no longer modify
+ the current result. */
+ if (cn < 0 && lastc == 0)
+ goto end;
+
+ for (; bn >= 0 ; bn--)
{
if (diff_exp >= GMP_NUMB_BITS)
- diff_exp -= GMP_NUMB_BITS;
- else
+ {
+ diff_exp -= GMP_NUMB_BITS;
+ MPFR_ASSERTD (cc == 0);
+ }
+ else if (MPFR_UNLIKELY (cn < 0))
{
cc = lastc;
- if (cn >= 0)
- {
- cc += cp[cn] >> diff_exp;
- if (diff_exp != 0)
- lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp);
- }
- else
- lastc = 0;
- cn--;
+ lastc = 0;
+ }
+ else if (diff_exp == 0)
+ {
+ cc = cp[cn--];
+ }
+ else
+ {
+ MPFR_ASSERTD (diff_exp >= 1 && diff_exp < GMP_NUMB_BITS);
+ cc = lastc + (cp[cn] >> diff_exp);
+ lastc = cp[cn--] << (GMP_NUMB_BITS - diff_exp);
}
+
if (bp[bn] != cc)
{
- *cancel = res + (bp[bn] < cc);
- return sign;
+ res += bp[bn] < cc;
+ goto end;
}
- bn--;
}
- if (bn < 0)
+ /* b has entirely been read. Determine whether the trailing part of c
+ is non-zero. */
+
+ if (lastc != 0)
+ res++;
+ else
{
- if (lastc != 0)
+ while (cn >= 0 && cp[cn] == 0)
+ cn--;
+ if (cn >= 0)
res++;
- else
- {
- while (cn >= 0 && cp[cn] == 0)
- cn--;
- if (cn >= 0)
- res++;
- }
}
+ end:
*cancel = res;
return sign;
}