1 /* mpfr_round_raw_generic -- Generic rounding function
2
3 Copyright 1999-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 #ifndef flag
24 # error "ERROR: flag must be defined (0 / 1)"
25 #endif
26 #ifndef use_inexp
27 # error "ERROR: use_enexp must be defined (0 / 1)"
28 #endif
29 #ifndef mpfr_round_raw_generic
30 # error "ERROR: mpfr_round_raw_generic must be defined"
31 #endif
32
33 /*
34 * If flag = 0, puts in y the value of xp (with precision xprec and
35 * sign 1 if negative=0, -1 otherwise) rounded to precision yprec and
36 * direction rnd_mode. Supposes x is not zero nor NaN nor +/- Infinity
37 * (i.e. *xp != 0). In that case, the return value is a possible carry
38 * (0 or 1) that may happen during the rounding, in which case the result
39 * is a power of two.
40 *
41 * If inexp != NULL, put in *inexp the inexact flag of the rounding (0, 1, -1).
42 * In case of even rounding when rnd = MPFR_RNDN, put MPFR_EVEN_INEX (2) or
43 * -MPFR_EVEN_INEX (-2) in *inexp.
44 *
45 * If flag = 1, just returns whether one should add 1 or not for rounding.
46 *
47 * Note: yprec may be < MPFR_PREC_MIN; in particular, it may be equal
48 * to 1. In this case, the even rounding is done away from 0, which is
49 * a natural generalization. Indeed, a number with 1-bit precision can
50 * be seen as a subnormal number with more precision.
51 *
52 * MPFR_RNDNA is now supported, but needs to be tested [TODO] and is
53 * still not part of the API. In particular, the MPFR_RNDNA value (-1)
54 * may change in the future without notice.
55 */
56
57 #if !(flag == 0 || flag == 1)
58 #error "flag must be 0 or 1"
59 #endif
60
61 int
mpfr_round_raw_generic(mp_limb_t * yp,const mp_limb_t * xp,mpfr_prec_t xprec,int neg,mpfr_prec_t yprec,mpfr_rnd_t rnd_mode,int * inexp)62 mpfr_round_raw_generic(
63 #if flag == 0
64 mp_limb_t *yp,
65 #endif
66 const mp_limb_t *xp, mpfr_prec_t xprec,
67 int neg, mpfr_prec_t yprec, mpfr_rnd_t rnd_mode
68 #if use_inexp != 0
69 , int *inexp
70 #endif
71 )
72 {
73 mp_size_t xsize, nw;
74 mp_limb_t himask, lomask, sb;
75 int rw, new_use_inexp;
76 #if flag == 0
77 int carry;
78 #endif
79
80 #if use_inexp != 0
81 MPFR_ASSERTD (inexp != ((int*) 0));
82 #endif
83 MPFR_ASSERTD (neg == 0 || neg == 1);
84 #if flag == 1
85 /* rnd_mode = RNDF is only possible for flag = 0. */
86 MPFR_ASSERTD (rnd_mode != MPFR_RNDF);
87 #endif
88
89 if (rnd_mode == MPFR_RNDF)
90 {
91 #if use_inexp != 0
92 *inexp = 0; /* make sure it has a valid value */
93 #endif
94 rnd_mode = MPFR_RNDZ; /* faster */
95 new_use_inexp = 0;
96 }
97 else
98 {
99 if (flag && !use_inexp &&
100 (xprec <= yprec || MPFR_IS_LIKE_RNDZ (rnd_mode, neg)))
101 return 0;
102 new_use_inexp = use_inexp;
103 }
104
105 xsize = MPFR_PREC2LIMBS (xprec);
106 nw = yprec / GMP_NUMB_BITS;
107 rw = yprec & (GMP_NUMB_BITS - 1);
108
109 if (MPFR_UNLIKELY(xprec <= yprec))
110 { /* No rounding is necessary. */
111 /* if yp=xp, maybe an overlap: mpn_copyd is OK when src <= dst */
112 if (MPFR_LIKELY(rw))
113 nw++;
114 MPFR_ASSERTD(nw >= 1);
115 MPFR_ASSERTD(nw >= xsize);
116 #if use_inexp != 0
117 *inexp = 0;
118 #endif
119 #if flag == 0
120 mpn_copyd (yp + (nw - xsize), xp, xsize);
121 MPN_ZERO(yp, nw - xsize);
122 #endif
123 return 0;
124 }
125
126 if (new_use_inexp || !MPFR_IS_LIKE_RNDZ(rnd_mode, neg))
127 {
128 mp_size_t k = xsize - nw - 1;
129
130 if (MPFR_LIKELY(rw))
131 {
132 nw++;
133 lomask = MPFR_LIMB_MASK (GMP_NUMB_BITS - rw);
134 himask = ~lomask;
135 }
136 else
137 {
138 lomask = MPFR_LIMB_MAX;
139 himask = MPFR_LIMB_MAX;
140 }
141 MPFR_ASSERTD(k >= 0);
142 sb = xp[k] & lomask; /* First non-significant bits */
143 /* Rounding to nearest? */
144 if (rnd_mode == MPFR_RNDN || rnd_mode == MPFR_RNDNA)
145 {
146 /* Rounding to nearest */
147 mp_limb_t rbmask = MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1 - rw);
148
149 if ((sb & rbmask) == 0) /* rounding bit = 0 ? */
150 goto rnd_RNDZ; /* yes, behave like rounding toward zero */
151 /* Rounding to nearest with rounding bit = 1 */
152 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDNA))
153 goto away_addone_ulp; /* like rounding away from zero */
154 sb &= ~rbmask; /* first bits after the rounding bit */
155 while (MPFR_UNLIKELY(sb == 0) && k > 0)
156 sb = xp[--k];
157 if (MPFR_UNLIKELY(sb == 0)) /* Even rounding. */
158 {
159 /* sb == 0 && rnd_mode == MPFR_RNDN */
160 sb = xp[xsize - nw] & (himask ^ (himask << 1));
161 if (sb == 0)
162 {
163 #if use_inexp != 0
164 *inexp = 2 * MPFR_EVEN_INEX * neg - MPFR_EVEN_INEX;
165 #endif
166 /* ((neg!=0)^(sb!=0)) ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX */
167 /* since neg = 0 or 1 and sb = 0 */
168 #if flag == 0
169 mpn_copyi (yp, xp + xsize - nw, nw);
170 yp[0] &= himask;
171 #endif
172 return 0; /* sb != 0 && rnd_mode != MPFR_RNDZ */
173 }
174 else
175 {
176 away_addone_ulp:
177 /* sb != 0 && rnd_mode == MPFR_RNDN */
178 #if use_inexp != 0
179 *inexp = MPFR_EVEN_INEX - 2 * MPFR_EVEN_INEX * neg;
180 #endif
181 /* ((neg!=0)^(sb!=0)) ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX */
182 /* since neg = 0 or 1 and sb != 0 */
183 goto rnd_RNDN_add_one_ulp;
184 }
185 }
186 else /* sb != 0 && rnd_mode == MPFR_RNDN */
187 {
188 #if use_inexp != 0
189 *inexp = 1 - 2 * neg; /* neg == 0 ? 1 : -1 */
190 #endif
191 rnd_RNDN_add_one_ulp:
192 #if flag == 1
193 return 1; /* sb != 0 && rnd_mode != MPFR_RNDZ */
194 #else
195 carry = mpn_add_1 (yp, xp + xsize - nw, nw,
196 rw ?
197 MPFR_LIMB_ONE << (GMP_NUMB_BITS - rw)
198 : MPFR_LIMB_ONE);
199 yp[0] &= himask;
200 return carry;
201 #endif
202 }
203 }
204 /* Rounding toward zero? */
205 else if (MPFR_IS_LIKE_RNDZ(rnd_mode, neg))
206 {
207 /* rnd_mode == MPFR_RNDZ */
208 rnd_RNDZ:
209 while (MPFR_UNLIKELY (sb == 0) && k > 0)
210 sb = xp[--k];
211 #if use_inexp != 0
212 /* rnd_mode == MPFR_RNDZ and neg = 0 or 1 */
213 /* ((neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1 */
214 *inexp = MPFR_UNLIKELY (sb == 0) ? 0 : 2 * neg - 1;
215 #endif
216 #if flag == 0
217 mpn_copyi (yp, xp + xsize - nw, nw);
218 yp[0] &= himask;
219 #endif
220 return 0; /* sb != 0 && rnd_mode != MPFR_RNDZ */
221 }
222 else
223 {
224 /* Rounding away from zero */
225 while (MPFR_UNLIKELY (sb == 0) && k > 0)
226 sb = xp[--k];
227 if (MPFR_UNLIKELY (sb == 0))
228 {
229 /* sb = 0 && rnd_mode != MPFR_RNDZ */
230 #if use_inexp != 0
231 /* ((neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1 */
232 *inexp = 0;
233 #endif
234 #if flag == 0
235 mpn_copyi (yp, xp + xsize - nw, nw);
236 yp[0] &= himask;
237 #endif
238 return 0;
239 }
240 else
241 {
242 /* sb != 0 && rnd_mode != MPFR_RNDZ */
243 #if use_inexp != 0
244 *inexp = 1 - 2 * neg; /* neg == 0 ? 1 : -1 */
245 #endif
246 #if flag == 1
247 return 1;
248 #else
249 carry = mpn_add_1(yp, xp + xsize - nw, nw,
250 rw ? MPFR_LIMB_ONE << (GMP_NUMB_BITS - rw)
251 : MPFR_LIMB_ONE);
252 yp[0] &= himask;
253 return carry;
254 #endif
255 }
256 }
257 }
258 else
259 {
260 /* Rounding toward zero / no inexact flag */
261 #if flag == 0
262 if (MPFR_LIKELY(rw))
263 {
264 nw++;
265 himask = ~MPFR_LIMB_MASK (GMP_NUMB_BITS - rw);
266 }
267 else
268 himask = MPFR_LIMB_MAX;
269 mpn_copyi (yp, xp + xsize - nw, nw);
270 yp[0] &= himask;
271 #endif
272 return 0;
273 }
274 }
275
276 #undef flag
277 #undef use_inexp
278 #undef mpfr_round_raw_generic
279