contrib/mpfr/pow.c

   1 /* mpfr_pow -- power function x^y
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* return non zero iff x^y is exact.
  27    Assumes x and y are ordinary numbers,
  28    y is not an integer, x is not a power of 2 and x is positive
  29
  30    If x^y is exact, it computes it and sets *inexact.
  31 */
  32 static int
  33 mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
  34                    mp_rnd_t rnd_mode, int *inexact)
  35 {
  36   mpz_t a, c;
  37   mp_exp_t d, b;
  38   unsigned long i;
  39   int res;
  40
  41   MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));
  42   MPFR_ASSERTD (!MPFR_IS_SINGULAR (x));
  43   MPFR_ASSERTD (!mpfr_integer_p (y));
  44   MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_INT_SIGN (x),
  45                                   MPFR_GET_EXP (x) - 1) != 0);
  46   MPFR_ASSERTD (MPFR_IS_POS (x));
  47
  48   if (MPFR_IS_NEG (y))
  49     return 0; /* x is not a power of two => x^-y is not exact */
  50
  51   /* compute d such that y = c*2^d with c odd integer */
  52   mpz_init (c);
  53   d = mpfr_get_z_exp (c, y);
  54   i = mpz_scan1 (c, 0);
  55   mpz_div_2exp (c, c, i);
  56   d += i;
  57   /* now y=c*2^d with c odd */
  58   /* Since y is not an integer, d is necessarily < 0 */
  59   MPFR_ASSERTD (d < 0);
  60
  61   /* Compute a,b such that x=a*2^b */
  62   mpz_init (a);
  63   b = mpfr_get_z_exp (a, x);
  64   i = mpz_scan1 (a, 0);
  65   mpz_div_2exp (a, a, i);
  66   b += i;
  67   /* now x=a*2^b with a is odd */
  68
  69   for (res = 1 ; d != 0 ; d++)
  70     {
  71       /* a*2^b is a square iff
  72             (i)  a is a square when b is even
  73             (ii) 2*a is a square when b is odd */
  74       if (b % 2 != 0)
  75         {
  76           mpz_mul_2exp (a, a, 1); /* 2*a */
  77           b --;
  78         }
  79       MPFR_ASSERTD ((b % 2) == 0);
  80       if (!mpz_perfect_square_p (a))
  81         {
  82           res = 0;
  83           goto end;
  84         }
  85       mpz_sqrt (a, a);
  86       b = b / 2;
  87     }
  88   /* Now x = (a'*2^b')^(2^-d) with d < 0
  89      so x^y = ((a'*2^b')^(2^-d))^(c*2^d)
  90             = ((a'*2^b')^c with c odd integer */
  91   {
  92     mpfr_t tmp;
  93     mp_prec_t p;
  94     MPFR_MPZ_SIZEINBASE2 (p, a);
  95     mpfr_init2 (tmp, p); /* prec = 1 should not be possible */
  96     res = mpfr_set_z (tmp, a, GMP_RNDN);
  97     MPFR_ASSERTD (res == 0);
  98     res = mpfr_mul_2si (tmp, tmp, b, GMP_RNDN);
  99     MPFR_ASSERTD (res == 0);
 100     *inexact = mpfr_pow_z (z, tmp, c, rnd_mode);
 101     mpfr_clear (tmp);
 102     res = 1;
 103   }
 104  end:
 105   mpz_clear (a);
 106   mpz_clear (c);
 107   return res;
 108 }
 109
 110 /* Return 1 if y is an odd integer, 0 otherwise. */
 111 static int
 112 is_odd (mpfr_srcptr y)
 113 {
 114   mp_exp_t expo;
 115   mp_prec_t prec;
 116   mp_size_t yn;
 117   mp_limb_t *yp;
 118
 119   /* NAN, INF or ZERO are not allowed */
 120   MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));
 121
 122   expo = MPFR_GET_EXP (y);
 123   if (expo <= 0)
 124     return 0;  /* |y| < 1 and not 0 */
 125
 126   prec = MPFR_PREC(y);
 127   if ((mpfr_prec_t) expo > prec)
 128     return 0;  /* y is a multiple of 2^(expo-prec), thus not odd */
 129
 130   /* 0 < expo <= prec:
 131      y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
 132           expo bits   (prec-expo) bits
 133
 134      We have to check that:
 135      (a) the bit 't' is set
 136      (b) all the 'z' bits are zero
 137   */
 138
 139   prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
 140   /* number of z+0 bits */
 141
 142   yn = prec / BITS_PER_MP_LIMB;
 143   MPFR_ASSERTN(yn >= 0);
 144   /* yn is the index of limb containing the 't' bit */
 145
 146   yp = MPFR_MANT(y);
 147   /* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */
 148   if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0
 149       : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT)
 150     return 0;
 151   while (--yn >= 0)
 152     if (yp[yn] != 0)
 153       return 0;
 154   return 1;
 155 }
 156
 157 /* Assumes that the exponent range has already been extended and if y is
 158    an integer, then the result is not exact in unbounded exponent range. */
 159 int
 160 mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
 161                   mp_rnd_t rnd_mode, int y_is_integer, mpfr_save_expo_t *expo)
 162 {
 163   mpfr_t t, u, k, absx;
 164   int k_non_zero = 0;
 165   int check_exact_case = 0;
 166   int inexact;
 167   /* Declaration of the size variable */
 168   mp_prec_t Nz = MPFR_PREC(z);               /* target precision */
 169   mp_prec_t Nt;                              /* working precision */
 170   mp_exp_t err;                              /* error */
 171   MPFR_ZIV_DECL (ziv_loop);
 172
 173
 174   MPFR_LOG_FUNC (("x[%#R]=%R y[%#R]=%R rnd=%d", x, x, y, y, rnd_mode),
 175                  ("z[%#R]=%R inexact=%d", z, z, inexact));
 176
 177   /* We put the absolute value of x in absx, pointing to the significand
 178      of x to avoid allocating memory for the significand of absx. */
 179   MPFR_ALIAS(absx, x, /*sign=*/ 1, /*EXP=*/ MPFR_EXP(x));
 180
 181   /* We will compute the absolute value of the result. So, let's
 182      invert the rounding mode if the result is negative. */
 183   if (MPFR_IS_NEG (x) && is_odd (y))
 184     rnd_mode = MPFR_INVERT_RND (rnd_mode);
 185
 186   /* compute the precision of intermediary variable */
 187   /* the optimal number of bits : see algorithms.tex */
 188   Nt = Nz + 5 + MPFR_INT_CEIL_LOG2 (Nz);
 189
 190   /* initialise of intermediary variable */
 191   mpfr_init2 (t, Nt);
 192
 193   MPFR_ZIV_INIT (ziv_loop, Nt);
 194   for (;;)
 195     {
 196       MPFR_BLOCK_DECL (flags1);
 197
 198       /* compute exp(y*ln|x|), using GMP_RNDU to get an upper bound, so
 199          that we can detect underflows. */
 200       mpfr_log (t, absx, MPFR_IS_NEG (y) ? GMP_RNDD : GMP_RNDU); /* ln|x| */
 201       mpfr_mul (t, y, t, GMP_RNDU);                              /* y*ln|x| */
 202       if (k_non_zero)
 203         {
 204           MPFR_LOG_MSG (("subtract k * ln(2)\n", 0));
 205           mpfr_const_log2 (u, GMP_RNDD);
 206           mpfr_mul (u, u, k, GMP_RNDD);
 207           /* Error on u = k * log(2): < k * 2^(-Nt) < 1. */
 208           mpfr_sub (t, t, u, GMP_RNDU);
 209           MPFR_LOG_MSG (("t = y * ln|x| - k * ln(2)\n", 0));
 210           MPFR_LOG_VAR (t);
 211         }
 212       /* estimate of the error -- see pow function in algorithms.tex.
 213          The error on t is at most 1/2 + 3*2^(EXP(t)+1) ulps, which is
 214          <= 2^(EXP(t)+3) for EXP(t) >= -1, and <= 2 ulps for EXP(t) <= -2.
 215          Additional error if k_no_zero: treal = t * errk, with
 216          1 - |k| * 2^(-Nt) <= exp(-|k| * 2^(-Nt)) <= errk <= 1,
 217          i.e., additional absolute error <= 2^(EXP(k)+EXP(t)-Nt).
 218          Total error <= 2^err1 + 2^err2 <= 2^(max(err1,err2)+1). */
 219       err = MPFR_NOTZERO (t) && MPFR_GET_EXP (t) >= -1 ?
 220         MPFR_GET_EXP (t) + 3 : 1;
 221       if (k_non_zero)
 222         {
 223           if (MPFR_GET_EXP (k) > err)
 224             err = MPFR_GET_EXP (k);
 225           err++;
 226         }
 227       MPFR_BLOCK (flags1, mpfr_exp (t, t, GMP_RNDN));  /* exp(y*ln|x|)*/
 228       /* We need to test */
 229       if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (t) || MPFR_UNDERFLOW (flags1)))
 230         {
 231           mp_prec_t Ntmin;
 232           MPFR_BLOCK_DECL (flags2);
 233
 234           MPFR_ASSERTN (!k_non_zero);
 235           MPFR_ASSERTN (!MPFR_IS_NAN (t));
 236
 237           /* Real underflow? */
 238           if (MPFR_IS_ZERO (t))
 239             {
 240               /* Underflow. We computed rndn(exp(t)), where t >= y*ln|x|.
 241                  Therefore rndn(|x|^y) = 0, and we have a real underflow on
 242                  |x|^y. */
 243               inexact = mpfr_underflow (z, rnd_mode == GMP_RNDN ? GMP_RNDZ
 244                                         : rnd_mode, MPFR_SIGN_POS);
 245               if (expo != NULL)
 246                 MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
 247                                              | MPFR_FLAGS_UNDERFLOW);
 248               break;
 249             }
 250
 251           /* Real overflow? */
 252           if (MPFR_IS_INF (t))
 253             {
 254               /* Note: we can probably use a low precision for this test. */
 255               mpfr_log (t, absx, MPFR_IS_NEG (y) ? GMP_RNDU : GMP_RNDD);
 256               mpfr_mul (t, y, t, GMP_RNDD);            /* y * ln|x| */
 257               MPFR_BLOCK (flags2, mpfr_exp (t, t, GMP_RNDD));
 258               /* t = lower bound on exp(y * ln|x|) */
 259               if (MPFR_OVERFLOW (flags2))
 260                 {
 261                   /* We have computed a lower bound on |x|^y, and it
 262                      overflowed. Therefore we have a real overflow
 263                      on |x|^y. */
 264                   inexact = mpfr_overflow (z, rnd_mode, MPFR_SIGN_POS);
 265                   if (expo != NULL)
 266                     MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
 267                                                  | MPFR_FLAGS_OVERFLOW);
 268                   break;
 269                 }
 270             }
 271
 272           k_non_zero = 1;
 273           Ntmin = sizeof(mp_exp_t) * CHAR_BIT;
 274           if (Ntmin > Nt)
 275             {
 276               Nt = Ntmin;
 277               mpfr_set_prec (t, Nt);
 278             }
 279           mpfr_init2 (u, Nt);
 280           mpfr_init2 (k, Ntmin);
 281           mpfr_log2 (k, absx, GMP_RNDN);
 282           mpfr_mul (k, y, k, GMP_RNDN);
 283           mpfr_round (k, k);
 284           MPFR_LOG_VAR (k);
 285           /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */
 286           continue;
 287         }
 288       if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Nz, rnd_mode)))
 289         {
 290           inexact = mpfr_set (z, t, rnd_mode);
 291           break;
 292         }
 293
 294       /* check exact power, except when y is an integer (since the
 295          exact cases for y integer have already been filtered out) */
 296       if (check_exact_case == 0 && ! y_is_integer)
 297         {
 298           if (mpfr_pow_is_exact (z, absx, y, rnd_mode, &inexact))
 299             break;
 300           check_exact_case = 1;
 301         }
 302
 303       /* reactualisation of the precision */
 304       MPFR_ZIV_NEXT (ziv_loop, Nt);
 305       mpfr_set_prec (t, Nt);
 306       if (k_non_zero)
 307         mpfr_set_prec (u, Nt);
 308     }
 309   MPFR_ZIV_FREE (ziv_loop);
 310
 311   if (k_non_zero)
 312     {
 313       int inex2;
 314       long lk;
 315
 316       /* The rounded result in an unbounded exponent range is z * 2^k. As
 317        * MPFR chooses underflow after rounding, the mpfr_mul_2si below will
 318        * correctly detect underflows and overflows. However, in rounding to
 319        * nearest, if z * 2^k = 2^(emin - 2), then the double rounding may
 320        * affect the result. We need to cope with that before overwriting z.
 321        * If inexact >= 0, then the real result is <= 2^(emin - 2), so that
 322        * o(2^(emin - 2)) = +0 is correct. If inexact < 0, then the real
 323        * result is > 2^(emin - 2) and we need to round to 2^(emin - 1).
 324        */
 325       MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
 326       lk = mpfr_get_si (k, GMP_RNDN);
 327       if (rnd_mode == GMP_RNDN && inexact < 0 &&
 328           MPFR_GET_EXP (z) + lk == __gmpfr_emin - 1 && mpfr_powerof2_raw (z))
 329         {
 330           /* Rounding to nearest, real result > z * 2^k = 2^(emin - 2),
 331            * underflow case: as the minimum precision is > 1, we will
 332            * obtain the correct result and exceptions by replacing z by
 333            * nextabove(z).
 334            */
 335           MPFR_ASSERTN (MPFR_PREC_MIN > 1);
 336           mpfr_nextabove (z);
 337         }
 338       mpfr_clear_flags ();
 339       inex2 = mpfr_mul_2si (z, z, lk, rnd_mode);
 340       if (inex2)  /* underflow or overflow */
 341         {
 342           inexact = inex2;
 343           if (expo != NULL)
 344             MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, __gmpfr_flags);
 345         }
 346       mpfr_clears (u, k, (mpfr_ptr) 0);
 347     }
 348   mpfr_clear (t);
 349
 350   /* update the sign of the result if x was negative */
 351   if (MPFR_IS_NEG (x) && is_odd (y))
 352     {
 353       MPFR_SET_NEG(z);
 354       inexact = -inexact;
 355     }
 356
 357   return inexact;
 358 }
 359
 360 /* The computation of z = pow(x,y) is done by
 361    z = exp(y * log(x)) = x^y
 362    For the special cases, see Section F.9.4.4 of the C standard:
 363      _ pow(±0, y) = ±inf for y an odd integer < 0.
 364      _ pow(±0, y) = +inf for y < 0 and not an odd integer.
 365      _ pow(±0, y) = ±0 for y an odd integer > 0.
 366      _ pow(±0, y) = +0 for y > 0 and not an odd integer.
 367      _ pow(-1, ±inf) = 1.
 368      _ pow(+1, y) = 1 for any y, even a NaN.
 369      _ pow(x, ±0) = 1 for any x, even a NaN.
 370      _ pow(x, y) = NaN for finite x < 0 and finite non-integer y.
 371      _ pow(x, -inf) = +inf for |x| < 1.
 372      _ pow(x, -inf) = +0 for |x| > 1.
 373      _ pow(x, +inf) = +0 for |x| < 1.
 374      _ pow(x, +inf) = +inf for |x| > 1.
 375      _ pow(-inf, y) = -0 for y an odd integer < 0.
 376      _ pow(-inf, y) = +0 for y < 0 and not an odd integer.
 377      _ pow(-inf, y) = -inf for y an odd integer > 0.
 378      _ pow(-inf, y) = +inf for y > 0 and not an odd integer.
 379      _ pow(+inf, y) = +0 for y < 0.
 380      _ pow(+inf, y) = +inf for y > 0. */
 381 int
 382 mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 383 {
 384   int inexact;
 385   int cmp_x_1;
 386   int y_is_integer;
 387   MPFR_SAVE_EXPO_DECL (expo);
 388
 389   MPFR_LOG_FUNC (("x[%#R]=%R y[%#R]=%R rnd=%d", x, x, y, y, rnd_mode),
 390                  ("z[%#R]=%R inexact=%d", z, z, inexact));
 391
 392   if (MPFR_ARE_SINGULAR (x, y))
 393     {
 394       /* pow(x, 0) returns 1 for any x, even a NaN. */
 395       if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
 396         return mpfr_set_ui (z, 1, rnd_mode);
 397       else if (MPFR_IS_NAN (x))
 398         {
 399           MPFR_SET_NAN (z);
 400           MPFR_RET_NAN;
 401         }
 402       else if (MPFR_IS_NAN (y))
 403         {
 404           /* pow(+1, NaN) returns 1. */
 405           if (mpfr_cmp_ui (x, 1) == 0)
 406             return mpfr_set_ui (z, 1, rnd_mode);
 407           MPFR_SET_NAN (z);
 408           MPFR_RET_NAN;
 409         }
 410       else if (MPFR_IS_INF (y))
 411         {
 412           if (MPFR_IS_INF (x))
 413             {
 414               if (MPFR_IS_POS (y))
 415                 MPFR_SET_INF (z);
 416               else
 417                 MPFR_SET_ZERO (z);
 418               MPFR_SET_POS (z);
 419               MPFR_RET (0);
 420             }
 421           else
 422             {
 423               int cmp;
 424               cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y);
 425               MPFR_SET_POS (z);
 426               if (cmp > 0)
 427                 {
 428                   /* Return +inf. */
 429                   MPFR_SET_INF (z);
 430                   MPFR_RET (0);
 431                 }
 432               else if (cmp < 0)
 433                 {
 434                   /* Return +0. */
 435                   MPFR_SET_ZERO (z);
 436                   MPFR_RET (0);
 437                 }
 438               else
 439                 {
 440                   /* Return 1. */
 441                   return mpfr_set_ui (z, 1, rnd_mode);
 442                 }
 443             }
 444         }
 445       else if (MPFR_IS_INF (x))
 446         {
 447           int negative;
 448           /* Determine the sign now, in case y and z are the same object */
 449           negative = MPFR_IS_NEG (x) && is_odd (y);
 450           if (MPFR_IS_POS (y))
 451             MPFR_SET_INF (z);
 452           else
 453             MPFR_SET_ZERO (z);
 454           if (negative)
 455             MPFR_SET_NEG (z);
 456           else
 457             MPFR_SET_POS (z);
 458           MPFR_RET (0);
 459         }
 460       else
 461         {
 462           int negative;
 463           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 464           /* Determine the sign now, in case y and z are the same object */
 465           negative = MPFR_IS_NEG(x) && is_odd (y);
 466           if (MPFR_IS_NEG (y))
 467             MPFR_SET_INF (z);
 468           else
 469             MPFR_SET_ZERO (z);
 470           if (negative)
 471             MPFR_SET_NEG (z);
 472           else
 473             MPFR_SET_POS (z);
 474           MPFR_RET (0);
 475         }
 476     }
 477
 478   /* x^y for x < 0 and y not an integer is not defined */
 479   y_is_integer = mpfr_integer_p (y);
 480   if (MPFR_IS_NEG (x) && ! y_is_integer)
 481     {
 482       MPFR_SET_NAN (z);
 483       MPFR_RET_NAN;
 484     }
 485
 486   /* now the result cannot be NaN:
 487      (1) either x > 0
 488      (2) or x < 0 and y is an integer */
 489
 490   cmp_x_1 = mpfr_cmpabs (x, __gmpfr_one);
 491   if (cmp_x_1 == 0)
 492     return mpfr_set_si (z, MPFR_IS_NEG (x) && is_odd (y) ? -1 : 1, rnd_mode);
 493
 494   /* now we have:
 495      (1) either x > 0
 496      (2) or x < 0 and y is an integer
 497      and in addition |x| <> 1.
 498   */
 499
 500   /* detect overflow: an overflow is possible if
 501      (a) |x| > 1 and y > 0
 502      (b) |x| < 1 and y < 0.
 503      FIXME: this assumes 1 is always representable.
 504
 505      FIXME2: maybe we can test overflow and underflow simultaneously.
 506      The idea is the following: first compute an approximation to
 507      y * log2|x|, using rounding to nearest. If |x| is not too near from 1,
 508      this approximation should be accurate enough, and in most cases enable
 509      one to prove that there is no underflow nor overflow.
 510      Otherwise, it should enable one to check only underflow or overflow,
 511      instead of both cases as in the present case.
 512   */
 513   if (cmp_x_1 * MPFR_SIGN (y) > 0)
 514     {
 515       mpfr_t t;
 516       int negative, overflow;
 517
 518       MPFR_SAVE_EXPO_MARK (expo);
 519       mpfr_init2 (t, 53);
 520       /* we want a lower bound on y*log2|x|:
 521          (i) if x > 0, it suffices to round log2(x) towards zero, and
 522              to round y*o(log2(x)) towards zero too;
 523          (ii) if x < 0, we first compute t = o(-x), with rounding towards 1,
 524               and then follow as in case (1). */
 525       if (MPFR_SIGN (x) > 0)
 526         mpfr_log2 (t, x, GMP_RNDZ);
 527       else
 528         {
 529           mpfr_neg (t, x, (cmp_x_1 > 0) ? GMP_RNDZ : GMP_RNDU);
 530           mpfr_log2 (t, t, GMP_RNDZ);
 531         }
 532       mpfr_mul (t, t, y, GMP_RNDZ);
 533       overflow = mpfr_cmp_si (t, __gmpfr_emax) > 0;
 534       mpfr_clear (t);
 535       MPFR_SAVE_EXPO_FREE (expo);
 536       if (overflow)
 537         {
 538           MPFR_LOG_MSG (("early overflow detection\n", 0));
 539           negative = MPFR_SIGN(x) < 0 && is_odd (y);
 540           return mpfr_overflow (z, rnd_mode, negative ? -1 : 1);
 541         }
 542     }
 543
 544   /* Basic underflow checking. One has:
 545    *   - if y > 0, |x^y| < 2^(EXP(x) * y);
 546    *   - if y < 0, |x^y| <= 2^((EXP(x) - 1) * y);
 547    * so that one can compute a value ebound such that |x^y| < 2^ebound.
 548    * If we have ebound <= emin - 2 (emin - 1 in directed rounding modes),
 549    * then there is an underflow and we can decide the return value.
 550    */
 551   if (MPFR_IS_NEG (y) ? (MPFR_GET_EXP (x) > 1) : (MPFR_GET_EXP (x) < 0))
 552     {
 553       mpfr_t tmp;
 554       mpfr_eexp_t ebound;
 555       int inex2;
 556
 557       /* We must restore the flags. */
 558       MPFR_SAVE_EXPO_MARK (expo);
 559       mpfr_init2 (tmp, sizeof (mp_exp_t) * CHAR_BIT);
 560       inex2 = mpfr_set_exp_t (tmp, MPFR_GET_EXP (x), GMP_RNDN);
 561       MPFR_ASSERTN (inex2 == 0);
 562       if (MPFR_IS_NEG (y))
 563         {
 564           inex2 = mpfr_sub_ui (tmp, tmp, 1, GMP_RNDN);
 565           MPFR_ASSERTN (inex2 == 0);
 566         }
 567       mpfr_mul (tmp, tmp, y, GMP_RNDU);
 568       if (MPFR_IS_NEG (y))
 569         mpfr_nextabove (tmp);
 570       /* tmp doesn't necessarily fit in ebound, but that doesn't matter
 571          since we get the minimum value in such a case. */
 572       ebound = mpfr_get_exp_t (tmp, GMP_RNDU);
 573       mpfr_clear (tmp);
 574       MPFR_SAVE_EXPO_FREE (expo);
 575       if (MPFR_UNLIKELY (ebound <=
 576                          __gmpfr_emin - (rnd_mode == GMP_RNDN ? 2 : 1)))
 577         {
 578           /* warning: mpfr_underflow rounds away from 0 for GMP_RNDN */
 579           MPFR_LOG_MSG (("early underflow detection\n", 0));
 580           return mpfr_underflow (z,
 581                                  rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,
 582                                  MPFR_SIGN (x) < 0 && is_odd (y) ? -1 : 1);
 583         }
 584     }
 585
 586   /* If y is an integer, we can use mpfr_pow_z (based on multiplications),
 587      but if y is very large (I'm not sure about the best threshold -- VL),
 588      we shouldn't use it, as it can be very slow and take a lot of memory
 589      (and even crash or make other programs crash, as several hundred of
 590      MBs may be necessary). Note that in such a case, either x = +/-2^b
 591      (this case is handled below) or x^y cannot be represented exactly in
 592      any precision supported by MPFR (the general case uses this property).
 593   */
 594   if (y_is_integer && (MPFR_GET_EXP (y) <= 256))
 595     {
 596       mpz_t zi;
 597
 598       MPFR_LOG_MSG (("special code for y not too large integer\n", 0));
 599       mpz_init (zi);
 600       mpfr_get_z (zi, y, GMP_RNDN);
 601       inexact = mpfr_pow_z (z, x, zi, rnd_mode);
 602       mpz_clear (zi);
 603       return inexact;
 604     }
 605
 606   /* Special case (+/-2^b)^Y which could be exact. If x is negative, then
 607      necessarily y is a large integer. */
 608   {
 609     mp_exp_t b = MPFR_GET_EXP (x) - 1;
 610
 611     MPFR_ASSERTN (b >= LONG_MIN && b <= LONG_MAX);  /* FIXME... */
 612     if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), b) == 0)
 613       {
 614         mpfr_t tmp;
 615         int sgnx = MPFR_SIGN (x);
 616
 617         MPFR_LOG_MSG (("special case (+/-2^b)^Y\n", 0));
 618         /* now x = +/-2^b, so x^y = (+/-1)^y*2^(b*y) is exact whenever b*y is
 619            an integer */
 620         MPFR_SAVE_EXPO_MARK (expo);
 621         mpfr_init2 (tmp, MPFR_PREC (y) + sizeof (long) * CHAR_BIT);
 622         inexact = mpfr_mul_si (tmp, y, b, GMP_RNDN); /* exact */
 623         MPFR_ASSERTN (inexact == 0);
 624         /* Note: as the exponent range has been extended, an overflow is not
 625            possible (due to basic overflow and underflow checking above, as
 626            the result is ~ 2^tmp), and an underflow is not possible either
 627            because b is an integer (thus either 0 or >= 1). */
 628         mpfr_clear_flags ();
 629         inexact = mpfr_exp2 (z, tmp, rnd_mode);
 630         mpfr_clear (tmp);
 631         if (sgnx < 0 && is_odd (y))
 632           {
 633             mpfr_neg (z, z, rnd_mode);
 634             inexact = -inexact;
 635           }
 636         /* Without the following, the overflows3 test in tpow.c fails. */
 637         MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 638         MPFR_SAVE_EXPO_FREE (expo);
 639         return mpfr_check_range (z, inexact, rnd_mode);
 640       }
 641   }
 642
 643   MPFR_SAVE_EXPO_MARK (expo);
 644
 645   /* Case where |y * log(x)| is very small. Warning: x can be negative, in
 646      that case y is a large integer. */
 647   {
 648     mpfr_t t;
 649     mp_exp_t err;
 650
 651     /* We need an upper bound on the exponent of y * log(x). */
 652     mpfr_init2 (t, 16);
 653     if (MPFR_IS_POS(x))
 654       mpfr_log (t, x, cmp_x_1 < 0 ? GMP_RNDD : GMP_RNDU); /* away from 0 */
 655     else
 656       {
 657         /* if x < -1, round to +Inf, else round to zero */
 658         mpfr_neg (t, x, (mpfr_cmp_si (x, -1) < 0) ? GMP_RNDU : GMP_RNDD);
 659         mpfr_log (t, t, (mpfr_cmp_ui (t, 1) < 0) ? GMP_RNDD : GMP_RNDU);
 660       }
 661     MPFR_ASSERTN (MPFR_IS_PURE_FP (t));
 662     err = MPFR_GET_EXP (y) + MPFR_GET_EXP (t);
 663     mpfr_clear (t);
 664     mpfr_clear_flags ();
 665     MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (z, __gmpfr_one, - err, 0,
 666                                       (MPFR_SIGN (y) > 0) ^ (cmp_x_1 < 0),
 667                                       rnd_mode, expo, {});
 668   }
 669
 670   /* General case */
 671   inexact = mpfr_pow_general (z, x, y, rnd_mode, y_is_integer, &expo);
 672
 673   MPFR_SAVE_EXPO_FREE (expo);
 674   return mpfr_check_range (z, inexact, rnd_mode);
 675 }