/
niceSOFT
/
mpfr
Обзор
Документация
Войти
/
niceSOFT
/
mpfr
Код
Задачи
Вики
Пакеты
0
Релизы
0
Аналитика
Безопасность
emptycommit-12104
src/round_prec.c
532 строки
20 KB
Vincent Lefèvre
Copyright notice update: added 2018 with
03 янв 2018, 13:53
03 янв 2018, 13:53
2fce388
Код
Авторство
О чём код?
/* 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-2018 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 http://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" #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); /* Transform RNDD and RNDU to Zero / Away */ MPFR_ASSERTD (neg == 0 || neg == 1); /* transform RNDF to RNDN */ if (rnd1 == MPFR_RNDF) rnd1 = MPFR_RNDN; if (rnd1 != MPFR_RNDN) rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA; if (rnd2 == MPFR_RNDF) rnd2 = (rnd1 == MPFR_RNDN) ? MPFR_RNDN : MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDA : MPFR_RNDZ; if (rnd2 != MPFR_RNDN) rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : 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) { if (rnd1 == MPFR_RNDN && prec < (mpfr_prec_t) bn * GMP_NUMB_BITS) { /* then rnd2 = RNDN, and for prec = bn * GMP_NUMB_BITS we cannot have b the middle of two representable numbers */ 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 representable in precision prec+1 and ends with a 1 */ return 0; else return 1; } else 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); MPFR_ASSERTD (k > 0); switch (rnd1) { case MPFR_RNDZ: /* rnd1 = Round to Zero */ cc = (bp[bn - 1] >> s1) & 1; /* 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] == 0; 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] == 0; 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; }