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