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