1/* mpfr_div_ui -- divide a floating-point number by a machine integer
2
3Copyright 1999-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#define MPFR_NEED_LONGLONG_H
24#include "mpfr-impl.h"
25
26#ifdef MPFR_COV_CHECK
27int __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
32MPFR_HOT_FUNCTION_ATTR int
33mpfr_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[%Pu]=%.*Rg u=%lu rnd=%d",
50      mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
51     ("y[%Pu]=%.*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}
317