1 /* mpfr_div_ui -- divide a floating-point number by a machine integer
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 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 #ifdef MPFR_COV_CHECK
27 int __gmpfr_cov_div_ui_sb[10][2] = { 0 };
28 #endif
29
30 /* returns 0 if result exact, non-zero otherwise */
31 #undef mpfr_div_ui
32 MPFR_HOT_FUNCTION_ATTR int
33 mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u,
34 mpfr_rnd_t rnd_mode)
35 {
36 int inexact;
37
38 #ifdef MPFR_LONG_WITHIN_LIMB
39
40 int sh;
41 mp_size_t i, xn, yn, dif;
42 mp_limb_t *xp, *yp, *tmp, c, d;
43 mpfr_exp_t exp;
44 mp_limb_t rb; /* round bit */
45 mp_limb_t sb; /* sticky bit */
46 MPFR_TMP_DECL(marker);
47
48 MPFR_LOG_FUNC
49 (("x[%Pd]=%.*Rg u=%lu rnd=%d",
50 mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
51 ("y[%Pd]=%.*Rg inexact=%d",
52 mpfr_get_prec(y), mpfr_log_prec, y, inexact));
53
54 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
55 {
56 if (MPFR_IS_NAN (x))
57 {
58 MPFR_SET_NAN (y);
59 MPFR_RET_NAN;
60 }
61 else if (MPFR_IS_INF (x))
62 {
63 MPFR_SET_INF (y);
64 MPFR_SET_SAME_SIGN (y, x);
65 MPFR_RET (0);
66 }
67 else
68 {
69 MPFR_ASSERTD (MPFR_IS_ZERO (x));
70 if (u == 0) /* 0/0 is NaN */
71 {
72 MPFR_SET_NAN (y);
73 MPFR_RET_NAN;
74 }
75 else
76 {
77 MPFR_SET_ZERO(y);
78 MPFR_SET_SAME_SIGN (y, x);
79 MPFR_RET(0);
80 }
81 }
82 }
83 else if (MPFR_UNLIKELY (u <= 1))
84 {
85 if (u < 1)
86 {
87 /* x/0 is Inf since x != 0 */
88 MPFR_SET_INF (y);
89 MPFR_SET_SAME_SIGN (y, x);
90 MPFR_SET_DIVBY0 ();
91 MPFR_RET (0);
92 }
93 else /* y = x/1 = x */
94 return mpfr_set (y, x, rnd_mode);
95 }
96 else if (MPFR_UNLIKELY (IS_POW2 (u)))
97 return mpfr_div_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode);
98
99 MPFR_SET_SAME_SIGN (y, x);
100
101 MPFR_TMP_MARK (marker);
102
103 xn = MPFR_LIMB_SIZE (x);
104 yn = MPFR_LIMB_SIZE (y);
105
106 xp = MPFR_MANT (x);
107 yp = MPFR_MANT (y);
108 exp = MPFR_GET_EXP (x);
109
110 dif = yn + 1 - xn;
111
112 /* we need to store yn + 1 = xn + dif limbs of the quotient */
113 tmp = MPFR_TMP_LIMBS_ALLOC (yn + 1);
114
115 /* Notation: {p, n} denotes the integer formed by the n limbs
116 from p[0] to p[n-1]. Let B = 2^GMP_NUMB_BITS.
117 One has: 0 <= {p, n} < B^n. */
118
119 if (dif >= 0)
120 {
121 c = mpn_divrem_1 (tmp, dif, xp, xn, u); /* used all the dividend */
122 /* {xp, xn} = ({tmp, xn+dif} * u + c) * B^(-dif)
123 = ({tmp, yn+1} * u + c) * B^(-dif) */
124 }
125 else /* dif < 0, i.e. xn > yn+1; ignore the (-dif) low limbs from x */
126 {
127 c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, u);
128 /* {xp-dif, yn+1} = {tmp, yn+1} * u + c
129 thus
130 {xp, xn} = {xp, -dif} + {xp-dif, yn+1} * B^(-dif)
131 = {xp, -dif} + ({tmp, yn+1} * u + c) * B^(-dif) */
132 }
133
134 /* Let r = {xp, -dif} / B^(-dif) if dif < 0, r = 0 otherwise; 0 <= r < 1.
135 Then {xp, xn} = ({tmp, yn+1} * u + c + r) * B^(-dif).
136 x / u = ({xp, xn} / u) * B^(-xn) * 2^exp
137 = ({tmp, yn+1} + (c + r) / u) * B^(-(yn+1)) * 2^exp
138 where 0 <= (c + r) / u < 1. */
139
140 for (sb = 0, i = 0; sb == 0 && i < -dif; i++)
141 if (xp[i])
142 sb = 1;
143 /* sb != 0 iff r != 0 */
144
145 /*
146 If the highest limb of the result is 0 (xp[xn-1] < u), remove it.
147 Otherwise, compute the left shift to be performed to normalize.
148 In the latter case, we discard some low bits computed. They
149 contain information useful for the rounding, hence the updating
150 of middle and inexact.
151 */
152
153 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (y));
154 /* sh: number of the trailing bits of y */
155
156 if (tmp[yn] == 0)
157 {
158 MPN_COPY(yp, tmp, yn);
159 exp -= GMP_NUMB_BITS;
160 if (sh == 0) /* round bit is 1 iff (c + r) / u >= 1/2 */
161 {
162 /* In this case tmp[yn]=0 and sh=0, the round bit is not in
163 {tmp,yn+1}. It is 1 iff 2*(c+r) - u >= 0. This means that in
164 some cases, we should look at the most significant bit of r. */
165 if (c >= u - c) /* i.e. 2c >= u: round bit is always 1 */
166 {
167 rb = 1;
168 /* The sticky bit is 1 unless 2c-u = 0 and r = 0. */
169 sb |= 2 * c - u;
170 MPFR_COV_SET (div_ui_sb[0][!!sb]);
171 }
172 else /* 2*c < u */
173 {
174 /* The round bit is 1 iff r >= 1/2 and 2*(c+1/2) = u. */
175 rb = (c == u/2) && (dif < 0) && (xp[-dif-1] & MPFR_LIMB_HIGHBIT);
176 /* If rb is set, we need to recompute sb, since it might have
177 taken into account the msb of xp[-dif-1]. */
178 if (rb)
179 {
180 sb = xp[-dif-1] << 1; /* discard the most significant bit */
181 for (i = 0; sb == 0 && i < -dif-1; i++)
182 if (xp[i])
183 sb = 1;
184 /* The dif < -1 case with sb = 0, i.e. [2][0], will
185 ensure that the body of the loop is covered. */
186 MPFR_COV_SET (div_ui_sb[1 + (dif < -1)][!!sb]);
187 }
188 else
189 {
190 sb |= c;
191 MPFR_COV_SET (div_ui_sb[3][!!sb]);
192 }
193 }
194 }
195 else
196 {
197 /* round bit is in tmp[0] */
198 rb = tmp[0] & (MPFR_LIMB_ONE << (sh - 1));
199 sb |= (tmp[0] & MPFR_LIMB_MASK(sh - 1)) | c;
200 MPFR_COV_SET (div_ui_sb[4+!!rb][!!sb]);
201 }
202 }
203 else /* tmp[yn] != 0 */
204 {
205 int shlz;
206 mp_limb_t w;
207
208 MPFR_ASSERTD (tmp[yn] != 0);
209 count_leading_zeros (shlz, tmp[yn]);
210
211 MPFR_ASSERTD (u >= 2); /* see special cases at the beginning */
212 MPFR_ASSERTD (shlz > 0); /* since u >= 2 */
213
214 /* shift left to normalize */
215 w = tmp[0] << shlz;
216 mpn_lshift (yp, tmp + 1, yn, shlz);
217 yp[0] |= tmp[0] >> (GMP_NUMB_BITS - shlz);
218 /* now {yp, yn} is the approximate quotient, w is the next limb */
219
220 if (sh == 0) /* round bit is upper bit from w */
221 {
222 rb = w & MPFR_LIMB_HIGHBIT;
223 sb |= (w - rb) | c;
224 MPFR_COV_SET (div_ui_sb[6+!!rb][!!sb]);
225 }
226 else
227 {
228 rb = yp[0] & (MPFR_LIMB_ONE << (sh - 1));
229 sb |= (yp[0] & MPFR_LIMB_MASK(sh - 1)) | w | c;
230 MPFR_COV_SET (div_ui_sb[8+!!rb][!!sb]);
231 }
232
233 exp -= shlz;
234 }
235
236 d = yp[0] & MPFR_LIMB_MASK (sh);
237 yp[0] ^= d; /* clear the lowest sh bits */
238
239 MPFR_TMP_FREE (marker);
240
241 if (MPFR_UNLIKELY (exp < __gmpfr_emin - 1))
242 return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
243 MPFR_SIGN (y));
244
245 if (MPFR_UNLIKELY (rb == 0 && sb == 0))
246 inexact = 0; /* result is exact */
247 else
248 {
249 int nexttoinf;
250
251 MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (y));
252 switch (rnd_mode)
253 {
254 case MPFR_RNDZ:
255 case MPFR_RNDF:
256 inexact = - MPFR_INT_SIGN (y); /* result is inexact */
257 nexttoinf = 0;
258 break;
259
260 case MPFR_RNDA:
261 inexact = MPFR_INT_SIGN (y);
262 nexttoinf = 1;
263 break;
264
265 default: /* should be MPFR_RNDN */
266 MPFR_ASSERTD (rnd_mode == MPFR_RNDN);
267 /* We have one more significant bit in yn. */
268 if (rb == 0)
269 {
270 inexact = - MPFR_INT_SIGN (y);
271 nexttoinf = 0;
272 }
273 else if (sb != 0) /* necessarily rb != 0 */
274 {
275 inexact = MPFR_INT_SIGN (y);
276 nexttoinf = 1;
277 }
278 else /* middle case */
279 {
280 if (yp[0] & (MPFR_LIMB_ONE << sh))
281 {
282 inexact = MPFR_INT_SIGN (y);
283 nexttoinf = 1;
284 }
285 else
286 {
287 inexact = - MPFR_INT_SIGN (y);
288 nexttoinf = 0;
289 }
290 }
291 }
292 if (nexttoinf &&
293 MPFR_UNLIKELY (mpn_add_1 (yp, yp, yn, MPFR_LIMB_ONE << sh)))
294 {
295 exp++;
296 yp[yn-1] = MPFR_LIMB_HIGHBIT;
297 }
298 }
299
300 /* Set the exponent. Warning! One may still have an underflow. */
301 MPFR_EXP (y) = exp;
302 #else /* MPFR_LONG_WITHIN_LIMB */
303 mpfr_t uu;
304 MPFR_SAVE_EXPO_DECL (expo);
305
306 MPFR_SAVE_EXPO_MARK (expo);
307 mpfr_init2 (uu, sizeof (unsigned long) * CHAR_BIT);
308 mpfr_set_ui (uu, u, MPFR_RNDZ);
309 inexact = mpfr_div (y, x, uu, rnd_mode);
310 mpfr_clear (uu);
311 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
312 MPFR_SAVE_EXPO_FREE (expo);
313 #endif
314
315 return mpfr_check_range (y, inexact, rnd_mode);
316 }