contrib/mpfr/src/pow_ui.c

   1 /* mpfr_pow_ui-- compute the power of a floating-point
   2                                   by a machine integer
   3
   4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Caramel projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 3 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  23
  24 #define MPFR_NEED_LONGLONG_H
  25 #include "mpfr-impl.h"
  26
  27 /* sets y to x^n, and return 0 if exact, non-zero otherwise */
  28 int
  29 mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mpfr_rnd_t rnd)
  30 {
  31   unsigned long m;
  32   mpfr_t res;
  33   mpfr_prec_t prec, err;
  34   int inexact;
  35   mpfr_rnd_t rnd1;
  36   MPFR_SAVE_EXPO_DECL (expo);
  37   MPFR_ZIV_DECL (loop);
  38   MPFR_BLOCK_DECL (flags);
  39
  40   MPFR_LOG_FUNC
  41     (("x[%Pu]=%.*Rg n=%lu rnd=%d",
  42       mpfr_get_prec (x), mpfr_log_prec, x, n, rnd),
  43      ("y[%Pu]=%.*Rg inexact=%d",
  44       mpfr_get_prec (y), mpfr_log_prec, y, inexact));
  45
  46   /* x^0 = 1 for any x, even a NaN */
  47   if (MPFR_UNLIKELY (n == 0))
  48     return mpfr_set_ui (y, 1, rnd);
  49
  50   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  51     {
  52       if (MPFR_IS_NAN (x))
  53         {
  54           MPFR_SET_NAN (y);
  55           MPFR_RET_NAN;
  56         }
  57       else if (MPFR_IS_INF (x))
  58         {
  59           /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
  60           if (MPFR_IS_NEG (x) && (n & 1) == 1)
  61             MPFR_SET_NEG (y);
  62           else
  63             MPFR_SET_POS (y);
  64           MPFR_SET_INF (y);
  65           MPFR_RET (0);
  66         }
  67       else /* x is zero */
  68         {
  69           MPFR_ASSERTD (MPFR_IS_ZERO (x));
  70           /* 0^n = 0 for any n */
  71           MPFR_SET_ZERO (y);
  72           if (MPFR_IS_POS (x) || (n & 1) == 0)
  73             MPFR_SET_POS (y);
  74           else
  75             MPFR_SET_NEG (y);
  76           MPFR_RET (0);
  77         }
  78     }
  79   else if (MPFR_UNLIKELY (n <= 2))
  80     {
  81       if (n < 2)
  82         /* x^1 = x */
  83         return mpfr_set (y, x, rnd);
  84       else
  85         /* x^2 = sqr(x) */
  86         return mpfr_sqr (y, x, rnd);
  87     }
  88
  89   /* Augment exponent range */
  90   MPFR_SAVE_EXPO_MARK (expo);
  91
  92   /* setup initial precision */
  93   prec = MPFR_PREC (y) + 3 + GMP_NUMB_BITS
  94     + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
  95   mpfr_init2 (res, prec);
  96
  97   rnd1 = MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD; /* away */
  98
  99   MPFR_ZIV_INIT (loop, prec);
 100   for (;;)
 101     {
 102       int i;
 103
 104       for (m = n, i = 0; m; i++, m >>= 1)
 105         ;
 106       /* now 2^(i-1) <= n < 2^i */
 107       MPFR_ASSERTD (prec > (mpfr_prec_t) i);
 108       err = prec - 1 - (mpfr_prec_t) i;
 109       /* First step: compute square from x */
 110       MPFR_BLOCK (flags,
 111                   inexact = mpfr_mul (res, x, x, MPFR_RNDU);
 112                   MPFR_ASSERTD (i >= 2);
 113                   if (n & (1UL << (i-2)))
 114                     inexact |= mpfr_mul (res, res, x, rnd1);
 115                   for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
 116                     {
 117                       inexact |= mpfr_mul (res, res, res, MPFR_RNDU);
 118                       if (n & (1UL << i))
 119                         inexact |= mpfr_mul (res, res, x, rnd1);
 120                     });
 121       /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2,
 122          and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1.
 123          Using Higham's method, to each rounding corresponds a factor
 124          (1-theta) with 0 <= theta <= 2^(1-p), thus at the end the
 125          absolute error is bounded by (n-1)*2^(1-p)*res <= 2*(n-1)*ulp(res)
 126          since 2^(-p)*x <= ulp(x). Since n < 2^i, this gives a maximal
 127          error of 2^(1+i)*ulp(res).
 128       */
 129       if (MPFR_LIKELY (inexact == 0
 130                        || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
 131                        || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
 132         break;
 133       /* Actualisation of the precision */
 134       MPFR_ZIV_NEXT (loop, prec);
 135       mpfr_set_prec (res, prec);
 136     }
 137   MPFR_ZIV_FREE (loop);
 138
 139   if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)))
 140     {
 141       mpz_t z;
 142
 143       /* Internal overflow or underflow. However the approximation error has
 144        * not been taken into account. So, let's solve this problem by using
 145        * mpfr_pow_z, which can handle it. This case could be improved in the
 146        * future, without having to use mpfr_pow_z.
 147        */
 148       MPFR_LOG_MSG (("Internal overflow or underflow,"
 149                      " let's use mpfr_pow_z.\n", 0));
 150       mpfr_clear (res);
 151       MPFR_SAVE_EXPO_FREE (expo);
 152       mpz_init (z);
 153       mpz_set_ui (z, n);
 154       inexact = mpfr_pow_z (y, x, z, rnd);
 155       mpz_clear (z);
 156       return inexact;
 157     }
 158
 159   inexact = mpfr_set (y, res, rnd);
 160   mpfr_clear (res);
 161
 162   MPFR_SAVE_EXPO_FREE (expo);
 163   return mpfr_check_range (y, inexact, rnd);
 164 }