1 /* mpfr_expm1 -- Compute exp(x)-1
3 Copyright 2001-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of expm1 is done by
31 mpfr_expm1 (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
35 MPFR_SAVE_EXPO_DECL (expo
);
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
,
42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
49 /* check for inf or -inf (expm1(-inf)=-1) */
50 else if (MPFR_IS_INF (x
))
59 return mpfr_set_si (y
, -1, rnd_mode
);
63 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
64 MPFR_SET_ZERO (y
); /* expm1(+/- 0) = +/- 0 */
65 MPFR_SET_SAME_SIGN (y
, x
);
70 ex
= MPFR_GET_EXP (x
);
73 /* For -1 < x < 0, abs(expm1(x)-x) < x^2/2.
74 For 0 < x < 1, abs(expm1(x)-x) < x^2. */
76 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y
, x
, - ex
, 0, 1, rnd_mode
, {});
78 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y
, x
, - ex
, 1, 0, rnd_mode
, {});
81 MPFR_SAVE_EXPO_MARK (expo
);
83 if (MPFR_IS_NEG (x
) && ex
> 5) /* x <= -32 */
88 mpfr_init2 (minus_one
, 2);
90 mpfr_set_si (minus_one
, -1, MPFR_RNDN
);
91 mpfr_const_log2 (t
, MPFR_RNDU
); /* round upward since x is negative */
92 mpfr_div (t
, x
, t
, MPFR_RNDU
); /* > x / ln(2) */
93 err
= mpfr_cmp_si (t
, MPFR_EMIN_MIN
>= -LONG_MAX
?
94 MPFR_EMIN_MIN
: -LONG_MAX
) <= 0 ?
95 - (MPFR_EMIN_MIN
>= -LONG_MAX
? MPFR_EMIN_MIN
: -LONG_MAX
) :
96 - mpfr_get_si (t
, MPFR_RNDU
);
97 /* exp(x) = 2^(x/ln(2))
98 <= 2^max(MPFR_EMIN_MIN,-LONG_MAX,ceil(x/ln(2)+epsilon))
101 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y
, minus_one
, err
, 0, 0, rnd_mode
,
102 expo
, { mpfr_clear (minus_one
); });
103 mpfr_clear (minus_one
);
108 /* Declaration of the intermediary variable */
110 /* Declaration of the size variable */
111 mpfr_prec_t Ny
= MPFR_PREC(y
); /* target precision */
112 mpfr_prec_t Nt
; /* working precision */
113 mpfr_exp_t err
, exp_te
; /* error */
114 MPFR_ZIV_DECL (loop
);
116 /* compute the precision of intermediary variable */
117 /* the optimal number of bits : see algorithms.tex */
118 Nt
= Ny
+ MPFR_INT_CEIL_LOG2 (Ny
) + 6;
120 /* if |x| is smaller than 2^(-e), we will loose about e bits in the
121 subtraction exp(x) - 1 */
125 /* initialize auxiliary variable */
128 /* First computation of expm1 */
129 MPFR_ZIV_INIT (loop
, Nt
);
132 MPFR_BLOCK_DECL (flags
);
134 /* exp(x) may overflow and underflow */
135 MPFR_BLOCK (flags
, mpfr_exp (t
, x
, MPFR_RNDN
));
136 if (MPFR_OVERFLOW (flags
))
138 inexact
= mpfr_overflow (y
, rnd_mode
, MPFR_SIGN_POS
);
139 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
142 else if (MPFR_UNDERFLOW (flags
))
144 inexact
= mpfr_set_si (y
, -1, rnd_mode
);
145 MPFR_ASSERTD (inexact
== 0);
147 if (MPFR_IS_LIKE_RNDZ (rnd_mode
, 1))
155 exp_te
= MPFR_GET_EXP (t
); /* FIXME: exp(x) may overflow! */
156 mpfr_sub_ui (t
, t
, 1, MPFR_RNDN
); /* exp(x)-1 */
159 /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))));*/
160 err
= Nt
- (MAX (exp_te
- MPFR_GET_EXP (t
), 0) + 1);
162 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, err
, Ny
, rnd_mode
)))
164 inexact
= mpfr_set (y
, t
, rnd_mode
);
168 /* increase the precision */
169 MPFR_ZIV_NEXT (loop
, Nt
);
170 mpfr_set_prec (t
, Nt
);
172 MPFR_ZIV_FREE (loop
);
177 MPFR_SAVE_EXPO_FREE (expo
);
178 return mpfr_check_range (y
, inexact
, rnd_mode
);