summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/sum.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sum.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sum.c1395
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);
- }
-}