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