diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sum.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/sum.c | 1395 |
1 files changed, 0 insertions, 1395 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sum.c b/Build/source/libs/mpfr/mpfr-src/src/sum.c deleted file mode 100644 index ca3cba48c74..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/sum.c +++ /dev/null @@ -1,1395 +0,0 @@ -/* Sum -- efficiently sum a list of floating-point numbers - -Copyright 2014-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. */ - -#define MPFR_NEED_LONGLONG_H -#include "mpfr-impl.h" - -/* Note: In the prototypes, one uses - * - * const mpfr_ptr *x i.e.: __mpfr_struct *const *x - * - * instead of - * - * const mpfr_srcptr *x i.e.: const __mpfr_struct *const *x - * - * because here one has a double indirection and the type matching rules - * from the C standard in such a case are stricter and they would yield - * annoying errors for the user in practice. See: - * - * Why can't I pass a char ** to a function which expects a const char **? - * - * in the comp.lang.c FAQ: - * - * http://c-faq.com/ansi/constmismatch.html - */ - -/* See the doc/sum.txt file for the algorithm and a part of its proof -(this will later go into algorithms.tex). - -TODO [VL, after a discussion with James Demmel]: Compared to - James Demmel and Yozo Hida, Fast and accurate floating-point summation - with application to computational geometry, Numerical Algorithms, - volume 37, number 1-4, pages 101--112, 2004. -sorting is not necessary here. It is not done because in the most common -cases (where big cancellations are rare), it would take time and be -useless. However the lack of sorting increases the worst case complexity. -For instance, consider many inputs that cancel one another (two by two). -One would need n/2 iterations, where each iteration reads the exponent -of each input, therefore n*n/2 read operations. Using a worst-case sort -in O(n log n) could give a O(n log n) worst-case complexity. As we don't -want to slow down the most common cases, this could be done at the 3rd -iteration. But are there practical applications which would be used as -tests? - -Note: see the following paper and its references: -http://www.eecs.berkeley.edu/~hdnguyen/public/papers/ARITH21_Fast_Sum.pdf -VL: This is very different: - In MPFR In the paper & references - arbitrary precision fixed precision - correct rounding just reproducible rounding - integer operations floating-point operations - sequential parallel (& sequential) -*/ - -#ifdef MPFR_COV_CHECK -int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2] = { 0 }; -#endif - -/* Update minexp (V) after detecting a potential integer overflow in - extreme cases (only a 32-bit ABI may be concerned in practice). - Instead of an assertion failure below, we could - 1. check that the ulp of each regular input has an exponent >= MPFR_EXP_MIN - (with an assertion failure if this is not the case); - 2. set minexp to MPFR_EXP_MIN and shift the accumulator accordingly - (the sum will then be exact). - However, such cases, which involve huge precisions, will probably - never occur in practice (at least with a 64-bit ABI) and are not - easily testable due to these huge precisions. Moreover, switching - to a 64-bit ABI would be a better solution for such computations. - So, let's leave this unimplemented. */ -#define SAFE_SUB(V,E,SH) \ - do \ - { \ - mpfr_prec_t sh = (SH); \ - MPFR_ASSERTN ((E) >= MPFR_EXP_MIN + sh); \ - V = (E) - sh; \ - } \ - while (0) - -/* Function sum_raw - * ================ - * - * Accumulate a new [minexp,maxexp[ block into (wp,ws). If e and err denote - * the exponents of the computed result and of the error bound respectively, - * while e - err is less than some given bound (due to cancellation), shift - * the accumulator and reiterate. - * - * Inputs: - * wp: pointer to the accumulator (least significant limb first). - * ws: size of the accumulator (in limbs). - * wq: precision of the accumulator (ws * GMP_NUMB_BITS). - * x: array of the input numbers. - * n: size of this array (number of inputs, regular or not). - * minexp: exponent of the least significant bit of the first block. - * maxexp: exponent of the first block (exponent of its MSB + 1). - * tp: pointer to a temporary area (pre-allocated). - * ts: size of this temporary area. - * logn: ceil(log2(rn)), where rn is the number of regular inputs. - * prec: lower bound for e - err (as described above). - * ep: pointer to mpfr_exp_t (see below), or a null pointer. - * minexpp: pointer to mpfr_exp_t (see below), or a null pointer. - * maxexpp: pointer to mpfr_exp_t (see below), or a null pointer. - * - * Preconditions: - * prec >= 1 - * wq >= logn + prec + 2 - * - * This function returns 0 if the accumulator is 0 (which implies that - * the exact sum for this sum_raw invocation is 0), otherwise the number - * of cancelled bits (>= 1), defined as the number of identical bits on - * the most significant part of the accumulator. In the latter case, it - * also returns the following data in variables passed by reference, if - * the pointers are not NULL: - * - in ep: the exponent e of the computed result; - * - in minexpp: the last value of minexp; - * - in maxexpp: the new value of maxexp (for the next iteration after - * the first invocation of sum_raw in the main code). - * - * Notes: - * - minexp is also the exponent of the least significant bit of the - * accumulator; - * - the temporary area must be large enough to hold a shifted input - * block, and the value of ts is used only when the full assertions - * are checked (i.e. with the --enable-assert configure option), to - * check that a buffer overflow doesn't occur; - * - contrary to the returned value of minexp (the value in the last - * iteration), the returned value of maxexp is the one for the next - * iteration (= maxexp2 of the last iteration). - */ -static mpfr_prec_t -sum_raw (mp_limb_t *wp, mp_size_t ws, mpfr_prec_t wq, const mpfr_ptr *x, - unsigned long n, mpfr_exp_t minexp, mpfr_exp_t maxexp, - mp_limb_t *tp, mp_size_t ts, int logn, mpfr_prec_t prec, - mpfr_exp_t *ep, mpfr_exp_t *minexpp, mpfr_exp_t *maxexpp) -{ - MPFR_LOG_FUNC - (("ws=%Pd ts=%Pd prec=%Pd", (mpfr_prec_t) ws, (mpfr_prec_t) ts, prec), - ("", 0)); - - /* The C code below requires prec >= 0 due to the use of unsigned - integer arithmetic on it. Actually the computation makes sense - only with prec >= 1 (otherwise one can't even know the sign of - the result), hence the following assertion. */ - MPFR_ASSERTD (prec >= 1); - - /* Consistency check. */ - MPFR_ASSERTD (wq == (mpfr_prec_t) ws * GMP_NUMB_BITS); - - /* The following precondition together with prec >= 1 will imply: - minexp - shiftq < maxexp2, as required by the algorithm. */ - MPFR_ASSERTD (wq >= logn + prec + 2); - - while (1) - { - mpfr_exp_t maxexp2 = MPFR_EXP_MIN; - unsigned long i; - - MPFR_LOG_MSG (("sum_raw loop: " - "maxexp=%" MPFR_EXP_FSPEC "d " - "minexp=%" MPFR_EXP_FSPEC "d\n", - (mpfr_eexp_t) maxexp, (mpfr_eexp_t) minexp)); - - MPFR_ASSERTD (maxexp > minexp); - - for (i = 0; i < n; i++) - if (! MPFR_IS_SINGULAR (x[i])) /* Step 1 (see sum_raw in sum.txt) */ - { - mp_limb_t *dp, *vp; - mp_size_t ds, vs, vds; - mpfr_exp_t xe, vd; - mpfr_prec_t xq; - int tr; - - xe = MPFR_GET_EXP (x[i]); - xq = MPFR_GET_PREC (x[i]); - - vp = MPFR_MANT (x[i]); - vs = MPFR_PREC2LIMBS (xq); - vd = xe - vs * GMP_NUMB_BITS - minexp; - /* vd is the exponent of the least significant represented bit of - x[i] (including the trailing bits, whose value is 0) minus the - exponent of the least significant bit of the accumulator. To - make the code simpler, we won't try to filter out the trailing - bits of x[i]. */ - - /* Steps 2, 3, 4 (see sum_raw in sum.txt) */ - - if (vd < 0) - { - /* This covers the following cases: - * [-+- accumulator ---] - * [---|----- x[i] ------|--] - * | [----- x[i] --|--] - * | |[----- x[i] -----] - * | | [----- x[i] -----] - * maxexp minexp - */ - - /* Step 2 for subcase vd < 0 */ - - if (xe <= minexp) - { - /* x[i] is entirely after the LSB of the accumulator, - so that it will be ignored at this iteration. */ - if (xe > maxexp2) - { - maxexp2 = xe; - /* And since the exponent of x[i] is valid... */ - MPFR_ASSERTD (maxexp2 >= MPFR_EMIN_MIN); - } - continue; - } - - /* Step 3 for subcase vd < 0 */ - - /* If some significant bits of x[i] are after the LSB of the - accumulator, then maxexp2 will necessarily be minexp. */ - if (MPFR_LIKELY (xe - xq < minexp)) - maxexp2 = minexp; - - /* Step 4 for subcase vd < 0 */ - - /* We need to ignore the least |vd| significant bits of x[i]. - First, let's ignore the least vds = |vd| / GMP_NUMB_BITS - limbs. */ - vd = - vd; - vds = vd / GMP_NUMB_BITS; - vs -= vds; - MPFR_ASSERTD (vs > 0); /* see xe <= minexp test above */ - vp += vds; - vd -= vds * GMP_NUMB_BITS; - MPFR_ASSERTD (vd >= 0 && vd < GMP_NUMB_BITS); - - if (xe > maxexp) - { - vs -= (xe - maxexp) / GMP_NUMB_BITS; - MPFR_ASSERTD (vs > 0); - tr = (xe - maxexp) % GMP_NUMB_BITS; - } - else - tr = 0; - - if (vd != 0) - { - MPFR_ASSERTD (vs <= ts); - mpn_rshift (tp, vp, vs, vd); - vp = tp; - tr += vd; - if (tr >= GMP_NUMB_BITS) - { - vs--; - tr -= GMP_NUMB_BITS; - } - MPFR_ASSERTD (vs >= 1); - MPFR_ASSERTD (tr >= 0 && tr < GMP_NUMB_BITS); - if (tr != 0) - { - tp[vs-1] &= MPFR_LIMB_MASK (GMP_NUMB_BITS - tr); - tr = 0; - } - /* Truncation has now been taken into account. */ - MPFR_ASSERTD (tr == 0); - } - - dp = wp; - ds = ws; - } - else /* vd >= 0 */ - { - /* This covers the following cases: - * [-+- accumulator ---] - * [- x[i] -] | | - * [---|-- x[i] ------] | - * [------|-- x[i] ---------] - * | [- x[i] -] | - * maxexp minexp - */ - - /* Steps 2 and 3 for subcase vd >= 0 */ - - MPFR_ASSERTD (xe - xq >= minexp); /* see definition of vd */ - - /* Step 4 for subcase vd >= 0 */ - - /* We need to ignore the least vd significant bits - of the accumulator. First, let's ignore the least - vds = vd / GMP_NUMB_BITS limbs. -> (dp,ds) */ - vds = vd / GMP_NUMB_BITS; - ds = ws - vds; - if (ds <= 0) - continue; - dp = wp + vds; - vd -= vds * GMP_NUMB_BITS; - MPFR_ASSERTD (vd >= 0 && vd < GMP_NUMB_BITS); - - /* The low part of x[i] (to be determined) will have to be - shifted vd bits to the left if vd != 0. */ - - if (xe > maxexp) - { - vs -= (xe - maxexp) / GMP_NUMB_BITS; - if (vs <= 0) - continue; - tr = (xe - maxexp) % GMP_NUMB_BITS; - } - else - tr = 0; - - MPFR_ASSERTD (tr >= 0 && tr < GMP_NUMB_BITS && vs > 0); - - /* We need to consider the least significant vs limbs of x[i] - except the most significant tr bits. */ - - if (vd != 0) - { - mp_limb_t carry; - - MPFR_ASSERTD (vs <= ts); - carry = mpn_lshift (tp, vp, vs, vd); - tr -= vd; - if (tr < 0) - { - tr += GMP_NUMB_BITS; - MPFR_ASSERTD (vs + 1 <= ts); - tp[vs++] = carry; - } - MPFR_ASSERTD (tr >= 0 && tr < GMP_NUMB_BITS); - vp = tp; - } - } /* vd >= 0 */ - - MPFR_ASSERTD (vs > 0 && vs <= ds); - - /* We can't truncate the most significant limb of the input - (in case it hasn't been shifted to the temporary area). - So, let's ignore it now. It will be taken into account - via carry propagation after the addition. */ - if (tr != 0) - vs--; - - /* Step 5 (see sum_raw in sum.txt) */ - - if (MPFR_IS_POS (x[i])) - { - mp_limb_t carry; - - carry = vs > 0 ? mpn_add_n (dp, dp, vp, vs) : 0; - MPFR_ASSERTD (carry <= 1); - if (tr != 0) - carry += vp[vs] & MPFR_LIMB_MASK (GMP_NUMB_BITS - tr); - if (ds > vs) - mpn_add_1 (dp + vs, dp + vs, ds - vs, carry); - } - else - { - mp_limb_t borrow; - - borrow = vs > 0 ? mpn_sub_n (dp, dp, vp, vs) : 0; - MPFR_ASSERTD (borrow <= 1); - if (tr != 0) - borrow += vp[vs] & MPFR_LIMB_MASK (GMP_NUMB_BITS - tr); - if (ds > vs) - mpn_sub_1 (dp + vs, dp + vs, ds - vs, borrow); - } - } - - { - mpfr_prec_t cancel; /* number of cancelled bits */ - mp_size_t wi; /* index in the accumulator */ - mp_limb_t a, b; - int cnt; - - cancel = 0; - wi = ws - 1; - MPFR_ASSERTD (wi >= 0); - a = wp[wi] >> (GMP_NUMB_BITS - 1) ? MPFR_LIMB_MAX : MPFR_LIMB_ZERO; - - while (wi >= 0) - if ((b = wp[wi]) == a) - { - cancel += GMP_NUMB_BITS; - wi--; - } - else - { - b ^= a; - MPFR_ASSERTD (b != 0); - count_leading_zeros (cnt, b); - cancel += cnt; - break; - } - - if (wi >= 0 || a != MPFR_LIMB_ZERO) /* accumulator != 0 */ - { - mpfr_exp_t e; /* exponent of the computed result */ - mpfr_exp_t err; /* exponent of the error bound */ - - MPFR_LOG_MSG (("accumulator %s 0, cancel=%Pd\n", - a != MPFR_LIMB_ZERO ? "<" : ">", cancel)); - - MPFR_ASSERTD (cancel > 0); - e = minexp + wq - cancel; - MPFR_ASSERTD (e >= minexp); - err = maxexp2 + logn; /* OK even if maxexp2 == MPFR_EXP_MIN */ - - /* The absolute value of the truncated sum is in the binade - [2^(e-1),2^e] (closed on both ends due to two's complement). - The error is strictly less than 2^err (and is 0 if - maxexp2 == MPFR_EXP_MIN). */ - - /* This basically tests whether err <= e - prec without - potential integer overflow (since prec >= 0)... - Note that the maxexp2 == MPFR_EXP_MIN test is there just for - the potential corner case e - prec < MPFR_EXP_MIN + logn. - Such corner cases, involving specific huge-precision numbers, - are probably not supported in many places in MPFR, but this - test doesn't hurt... */ - if (maxexp2 == MPFR_EXP_MIN || - (err <= e && SAFE_DIFF (mpfr_uexp_t, e, err) >= prec)) - { - MPFR_LOG_MSG (("(err=%" MPFR_EXP_FSPEC "d) <= (e=%" - MPFR_EXP_FSPEC "d) - (prec=%Pd)\n", - (mpfr_eexp_t) err, (mpfr_eexp_t) e, prec)); - /* To avoid tests or copies, we consider the only two cases - that will occur in sum_aux. */ - MPFR_ASSERTD ((ep != NULL && - minexpp != NULL && - maxexpp != NULL) || - (ep == NULL && - minexpp == NULL && - maxexpp == NULL)); - if (ep != NULL) - { - *ep = e; - *minexpp = minexp; - *maxexpp = maxexp2; - } - MPFR_LOG_MSG (("return with minexp=%" MPFR_EXP_FSPEC - "d maxexp2=%" MPFR_EXP_FSPEC "d%s\n", - (mpfr_eexp_t) minexp, (mpfr_eexp_t) maxexp2, - maxexp2 == MPFR_EXP_MIN ? - " (MPFR_EXP_MIN)" : "")); - return cancel; - } - else - { - mpfr_exp_t diffexp; - mpfr_prec_t shiftq; - mpfr_size_t shifts; - int shiftc; - - MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d err=%" MPFR_EXP_FSPEC - "d maxexp2=%" MPFR_EXP_FSPEC "d%s\n", - (mpfr_eexp_t) e, (mpfr_eexp_t) err, - (mpfr_eexp_t) maxexp2, - maxexp2 == MPFR_EXP_MIN ? - " (MPFR_EXP_MIN)" : "")); - - diffexp = err - e; - if (diffexp < 0) - diffexp = 0; - /* diffexp = max(0, err - e) */ - - MPFR_LOG_MSG (("diffexp=%" MPFR_EXP_FSPEC "d\n", - (mpfr_eexp_t) diffexp)); - - MPFR_ASSERTD (diffexp < cancel - 2); - shiftq = cancel - 2 - (mpfr_prec_t) diffexp; - /* equivalent to: minexp + wq - 2 - max(e,err) */ - MPFR_ASSERTD (shiftq > 0); - shifts = shiftq / GMP_NUMB_BITS; - shiftc = shiftq % GMP_NUMB_BITS; - MPFR_LOG_MSG (("shiftq = %Pd = %Pd * GMP_NUMB_BITS + %d\n", - shiftq, (mpfr_prec_t) shifts, shiftc)); - if (MPFR_LIKELY (shiftc != 0)) - mpn_lshift (wp + shifts, wp, ws - shifts, shiftc); - else - mpn_copyd (wp + shifts, wp, ws - shifts); - MPN_ZERO (wp, shifts); - /* Compute minexp = minexp - shiftq safely. */ - SAFE_SUB (minexp, minexp, shiftq); - MPFR_ASSERTD (minexp < maxexp2); - } - } - else if (maxexp2 == MPFR_EXP_MIN) - { - MPFR_LOG_MSG (("accumulator = 0, maxexp2 = MPFR_EXP_MIN\n", 0)); - return 0; - } - else - { - MPFR_LOG_MSG (("accumulator = 0, reiterate\n", 0)); - /* Compute minexp = maxexp2 - (wq - (logn + 1)) safely. */ - SAFE_SUB (minexp, maxexp2, wq - (logn + 1)); - /* Note: the logn + 1 corresponds to cq in the main code. */ - } - } - - maxexp = maxexp2; - } -} - -/**********************************************************************/ - -/* Generic case: all the inputs are finite numbers, - with at least 3 regular numbers. */ -static int -sum_aux (mpfr_ptr sum, const mpfr_ptr *x, unsigned long n, mpfr_rnd_t rnd, - mpfr_exp_t maxexp, unsigned long rn) -{ - mp_limb_t *sump; - mp_limb_t *tp; /* pointer to a temporary area */ - mp_limb_t *wp; /* pointer to the accumulator */ - mp_size_t ts; /* size of the temporary area, in limbs */ - mp_size_t ws; /* size of the accumulator, in limbs */ - mp_size_t zs; /* size of the TMD accumulator, in limbs */ - mpfr_prec_t wq; /* size of the accumulator, in bits */ - int logn; /* ceil(log2(rn)) */ - int cq; - mpfr_prec_t sq; - int inex; - MPFR_TMP_DECL (marker); - - MPFR_LOG_FUNC - (("n=%lu rnd=%d maxexp=%" MPFR_EXP_FSPEC "d rn=%lu", - n, rnd, (mpfr_eexp_t) maxexp, rn), - ("sum[%Pu]=%.*Rg", mpfr_get_prec (sum), mpfr_log_prec, sum)); - - MPFR_ASSERTD (rn >= 3 && rn <= n); - - /* In practice, no integer overflow on the exponent. */ - MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX - MPFR_EMAX_MAX >= - sizeof (unsigned long) * CHAR_BIT); - - /* Set up some variables and the accumulator. */ - - sump = MPFR_MANT (sum); - - /* rn is the number of regular inputs (the singular ones will be - ignored). Compute logn = ceil(log2(rn)). */ - logn = MPFR_INT_CEIL_LOG2 (rn); - MPFR_ASSERTD (logn >= 2); - - MPFR_LOG_MSG (("logn=%d maxexp=%" MPFR_EXP_FSPEC "d\n", - logn, (mpfr_eexp_t) maxexp)); - - sq = MPFR_GET_PREC (sum); - cq = logn + 1; - - /* First determine the size of the accumulator. - * cq + sq + logn + 2 >= logn + sq + 5, which will be used later. - * The assertion wq - cq - sq >= 4 is another way to check that. - */ - ws = MPFR_PREC2LIMBS (cq + sq + logn + 2); - wq = (mpfr_prec_t) ws * GMP_NUMB_BITS; - MPFR_ASSERTD (wq - cq - sq >= 4); - - /* TODO: timings, comparing with a larger zs. */ - zs = MPFR_PREC2LIMBS (wq - sq); - - MPFR_LOG_MSG (("cq=%d sq=%Pd logn=%d wq=%Pd\n", cq, sq, logn, wq)); - - /* An input block will have up to wq - cq bits, and its shifted value - (to be correctly aligned) may take GMP_NUMB_BITS - 1 additional bits. */ - ts = MPFR_PREC2LIMBS (wq - cq + GMP_NUMB_BITS - 1); - - MPFR_TMP_MARK (marker); - - /* Note: If the TMD does not occur, which should be the case for most - sums, allocating zs limbs is not necessary. However, we choose to - do this now (thus in all cases) because zs is very small, so that - the difference on the memory footprint will not be noticeable. - More precisely, zs is at most 2 in practice with the current code; - we may want to increase it in order to avoid performance issues in - some unlikely corner cases, but even in this case, it will remain - small. - One will have: - [------ ts ------][------ ws ------][- zs -] - The following would probably be better: - [------ ts ------] [------ ws ------] - [- zs -] - i.e. where the TMD accumulator (partially or completely) takes - some unneeded part of the temporary area in order to improve - data locality. But - * in low precision, data locality is regarded as ensured even - with the actual choice; - * in high precision, data locality for TMD resolution may not - be that important. - */ - tp = MPFR_TMP_LIMBS_ALLOC (ts + ws + zs); - wp = tp + ts; - - MPN_ZERO (wp, ws); /* zero the accumulator */ - - { - mpfr_exp_t minexp; /* exponent of the LSB of the block for sum_raw */ - mpfr_prec_t cancel; /* number of cancelled bits */ - mpfr_exp_t e; /* temporary exponent of the result */ - mpfr_exp_t u; /* temporary exponent of the ulp (quantum) */ - mp_limb_t lbit; /* last bit (useful if even rounding) */ - mp_limb_t rbit; /* rounding bit (corrected in halfway case) */ - int corr; /* correction term (from -1 to 2) */ - int sd, sh; /* shift counts */ - mp_size_t sn; /* size of the output number */ - int tmd; /* 0: the TMD does not occur - 1: the TMD occurs on a machine number - 2: the TMD occurs on a midpoint */ - int neg; /* 0 if positive sum, 1 if negative */ - int sgn; /* +1 if positive sum, -1 if negative */ - - MPFR_LOG_MSG (("Compute an approximation with sum_raw...\n", 0)); - - /* Compute minexp = maxexp - (wq - cq) safely. */ - SAFE_SUB (minexp, maxexp, wq - cq); - MPFR_ASSERTD (wq >= logn + sq + 5); - cancel = sum_raw (wp, ws, wq, x, n, minexp, maxexp, tp, ts, - logn, sq + 3, &e, &minexp, &maxexp); - - if (MPFR_UNLIKELY (cancel == 0)) - { - /* The exact sum is zero. Since not all inputs are 0, the sum - * is +0 except in MPFR_RNDD, as specified according to the - * IEEE 754 rules for the addition of two numbers. - */ - MPFR_SET_SIGN (sum, (rnd != MPFR_RNDD ? - MPFR_SIGN_POS : MPFR_SIGN_NEG)); - MPFR_SET_ZERO (sum); - MPFR_TMP_FREE (marker); - MPFR_RET (0); - } - - /* The absolute value of the truncated sum is in the binade - [2^(e-1),2^e] (closed on both ends due to two's complement). - The error is strictly less than 2^(maxexp + logn) (and is 0 - if maxexp == MPFR_EXP_MIN). */ - - u = e - sq; /* e being the exponent, u is the ulp of the target */ - - /* neg = 1 if negative, 0 if positive. */ - neg = wp[ws-1] >> (GMP_NUMB_BITS - 1); - MPFR_ASSERTD (neg == 0 || neg == 1); - - sgn = neg ? -1 : 1; - MPFR_ASSERTN (sgn == (neg ? MPFR_SIGN_NEG : MPFR_SIGN_POS)); - - MPFR_LOG_MSG (("neg=%d sgn=%d cancel=%Pd" - " e=%" MPFR_EXP_FSPEC "d" - " u=%" MPFR_EXP_FSPEC "d" - " maxexp=%" MPFR_EXP_FSPEC "d%s\n", - neg, sgn, cancel, (mpfr_eexp_t) e, (mpfr_eexp_t) u, - (mpfr_eexp_t) maxexp, - maxexp == MPFR_EXP_MIN ? " (MPFR_EXP_MIN)" : "")); - - if (rnd == MPFR_RNDF) - { - /* Rounding the approximate value to nearest (ties don't matter) is - sufficient. We need to get the rounding bit; the code is similar - to a part from the generic code (here, corr = rbit). */ - if (MPFR_LIKELY (u > minexp)) - { - mpfr_prec_t tq; - mp_size_t wi; - int td; - - tq = u - minexp; - MPFR_ASSERTD (tq > 0); /* number of trailing bits */ - MPFR_LOG_MSG (("tq=%Pd\n", tq)); - - wi = tq / GMP_NUMB_BITS; - td = tq % GMP_NUMB_BITS; - corr = td >= 1 ? ((wp[wi] >> (td - 1)) & MPFR_LIMB_ONE) : - (MPFR_ASSERTD (wi >= 1), wp[wi-1] >> (GMP_NUMB_BITS - 1)); - } - else - corr = 0; - inex = 0; /* not meaningful, but needs to have a value */ - } - else /* rnd != MPFR_RNDF */ - { - if (MPFR_LIKELY (u > minexp)) - { - mpfr_prec_t tq; - mp_size_t wi; - int td; - - tq = u - minexp; - MPFR_ASSERTD (tq > 0); /* number of trailing bits */ - MPFR_LOG_MSG (("tq=%Pd\n", tq)); - - wi = tq / GMP_NUMB_BITS; - - /* Determine the rounding bit, which is represented. */ - td = tq % GMP_NUMB_BITS; - lbit = (wp[wi] >> td) & MPFR_LIMB_ONE; - rbit = td >= 1 ? ((wp[wi] >> (td - 1)) & MPFR_LIMB_ONE) : - (MPFR_ASSERTD (wi >= 1), wp[wi-1] >> (GMP_NUMB_BITS - 1)); - MPFR_ASSERTD (rbit == 0 || rbit == 1); - MPFR_LOG_MSG (("rbit=%d\n", (int) rbit)); - - if (maxexp == MPFR_EXP_MIN) - { - /* The sum in the accumulator is exact. Determine inex: - inex = 0 if the final sum is exact, else 1, i.e. - inex = rounding bit || sticky bit. In round to nearest, - also determine the rounding direction: obtained from - the rounding bit possibly except in halfway cases. - Halfway cases are rounded toward -inf iff the last bit - of the truncated significand in two's complement is 0 - (in precision > 1, because the parity after rounding is - the same in two's complement and sign + magnitude; in - precision 1, one checks that the rule works for both - positive (lbit == 1) and negative (lbit == 0) numbers, - rounding halfway cases away from zero). */ - if (MPFR_LIKELY (rbit == 0 || (rnd == MPFR_RNDN && lbit == 0))) - { - /* We need to determine the sticky bit, either to set inex - (if the rounding bit is 0) or to possibly "correct" rbit - (round to nearest, halfway case rounded downward) from - which the rounding direction will be determined. */ - MPFR_LOG_MSG (("Determine the sticky bit...\n", 0)); - - inex = td >= 2 ? (wp[wi] & MPFR_LIMB_MASK (td - 1)) != 0 - : td == 0 ? - (MPFR_ASSERTD (wi >= 1), - (wp[--wi] & MPFR_LIMB_MASK (GMP_NUMB_BITS - 1)) != 0) - : 0; - - if (!inex) - { - while (!inex && wi > 0) - inex = wp[--wi] != 0; - if (!inex && rbit != 0) - { - /* sticky bit = 0, rounding bit = 1, - i.e. halfway case, which will be - rounded downward (see earlier if). */ - MPFR_ASSERTD (rnd == MPFR_RNDN); - inex = 1; - rbit = 0; /* even rounding downward */ - MPFR_LOG_MSG (("Halfway case rounded downward;" - " set inex=1 rbit=0\n", 0)); - } - } - } - else - inex = 1; - tmd = 0; /* We can round correctly -> no TMD. */ - } - else /* maxexp > MPFR_EXP_MIN */ - { - mpfr_exp_t d; - mp_limb_t limb, mask; - int nbits; - - /* Since maxexp was set to either the exponent of a x[i] or - to minexp... */ - MPFR_ASSERTD (maxexp >= MPFR_EMIN_MIN || maxexp == minexp); - - inex = 1; /* We do not know whether the sum is exact. */ - - MPFR_ASSERTD (u <= MPFR_EMAX_MAX && u <= minexp + wq); - d = u - (maxexp + logn); /* representable */ - MPFR_ASSERTD (d >= 3); /* due to prec = sq + 3 in sum_raw */ - - /* Let's see whether the TMD occurs by looking at the d bits - following the ulp bit, or the d-1 bits after the rounding - bit. */ - - /* First chunk after the rounding bit... It starts at: - (wi,td-2) if td >= 2, - (wi-1,td-2+GMP_NUMB_BITS) if td < 2. */ - if (td == 0) - { - MPFR_ASSERTD (wi >= 1); - limb = wp[--wi]; - mask = MPFR_LIMB_MASK (GMP_NUMB_BITS - 1); - nbits = GMP_NUMB_BITS; - } - else if (td == 1) - { - limb = wi >= 1 ? wp[--wi] : MPFR_LIMB_ZERO; - mask = MPFR_LIMB_MAX; - nbits = GMP_NUMB_BITS + 1; - } - else /* td >= 2 */ - { - MPFR_ASSERTD (td >= 2); - limb = wp[wi]; - mask = MPFR_LIMB_MASK (td - 1); - nbits = td; - } - - /* nbits: number of bits of the first chunk + 1 - (the +1 is for the rounding bit). */ - - if (nbits > d) - { - /* Some low significant bits must be ignored. */ - limb >>= nbits - d; - mask >>= nbits - d; - d = 0; - } - else - { - d -= nbits; - MPFR_ASSERTD (d >= 0); - } - - limb &= mask; - tmd = - limb == MPFR_LIMB_ZERO ? - (rbit == 0 ? 1 : rnd == MPFR_RNDN ? 2 : 0) : - limb == mask ? - (limb = MPFR_LIMB_MAX, - rbit != 0 ? 1 : rnd == MPFR_RNDN ? 2 : 0) : 0; - - while (tmd != 0 && d != 0) - { - mp_limb_t limb2; - - MPFR_ASSERTD (d > 0); - if (wi == 0) - { - /* The non-represented bits are 0's. */ - if (limb != MPFR_LIMB_ZERO) - tmd = 0; - break; - } - MPFR_ASSERTD (wi > 0); - limb2 = wp[--wi]; - if (d < GMP_NUMB_BITS) - { - int c = GMP_NUMB_BITS - d; - MPFR_ASSERTD (c > 0 && c < GMP_NUMB_BITS); - if ((limb2 >> c) != (limb >> c)) - tmd = 0; - break; - } - if (limb2 != limb) - tmd = 0; - d -= GMP_NUMB_BITS; - } - } - } - else /* u <= minexp */ - { - /* The exact value of the accumulator will be copied. - * The TMD occurs if and only if there are bits still - * not taken into account, and if it occurs, this is - * necessarily on a machine number (-> tmd = 1). - */ - lbit = u == minexp ? wp[0] & MPFR_LIMB_ONE : 0; - rbit = 0; - inex = tmd = maxexp != MPFR_EXP_MIN; - } - - MPFR_ASSERTD (rbit == 0 || rbit == 1); - - MPFR_LOG_MSG (("tmd=%d lbit=%d rbit=%d inex=%d neg=%d\n", - tmd, (int) lbit, (int) rbit, inex, neg)); - - /* Here, if the final sum is known to be exact, inex = 0, otherwise - * inex = 1. We have a truncated significand, a trailing term t such - * that 0 <= t < 1 ulp, and an error on the trailing term bounded by - * t' in absolute value. Thus the error e on the truncated significand - * satisfies -t' <= e < 1 ulp + t'. Thus one has 4 correction cases - * denoted by a corr value between -1 and 2 depending on e, neg, rbit, - * and the rounding mode: - * -1: equivalent to nextbelow; - * 0: the truncated significand is not corrected; - * 1: add 1 ulp; - * 2: add 1 ulp, then nextabove. - * The nextbelow and nextabove are used here since there may be a - * change of the binade. - */ - - if (tmd == 0) /* no TMD */ - { - switch (rnd) - { - case MPFR_RNDD: - corr = 0; - break; - case MPFR_RNDU: - corr = inex; - break; - case MPFR_RNDZ: - corr = inex && neg; - break; - case MPFR_RNDA: - corr = inex && !neg; - break; - default: - MPFR_ASSERTN (rnd == MPFR_RNDN); - /* Note: for halfway cases (maxexp == MPFR_EXP_MIN) that are - rounded downward, rbit has been changed to 0 so that corr - is set correctly. */ - corr = rbit; - } - MPFR_ASSERTD (corr == 0 || corr == 1); - if (inex && - corr == 0) /* two's complement significand decreased */ - inex = -1; - } - else /* tmd */ - { - mpfr_exp_t minexp2; - mpfr_prec_t cancel2; - mpfr_exp_t err; /* exponent of the error bound */ - mp_size_t zz; /* nb of limbs to zero in the TMD accumulator */ - mp_limb_t *zp; /* pointer to the TMD accumulator */ - mpfr_prec_t zq; /* size of the TMD accumulator, in bits */ - int sst; /* sign of the secondary term */ - - /* TMD case. Here we use a new variable minexp2, with the same - meaning as minexp, as we want to keep the minexp value for - the copy to the destination. */ - - MPFR_ASSERTD (maxexp > MPFR_EXP_MIN); - MPFR_ASSERTD (tmd == 1 || tmd == 2); - - /* TMD accumulator */ - zp = wp + ws; - zq = (mpfr_prec_t) zs * GMP_NUMB_BITS; - - err = maxexp + logn; - - MPFR_LOG_MSG (("TMD with" - " maxexp=%" MPFR_EXP_FSPEC "d" - " err=%" MPFR_EXP_FSPEC "d" - " zs=%Pd" - " zq=%Pd\n", - (mpfr_eexp_t) maxexp, (mpfr_eexp_t) err, - (mpfr_prec_t) zs, zq)); - - /* The d-1 bits from u-2 to u-d (= err) are identical. */ - - if (err >= minexp) - { - mpfr_prec_t tq; - mp_size_t wi; - int td; - - /* Let's keep the last 2 over the d-1 identical bits and the - following bits, i.e. the bits from err+1 to minexp. */ - tq = err - minexp + 2; /* tq = number of such bits */ - MPFR_LOG_MSG (("[TMD] tq=%Pd\n", tq)); - MPFR_ASSERTD (tq >= 2); - - wi = tq / GMP_NUMB_BITS; - td = tq % GMP_NUMB_BITS; - - /* Note: The "else" (td == 0) branch below can be executed - only if tq >= GMP_NUMB_BITS, which is possible only when - logn is large enough. Indeed, if tq > logn + some constant, - this means that the TMD did not occur. - TODO: Find an upper bound on tq, and add a corresponding - MPFR_ASSERTD assertion / hint. On some platforms, this - branch might be dead code, and such information would - allow the compiler to remove it. - It seems that this branch is never tested (r12754). */ - - if (td != 0) - { - wi++; /* number of words with represented bits */ - td = GMP_NUMB_BITS - td; - zz = zs - wi; - MPFR_ASSERTD (zz >= 0 && zz < zs); - mpn_lshift (zp + zz, wp, wi, td); - } - else - { - MPFR_ASSERTD (wi > 0); - zz = zs - wi; - MPFR_ASSERTD (zz >= 0 && zz < zs); - if (zz > 0) - MPN_COPY (zp + zz, wp, wi); - } - - /* Compute minexp2 = minexp - (zs * GMP_NUMB_BITS + td) - safely. */ - SAFE_SUB (minexp2, minexp, zz * GMP_NUMB_BITS + td); - MPFR_ASSERTD (minexp2 == err + 2 - zq); - } - else /* err < minexp */ - { - /* At least one of the identical bits is not represented, - meaning that it is 0 and all these bits are 0's. Thus - the accumulator will be 0. The new minexp is determined - from maxexp, with cq bits reserved to avoid an overflow - (as in the early steps). */ - MPFR_LOG_MSG (("[TMD] err < minexp\n", 0)); - zz = zs; - - /* Compute minexp2 = maxexp - (zq - cq) safely. */ - SAFE_SUB (minexp2, maxexp, zq - cq); - MPFR_ASSERTD (minexp2 == err + 1 - zq); - } - - MPN_ZERO (zp, zz); - - /* We need to determine the sign sst of the secondary term. - In sum_raw, since the truncated sum corresponding to this - secondary term will be in [2^(e-1),2^e] and the error - strictly less than 2^err, we can stop the iterations when - e - err >= 1 (this bound is the 11th argument of sum_raw). */ - cancel2 = sum_raw (zp, zs, zq, x, n, minexp2, maxexp, tp, ts, - logn, 1, NULL, NULL, NULL); - - if (cancel2 != 0) - sst = MPFR_LIMB_MSB (zp[zs-1]) == 0 ? 1 : -1; - else if (tmd == 1) - sst = 0; - else - { - /* For halfway cases, let's virtually eliminate them - by setting a sst equivalent to a non-halfway case, - which depends on the last bit of the pre-rounded - result. */ - MPFR_ASSERTD (rnd == MPFR_RNDN && tmd == 2); - sst = lbit != 0 ? 1 : -1; - } - - MPFR_LOG_MSG (("[TMD] tmd=%d rbit=%d sst=%d\n", - tmd, (int) rbit, sst)); - - /* Do not consider the corrected sst for MPFR_COV_SET */ - MPFR_COV_SET (sum_tmd[(int) rnd][tmd-1][rbit] - [cancel2 == 0 ? 1 : sst+1][neg][sq > MPFR_PREC_MIN]); - - inex = - MPFR_IS_LIKE_RNDD (rnd, sgn) ? (sst ? -1 : 0) : - MPFR_IS_LIKE_RNDU (rnd, sgn) ? (sst ? 1 : 0) : - (MPFR_ASSERTD (rnd == MPFR_RNDN), - tmd == 1 ? - sst : sst); - - if (tmd == 2 && sst == (rbit != 0 ? -1 : 1)) - corr = 1 - (int) rbit; - else if (MPFR_IS_LIKE_RNDD (rnd, sgn) && sst == -1) - corr = (int) rbit - 1; - else if (MPFR_IS_LIKE_RNDU (rnd, sgn) && sst == +1) - corr = (int) rbit + 1; - else - corr = (int) rbit; - } /* tmd */ - } /* rnd != MPFR_RNDF */ - - MPFR_LOG_MSG (("neg=%d corr=%d inex=%d\n", neg, corr, inex)); - - /* Sign handling (-> absolute value and sign), together with - rounding. The most common cases are corr = 0 and corr = 1 - as this is necessarily the case when the TMD did not occur. */ - - MPFR_ASSERTD (corr >= -1 && corr <= 2); - - MPFR_SIGN (sum) = sgn; - - /* Let's copy/shift the bits [max(u,minexp),e) to the - most significant part of the destination, and zero - the least significant part (there can be one only if - u < minexp). The trailing bits of the destination may - contain garbage at this point. */ - - sn = MPFR_PREC2LIMBS (sq); - sd = (mpfr_prec_t) sn * GMP_NUMB_BITS - sq; - sh = cancel % GMP_NUMB_BITS; - - MPFR_ASSERTD (sd >= 0 && sd < GMP_NUMB_BITS); - - if (MPFR_LIKELY (u > minexp)) - { - mp_size_t wi; - - /* Recompute the initial value of wi. */ - wi = (u - minexp) / GMP_NUMB_BITS; - if (MPFR_LIKELY (sh != 0)) - { - mp_size_t fi; - - fi = (e - minexp) / GMP_NUMB_BITS - (sn - 1); - MPFR_ASSERTD (fi == wi || fi == wi + 1); - mpn_lshift (sump, wp + fi, sn, sh); - if (fi != wi) - sump[0] |= wp[wi] >> (GMP_NUMB_BITS - sh); - } - else - { - MPFR_ASSERTD ((mpfr_prec_t) (ws - (wi + sn)) * GMP_NUMB_BITS - == cancel); - MPN_COPY (sump, wp + wi, sn); - } - } - else /* u <= minexp */ - { - mp_size_t en; - - en = (e - minexp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS; - if (MPFR_LIKELY (sh != 0)) - mpn_lshift (sump + sn - en, wp, en, sh); - else if (MPFR_UNLIKELY (en > 0)) - MPN_COPY (sump + sn - en, wp, en); - if (sn > en) - MPN_ZERO (sump, sn - en); - } - - /* Let's take the complement if the result is negative, and at - the same time, do the rounding and zero the trailing bits. - As this is valid only for precisions >= 2, there is special - code for precision 1 first. */ - - if (MPFR_UNLIKELY (sq == 1)) /* precision 1 */ - { - sump[0] = MPFR_LIMB_HIGHBIT; - e += neg ? 1 - corr : corr; - } - else if (neg) /* negative result with sq > 1 */ - { - MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) == 0); - - /* abs(x + corr) = - (x + corr) = com(x) + (1 - corr) */ - - /* We want to avoid separate mpn_com (or mpn_neg) and mpn_add_1 - (or mpn_sub_1) operations, as they could yield two loops in - some particular cases involving a long sequence of 0's in - the low significant bits (except the least significant bit, - which doesn't matter). */ - - if (corr <= 1) - { - mp_limb_t corr2; - - /* Here we can just do the correction operation on the - least significant limb, then do either a mpn_com or - a mpn_neg on the remaining limbs, depending on the - carry (BTW, mpn_neg is just a mpn_com with an initial - carry propagation: after some point, mpn_neg does a - complement). */ - - corr2 = (mp_limb_t) (1 - corr) << sd; - /* Note: If corr = -1, this can overflow to corr2 = 0. - This case is taken into account below. */ - - sump[0] = (~ (sump[0] | MPFR_LIMB_MASK (sd))) + corr2; - - if (sump[0] < corr2 || (corr2 == 0 && corr < 0)) - { - if (sn == 1 || ! mpn_neg (sump + 1, sump + 1, sn - 1)) - { - /* Note: The | is important when sump[sn-1] is not 0 - (this can occur with sn = 1 and corr = -1). TODO: - Add something to make sure that this is tested. */ - sump[sn-1] |= MPFR_LIMB_HIGHBIT; - e++; - } - } - else if (sn > 1) - mpn_com (sump + 1, sump + 1, sn - 1); - } - else /* corr == 2 */ - { - mp_limb_t corr2, c; - mp_size_t i = 1; - - /* We want to compute com(x) - 1, but GMP doesn't have an - operation for that. The fact is that a sequence of low - significant bits 1 is invariant. Starting at the first - low significant bit 0, we can do the complement with - mpn_com. */ - - corr2 = MPFR_LIMB_ONE << sd; - c = ~ (sump[0] | MPFR_LIMB_MASK (sd)); - sump[0] = c - corr2; - - if (c == 0) - { - while (MPFR_ASSERTD (i < sn), sump[i] == MPFR_LIMB_MAX) - i++; - sump[i] = (~ sump[i]) - 1; - i++; - } - - if (i < sn) - mpn_com (sump + i, sump + i, sn - i); - else if (MPFR_UNLIKELY (MPFR_LIMB_MSB (sump[sn-1]) == 0)) - { - /* Happens on 01111...111, whose complement is - 10000...000, and com(x) - 1 is 01111...111. */ - sump[sn-1] |= MPFR_LIMB_HIGHBIT; - e--; - } - } - } - else /* positive result with sq > 1 */ - { - MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) != 0); - sump[0] &= ~ MPFR_LIMB_MASK (sd); - - if (corr > 0) - { - mp_limb_t corr2, carry_out; - - corr2 = (mp_limb_t) corr << sd; - /* If corr == 2 && sd == GMP_NUMB_BITS - 1, this overflows - to corr2 = 0. This case is taken into account below. */ - - carry_out = corr2 != 0 ? mpn_add_1 (sump, sump, sn, corr2) : - (MPFR_ASSERTD (sn > 1), - mpn_add_1 (sump + 1, sump + 1, sn - 1, MPFR_LIMB_ONE)); - - MPFR_ASSERTD (sump[sn-1] >> (GMP_NUMB_BITS - 1) == !carry_out); - - if (MPFR_UNLIKELY (carry_out)) - { - /* Note: The | is important when sump[sn-1] is not 0 - (this can occur with sn = 1 and corr = 2). TODO: - Add something to make sure that this is tested. */ - sump[sn-1] |= MPFR_LIMB_HIGHBIT; - e++; - } - } - - if (corr < 0) - { - mpn_sub_1 (sump, sump, sn, MPFR_LIMB_ONE << sd); - - if (MPFR_UNLIKELY (MPFR_LIMB_MSB (sump[sn-1]) == 0)) - { - sump[sn-1] |= MPFR_LIMB_HIGHBIT; - e--; - } - } - } - - MPFR_ASSERTD (MPFR_LIMB_MSB (sump[sn-1]) != 0); - MPFR_LOG_MSG (("Set exponent e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e)); - /* e may be outside the current exponent range, but this will be checked - with mpfr_check_range below. */ - MPFR_EXP (sum) = e; - } /* main block */ - - MPFR_TMP_FREE (marker); - return mpfr_check_range (sum, inex, rnd); -} - -/**********************************************************************/ - -int -mpfr_sum (mpfr_ptr sum, const mpfr_ptr *x, unsigned long n, mpfr_rnd_t rnd) -{ - MPFR_LOG_FUNC - (("n=%lu rnd=%d", n, rnd), - ("sum[%Pu]=%.*Rg", mpfr_get_prec (sum), mpfr_log_prec, sum)); - - if (MPFR_UNLIKELY (n <= 2)) - { - if (n == 0) - { - MPFR_SET_ZERO (sum); - MPFR_SET_POS (sum); - MPFR_RET (0); - } - else if (n == 1) - return mpfr_set (sum, x[0], rnd); - else - return mpfr_add (sum, x[0], x[1], rnd); - } - else - { - mpfr_exp_t maxexp = MPFR_EXP_MIN; /* max(Empty) */ - unsigned long i; - unsigned long rn = 0; /* will be the number of regular inputs */ - /* sign of infinities and zeros (0: currently unknown) */ - int sign_inf = 0, sign_zero = 0; - - MPFR_LOG_MSG (("Check for special inputs (n = %lu >= 3)\n", n)); - - for (i = 0; i < n; i++) - { - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x[i]))) - { - if (MPFR_IS_NAN (x[i])) - { - /* The current value x[i] is NaN. Then the sum is NaN. */ - nan: - MPFR_SET_NAN (sum); - MPFR_RET_NAN; - } - else if (MPFR_IS_INF (x[i])) - { - /* The current value x[i] is an infinity. - There are two cases: - 1. This is the first infinity value (sign_inf == 0). - Then set sign_inf to its sign, and go on. - 2. All the infinities found until now have the same - sign sign_inf. If this new infinity has a different - sign, then return NaN immediately, else go on. */ - if (sign_inf == 0) - sign_inf = MPFR_SIGN (x[i]); - else if (MPFR_SIGN (x[i]) != sign_inf) - goto nan; - } - else if (MPFR_UNLIKELY (rn == 0)) - { - /* The current value x[i] is a zero. The code below matters - only when all values found until now are zeros, otherwise - it is harmless (the test rn == 0 above is just a minor - optimization). - Here we track the sign of the zero result when all inputs - are zeros: if all zeros have the same sign, the result - will have this sign, otherwise (i.e. if there is at least - a zero of each sign), the sign of the zero result depends - only on the rounding mode (note that this choice is - sticky when new zeros are considered). */ - MPFR_ASSERTD (MPFR_IS_ZERO (x[i])); - if (sign_zero == 0) - sign_zero = MPFR_SIGN (x[i]); - else if (MPFR_SIGN (x[i]) != sign_zero) - sign_zero = rnd == MPFR_RNDD ? -1 : 1; - } - } - else - { - /* The current value x[i] is a regular number. */ - mpfr_exp_t e = MPFR_GET_EXP (x[i]); - if (e > maxexp) - maxexp = e; /* maximum exponent found until now */ - rn++; /* current number of regular inputs */ - } - } - - MPFR_LOG_MSG (("rn=%lu sign_inf=%d sign_zero=%d\n", - rn, sign_inf, sign_zero)); - - /* At this point the result cannot be NaN (this case has already - been filtered out). */ - - if (MPFR_UNLIKELY (sign_inf != 0)) - { - /* At least one infinity, and all of them have the same sign - sign_inf. The sum is the infinity of this sign. */ - MPFR_SET_INF (sum); - MPFR_SET_SIGN (sum, sign_inf); - MPFR_RET (0); - } - - /* At this point, all the inputs are finite numbers. */ - - if (MPFR_UNLIKELY (rn == 0)) - { - /* All the numbers were zeros (and there is at least one). - The sum is zero with sign sign_zero. */ - MPFR_ASSERTD (sign_zero != 0); - MPFR_SET_ZERO (sum); - MPFR_SET_SIGN (sum, sign_zero); - MPFR_RET (0); - } - - /* Optimize the case where there are only two regular numbers. */ - if (MPFR_UNLIKELY (rn <= 2)) - { - unsigned long h = ULONG_MAX; - - for (i = 0; i < n; i++) - if (! MPFR_IS_SINGULAR (x[i])) - { - if (rn == 1) - return mpfr_set (sum, x[i], rnd); - if (h != ULONG_MAX) - return mpfr_add (sum, x[h], x[i], rnd); - h = i; - } - MPFR_RET_NEVER_GO_HERE(); - } - - return sum_aux (sum, x, n, rnd, maxexp, rn); - } -} |