1 /* Implementations of operations between mpfr and mpz/mpq data
3 Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 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 */
28 init_set_z (mpfr_ptr t
, mpz_srcptr z
)
33 if (mpz_size (z
) <= 1)
36 MPFR_MPZ_SIZEINBASE2 (p
, z
);
38 i
= mpfr_set_z (t
, z
, GMP_RNDN
);
39 MPFR_ASSERTD (i
== 0);
42 /* Init, set a mpfr_t with enough precision to store a mpz_t without round,
43 call the function, and clear the allocated mpfr_t */
45 foo (mpfr_ptr x
, mpfr_srcptr y
, mpz_srcptr z
, mp_rnd_t r
,
46 int (*f
)(mpfr_ptr
, mpfr_srcptr
, mpfr_srcptr
, mp_rnd_t
))
51 i
= (*f
) (x
, y
, t
, r
);
57 mpfr_mul_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mp_rnd_t r
)
59 return foo (y
, x
, z
, r
, mpfr_mul
);
63 mpfr_div_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mp_rnd_t r
)
65 return foo (y
, x
, z
, r
, mpfr_div
);
69 mpfr_add_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mp_rnd_t r
)
71 /* Mpz 0 is unsigned */
72 if (MPFR_UNLIKELY (mpz_sgn (z
) == 0))
73 return mpfr_set (y
, x
, r
);
75 return foo (y
, x
, z
, r
, mpfr_add
);
79 mpfr_sub_z (mpfr_ptr y
, mpfr_srcptr x
, mpz_srcptr z
, mp_rnd_t r
)
81 /* Mpz 0 is unsigned */
82 if (MPFR_UNLIKELY (mpz_sgn (z
) == 0))
83 return mpfr_set (y
, x
, r
);
85 return foo (y
, x
, z
, r
, mpfr_sub
);
89 mpfr_cmp_z (mpfr_srcptr x
, mpz_srcptr z
)
94 res
= mpfr_cmp (x
, t
);
100 mpfr_mul_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
, mp_rnd_t rnd_mode
)
106 if (MPFR_UNLIKELY (mpq_sgn (z
) == 0))
107 return mpfr_mul_ui (y
, x
, 0, rnd_mode
);
110 MPFR_MPZ_SIZEINBASE2 (p
, mpq_numref (z
));
111 mpfr_init2 (tmp
, MPFR_PREC (x
) + p
);
112 res
= mpfr_mul_z (tmp
, x
, mpq_numref(z
), GMP_RNDN
);
113 MPFR_ASSERTD (res
== 0);
114 res
= mpfr_div_z (y
, tmp
, mpq_denref(z
), rnd_mode
);
121 mpfr_div_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
, mp_rnd_t rnd_mode
)
127 if (MPFR_UNLIKELY (mpq_sgn (z
) == 0))
128 return mpfr_div_ui (y
, x
, 0, rnd_mode
);
129 else if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z
)) == 0))
132 MPFR_MPZ_SIZEINBASE2 (p
, mpq_denref (z
));
133 mpfr_init2 (tmp
, MPFR_PREC(x
) + p
);
134 res
= mpfr_mul_z (tmp
, x
, mpq_denref(z
), GMP_RNDN
);
135 MPFR_ASSERTD( res
== 0 );
136 res
= mpfr_div_z (y
, tmp
, mpq_numref(z
), rnd_mode
);
142 mpfr_add_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
, mp_rnd_t rnd_mode
)
148 MPFR_ZIV_DECL (loop
);
150 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
157 else if (MPFR_IS_INF (x
))
159 MPFR_ASSERTD (mpz_sgn (mpq_denref (z
)) != 0);
161 MPFR_SET_SAME_SIGN (y
, x
);
166 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
167 if (MPFR_UNLIKELY (mpq_sgn (z
) == 0))
168 return mpfr_set (y
, x
, rnd_mode
); /* signed 0 - Unsigned 0 */
170 return mpfr_set_q (y
, z
, rnd_mode
);
174 p
= MPFR_PREC (y
) + 10;
178 MPFR_ZIV_INIT (loop
, p
);
181 res
= mpfr_set_q (q
, z
, GMP_RNDN
); /* Error <= 1/2 ulp(q) */
182 /* If z if @INF@ (1/0), res = 0, so it quits immediately */
183 if (MPFR_UNLIKELY (res
== 0))
184 /* Result is exact so we can add it directly! */
186 res
= mpfr_add (y
, x
, q
, rnd_mode
);
189 mpfr_add (t
, x
, q
, GMP_RNDN
); /* Error <= 1/2 ulp(t) */
190 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
191 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
193 If EXP(q)-EXP(t)<0, <= 2^0 */
194 /* We can get 0, but we can't round since q is inexact */
195 if (MPFR_LIKELY (!MPFR_IS_ZERO (t
)))
197 err
= (mp_exp_t
) p
- 1 - MAX (MPFR_GET_EXP(q
)-MPFR_GET_EXP(t
), 0);
198 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, err
, MPFR_PREC (y
), rnd_mode
)))
200 res
= mpfr_set (y
, t
, rnd_mode
);
204 MPFR_ZIV_NEXT (loop
, p
);
205 mpfr_set_prec (t
, p
);
206 mpfr_set_prec (q
, p
);
208 MPFR_ZIV_FREE (loop
);
215 mpfr_sub_q (mpfr_ptr y
, mpfr_srcptr x
, mpq_srcptr z
,mp_rnd_t rnd_mode
)
221 MPFR_ZIV_DECL (loop
);
223 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
230 else if (MPFR_IS_INF (x
))
232 MPFR_ASSERTD (mpz_sgn (mpq_denref (z
)) != 0);
234 MPFR_SET_SAME_SIGN (y
, x
);
239 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
241 if (MPFR_UNLIKELY (mpq_sgn (z
) == 0))
242 return mpfr_set (y
, x
, rnd_mode
); /* signed 0 - Unsigned 0 */
245 res
= mpfr_set_q (y
, z
, MPFR_INVERT_RND (rnd_mode
));
246 MPFR_CHANGE_SIGN (y
);
252 p
= MPFR_PREC (y
) + 10;
256 MPFR_ZIV_INIT (loop
, p
);
259 res
= mpfr_set_q(q
, z
, GMP_RNDN
); /* Error <= 1/2 ulp(q) */
260 /* If z if @INF@ (1/0), res = 0, so it quits immediately */
261 if (MPFR_UNLIKELY (res
== 0))
262 /* Result is exact so we can add it directly!*/
264 res
= mpfr_sub (y
, x
, q
, rnd_mode
);
267 mpfr_sub (t
, x
, q
, GMP_RNDN
); /* Error <= 1/2 ulp(t) */
268 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
269 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
271 If EXP(q)-EXP(t)<0, <= 2^0 */
272 /* We can get 0, but we can't round since q is inexact */
273 if (MPFR_LIKELY (!MPFR_IS_ZERO (t
)))
275 err
= (mp_exp_t
) p
- 1 - MAX (MPFR_GET_EXP(q
)-MPFR_GET_EXP(t
), 0);
276 res
= MPFR_CAN_ROUND (t
, err
, MPFR_PREC (y
), rnd_mode
);
277 if (MPFR_LIKELY (res
!= 0)) /* We can round! */
279 res
= mpfr_set (y
, t
, rnd_mode
);
283 MPFR_ZIV_NEXT (loop
, p
);
284 mpfr_set_prec (t
, p
);
285 mpfr_set_prec (q
, p
);
287 MPFR_ZIV_FREE (loop
);
294 mpfr_cmp_q (mpfr_srcptr x
, mpq_srcptr z
)
299 /* x < a/b ? <=> x*b < a */
300 MPFR_ASSERTD (mpz_sgn (mpq_denref (z
)) != 0);
301 MPFR_MPZ_SIZEINBASE2 (p
, mpq_denref (z
));
302 mpfr_init2 (t
, MPFR_PREC(x
) + p
);
303 res
= mpfr_mul_z (t
, x
, mpq_denref (z
), GMP_RNDN
);
304 MPFR_ASSERTD (res
== 0);
305 res
= mpfr_cmp_z (t
, mpq_numref (z
) );
311 mpfr_cmp_f (mpfr_srcptr x
, mpf_srcptr z
)
316 mpfr_init2 (t
, MPFR_PREC_MIN
+ ABS(SIZ(z
)) * BITS_PER_MP_LIMB
);
317 res
= mpfr_set_f (t
, z
, GMP_RNDN
);
318 MPFR_ASSERTD (res
== 0);
319 res
= mpfr_cmp (x
, t
);