diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/rndna.c')
| -rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/rndna.c | 199 |
1 files changed, 199 insertions, 0 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/rndna.c b/Build/source/libs/mpfr/mpfr-src/src/rndna.c new file mode 100644 index 00000000000..23997bf0ee7 --- /dev/null +++ b/Build/source/libs/mpfr/mpfr-src/src/rndna.c @@ -0,0 +1,199 @@ +/* mpfr_round_nearest_away -- round to nearest away + +Copyright 2012-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" + +/* 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. + Note: The cast to void * prevents the compiler from emitting a warning + (or error), such as: + cast increases required alignment of target type + with the -Wcast-align GCC option. Correct alignment is a consequence + of the code that built rop. + */ + ext = ((mpfr_size_limb_extended_t *) (void *) 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); +} |