1/* mpfr_div_{ui,si} -- divide a floating-point number by a machine integer
2
3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel 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
20http://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/* returns 0 if result exact, non-zero otherwise */
27int
28mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode)
29{
30  long i;
31  int sh;
32  mp_size_t xn, yn, dif;
33  mp_limb_t *xp, *yp, *tmp, c, d;
34  mpfr_exp_t exp;
35  int inexact, middle = 1, nexttoinf;
36  MPFR_TMP_DECL(marker);
37
38  MPFR_LOG_FUNC
39    (("x[%Pu]=%.*Rg u=%lu rnd=%d",
40      mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
41     ("y[%Pu]=%.*Rg inexact=%d",
42      mpfr_get_prec(y), mpfr_log_prec, y, inexact));
43
44  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
45    {
46      if (MPFR_IS_NAN (x))
47        {
48          MPFR_SET_NAN (y);
49          MPFR_RET_NAN;
50        }
51      else if (MPFR_IS_INF (x))
52        {
53          MPFR_SET_INF (y);
54          MPFR_SET_SAME_SIGN (y, x);
55          MPFR_RET (0);
56        }
57      else
58        {
59          MPFR_ASSERTD (MPFR_IS_ZERO(x));
60          if (u == 0) /* 0/0 is NaN */
61            {
62              MPFR_SET_NAN(y);
63              MPFR_RET_NAN;
64            }
65          else
66            {
67              MPFR_SET_ZERO(y);
68              MPFR_SET_SAME_SIGN (y, x);
69              MPFR_RET(0);
70            }
71        }
72    }
73  else if (MPFR_UNLIKELY (u <= 1))
74    {
75      if (u < 1)
76        {
77          /* x/0 is Inf since x != 0*/
78          MPFR_SET_INF (y);
79          MPFR_SET_SAME_SIGN (y, x);
80          mpfr_set_divby0 ();
81          MPFR_RET (0);
82        }
83      else /* y = x/1 = x */
84        return mpfr_set (y, x, rnd_mode);
85    }
86  else if (MPFR_UNLIKELY (IS_POW2 (u)))
87    return mpfr_div_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode);
88
89  MPFR_SET_SAME_SIGN (y, x);
90
91  MPFR_TMP_MARK (marker);
92  xn = MPFR_LIMB_SIZE (x);
93  yn = MPFR_LIMB_SIZE (y);
94
95  xp = MPFR_MANT (x);
96  yp = MPFR_MANT (y);
97  exp = MPFR_GET_EXP (x);
98
99  dif = yn + 1 - xn;
100
101  /* we need to store yn+1 = xn + dif limbs of the quotient */
102  /* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */
103  tmp = MPFR_TMP_LIMBS_ALLOC (yn + 1);
104
105  c = (mp_limb_t) u;
106  MPFR_ASSERTN (u == c);
107  if (dif >= 0)
108    c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */
109  else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */
110    c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c);
111
112  inexact = (c != 0);
113
114  /* First pass in estimating next bit of the quotient, in case of RNDN    *
115   * In case we just have the right number of bits (postpone this ?),      *
116   * we need to check whether the remainder is more or less than half      *
117   * the divisor. The test must be performed with a subtraction, so as     *
118   * to prevent carries.                                                   */
119
120  if (MPFR_LIKELY (rnd_mode == MPFR_RNDN))
121    {
122      if (c < (mp_limb_t) u - c) /* We have u > c */
123        middle = -1;
124      else if (c > (mp_limb_t) u - c)
125        middle = 1;
126      else
127        middle = 0; /* exactly in the middle */
128    }
129
130  /* If we believe that we are right in the middle or exact, we should check
131     that we did not neglect any word of x (division large / 1 -> small). */
132
133  for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++)
134    if (xp[i])
135      inexact = middle = 1; /* larger than middle */
136
137  /*
138     If the high limb of the result is 0 (xp[xn-1] < u), remove it.
139     Otherwise, compute the left shift to be performed to normalize.
140     In the latter case, we discard some low bits computed. They
141     contain information useful for the rounding, hence the updating
142     of middle and inexact.
143  */
144
145  if (tmp[yn] == 0)
146    {
147      MPN_COPY(yp, tmp, yn);
148      exp -= GMP_NUMB_BITS;
149    }
150  else
151    {
152      int shlz;
153
154      count_leading_zeros (shlz, tmp[yn]);
155
156      /* shift left to normalize */
157      if (MPFR_LIKELY (shlz != 0))
158        {
159          mp_limb_t w = tmp[0] << shlz;
160
161          mpn_lshift (yp, tmp + 1, yn, shlz);
162          yp[0] += tmp[0] >> (GMP_NUMB_BITS - shlz);
163
164          if (w > (MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1)))
165            { middle = 1; }
166          else if (w < (MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1)))
167            { middle = -1; }
168          else
169            { middle = (c != 0); }
170
171          inexact = inexact || (w != 0);
172          exp -= shlz;
173        }
174      else
175        { /* this happens only if u == 1 and xp[xn-1] >=
176             1<<(GMP_NUMB_BITS-1). It might be better to handle the
177             u == 1 case seperately ?
178          */
179
180             MPN_COPY (yp, tmp + 1, yn);
181        }
182    }
183
184  MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (y));
185  /* it remains sh bits in less significant limb of y */
186
187  d = *yp & MPFR_LIMB_MASK (sh);
188  *yp ^= d; /* set to zero lowest sh bits */
189
190  MPFR_TMP_FREE (marker);
191
192  if (exp < __gmpfr_emin - 1)
193    return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
194                           MPFR_SIGN (y));
195
196  if (MPFR_UNLIKELY (d == 0 && inexact == 0))
197    nexttoinf = 0;  /* result is exact */
198  else
199    {
200      MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN (y));
201      switch (rnd_mode)
202        {
203        case MPFR_RNDZ:
204          inexact = - MPFR_INT_SIGN (y);  /* result is inexact */
205          nexttoinf = 0;
206          break;
207
208        case MPFR_RNDA:
209          inexact = MPFR_INT_SIGN (y);
210          nexttoinf = 1;
211          break;
212
213        default: /* should be MPFR_RNDN */
214          MPFR_ASSERTD (rnd_mode == MPFR_RNDN);
215          /* We have one more significant bit in yn. */
216          if (sh && d < (MPFR_LIMB_ONE << (sh - 1)))
217            {
218              inexact = - MPFR_INT_SIGN (y);
219              nexttoinf = 0;
220            }
221          else if (sh && d > (MPFR_LIMB_ONE << (sh - 1)))
222            {
223              inexact = MPFR_INT_SIGN (y);
224              nexttoinf = 1;
225            }
226          else /* sh = 0 or d = 1 << (sh-1) */
227            {
228              /* The first case is "false" even rounding (significant bits
229                 indicate even rounding, but the result is inexact, so up) ;
230                 The second case is the case where middle should be used to
231                 decide the direction of rounding (no further bit computed) ;
232                 The third is the true even rounding.
233              */
234              if ((sh && inexact) || (!sh && middle > 0) ||
235                  (!inexact && *yp & (MPFR_LIMB_ONE << sh)))
236                {
237                  inexact = MPFR_INT_SIGN (y);
238                  nexttoinf = 1;
239                }
240              else
241                {
242                  inexact = - MPFR_INT_SIGN (y);
243                  nexttoinf = 0;
244                }
245            }
246        }
247    }
248
249  if (nexttoinf &&
250      MPFR_UNLIKELY (mpn_add_1 (yp, yp, yn, MPFR_LIMB_ONE << sh)))
251    {
252      exp++;
253      yp[yn-1] = MPFR_LIMB_HIGHBIT;
254    }
255
256  /* Set the exponent. Warning! One may still have an underflow. */
257  MPFR_EXP (y) = exp;
258
259  return mpfr_check_range (y, inexact, rnd_mode);
260}
261
262int
263mpfr_div_si (mpfr_ptr y, mpfr_srcptr x, long int u, mpfr_rnd_t rnd_mode)
264{
265  int res;
266
267  MPFR_LOG_FUNC
268    (("x[%Pu]=%.*Rg u=%ld rnd=%d",
269      mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
270     ("y[%Pu]=%.*Rg inexact=%d",
271      mpfr_get_prec(y), mpfr_log_prec, y, res));
272
273  if (u >= 0)
274    res = mpfr_div_ui (y, x, u, rnd_mode);
275  else
276    {
277      res = -mpfr_div_ui (y, x, -u, MPFR_INVERT_RND (rnd_mode));
278      MPFR_CHANGE_SIGN (y);
279    }
280  return res;
281}
282