1/* mpfr_round_nearest_away -- round to nearest away 2 3Copyright 2012-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-impl.h" 24 25/* Note: this doesn't work for 2^(emin-2). Currently, the best that can be 26 done is to extend the exponent range internally, and do not support the 27 case emin = MPFR_EMIN_MIN from the caller. */ 28 29/* In order to work, we'll save the context within the mantissa of 'rop'. 30 For that, we need to do some low level stuff like for init2/clear to create 31 a mantissa of bigger size, which contains the extra context. 32 To add a new field of type T in the context, add its type in 33 mpfr_size_limb_extended_t (if it is not already present) 34 and add a new value for the enum mpfr_index_extended_t before MANTISSA. */ 35typedef union { 36 mp_size_t si; 37 mp_limb_t li; 38 mpfr_exp_t ex; 39 mpfr_prec_t pr; 40 mpfr_sign_t sg; 41 mpfr_flags_t fl; 42 mpfr_limb_ptr pi; 43} mpfr_size_limb_extended_t; 44typedef enum { 45 ALLOC_SIZE = 0, 46 OLD_MANTISSA, 47 OLD_EXPONENT, 48 OLD_SIGN, 49 OLD_PREC, 50 OLD_FLAGS, 51 OLD_EXP_MIN, 52 OLD_EXP_MAX, 53 MANTISSA 54} mpfr_index_extended_t ; 55 56#define MPFR_MALLOC_EXTENDED_SIZE(s) \ 57 (MANTISSA * sizeof(mpfr_size_limb_extended_t) + \ 58 MPFR_BYTES_PER_MP_LIMB * (size_t) (s)) 59 60/* This function is called before the applied function 61 and prepares rop to give it one more bit of precision 62 and to save its old value within it. */ 63void 64mpfr_round_nearest_away_begin (mpfr_ptr rop) 65{ 66 mpfr_t tmp; 67 mp_size_t xsize; 68 mpfr_size_limb_extended_t *ext; 69 mpfr_prec_t p; 70 int inexact; 71 MPFR_SAVE_EXPO_DECL (expo); 72 73 /* when emin is the smallest possible value, we cannot 74 determine the correct round-to-nearest-away rounding for 75 0.25*2^emin, which gets rounded to 0 with nearest-even, 76 instead of 0.5^2^emin. */ 77 MPFR_ASSERTN(__gmpfr_emin > MPFR_EMIN_MIN); 78 79 p = MPFR_PREC (rop) + 1; 80 MPFR_SAVE_EXPO_MARK (expo); 81 82 /* We can't use mpfr_init2 since we need to store an additional context 83 within the mantissa. */ 84 MPFR_ASSERTN(p <= MPFR_PREC_MAX); 85 /* Allocate the context within the needed mantissa. */ 86 xsize = MPFR_PREC2LIMBS (p); 87 ext = (mpfr_size_limb_extended_t *) 88 mpfr_allocate_func (MPFR_MALLOC_EXTENDED_SIZE(xsize)); 89 90 /* Save the context first. */ 91 ext[ALLOC_SIZE].si = xsize; 92 ext[OLD_MANTISSA].pi = MPFR_MANT(rop); 93 ext[OLD_EXPONENT].ex = MPFR_EXP(rop); 94 ext[OLD_SIGN].sg = MPFR_SIGN(rop); 95 ext[OLD_PREC].pr = MPFR_PREC(rop); 96 ext[OLD_FLAGS].fl = expo.saved_flags; 97 ext[OLD_EXP_MIN].ex = expo.saved_emin; 98 ext[OLD_EXP_MAX].ex = expo.saved_emax; 99 100 /* Create tmp as a proper NAN. */ 101 MPFR_PREC(tmp) = p; /* Set prec */ 102 MPFR_SET_POS(tmp); /* Set a sign */ 103 MPFR_MANT(tmp) = (mp_limb_t*)(ext+MANTISSA); /* Set Mantissa ptr */ 104 MPFR_SET_NAN(tmp); /* initializes to NaN */ 105 106 /* Note: This full initialization to NaN may be unnecessary because of 107 the mpfr_set just below, but this is cleaner in case internals would 108 change in the future (for instance, some fields of tmp could be read 109 before being set, and reading an uninitialized value is UB here). */ 110 111 /* Copy rop into tmp now (rop may be used as input later). */ 112 MPFR_DBGRES (inexact = mpfr_set(tmp, rop, MPFR_RNDN)); 113 MPFR_ASSERTD(inexact == 0); /* Shall be exact since prec(tmp) > prec(rop) */ 114 115 /* Overwrite rop with tmp. */ 116 rop[0] = tmp[0]; 117 118 /* The new rop now contains a copy of the old rop, with one extra bit of 119 precision while keeping the old rop "hidden" within its mantissa. */ 120 121 return; 122} 123 124/* This function is called after the applied function. 125 rop is the one prepared in the function above, 126 and contains the result of the applied function. 127 This function restores rop like it was before the 128 calls to mpfr_round_nearest_away_begin while 129 copying it back the result of the applied function 130 and performing additional roundings. */ 131int 132mpfr_round_nearest_away_end (mpfr_ptr rop, int inex) 133{ 134 mpfr_t tmp; 135 mp_size_t xsize; 136 mpfr_size_limb_extended_t *ext; 137 mpfr_prec_t n; 138 MPFR_SAVE_EXPO_DECL (expo); 139 140 /* Get back the hidden context. 141 Note: The cast to void * prevents the compiler from emitting a warning 142 (or error), such as: 143 cast increases required alignment of target type 144 with the -Wcast-align GCC option. Correct alignment is a consequence 145 of the code that built rop. 146 */ 147 ext = ((mpfr_size_limb_extended_t *) (void *) MPFR_MANT(rop)) - MANTISSA; 148 149 /* Create tmp with the result of the function. */ 150 tmp[0] = rop[0]; 151 152 /* Revert rop back to what it was before calling 153 mpfr_round_neareset_away_begin. */ 154 MPFR_PREC(rop) = ext[OLD_PREC].pr; 155 MPFR_SIGN(rop) = ext[OLD_SIGN].sg; 156 MPFR_EXP(rop) = ext[OLD_EXPONENT].ex; 157 MPFR_MANT(rop) = ext[OLD_MANTISSA].pi; 158 MPFR_ASSERTD(MPFR_PREC(tmp) == MPFR_PREC(rop)+1); 159 160 /* Restore the saved context. */ 161 expo.saved_flags = ext[OLD_FLAGS].fl; 162 expo.saved_emin = ext[OLD_EXP_MIN].ex; 163 expo.saved_emax = ext[OLD_EXP_MAX].ex; 164 xsize = ext[ALLOC_SIZE].si; 165 166 /* Perform RNDNA. */ 167 n = MPFR_PREC(rop); 168 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp))) 169 mpfr_set (rop, tmp, MPFR_RNDN); /* inex unchanged */ 170 else 171 { 172 int lastbit, sh; 173 174 MPFR_UNSIGNED_MINUS_MODULO(sh, n + 1); 175 lastbit = (MPFR_MANT(tmp)[0] >> sh) & 1; 176 177 if (lastbit == 0) 178 mpfr_set (rop, tmp, MPFR_RNDN); /* exact, inex unchanged */ 179 else if (inex == 0) /* midpoint: round away from zero */ 180 inex = mpfr_set (rop, tmp, MPFR_RNDA); 181 else /* lastbit == 1, inex != 0: double rounding */ 182 inex = mpfr_set (rop, tmp, (inex > 0) ? MPFR_RNDD : MPFR_RNDU); 183 } 184 185 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 186 MPFR_SAVE_EXPO_FREE (expo); 187 188 /* special treatment for the case rop = +/- 2^(emin-2), which should be 189 rounded to +/- 2^(emin-1). We do as if it was rounded to zero, thus the 190 mpfr_check_range() call will round it to +/- 2^(emin-1). */ 191 if (inex == 0 && mpfr_cmp_si_2exp (rop, (mpfr_sgn (rop) > 0) ? 1 : -1, 192 __gmpfr_emin - 2) == 0) 193 inex = -mpfr_sgn (rop); 194 195 /* Free tmp (cannot call mpfr_clear): free the associated context. */ 196 mpfr_free_func(ext, MPFR_MALLOC_EXTENDED_SIZE(xsize)); 197 198 return mpfr_check_range (rop, inex, MPFR_RNDN); 199} 200