MPFR v4.0.2
[mpfr.git] / src / expm1.c
1 /* mpfr_expm1 -- Compute exp(x)-1
2
3 Copyright 2001-2019 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 expm1 is done by
27     expm1(x)=exp(x)-1
28  */
29
30 int
31 mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mpfr_rnd_t rnd_mode)
32 {
33   int inexact;
34   mpfr_exp_t ex;
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       /* check for inf or -inf (expm1(-inf)=-1) */
50       else if (MPFR_IS_INF (x))
51         {
52           if (MPFR_IS_POS (x))
53             {
54               MPFR_SET_INF (y);
55               MPFR_SET_POS (y);
56               MPFR_RET (0);
57             }
58           else
59             return mpfr_set_si (y, -1, rnd_mode);
60         }
61       else
62         {
63           MPFR_ASSERTD (MPFR_IS_ZERO (x));
64           MPFR_SET_ZERO (y);   /* expm1(+/- 0) = +/- 0 */
65           MPFR_SET_SAME_SIGN (y, x);
66           MPFR_RET (0);
67         }
68     }
69
70   ex = MPFR_GET_EXP (x);
71   if (ex < 0)
72     {
73       /* For -1 < x < 0, abs(expm1(x)-x) < x^2/2.
74          For 0 < x < 1,  abs(expm1(x)-x) < x^2. */
75       if (MPFR_IS_POS (x))
76         MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {});
77       else
78         MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, 0, rnd_mode, {});
79     }
80
81   MPFR_SAVE_EXPO_MARK (expo);
82
83   if (MPFR_IS_NEG (x) && ex > 5)  /* x <= -32 */
84     {
85       mp_limb_t t_limb[(64 - 1) / GMP_NUMB_BITS + 1];
86       mpfr_t t;
87       mpfr_exp_t err;
88
89       MPFR_TMP_INIT1(t_limb, t, 64);
90       mpfr_div (t, x, __gmpfr_const_log2_RNDU, MPFR_RNDU); /* > x / ln(2) */
91       err = mpfr_cmp_si (t, MPFR_EMIN_MIN >= -LONG_MAX ?
92                          MPFR_EMIN_MIN : -LONG_MAX) <= 0 ?
93         - (MPFR_EMIN_MIN >= -LONG_MAX ? MPFR_EMIN_MIN : -LONG_MAX) :
94         - mpfr_get_si (t, MPFR_RNDU);
95       /* exp(x) = 2^(x/ln(2))
96                <= 2^max(MPFR_EMIN_MIN,-LONG_MAX,ceil(x/ln(2)+epsilon))
97          with epsilon > 0 */
98       MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_mone, err, 0, 0,
99                                         rnd_mode, expo, {});
100     }
101
102   /* General case */
103   {
104     /* Declaration of the intermediary variable */
105     mpfr_t t;
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, exp_te;          /* 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 + MPFR_INT_CEIL_LOG2 (Ny) + 6;
115
116     /* if |x| is smaller than 2^(-e), we will loose about e bits in the
117        subtraction exp(x) - 1 */
118     if (ex < 0)
119       Nt += - ex;
120
121     /* initialize auxiliary variable */
122     mpfr_init2 (t, Nt);
123
124     /* First computation of expm1 */
125     MPFR_ZIV_INIT (loop, Nt);
126     for (;;)
127       {
128         MPFR_BLOCK_DECL (flags);
129
130         /* exp(x) may overflow and underflow */
131         MPFR_BLOCK (flags, mpfr_exp (t, x, MPFR_RNDN));
132         if (MPFR_OVERFLOW (flags))
133           {
134             inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN_POS);
135             MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
136             break;
137           }
138         else if (MPFR_UNDERFLOW (flags))
139           {
140             inexact = mpfr_set_si (y, -1, rnd_mode);
141             MPFR_ASSERTD (inexact == 0);
142             inexact = -1;
143             if (MPFR_IS_LIKE_RNDZ (rnd_mode, 1))
144               {
145                 inexact = 1;
146                 mpfr_nexttozero (y);
147               }
148             break;
149           }
150
151         exp_te = MPFR_GET_EXP (t);         /* FIXME: exp(x) may overflow! */
152         mpfr_sub_ui (t, t, 1, MPFR_RNDN);   /* exp(x)-1 */
153
154         /* error estimate */
155         /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))));*/
156         err = Nt - (MAX (exp_te - MPFR_GET_EXP (t), 0) + 1);
157
158         if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
159           {
160             inexact = mpfr_set (y, t, rnd_mode);
161             break;
162           }
163
164         /* increase the precision */
165         MPFR_ZIV_NEXT (loop, Nt);
166         mpfr_set_prec (t, Nt);
167       }
168     MPFR_ZIV_FREE (loop);
169
170     mpfr_clear (t);
171   }
172
173   MPFR_SAVE_EXPO_FREE (expo);
174   return mpfr_check_range (y, inexact, rnd_mode);
175 }