exp.c

   1 /* mpfr_exp -- exponential of a floating-point number
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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 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 #include "mpfr-impl.h"
  24
  25 /* #define DEBUG */
  26
  27 int
  28 mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  29 {
  30   mpfr_exp_t expx;
  31   mpfr_prec_t precy;
  32   int inexact;
  33   MPFR_SAVE_EXPO_DECL (expo);
  34
  35   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
  36                  ("y[%#R]=%R inexact=%d", y, y, inexact));
  37
  38   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
  39     {
  40       if (MPFR_IS_NAN(x))
  41         {
  42           MPFR_SET_NAN(y);
  43           MPFR_RET_NAN;
  44         }
  45       else if (MPFR_IS_INF(x))
  46         {
  47           if (MPFR_IS_POS(x))
  48             MPFR_SET_INF(y);
  49           else
  50             MPFR_SET_ZERO(y);
  51           MPFR_SET_POS(y);
  52           MPFR_RET(0);
  53         }
  54       else
  55         {
  56           MPFR_ASSERTD(MPFR_IS_ZERO(x));
  57           return mpfr_set_ui (y, 1, rnd_mode);
  58         }
  59     }
  60
  61   /* First, let's detect most overflow and underflow cases. */
  62   {
  63     mpfr_t e, bound;
  64
  65     /* We must extended the exponent range and save the flags now. */
  66     MPFR_SAVE_EXPO_MARK (expo);
  67
  68     mpfr_init2 (e, sizeof (mpfr_exp_t) * CHAR_BIT);
  69     mpfr_init2 (bound, 32);
  70
  71     inexact = mpfr_set_exp_t (e, expo.saved_emax, MPFR_RNDN);
  72     MPFR_ASSERTD (inexact == 0);
  73     mpfr_const_log2 (bound, expo.saved_emax < 0 ? MPFR_RNDD : MPFR_RNDU);
  74     mpfr_mul (bound, bound, e, MPFR_RNDU);
  75     if (MPFR_UNLIKELY (mpfr_cmp (x, bound) >= 0))
  76       {
  77         /* x > log(2^emax), thus exp(x) > 2^emax */
  78         mpfr_clears (e, bound, (mpfr_ptr) 0);
  79         MPFR_SAVE_EXPO_FREE (expo);
  80         return mpfr_overflow (y, rnd_mode, 1);
  81       }
  82
  83     inexact = mpfr_set_exp_t (e, expo.saved_emin, MPFR_RNDN);
  84     MPFR_ASSERTD (inexact == 0);
  85     inexact = mpfr_sub_ui (e, e, 2, MPFR_RNDN);
  86     MPFR_ASSERTD (inexact == 0);
  87     mpfr_const_log2 (bound, expo.saved_emin < 0 ? MPFR_RNDU : MPFR_RNDD);
  88     mpfr_mul (bound, bound, e, MPFR_RNDD);
  89     if (MPFR_UNLIKELY (mpfr_cmp (x, bound) <= 0))
  90       {
  91         /* x < log(2^(emin - 2)), thus exp(x) < 2^(emin - 2) */
  92         mpfr_clears (e, bound, (mpfr_ptr) 0);
  93         MPFR_SAVE_EXPO_FREE (expo);
  94         return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
  95                                1);
  96       }
  97
  98     /* Other overflow/underflow cases must be detected
  99        by the generic routines. */
 100     mpfr_clears (e, bound, (mpfr_ptr) 0);
 101     MPFR_SAVE_EXPO_FREE (expo);
 102   }
 103
 104   expx  = MPFR_GET_EXP (x);
 105   precy = MPFR_PREC (y);
 106
 107   /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */
 108   if (MPFR_UNLIKELY (expx < 0 && (mpfr_uexp_t) (-expx) > precy))
 109     {
 110       mpfr_exp_t emin = __gmpfr_emin;
 111       mpfr_exp_t emax = __gmpfr_emax;
 112       int signx = MPFR_SIGN (x);
 113
 114       MPFR_SET_POS (y);
 115       if (MPFR_IS_NEG_SIGN (signx) && (rnd_mode == MPFR_RNDD ||
 116                                        rnd_mode == MPFR_RNDZ))
 117         {
 118           __gmpfr_emin = 0;
 119           __gmpfr_emax = 0;
 120           mpfr_setmax (y, 0);  /* y = 1 - epsilon */
 121           inexact = -1;
 122         }
 123       else
 124         {
 125           __gmpfr_emin = 1;
 126           __gmpfr_emax = 1;
 127           mpfr_setmin (y, 1);  /* y = 1 */
 128           if (MPFR_IS_POS_SIGN (signx) && (rnd_mode == MPFR_RNDU ||
 129                                            rnd_mode == MPFR_RNDA))
 130             {
 131               mp_size_t yn;
 132               int sh;
 133
 134               yn = 1 + (MPFR_PREC(y) - 1) / GMP_NUMB_BITS;
 135               sh = (mpfr_prec_t) yn * GMP_NUMB_BITS - MPFR_PREC(y);
 136               MPFR_MANT(y)[0] += MPFR_LIMB_ONE << sh;
 137               inexact = 1;
 138             }
 139           else
 140             inexact = -MPFR_FROM_SIGN_TO_INT(signx);
 141         }
 142
 143       __gmpfr_emin = emin;
 144       __gmpfr_emax = emax;
 145     }
 146   else  /* General case */
 147     {
 148       if (MPFR_UNLIKELY (precy >= MPFR_EXP_THRESHOLD))
 149         /* mpfr_exp_3 saves the exponent range and flags itself, otherwise
 150            the flag changes in mpfr_exp_3 are lost */
 151         inexact = mpfr_exp_3 (y, x, rnd_mode); /* O(M(n) log(n)^2) */
 152       else
 153         {
 154           MPFR_SAVE_EXPO_MARK (expo);
 155           inexact = mpfr_exp_2 (y, x, rnd_mode); /* O(n^(1/3) M(n)) */
 156           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 157           MPFR_SAVE_EXPO_FREE (expo);
 158         }
 159     }
 160
 161   return mpfr_check_range (y, inexact, rnd_mode);
 162 }