contrib/mpfr/pow_si.c

   1 /* mpfr_pow_si -- power function x^y with y a signed int
   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 /* The computation of y = pow_si(x,n) is done by
  27  *    y = pow_ui(x,n)       if n >= 0
  28  *    y = 1 / pow_ui(x,-n)  if n < 0
  29  */
  30
  31 int
  32 mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
  33 {
  34   MPFR_LOG_FUNC (("x[%#R]=%R n=%ld rnd=%d", x, x, n, rnd),
  35                  ("y[%#R]=%R", y, y));
  36
  37   if (n >= 0)
  38     return mpfr_pow_ui (y, x, n, rnd);
  39   else
  40     {
  41       if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  42         {
  43           if (MPFR_IS_NAN (x))
  44             {
  45               MPFR_SET_NAN (y);
  46               MPFR_RET_NAN;
  47             }
  48           else if (MPFR_IS_INF (x))
  49             {
  50               MPFR_SET_ZERO (y);
  51               if (MPFR_IS_POS (x) || ((unsigned) n & 1) == 0)
  52                 MPFR_SET_POS (y);
  53               else
  54                 MPFR_SET_NEG (y);
  55               MPFR_RET (0);
  56             }
  57           else /* x is zero */
  58             {
  59               MPFR_ASSERTD (MPFR_IS_ZERO (x));
  60               MPFR_SET_INF(y);
  61               if (MPFR_IS_POS (x) || ((unsigned) n & 1) == 0)
  62                 MPFR_SET_POS (y);
  63               else
  64                 MPFR_SET_NEG (y);
  65               MPFR_RET(0);
  66             }
  67         }
  68       MPFR_CLEAR_FLAGS (y);
  69
  70       /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
  71       if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0)
  72         {
  73           mp_exp_t expx = MPFR_EXP (x) - 1, expy;
  74           MPFR_ASSERTD (n < 0);
  75           /* Warning: n * expx may overflow!
  76            * Some systems (apparently alpha-freebsd) abort with
  77            * LONG_MIN / 1, and LONG_MIN / -1 is undefined.
  78            * Proof of the overflow checking. The expressions below are
  79            * assumed to be on the rational numbers, but the word "overflow"
  80            * still has its own meaning in the C context. / still denotes
  81            * the integer (truncated) division, and // denotes the exact
  82            * division.
  83            * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n
  84            *   cannot overflow due to the constraints on the exponents of
  85            *   MPFR numbers.
  86            * - If n = -1, then n * expx = - expx, which is representable
  87            *   because of the constraints on the exponents of MPFR numbers.
  88            * - If expx = 0, then n * expx = 0, which is representable.
  89            * - If n < -1 and expx > 0:
  90            *   + If expx > (__gmpfr_emin - 1) / n, then
  91            *           expx >= (__gmpfr_emin - 1) / n + 1
  92            *                > (__gmpfr_emin - 1) // n,
  93            *     and
  94            *           n * expx < __gmpfr_emin - 1,
  95            *     i.e.
  96            *           n * expx <= __gmpfr_emin - 2.
  97            *     This corresponds to an underflow, with a null result in
  98            *     the rounding-to-nearest mode.
  99            *   + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot
 100            *     overflow since 0 < expx <= (__gmpfr_emin - 1) / n and
 101            *           0 > n * expx >= n * ((__gmpfr_emin - 1) / n)
 102            *                        >= __gmpfr_emin - 1.
 103            * - If n < -1 and expx < 0:
 104            *   + If expx < (__gmpfr_emax - 1) / n, then
 105            *           expx <= (__gmpfr_emax - 1) / n - 1
 106            *                < (__gmpfr_emax - 1) // n,
 107            *     and
 108            *           n * expx > __gmpfr_emax - 1,
 109            *     i.e.
 110            *           n * expx >= __gmpfr_emax.
 111            *     This corresponds to an overflow (2^(n * expx) has an
 112            *     exponent > __gmpfr_emax).
 113            *   + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot
 114            *     overflow since 0 > expx >= (__gmpfr_emax - 1) / n and
 115            *           0 < n * expx <= n * ((__gmpfr_emax - 1) / n)
 116            *                        <= __gmpfr_emax - 1.
 117            * Note: one could use expx bounds based on MPFR_EXP_MIN and
 118            * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The
 119            * current bounds do not lead to noticeably slower code and
 120            * allow us to avoid a bug in Sun's compiler for Solaris/x86
 121            * (when optimizations are enabled).
 122            */
 123           expy =
 124             n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
 125             MPFR_EMIN_MIN - 2 /* Underflow */ :
 126             n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ?
 127             MPFR_EMAX_MAX /* Overflow */ : n * expx;
 128           return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1,
 129                                    expy, rnd);
 130         }
 131
 132       /* General case */
 133       {
 134         /* Declaration of the intermediary variable */
 135         mpfr_t t;
 136         /* Declaration of the size variable */
 137         mp_prec_t Ny;                              /* target precision */
 138         mp_prec_t Nt;                              /* working precision */
 139         mp_rnd_t rnd1;
 140         int size_n;
 141         int inexact;
 142         unsigned long abs_n;
 143         MPFR_SAVE_EXPO_DECL (expo);
 144         MPFR_ZIV_DECL (loop);
 145
 146         abs_n = - (unsigned long) n;
 147         count_leading_zeros (size_n, (mp_limb_t) abs_n);
 148         size_n = BITS_PER_MP_LIMB - size_n;
 149
 150         /* initial working precision */
 151         Ny = MPFR_PREC (y);
 152         Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny);
 153
 154         MPFR_SAVE_EXPO_MARK (expo);
 155
 156         /* initialise of intermediary   variable */
 157         mpfr_init2 (t, Nt);
 158
 159         /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
 160            toward sign(x), to avoid spurious overflow or underflow, as in
 161            mpfr_pow_z. */
 162         rnd1 = MPFR_EXP (x) < 1 ? GMP_RNDZ :
 163           (MPFR_SIGN (x) > 0 ? GMP_RNDU : GMP_RNDD);
 164
 165         MPFR_ZIV_INIT (loop, Nt);
 166         for (;;)
 167           {
 168             MPFR_BLOCK_DECL (flags);
 169
 170             /* compute (1/x)^|n| */
 171             MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
 172             MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
 173             /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
 174             if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
 175               goto overflow;
 176             MPFR_BLOCK (flags, mpfr_pow_ui (t, t, abs_n, rnd));
 177             /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
 178                If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
 179                Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
 180                from algorithms.tex, which yields x^n*(1+theta) with
 181                |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
 182                2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
 183             if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
 184               {
 185               overflow:
 186                 MPFR_ZIV_FREE (loop);
 187                 mpfr_clear (t);
 188                 MPFR_SAVE_EXPO_FREE (expo);
 189                 MPFR_LOG_MSG (("overflow\n", 0));
 190                 return mpfr_overflow (y, rnd, abs_n & 1 ?
 191                                       MPFR_SIGN (x) : MPFR_SIGN_POS);
 192               }
 193             if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
 194               {
 195                 MPFR_ZIV_FREE (loop);
 196                 mpfr_clear (t);
 197                 MPFR_LOG_MSG (("underflow\n", 0));
 198                 if (rnd == GMP_RNDN)
 199                   {
 200                     mpfr_t y2, nn;
 201
 202                     /* We cannot decide now whether the result should be
 203                        rounded toward zero or away from zero. So, like
 204                        in mpfr_pow_pos_z, let's use the general case of
 205                        mpfr_pow in precision 2. */
 206                     MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
 207                                                     MPFR_EXP (x) - 1) != 0);
 208                     mpfr_init2 (y2, 2);
 209                     mpfr_init2 (nn, sizeof (long) * CHAR_BIT);
 210                     inexact = mpfr_set_si (nn, n, GMP_RNDN);
 211                     MPFR_ASSERTN (inexact == 0);
 212                     inexact = mpfr_pow_general (y2, x, nn, rnd, 1,
 213                                                 (mpfr_save_expo_t *) NULL);
 214                     mpfr_clear (nn);
 215                     mpfr_set (y, y2, GMP_RNDN);
 216                     mpfr_clear (y2);
 217                     MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
 218                     goto end;
 219                   }
 220                 else
 221                   {
 222                     MPFR_SAVE_EXPO_FREE (expo);
 223                     return mpfr_underflow (y, rnd, abs_n & 1 ?
 224                                            MPFR_SIGN (x) : MPFR_SIGN_POS);
 225                   }
 226               }
 227             /* error estimate -- see pow function in algorithms.ps */
 228             if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_n - 2, Ny, rnd)))
 229               break;
 230
 231             /* actualisation of the precision */
 232             MPFR_ZIV_NEXT (loop, Nt);
 233             mpfr_set_prec (t, Nt);
 234           }
 235         MPFR_ZIV_FREE (loop);
 236
 237         inexact = mpfr_set (y, t, rnd);
 238         mpfr_clear (t);
 239
 240       end:
 241         MPFR_SAVE_EXPO_FREE (expo);
 242         return mpfr_check_range (y, inexact, rnd);
 243       }
 244     }
 245 }