1 /* mpfr_round_nearest_away -- round to nearest away
2
3 Copyright 2012-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 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. */
35 typedef 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;
44 typedef 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. */
63 void
64 mpfr_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. */
131 int
132 mpfr_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 }