math/special_functions/log1p.hpp

0001 //  (C) Copyright John Maddock 2005-2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_LOG1P_INCLUDED
0007 #define BOOST_MATH_LOG1P_INCLUDED
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #pragma warning(push)
0012 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
0013 #endif
0014
0015 #include <boost/math/tools/config.hpp>
0016 #include <boost/math/tools/series.hpp>
0017 #include <boost/math/tools/rational.hpp>
0018 #include <boost/math/tools/big_constant.hpp>
0019 #include <boost/math/tools/numeric_limits.hpp>
0020 #include <boost/math/tools/cstdint.hpp>
0021 #include <boost/math/tools/promotion.hpp>
0022 #include <boost/math/tools/precision.hpp>
0023 #include <boost/math/policies/error_handling.hpp>
0024 #include <boost/math/special_functions/math_fwd.hpp>
0025 #include <boost/math/tools/assert.hpp>
0026 #include <boost/math/special_functions/fpclassify.hpp>
0027
0028 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0029 //
0030 // This is the only way we can avoid
0031 // warning: non-standard suffix on floating constant [-Wpedantic]
0032 // when building with -Wall -pedantic.  Neither __extension__
0033 // nor #pragma diagnostic ignored work :(
0034 //
0035 #pragma GCC system_header
0036 #endif
0037
0038 namespace boost{ namespace math{
0039
0040 namespace detail
0041 {
0042   // Functor log1p_series returns the next term in the Taylor series
0043   //   pow(-1, k-1)*pow(x, k) / k
0044   // each time that operator() is invoked.
0045   //
0046   template <class T>
0047   struct log1p_series
0048   {
0049      typedef T result_type;
0050
0051      BOOST_MATH_GPU_ENABLED log1p_series(T x)
0052         : k(0), m_mult(-x), m_prod(-1){}
0053
0054      BOOST_MATH_GPU_ENABLED T operator()()
0055      {
0056         m_prod *= m_mult;
0057         return m_prod / ++k;
0058      }
0059
0060      BOOST_MATH_GPU_ENABLED int count()const
0061      {
0062         return k;
0063      }
0064
0065   private:
0066      int k;
0067      const T m_mult;
0068      T m_prod;
0069      log1p_series(const log1p_series&) = delete;
0070      log1p_series& operator=(const log1p_series&) = delete;
0071   };
0072
0073 // Algorithm log1p is part of C99, but is not yet provided by many compilers.
0074 //
0075 // This version uses a Taylor series expansion for 0.5 > x > epsilon, which may
0076 // require up to std::numeric_limits<T>::digits+1 terms to be calculated.
0077 // It would be much more efficient to use the equivalence:
0078 //   log(1+x) == (log(1+x) * x) / ((1-x) - 1)
0079 // Unfortunately many optimizing compilers make such a mess of this, that
0080 // it performs no better than log(1+x): which is to say not very well at all.
0081 //
0082 template <class T, class Policy>
0083 BOOST_MATH_GPU_ENABLED T log1p_imp(T const & x, const Policy& pol, const boost::math::integral_constant<int, 0>&)
0084 { // The function returns the natural logarithm of 1 + x.
0085    typedef typename tools::promote_args<T>::type result_type;
0086    BOOST_MATH_STD_USING
0087
0088    constexpr auto function = "boost::math::log1p<%1%>(%1%)";
0089
0090    if((x < -1) || (boost::math::isnan)(x))
0091       return policies::raise_domain_error<T>(
0092          function, "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0093    if(x == -1)
0094       return -policies::raise_overflow_error<T>(
0095          function, nullptr, pol);
0096
0097    result_type a = abs(result_type(x));
0098    if(a > result_type(0.5f))
0099       return log(1 + result_type(x));
0100    // Note that without numeric_limits specialisation support,
0101    // epsilon just returns zero, and our "optimisation" will always fail:
0102    if(a < tools::epsilon<result_type>())
0103       return x;
0104    detail::log1p_series<result_type> s(x);
0105    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0106
0107    result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter);
0108
0109    policies::check_series_iterations<T>(function, max_iter, pol);
0110    return result;
0111 }
0112
0113 template <class T, class Policy>
0114 BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const boost::math::integral_constant<int, 53>&)
0115 { // The function returns the natural logarithm of 1 + x.
0116    BOOST_MATH_STD_USING
0117
0118    constexpr auto function = "boost::math::log1p<%1%>(%1%)";
0119
0120    if(x < -1)
0121       return policies::raise_domain_error<T>(
0122          function, "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0123    if(x == -1)
0124       return -policies::raise_overflow_error<T>(
0125          function, nullptr, pol);
0126
0127    T a = fabs(x);
0128    if(a > 0.5f)
0129       return log(1 + x);
0130    // Note that without numeric_limits specialisation support,
0131    // epsilon just returns zero, and our "optimisation" will always fail:
0132    if(a < tools::epsilon<T>())
0133       return x;
0134
0135    // Maximum Deviation Found:                     1.846e-017
0136    // Expected Error Term:                         1.843e-017
0137    // Maximum Relative Change in Control Points:   8.138e-004
0138    // Max Error found at double precision =        3.250766e-016
0139    BOOST_MATH_STATIC const T P[] = {
0140        static_cast<T>(0.15141069795941984e-16L),
0141        static_cast<T>(0.35495104378055055e-15L),
0142        static_cast<T>(0.33333333333332835L),
0143        static_cast<T>(0.99249063543365859L),
0144        static_cast<T>(1.1143969784156509L),
0145        static_cast<T>(0.58052937949269651L),
0146        static_cast<T>(0.13703234928513215L),
0147        static_cast<T>(0.011294864812099712L)
0148      };
0149    BOOST_MATH_STATIC const T Q[] = {
0150        static_cast<T>(1L),
0151        static_cast<T>(3.7274719063011499L),
0152        static_cast<T>(5.5387948649720334L),
0153        static_cast<T>(4.159201143419005L),
0154        static_cast<T>(1.6423855110312755L),
0155        static_cast<T>(0.31706251443180914L),
0156        static_cast<T>(0.022665554431410243L),
0157        static_cast<T>(-0.29252538135177773e-5L)
0158      };
0159
0160    T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x);
0161    result *= x;
0162
0163    return result;
0164 }
0165
0166 template <class T, class Policy>
0167 BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const boost::math::integral_constant<int, 64>&)
0168 { // The function returns the natural logarithm of 1 + x.
0169    BOOST_MATH_STD_USING
0170
0171    constexpr auto function = "boost::math::log1p<%1%>(%1%)";
0172
0173    if(x < -1)
0174       return policies::raise_domain_error<T>(
0175          function, "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0176    if(x == -1)
0177       return -policies::raise_overflow_error<T>(
0178          function, nullptr, pol);
0179
0180    T a = fabs(x);
0181    if(a > 0.5f)
0182       return log(1 + x);
0183    // Note that without numeric_limits specialisation support,
0184    // epsilon just returns zero, and our "optimisation" will always fail:
0185    if(a < tools::epsilon<T>())
0186       return x;
0187
0188    // Maximum Deviation Found:                     8.089e-20
0189    // Expected Error Term:                         8.088e-20
0190    // Maximum Relative Change in Control Points:   9.648e-05
0191    // Max Error found at long double precision =   2.242324e-19
0192    BOOST_MATH_STATIC const T P[] = {
0193       BOOST_MATH_BIG_CONSTANT(T, 64, -0.807533446680736736712e-19),
0194       BOOST_MATH_BIG_CONSTANT(T, 64, -0.490881544804798926426e-18),
0195       BOOST_MATH_BIG_CONSTANT(T, 64, 0.333333333333333373941),
0196       BOOST_MATH_BIG_CONSTANT(T, 64, 1.17141290782087994162),
0197       BOOST_MATH_BIG_CONSTANT(T, 64, 1.62790522814926264694),
0198       BOOST_MATH_BIG_CONSTANT(T, 64, 1.13156411870766876113),
0199       BOOST_MATH_BIG_CONSTANT(T, 64, 0.408087379932853785336),
0200       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0706537026422828914622),
0201       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00441709903782239229447)
0202    };
0203    BOOST_MATH_STATIC const T Q[] = {
0204       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0205       BOOST_MATH_BIG_CONSTANT(T, 64, 4.26423872346263928361),
0206       BOOST_MATH_BIG_CONSTANT(T, 64, 7.48189472704477708962),
0207       BOOST_MATH_BIG_CONSTANT(T, 64, 6.94757016732904280913),
0208       BOOST_MATH_BIG_CONSTANT(T, 64, 3.6493508622280767304),
0209       BOOST_MATH_BIG_CONSTANT(T, 64, 1.06884863623790638317),
0210       BOOST_MATH_BIG_CONSTANT(T, 64, 0.158292216998514145947),
0211       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00885295524069924328658),
0212       BOOST_MATH_BIG_CONSTANT(T, 64, -0.560026216133415663808e-6)
0213    };
0214
0215    T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x);
0216    result *= x;
0217
0218    return result;
0219 }
0220
0221 template <class T, class Policy>
0222 BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const boost::math::integral_constant<int, 24>&)
0223 { // The function returns the natural logarithm of 1 + x.
0224    BOOST_MATH_STD_USING
0225
0226    constexpr auto function = "boost::math::log1p<%1%>(%1%)";
0227
0228    if(x < -1)
0229       return policies::raise_domain_error<T>(
0230          function, "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0231    if(x == -1)
0232       return -policies::raise_overflow_error<T>(
0233          function, nullptr, pol);
0234
0235    T a = fabs(x);
0236    if(a > 0.5f)
0237       return log(1 + x);
0238    // Note that without numeric_limits specialisation support,
0239    // epsilon just returns zero, and our "optimisation" will always fail:
0240    if(a < tools::epsilon<T>())
0241       return x;
0242
0243    // Maximum Deviation Found:                     6.910e-08
0244    // Expected Error Term:                         6.910e-08
0245    // Maximum Relative Change in Control Points:   2.509e-04
0246    // Max Error found at double precision =        6.910422e-08
0247    // Max Error found at float precision =         8.357242e-08
0248    BOOST_MATH_STATIC const T P[] = {
0249       -0.671192866803148236519e-7L,
0250       0.119670999140731844725e-6L,
0251       0.333339469182083148598L,
0252       0.237827183019664122066L
0253    };
0254    BOOST_MATH_STATIC const T Q[] = {
0255       1L,
0256       1.46348272586988539733L,
0257       0.497859871350117338894L,
0258       -0.00471666268910169651936L
0259    };
0260
0261    T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x);
0262    result *= x;
0263
0264    return result;
0265 }
0266
0267 template <class T, class Policy, class tag>
0268 struct log1p_initializer
0269 {
0270    struct init
0271    {
0272       BOOST_MATH_GPU_ENABLED init()
0273       {
0274          do_init(tag());
0275       }
0276       template <int N>
0277       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, N>&){}
0278       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 64>&)
0279       {
0280          boost::math::log1p(static_cast<T>(0.25), Policy());
0281       }
0282       BOOST_MATH_GPU_ENABLED void force_instantiate()const{}
0283    };
0284    BOOST_MATH_STATIC const init initializer;
0285    BOOST_MATH_GPU_ENABLED static void force_instantiate()
0286    {
0287       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0288       initializer.force_instantiate();
0289       #endif
0290    }
0291 };
0292
0293 template <class T, class Policy, class tag>
0294 const typename log1p_initializer<T, Policy, tag>::init log1p_initializer<T, Policy, tag>::initializer;
0295
0296
0297 } // namespace detail
0298
0299 template <class T, class Policy>
0300 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1p(T x, const Policy&)
0301 {
0302    typedef typename tools::promote_args<T>::type result_type;
0303    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0304    typedef typename policies::precision<result_type, Policy>::type precision_type;
0305    typedef typename policies::normalise<
0306       Policy,
0307       policies::promote_float<false>,
0308       policies::promote_double<false>,
0309       policies::discrete_quantile<>,
0310       policies::assert_undefined<> >::type forwarding_policy;
0311
0312    typedef boost::math::integral_constant<int,
0313       precision_type::value <= 0 ? 0 :
0314       precision_type::value <= 53 ? 53 :
0315       precision_type::value <= 64 ? 64 : 0
0316    > tag_type;
0317
0318    detail::log1p_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
0319
0320    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0321       detail::log1p_imp(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)");
0322 }
0323
0324 #ifdef log1p
0325 #  ifndef BOOST_HAS_LOG1P
0326 #     define BOOST_HAS_LOG1P
0327 #  endif
0328 #  undef log1p
0329 #endif
0330
0331 #if defined(BOOST_HAS_LOG1P) && !(defined(__osf__) && defined(__DECCXX_VER))
0332 #  ifdef BOOST_MATH_USE_C99
0333 template <class Policy>
0334 BOOST_MATH_GPU_ENABLED inline float log1p(float x, const Policy& pol)
0335 {
0336    if(x < -1)
0337       return policies::raise_domain_error<float>(
0338          "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0339    if(x == -1)
0340       return -policies::raise_overflow_error<float>(
0341          "log1p<%1%>(%1%)", nullptr, pol);
0342    return ::log1pf(x);
0343 }
0344 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0345 template <class Policy>
0346 BOOST_MATH_GPU_ENABLED inline long double log1p(long double x, const Policy& pol)
0347 {
0348    if(x < -1)
0349       return policies::raise_domain_error<long double>(
0350          "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0351    if(x == -1)
0352       return -policies::raise_overflow_error<long double>(
0353          "log1p<%1%>(%1%)", nullptr, pol);
0354    return ::log1pl(x);
0355 }
0356 #endif
0357 #else
0358 template <class Policy>
0359 inline float log1p(float x, const Policy& pol)
0360 {
0361    if(x < -1)
0362       return policies::raise_domain_error<float>(
0363          "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0364    if(x == -1)
0365       return -policies::raise_overflow_error<float>(
0366          "log1p<%1%>(%1%)", nullptr, pol);
0367    return ::log1p(x);
0368 }
0369 #endif
0370 template <class Policy>
0371 BOOST_MATH_GPU_ENABLED inline double log1p(double x, const Policy& pol)
0372 {
0373    if(x < -1)
0374       return policies::raise_domain_error<double>(
0375          "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0376    if(x == -1)
0377       return -policies::raise_overflow_error<double>(
0378          "log1p<%1%>(%1%)", nullptr, pol);
0379    return ::log1p(x);
0380 }
0381 #elif defined(_MSC_VER) && (BOOST_MSVC >= 1400)
0382 //
0383 // You should only enable this branch if you are absolutely sure
0384 // that your compilers optimizer won't mess this code up!!
0385 // Currently tested with VC8 and Intel 9.1.
0386 //
0387 template <class Policy>
0388 inline double log1p(double x, const Policy& pol)
0389 {
0390    if(x < -1)
0391       return policies::raise_domain_error<double>(
0392          "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0393    if(x == -1)
0394       return -policies::raise_overflow_error<double>(
0395          "log1p<%1%>(%1%)", nullptr, pol);
0396    double u = 1+x;
0397    if(u == 1.0)
0398       return x;
0399    else
0400       return ::log(u)*(x/(u-1.0));
0401 }
0402 template <class Policy>
0403 inline float log1p(float x, const Policy& pol)
0404 {
0405    return static_cast<float>(boost::math::log1p(static_cast<double>(x), pol));
0406 }
0407 #ifndef _WIN32_WCE
0408 //
0409 // For some reason this fails to compile under WinCE...
0410 // Needs more investigation.
0411 //
0412 template <class Policy>
0413 inline long double log1p(long double x, const Policy& pol)
0414 {
0415    if(x < -1)
0416       return policies::raise_domain_error<long double>(
0417          "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
0418    if(x == -1)
0419       return -policies::raise_overflow_error<long double>(
0420          "log1p<%1%>(%1%)", nullptr, pol);
0421    long double u = 1+x;
0422    if(u == 1.0)
0423       return x;
0424    else
0425       return ::logl(u)*(x/(u-1.0));
0426 }
0427 #endif
0428 #endif
0429
0430 template <class T>
0431 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1p(T x)
0432 {
0433    return boost::math::log1p(x, policies::policy<>());
0434 }
0435 //
0436 // Compute log(1+x)-x:
0437 //
0438 template <class T, class Policy>
0439 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
0440    log1pmx(T x, const Policy& pol)
0441 {
0442    typedef typename tools::promote_args<T>::type result_type;
0443    BOOST_MATH_STD_USING
0444    constexpr auto function = "boost::math::log1pmx<%1%>(%1%)";
0445
0446    if(x < -1)
0447       return policies::raise_domain_error<T>(
0448          function, "log1pmx(x) requires x > -1, but got x = %1%.", x, pol);
0449    if(x == -1)
0450       return -policies::raise_overflow_error<T>(
0451          function, nullptr, pol);
0452
0453    result_type a = abs(result_type(x));
0454    if(a > result_type(0.95f))
0455       return log(1 + result_type(x)) - result_type(x);
0456    // Note that without numeric_limits specialisation support,
0457    // epsilon just returns zero, and our "optimisation" will always fail:
0458    if(a < tools::epsilon<result_type>())
0459       return -x * x / 2;
0460    boost::math::detail::log1p_series<T> s(x);
0461    s();
0462    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0463
0464    T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
0465
0466    policies::check_series_iterations<T>(function, max_iter, pol);
0467    return result;
0468 }
0469
0470 template <class T>
0471 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1pmx(T x)
0472 {
0473    return log1pmx(x, policies::policy<>());
0474 }
0475
0476 } // namespace math
0477 } // namespace boost
0478
0479 #ifdef _MSC_VER
0480 #pragma warning(pop)
0481 #endif
0482
0483 #endif // BOOST_MATH_LOG1P_INCLUDED
0484
0485
0486