contrib/mpfr/expm1.c

   1 /* mpfr_expm1 -- Compute exp(x)-1
   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 expm1 is done by
  27     expm1(x)=exp(x)-1
  28  */
  29
  30 int
  31 mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode)
  32 {
  33   int inexact;
  34   mp_exp_t ex;
  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       /* check for inf or -inf (expm1(-inf)=-1) */
  45       else if (MPFR_IS_INF (x))
  46         {
  47           if (MPFR_IS_POS (x))
  48             {
  49               MPFR_SET_INF (y);
  50               MPFR_SET_POS (y);
  51               MPFR_RET (0);
  52             }
  53           else
  54             return mpfr_set_si (y, -1, rnd_mode);
  55         }
  56       else
  57         {
  58           MPFR_ASSERTD (MPFR_IS_ZERO (x));
  59           MPFR_SET_ZERO (y);   /* expm1(+/- 0) = +/- 0 */
  60           MPFR_SET_SAME_SIGN (y, x);
  61           MPFR_RET (0);
  62         }
  63     }
  64
  65   ex = MPFR_GET_EXP (x);
  66   if (ex < 0)
  67     {
  68       /* For -1 < x < 0, abs(expm1(x)-x) < x^2/2.
  69          For 0 < x < 1,  abs(expm1(x)-x) < x^2. */
  70       if (MPFR_IS_POS (x))
  71         MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {});
  72       else
  73         MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, 0, rnd_mode, {});
  74     }
  75
  76   MPFR_SAVE_EXPO_MARK (expo);
  77
  78   if (MPFR_IS_NEG (x) && ex > 5)  /* x <= -32 */
  79     {
  80       mpfr_t minus_one, t;
  81       mp_exp_t err;
  82
  83       mpfr_init2 (minus_one, 2);
  84       mpfr_init2 (t, 64);
  85       mpfr_set_si (minus_one, -1, GMP_RNDN);
  86       mpfr_const_log2 (t, GMP_RNDU); /* round upward since x is negative */
  87       mpfr_div (t, x, t, GMP_RNDU); /* > x / ln(2) */
  88       err = mpfr_cmp_si (t, MPFR_EMIN_MIN >= -LONG_MAX ?
  89                          MPFR_EMIN_MIN : -LONG_MAX) <= 0 ?
  90         - (MPFR_EMIN_MIN >= -LONG_MAX ? MPFR_EMIN_MIN : -LONG_MAX) :
  91         - mpfr_get_si (t, GMP_RNDU);
  92       /* exp(x) = 2^(x/ln(2))
  93                <= 2^max(MPFR_EMIN_MIN,-LONG_MAX,ceil(x/ln(2)+epsilon))
  94          with epsilon > 0 */
  95       mpfr_clear (t);
  96       MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, minus_one, err, 0, 0, rnd_mode,
  97                                         expo, { mpfr_clear (minus_one); });
  98       mpfr_clear (minus_one);
  99     }
 100
 101   /* General case */
 102   {
 103     /* Declaration of the intermediary variable */
 104     mpfr_t t;
 105     /* Declaration of the size variable */
 106     mp_prec_t Ny = MPFR_PREC(y);   /* target precision */
 107     mp_prec_t Nt;                  /* working precision */
 108     mp_exp_t err, exp_te;          /* error */
 109     MPFR_ZIV_DECL (loop);
 110
 111     /* compute the precision of intermediary variable */
 112     /* the optimal number of bits : see algorithms.tex */
 113     Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6;
 114
 115     /* if |x| is smaller than 2^(-e), we will loose about e bits in the
 116        subtraction exp(x) - 1 */
 117     if (ex < 0)
 118       Nt += - ex;
 119
 120     /* initialize auxiliary variable */
 121     mpfr_init2 (t, Nt);
 122
 123     /* First computation of expm1 */
 124     MPFR_ZIV_INIT (loop, Nt);
 125     for (;;)
 126       {
 127         MPFR_BLOCK_DECL (flags);
 128
 129         /* exp(x) may overflow and underflow */
 130         MPFR_BLOCK (flags, mpfr_exp (t, x, GMP_RNDN));
 131         if (MPFR_OVERFLOW (flags))
 132           {
 133             inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN_POS);
 134             MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 135             break;
 136           }
 137         else if (MPFR_UNDERFLOW (flags))
 138           {
 139             inexact = mpfr_set_si (y, -1, rnd_mode);
 140             MPFR_ASSERTD (inexact == 0);
 141             inexact = -1;
 142             if (MPFR_IS_LIKE_RNDZ (rnd_mode, 1))
 143               {
 144                 inexact = 1;
 145                 mpfr_nexttozero (y);
 146               }
 147             break;
 148           }
 149
 150         exp_te = MPFR_GET_EXP (t);         /* FIXME: exp(x) may overflow! */
 151         mpfr_sub_ui (t, t, 1, GMP_RNDN);   /* exp(x)-1 */
 152
 153         /* error estimate */
 154         /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))));*/
 155         err = Nt - (MAX (exp_te - MPFR_GET_EXP (t), 0) + 1);
 156
 157         if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
 158           {
 159             inexact = mpfr_set (y, t, rnd_mode);
 160             break;
 161           }
 162
 163         /* increase the precision */
 164         MPFR_ZIV_NEXT (loop, Nt);
 165         mpfr_set_prec (t, Nt);
 166       }
 167     MPFR_ZIV_FREE (loop);
 168
 169     mpfr_clear (t);
 170   }
 171
 172   MPFR_SAVE_EXPO_FREE (expo);
 173   return mpfr_check_range (y, inexact, rnd_mode);
 174 }