1 /* mpfr_div_{ui,si} -- divide a floating-point number by a machine integer
2
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel 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 http://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 /* returns 0 if result exact, non-zero otherwise */
27 int
mpfr_div_ui(mpfr_ptr y,mpfr_srcptr x,unsigned long int u,mpfr_rnd_t rnd_mode)28 mpfr_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
262 int
mpfr_div_si(mpfr_ptr y,mpfr_srcptr x,long int u,mpfr_rnd_t rnd_mode)263 mpfr_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