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