/* mpfr_round_nearest_away -- round to nearest away Copyright 2012-2017 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" /* Note: this doesn't work for 2^(emin-2). Currently, the best that can be done is to extend the exponent range internally, and do not support the case emin = MPFR_EMIN_MIN from the caller. */ /* In order to work, we'll save the context within the mantissa of 'rop'. For that, we need to do some low level stuff like for init2/clear to create a mantissa of bigger size, which contains the extra context. To add a new field of type T in the context, add its type in mpfr_size_limb_extended_t (if it is not already present) and add a new value for the enum mpfr_index_extended_t before MANTISSA. */ typedef union { mp_size_t si; mp_limb_t li; mpfr_exp_t ex; mpfr_prec_t pr; mpfr_sign_t sg; mpfr_flags_t fl; mpfr_limb_ptr pi; } mpfr_size_limb_extended_t; typedef enum { ALLOC_SIZE = 0, OLD_MANTISSA, OLD_EXPONENT, OLD_SIGN, OLD_PREC, OLD_FLAGS, OLD_EXP_MIN, OLD_EXP_MAX, MANTISSA } mpfr_index_extended_t ; #define MPFR_MALLOC_EXTENDED_SIZE(s) \ (MANTISSA * sizeof(mpfr_size_limb_extended_t) + \ MPFR_BYTES_PER_MP_LIMB * (size_t) (s)) /* This function is called before the applied function and prepares rop to give it one more bit of precision and to save its old value within it. */ void mpfr_round_nearest_away_begin (mpfr_t rop) { mpfr_t tmp; mp_size_t xsize; mpfr_size_limb_extended_t *ext; mpfr_prec_t p; int inexact; MPFR_SAVE_EXPO_DECL (expo); /* when emin is the smallest possible value, we cannot determine the correct round-to-nearest-away rounding for 0.25*2^emin, which gets rounded to 0 with nearest-even, instead of 0.5^2^emin. */ MPFR_ASSERTN(__gmpfr_emin > MPFR_EMIN_MIN); p = MPFR_PREC (rop) + 1; MPFR_SAVE_EXPO_MARK (expo); /* We can't use mpfr_init2 since we need to store an additional context within the mantissa. */ MPFR_ASSERTN(p <= MPFR_PREC_MAX); /* Allocate the context within the needed mantissa. */ xsize = MPFR_PREC2LIMBS (p); ext = (mpfr_size_limb_extended_t *) mpfr_allocate_func (MPFR_MALLOC_EXTENDED_SIZE(xsize)); /* Save the context first. */ ext[ALLOC_SIZE].si = xsize; ext[OLD_MANTISSA].pi = MPFR_MANT(rop); ext[OLD_EXPONENT].ex = MPFR_EXP(rop); ext[OLD_SIGN].sg = MPFR_SIGN(rop); ext[OLD_PREC].pr = MPFR_PREC(rop); ext[OLD_FLAGS].fl = expo.saved_flags; ext[OLD_EXP_MIN].ex = expo.saved_emin; ext[OLD_EXP_MAX].ex = expo.saved_emax; /* Create tmp as a proper NAN. */ MPFR_PREC(tmp) = p; /* Set prec */ MPFR_SET_POS(tmp); /* Set a sign */ MPFR_MANT(tmp) = (mp_limb_t*)(ext+MANTISSA); /* Set Mantissa ptr */ MPFR_SET_NAN(tmp); /* initializes to NaN */ /* Note: This full initialization to NaN may be unnecessary because of the mpfr_set just below, but this is cleaner in case internals would change in the future (for instance, some fields of tmp could be read before being set, and reading an uninitialized value is UB here). */ /* Copy rop into tmp now (rop may be used as input later). */ MPFR_DBGRES (inexact = mpfr_set(tmp, rop, MPFR_RNDN)); MPFR_ASSERTD(inexact == 0); /* Shall be exact since prec(tmp) > prec(rop) */ /* Overwrite rop with tmp. */ rop[0] = tmp[0]; /* The new rop now contains a copy of the old rop, with one extra bit of precision while keeping the old rop "hidden" within its mantissa. */ return; } /* This function is called after the applied function. rop is the one prepared in the function above, and contains the result of the applied function. This function restores rop like it was before the calls to mpfr_round_nearest_away_begin while copying it back the result of the applied function and performing additional roundings. */ int mpfr_round_nearest_away_end (mpfr_t rop, int inex) { mpfr_t tmp; mp_size_t xsize; mpfr_size_limb_extended_t *ext; mpfr_prec_t n; MPFR_SAVE_EXPO_DECL (expo); /* Get back the hidden context. */ ext = ((mpfr_size_limb_extended_t *) MPFR_MANT(rop)) - MANTISSA; /* Create tmp with the result of the function. */ tmp[0] = rop[0]; /* Revert rop back to what it was before calling mpfr_round_neareset_away_begin. */ MPFR_PREC(rop) = ext[OLD_PREC].pr; MPFR_SIGN(rop) = ext[OLD_SIGN].sg; MPFR_EXP(rop) = ext[OLD_EXPONENT].ex; MPFR_MANT(rop) = ext[OLD_MANTISSA].pi; MPFR_ASSERTD(MPFR_PREC(tmp) == MPFR_PREC(rop)+1); /* Restore the saved context. */ expo.saved_flags = ext[OLD_FLAGS].fl; expo.saved_emin = ext[OLD_EXP_MIN].ex; expo.saved_emax = ext[OLD_EXP_MAX].ex; xsize = ext[ALLOC_SIZE].si; /* Perform RNDNA. */ n = MPFR_PREC(rop); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp))) mpfr_set (rop, tmp, MPFR_RNDN); /* inex unchanged */ else { int lastbit, sh; MPFR_UNSIGNED_MINUS_MODULO(sh, n + 1); lastbit = (MPFR_MANT(tmp)[0] >> sh) & 1; if (lastbit == 0) mpfr_set (rop, tmp, MPFR_RNDN); /* exact, inex unchanged */ else if (inex == 0) /* midpoint: round away from zero */ inex = mpfr_set (rop, tmp, MPFR_RNDA); else /* lastbit == 1, inex != 0: double rounding */ inex = mpfr_set (rop, tmp, (inex > 0) ? MPFR_RNDD : MPFR_RNDU); } MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); MPFR_SAVE_EXPO_FREE (expo); /* special treatment for the case rop = +/- 2^(emin-2), which should be rounded to +/- 2^(emin-1). We do as if it was rounded to zero, thus the mpfr_check_range() call will round it to +/- 2^(emin-1). */ if (inex == 0 && mpfr_cmp_si_2exp (rop, (mpfr_sgn (rop) > 0) ? 1 : -1, __gmpfr_emin - 2) == 0) inex = -mpfr_sgn (rop); /* Free tmp (cannot call mpfr_clear): free the associated context. */ mpfr_free_func(ext, MPFR_MALLOC_EXTENDED_SIZE(xsize)); return mpfr_check_range (rop, inex, MPFR_RNDN); }