]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2005-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_LOG1P_INCLUDED | |
7 | #define BOOST_MATH_LOG1P_INCLUDED | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #pragma warning(push) | |
12 | #pragma warning(disable:4702) // Unreachable code (release mode only warning) | |
13 | #endif | |
14 | ||
15 | #include <boost/config/no_tr1/cmath.hpp> | |
16 | #include <math.h> // platform's ::log1p | |
17 | #include <boost/limits.hpp> | |
18 | #include <boost/math/tools/config.hpp> | |
19 | #include <boost/math/tools/series.hpp> | |
20 | #include <boost/math/tools/rational.hpp> | |
21 | #include <boost/math/tools/big_constant.hpp> | |
22 | #include <boost/math/policies/error_handling.hpp> | |
23 | #include <boost/math/special_functions/math_fwd.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 log1p_series returns the next term in the Taylor series | |
36 | // pow(-1, k-1)*pow(x, k) / k | |
37 | // each time that operator() is invoked. | |
38 | // | |
39 | template <class T> | |
40 | struct log1p_series | |
41 | { | |
42 | typedef T result_type; | |
43 | ||
44 | log1p_series(T x) | |
45 | : k(0), m_mult(-x), m_prod(-1){} | |
46 | ||
47 | T operator()() | |
48 | { | |
49 | m_prod *= m_mult; | |
50 | return m_prod / ++k; | |
51 | } | |
52 | ||
53 | int count()const | |
54 | { | |
55 | return k; | |
56 | } | |
57 | ||
58 | private: | |
59 | int k; | |
60 | const T m_mult; | |
61 | T m_prod; | |
62 | log1p_series(const log1p_series&); | |
63 | log1p_series& operator=(const log1p_series&); | |
64 | }; | |
65 | ||
66 | // Algorithm log1p is part of C99, but is not yet provided by many compilers. | |
67 | // | |
68 | // This version uses a Taylor series expansion for 0.5 > x > epsilon, which may | |
69 | // require up to std::numeric_limits<T>::digits+1 terms to be calculated. | |
70 | // It would be much more efficient to use the equivalence: | |
71 | // log(1+x) == (log(1+x) * x) / ((1-x) - 1) | |
72 | // Unfortunately many optimizing compilers make such a mess of this, that | |
73 | // it performs no better than log(1+x): which is to say not very well at all. | |
74 | // | |
75 | template <class T, class Policy> | |
76 | T log1p_imp(T const & x, const Policy& pol, const mpl::int_<0>&) | |
77 | { // The function returns the natural logarithm of 1 + x. | |
78 | typedef typename tools::promote_args<T>::type result_type; | |
79 | BOOST_MATH_STD_USING | |
80 | ||
81 | static const char* function = "boost::math::log1p<%1%>(%1%)"; | |
82 | ||
b32b8144 | 83 | if((x < -1) || (boost::math::isnan)(x)) |
7c673cae FG |
84 | return policies::raise_domain_error<T>( |
85 | function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
86 | if(x == -1) | |
87 | return -policies::raise_overflow_error<T>( | |
88 | function, 0, pol); | |
89 | ||
90 | result_type a = abs(result_type(x)); | |
91 | if(a > result_type(0.5f)) | |
92 | return log(1 + result_type(x)); | |
93 | // Note that without numeric_limits specialisation support, | |
94 | // epsilon just returns zero, and our "optimisation" will always fail: | |
95 | if(a < tools::epsilon<result_type>()) | |
96 | return x; | |
97 | detail::log1p_series<result_type> s(x); | |
98 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
99 | #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245) | |
100 | result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter); | |
101 | #else | |
102 | result_type zero = 0; | |
103 | result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter, zero); | |
104 | #endif | |
105 | policies::check_series_iterations<T>(function, max_iter, pol); | |
106 | return result; | |
107 | } | |
108 | ||
109 | template <class T, class Policy> | |
110 | T log1p_imp(T const& x, const Policy& pol, const mpl::int_<53>&) | |
111 | { // The function returns the natural logarithm of 1 + x. | |
112 | BOOST_MATH_STD_USING | |
113 | ||
114 | static const char* function = "boost::math::log1p<%1%>(%1%)"; | |
115 | ||
116 | if(x < -1) | |
117 | return policies::raise_domain_error<T>( | |
118 | function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
119 | if(x == -1) | |
120 | return -policies::raise_overflow_error<T>( | |
121 | function, 0, pol); | |
122 | ||
123 | T a = fabs(x); | |
124 | if(a > 0.5f) | |
125 | return log(1 + x); | |
126 | // Note that without numeric_limits specialisation support, | |
127 | // epsilon just returns zero, and our "optimisation" will always fail: | |
128 | if(a < tools::epsilon<T>()) | |
129 | return x; | |
130 | ||
131 | // Maximum Deviation Found: 1.846e-017 | |
132 | // Expected Error Term: 1.843e-017 | |
133 | // Maximum Relative Change in Control Points: 8.138e-004 | |
134 | // Max Error found at double precision = 3.250766e-016 | |
135 | static const T P[] = { | |
136 | 0.15141069795941984e-16L, | |
137 | 0.35495104378055055e-15L, | |
138 | 0.33333333333332835L, | |
139 | 0.99249063543365859L, | |
140 | 1.1143969784156509L, | |
141 | 0.58052937949269651L, | |
142 | 0.13703234928513215L, | |
143 | 0.011294864812099712L | |
144 | }; | |
145 | static const T Q[] = { | |
146 | 1L, | |
147 | 3.7274719063011499L, | |
148 | 5.5387948649720334L, | |
149 | 4.159201143419005L, | |
150 | 1.6423855110312755L, | |
151 | 0.31706251443180914L, | |
152 | 0.022665554431410243L, | |
153 | -0.29252538135177773e-5L | |
154 | }; | |
155 | ||
156 | T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x); | |
157 | result *= x; | |
158 | ||
159 | return result; | |
160 | } | |
161 | ||
162 | template <class T, class Policy> | |
163 | T log1p_imp(T const& x, const Policy& pol, const mpl::int_<64>&) | |
164 | { // The function returns the natural logarithm of 1 + x. | |
165 | BOOST_MATH_STD_USING | |
166 | ||
167 | static const char* function = "boost::math::log1p<%1%>(%1%)"; | |
168 | ||
169 | if(x < -1) | |
170 | return policies::raise_domain_error<T>( | |
171 | function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
172 | if(x == -1) | |
173 | return -policies::raise_overflow_error<T>( | |
174 | function, 0, pol); | |
175 | ||
176 | T a = fabs(x); | |
177 | if(a > 0.5f) | |
178 | return log(1 + x); | |
179 | // Note that without numeric_limits specialisation support, | |
180 | // epsilon just returns zero, and our "optimisation" will always fail: | |
181 | if(a < tools::epsilon<T>()) | |
182 | return x; | |
183 | ||
184 | // Maximum Deviation Found: 8.089e-20 | |
185 | // Expected Error Term: 8.088e-20 | |
186 | // Maximum Relative Change in Control Points: 9.648e-05 | |
187 | // Max Error found at long double precision = 2.242324e-19 | |
188 | static const T P[] = { | |
189 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.807533446680736736712e-19), | |
190 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.490881544804798926426e-18), | |
191 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.333333333333333373941), | |
192 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.17141290782087994162), | |
193 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.62790522814926264694), | |
194 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.13156411870766876113), | |
195 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.408087379932853785336), | |
196 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0706537026422828914622), | |
197 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00441709903782239229447) | |
198 | }; | |
199 | static const T Q[] = { | |
200 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
201 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.26423872346263928361), | |
202 | BOOST_MATH_BIG_CONSTANT(T, 64, 7.48189472704477708962), | |
203 | BOOST_MATH_BIG_CONSTANT(T, 64, 6.94757016732904280913), | |
204 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.6493508622280767304), | |
205 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.06884863623790638317), | |
206 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.158292216998514145947), | |
207 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00885295524069924328658), | |
208 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.560026216133415663808e-6) | |
209 | }; | |
210 | ||
211 | T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x); | |
212 | result *= x; | |
213 | ||
214 | return result; | |
215 | } | |
216 | ||
217 | template <class T, class Policy> | |
218 | T log1p_imp(T const& x, const Policy& pol, const mpl::int_<24>&) | |
219 | { // The function returns the natural logarithm of 1 + x. | |
220 | BOOST_MATH_STD_USING | |
221 | ||
222 | static const char* function = "boost::math::log1p<%1%>(%1%)"; | |
223 | ||
224 | if(x < -1) | |
225 | return policies::raise_domain_error<T>( | |
226 | function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
227 | if(x == -1) | |
228 | return -policies::raise_overflow_error<T>( | |
229 | function, 0, pol); | |
230 | ||
231 | T a = fabs(x); | |
232 | if(a > 0.5f) | |
233 | return log(1 + x); | |
234 | // Note that without numeric_limits specialisation support, | |
235 | // epsilon just returns zero, and our "optimisation" will always fail: | |
236 | if(a < tools::epsilon<T>()) | |
237 | return x; | |
238 | ||
239 | // Maximum Deviation Found: 6.910e-08 | |
240 | // Expected Error Term: 6.910e-08 | |
241 | // Maximum Relative Change in Control Points: 2.509e-04 | |
242 | // Max Error found at double precision = 6.910422e-08 | |
243 | // Max Error found at float precision = 8.357242e-08 | |
244 | static const T P[] = { | |
245 | -0.671192866803148236519e-7L, | |
246 | 0.119670999140731844725e-6L, | |
247 | 0.333339469182083148598L, | |
248 | 0.237827183019664122066L | |
249 | }; | |
250 | static const T Q[] = { | |
251 | 1L, | |
252 | 1.46348272586988539733L, | |
253 | 0.497859871350117338894L, | |
254 | -0.00471666268910169651936L | |
255 | }; | |
256 | ||
257 | T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x); | |
258 | result *= x; | |
259 | ||
260 | return result; | |
261 | } | |
262 | ||
263 | template <class T, class Policy, class tag> | |
264 | struct log1p_initializer | |
265 | { | |
266 | struct init | |
267 | { | |
268 | init() | |
269 | { | |
270 | do_init(tag()); | |
271 | } | |
272 | template <int N> | |
273 | static void do_init(const mpl::int_<N>&){} | |
274 | static void do_init(const mpl::int_<64>&) | |
275 | { | |
276 | boost::math::log1p(static_cast<T>(0.25), Policy()); | |
277 | } | |
278 | void force_instantiate()const{} | |
279 | }; | |
280 | static const init initializer; | |
281 | static void force_instantiate() | |
282 | { | |
283 | initializer.force_instantiate(); | |
284 | } | |
285 | }; | |
286 | ||
287 | template <class T, class Policy, class tag> | |
288 | const typename log1p_initializer<T, Policy, tag>::init log1p_initializer<T, Policy, tag>::initializer; | |
289 | ||
290 | ||
291 | } // namespace detail | |
292 | ||
293 | template <class T, class Policy> | |
294 | inline typename tools::promote_args<T>::type log1p(T x, const Policy&) | |
295 | { | |
296 | typedef typename tools::promote_args<T>::type result_type; | |
297 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
298 | typedef typename policies::precision<result_type, Policy>::type precision_type; | |
299 | typedef typename policies::normalise< | |
300 | Policy, | |
301 | policies::promote_float<false>, | |
302 | policies::promote_double<false>, | |
303 | policies::discrete_quantile<>, | |
304 | policies::assert_undefined<> >::type forwarding_policy; | |
305 | ||
306 | typedef typename mpl::if_< | |
307 | mpl::less_equal<precision_type, mpl::int_<0> >, | |
308 | mpl::int_<0>, | |
309 | typename mpl::if_< | |
310 | mpl::less_equal<precision_type, mpl::int_<53> >, | |
311 | mpl::int_<53>, // double | |
312 | typename mpl::if_< | |
313 | mpl::less_equal<precision_type, mpl::int_<64> >, | |
314 | mpl::int_<64>, // 80-bit long double | |
315 | mpl::int_<0> // too many bits, use generic version. | |
316 | >::type | |
317 | >::type | |
318 | >::type tag_type; | |
319 | ||
320 | detail::log1p_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); | |
321 | ||
322 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
323 | detail::log1p_imp(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)"); | |
324 | } | |
325 | ||
326 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564)) | |
327 | // These overloads work around a type deduction bug: | |
328 | inline float log1p(float z) | |
329 | { | |
330 | return log1p<float>(z); | |
331 | } | |
332 | inline double log1p(double z) | |
333 | { | |
334 | return log1p<double>(z); | |
335 | } | |
336 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
337 | inline long double log1p(long double z) | |
338 | { | |
339 | return log1p<long double>(z); | |
340 | } | |
341 | #endif | |
342 | #endif | |
343 | ||
344 | #ifdef log1p | |
345 | # ifndef BOOST_HAS_LOG1P | |
346 | # define BOOST_HAS_LOG1P | |
347 | # endif | |
348 | # undef log1p | |
349 | #endif | |
350 | ||
351 | #if defined(BOOST_HAS_LOG1P) && !(defined(__osf__) && defined(__DECCXX_VER)) | |
352 | # ifdef BOOST_MATH_USE_C99 | |
353 | template <class Policy> | |
354 | inline float log1p(float x, const Policy& pol) | |
355 | { | |
356 | if(x < -1) | |
357 | return policies::raise_domain_error<float>( | |
358 | "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
359 | if(x == -1) | |
360 | return -policies::raise_overflow_error<float>( | |
361 | "log1p<%1%>(%1%)", 0, pol); | |
362 | return ::log1pf(x); | |
363 | } | |
364 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
365 | template <class Policy> | |
366 | inline long double log1p(long double x, const Policy& pol) | |
367 | { | |
368 | if(x < -1) | |
369 | return policies::raise_domain_error<long double>( | |
370 | "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
371 | if(x == -1) | |
372 | return -policies::raise_overflow_error<long double>( | |
373 | "log1p<%1%>(%1%)", 0, pol); | |
374 | return ::log1pl(x); | |
375 | } | |
376 | #endif | |
377 | #else | |
378 | template <class Policy> | |
379 | inline float log1p(float x, const Policy& pol) | |
380 | { | |
381 | if(x < -1) | |
382 | return policies::raise_domain_error<float>( | |
383 | "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
384 | if(x == -1) | |
385 | return -policies::raise_overflow_error<float>( | |
386 | "log1p<%1%>(%1%)", 0, pol); | |
387 | return ::log1p(x); | |
388 | } | |
389 | #endif | |
390 | template <class Policy> | |
391 | inline double log1p(double x, const Policy& pol) | |
392 | { | |
393 | if(x < -1) | |
394 | return policies::raise_domain_error<double>( | |
395 | "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
396 | if(x == -1) | |
397 | return -policies::raise_overflow_error<double>( | |
398 | "log1p<%1%>(%1%)", 0, pol); | |
399 | return ::log1p(x); | |
400 | } | |
401 | #elif defined(_MSC_VER) && (BOOST_MSVC >= 1400) | |
402 | // | |
403 | // You should only enable this branch if you are absolutely sure | |
404 | // that your compilers optimizer won't mess this code up!! | |
405 | // Currently tested with VC8 and Intel 9.1. | |
406 | // | |
407 | template <class Policy> | |
408 | inline double log1p(double x, const Policy& pol) | |
409 | { | |
410 | if(x < -1) | |
411 | return policies::raise_domain_error<double>( | |
412 | "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
413 | if(x == -1) | |
414 | return -policies::raise_overflow_error<double>( | |
415 | "log1p<%1%>(%1%)", 0, pol); | |
416 | double u = 1+x; | |
417 | if(u == 1.0) | |
418 | return x; | |
419 | else | |
420 | return ::log(u)*(x/(u-1.0)); | |
421 | } | |
422 | template <class Policy> | |
423 | inline float log1p(float x, const Policy& pol) | |
424 | { | |
425 | return static_cast<float>(boost::math::log1p(static_cast<double>(x), pol)); | |
426 | } | |
427 | #ifndef _WIN32_WCE | |
428 | // | |
429 | // For some reason this fails to compile under WinCE... | |
430 | // Needs more investigation. | |
431 | // | |
432 | template <class Policy> | |
433 | inline long double log1p(long double x, const Policy& pol) | |
434 | { | |
435 | if(x < -1) | |
436 | return policies::raise_domain_error<long double>( | |
437 | "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); | |
438 | if(x == -1) | |
439 | return -policies::raise_overflow_error<long double>( | |
440 | "log1p<%1%>(%1%)", 0, pol); | |
441 | long double u = 1+x; | |
442 | if(u == 1.0) | |
443 | return x; | |
444 | else | |
445 | return ::logl(u)*(x/(u-1.0)); | |
446 | } | |
447 | #endif | |
448 | #endif | |
449 | ||
450 | template <class T> | |
451 | inline typename tools::promote_args<T>::type log1p(T x) | |
452 | { | |
453 | return boost::math::log1p(x, policies::policy<>()); | |
454 | } | |
455 | // | |
456 | // Compute log(1+x)-x: | |
457 | // | |
458 | template <class T, class Policy> | |
459 | inline typename tools::promote_args<T>::type | |
460 | log1pmx(T x, const Policy& pol) | |
461 | { | |
462 | typedef typename tools::promote_args<T>::type result_type; | |
463 | BOOST_MATH_STD_USING | |
464 | static const char* function = "boost::math::log1pmx<%1%>(%1%)"; | |
465 | ||
466 | if(x < -1) | |
467 | return policies::raise_domain_error<T>( | |
468 | function, "log1pmx(x) requires x > -1, but got x = %1%.", x, pol); | |
469 | if(x == -1) | |
470 | return -policies::raise_overflow_error<T>( | |
471 | function, 0, pol); | |
472 | ||
473 | result_type a = abs(result_type(x)); | |
474 | if(a > result_type(0.95f)) | |
475 | return log(1 + result_type(x)) - result_type(x); | |
476 | // Note that without numeric_limits specialisation support, | |
477 | // epsilon just returns zero, and our "optimisation" will always fail: | |
478 | if(a < tools::epsilon<result_type>()) | |
479 | return -x * x / 2; | |
480 | boost::math::detail::log1p_series<T> s(x); | |
481 | s(); | |
482 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
483 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
484 | T zero = 0; | |
485 | T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero); | |
486 | #else | |
487 | T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter); | |
488 | #endif | |
489 | policies::check_series_iterations<T>(function, max_iter, pol); | |
490 | return result; | |
491 | } | |
492 | ||
493 | template <class T> | |
494 | inline typename tools::promote_args<T>::type log1pmx(T x) | |
495 | { | |
496 | return log1pmx(x, policies::policy<>()); | |
497 | } | |
498 | ||
499 | } // namespace math | |
500 | } // namespace boost | |
501 | ||
502 | #ifdef _MSC_VER | |
503 | #pragma warning(pop) | |
504 | #endif | |
505 | ||
506 | #endif // BOOST_MATH_LOG1P_INCLUDED | |
507 | ||
508 | ||
509 |