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