Chris@16: // (C) Copyright John Maddock 2005-2006. Chris@16: // Use, modification and distribution are subject to the Chris@16: // Boost Software License, Version 1.0. (See accompanying file Chris@16: // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) Chris@16: Chris@16: #ifndef BOOST_MATH_LOG1P_INCLUDED Chris@16: #define BOOST_MATH_LOG1P_INCLUDED Chris@16: Chris@16: #ifdef _MSC_VER Chris@16: #pragma once Chris@16: #endif Chris@16: Chris@16: #include Chris@16: #include // platform's ::log1p Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS Chris@16: # include Chris@16: #else Chris@16: # include Chris@16: #endif Chris@16: Chris@16: namespace boost{ namespace math{ Chris@16: Chris@16: namespace detail Chris@16: { Chris@16: // Functor log1p_series returns the next term in the Taylor series Chris@16: // pow(-1, k-1)*pow(x, k) / k Chris@16: // each time that operator() is invoked. Chris@16: // Chris@16: template Chris@16: struct log1p_series Chris@16: { Chris@16: typedef T result_type; Chris@16: Chris@16: log1p_series(T x) Chris@16: : k(0), m_mult(-x), m_prod(-1){} Chris@16: Chris@16: T operator()() Chris@16: { Chris@16: m_prod *= m_mult; Chris@16: return m_prod / ++k; Chris@16: } Chris@16: Chris@16: int count()const Chris@16: { Chris@16: return k; Chris@16: } Chris@16: Chris@16: private: Chris@16: int k; Chris@16: const T m_mult; Chris@16: T m_prod; Chris@16: log1p_series(const log1p_series&); Chris@16: log1p_series& operator=(const log1p_series&); Chris@16: }; Chris@16: Chris@16: // Algorithm log1p is part of C99, but is not yet provided by many compilers. Chris@16: // Chris@16: // This version uses a Taylor series expansion for 0.5 > x > epsilon, which may Chris@16: // require up to std::numeric_limits::digits+1 terms to be calculated. Chris@16: // It would be much more efficient to use the equivalence: Chris@16: // log(1+x) == (log(1+x) * x) / ((1-x) - 1) Chris@16: // Unfortunately many optimizing compilers make such a mess of this, that Chris@16: // it performs no better than log(1+x): which is to say not very well at all. Chris@16: // Chris@16: template Chris@16: T log1p_imp(T const & x, const Policy& pol, const mpl::int_<0>&) Chris@16: { // The function returns the natural logarithm of 1 + x. Chris@16: typedef typename tools::promote_args::type result_type; Chris@16: BOOST_MATH_STD_USING Chris@16: Chris@16: static const char* function = "boost::math::log1p<%1%>(%1%)"; Chris@16: Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: function, 0, pol); Chris@16: Chris@16: result_type a = abs(result_type(x)); Chris@16: if(a > result_type(0.5f)) Chris@16: return log(1 + result_type(x)); Chris@16: // Note that without numeric_limits specialisation support, Chris@16: // epsilon just returns zero, and our "optimisation" will always fail: Chris@16: if(a < tools::epsilon()) Chris@16: return x; Chris@16: detail::log1p_series s(x); Chris@16: boost::uintmax_t max_iter = policies::get_max_series_iterations(); Chris@16: #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245) Chris@16: result_type result = tools::sum_series(s, policies::get_epsilon(), max_iter); Chris@16: #else Chris@16: result_type zero = 0; Chris@16: result_type result = tools::sum_series(s, policies::get_epsilon(), max_iter, zero); Chris@16: #endif Chris@16: policies::check_series_iterations(function, max_iter, pol); Chris@16: return result; Chris@16: } Chris@16: Chris@16: template Chris@16: T log1p_imp(T const& x, const Policy& pol, const mpl::int_<53>&) Chris@16: { // The function returns the natural logarithm of 1 + x. Chris@16: BOOST_MATH_STD_USING Chris@16: Chris@16: static const char* function = "boost::math::log1p<%1%>(%1%)"; Chris@16: Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: function, 0, pol); Chris@16: Chris@16: T a = fabs(x); Chris@16: if(a > 0.5f) Chris@16: return log(1 + x); Chris@16: // Note that without numeric_limits specialisation support, Chris@16: // epsilon just returns zero, and our "optimisation" will always fail: Chris@16: if(a < tools::epsilon()) Chris@16: return x; Chris@16: Chris@16: // Maximum Deviation Found: 1.846e-017 Chris@16: // Expected Error Term: 1.843e-017 Chris@16: // Maximum Relative Change in Control Points: 8.138e-004 Chris@16: // Max Error found at double precision = 3.250766e-016 Chris@16: static const T P[] = { Chris@16: 0.15141069795941984e-16L, Chris@16: 0.35495104378055055e-15L, Chris@16: 0.33333333333332835L, Chris@16: 0.99249063543365859L, Chris@16: 1.1143969784156509L, Chris@16: 0.58052937949269651L, Chris@16: 0.13703234928513215L, Chris@16: 0.011294864812099712L Chris@16: }; Chris@16: static const T Q[] = { Chris@16: 1L, Chris@16: 3.7274719063011499L, Chris@16: 5.5387948649720334L, Chris@16: 4.159201143419005L, Chris@16: 1.6423855110312755L, Chris@16: 0.31706251443180914L, Chris@16: 0.022665554431410243L, Chris@16: -0.29252538135177773e-5L Chris@16: }; Chris@16: Chris@16: T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x); Chris@16: result *= x; Chris@16: Chris@16: return result; Chris@16: } Chris@16: Chris@16: template Chris@16: T log1p_imp(T const& x, const Policy& pol, const mpl::int_<64>&) Chris@16: { // The function returns the natural logarithm of 1 + x. Chris@16: BOOST_MATH_STD_USING Chris@16: Chris@16: static const char* function = "boost::math::log1p<%1%>(%1%)"; Chris@16: Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: function, 0, pol); Chris@16: Chris@16: T a = fabs(x); Chris@16: if(a > 0.5f) Chris@16: return log(1 + x); Chris@16: // Note that without numeric_limits specialisation support, Chris@16: // epsilon just returns zero, and our "optimisation" will always fail: Chris@16: if(a < tools::epsilon()) Chris@16: return x; Chris@16: Chris@16: // Maximum Deviation Found: 8.089e-20 Chris@16: // Expected Error Term: 8.088e-20 Chris@16: // Maximum Relative Change in Control Points: 9.648e-05 Chris@16: // Max Error found at long double precision = 2.242324e-19 Chris@16: static const T P[] = { Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, -0.807533446680736736712e-19), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, -0.490881544804798926426e-18), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 0.333333333333333373941), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 1.17141290782087994162), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 1.62790522814926264694), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 1.13156411870766876113), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 0.408087379932853785336), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 0.0706537026422828914622), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 0.00441709903782239229447) Chris@16: }; Chris@16: static const T Q[] = { Chris@101: BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 4.26423872346263928361), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 7.48189472704477708962), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 6.94757016732904280913), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 3.6493508622280767304), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 1.06884863623790638317), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 0.158292216998514145947), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, 0.00885295524069924328658), Chris@16: BOOST_MATH_BIG_CONSTANT(T, 64, -0.560026216133415663808e-6) Chris@16: }; Chris@16: Chris@16: T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x); Chris@16: result *= x; Chris@16: Chris@16: return result; Chris@16: } Chris@16: Chris@16: template Chris@16: T log1p_imp(T const& x, const Policy& pol, const mpl::int_<24>&) Chris@16: { // The function returns the natural logarithm of 1 + x. Chris@16: BOOST_MATH_STD_USING Chris@16: Chris@16: static const char* function = "boost::math::log1p<%1%>(%1%)"; Chris@16: Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: function, "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: function, 0, pol); Chris@16: Chris@16: T a = fabs(x); Chris@16: if(a > 0.5f) Chris@16: return log(1 + x); Chris@16: // Note that without numeric_limits specialisation support, Chris@16: // epsilon just returns zero, and our "optimisation" will always fail: Chris@16: if(a < tools::epsilon()) Chris@16: return x; Chris@16: Chris@16: // Maximum Deviation Found: 6.910e-08 Chris@16: // Expected Error Term: 6.910e-08 Chris@16: // Maximum Relative Change in Control Points: 2.509e-04 Chris@16: // Max Error found at double precision = 6.910422e-08 Chris@16: // Max Error found at float precision = 8.357242e-08 Chris@16: static const T P[] = { Chris@16: -0.671192866803148236519e-7L, Chris@16: 0.119670999140731844725e-6L, Chris@16: 0.333339469182083148598L, Chris@16: 0.237827183019664122066L Chris@16: }; Chris@16: static const T Q[] = { Chris@16: 1L, Chris@16: 1.46348272586988539733L, Chris@16: 0.497859871350117338894L, Chris@16: -0.00471666268910169651936L Chris@16: }; Chris@16: Chris@16: T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x); Chris@16: result *= x; Chris@16: Chris@16: return result; Chris@16: } Chris@16: Chris@16: template Chris@16: struct log1p_initializer Chris@16: { Chris@16: struct init Chris@16: { Chris@16: init() Chris@16: { Chris@16: do_init(tag()); Chris@16: } Chris@16: template Chris@16: static void do_init(const mpl::int_&){} Chris@16: static void do_init(const mpl::int_<64>&) Chris@16: { Chris@16: boost::math::log1p(static_cast(0.25), Policy()); Chris@16: } Chris@16: void force_instantiate()const{} Chris@16: }; Chris@16: static const init initializer; Chris@16: static void force_instantiate() Chris@16: { Chris@16: initializer.force_instantiate(); Chris@16: } Chris@16: }; Chris@16: Chris@16: template Chris@16: const typename log1p_initializer::init log1p_initializer::initializer; Chris@16: Chris@16: Chris@16: } // namespace detail Chris@16: Chris@16: template Chris@16: inline typename tools::promote_args::type log1p(T x, const Policy&) Chris@16: { Chris@16: typedef typename tools::promote_args::type result_type; Chris@16: typedef typename policies::evaluation::type value_type; Chris@16: typedef typename policies::precision::type precision_type; Chris@16: typedef typename policies::normalise< Chris@16: Policy, Chris@16: policies::promote_float, Chris@16: policies::promote_double, Chris@16: policies::discrete_quantile<>, Chris@16: policies::assert_undefined<> >::type forwarding_policy; Chris@16: Chris@16: typedef typename mpl::if_< Chris@16: mpl::less_equal >, Chris@16: mpl::int_<0>, Chris@16: typename mpl::if_< Chris@16: mpl::less_equal >, Chris@16: mpl::int_<53>, // double Chris@16: typename mpl::if_< Chris@16: mpl::less_equal >, Chris@16: mpl::int_<64>, // 80-bit long double Chris@16: mpl::int_<0> // too many bits, use generic version. Chris@16: >::type Chris@16: >::type Chris@16: >::type tag_type; Chris@16: Chris@16: detail::log1p_initializer::force_instantiate(); Chris@16: Chris@16: return policies::checked_narrowing_cast( Chris@16: detail::log1p_imp(static_cast(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)"); Chris@16: } Chris@16: Chris@16: #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564)) Chris@16: // These overloads work around a type deduction bug: Chris@16: inline float log1p(float z) Chris@16: { Chris@16: return log1p(z); Chris@16: } Chris@16: inline double log1p(double z) Chris@16: { Chris@16: return log1p(z); Chris@16: } Chris@16: #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS Chris@16: inline long double log1p(long double z) Chris@16: { Chris@16: return log1p(z); Chris@16: } Chris@16: #endif Chris@16: #endif Chris@16: Chris@16: #ifdef log1p Chris@16: # ifndef BOOST_HAS_LOG1P Chris@16: # define BOOST_HAS_LOG1P Chris@16: # endif Chris@16: # undef log1p Chris@16: #endif Chris@16: Chris@16: #if defined(BOOST_HAS_LOG1P) && !(defined(__osf__) && defined(__DECCXX_VER)) Chris@16: # ifdef BOOST_MATH_USE_C99 Chris@16: template Chris@16: inline float log1p(float x, const Policy& pol) Chris@16: { Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: "log1p<%1%>(%1%)", 0, pol); Chris@16: return ::log1pf(x); Chris@16: } Chris@16: #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS Chris@16: template Chris@16: inline long double log1p(long double x, const Policy& pol) Chris@16: { Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: "log1p<%1%>(%1%)", 0, pol); Chris@16: return ::log1pl(x); Chris@16: } Chris@16: #endif Chris@16: #else Chris@16: template Chris@16: inline float log1p(float x, const Policy& pol) Chris@16: { Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: "log1p<%1%>(%1%)", 0, pol); Chris@16: return ::log1p(x); Chris@16: } Chris@16: #endif Chris@16: template Chris@16: inline double log1p(double x, const Policy& pol) Chris@16: { Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: "log1p<%1%>(%1%)", 0, pol); Chris@16: return ::log1p(x); Chris@16: } Chris@16: #elif defined(_MSC_VER) && (BOOST_MSVC >= 1400) Chris@16: // Chris@16: // You should only enable this branch if you are absolutely sure Chris@16: // that your compilers optimizer won't mess this code up!! Chris@16: // Currently tested with VC8 and Intel 9.1. Chris@16: // Chris@16: template Chris@16: inline double log1p(double x, const Policy& pol) Chris@16: { Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: "log1p<%1%>(%1%)", 0, pol); Chris@16: double u = 1+x; Chris@16: if(u == 1.0) Chris@16: return x; Chris@16: else Chris@16: return ::log(u)*(x/(u-1.0)); Chris@16: } Chris@16: template Chris@16: inline float log1p(float x, const Policy& pol) Chris@16: { Chris@16: return static_cast(boost::math::log1p(static_cast(x), pol)); Chris@16: } Chris@16: #ifndef _WIN32_WCE Chris@16: // Chris@16: // For some reason this fails to compile under WinCE... Chris@16: // Needs more investigation. Chris@16: // Chris@16: template Chris@16: inline long double log1p(long double x, const Policy& pol) Chris@16: { Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: "log1p<%1%>(%1%)", 0, pol); Chris@16: long double u = 1+x; Chris@16: if(u == 1.0) Chris@16: return x; Chris@16: else Chris@16: return ::logl(u)*(x/(u-1.0)); Chris@16: } Chris@16: #endif Chris@16: #endif Chris@16: Chris@16: template Chris@16: inline typename tools::promote_args::type log1p(T x) Chris@16: { Chris@16: return boost::math::log1p(x, policies::policy<>()); Chris@16: } Chris@16: // Chris@16: // Compute log(1+x)-x: Chris@16: // Chris@16: template Chris@16: inline typename tools::promote_args::type Chris@16: log1pmx(T x, const Policy& pol) Chris@16: { Chris@16: typedef typename tools::promote_args::type result_type; Chris@16: BOOST_MATH_STD_USING Chris@16: static const char* function = "boost::math::log1pmx<%1%>(%1%)"; Chris@16: Chris@16: if(x < -1) Chris@16: return policies::raise_domain_error( Chris@16: function, "log1pmx(x) requires x > -1, but got x = %1%.", x, pol); Chris@16: if(x == -1) Chris@16: return -policies::raise_overflow_error( Chris@16: function, 0, pol); Chris@16: Chris@16: result_type a = abs(result_type(x)); Chris@16: if(a > result_type(0.95f)) Chris@16: return log(1 + result_type(x)) - result_type(x); Chris@16: // Note that without numeric_limits specialisation support, Chris@16: // epsilon just returns zero, and our "optimisation" will always fail: Chris@16: if(a < tools::epsilon()) Chris@16: return -x * x / 2; Chris@16: boost::math::detail::log1p_series s(x); Chris@16: s(); Chris@16: boost::uintmax_t max_iter = policies::get_max_series_iterations(); Chris@16: #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) Chris@16: T zero = 0; Chris@16: T result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter, zero); Chris@16: #else Chris@16: T result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter); Chris@16: #endif Chris@16: policies::check_series_iterations(function, max_iter, pol); Chris@16: return result; Chris@16: } Chris@16: Chris@16: template Chris@16: inline typename tools::promote_args::type log1pmx(T x) Chris@16: { Chris@16: return log1pmx(x, policies::policy<>()); Chris@16: } Chris@16: Chris@16: } // namespace math Chris@16: } // namespace boost Chris@16: Chris@16: #endif // BOOST_MATH_LOG1P_INCLUDED Chris@16: Chris@16: Chris@16: