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/policies/error_handling.hpp>
20 #include <boost/math/tools/rational.hpp>
21 #include <boost/math/special_functions/math_fwd.hpp>
22 #include <boost/mpl/less_equal.hpp>
23
24 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
25 # include <boost/static_assert.hpp>
26 #else
27 # include <boost/assert.hpp>
28 #endif
29
30 namespace boost{ namespace math{
31
32 namespace detail
33 {
34 // Functor expm1_series returns the next term in the Taylor series
35 // x^k / k!
36 // each time that operator() is invoked.
37 //
38 template <class T>
39 struct expm1_series
40 {
41 typedef T result_type;
42
expm1_seriesboost::math::detail::expm1_series43 expm1_series(T x)
44 : k(0), m_x(x), m_term(1) {}
45
operator ()boost::math::detail::expm1_series46 T operator()()
47 {
48 ++k;
49 m_term *= m_x;
50 m_term /= k;
51 return m_term;
52 }
53
countboost::math::detail::expm1_series54 int count()const
55 {
56 return k;
57 }
58
59 private:
60 int k;
61 const T m_x;
62 T m_term;
63 expm1_series(const expm1_series&);
64 expm1_series& operator=(const expm1_series&);
65 };
66
67 //
68 // Algorithm expm1 is part of C99, but is not yet provided by many compilers.
69 //
70 // This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
71 //
72 template <class T, class Policy>
73 T expm1_imp(T x, const mpl::int_<0>&, const Policy& pol)
74 {
75 BOOST_MATH_STD_USING
76
77 T a = fabs(x);
78 if(a > T(0.5f))
79 {
80 if(a >= tools::log_max_value<T>())
81 {
82 if(x > 0)
83 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
84 return -1;
85 }
86 return exp(x) - T(1);
87 }
88 if(a < tools::epsilon<T>())
89 return x;
90 detail::expm1_series<T> s(x);
91 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
92 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
93 T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
94 #else
95 T zero = 0;
96 T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
97 #endif
98 policies::check_series_iterations("boost::math::expm1<%1%>(%1%)", max_iter, pol);
99 return result;
100 }
101
102 template <class T, class P>
103 T expm1_imp(T x, const mpl::int_<53>&, const P& pol)
104 {
105 BOOST_MATH_STD_USING
106
107 T a = fabs(x);
108 if(a > T(0.5L))
109 {
110 if(a >= tools::log_max_value<T>())
111 {
112 if(x > 0)
113 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
114 return -1;
115 }
116 return exp(x) - T(1);
117 }
118 if(a < tools::epsilon<T>())
119 return x;
120
121 static const float Y = 0.10281276702880859e1f;
122 static const T n[] = { -0.28127670288085937e-1, 0.51278186299064534e0, -0.6310029069350198e-1, 0.11638457975729296e-1, -0.52143390687521003e-3, 0.21491399776965688e-4 };
123 static const T d[] = { 1, -0.45442309511354755e0, 0.90850389570911714e-1, -0.10088963629815502e-1, 0.63003407478692265e-3, -0.17976570003654402e-4 };
124
125 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
126 return result;
127 }
128
129 template <class T, class P>
130 T expm1_imp(T x, const mpl::int_<64>&, const P& pol)
131 {
132 BOOST_MATH_STD_USING
133
134 T a = fabs(x);
135 if(a > T(0.5L))
136 {
137 if(a >= tools::log_max_value<T>())
138 {
139 if(x > 0)
140 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
141 return -1;
142 }
143 return exp(x) - T(1);
144 }
145 if(a < tools::epsilon<T>())
146 return x;
147
148 static const float Y = 0.10281276702880859375e1f;
149 static const T n[] = {
150 -0.281276702880859375e-1L,
151 0.512980290285154286358e0L,
152 -0.667758794592881019644e-1L,
153 0.131432469658444745835e-1L,
154 -0.72303795326880286965e-3L,
155 0.447441185192951335042e-4L,
156 -0.714539134024984593011e-6L
157 };
158 static const T d[] = {
159 1,
160 -0.461477618025562520389e0L,
161 0.961237488025708540713e-1L,
162 -0.116483957658204450739e-1L,
163 0.873308008461557544458e-3L,
164 -0.387922804997682392562e-4L,
165 0.807473180049193557294e-6L
166 };
167
168 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
169 return result;
170 }
171
172 template <class T, class P>
173 T expm1_imp(T x, const mpl::int_<113>&, const P& pol)
174 {
175 BOOST_MATH_STD_USING
176
177 T a = fabs(x);
178 if(a > T(0.5L))
179 {
180 if(a >= tools::log_max_value<T>())
181 {
182 if(x > 0)
183 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
184 return -1;
185 }
186 return exp(x) - T(1);
187 }
188 if(a < tools::epsilon<T>())
189 return x;
190
191 static const float Y = 0.10281276702880859375e1f;
192 static const T n[] = {
193 -0.28127670288085937499999999999999999854e-1L,
194 0.51278156911210477556524452177540792214e0L,
195 -0.63263178520747096729500254678819588223e-1L,
196 0.14703285606874250425508446801230572252e-1L,
197 -0.8675686051689527802425310407898459386e-3L,
198 0.88126359618291165384647080266133492399e-4L,
199 -0.25963087867706310844432390015463138953e-5L,
200 0.14226691087800461778631773363204081194e-6L,
201 -0.15995603306536496772374181066765665596e-8L,
202 0.45261820069007790520447958280473183582e-10L
203 };
204 static const T d[] = {
205 1,
206 -0.45441264709074310514348137469214538853e0L,
207 0.96827131936192217313133611655555298106e-1L,
208 -0.12745248725908178612540554584374876219e-1L,
209 0.11473613871583259821612766907781095472e-2L,
210 -0.73704168477258911962046591907690764416e-4L,
211 0.34087499397791555759285503797256103259e-5L,
212 -0.11114024704296196166272091230695179724e-6L,
213 0.23987051614110848595909588343223896577e-8L,
214 -0.29477341859111589208776402638429026517e-10L,
215 0.13222065991022301420255904060628100924e-12L
216 };
217
218 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
219 return result;
220 }
221
222 } // namespace detail
223
224 template <class T, class Policy>
expm1(T x,const Policy &)225 inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
226 {
227 typedef typename tools::promote_args<T>::type result_type;
228 typedef typename policies::evaluation<result_type, Policy>::type value_type;
229 typedef typename policies::precision<result_type, Policy>::type precision_type;
230 typedef typename policies::normalise<
231 Policy,
232 policies::promote_float<false>,
233 policies::promote_double<false>,
234 policies::discrete_quantile<>,
235 policies::assert_undefined<> >::type forwarding_policy;
236
237 typedef typename mpl::if_c<
238 ::std::numeric_limits<result_type>::is_specialized == 0,
239 mpl::int_<0>, // no numeric_limits, use generic solution
240 typename mpl::if_<
241 typename mpl::less_equal<precision_type, mpl::int_<53> >::type,
242 mpl::int_<53>, // double
243 typename mpl::if_<
244 typename mpl::less_equal<precision_type, mpl::int_<64> >::type,
245 mpl::int_<64>, // 80-bit long double
246 typename mpl::if_<
247 typename mpl::less_equal<precision_type, mpl::int_<113> >::type,
248 mpl::int_<113>, // 128-bit long double
249 mpl::int_<0> // too many bits, use generic version.
250 >::type
251 >::type
252 >::type
253 >::type tag_type;
254
255 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
256 static_cast<value_type>(x),
257 tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
258 }
259
260 #ifdef expm1
261 # ifndef BOOST_HAS_expm1
262 # define BOOST_HAS_expm1
263 # endif
264 # undef expm1
265 #endif
266
267 #if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER))
268 # ifdef BOOST_MATH_USE_C99
expm1(float x,const policies::policy<> &)269 inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x); }
270 # ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
expm1(long double x,const policies::policy<> &)271 inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x); }
272 # endif
273 # else
expm1(float x,const policies::policy<> &)274 inline float expm1(float x, const policies::policy<>&){ return ::expm1(x); }
275 # endif
expm1(double x,const policies::policy<> &)276 inline double expm1(double x, const policies::policy<>&){ return ::expm1(x); }
277 #endif
278
279 template <class T>
expm1(T x)280 inline typename tools::promote_args<T>::type expm1(T x)
281 {
282 return expm1(x, policies::policy<>());
283 }
284
285 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))
expm1(float z)286 inline float expm1(float z)
287 {
288 return expm1<float>(z);
289 }
expm1(double z)290 inline double expm1(double z)
291 {
292 return expm1<double>(z);
293 }
294 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
expm1(long double z)295 inline long double expm1(long double z)
296 {
297 return expm1<long double>(z);
298 }
299 #endif
300 #endif
301
302 } // namespace math
303 } // namespace boost
304
305 #endif // BOOST_MATH_HYPOT_INCLUDED
306
307
308
309
310