contrib/mpfr/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 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 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, mp_rnd_t rnd)
  30 {
  31   unsigned long m;
  32   mpfr_t res;
  33   mp_prec_t prec, err;
  34   int inexact;
  35   mp_rnd_t rnd1;
  36   MPFR_SAVE_EXPO_DECL (expo);
  37   MPFR_ZIV_DECL (loop);
  38   MPFR_BLOCK_DECL (flags);
  39
  40   MPFR_LOG_FUNC (("x[%#R]=%R n=%lu rnd=%d", x, x, n, rnd),
  41                  ("y[%#R]=%R inexact=%d", y, y, inexact));
  42
  43   /* x^0 = 1 for any x, even a NaN */
  44   if (MPFR_UNLIKELY (n == 0))
  45     return mpfr_set_ui (y, 1, rnd);
  46
  47   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  48     {
  49       if (MPFR_IS_NAN (x))
  50         {
  51           MPFR_SET_NAN (y);
  52           MPFR_RET_NAN;
  53         }
  54       else if (MPFR_IS_INF (x))
  55         {
  56           /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
  57           if (MPFR_IS_NEG (x) && (n & 1) == 1)
  58             MPFR_SET_NEG (y);
  59           else
  60             MPFR_SET_POS (y);
  61           MPFR_SET_INF (y);
  62           MPFR_RET (0);
  63         }
  64       else /* x is zero */
  65         {
  66           MPFR_ASSERTD (MPFR_IS_ZERO (x));
  67           /* 0^n = 0 for any n */
  68           MPFR_SET_ZERO (y);
  69           if (MPFR_IS_POS (x) || (n & 1) == 0)
  70             MPFR_SET_POS (y);
  71           else
  72             MPFR_SET_NEG (y);
  73           MPFR_RET (0);
  74         }
  75     }
  76   else if (MPFR_UNLIKELY (n <= 2))
  77     {
  78       if (n < 2)
  79         /* x^1 = x */
  80         return mpfr_set (y, x, rnd);
  81       else
  82         /* x^2 = sqr(x) */
  83         return mpfr_sqr (y, x, rnd);
  84     }
  85
  86   /* Augment exponent range */
  87   MPFR_SAVE_EXPO_MARK (expo);
  88
  89   /* setup initial precision */
  90   prec = MPFR_PREC (y) + 3 + BITS_PER_MP_LIMB
  91     + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
  92   mpfr_init2 (res, prec);
  93
  94   rnd1 = MPFR_IS_POS (x) ? GMP_RNDU : GMP_RNDD; /* away */
  95
  96   MPFR_ZIV_INIT (loop, prec);
  97   for (;;)
  98     {
  99       int i;
 100
 101       for (m = n, i = 0; m; i++, m >>= 1)
 102         ;
 103       /* now 2^(i-1) <= n < 2^i */
 104       MPFR_ASSERTD (prec > (mpfr_prec_t) i);
 105       err = prec - 1 - (mpfr_prec_t) i;
 106       /* First step: compute square from x */
 107       MPFR_BLOCK (flags,
 108                   inexact = mpfr_mul (res, x, x, GMP_RNDU);
 109                   MPFR_ASSERTD (i >= 2);
 110                   if (n & (1UL << (i-2)))
 111                     inexact |= mpfr_mul (res, res, x, rnd1);
 112                   for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
 113                     {
 114                       inexact |= mpfr_mul (res, res, res, GMP_RNDU);
 115                       if (n & (1UL << i))
 116                         inexact |= mpfr_mul (res, res, x, rnd1);
 117                     });
 118       /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2,
 119          and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1.
 120          Using Higham's method, to each rounding corresponds a factor
 121          (1-theta) with 0 <= theta <= 2^(1-p), thus at the end the
 122          absolute error is bounded by (n-1)*2^(1-p)*res <= 2*(n-1)*ulp(res)
 123          since 2^(-p)*x <= ulp(x). Since n < 2^i, this gives a maximal
 124          error of 2^(1+i)*ulp(res).
 125       */
 126       if (MPFR_LIKELY (inexact == 0
 127                        || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
 128                        || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
 129         break;
 130       /* Actualisation of the precision */
 131       MPFR_ZIV_NEXT (loop, prec);
 132       mpfr_set_prec (res, prec);
 133     }
 134   MPFR_ZIV_FREE (loop);
 135
 136   if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)))
 137     {
 138       mpz_t z;
 139
 140       /* Internal overflow or underflow. However the approximation error has
 141        * not been taken into account. So, let's solve this problem by using
 142        * mpfr_pow_z, which can handle it. This case could be improved in the
 143        * future, without having to use mpfr_pow_z.
 144        */
 145       MPFR_LOG_MSG (("Internal overflow or underflow,"
 146                      " let's use mpfr_pow_z.\n", 0));
 147       mpfr_clear (res);
 148       MPFR_SAVE_EXPO_FREE (expo);
 149       mpz_init (z);
 150       mpz_set_ui (z, n);
 151       inexact = mpfr_pow_z (y, x, z, rnd);
 152       mpz_clear (z);
 153       return inexact;
 154     }
 155
 156   inexact = mpfr_set (y, res, rnd);
 157   mpfr_clear (res);
 158
 159   MPFR_SAVE_EXPO_FREE (expo);
 160   return mpfr_check_range (y, inexact, rnd);
 161 }