1/* Implementations of operations between mpfr and mpz/mpq data
2
3Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao 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/* Init and set a mpfr_t with enough precision to store a mpz */
27static void
28init_set_z (mpfr_ptr t, mpz_srcptr z)
29{
30  mpfr_prec_t p;
31  int i;
32
33  if (mpz_size (z) <= 1)
34    p = GMP_NUMB_BITS;
35  else
36    MPFR_MPZ_SIZEINBASE2 (p, z);
37  mpfr_init2 (t, p);
38  i = mpfr_set_z (t, z, MPFR_RNDN);
39  MPFR_ASSERTD (i == 0);
40}
41
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  */
44static int
45foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mpfr_rnd_t r,
46     int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t))
47{
48  mpfr_t t;
49  int i;
50  init_set_z (t, z);
51  i = (*f) (x, y, t, r);
52  mpfr_clear (t);
53  return i;
54}
55
56int
57mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
58{
59  return foo (y, x, z, r, mpfr_mul);
60}
61
62int
63mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
64{
65  return foo (y, x, z, r, mpfr_div);
66}
67
68int
69mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
70{
71  /* Mpz 0 is unsigned */
72  if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
73    return mpfr_set (y, x, r);
74  else
75    return foo (y, x, z, r, mpfr_add);
76}
77
78int
79mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r)
80{
81  /* Mpz 0 is unsigned */
82  if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
83    return mpfr_set (y, x, r);
84  else
85    return foo (y, x, z, r, mpfr_sub);
86}
87
88int
89mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z)
90{
91  mpfr_t t;
92  int res;
93  init_set_z (t, z);
94  res = mpfr_cmp (x, t);
95  mpfr_clear (t);
96  return res;
97}
98
99int
100mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode)
101{
102  mpfr_t tmp;
103  int res;
104  mpfr_prec_t p;
105
106  if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
107    return mpfr_mul_ui (y, x, 0, rnd_mode);
108  else
109    {
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), MPFR_RNDN );
113      MPFR_ASSERTD (res == 0);
114      res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode);
115      mpfr_clear (tmp);
116      return res;
117    }
118}
119
120int
121mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode)
122{
123  mpfr_t tmp;
124  int res;
125  mpfr_prec_t p;
126
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))
130    p = 0;
131  else
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), MPFR_RNDN );
135  MPFR_ASSERTD( res == 0 );
136  res = mpfr_div_z (y, tmp, mpq_numref(z), rnd_mode);
137  mpfr_clear (tmp);
138  return res;
139}
140
141int
142mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode)
143{
144  mpfr_t      t,q;
145  mpfr_prec_t p;
146  mpfr_exp_t  err;
147  int res;
148  MPFR_ZIV_DECL (loop);
149
150  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
151    {
152      if (MPFR_IS_NAN (x))
153        {
154          MPFR_SET_NAN (y);
155          MPFR_RET_NAN;
156        }
157      else if (MPFR_IS_INF (x))
158        {
159          MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
160          MPFR_SET_INF (y);
161          MPFR_SET_SAME_SIGN (y, x);
162          MPFR_RET (0);
163        }
164      else
165        {
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 */
169          else
170            return mpfr_set_q (y, z, rnd_mode);
171        }
172    }
173
174  p = MPFR_PREC (y) + 10;
175  mpfr_init2 (t, p);
176  mpfr_init2 (q, p);
177
178  MPFR_ZIV_INIT (loop, p);
179  for (;;)
180    {
181      res = mpfr_set_q (q, z, MPFR_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! */
185        {
186          res = mpfr_add (y, x, q, rnd_mode);
187          break;
188        }
189      mpfr_add (t, x, q, MPFR_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)))
192                             <= 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)))
196        {
197          err = (mpfr_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)))
199            {
200              res = mpfr_set (y, t, rnd_mode);
201              break;
202            }
203        }
204      MPFR_ZIV_NEXT (loop, p);
205      mpfr_set_prec (t, p);
206      mpfr_set_prec (q, p);
207    }
208  MPFR_ZIV_FREE (loop);
209  mpfr_clear (t);
210  mpfr_clear (q);
211  return res;
212}
213
214int
215mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mpfr_rnd_t rnd_mode)
216{
217  mpfr_t t,q;
218  mpfr_prec_t p;
219  int res;
220  mpfr_exp_t err;
221  MPFR_ZIV_DECL (loop);
222
223  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
224    {
225      if (MPFR_IS_NAN (x))
226        {
227          MPFR_SET_NAN (y);
228          MPFR_RET_NAN;
229        }
230      else if (MPFR_IS_INF (x))
231        {
232          MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
233          MPFR_SET_INF (y);
234          MPFR_SET_SAME_SIGN (y, x);
235          MPFR_RET (0);
236        }
237      else
238        {
239          MPFR_ASSERTD (MPFR_IS_ZERO (x));
240
241          if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
242            return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
243          else
244            {
245              res =  mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode));
246              MPFR_CHANGE_SIGN (y);
247              return -res;
248            }
249        }
250    }
251
252  p = MPFR_PREC (y) + 10;
253  mpfr_init2 (t, p);
254  mpfr_init2 (q, p);
255
256  MPFR_ZIV_INIT (loop, p);
257  for(;;)
258    {
259      res = mpfr_set_q(q, z, MPFR_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!*/
263        {
264          res = mpfr_sub (y, x, q, rnd_mode);
265          break;
266        }
267      mpfr_sub (t, x, q, MPFR_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)))
270                             <= 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)))
274        {
275          err = (mpfr_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! */
278            {
279              res = mpfr_set (y, t, rnd_mode);
280              break;
281            }
282        }
283      MPFR_ZIV_NEXT (loop, p);
284      mpfr_set_prec (t, p);
285      mpfr_set_prec (q, p);
286    }
287  MPFR_ZIV_FREE (loop);
288  mpfr_clear (t);
289  mpfr_clear (q);
290  return res;
291}
292
293int
294mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr z)
295{
296  mpfr_t t;
297  int res;
298  mpfr_prec_t p;
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), MPFR_RNDN );
304  MPFR_ASSERTD (res == 0);
305  res = mpfr_cmp_z (t, mpq_numref (z) );
306  mpfr_clear (t);
307  return res;
308}
309
310int
311mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z)
312{
313  mpfr_t t;
314  int res;
315
316  mpfr_init2 (t, MPFR_PREC_MIN + ABS(SIZ(z)) * GMP_NUMB_BITS );
317  res = mpfr_set_f (t, z, MPFR_RNDN);
318  MPFR_ASSERTD (res == 0);
319  res = mpfr_cmp (x, t);
320  mpfr_clear (t);
321  return res;
322}
323