contrib/mpfr/exp2.c

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