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