1 /* Implementations of operations between mpfr and mpz/mpq data
3 Copyright 2001, 2003-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
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. */
29 init_set_z (mpfr_ptr t
, mpz_srcptr z
)
34 if (mpz_size (z
) <= 1)
37 MPFR_MPZ_SIZEINBASE2 (p
, z
);
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 */
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 */
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
))
55 MPFR_SAVE_EXPO_DECL (expo
);
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
);
62 MPFR_SAVE_EXPO_FREE (expo
);
63 return mpfr_check_range (x
, i
, r
);
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
))
72 MPFR_SAVE_EXPO_DECL (expo
);
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
);
79 MPFR_SAVE_EXPO_FREE (expo
);
80 return mpfr_check_range (x
, i
, r
);
84 mpfr_mul_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mpfr_rnd_t r
)
86 return foo (y
, x
, z
, r
, mpfr_mul
);
90 mpfr_div_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mpfr_rnd_t r
)
92 return foo (y
, x
, z
, r
, mpfr_div
);
96 mpfr_add_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mpfr_rnd_t r
)
98 /* Mpz 0 is unsigned */
99 if (MPFR_UNLIKELY (mpz_sgn (z
) == 0))
100 return mpfr_set (y
, x
, r
);
102 return foo (y
, x
, z
, r
, mpfr_add
);
106 mpfr_sub_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mpfr_rnd_t r
)
108 /* Mpz 0 is unsigned */
109 if (MPFR_UNLIKELY (mpz_sgn (z
) == 0))
110 return mpfr_set (y
, x
, r
);
112 return foo (y
, x
, z
, r
, mpfr_sub
);
116 mpfr_z_sub (mpfr_ptr y
, mpz_srcptr x
, mpfr_srcptr z
, mpfr_rnd_t r
)
118 /* Mpz 0 is unsigned */
119 if (MPFR_UNLIKELY (mpz_sgn (x
) == 0))
120 return mpfr_neg (y
, z
, r
);
122 return foo2 (y
, x
, z
, r
, mpfr_sub
);
126 mpfr_cmp_z (mpfr_srcptr x
, mpz_srcptr z
)
133 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
134 return mpfr_cmp_si (x
, mpz_sgn (z
));
136 if (mpz_size (z
) <= 1)
139 MPFR_MPZ_SIZEINBASE2 (p
, z
);
141 flags
= __gmpfr_flags
;
142 if (mpfr_set_z (t
, z
, MPFR_RNDN
))
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
151 res
= mpfr_cmp (x
, t
);
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. */
164 mpfr_muldiv_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr n
, mpz_srcptr d
,
167 if (MPFR_UNLIKELY (mpz_sgn (n
) == 0))
169 if (MPFR_UNLIKELY (mpz_sgn (d
) == 0))
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
);
179 else if (MPFR_UNLIKELY (mpz_sgn (d
) == 0))
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
);
191 MPFR_SAVE_EXPO_DECL (expo
);
193 MPFR_SAVE_EXPO_MARK (expo
);
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))
208 MPFR_BLOCK_DECL (flags
);
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 */
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
));
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. */
227 inexact
= mpfr_div_z (y
, tmp
, d
, rnd_mode
);
231 MPFR_SAVE_EXPO_FREE (expo
);
232 return mpfr_check_range (y
, inexact
, rnd_mode
);
237 mpfr_mul_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
, mpfr_rnd_t rnd_mode
)
239 return mpfr_muldiv_z (y
, x
, mpq_numref (z
), mpq_denref (z
), rnd_mode
);
243 mpfr_div_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
, mpfr_rnd_t rnd_mode
)
245 return mpfr_muldiv_z (y
, x
, mpq_denref (z
), mpq_numref (z
), rnd_mode
);
249 mpfr_add_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
, mpfr_rnd_t rnd_mode
)
255 MPFR_SAVE_EXPO_DECL (expo
);
256 MPFR_ZIV_DECL (loop
);
258 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
265 else if (MPFR_IS_INF (x
))
267 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z
)) == 0 &&
268 MPFR_MULT_SIGN (mpz_sgn (mpq_numref (z
)),
269 MPFR_SIGN (x
)) <= 0))
275 MPFR_SET_SAME_SIGN (y
, x
);
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 */
284 return mpfr_set_q (y
, z
, rnd_mode
);
288 MPFR_SAVE_EXPO_MARK (expo
);
290 p
= MPFR_PREC (y
) + 10;
294 MPFR_ZIV_INIT (loop
, p
);
297 MPFR_BLOCK_DECL (flags
);
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! */
304 res
= mpfr_add (y
, x
, q
, rnd_mode
);
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)))
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
)))
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
)))
322 res
= mpfr_set (y
, t
, rnd_mode
);
326 MPFR_ZIV_NEXT (loop
, p
);
327 mpfr_set_prec (t
, p
);
328 mpfr_set_prec (q
, p
);
330 MPFR_ZIV_FREE (loop
);
334 MPFR_SAVE_EXPO_FREE (expo
);
335 return mpfr_check_range (y
, res
, rnd_mode
);
339 mpfr_sub_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
,mpfr_rnd_t rnd_mode
)
345 MPFR_SAVE_EXPO_DECL (expo
);
346 MPFR_ZIV_DECL (loop
);
348 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
355 else if (MPFR_IS_INF (x
))
357 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z
)) == 0 &&
358 MPFR_MULT_SIGN (mpz_sgn (mpq_numref (z
)),
359 MPFR_SIGN (x
)) >= 0))
365 MPFR_SET_SAME_SIGN (y
, x
);
370 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
372 if (MPFR_UNLIKELY (mpq_sgn (z
) == 0))
373 return mpfr_set (y
, x
, rnd_mode
); /* signed 0 - Unsigned 0 */
376 res
= mpfr_set_q (y
, z
, MPFR_INVERT_RND (rnd_mode
));
377 MPFR_CHANGE_SIGN (y
);
383 MPFR_SAVE_EXPO_MARK (expo
);
385 p
= MPFR_PREC (y
) + 10;
389 MPFR_ZIV_INIT (loop
, p
);
392 MPFR_BLOCK_DECL (flags
);
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!*/
399 res
= mpfr_sub (y
, x
, q
, rnd_mode
);
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)))
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
)))
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! */
418 res
= mpfr_set (y
, t
, rnd_mode
);
422 MPFR_ZIV_NEXT (loop
, p
);
423 mpfr_set_prec (t
, p
);
424 mpfr_set_prec (q
, p
);
426 MPFR_ZIV_FREE (loop
);
430 MPFR_SAVE_EXPO_FREE (expo
);
431 return mpfr_check_range (y
, res
, rnd_mode
);
435 mpfr_cmp_q (mpfr_srcptr x
, mpq_srcptr q
)
440 MPFR_SAVE_EXPO_DECL (expo
);
442 if (MPFR_UNLIKELY (mpq_denref (q
) == 0))
444 /* q is an infinity or NaN */
446 mpfr_set_q (t
, q
, MPFR_RNDN
);
447 res
= mpfr_cmp (x
, t
);
452 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
453 return mpfr_cmp_si (x
, mpq_sgn (q
));
455 MPFR_SAVE_EXPO_MARK (expo
);
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
));
465 MPFR_SAVE_EXPO_FREE (expo
);
470 mpfr_cmp_f (mpfr_srcptr x
, mpf_srcptr z
)
474 MPFR_SAVE_EXPO_DECL (expo
);
476 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
477 return mpfr_cmp_si (x
, mpf_sgn (z
));
479 MPFR_SAVE_EXPO_MARK (expo
);
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
);
487 MPFR_SAVE_EXPO_FREE (expo
);