detail/functions/pow.hpp

0001
0002 // Copyright Christopher Kormanyos 2002 - 2013.
0003 // Copyright 2011 - 2013 John Maddock.
0004 // Distributed under the Boost Software License, Version 1.0.
0005 //    (See accompanying file LICENSE_1_0.txt or copy at
0006 //          http://www.boost.org/LICENSE_1_0.txt)
0007
0008 // This work is based on an earlier work:
0009 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
0010 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
0011 //
0012 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
0013 //
0014
0015 #ifdef BOOST_MSVC
0016 #pragma warning(push)
0017 #pragma warning(disable : 6326) // comparison of two constants
0018 #pragma warning(disable : 4127) // conditional expression is constant
0019 #endif
0020
0021 #include <boost/multiprecision/detail/standalone_config.hpp>
0022 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0023 #include <boost/multiprecision/detail/assert.hpp>
0024
0025 namespace detail {
0026
0027 template <typename T, typename U>
0028 inline void pow_imp(T& result, const T& t, const U& p, const std::integral_constant<bool, false>&)
0029 {
0030    // Compute the pure power of typename T t^p.
0031    // Use the S-and-X binary method, as described in
0032    // D. E. Knuth, "The Art of Computer Programming", Vol. 2,
0033    // Section 4.6.3 . The resulting computational complexity
0034    // is order log2[abs(p)].
0035
0036    using int_type = typename boost::multiprecision::detail::canonical<U, T>::type;
0037
0038    if (&result == &t)
0039    {
0040       T temp;
0041       pow_imp(temp, t, p, std::integral_constant<bool, false>());
0042       result = temp;
0043       return;
0044    }
0045
0046    // This will store the result.
0047    if (U(p % U(2)) != U(0))
0048    {
0049       result = t;
0050    }
0051    else
0052       result = int_type(1);
0053
0054    U p2(p);
0055
0056    // The variable x stores the binary powers of t.
0057    T x(t);
0058
0059    while (U(p2 /= 2) != U(0))
0060    {
0061       // Square x for each binary power.
0062       eval_multiply(x, x);
0063
0064       const bool has_binary_power = (U(p2 % U(2)) != U(0));
0065
0066       if (has_binary_power)
0067       {
0068          // Multiply the result with each binary power contained in the exponent.
0069          eval_multiply(result, x);
0070       }
0071    }
0072 }
0073
0074 template <typename T, typename U>
0075 inline void pow_imp(T& result, const T& t, const U& p, const std::integral_constant<bool, true>&)
0076 {
0077    // Signed integer power, just take care of the sign then call the unsigned version:
0078    using int_type = typename boost::multiprecision::detail::canonical<U, T>::type;
0079    using ui_type = typename boost::multiprecision::detail::make_unsigned<U>::type                         ;
0080
0081    if (p < 0)
0082    {
0083       T temp;
0084       temp = static_cast<int_type>(1);
0085       T denom;
0086       pow_imp(denom, t, static_cast<ui_type>(-p), std::integral_constant<bool, false>());
0087       eval_divide(result, temp, denom);
0088       return;
0089    }
0090    pow_imp(result, t, static_cast<ui_type>(p), std::integral_constant<bool, false>());
0091 }
0092
0093 } // namespace detail
0094
0095 template <typename T, typename U>
0096 inline typename std::enable_if<boost::multiprecision::detail::is_integral<U>::value>::type eval_pow(T& result, const T& t, const U& p)
0097 {
0098    detail::pow_imp(result, t, p, boost::multiprecision::detail::is_signed<U>());
0099 }
0100
0101 template <class T>
0102 void hyp0F0(T& H0F0, const T& x)
0103 {
0104    // Compute the series representation of Hypergeometric0F0 taken from
0105    // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F0/06/01/
0106    // There are no checks on input range or parameter boundaries.
0107
0108    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
0109
0110    BOOST_MP_ASSERT(&H0F0 != &x);
0111    long tol = boost::multiprecision::detail::digits2<number<T, et_on> >::value();
0112    T    t;
0113
0114    T x_pow_n_div_n_fact(x);
0115
0116    eval_add(H0F0, x_pow_n_div_n_fact, ui_type(1));
0117
0118    T lim;
0119    eval_ldexp(lim, H0F0, static_cast<int>(1L - tol));
0120    if (eval_get_sign(lim) < 0)
0121       lim.negate();
0122
0123    ui_type n;
0124
0125    const unsigned series_limit =
0126        boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
0127            ? 100
0128            : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
0129    // Series expansion of hyperg_0f0(; ; x).
0130    for (n = 2; n < series_limit; ++n)
0131    {
0132       eval_multiply(x_pow_n_div_n_fact, x);
0133       eval_divide(x_pow_n_div_n_fact, n);
0134       eval_add(H0F0, x_pow_n_div_n_fact);
0135       bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0;
0136       if (neg)
0137          x_pow_n_div_n_fact.negate();
0138       if (lim.compare(x_pow_n_div_n_fact) > 0)
0139          break;
0140       if (neg)
0141          x_pow_n_div_n_fact.negate();
0142    }
0143    if (n >= series_limit)
0144       BOOST_MP_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge"));
0145 }
0146
0147 template <class T>
0148 void hyp1F0(T& H1F0, const T& a, const T& x)
0149 {
0150    // Compute the series representation of Hypergeometric1F0 taken from
0151    // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F0/06/01/01/
0152    // and also see the corresponding section for the power function (i.e. x^a).
0153    // There are no checks on input range or parameter boundaries.
0154
0155    using si_type = typename boost::multiprecision::detail::canonical<int, T>::type;
0156
0157    BOOST_MP_ASSERT(&H1F0 != &x);
0158    BOOST_MP_ASSERT(&H1F0 != &a);
0159
0160    T x_pow_n_div_n_fact(x);
0161    T pochham_a(a);
0162    T ap(a);
0163
0164    eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact);
0165    eval_add(H1F0, si_type(1));
0166    T lim;
0167    eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0168    if (eval_get_sign(lim) < 0)
0169       lim.negate();
0170
0171    si_type n;
0172    T       term, part;
0173
0174    const si_type series_limit =
0175        boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
0176            ? 100
0177            : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
0178    // Series expansion of hyperg_1f0(a; ; x).
0179    for (n = 2; n < series_limit; n++)
0180    {
0181       eval_multiply(x_pow_n_div_n_fact, x);
0182       eval_divide(x_pow_n_div_n_fact, n);
0183       eval_increment(ap);
0184       eval_multiply(pochham_a, ap);
0185       eval_multiply(term, pochham_a, x_pow_n_div_n_fact);
0186       eval_add(H1F0, term);
0187       if (eval_get_sign(term) < 0)
0188          term.negate();
0189       if (lim.compare(term) >= 0)
0190          break;
0191    }
0192    if (n >= series_limit)
0193       BOOST_MP_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge"));
0194 }
0195
0196 template <class T>
0197 void eval_exp(T& result, const T& x)
0198 {
0199    static_assert(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types.");
0200    if (&x == &result)
0201    {
0202       T temp;
0203       eval_exp(temp, x);
0204       result = temp;
0205       return;
0206    }
0207    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
0208    using si_type = typename boost::multiprecision::detail::canonical<int, T>::type     ;
0209    using exp_type = typename T::exponent_type                                           ;
0210    using canonical_exp_type = typename boost::multiprecision::detail::canonical<exp_type, T>::type;
0211
0212    // Handle special arguments.
0213    int  type  = eval_fpclassify(x);
0214    bool isneg = eval_get_sign(x) < 0;
0215    if (type == static_cast<int>(FP_NAN))
0216    {
0217       result = x;
0218       errno  = EDOM;
0219       return;
0220    }
0221    else if (type == static_cast<int>(FP_INFINITE))
0222    {
0223       if (isneg)
0224          result = ui_type(0u);
0225       else
0226          result = x;
0227       return;
0228    }
0229    else if (type == static_cast<int>(FP_ZERO))
0230    {
0231       result = ui_type(1);
0232       return;
0233    }
0234
0235    // Get local copy of argument and force it to be positive.
0236    T xx = x;
0237    T exp_series;
0238    if (isneg)
0239       xx.negate();
0240
0241    // Check the range of the argument.
0242    if (xx.compare(si_type(1)) <= 0)
0243    {
0244       //
0245       // Use series for exp(x) - 1:
0246       //
0247       T lim;
0248       BOOST_IF_CONSTEXPR(std::numeric_limits<number<T, et_on> >::is_specialized)
0249          lim = std::numeric_limits<number<T, et_on> >::epsilon().backend();
0250       else
0251       {
0252          result = ui_type(1);
0253          eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0254       }
0255       unsigned k = 2;
0256       exp_series = xx;
0257       result     = si_type(1);
0258       if (isneg)
0259          eval_subtract(result, exp_series);
0260       else
0261          eval_add(result, exp_series);
0262       eval_multiply(exp_series, xx);
0263       eval_divide(exp_series, ui_type(k));
0264       eval_add(result, exp_series);
0265       while (exp_series.compare(lim) > 0)
0266       {
0267          ++k;
0268          eval_multiply(exp_series, xx);
0269          eval_divide(exp_series, ui_type(k));
0270          if (isneg && (k & 1))
0271             eval_subtract(result, exp_series);
0272          else
0273             eval_add(result, exp_series);
0274       }
0275       return;
0276    }
0277
0278    // Check for pure-integer arguments which can be either signed or unsigned.
0279    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type ll;
0280    eval_trunc(exp_series, x);
0281    eval_convert_to(&ll, exp_series);
0282    if (x.compare(ll) == 0)
0283    {
0284       detail::pow_imp(result, get_constant_e<T>(), ll, std::integral_constant<bool, true>());
0285       return;
0286    }
0287    else if (exp_series.compare(x) == 0)
0288    {
0289       // We have a value that has no fractional part, but is too large to fit
0290       // in a long long, in this situation the code below will fail, so
0291       // we're just going to assume that this will overflow:
0292       if (isneg)
0293          result = ui_type(0);
0294       else
0295          result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
0296       return;
0297    }
0298
0299    // The algorithm for exp has been taken from MPFUN.
0300    // exp(t) = [ (1 + r + r^2/2! + r^3/3! + r^4/4! ...)^p2 ] * 2^n
0301    // where p2 is a power of 2 such as 2048, r = t_prime / p2, and
0302    // t_prime = t - n*ln2, with n chosen to minimize the absolute
0303    // value of t_prime. In the resulting Taylor series, which is
0304    // implemented as a hypergeometric function, |r| is bounded by
0305    // ln2 / p2. For small arguments, no scaling is done.
0306
0307    // Compute the exponential series of the (possibly) scaled argument.
0308
0309    eval_divide(result, xx, get_constant_ln2<T>());
0310    exp_type n;
0311    eval_convert_to(&n, result);
0312
0313    if (n == (std::numeric_limits<exp_type>::max)())
0314    {
0315       // Exponent is too large to fit in our exponent type:
0316       if (isneg)
0317          result = ui_type(0);
0318       else
0319          result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
0320       return;
0321    }
0322
0323    // The scaling is 2^11 = 2048.
0324    const si_type p2 = static_cast<si_type>(si_type(1) << 11);
0325
0326    eval_multiply(exp_series, get_constant_ln2<T>(), static_cast<canonical_exp_type>(n));
0327    eval_subtract(exp_series, xx);
0328    eval_divide(exp_series, p2);
0329    exp_series.negate();
0330    hyp0F0(result, exp_series);
0331
0332    detail::pow_imp(exp_series, result, p2, std::integral_constant<bool, true>());
0333    result = ui_type(1);
0334    eval_ldexp(result, result, n);
0335    eval_multiply(exp_series, result);
0336
0337    if (isneg)
0338       eval_divide(result, ui_type(1), exp_series);
0339    else
0340       result = exp_series;
0341 }
0342
0343 template <class T>
0344 void eval_log(T& result, const T& arg)
0345 {
0346    static_assert(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
0347    //
0348    // We use a variation of http://dlmf.nist.gov/4.45#i
0349    // using frexp to reduce the argument to x * 2^n,
0350    // then let y = x - 1 and compute:
0351    // log(x) = log(2) * n + log1p(1 + y)
0352    //
0353    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
0354    using exp_type = typename T::exponent_type                                           ;
0355    using canonical_exp_type = typename boost::multiprecision::detail::canonical<exp_type, T>::type;
0356    using fp_type = typename std::tuple_element<0, typename T::float_types>::type                  ;
0357    int                                                                          s = eval_signbit(arg);
0358    switch (eval_fpclassify(arg))
0359    {
0360    case FP_NAN:
0361       result = arg;
0362       errno  = EDOM;
0363       return;
0364    case FP_INFINITE:
0365       if (s)
0366          break;
0367       result = arg;
0368       return;
0369    case FP_ZERO:
0370       result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
0371       result.negate();
0372       errno = ERANGE;
0373       return;
0374    }
0375    if (s)
0376    {
0377       result = std::numeric_limits<number<T> >::quiet_NaN().backend();
0378       errno  = EDOM;
0379       return;
0380    }
0381
0382    exp_type e;
0383    T        t;
0384    eval_frexp(t, arg, &e);
0385    bool alternate = false;
0386
0387    if (t.compare(fp_type(2) / fp_type(3)) <= 0)
0388    {
0389       alternate = true;
0390       eval_ldexp(t, t, 1);
0391       --e;
0392    }
0393
0394    eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e));
0395    INSTRUMENT_BACKEND(result);
0396    eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */
0397    if (!alternate)
0398       t.negate(); /* 0 <= t <= 0.33333 */
0399    T pow = t;
0400    T lim;
0401    T t2;
0402
0403    if (alternate)
0404       eval_add(result, t);
0405    else
0406       eval_subtract(result, t);
0407
0408    BOOST_IF_CONSTEXPR(std::numeric_limits<number<T, et_on> >::is_specialized)
0409       eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend());
0410    else
0411       eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0412    if (eval_get_sign(lim) < 0)
0413       lim.negate();
0414    INSTRUMENT_BACKEND(lim);
0415
0416    ui_type k = 1;
0417    do
0418    {
0419       ++k;
0420       eval_multiply(pow, t);
0421       eval_divide(t2, pow, k);
0422       INSTRUMENT_BACKEND(t2);
0423       if (alternate && ((k & 1) != 0))
0424          eval_add(result, t2);
0425       else
0426          eval_subtract(result, t2);
0427       INSTRUMENT_BACKEND(result);
0428    } while (lim.compare(t2) < 0);
0429 }
0430
0431 template <class T>
0432 const T& get_constant_log10()
0433 {
0434    static BOOST_MP_THREAD_LOCAL T             result;
0435    static BOOST_MP_THREAD_LOCAL long digits = 0;
0436    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
0437    {
0438       using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
0439       T                                                                            ten;
0440       ten = ui_type(10u);
0441       eval_log(result, ten);
0442       digits = boost::multiprecision::detail::digits2<number<T> >::value();
0443    }
0444
0445    return result;
0446 }
0447
0448 template <class T>
0449 void eval_log10(T& result, const T& arg)
0450 {
0451    static_assert(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types.");
0452    eval_log(result, arg);
0453    eval_divide(result, get_constant_log10<T>());
0454 }
0455
0456 template <class R, class T>
0457 inline void eval_log2(R& result, const T& a)
0458 {
0459    eval_log(result, a);
0460    eval_divide(result, get_constant_ln2<R>());
0461 }
0462
0463 template <typename T>
0464 inline void eval_pow(T& result, const T& x, const T& a)
0465 {
0466    static_assert(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types.");
0467    using si_type = typename boost::multiprecision::detail::canonical<int, T>::type;
0468    using fp_type = typename std::tuple_element<0, typename T::float_types>::type             ;
0469
0470    if ((&result == &x) || (&result == &a))
0471    {
0472       T t;
0473       eval_pow(t, x, a);
0474       result = t;
0475       return;
0476    }
0477
0478    if ((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
0479    {
0480       result = x;
0481       return;
0482    }
0483    if (a.compare(si_type(0)) == 0)
0484    {
0485       result = si_type(1);
0486       return;
0487    }
0488
0489    int type = eval_fpclassify(x);
0490
0491    switch (type)
0492    {
0493    case FP_ZERO:
0494       switch (eval_fpclassify(a))
0495       {
0496       case FP_ZERO:
0497          result = si_type(1);
0498          break;
0499       case FP_NAN:
0500          result = a;
0501          break;
0502       case FP_NORMAL: {
0503          // Need to check for a an odd integer as a special case:
0504          BOOST_MP_TRY
0505          {
0506             typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type i;
0507             eval_convert_to(&i, a);
0508             if (a.compare(i) == 0)
0509             {
0510                if (eval_signbit(a))
0511                {
0512                   if (i & 1)
0513                   {
0514                      result = std::numeric_limits<number<T> >::infinity().backend();
0515                      if (eval_signbit(x))
0516                         result.negate();
0517                      errno = ERANGE;
0518                   }
0519                   else
0520                   {
0521                      result = std::numeric_limits<number<T> >::infinity().backend();
0522                      errno  = ERANGE;
0523                   }
0524                }
0525                else if (i & 1)
0526                {
0527                   result = x;
0528                }
0529                else
0530                   result = si_type(0);
0531                return;
0532             }
0533          }
0534          BOOST_MP_CATCH(const std::exception&)
0535          {
0536             // fallthrough..
0537          }
0538          BOOST_MP_CATCH_END
0539          BOOST_FALLTHROUGH;
0540       }
0541       default:
0542          if (eval_signbit(a))
0543          {
0544             result = std::numeric_limits<number<T> >::infinity().backend();
0545             errno  = ERANGE;
0546          }
0547          else
0548             result = x;
0549          break;
0550       }
0551       return;
0552    case FP_NAN:
0553       result = x;
0554       errno  = ERANGE;
0555       return;
0556    default:;
0557    }
0558
0559    int s = eval_get_sign(a);
0560    if (s == 0)
0561    {
0562       result = si_type(1);
0563       return;
0564    }
0565
0566    if (s < 0)
0567    {
0568       T t, da;
0569       t = a;
0570       t.negate();
0571       eval_pow(da, x, t);
0572       eval_divide(result, si_type(1), da);
0573       return;
0574    }
0575
0576    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type an;
0577    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type max_an =
0578        std::numeric_limits<typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type>::max)() : static_cast<typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type>(1) << (sizeof(typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type) * CHAR_BIT - 2);
0579    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type min_an =
0580        std::numeric_limits<typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type>::min)() : -min_an;
0581
0582    T fa;
0583    BOOST_MP_TRY
0584    {
0585       eval_convert_to(&an, a);
0586       if (a.compare(an) == 0)
0587       {
0588          detail::pow_imp(result, x, an, std::integral_constant<bool, true>());
0589          return;
0590       }
0591    }
0592    BOOST_MP_CATCH(const std::exception&)
0593    {
0594       // conversion failed, just fall through, value is not an integer.
0595       an = (std::numeric_limits<std::intmax_t>::max)();
0596    }
0597    BOOST_MP_CATCH_END
0598    if ((eval_get_sign(x) < 0))
0599    {
0600       typename boost::multiprecision::detail::canonical<std::uintmax_t, T>::type aun;
0601       BOOST_MP_TRY
0602       {
0603          eval_convert_to(&aun, a);
0604          if (a.compare(aun) == 0)
0605          {
0606             fa = x;
0607             fa.negate();
0608             eval_pow(result, fa, a);
0609             if (aun & 1u)
0610                result.negate();
0611             return;
0612          }
0613       }
0614       BOOST_MP_CATCH(const std::exception&)
0615       {
0616          // conversion failed, just fall through, value is not an integer.
0617       }
0618       BOOST_MP_CATCH_END
0619
0620       eval_floor(result, a);
0621       // -1^INF is a special case in C99:
0622       if ((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
0623       {
0624          result = si_type(1);
0625       }
0626       else if (a.compare(result) == 0)
0627       {
0628          // exponent is so large we have no fractional part:
0629          if (x.compare(si_type(-1)) < 0)
0630          {
0631             result = std::numeric_limits<number<T, et_on> >::infinity().backend();
0632          }
0633          else
0634          {
0635             result = si_type(0);
0636          }
0637       }
0638       else if (type == FP_INFINITE)
0639       {
0640          result = std::numeric_limits<number<T, et_on> >::infinity().backend();
0641       }
0642       else BOOST_IF_CONSTEXPR (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
0643       {
0644          result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
0645          errno  = EDOM;
0646       }
0647       else
0648       {
0649          BOOST_MP_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type."));
0650       }
0651       return;
0652    }
0653
0654    T t, da;
0655
0656    eval_subtract(da, a, an);
0657
0658    if ((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
0659    {
0660       if (a.compare(fp_type(1e-5f)) <= 0)
0661       {
0662          // Series expansion for small a.
0663          eval_log(t, x);
0664          eval_multiply(t, a);
0665          hyp0F0(result, t);
0666          return;
0667       }
0668       else
0669       {
0670          // Series expansion for moderately sized x. Note that for large power of a,
0671          // the power of the integer part of a is calculated using the pown function.
0672          if (an)
0673          {
0674             da.negate();
0675             t = si_type(1);
0676             eval_subtract(t, x);
0677             hyp1F0(result, da, t);
0678             detail::pow_imp(t, x, an, std::integral_constant<bool, true>());
0679             eval_multiply(result, t);
0680          }
0681          else
0682          {
0683             da = a;
0684             da.negate();
0685             t = si_type(1);
0686             eval_subtract(t, x);
0687             hyp1F0(result, da, t);
0688          }
0689       }
0690    }
0691    else
0692    {
0693       // Series expansion for pow(x, a). Note that for large power of a, the power
0694       // of the integer part of a is calculated using the pown function.
0695       if (an)
0696       {
0697          eval_log(t, x);
0698          eval_multiply(t, da);
0699          eval_exp(result, t);
0700          detail::pow_imp(t, x, an, std::integral_constant<bool, true>());
0701          eval_multiply(result, t);
0702       }
0703       else
0704       {
0705          eval_log(t, x);
0706          eval_multiply(t, a);
0707          eval_exp(result, t);
0708       }
0709    }
0710 }
0711
0712 template <class T, class A>
0713 #if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
0714 inline typename std::enable_if<!boost::multiprecision::detail::is_integral<A>::value, void>::type
0715 #else
0716 inline typename std::enable_if<is_compatible_arithmetic_type<A, number<T> >::value && !boost::multiprecision::detail::is_integral<A>::value, void>::type
0717 #endif
0718 eval_pow(T& result, const T& x, const A& a)
0719 {
0720    // Note this one is restricted to float arguments since pow.hpp already has a version for
0721    // integer powers....
0722    using canonical_type = typename boost::multiprecision::detail::canonical<A, T>::type         ;
0723    using cast_type = typename std::conditional<std::is_same<A, canonical_type>::value, T, canonical_type>::type;
0724    cast_type                                                                      c;
0725    c = a;
0726    eval_pow(result, x, c);
0727 }
0728
0729 template <class T, class A>
0730 #if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
0731 inline void
0732 #else
0733 inline typename std::enable_if<is_compatible_arithmetic_type<A, number<T> >::value, void>::type
0734 #endif
0735 eval_pow(T& result, const A& x, const T& a)
0736 {
0737    using canonical_type = typename boost::multiprecision::detail::canonical<A, T>::type         ;
0738    using cast_type = typename std::conditional<std::is_same<A, canonical_type>::value, T, canonical_type>::type;
0739    cast_type                                                                      c;
0740    c = x;
0741    eval_pow(result, c, a);
0742 }
0743
0744 template <class T>
0745 void eval_exp2(T& result, const T& arg)
0746 {
0747    static_assert(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
0748
0749    // Check for pure-integer arguments which can be either signed or unsigned.
0750    typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
0751    T                                                                                     temp;
0752    BOOST_MP_TRY
0753    {
0754       eval_trunc(temp, arg);
0755       eval_convert_to(&i, temp);
0756       if (arg.compare(i) == 0)
0757       {
0758          temp = static_cast<typename std::tuple_element<0, typename T::unsigned_types>::type>(1u);
0759          eval_ldexp(result, temp, i);
0760          return;
0761       }
0762    }
0763    #ifdef BOOST_MP_MATH_AVAILABLE
0764    BOOST_MP_CATCH(const boost::math::rounding_error&)
0765    { /* Fallthrough */
0766    }
0767    #endif
0768    BOOST_MP_CATCH(const std::runtime_error&)
0769    { /* Fallthrough */
0770    }
0771    BOOST_MP_CATCH_END
0772
0773    temp = static_cast<typename std::tuple_element<0, typename T::unsigned_types>::type>(2u);
0774    eval_pow(result, temp, arg);
0775 }
0776
0777 namespace detail {
0778
0779 template <class T>
0780 void small_sinh_series(T x, T& result)
0781 {
0782    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
0783    bool                                                                         neg = eval_get_sign(x) < 0;
0784    if (neg)
0785       x.negate();
0786    T p(x);
0787    T mult(x);
0788    eval_multiply(mult, x);
0789    result    = x;
0790    ui_type k = 1;
0791
0792    T lim(x);
0793    eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0794
0795    do
0796    {
0797       eval_multiply(p, mult);
0798       eval_divide(p, ++k);
0799       eval_divide(p, ++k);
0800       eval_add(result, p);
0801    } while (p.compare(lim) >= 0);
0802    if (neg)
0803       result.negate();
0804 }
0805
0806 template <class T>
0807 void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
0808 {
0809    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
0810    using fp_type = typename std::tuple_element<0, typename T::float_types>::type                  ;
0811
0812    switch (eval_fpclassify(x))
0813    {
0814    case FP_NAN:
0815       errno = EDOM;
0816       // fallthrough...
0817    case FP_INFINITE:
0818       if (p_sinh)
0819          *p_sinh = x;
0820       if (p_cosh)
0821       {
0822          *p_cosh = x;
0823          if (eval_get_sign(x) < 0)
0824             p_cosh->negate();
0825       }
0826       return;
0827    case FP_ZERO:
0828       if (p_sinh)
0829          *p_sinh = x;
0830       if (p_cosh)
0831          *p_cosh = ui_type(1);
0832       return;
0833    default:;
0834    }
0835
0836    bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;
0837
0838    if (p_cosh || !small_sinh)
0839    {
0840       T e_px, e_mx;
0841       eval_exp(e_px, x);
0842       eval_divide(e_mx, ui_type(1), e_px);
0843       if (eval_signbit(e_mx) != eval_signbit(e_px))
0844          e_mx.negate(); // Handles lack of signed zero in some types
0845
0846       if (p_sinh)
0847       {
0848          if (small_sinh)
0849          {
0850             small_sinh_series(x, *p_sinh);
0851          }
0852          else
0853          {
0854             eval_subtract(*p_sinh, e_px, e_mx);
0855             eval_ldexp(*p_sinh, *p_sinh, -1);
0856          }
0857       }
0858       if (p_cosh)
0859       {
0860          eval_add(*p_cosh, e_px, e_mx);
0861          eval_ldexp(*p_cosh, *p_cosh, -1);
0862       }
0863    }
0864    else
0865    {
0866       small_sinh_series(x, *p_sinh);
0867    }
0868 }
0869
0870 } // namespace detail
0871
0872 template <class T>
0873 inline void eval_sinh(T& result, const T& x)
0874 {
0875    static_assert(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types.");
0876    detail::sinhcosh(x, &result, static_cast<T*>(0));
0877 }
0878
0879 template <class T>
0880 inline void eval_cosh(T& result, const T& x)
0881 {
0882    static_assert(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types.");
0883    detail::sinhcosh(x, static_cast<T*>(0), &result);
0884 }
0885
0886 template <class T>
0887 inline void eval_tanh(T& result, const T& x)
0888 {
0889    static_assert(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
0890    T c;
0891    detail::sinhcosh(x, &result, &c);
0892    if ((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE))
0893    {
0894       bool s = eval_signbit(result) != eval_signbit(c);
0895       result = static_cast<typename std::tuple_element<0, typename T::unsigned_types>::type>(1u);
0896       if (s)
0897          result.negate();
0898       return;
0899    }
0900    eval_divide(result, c);
0901 }
0902
0903 #ifdef BOOST_MSVC
0904 #pragma warning(pop)
0905 #endif