1// (C) Copyright John Maddock 2006.
2// Use, modification and distribution are subject to the
3// Boost Software License, Version 1.0. (See accompanying file
4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6#ifndef BOOST_MATH_EXPM1_INCLUDED
7#define BOOST_MATH_EXPM1_INCLUDED
8
9#ifdef _MSC_VER
10#pragma once
11#endif
12
13#include <boost/config/no_tr1/cmath.hpp>
14#include <math.h> // platform's ::expm1
15#include <boost/limits.hpp>
16#include <boost/math/tools/config.hpp>
17#include <boost/math/tools/series.hpp>
18#include <boost/math/tools/precision.hpp>
19#include <boost/math/tools/big_constant.hpp>
20#include <boost/math/policies/error_handling.hpp>
21#include <boost/math/tools/rational.hpp>
22#include <boost/math/special_functions/math_fwd.hpp>
23#include <boost/mpl/less_equal.hpp>
24
25#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
26# include <boost/static_assert.hpp>
27#else
28# include <boost/assert.hpp>
29#endif
30
31#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
32//
33// This is the only way we can avoid
34// warning: non-standard suffix on floating constant [-Wpedantic]
35// when building with -Wall -pedantic. Neither __extension__
36// nor #pragma diagnostic ignored work :(
37//
38#pragma GCC system_header
39#endif
40
41namespace boost{ namespace math{
42
43namespace detail
44{
45 // Functor expm1_series returns the next term in the Taylor series
46 // x^k / k!
47 // each time that operator() is invoked.
48 //
49 template <class T>
50 struct expm1_series
51 {
52 typedef T result_type;
53
54 expm1_series(T x)
55 : k(0), m_x(x), m_term(1) {}
56
57 T operator()()
58 {
59 ++k;
60 m_term *= m_x;
61 m_term /= k;
62 return m_term;
63 }
64
65 int count()const
66 {
67 return k;
68 }
69
70 private:
71 int k;
72 const T m_x;
73 T m_term;
74 expm1_series(const expm1_series&);
75 expm1_series& operator=(const expm1_series&);
76 };
77
78template <class T, class Policy, class tag>
79struct expm1_initializer
80{
81 struct init
82 {
83 init()
84 {
85 do_init(tag());
86 }
87 template <int N>
88 static void do_init(const boost::integral_constant<int, N>&){}
89 static void do_init(const boost::integral_constant<int, 64>&)
90 {
91 expm1(T(0.5));
92 }
93 static void do_init(const boost::integral_constant<int, 113>&)
94 {
95 expm1(T(0.5));
96 }
97 void force_instantiate()const{}
98 };
99 static const init initializer;
100 static void force_instantiate()
101 {
102 initializer.force_instantiate();
103 }
104};
105
106template <class T, class Policy, class tag>
107const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer;
108
109//
110// Algorithm expm1 is part of C99, but is not yet provided by many compilers.
111//
112// This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
113//
114template <class T, class Policy>
115T expm1_imp(T x, const boost::integral_constant<int, 0>&, const Policy& pol)
116{
117 BOOST_MATH_STD_USING
118
119 T a = fabs(x);
120 if((boost::math::isnan)(a))
121 {
122 return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
123 }
124 if(a > T(0.5f))
125 {
126 if(a >= tools::log_max_value<T>())
127 {
128 if(x > 0)
129 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
130 return -1;
131 }
132 return exp(x) - T(1);
133 }
134 if(a < tools::epsilon<T>())
135 return x;
136 detail::expm1_series<T> s(x);
137 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
138#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
139 T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
140#else
141 T zero = 0;
142 T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
143#endif
144 policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
145 return result;
146}
147
148template <class T, class P>
149T expm1_imp(T x, const boost::integral_constant<int, 53>&, const P& pol)
150{
151 BOOST_MATH_STD_USING
152
153 T a = fabs(x);
154 if(a > T(0.5L))
155 {
156 if(a >= tools::log_max_value<T>())
157 {
158 if(x > 0)
159 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
160 return -1;
161 }
162 return exp(x) - T(1);
163 }
164 if(a < tools::epsilon<T>())
165 return x;
166
167 static const float Y = 0.10281276702880859e1f;
168 static const T n[] = { static_cast<T>(-0.28127670288085937e-1), static_cast<T>(0.51278186299064534e0), static_cast<T>(-0.6310029069350198e-1), static_cast<T>(0.11638457975729296e-1), static_cast<T>(-0.52143390687521003e-3), static_cast<T>(0.21491399776965688e-4) };
169 static const T d[] = { 1, static_cast<T>(-0.45442309511354755e0), static_cast<T>(0.90850389570911714e-1), static_cast<T>(-0.10088963629815502e-1), static_cast<T>(0.63003407478692265e-3), static_cast<T>(-0.17976570003654402e-4) };
170
171 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
172 return result;
173}
174
175template <class T, class P>
176T expm1_imp(T x, const boost::integral_constant<int, 64>&, const P& pol)
177{
178 BOOST_MATH_STD_USING
179
180 T a = fabs(x);
181 if(a > T(0.5L))
182 {
183 if(a >= tools::log_max_value<T>())
184 {
185 if(x > 0)
186 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
187 return -1;
188 }
189 return exp(x) - T(1);
190 }
191 if(a < tools::epsilon<T>())
192 return x;
193
194 static const float Y = 0.10281276702880859375e1f;
195 static const T n[] = {
196 BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),
197 BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),
198 BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
199 BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
200 BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
201 BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
202 BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
203 };
204 static const T d[] = {
205 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
206 BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
207 BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
208 BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
209 BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
210 BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
211 BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
212 };
213
214 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
215 return result;
216}
217
218template <class T, class P>
219T expm1_imp(T x, const boost::integral_constant<int, 113>&, const P& pol)
220{
221 BOOST_MATH_STD_USING
222
223 T a = fabs(x);
224 if(a > T(0.5L))
225 {
226 if(a >= tools::log_max_value<T>())
227 {
228 if(x > 0)
229 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
230 return -1;
231 }
232 return exp(x) - T(1);
233 }
234 if(a < tools::epsilon<T>())
235 return x;
236
237 static const float Y = 0.10281276702880859375e1f;
238 static const T n[] = {
239 BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
240 BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
241 BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
242 BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
243 BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
244 BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
245 BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
246 BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
247 BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
248 BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
249 };
250 static const T d[] = {
251 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
252 BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
253 BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
254 BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
255 BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
256 BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
257 BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
258 BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
259 BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
260 BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
261 BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
262 };
263
264 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
265 return result;
266}
267
268} // namespace detail
269
270template <class T, class Policy>
271inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
272{
273 typedef typename tools::promote_args<T>::type result_type;
274 typedef typename policies::evaluation<result_type, Policy>::type value_type;
275 typedef typename policies::precision<result_type, Policy>::type precision_type;
276 typedef typename policies::normalise<
277 Policy,
278 policies::promote_float<false>,
279 policies::promote_double<false>,
280 policies::discrete_quantile<>,
281 policies::assert_undefined<> >::type forwarding_policy;
282
283 typedef boost::integral_constant<int,
284 precision_type::value <= 0 ? 0 :
285 precision_type::value <= 53 ? 53 :
286 precision_type::value <= 64 ? 64 :
287 precision_type::value <= 113 ? 113 : 0
288 > tag_type;
289
290 detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
291
292 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
293 static_cast<value_type>(x),
294 tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
295}
296
297#ifdef expm1
298# ifndef BOOST_HAS_expm1
299# define BOOST_HAS_expm1
300# endif
301# undef expm1
302#endif
303
304#if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER))
305# ifdef BOOST_MATH_USE_C99
306inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x: x); }
307# ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
308inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x: x); }
309# endif
310# else
311inline float expm1(float x, const policies::policy<>&){ return static_cast<float>(::expm1(x)); }
312# endif
313inline double expm1(double x, const policies::policy<>&){ return ::expm1(x: x); }
314#endif
315
316template <class T>
317inline typename tools::promote_args<T>::type expm1(T x)
318{
319 return expm1(x, policies::policy<>());
320}
321
322#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))
323inline float expm1(float z)
324{
325 return expm1<float>(z);
326}
327inline double expm1(double z)
328{
329 return expm1<double>(z);
330}
331#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
332inline long double expm1(long double z)
333{
334 return expm1<long double>(z);
335}
336#endif
337#endif
338
339} // namespace math
340} // namespace boost
341
342#endif // BOOST_MATH_HYPOT_INCLUDED
343
344
345
346
347

source code of include/boost/math/special_functions/expm1.hpp