boost/math/special_functions/log1p.hpp

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