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