contrib/mpfr/src/pow.c

   1 /* mpfr_pow -- power function x^y
   2
   3 Copyright 2001, 2002, 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 /* 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                    mpfr_rnd_t rnd_mode, int *inexact)
  35 {
  36   mpz_t a, c;
  37   mpfr_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_2exp (c, y);
  54   i = mpz_scan1 (c, 0);
  55   mpz_fdiv_q_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_2exp (a, x);
  64   i = mpz_scan1 (a, 0);
  65   mpz_fdiv_q_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     mpfr_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, MPFR_RNDN);
  97     MPFR_ASSERTD (res == 0);
  98     res = mpfr_mul_2si (tmp, tmp, b, MPFR_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   mpfr_exp_t expo;
 115   mpfr_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) / GMP_NUMB_BITS + 1) * GMP_NUMB_BITS - expo;
 140   /* number of z+0 bits */
 141
 142   yn = prec / GMP_NUMB_BITS;
 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 GMP_NUMB_BITS, t is bit 0 */
 148   if (expo % GMP_NUMB_BITS == 0 ? (yp[yn] & 1) == 0
 149       : yp[yn] << ((expo % GMP_NUMB_BITS) - 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                   mpfr_rnd_t rnd_mode, int y_is_integer, mpfr_save_expo_t *expo)
 162 {
 163   mpfr_t t, u, k, absx;
 164   int neg_result = 0;
 165   int k_non_zero = 0;
 166   int check_exact_case = 0;
 167   int inexact;
 168   /* Declaration of the size variable */
 169   mpfr_prec_t Nz = MPFR_PREC(z);               /* target precision */
 170   mpfr_prec_t Nt;                              /* working precision */
 171   mpfr_exp_t err;                              /* error */
 172   MPFR_ZIV_DECL (ziv_loop);
 173
 174
 175   MPFR_LOG_FUNC
 176     (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d",
 177       mpfr_get_prec (x), mpfr_log_prec, x,
 178       mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode),
 179      ("z[%Pu]=%.*Rg inexact=%d",
 180       mpfr_get_prec (z), mpfr_log_prec, z, inexact));
 181
 182   /* We put the absolute value of x in absx, pointing to the significand
 183      of x to avoid allocating memory for the significand of absx. */
 184   MPFR_ALIAS(absx, x, /*sign=*/ 1, /*EXP=*/ MPFR_EXP(x));
 185
 186   /* We will compute the absolute value of the result. So, let's
 187      invert the rounding mode if the result is negative. */
 188   if (MPFR_IS_NEG (x) && is_odd (y))
 189     {
 190       neg_result = 1;
 191       rnd_mode = MPFR_INVERT_RND (rnd_mode);
 192     }
 193
 194   /* compute the precision of intermediary variable */
 195   /* the optimal number of bits : see algorithms.tex */
 196   Nt = Nz + 5 + MPFR_INT_CEIL_LOG2 (Nz);
 197
 198   /* initialise of intermediary variable */
 199   mpfr_init2 (t, Nt);
 200
 201   MPFR_ZIV_INIT (ziv_loop, Nt);
 202   for (;;)
 203     {
 204       MPFR_BLOCK_DECL (flags1);
 205
 206       /* compute exp(y*ln|x|), using MPFR_RNDU to get an upper bound, so
 207          that we can detect underflows. */
 208       mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDD : MPFR_RNDU); /* ln|x| */
 209       mpfr_mul (t, y, t, MPFR_RNDU);                              /* y*ln|x| */
 210       if (k_non_zero)
 211         {
 212           MPFR_LOG_MSG (("subtract k * ln(2)\n", 0));
 213           mpfr_const_log2 (u, MPFR_RNDD);
 214           mpfr_mul (u, u, k, MPFR_RNDD);
 215           /* Error on u = k * log(2): < k * 2^(-Nt) < 1. */
 216           mpfr_sub (t, t, u, MPFR_RNDU);
 217           MPFR_LOG_MSG (("t = y * ln|x| - k * ln(2)\n", 0));
 218           MPFR_LOG_VAR (t);
 219         }
 220       /* estimate of the error -- see pow function in algorithms.tex.
 221          The error on t is at most 1/2 + 3*2^(EXP(t)+1) ulps, which is
 222          <= 2^(EXP(t)+3) for EXP(t) >= -1, and <= 2 ulps for EXP(t) <= -2.
 223          Additional error if k_no_zero: treal = t * errk, with
 224          1 - |k| * 2^(-Nt) <= exp(-|k| * 2^(-Nt)) <= errk <= 1,
 225          i.e., additional absolute error <= 2^(EXP(k)+EXP(t)-Nt).
 226          Total error <= 2^err1 + 2^err2 <= 2^(max(err1,err2)+1). */
 227       err = MPFR_NOTZERO (t) && MPFR_GET_EXP (t) >= -1 ?
 228         MPFR_GET_EXP (t) + 3 : 1;
 229       if (k_non_zero)
 230         {
 231           if (MPFR_GET_EXP (k) > err)
 232             err = MPFR_GET_EXP (k);
 233           err++;
 234         }
 235       MPFR_BLOCK (flags1, mpfr_exp (t, t, MPFR_RNDN));  /* exp(y*ln|x|)*/
 236       /* We need to test */
 237       if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (t) || MPFR_UNDERFLOW (flags1)))
 238         {
 239           mpfr_prec_t Ntmin;
 240           MPFR_BLOCK_DECL (flags2);
 241
 242           MPFR_ASSERTN (!k_non_zero);
 243           MPFR_ASSERTN (!MPFR_IS_NAN (t));
 244
 245           /* Real underflow? */
 246           if (MPFR_IS_ZERO (t))
 247             {
 248               /* Underflow. We computed rndn(exp(t)), where t >= y*ln|x|.
 249                  Therefore rndn(|x|^y) = 0, and we have a real underflow on
 250                  |x|^y. */
 251               inexact = mpfr_underflow (z, rnd_mode == MPFR_RNDN ? MPFR_RNDZ
 252                                         : rnd_mode, MPFR_SIGN_POS);
 253               if (expo != NULL)
 254                 MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
 255                                              | MPFR_FLAGS_UNDERFLOW);
 256               break;
 257             }
 258
 259           /* Real overflow? */
 260           if (MPFR_IS_INF (t))
 261             {
 262               /* Note: we can probably use a low precision for this test. */
 263               mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDU : MPFR_RNDD);
 264               mpfr_mul (t, y, t, MPFR_RNDD);            /* y * ln|x| */
 265               MPFR_BLOCK (flags2, mpfr_exp (t, t, MPFR_RNDD));
 266               /* t = lower bound on exp(y * ln|x|) */
 267               if (MPFR_OVERFLOW (flags2))
 268                 {
 269                   /* We have computed a lower bound on |x|^y, and it
 270                      overflowed. Therefore we have a real overflow
 271                      on |x|^y. */
 272                   inexact = mpfr_overflow (z, rnd_mode, MPFR_SIGN_POS);
 273                   if (expo != NULL)
 274                     MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
 275                                                  | MPFR_FLAGS_OVERFLOW);
 276                   break;
 277                 }
 278             }
 279
 280           k_non_zero = 1;
 281           Ntmin = sizeof(mpfr_exp_t) * CHAR_BIT;
 282           if (Ntmin > Nt)
 283             {
 284               Nt = Ntmin;
 285               mpfr_set_prec (t, Nt);
 286             }
 287           mpfr_init2 (u, Nt);
 288           mpfr_init2 (k, Ntmin);
 289           mpfr_log2 (k, absx, MPFR_RNDN);
 290           mpfr_mul (k, y, k, MPFR_RNDN);
 291           mpfr_round (k, k);
 292           MPFR_LOG_VAR (k);
 293           /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */
 294           continue;
 295         }
 296       if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Nz, rnd_mode)))
 297         {
 298           inexact = mpfr_set (z, t, rnd_mode);
 299           break;
 300         }
 301
 302       /* check exact power, except when y is an integer (since the
 303          exact cases for y integer have already been filtered out) */
 304       if (check_exact_case == 0 && ! y_is_integer)
 305         {
 306           if (mpfr_pow_is_exact (z, absx, y, rnd_mode, &inexact))
 307             break;
 308           check_exact_case = 1;
 309         }
 310
 311       /* reactualisation of the precision */
 312       MPFR_ZIV_NEXT (ziv_loop, Nt);
 313       mpfr_set_prec (t, Nt);
 314       if (k_non_zero)
 315         mpfr_set_prec (u, Nt);
 316     }
 317   MPFR_ZIV_FREE (ziv_loop);
 318
 319   if (k_non_zero)
 320     {
 321       int inex2;
 322       long lk;
 323
 324       /* The rounded result in an unbounded exponent range is z * 2^k. As
 325        * MPFR chooses underflow after rounding, the mpfr_mul_2si below will
 326        * correctly detect underflows and overflows. However, in rounding to
 327        * nearest, if z * 2^k = 2^(emin - 2), then the double rounding may
 328        * affect the result. We need to cope with that before overwriting z.
 329        * This can occur only if k < 0 (this test is necessary to avoid a
 330        * potential integer overflow).
 331        * If inexact >= 0, then the real result is <= 2^(emin - 2), so that
 332        * o(2^(emin - 2)) = +0 is correct. If inexact < 0, then the real
 333        * result is > 2^(emin - 2) and we need to round to 2^(emin - 1).
 334        */
 335       MPFR_ASSERTN (MPFR_EXP_MAX <= LONG_MAX);
 336       lk = mpfr_get_si (k, MPFR_RNDN);
 337       /* Due to early overflow detection, |k| should not be much larger than
 338        * MPFR_EMAX_MAX, and as MPFR_EMAX_MAX <= MPFR_EXP_MAX/2 <= LONG_MAX/2,
 339        * an overflow should not be possible in mpfr_get_si (and lk is exact).
 340        * And one even has the following assertion. TODO: complete proof.
 341        */
 342       MPFR_ASSERTD (lk > LONG_MIN && lk < LONG_MAX);
 343       /* Note: even in case of overflow (lk inexact), the code is correct.
 344        * Indeed, for the 3 occurrences of lk:
 345        *   - The test lk < 0 is correct as sign(lk) = sign(k).
 346        *   - In the test MPFR_GET_EXP (z) == __gmpfr_emin - 1 - lk,
 347        *     if lk is inexact, then lk = LONG_MIN <= MPFR_EXP_MIN
 348        *     (the minimum value of the mpfr_exp_t type), and
 349        *     __gmpfr_emin - 1 - lk >= MPFR_EMIN_MIN - 1 - 2 * MPFR_EMIN_MIN
 350        *     >= - MPFR_EMIN_MIN - 1 = MPFR_EMAX_MAX - 1. However, from the
 351        *     choice of k, z has been chosen to be around 1, so that the
 352        *     result of the test is false, as if lk were exact.
 353        *   - In the mpfr_mul_2si (z, z, lk, rnd_mode), if lk is inexact,
 354        *     then |lk| >= LONG_MAX >= MPFR_EXP_MAX, and as z is around 1,
 355        *     mpfr_mul_2si underflows or overflows in the same way as if
 356        *     lk were exact.
 357        * TODO: give a bound on |t|, then on |EXP(z)|.
 358        */
 359       if (rnd_mode == MPFR_RNDN && inexact < 0 && lk < 0 &&
 360           MPFR_GET_EXP (z) == __gmpfr_emin - 1 - lk && mpfr_powerof2_raw (z))
 361         {
 362           /* Rounding to nearest, real result > z * 2^k = 2^(emin - 2),
 363            * underflow case: as the minimum precision is > 1, we will
 364            * obtain the correct result and exceptions by replacing z by
 365            * nextabove(z).
 366            */
 367           MPFR_ASSERTN (MPFR_PREC_MIN > 1);
 368           mpfr_nextabove (z);
 369         }
 370       mpfr_clear_flags ();
 371       inex2 = mpfr_mul_2si (z, z, lk, rnd_mode);
 372       if (inex2)  /* underflow or overflow */
 373         {
 374           inexact = inex2;
 375           if (expo != NULL)
 376             MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, __gmpfr_flags);
 377         }
 378       mpfr_clears (u, k, (mpfr_ptr) 0);
 379     }
 380   mpfr_clear (t);
 381
 382   /* update the sign of the result if x was negative */
 383   if (neg_result)
 384     {
 385       MPFR_SET_NEG(z);
 386       inexact = -inexact;
 387     }
 388
 389   return inexact;
 390 }
 391
 392 /* The computation of z = pow(x,y) is done by
 393    z = exp(y * log(x)) = x^y
 394    For the special cases, see Section F.9.4.4 of the C standard:
 395      _ pow(±0, y) = ±inf for y an odd integer < 0.
 396      _ pow(±0, y) = +inf for y < 0 and not an odd integer.
 397      _ pow(±0, y) = ±0 for y an odd integer > 0.
 398      _ pow(±0, y) = +0 for y > 0 and not an odd integer.
 399      _ pow(-1, ±inf) = 1.
 400      _ pow(+1, y) = 1 for any y, even a NaN.
 401      _ pow(x, ±0) = 1 for any x, even a NaN.
 402      _ pow(x, y) = NaN for finite x < 0 and finite non-integer y.
 403      _ pow(x, -inf) = +inf for |x| < 1.
 404      _ pow(x, -inf) = +0 for |x| > 1.
 405      _ pow(x, +inf) = +0 for |x| < 1.
 406      _ pow(x, +inf) = +inf for |x| > 1.
 407      _ pow(-inf, y) = -0 for y an odd integer < 0.
 408      _ pow(-inf, y) = +0 for y < 0 and not an odd integer.
 409      _ pow(-inf, y) = -inf for y an odd integer > 0.
 410      _ pow(-inf, y) = +inf for y > 0 and not an odd integer.
 411      _ pow(+inf, y) = +0 for y < 0.
 412      _ pow(+inf, y) = +inf for y > 0. */
 413 int
 414 mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
 415 {
 416   int inexact;
 417   int cmp_x_1;
 418   int y_is_integer;
 419   MPFR_SAVE_EXPO_DECL (expo);
 420
 421   MPFR_LOG_FUNC
 422     (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d",
 423       mpfr_get_prec (x), mpfr_log_prec, x,
 424       mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode),
 425      ("z[%Pu]=%.*Rg inexact=%d",
 426       mpfr_get_prec (z), mpfr_log_prec, z, inexact));
 427
 428   if (MPFR_ARE_SINGULAR (x, y))
 429     {
 430       /* pow(x, 0) returns 1 for any x, even a NaN. */
 431       if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
 432         return mpfr_set_ui (z, 1, rnd_mode);
 433       else if (MPFR_IS_NAN (x))
 434         {
 435           MPFR_SET_NAN (z);
 436           MPFR_RET_NAN;
 437         }
 438       else if (MPFR_IS_NAN (y))
 439         {
 440           /* pow(+1, NaN) returns 1. */
 441           if (mpfr_cmp_ui (x, 1) == 0)
 442             return mpfr_set_ui (z, 1, rnd_mode);
 443           MPFR_SET_NAN (z);
 444           MPFR_RET_NAN;
 445         }
 446       else if (MPFR_IS_INF (y))
 447         {
 448           if (MPFR_IS_INF (x))
 449             {
 450               if (MPFR_IS_POS (y))
 451                 MPFR_SET_INF (z);
 452               else
 453                 MPFR_SET_ZERO (z);
 454               MPFR_SET_POS (z);
 455               MPFR_RET (0);
 456             }
 457           else
 458             {
 459               int cmp;
 460               cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y);
 461               MPFR_SET_POS (z);
 462               if (cmp > 0)
 463                 {
 464                   /* Return +inf. */
 465                   MPFR_SET_INF (z);
 466                   MPFR_RET (0);
 467                 }
 468               else if (cmp < 0)
 469                 {
 470                   /* Return +0. */
 471                   MPFR_SET_ZERO (z);
 472                   MPFR_RET (0);
 473                 }
 474               else
 475                 {
 476                   /* Return 1. */
 477                   return mpfr_set_ui (z, 1, rnd_mode);
 478                 }
 479             }
 480         }
 481       else if (MPFR_IS_INF (x))
 482         {
 483           int negative;
 484           /* Determine the sign now, in case y and z are the same object */
 485           negative = MPFR_IS_NEG (x) && is_odd (y);
 486           if (MPFR_IS_POS (y))
 487             MPFR_SET_INF (z);
 488           else
 489             MPFR_SET_ZERO (z);
 490           if (negative)
 491             MPFR_SET_NEG (z);
 492           else
 493             MPFR_SET_POS (z);
 494           MPFR_RET (0);
 495         }
 496       else
 497         {
 498           int negative;
 499           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 500           /* Determine the sign now, in case y and z are the same object */
 501           negative = MPFR_IS_NEG(x) && is_odd (y);
 502           if (MPFR_IS_NEG (y))
 503             {
 504               MPFR_ASSERTD (! MPFR_IS_INF (y));
 505               MPFR_SET_INF (z);
 506               mpfr_set_divby0 ();
 507             }
 508           else
 509             MPFR_SET_ZERO (z);
 510           if (negative)
 511             MPFR_SET_NEG (z);
 512           else
 513             MPFR_SET_POS (z);
 514           MPFR_RET (0);
 515         }
 516     }
 517
 518   /* x^y for x < 0 and y not an integer is not defined */
 519   y_is_integer = mpfr_integer_p (y);
 520   if (MPFR_IS_NEG (x) && ! y_is_integer)
 521     {
 522       MPFR_SET_NAN (z);
 523       MPFR_RET_NAN;
 524     }
 525
 526   /* now the result cannot be NaN:
 527      (1) either x > 0
 528      (2) or x < 0 and y is an integer */
 529
 530   cmp_x_1 = mpfr_cmpabs (x, __gmpfr_one);
 531   if (cmp_x_1 == 0)
 532     return mpfr_set_si (z, MPFR_IS_NEG (x) && is_odd (y) ? -1 : 1, rnd_mode);
 533
 534   /* now we have:
 535      (1) either x > 0
 536      (2) or x < 0 and y is an integer
 537      and in addition |x| <> 1.
 538   */
 539
 540   /* detect overflow: an overflow is possible if
 541      (a) |x| > 1 and y > 0
 542      (b) |x| < 1 and y < 0.
 543      FIXME: this assumes 1 is always representable.
 544
 545      FIXME2: maybe we can test overflow and underflow simultaneously.
 546      The idea is the following: first compute an approximation to
 547      y * log2|x|, using rounding to nearest. If |x| is not too near from 1,
 548      this approximation should be accurate enough, and in most cases enable
 549      one to prove that there is no underflow nor overflow.
 550      Otherwise, it should enable one to check only underflow or overflow,
 551      instead of both cases as in the present case.
 552   */
 553   if (cmp_x_1 * MPFR_SIGN (y) > 0)
 554     {
 555       mpfr_t t;
 556       int negative, overflow;
 557
 558       MPFR_SAVE_EXPO_MARK (expo);
 559       mpfr_init2 (t, 53);
 560       /* we want a lower bound on y*log2|x|:
 561          (i) if x > 0, it suffices to round log2(x) toward zero, and
 562              to round y*o(log2(x)) toward zero too;
 563          (ii) if x < 0, we first compute t = o(-x), with rounding toward 1,
 564               and then follow as in case (1). */
 565       if (MPFR_SIGN (x) > 0)
 566         mpfr_log2 (t, x, MPFR_RNDZ);
 567       else
 568         {
 569           mpfr_neg (t, x, (cmp_x_1 > 0) ? MPFR_RNDZ : MPFR_RNDU);
 570           mpfr_log2 (t, t, MPFR_RNDZ);
 571         }
 572       mpfr_mul (t, t, y, MPFR_RNDZ);
 573       overflow = mpfr_cmp_si (t, __gmpfr_emax) > 0;
 574       mpfr_clear (t);
 575       MPFR_SAVE_EXPO_FREE (expo);
 576       if (overflow)
 577         {
 578           MPFR_LOG_MSG (("early overflow detection\n", 0));
 579           negative = MPFR_SIGN(x) < 0 && is_odd (y);
 580           return mpfr_overflow (z, rnd_mode, negative ? -1 : 1);
 581         }
 582     }
 583
 584   /* Basic underflow checking. One has:
 585    *   - if y > 0, |x^y| < 2^(EXP(x) * y);
 586    *   - if y < 0, |x^y| <= 2^((EXP(x) - 1) * y);
 587    * so that one can compute a value ebound such that |x^y| < 2^ebound.
 588    * If we have ebound <= emin - 2 (emin - 1 in directed rounding modes),
 589    * then there is an underflow and we can decide the return value.
 590    */
 591   if (MPFR_IS_NEG (y) ? (MPFR_GET_EXP (x) > 1) : (MPFR_GET_EXP (x) < 0))
 592     {
 593       mpfr_t tmp;
 594       mpfr_eexp_t ebound;
 595       int inex2;
 596
 597       /* We must restore the flags. */
 598       MPFR_SAVE_EXPO_MARK (expo);
 599       mpfr_init2 (tmp, sizeof (mpfr_exp_t) * CHAR_BIT);
 600       inex2 = mpfr_set_exp_t (tmp, MPFR_GET_EXP (x), MPFR_RNDN);
 601       MPFR_ASSERTN (inex2 == 0);
 602       if (MPFR_IS_NEG (y))
 603         {
 604           inex2 = mpfr_sub_ui (tmp, tmp, 1, MPFR_RNDN);
 605           MPFR_ASSERTN (inex2 == 0);
 606         }
 607       mpfr_mul (tmp, tmp, y, MPFR_RNDU);
 608       if (MPFR_IS_NEG (y))
 609         mpfr_nextabove (tmp);
 610       /* tmp doesn't necessarily fit in ebound, but that doesn't matter
 611          since we get the minimum value in such a case. */
 612       ebound = mpfr_get_exp_t (tmp, MPFR_RNDU);
 613       mpfr_clear (tmp);
 614       MPFR_SAVE_EXPO_FREE (expo);
 615       if (MPFR_UNLIKELY (ebound <=
 616                          __gmpfr_emin - (rnd_mode == MPFR_RNDN ? 2 : 1)))
 617         {
 618           /* warning: mpfr_underflow rounds away from 0 for MPFR_RNDN */
 619           MPFR_LOG_MSG (("early underflow detection\n", 0));
 620           return mpfr_underflow (z,
 621                                  rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
 622                                  MPFR_SIGN (x) < 0 && is_odd (y) ? -1 : 1);
 623         }
 624     }
 625
 626   /* If y is an integer, we can use mpfr_pow_z (based on multiplications),
 627      but if y is very large (I'm not sure about the best threshold -- VL),
 628      we shouldn't use it, as it can be very slow and take a lot of memory
 629      (and even crash or make other programs crash, as several hundred of
 630      MBs may be necessary). Note that in such a case, either x = +/-2^b
 631      (this case is handled below) or x^y cannot be represented exactly in
 632      any precision supported by MPFR (the general case uses this property).
 633   */
 634   if (y_is_integer && (MPFR_GET_EXP (y) <= 256))
 635     {
 636       mpz_t zi;
 637
 638       MPFR_LOG_MSG (("special code for y not too large integer\n", 0));
 639       mpz_init (zi);
 640       mpfr_get_z (zi, y, MPFR_RNDN);
 641       inexact = mpfr_pow_z (z, x, zi, rnd_mode);
 642       mpz_clear (zi);
 643       return inexact;
 644     }
 645
 646   /* Special case (+/-2^b)^Y which could be exact. If x is negative, then
 647      necessarily y is a large integer. */
 648   {
 649     mpfr_exp_t b = MPFR_GET_EXP (x) - 1;
 650
 651     MPFR_ASSERTN (b >= LONG_MIN && b <= LONG_MAX);  /* FIXME... */
 652     if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), b) == 0)
 653       {
 654         mpfr_t tmp;
 655         int sgnx = MPFR_SIGN (x);
 656
 657         MPFR_LOG_MSG (("special case (+/-2^b)^Y\n", 0));
 658         /* now x = +/-2^b, so x^y = (+/-1)^y*2^(b*y) is exact whenever b*y is
 659            an integer */
 660         MPFR_SAVE_EXPO_MARK (expo);
 661         mpfr_init2 (tmp, MPFR_PREC (y) + sizeof (long) * CHAR_BIT);
 662         inexact = mpfr_mul_si (tmp, y, b, MPFR_RNDN); /* exact */
 663         MPFR_ASSERTN (inexact == 0);
 664         /* Note: as the exponent range has been extended, an overflow is not
 665            possible (due to basic overflow and underflow checking above, as
 666            the result is ~ 2^tmp), and an underflow is not possible either
 667            because b is an integer (thus either 0 or >= 1). */
 668         mpfr_clear_flags ();
 669         inexact = mpfr_exp2 (z, tmp, rnd_mode);
 670         mpfr_clear (tmp);
 671         if (sgnx < 0 && is_odd (y))
 672           {
 673             mpfr_neg (z, z, rnd_mode);
 674             inexact = -inexact;
 675           }
 676         /* Without the following, the overflows3 test in tpow.c fails. */
 677         MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 678         MPFR_SAVE_EXPO_FREE (expo);
 679         return mpfr_check_range (z, inexact, rnd_mode);
 680       }
 681   }
 682
 683   MPFR_SAVE_EXPO_MARK (expo);
 684
 685   /* Case where |y * log(x)| is very small. Warning: x can be negative, in
 686      that case y is a large integer. */
 687   {
 688     mpfr_t t;
 689     mpfr_exp_t err;
 690
 691     /* We need an upper bound on the exponent of y * log(x). */
 692     mpfr_init2 (t, 16);
 693     if (MPFR_IS_POS(x))
 694       mpfr_log (t, x, cmp_x_1 < 0 ? MPFR_RNDD : MPFR_RNDU); /* away from 0 */
 695     else
 696       {
 697         /* if x < -1, round to +Inf, else round to zero */
 698         mpfr_neg (t, x, (mpfr_cmp_si (x, -1) < 0) ? MPFR_RNDU : MPFR_RNDD);
 699         mpfr_log (t, t, (mpfr_cmp_ui (t, 1) < 0) ? MPFR_RNDD : MPFR_RNDU);
 700       }
 701     MPFR_ASSERTN (MPFR_IS_PURE_FP (t));
 702     err = MPFR_GET_EXP (y) + MPFR_GET_EXP (t);
 703     mpfr_clear (t);
 704     mpfr_clear_flags ();
 705     MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (z, __gmpfr_one, - err, 0,
 706                                       (MPFR_SIGN (y) > 0) ^ (cmp_x_1 < 0),
 707                                       rnd_mode, expo, {});
 708   }
 709
 710   /* General case */
 711   inexact = mpfr_pow_general (z, x, y, rnd_mode, y_is_integer, &expo);
 712
 713   MPFR_SAVE_EXPO_FREE (expo);
 714   return mpfr_check_range (z, inexact, rnd_mode);
 715 }