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