contrib/mpfr/src/exp2.c

   1 /* mpfr_exp2 -- power of 2 function 2^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 /* The computation of y = 2^z is done by                           *
  27  *     y = exp(z*log(2)). The result is exact iff z is an integer. */
  28
  29 int
  30 mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  31 {
  32   int inexact;
  33   long xint;
  34   mpfr_t xfrac;
  35   MPFR_SAVE_EXPO_DECL (expo);
  36
  37   MPFR_LOG_FUNC
  38     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
  39      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y,
  40       inexact));
  41
  42   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  43     {
  44       if (MPFR_IS_NAN (x))
  45         {
  46           MPFR_SET_NAN (y);
  47           MPFR_RET_NAN;
  48         }
  49       else if (MPFR_IS_INF (x))
  50         {
  51           if (MPFR_IS_POS (x))
  52             MPFR_SET_INF (y);
  53           else
  54             MPFR_SET_ZERO (y);
  55           MPFR_SET_POS (y);
  56           MPFR_RET (0);
  57         }
  58       else /* 2^0 = 1 */
  59         {
  60           MPFR_ASSERTD (MPFR_IS_ZERO(x));
  61           return mpfr_set_ui (y, 1, rnd_mode);
  62         }
  63     }
  64
  65   /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin,
  66      if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */
  67   MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN + 2);
  68   if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 1) < 0))
  69     {
  70       mpfr_rnd_t rnd2 = rnd_mode;
  71       /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */
  72       if (rnd_mode == MPFR_RNDN &&
  73           mpfr_cmp_si_2exp (x, __gmpfr_emin - 2, 0) <= 0)
  74         rnd2 = MPFR_RNDZ;
  75       return mpfr_underflow (y, rnd2, 1);
  76     }
  77
  78   MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
  79   if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax) >= 0))
  80     return mpfr_overflow (y, rnd_mode, 1);
  81
  82   /* We now know that emin - 1 <= x < emax. */
  83
  84   MPFR_SAVE_EXPO_MARK (expo);
  85
  86   /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have
  87      |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1);
  88      if x < 0 we must round toward 0 (dir=0). */
  89   MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, - MPFR_GET_EXP (x), 0,
  90                                     MPFR_SIGN(x) > 0, rnd_mode, expo, {});
  91
  92   xint = mpfr_get_si (x, MPFR_RNDZ);
  93   mpfr_init2 (xfrac, MPFR_PREC (x));
  94   mpfr_sub_si (xfrac, x, xint, MPFR_RNDN); /* exact */
  95
  96   if (MPFR_IS_ZERO (xfrac))
  97     {
  98       mpfr_set_ui (y, 1, MPFR_RNDN);
  99       inexact = 0;
 100     }
 101   else
 102     {
 103       /* Declaration of the intermediary variable */
 104       mpfr_t t;
 105
 106       /* Declaration of the size variable */
 107       mpfr_prec_t Ny = MPFR_PREC(y);              /* target precision */
 108       mpfr_prec_t Nt;                             /* working precision */
 109       mpfr_exp_t err;                             /* error */
 110       MPFR_ZIV_DECL (loop);
 111
 112       /* compute the precision of intermediary variable */
 113       /* the optimal number of bits : see algorithms.tex */
 114       Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny);
 115
 116       /* initialise of intermediary variable */
 117       mpfr_init2 (t, Nt);
 118
 119       /* First computation */
 120       MPFR_ZIV_INIT (loop, Nt);
 121       for (;;)
 122         {
 123           /* compute exp(x*ln(2))*/
 124           mpfr_const_log2 (t, MPFR_RNDU);       /* ln(2) */
 125           mpfr_mul (t, xfrac, t, MPFR_RNDU);    /* xfrac * ln(2) */
 126           err = Nt - (MPFR_GET_EXP (t) + 2);   /* Estimate of the error */
 127           mpfr_exp (t, t, MPFR_RNDN);           /* exp(xfrac * ln(2)) */
 128
 129           if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
 130             break;
 131
 132           /* Actualisation of the precision */
 133           MPFR_ZIV_NEXT (loop, Nt);
 134           mpfr_set_prec (t, Nt);
 135         }
 136       MPFR_ZIV_FREE (loop);
 137
 138       inexact = mpfr_set (y, t, rnd_mode);
 139
 140       mpfr_clear (t);
 141     }
 142
 143   mpfr_clear (xfrac);
 144   mpfr_clear_flags ();
 145   mpfr_mul_2si (y, y, xint, MPFR_RNDN); /* exact or overflow */
 146   /* Note: We can have an overflow only when t was rounded up to 2. */
 147   MPFR_ASSERTD (MPFR_IS_PURE_FP (y) || inexact > 0);
 148   MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 149   MPFR_SAVE_EXPO_FREE (expo);
 150   return mpfr_check_range (y, inexact, rnd_mode);
 151 }