diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/round_prec.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/round_prec.c | 539 |
1 files changed, 0 insertions, 539 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/round_prec.c b/Build/source/libs/mpfr/mpfr-src/src/round_prec.c deleted file mode 100644 index f03fbac7be7..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/round_prec.c +++ /dev/null @@ -1,539 +0,0 @@ -/* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_prec_round, - mpfr_can_round, mpfr_can_round_raw -- various rounding functions - -Copyright 1999-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. */ - -#include "mpfr-impl.h" - -#define mpfr_round_raw_generic mpfr_round_raw -#define flag 0 -#define use_inexp 1 -#include "round_raw_generic.c" - -/* mpfr_round_raw_2 is called from mpfr_round_raw2 */ -#define mpfr_round_raw_generic mpfr_round_raw_2 -#define flag 1 -#define use_inexp 0 -#include "round_raw_generic.c" - -/* Seems to be unused. Remove comment to implement it. -#define mpfr_round_raw_generic mpfr_round_raw_3 -#define flag 1 -#define use_inexp 1 -#include "round_raw_generic.c" -*/ - -#define mpfr_round_raw_generic mpfr_round_raw_4 -#define flag 0 -#define use_inexp 0 -#include "round_raw_generic.c" - -/* Note: if the new prec is lower than the current one, a reallocation - must not be done (see exp_2.c). */ - -int -mpfr_prec_round (mpfr_ptr x, mpfr_prec_t prec, mpfr_rnd_t rnd_mode) -{ - mp_limb_t *tmp, *xp; - int carry, inexact; - mpfr_prec_t nw, ow; - MPFR_TMP_DECL(marker); - - MPFR_ASSERTN (MPFR_PREC_COND (prec)); - - nw = MPFR_PREC2LIMBS (prec); /* needed allocated limbs */ - - /* check if x has enough allocated space for the significand */ - /* Get the number of limbs from the precision. - (Compatible with all allocation methods) */ - ow = MPFR_LIMB_SIZE (x); - if (MPFR_UNLIKELY (nw > ow)) - { - /* FIXME: Variable can't be created using custom allocation, - MPFR_DECL_INIT or GROUP_ALLOC: How to detect? */ - ow = MPFR_GET_ALLOC_SIZE(x); - if (nw > ow) - { - mpfr_size_limb_t *tmpx; - - /* Realloc significand */ - tmpx = (mpfr_size_limb_t *) mpfr_reallocate_func - (MPFR_GET_REAL_PTR(x), MPFR_MALLOC_SIZE(ow), MPFR_MALLOC_SIZE(nw)); - MPFR_SET_MANT_PTR(x, tmpx); /* mant ptr must be set - before alloc size */ - MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */ - } - } - - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) - { - MPFR_PREC(x) = prec; /* Special value: need to set prec */ - if (MPFR_IS_NAN(x)) - MPFR_RET_NAN; - MPFR_ASSERTD(MPFR_IS_INF(x) || MPFR_IS_ZERO(x)); - return 0; /* infinity and zero are exact */ - } - - /* x is a non-zero real number */ - - MPFR_TMP_MARK(marker); - tmp = MPFR_TMP_LIMBS_ALLOC (nw); - xp = MPFR_MANT(x); - carry = mpfr_round_raw (tmp, xp, MPFR_PREC(x), MPFR_IS_NEG(x), - prec, rnd_mode, &inexact); - MPFR_PREC(x) = prec; - - if (MPFR_UNLIKELY(carry)) - { - mpfr_exp_t exp = MPFR_EXP (x); - - if (MPFR_UNLIKELY(exp == __gmpfr_emax)) - (void) mpfr_overflow(x, rnd_mode, MPFR_SIGN(x)); - else - { - MPFR_ASSERTD (exp < __gmpfr_emax); - MPFR_SET_EXP (x, exp + 1); - xp[nw - 1] = MPFR_LIMB_HIGHBIT; - if (nw - 1 > 0) - MPN_ZERO(xp, nw - 1); - } - } - else - MPN_COPY(xp, tmp, nw); - - MPFR_TMP_FREE(marker); - return inexact; -} - -/* assumption: GMP_NUMB_BITS is a power of 2 */ - -/* assuming b is an approximation to x in direction rnd1 with error at - most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly - x to precision prec with direction rnd2, and 0 otherwise. - Side effects: none. - - rnd1 = RNDN and RNDF are similar: the sign of the error is unknown. - - rnd2 = RNDF: assume that the user will round the approximation b - toward the direction of x, i.e. the opposite of rnd1 in directed - rounding modes, otherwise RNDN. Some details: - - u xinf v xsup w - -----|----+----------|--+------------|----- - [----- x -----] - rnd1 = RNDD b | - rnd1 = RNDU b - - where u, v and w are consecutive machine numbers. - - * If [xinf,xsup] contains no machine numbers, then return 1. - - * If [xinf,xsup] contains 2 machine numbers, then return 0. - - * If [xinf,xsup] contains a single machine number, then return 1 iff - the rounding of b is this machine number. - With the above choice for the rounding of b, this will always be - the case if rnd1 is a directed rounding mode; said otherwise, for - rnd2 = RNDF and rnd1 being a directed rounding mode, return 1 iff - [xinf,xsup] contains at most 1 machine number. -*/ - -int -mpfr_can_round (mpfr_srcptr b, mpfr_exp_t err, mpfr_rnd_t rnd1, - mpfr_rnd_t rnd2, mpfr_prec_t prec) -{ - if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b))) - return 0; /* We cannot round if Zero, Nan or Inf */ - else - return mpfr_can_round_raw (MPFR_MANT(b), MPFR_LIMB_SIZE(b), - MPFR_SIGN(b), err, rnd1, rnd2, prec); -} - -/* TODO: mpfr_can_round_raw currently does a memory allocation and some - mpn operations. A bit inspection like for mpfr_round_p (round_p.c) may - be sufficient, though this would be more complex than the one done in - mpfr_round_p, and in particular, for some rnd1/rnd2 combinations, one - needs to take care of changes of binade when the value is close to a - power of 2. */ - -int -mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err, - mpfr_rnd_t rnd1, mpfr_rnd_t rnd2, mpfr_prec_t prec) -{ - mpfr_prec_t prec2; - mp_size_t k, k1, tn; - int s, s1; - mp_limb_t cc, cc2; - mp_limb_t *tmp; - mp_limb_t cy = 0, tmp_hi; - int res; - MPFR_TMP_DECL(marker); - - /* Since mpfr_can_round is a function in the API, use MPFR_ASSERTN. - The specification makes sense only for prec >= 1. */ - MPFR_ASSERTN (prec >= 1); - - MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT); - - MPFR_ASSERT_SIGN(neg); - neg = MPFR_IS_NEG_SIGN(neg); - MPFR_ASSERTD (neg == 0 || neg == 1); - - /* For rnd1 and rnd2, transform RNDF / RNDD / RNDU to RNDN / RNDZ / RNDA - (with a special case for rnd1 directed rounding, rnd2 = RNDF). */ - - if (rnd1 == MPFR_RNDF) - rnd1 = MPFR_RNDN; /* transform RNDF to RNDN */ - else if (rnd1 != MPFR_RNDN) - rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA; - - MPFR_ASSERTD (rnd1 == MPFR_RNDN || - rnd1 == MPFR_RNDZ || - rnd1 == MPFR_RNDA); - - if (rnd2 == MPFR_RNDF) - { - if (rnd1 == MPFR_RNDN) - rnd2 = MPFR_RNDN; - else - { - rnd2 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDA : MPFR_RNDZ; - /* Warning: in this case (rnd1 directed rounding, rnd2 = RNDF), - the specification of mpfr_can_round says that we should - return non-zero (i.e., we can round) when {bp, bn} is - exactly representable in precision prec. */ - if (mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) == 0) - return 1; - } - } - else if (rnd2 != MPFR_RNDN) - rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : MPFR_RNDA; - - MPFR_ASSERTD (rnd2 == MPFR_RNDN || - rnd2 == MPFR_RNDZ || - rnd2 == MPFR_RNDA); - - /* For err < prec (+1 for rnd1=RNDN), we can never round correctly, since - the error is at least 2*ulp(b) >= ulp(round(b)). - However for err = prec (+1 for rnd1=RNDN), we can round correctly in some - rare cases where ulp(b) = 1/2*ulp(U) [see below for the definition of U], - which implies rnd1 = RNDZ or RNDN, and rnd2 = RNDA or RNDN. */ - - if (MPFR_UNLIKELY (err < prec + (rnd1 == MPFR_RNDN) || - (err == prec + (rnd1 == MPFR_RNDN) && - (rnd1 == MPFR_RNDA || - rnd2 == MPFR_RNDZ)))) - return 0; /* can't round */ - - /* As a consequence... */ - MPFR_ASSERTD (err >= prec); - - /* The bound c on the error |x-b| is: c = 2^(MPFR_EXP(b)-err) <= b/2. - * So, we now know that x and b have the same sign. By symmetry, - * assume x > 0 and b > 0. We have: L <= x <= U, where, depending - * on rnd1: - * MPFR_RNDN: L = b-c, U = b+c - * MPFR_RNDZ: L = b, U = b+c - * MPFR_RNDA: L = b-c, U = b - * - * We can round x iff round(L,prec,rnd2) = round(U,prec,rnd2). - */ - - if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS)) - { /* Then prec > PREC(b): we can round: - (i) in rounding to the nearest as long as err >= prec + 2. - When err = prec + 1 and b is not a power - of two (so that a change of binade cannot occur), then one - can round to nearest thanks to the even rounding rule (in the - target precision prec, the significand of b ends with a 0). - When err = prec + 1 and b is a power of two, when rnd1 = RNDZ one - can round too. - (ii) in directed rounding mode iff rnd1 is compatible with rnd2 - and err >= prec + 1, unless b = 2^k and rnd1 = RNDA or RNDN in - which case we need err >= prec + 2. - */ - if ((rnd1 == rnd2 || rnd2 == MPFR_RNDN) && err >= prec + 1) - { - if (rnd1 != MPFR_RNDZ && - err == prec + 1 && - mpfr_powerof2_raw2 (bp, bn)) - return 0; - else - return 1; - } - return 0; - } - - /* now prec <= bn * GMP_NUMB_BITS */ - - if (MPFR_UNLIKELY (err > (mpfr_prec_t) bn * GMP_NUMB_BITS)) - { - /* we distinguish the case where b is a power of two: - rnd1 rnd2 can round? - RNDZ RNDZ ok - RNDZ RNDA no - RNDZ RNDN ok - RNDA RNDZ no - RNDA RNDA ok except when err = prec + 1 - RNDA RNDN ok except when err = prec + 1 - RNDN RNDZ no - RNDN RNDA no - RNDN RNDN ok except when err = prec + 1 */ - if (mpfr_powerof2_raw2 (bp, bn)) - { - if ((rnd2 == MPFR_RNDZ || rnd2 == MPFR_RNDA) && rnd1 != rnd2) - return 0; - else if (rnd1 == MPFR_RNDZ) - return 1; /* RNDZ RNDZ and RNDZ RNDN */ - else - return err > prec + 1; - } - - /* now the general case where b is not a power of two: - rnd1 rnd2 can round? - RNDZ RNDZ ok - RNDZ RNDA except when b is representable in precision 'prec' - RNDZ RNDN except when b is the middle of two representable numbers in - precision 'prec' and b ends with 'xxx0[1]', - or b is representable in precision 'prec' - and err = prec + 1 and b ends with '1'. - RNDA RNDZ except when b is representable in precision 'prec' - RNDA RNDA ok - RNDA RNDN except when b is the middle of two representable numbers in - precision 'prec' and b ends with 'xxx1[1]', - or b is representable in precision 'prec' - and err = prec + 1 and b ends with '1'. - RNDN RNDZ except when b is representable in precision 'prec' - RNDN RNDA except when b is representable in precision 'prec' - RNDN RNDN except when b is the middle of two representable numbers in - precision 'prec', or b is representable in precision 'prec' - and err = prec + 1 and b ends with '1'. */ - if (rnd2 == MPFR_RNDN) - { - if (err == prec + 1 && (bp[0] & 1)) - return 0; /* err == prec + 1 implies prec = bn * GMP_NUMB_BITS */ - if (prec < (mpfr_prec_t) bn * GMP_NUMB_BITS) - { - k1 = MPFR_PREC2LIMBS (prec + 1); - MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1); - if (((bp[bn - k1] >> s1) & 1) && - mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0) - { /* b is the middle of two representable numbers */ - if (rnd1 == MPFR_RNDN) - return 0; - k1 = MPFR_PREC2LIMBS (prec); - MPFR_UNSIGNED_MINUS_MODULO(s1, prec); - return (rnd1 == MPFR_RNDZ) ^ - (((bp[bn - k1] >> s1) & 1) == 0); - } - } - return 1; - } - else if (rnd1 == rnd2) /* cases RNDZ RNDZ or RNDA RNDA: ok */ - return 1; - else - return mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) != 0; - } - - /* now err <= bn * GMP_NUMB_BITS */ - - /* warning: if k = m*GMP_NUMB_BITS, consider limb m-1 and not m */ - k = (err - 1) / GMP_NUMB_BITS; - MPFR_UNSIGNED_MINUS_MODULO(s, err); - /* the error corresponds to bit s in limb k, the most significant limb - being limb 0; in memory, limb k is bp[bn-1-k]. */ - - k1 = (prec - 1) / GMP_NUMB_BITS; - MPFR_UNSIGNED_MINUS_MODULO(s1, prec); - /* the least significant bit is bit s1 in limb k1 */ - - /* We don't need to consider the k1 most significant limbs. - They will be considered later only to detect when subtracting - the error bound yields a change of binade. - Warning! The number with updated bn may no longer be normalized. */ - k -= k1; - bn -= k1; - prec2 = prec - (mpfr_prec_t) k1 * GMP_NUMB_BITS; - - /* We can decide of the correct rounding if rnd2(b-eps) and rnd2(b+eps) - give the same result to the target precision 'prec', i.e., if when - adding or subtracting (1 << s) in bp[bn-1-k], it does not change the - rounding in direction 'rnd2' at ulp-position bp[bn-1] >> s1, taking also - into account the possible change of binade. */ - MPFR_TMP_MARK(marker); - tn = bn; - k++; /* since we work with k+1 everywhere */ - tmp = MPFR_TMP_LIMBS_ALLOC (tn); - if (bn > k) - MPN_COPY (tmp, bp, bn - k); /* copy low bn-k limbs of b into tmp */ - - MPFR_ASSERTD (k > 0); - - switch (rnd1) - { - case MPFR_RNDZ: - /* rnd1 = Round to Zero */ - cc = (bp[bn - 1] >> s1) & 1; /* cc is the least significant bit of b */ - /* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec), - and 0 otherwise */ - cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2); - /* cc is the new value of bit s1 in bp[bn-1] after rounding 'rnd2' */ - - /* now round b + 2^(MPFR_EXP(b)-err) */ - cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s); - /* propagate carry up to most significant limb */ - for (tn = 0; tn + 1 < k1 && cy != 0; tn ++) - cy = bp[bn + tn] == MPFR_LIMB_MAX; - if (cy == 0 && err == prec) - { - res = 0; - goto end; - } - if (MPFR_UNLIKELY(cy)) - { - /* when a carry occurs, we have b < 2^h <= b+c, we can round iff: - rnd2 = RNDZ: never, since b and b+c round to different values; - rnd2 = RNDA: when b+c is an exact power of two, and err > prec - (since for err = prec, b = 2^h - 1/2*ulp(2^h) is - exactly representable and thus rounds to itself); - rnd2 = RNDN: whenever cc = 0, since err >= prec implies - c <= ulp(b) = 1/2*ulp(2^h), thus b+c rounds to 2^h, - and b+c >= 2^h implies that bit 'prec' of b is 1, - thus cc = 0 means that b is rounded to 2^h too. */ - res = (rnd2 == MPFR_RNDZ) ? 0 - : (rnd2 == MPFR_RNDA) ? (err > prec && k == bn && tmp[0] == 0) - : cc == 0; - goto end; - } - break; - case MPFR_RNDN: - /* rnd1 = Round to nearest */ - - /* first round b+2^(MPFR_EXP(b)-err) */ - cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s); - /* propagate carry up to most significant limb */ - for (tn = 0; tn + 1 < k1 && cy != 0; tn ++) - cy = bp[bn + tn] == MPFR_LIMB_MAX; - cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */ - cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2); - /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */ - if (MPFR_UNLIKELY (cy != 0)) - { - /* when a carry occurs, we have b-c < b < 2^h <= b+c, we can round - iff: - rnd2 = RNDZ: never, since b-c and b+c round to different values; - rnd2 = RNDA: when b+c is an exact power of two, and - err > prec + 1 (since for err <= prec + 1, - b-c <= 2^h - 1/2*ulp(2^h) is exactly representable - and thus rounds to itself); - rnd2 = RNDN: whenever err > prec + 1, since for err = prec + 1, - b+c rounds to 2^h, and b-c rounds to nextbelow(2^h). - For err > prec + 1, c <= 1/4*ulp(b) <= 1/8*ulp(2^h), - thus - 2^h - 1/4*ulp(b) <= b-c < b+c <= 2^h + 1/8*ulp(2^h), - therefore both b-c and b+c round to 2^h. */ - res = (rnd2 == MPFR_RNDZ) ? 0 - : (rnd2 == MPFR_RNDA) ? (err > prec + 1 && k == bn && tmp[0] == 0) - : err > prec + 1; - goto end; - } - subtract_eps: - /* now round b-2^(MPFR_EXP(b)-err), this happens for - rnd1 = RNDN or RNDA */ - MPFR_ASSERTD(rnd1 == MPFR_RNDN || rnd1 == MPFR_RNDA); - cy = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s); - /* propagate the potential borrow up to the most significant limb - (it cannot propagate further since the most significant limb is - at least MPFR_LIMB_HIGHBIT). - Note: we use the same limb tmp[bn-1] to subtract. */ - tmp_hi = tmp[bn - 1]; - for (tn = 0; tn < k1 && cy != 0; tn ++) - cy = mpn_sub_1 (&tmp_hi, bp + bn + tn, 1, cy); - /* We have an exponent decrease when tn = k1 and - tmp[bn-1] < MPFR_LIMB_HIGHBIT: - b-c < 2^h <= b (for RNDA) or b+c (for RNDN). - Then we surely cannot round when rnd2 = RNDZ, since b or b+c round to - a value >= 2^h, and b-c rounds to a value < 2^h. - We also surely cannot round when (rnd1,rnd2) = (RNDN,RNDA), since - b-c rounds to a value <= 2^h, and b+c > 2^h rounds to a value > 2^h. - It thus remains: - (rnd1,rnd2) = (RNDA,RNDA), (RNDA,RNDN) and (RNDN,RNDN). - For (RNDA,RNDA) we can round only when b-c and b round to 2^h, which - implies b = 2^h and err > prec (which is true in that case): - a necessary condition is that cc = 0. - For (RNDA,RNDN) we can round only when b-c and b round to 2^h, which - implies b-c >= 2^h - 1/4*ulp(2^h), and b <= 2^h + 1/2*ulp(2^h); - since ulp(2^h) = ulp(b), this implies c <= 3/4*ulp(b), thus - err > prec. - For (RNDN,RNDN) we can round only when b-c and b+c round to 2^h, - which implies b-c >= 2^h - 1/4*ulp(2^h), and - b+c <= 2^h + 1/2*ulp(2^h); - since ulp(2^h) = ulp(b), this implies 2*c <= 3/4*ulp(b), thus - err > prec+1. - */ - if (tn == k1 && tmp_hi < MPFR_LIMB_HIGHBIT) /* exponent decrease */ - { - if (rnd2 == MPFR_RNDZ || (rnd1 == MPFR_RNDN && rnd2 == MPFR_RNDA) || - cc != 0 /* b or b+c does not round to 2^h */) - { - res = 0; - goto end; - } - /* in that case since the most significant bit of tmp is 0, we - should consider one more bit; res = 0 when b-c does not round - to 2^h. */ - res = mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2 + 1) != 0; - goto end; - } - if (err == prec + (rnd1 == MPFR_RNDN)) - { - /* No exponent increase nor decrease, thus we have |U-L| = ulp(b). - For rnd2 = RNDZ or RNDA, either [L,U] contains one representable - number in the target precision, and then L and U round - differently; or both L and U are representable: they round - differently too; thus in all cases we cannot round. - For rnd2 = RNDN, the only case where we can round is when the - middle of [L,U] (i.e. b) is representable, and ends with a 0. */ - res = (rnd2 == MPFR_RNDN && (((bp[bn - 1] >> s1) & 1) == 0) && - mpfr_round_raw2 (bp, bn, neg, MPFR_RNDZ, prec2) == - mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec2)); - goto end; - } - break; - default: - /* rnd1 = Round away */ - MPFR_ASSERTD (rnd1 == MPFR_RNDA); - cc = (bp[bn - 1] >> s1) & 1; - /* the mpfr_round_raw2() call below returns whether one should add 1 or - not for rounding */ - cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2); - /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */ - - goto subtract_eps; - } - - cc2 = (tmp[bn - 1] >> s1) & 1; - res = cc == (cc2 ^ mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2)); - - end: - MPFR_TMP_FREE(marker); - return res; -} |