diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/cmp2.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/cmp2.c | 254 |
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; } |