contrib/mpfr/src/gmp_op.c

   1 /* Implementations of operations between mpfr and mpz/mpq data
   2
   3 Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   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.
  12
  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.
  17
  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. */
  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. */
  28 static void
  29 init_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  */
  49 static int
  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))
  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
  66 static int
  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))
  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
  83 int
  84 mpfr_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
  89 int
  90 mpfr_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
  95 int
  96 mpfr_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
 105 int
 106 mpfr_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
 115 int
 116 mpfr_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
 125 int
 126 mpfr_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. */
 163 static int
 164 mpfr_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
 236 int
 237 mpfr_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
 242 int
 243 mpfr_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
 248 int
 249 mpfr_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
 338 int
 339 mpfr_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
 434 int
 435 mpfr_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
 469 int
 470 mpfr_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 }