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