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

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