1 /* Implementations of operations between mpfr and mpz/mpq data
2
3 Copyright 2001, 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 /* Init and set a mpfr_t with enough precision to store a mpz.
27 This function should be called in the extended exponent range. */
28 static void
init_set_z(mpfr_ptr t,mpz_srcptr z)29 init_set_z (mpfr_ptr t, mpz_srcptr z)
30 {
31 mpfr_prec_t p;
32 int i;
33
34 if (mpz_size (z) <= 1)
35 p = GMP_NUMB_BITS;
36 else
37 MPFR_MPZ_SIZEINBASE2 (p, z);
38 mpfr_init2 (t, p);
39 i = mpfr_set_z (t, z, MPFR_RNDN);
40 /* Possible assertion failure in case of overflow. Such cases,
41 which imply that z is huge (if the function is called in
42 the extended exponent range), are currently not supported,
43 just like precisions around MPFR_PREC_MAX. */
44 MPFR_ASSERTN (i == 0); (void) i; /* use i to avoid a warning */
45 }
46
47 /* Init, set a mpfr_t with enough precision to store a mpz_t without round,
48 call the function, and clear the allocated mpfr_t */
49 static int
foo(mpfr_ptr x,mpfr_srcptr y,mpz_srcptr z,mpfr_rnd_t r,int (* f)(mpfr_ptr,mpfr_srcptr,mpfr_srcptr,mpfr_rnd_t))50 foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mpfr_rnd_t r,
51 int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t))
52 {
53 mpfr_t t;
54 int i;
55 MPFR_SAVE_EXPO_DECL (expo);
56
57 MPFR_SAVE_EXPO_MARK (expo);
58 init_set_z (t, z); /* There should be no exceptions. */
59 i = (*f) (x, y, t, r);
60 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
61 mpfr_clear (t);
62 MPFR_SAVE_EXPO_FREE (expo);
63 return mpfr_check_range (x, i, r);
64 }
65
66 static int
foo2(mpfr_ptr x,mpz_srcptr y,mpfr_srcptr z,mpfr_rnd_t r,int (* f)(mpfr_ptr,mpfr_srcptr,mpfr_srcptr,mpfr_rnd_t))67 foo2 (mpfr_ptr x, mpz_srcptr y, mpfr_srcptr z, mpfr_rnd_t r,
68 int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t))
69 {
70 mpfr_t t;
71 int i;
72 MPFR_SAVE_EXPO_DECL (expo);
73
74 MPFR_SAVE_EXPO_MARK (expo);
75 init_set_z (t, y); /* There should be no exceptions. */
76 i = (*f) (x, t, z, r);
77 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
78 mpfr_clear (t);
79 MPFR_SAVE_EXPO_FREE (expo);
80 return mpfr_check_range (x, i, r);
81 }
82
83 int
mpfr_mul_z(mpfr_ptr y,mpfr_srcptr x,mpz_srcptr z,mpfr_rnd_t r)84 mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
85 {
86 return foo (y, x, z, r, mpfr_mul);
87 }
88
89 int
mpfr_div_z(mpfr_ptr y,mpfr_srcptr x,mpz_srcptr z,mpfr_rnd_t r)90 mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
91 {
92 return foo (y, x, z, r, mpfr_div);
93 }
94
95 int
mpfr_add_z(mpfr_ptr y,mpfr_srcptr x,mpz_srcptr z,mpfr_rnd_t r)96 mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
97 {
98 /* Mpz 0 is unsigned */
99 if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
100 return mpfr_set (y, x, r);
101 else
102 return foo (y, x, z, r, mpfr_add);
103 }
104
105 int
mpfr_sub_z(mpfr_ptr y,mpfr_srcptr x,mpz_srcptr z,mpfr_rnd_t r)106 mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
107 {
108 /* Mpz 0 is unsigned */
109 if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
110 return mpfr_set (y, x, r);
111 else
112 return foo (y, x, z, r, mpfr_sub);
113 }
114
115 int
mpfr_z_sub(mpfr_ptr y,mpz_srcptr x,mpfr_srcptr z,mpfr_rnd_t r)116 mpfr_z_sub (mpfr_ptr y, mpz_srcptr x, mpfr_srcptr z, mpfr_rnd_t r)
117 {
118 /* Mpz 0 is unsigned */
119 if (MPFR_UNLIKELY (mpz_sgn (x) == 0))
120 return mpfr_neg (y, z, r);
121 else
122 return foo2 (y, x, z, r, mpfr_sub);
123 }
124
125 int
mpfr_cmp_z(mpfr_srcptr x,mpz_srcptr z)126 mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z)
127 {
128 mpfr_t t;
129 int res;
130 mpfr_prec_t p;
131 unsigned int flags;
132
133 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
134 return mpfr_cmp_si (x, mpz_sgn (z));
135
136 if (mpz_size (z) <= 1)
137 p = GMP_NUMB_BITS;
138 else
139 MPFR_MPZ_SIZEINBASE2 (p, z);
140 mpfr_init2 (t, p);
141 flags = __gmpfr_flags;
142 if (mpfr_set_z (t, z, MPFR_RNDN))
143 {
144 /* overflow (t is an infinity) or underflow */
145 mpfr_div_2ui (t, t, 2, MPFR_RNDZ); /* if underflow, set t to zero */
146 __gmpfr_flags = flags; /* restore the flags */
147 /* The real value of t (= z), which falls outside the exponent range,
148 has been replaced by an equivalent value for the comparison: zero
149 or an infinity. */
150 }
151 res = mpfr_cmp (x, t);
152 mpfr_clear (t);
153 return res;
154 }
155
156 /* Compute y = RND(x*n/d), where n and d are mpz integers.
157 An integer 0 is assumed to have a positive sign.
158 This function is used by mpfr_mul_q and mpfr_div_q.
159 Note: the status of the rational 0/(-1) is not clear (if there is
160 a signed infinity, there should be a signed zero). But infinities
161 are not currently supported/documented in GMP, and if the rational
162 is canonicalized as it should be, the case 0/(-1) cannot occur. */
163 static int
mpfr_muldiv_z(mpfr_ptr y,mpfr_srcptr x,mpz_srcptr n,mpz_srcptr d,mpfr_rnd_t rnd_mode)164 mpfr_muldiv_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr n, mpz_srcptr d,
165 mpfr_rnd_t rnd_mode)
166 {
167 if (MPFR_UNLIKELY (mpz_sgn (n) == 0))
168 {
169 if (MPFR_UNLIKELY (mpz_sgn (d) == 0))
170 MPFR_SET_NAN (y);
171 else
172 {
173 mpfr_mul_ui (y, x, 0, MPFR_RNDN); /* exact: +0, -0 or NaN */
174 if (MPFR_UNLIKELY (mpz_sgn (d) < 0))
175 MPFR_CHANGE_SIGN (y);
176 }
177 return 0;
178 }
179 else if (MPFR_UNLIKELY (mpz_sgn (d) == 0))
180 {
181 mpfr_div_ui (y, x, 0, MPFR_RNDN); /* exact: +Inf, -Inf or NaN */
182 if (MPFR_UNLIKELY (mpz_sgn (n) < 0))
183 MPFR_CHANGE_SIGN (y);
184 return 0;
185 }
186 else
187 {
188 mpfr_prec_t p;
189 mpfr_t tmp;
190 int inexact;
191 MPFR_SAVE_EXPO_DECL (expo);
192
193 MPFR_SAVE_EXPO_MARK (expo);
194
195 /* With the current MPFR code, using mpfr_mul_z and mpfr_div_z
196 for the general case should be faster than doing everything
197 in mpn, mpz and/or mpq. MPFR_SAVE_EXPO_MARK could be avoided
198 here, but it would be more difficult to handle corner cases. */
199 MPFR_MPZ_SIZEINBASE2 (p, n);
200 mpfr_init2 (tmp, MPFR_PREC (x) + p);
201 inexact = mpfr_mul_z (tmp, x, n, MPFR_RNDN);
202 /* Since |n| >= 1, an underflow is not possible. And the precision of
203 tmp has been chosen so that inexact != 0 iff there's an overflow. */
204 if (MPFR_UNLIKELY (inexact != 0))
205 {
206 mpfr_t x0;
207 mpfr_exp_t ex;
208 MPFR_BLOCK_DECL (flags);
209
210 /* intermediate overflow case */
211 MPFR_ASSERTD (mpfr_inf_p (tmp));
212 ex = MPFR_GET_EXP (x); /* x is a pure FP number */
213 MPFR_ALIAS (x0, x, MPFR_SIGN(x), 0); /* x0 = x / 2^ex */
214 MPFR_BLOCK (flags,
215 inexact = mpfr_mul_z (tmp, x0, n, MPFR_RNDN);
216 MPFR_ASSERTD (inexact == 0);
217 inexact = mpfr_div_z (y, tmp, d, rnd_mode);
218 /* Just in case the division underflows
219 (highly unlikely, not supported)... */
220 MPFR_ASSERTN (!MPFR_BLOCK_EXCEP));
221 MPFR_EXP (y) += ex;
222 /* Detect highly unlikely, not supported corner cases... */
223 MPFR_ASSERTN (MPFR_EXP (y) >= __gmpfr_emin && MPFR_IS_PURE_FP (y));
224 /* The potential overflow will be detected by mpfr_check_range. */
225 }
226 else
227 inexact = mpfr_div_z (y, tmp, d, rnd_mode);
228
229 mpfr_clear (tmp);
230
231 MPFR_SAVE_EXPO_FREE (expo);
232 return mpfr_check_range (y, inexact, rnd_mode);
233 }
234 }
235
236 int
mpfr_mul_q(mpfr_ptr y,mpfr_srcptr x,mpq_srcptr z,mpfr_rnd_t rnd_mode)237 mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode)
238 {
239 return mpfr_muldiv_z (y, x, mpq_numref (z), mpq_denref (z), rnd_mode);
240 }
241
242 int
mpfr_div_q(mpfr_ptr y,mpfr_srcptr x,mpq_srcptr z,mpfr_rnd_t rnd_mode)243 mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode)
244 {
245 return mpfr_muldiv_z (y, x, mpq_denref (z), mpq_numref (z), rnd_mode);
246 }
247
248 int
mpfr_add_q(mpfr_ptr y,mpfr_srcptr x,mpq_srcptr z,mpfr_rnd_t rnd_mode)249 mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode)
250 {
251 mpfr_t t,q;
252 mpfr_prec_t p;
253 mpfr_exp_t err;
254 int res;
255 MPFR_SAVE_EXPO_DECL (expo);
256 MPFR_ZIV_DECL (loop);
257
258 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
259 {
260 if (MPFR_IS_NAN (x))
261 {
262 MPFR_SET_NAN (y);
263 MPFR_RET_NAN;
264 }
265 else if (MPFR_IS_INF (x))
266 {
267 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0 &&
268 MPFR_MULT_SIGN (mpz_sgn (mpq_numref (z)),
269 MPFR_SIGN (x)) <= 0))
270 {
271 MPFR_SET_NAN (y);
272 MPFR_RET_NAN;
273 }
274 MPFR_SET_INF (y);
275 MPFR_SET_SAME_SIGN (y, x);
276 MPFR_RET (0);
277 }
278 else
279 {
280 MPFR_ASSERTD (MPFR_IS_ZERO (x));
281 if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
282 return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
283 else
284 return mpfr_set_q (y, z, rnd_mode);
285 }
286 }
287
288 MPFR_SAVE_EXPO_MARK (expo);
289
290 p = MPFR_PREC (y) + 10;
291 mpfr_init2 (t, p);
292 mpfr_init2 (q, p);
293
294 MPFR_ZIV_INIT (loop, p);
295 for (;;)
296 {
297 MPFR_BLOCK_DECL (flags);
298
299 res = mpfr_set_q (q, z, MPFR_RNDN); /* Error <= 1/2 ulp(q) */
300 /* If z if @INF@ (1/0), res = 0, so it quits immediately */
301 if (MPFR_UNLIKELY (res == 0))
302 /* Result is exact so we can add it directly! */
303 {
304 res = mpfr_add (y, x, q, rnd_mode);
305 break;
306 }
307 MPFR_BLOCK (flags, mpfr_add (t, x, q, MPFR_RNDN));
308 /* Error on t is <= 1/2 ulp(t), except in case of overflow/underflow,
309 but such an exception is very unlikely as it would be possible
310 only if q has a huge numerator or denominator. Not supported! */
311 MPFR_ASSERTN (! (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)));
312 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
313 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
314 <= 2^(EXP(q)-EXP(t))
315 If EXP(q)-EXP(t)<0, <= 2^0 */
316 /* We can get 0, but we can't round since q is inexact */
317 if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
318 {
319 err = (mpfr_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
320 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode)))
321 {
322 res = mpfr_set (y, t, rnd_mode);
323 break;
324 }
325 }
326 MPFR_ZIV_NEXT (loop, p);
327 mpfr_set_prec (t, p);
328 mpfr_set_prec (q, p);
329 }
330 MPFR_ZIV_FREE (loop);
331 mpfr_clear (t);
332 mpfr_clear (q);
333
334 MPFR_SAVE_EXPO_FREE (expo);
335 return mpfr_check_range (y, res, rnd_mode);
336 }
337
338 int
mpfr_sub_q(mpfr_ptr y,mpfr_srcptr x,mpq_srcptr z,mpfr_rnd_t rnd_mode)339 mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mpfr_rnd_t rnd_mode)
340 {
341 mpfr_t t,q;
342 mpfr_prec_t p;
343 int res;
344 mpfr_exp_t err;
345 MPFR_SAVE_EXPO_DECL (expo);
346 MPFR_ZIV_DECL (loop);
347
348 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
349 {
350 if (MPFR_IS_NAN (x))
351 {
352 MPFR_SET_NAN (y);
353 MPFR_RET_NAN;
354 }
355 else if (MPFR_IS_INF (x))
356 {
357 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0 &&
358 MPFR_MULT_SIGN (mpz_sgn (mpq_numref (z)),
359 MPFR_SIGN (x)) >= 0))
360 {
361 MPFR_SET_NAN (y);
362 MPFR_RET_NAN;
363 }
364 MPFR_SET_INF (y);
365 MPFR_SET_SAME_SIGN (y, x);
366 MPFR_RET (0);
367 }
368 else
369 {
370 MPFR_ASSERTD (MPFR_IS_ZERO (x));
371
372 if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
373 return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
374 else
375 {
376 res = mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode));
377 MPFR_CHANGE_SIGN (y);
378 return -res;
379 }
380 }
381 }
382
383 MPFR_SAVE_EXPO_MARK (expo);
384
385 p = MPFR_PREC (y) + 10;
386 mpfr_init2 (t, p);
387 mpfr_init2 (q, p);
388
389 MPFR_ZIV_INIT (loop, p);
390 for(;;)
391 {
392 MPFR_BLOCK_DECL (flags);
393
394 res = mpfr_set_q(q, z, MPFR_RNDN); /* Error <= 1/2 ulp(q) */
395 /* If z if @INF@ (1/0), res = 0, so it quits immediately */
396 if (MPFR_UNLIKELY (res == 0))
397 /* Result is exact so we can add it directly!*/
398 {
399 res = mpfr_sub (y, x, q, rnd_mode);
400 break;
401 }
402 MPFR_BLOCK (flags, mpfr_sub (t, x, q, MPFR_RNDN));
403 /* Error on t is <= 1/2 ulp(t), except in case of overflow/underflow,
404 but such an exception is very unlikely as it would be possible
405 only if q has a huge numerator or denominator. Not supported! */
406 MPFR_ASSERTN (! (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)));
407 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
408 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
409 <= 2^(EXP(q)-EXP(t))
410 If EXP(q)-EXP(t)<0, <= 2^0 */
411 /* We can get 0, but we can't round since q is inexact */
412 if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
413 {
414 err = (mpfr_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
415 res = MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode);
416 if (MPFR_LIKELY (res != 0)) /* We can round! */
417 {
418 res = mpfr_set (y, t, rnd_mode);
419 break;
420 }
421 }
422 MPFR_ZIV_NEXT (loop, p);
423 mpfr_set_prec (t, p);
424 mpfr_set_prec (q, p);
425 }
426 MPFR_ZIV_FREE (loop);
427 mpfr_clear (t);
428 mpfr_clear (q);
429
430 MPFR_SAVE_EXPO_FREE (expo);
431 return mpfr_check_range (y, res, rnd_mode);
432 }
433
434 int
mpfr_cmp_q(mpfr_srcptr x,mpq_srcptr q)435 mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr q)
436 {
437 mpfr_t t;
438 int res;
439 mpfr_prec_t p;
440 MPFR_SAVE_EXPO_DECL (expo);
441
442 if (MPFR_UNLIKELY (mpq_denref (q) == 0))
443 {
444 /* q is an infinity or NaN */
445 mpfr_init2 (t, 2);
446 mpfr_set_q (t, q, MPFR_RNDN);
447 res = mpfr_cmp (x, t);
448 mpfr_clear (t);
449 return res;
450 }
451
452 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
453 return mpfr_cmp_si (x, mpq_sgn (q));
454
455 MPFR_SAVE_EXPO_MARK (expo);
456
457 /* x < a/b ? <=> x*b < a */
458 MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (q));
459 mpfr_init2 (t, MPFR_PREC(x) + p);
460 res = mpfr_mul_z (t, x, mpq_denref (q), MPFR_RNDN);
461 MPFR_ASSERTD (res == 0);
462 res = mpfr_cmp_z (t, mpq_numref (q));
463 mpfr_clear (t);
464
465 MPFR_SAVE_EXPO_FREE (expo);
466 return res;
467 }
468
469 int
mpfr_cmp_f(mpfr_srcptr x,mpf_srcptr z)470 mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z)
471 {
472 mpfr_t t;
473 int res;
474 MPFR_SAVE_EXPO_DECL (expo);
475
476 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
477 return mpfr_cmp_si (x, mpf_sgn (z));
478
479 MPFR_SAVE_EXPO_MARK (expo);
480
481 mpfr_init2 (t, MPFR_PREC_MIN + ABS(SIZ(z)) * GMP_NUMB_BITS );
482 res = mpfr_set_f (t, z, MPFR_RNDN);
483 MPFR_ASSERTD (res == 0);
484 res = mpfr_cmp (x, t);
485 mpfr_clear (t);
486
487 MPFR_SAVE_EXPO_FREE (expo);
488 return res;
489 }
490