1/* Implementations of operations between mpfr and mpz/mpq data
2
3Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel 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.
27   This function should be called in the extended exponent range. */
28static void
29init_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  */
49static int
50foo (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
66static int
67foo2 (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
83int
84mpfr_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
89int
90mpfr_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
95int
96mpfr_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
105int
106mpfr_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
115int
116mpfr_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
125int
126mpfr_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. */
163static int
164mpfr_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
236int
237mpfr_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
242int
243mpfr_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
248int
249mpfr_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
338int
339mpfr_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
434int
435mpfr_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
469int
470mpfr_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