ceph/src/boost/boost/math/special_functions/expm1.hpp

   1 //  (C) Copyright John Maddock 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_EXPM1_INCLUDED
   7 #define BOOST_MATH_EXPM1_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 ::expm1
  15 #include <boost/limits.hpp>
  16 #include <boost/math/tools/config.hpp>
  17 #include <boost/math/tools/series.hpp>
  18 #include <boost/math/tools/precision.hpp>
  19 #include <boost/math/tools/big_constant.hpp>
  20 #include <boost/math/policies/error_handling.hpp>
  21 #include <boost/math/tools/rational.hpp>
  22 #include <boost/math/special_functions/math_fwd.hpp>
  23 #include <boost/mpl/less_equal.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 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
  32 //
  33 // This is the only way we can avoid
  34 // warning: non-standard suffix on floating constant [-Wpedantic]
  35 // when building with -Wall -pedantic.  Neither __extension__
  36 // nor #pragma diagnostic ignored work :(
  37 //
  38 #pragma GCC system_header
  39 #endif
  40
  41 namespace boost{ namespace math{
  42
  43 namespace detail
  44 {
  45   // Functor expm1_series returns the next term in the Taylor series
  46   // x^k / k!
  47   // each time that operator() is invoked.
  48   //
  49   template <class T>
  50   struct expm1_series
  51   {
  52      typedef T result_type;
  53
  54      expm1_series(T x)
  55         : k(0), m_x(x), m_term(1) {}
  56
  57      T operator()()
  58      {
  59         ++k;
  60         m_term *= m_x;
  61         m_term /= k;
  62         return m_term;
  63      }
  64
  65      int count()const
  66      {
  67         return k;
  68      }
  69
  70   private:
  71      int k;
  72      const T m_x;
  73      T m_term;
  74      expm1_series(const expm1_series&);
  75      expm1_series& operator=(const expm1_series&);
  76   };
  77
  78 template <class T, class Policy, class tag>
  79 struct expm1_initializer
  80 {
  81    struct init
  82    {
  83       init()
  84       {
  85          do_init(tag());
  86       }
  87       template <int N>
  88       static void do_init(const boost::integral_constant<int, N>&){}
  89       static void do_init(const boost::integral_constant<int, 64>&)
  90       {
  91          expm1(T(0.5));
  92       }
  93       static void do_init(const boost::integral_constant<int, 113>&)
  94       {
  95          expm1(T(0.5));
  96       }
  97       void force_instantiate()const{}
  98    };
  99    static const init initializer;
 100    static void force_instantiate()
 101    {
 102       initializer.force_instantiate();
 103    }
 104 };
 105
 106 template <class T, class Policy, class tag>
 107 const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer;
 108
 109 //
 110 // Algorithm expm1 is part of C99, but is not yet provided by many compilers.
 111 //
 112 // This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
 113 //
 114 template <class T, class Policy>
 115 T expm1_imp(T x, const boost::integral_constant<int, 0>&, const Policy& pol)
 116 {
 117    BOOST_MATH_STD_USING
 118
 119    T a = fabs(x);
 120    if((boost::math::isnan)(a))
 121    {
 122       return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
 123    }
 124    if(a > T(0.5f))
 125    {
 126       if(a >= tools::log_max_value<T>())
 127       {
 128          if(x > 0)
 129             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
 130          return -1;
 131       }
 132       return exp(x) - T(1);
 133    }
 134    if(a < tools::epsilon<T>())
 135       return x;
 136    detail::expm1_series<T> s(x);
 137    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 138 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
 139    T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
 140 #else
 141    T zero = 0;
 142    T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
 143 #endif
 144    policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
 145    return result;
 146 }
 147
 148 template <class T, class P>
 149 T expm1_imp(T x, const boost::integral_constant<int, 53>&, const P& pol)
 150 {
 151    BOOST_MATH_STD_USING
 152
 153    T a = fabs(x);
 154    if(a > T(0.5L))
 155    {
 156       if(a >= tools::log_max_value<T>())
 157       {
 158          if(x > 0)
 159             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
 160          return -1;
 161       }
 162       return exp(x) - T(1);
 163    }
 164    if(a < tools::epsilon<T>())
 165       return x;
 166
 167    static const float Y = 0.10281276702880859e1f;
 168    static const T n[] = { static_cast<T>(-0.28127670288085937e-1), static_cast<T>(0.51278186299064534e0), static_cast<T>(-0.6310029069350198e-1), static_cast<T>(0.11638457975729296e-1), static_cast<T>(-0.52143390687521003e-3), static_cast<T>(0.21491399776965688e-4) };
 169    static const T d[] = { 1, static_cast<T>(-0.45442309511354755e0), static_cast<T>(0.90850389570911714e-1), static_cast<T>(-0.10088963629815502e-1), static_cast<T>(0.63003407478692265e-3), static_cast<T>(-0.17976570003654402e-4) };
 170
 171    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
 172    return result;
 173 }
 174
 175 template <class T, class P>
 176 T expm1_imp(T x, const boost::integral_constant<int, 64>&, const P& pol)
 177 {
 178    BOOST_MATH_STD_USING
 179
 180    T a = fabs(x);
 181    if(a > T(0.5L))
 182    {
 183       if(a >= tools::log_max_value<T>())
 184       {
 185          if(x > 0)
 186             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
 187          return -1;
 188       }
 189       return exp(x) - T(1);
 190    }
 191    if(a < tools::epsilon<T>())
 192       return x;
 193
 194    static const float Y = 0.10281276702880859375e1f;
 195    static const T n[] = {
 196       BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),
 197        BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),
 198        BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
 199        BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
 200        BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
 201        BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
 202        BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
 203    };
 204    static const T d[] = {
 205       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 206       BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
 207       BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
 208       BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
 209       BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
 210       BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
 211       BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
 212    };
 213
 214    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
 215    return result;
 216 }
 217
 218 template <class T, class P>
 219 T expm1_imp(T x, const boost::integral_constant<int, 113>&, const P& pol)
 220 {
 221    BOOST_MATH_STD_USING
 222
 223    T a = fabs(x);
 224    if(a > T(0.5L))
 225    {
 226       if(a >= tools::log_max_value<T>())
 227       {
 228          if(x > 0)
 229             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
 230          return -1;
 231       }
 232       return exp(x) - T(1);
 233    }
 234    if(a < tools::epsilon<T>())
 235       return x;
 236
 237    static const float Y = 0.10281276702880859375e1f;
 238    static const T n[] = {
 239       BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
 240       BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
 241       BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
 242       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
 243       BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
 244       BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
 245       BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
 246       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
 247       BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
 248       BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
 249    };
 250    static const T d[] = {
 251       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 252       BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
 253       BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
 254       BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
 255       BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
 256       BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
 257       BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
 258       BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
 259       BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
 260       BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
 261       BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
 262    };
 263
 264    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
 265    return result;
 266 }
 267
 268 } // namespace detail
 269
 270 template <class T, class Policy>
 271 inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
 272 {
 273    typedef typename tools::promote_args<T>::type result_type;
 274    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 275    typedef typename policies::precision<result_type, Policy>::type precision_type;
 276    typedef typename policies::normalise<
 277       Policy,
 278       policies::promote_float<false>,
 279       policies::promote_double<false>,
 280       policies::discrete_quantile<>,
 281       policies::assert_undefined<> >::type forwarding_policy;
 282
 283    typedef boost::integral_constant<int,
 284       precision_type::value <= 0 ? 0 :
 285       precision_type::value <= 53 ? 53 :
 286       precision_type::value <= 64 ? 64 :
 287       precision_type::value <= 113 ? 113 : 0
 288    > tag_type;
 289
 290    detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
 291
 292    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
 293       static_cast<value_type>(x),
 294       tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
 295 }
 296
 297 #ifdef expm1
 298 #  ifndef BOOST_HAS_expm1
 299 #     define BOOST_HAS_expm1
 300 #  endif
 301 #  undef expm1
 302 #endif
 303
 304 #if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER))
 305 #  ifdef BOOST_MATH_USE_C99
 306 inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x); }
 307 #     ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 308 inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x); }
 309 #     endif
 310 #  else
 311 inline float expm1(float x, const policies::policy<>&){ return static_cast<float>(::expm1(x)); }
 312 #  endif
 313 inline double expm1(double x, const policies::policy<>&){ return ::expm1(x); }
 314 #endif
 315
 316 template <class T>
 317 inline typename tools::promote_args<T>::type expm1(T x)
 318 {
 319    return expm1(x, policies::policy<>());
 320 }
 321
 322 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))
 323 inline float expm1(float z)
 324 {
 325    return expm1<float>(z);
 326 }
 327 inline double expm1(double z)
 328 {
 329    return expm1<double>(z);
 330 }
 331 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 332 inline long double expm1(long double z)
 333 {
 334    return expm1<long double>(z);
 335 }
 336 #endif
 337 #endif
 338
 339 } // namespace math
 340 } // namespace boost
 341
 342 #endif // BOOST_MATH_HYPOT_INCLUDED
 343
 344
 345
 346