ceph/src/boost/boost/multiprecision/detail/functions/pow.hpp

   1
   2 // Copyright Christopher Kormanyos 2002 - 2013.
   3 // Copyright 2011 - 2013 John Maddock.
   4 // Distributed under the Boost Software License, Version 1.0.
   5 //    (See accompanying file LICENSE_1_0.txt or copy at
   6 //          http://www.boost.org/LICENSE_1_0.txt)
   7
   8 // This work is based on an earlier work:
   9 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
  10 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
  11 //
  12 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
  13 //
  14
  15 #ifdef BOOST_MSVC
  16 #pragma warning(push)
  17 #pragma warning(disable : 6326) // comparison of two constants
  18 #pragma warning(disable : 4127) // conditional expression is constant
  19 #endif
  20
  21 #include <boost/multiprecision/detail/standalone_config.hpp>
  22 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
  23 #include <boost/multiprecision/detail/assert.hpp>
  24
  25 namespace detail {
  26
  27 template <typename T, typename U>
  28 inline void pow_imp(T& result, const T& t, const U& p, const std::integral_constant<bool, false>&)
  29 {
  30    // Compute the pure power of typename T t^p.
  31    // Use the S-and-X binary method, as described in
  32    // D. E. Knuth, "The Art of Computer Programming", Vol. 2,
  33    // Section 4.6.3 . The resulting computational complexity
  34    // is order log2[abs(p)].
  35
  36    using int_type = typename boost::multiprecision::detail::canonical<U, T>::type;
  37
  38    if (&result == &t)
  39    {
  40       T temp;
  41       pow_imp(temp, t, p, std::integral_constant<bool, false>());
  42       result = temp;
  43       return;
  44    }
  45
  46    // This will store the result.
  47    if (U(p % U(2)) != U(0))
  48    {
  49       result = t;
  50    }
  51    else
  52       result = int_type(1);
  53
  54    U p2(p);
  55
  56    // The variable x stores the binary powers of t.
  57    T x(t);
  58
  59    while (U(p2 /= 2) != U(0))
  60    {
  61       // Square x for each binary power.
  62       eval_multiply(x, x);
  63
  64       const bool has_binary_power = (U(p2 % U(2)) != U(0));
  65
  66       if (has_binary_power)
  67       {
  68          // Multiply the result with each binary power contained in the exponent.
  69          eval_multiply(result, x);
  70       }
  71    }
  72 }
  73
  74 template <typename T, typename U>
  75 inline void pow_imp(T& result, const T& t, const U& p, const std::integral_constant<bool, true>&)
  76 {
  77    // Signed integer power, just take care of the sign then call the unsigned version:
  78    using int_type = typename boost::multiprecision::detail::canonical<U, T>::type;
  79    using ui_type = typename boost::multiprecision::detail::make_unsigned<U>::type                         ;
  80
  81    if (p < 0)
  82    {
  83       T temp;
  84       temp = static_cast<int_type>(1);
  85       T denom;
  86       pow_imp(denom, t, static_cast<ui_type>(-p), std::integral_constant<bool, false>());
  87       eval_divide(result, temp, denom);
  88       return;
  89    }
  90    pow_imp(result, t, static_cast<ui_type>(p), std::integral_constant<bool, false>());
  91 }
  92
  93 } // namespace detail
  94
  95 template <typename T, typename U>
  96 inline typename std::enable_if<boost::multiprecision::detail::is_integral<U>::value>::type eval_pow(T& result, const T& t, const U& p)
  97 {
  98    detail::pow_imp(result, t, p, boost::multiprecision::detail::is_signed<U>());
  99 }
 100
 101 template <class T>
 102 void hyp0F0(T& H0F0, const T& x)
 103 {
 104    // Compute the series representation of Hypergeometric0F0 taken from
 105    // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F0/06/01/
 106    // There are no checks on input range or parameter boundaries.
 107
 108    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
 109
 110    BOOST_MP_ASSERT(&H0F0 != &x);
 111    long tol = boost::multiprecision::detail::digits2<number<T, et_on> >::value();
 112    T    t;
 113
 114    T x_pow_n_div_n_fact(x);
 115
 116    eval_add(H0F0, x_pow_n_div_n_fact, ui_type(1));
 117
 118    T lim;
 119    eval_ldexp(lim, H0F0, 1 - tol);
 120    if (eval_get_sign(lim) < 0)
 121       lim.negate();
 122
 123    ui_type n;
 124
 125    const unsigned series_limit =
 126        boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
 127            ? 100
 128            : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
 129    // Series expansion of hyperg_0f0(; ; x).
 130    for (n = 2; n < series_limit; ++n)
 131    {
 132       eval_multiply(x_pow_n_div_n_fact, x);
 133       eval_divide(x_pow_n_div_n_fact, n);
 134       eval_add(H0F0, x_pow_n_div_n_fact);
 135       bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0;
 136       if (neg)
 137          x_pow_n_div_n_fact.negate();
 138       if (lim.compare(x_pow_n_div_n_fact) > 0)
 139          break;
 140       if (neg)
 141          x_pow_n_div_n_fact.negate();
 142    }
 143    if (n >= series_limit)
 144       BOOST_MP_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge"));
 145 }
 146
 147 template <class T>
 148 void hyp1F0(T& H1F0, const T& a, const T& x)
 149 {
 150    // Compute the series representation of Hypergeometric1F0 taken from
 151    // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F0/06/01/01/
 152    // and also see the corresponding section for the power function (i.e. x^a).
 153    // There are no checks on input range or parameter boundaries.
 154
 155    using si_type = typename boost::multiprecision::detail::canonical<int, T>::type;
 156
 157    BOOST_MP_ASSERT(&H1F0 != &x);
 158    BOOST_MP_ASSERT(&H1F0 != &a);
 159
 160    T x_pow_n_div_n_fact(x);
 161    T pochham_a(a);
 162    T ap(a);
 163
 164    eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact);
 165    eval_add(H1F0, si_type(1));
 166    T lim;
 167    eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 168    if (eval_get_sign(lim) < 0)
 169       lim.negate();
 170
 171    si_type n;
 172    T       term, part;
 173
 174    const si_type series_limit =
 175        boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
 176            ? 100
 177            : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
 178    // Series expansion of hyperg_1f0(a; ; x).
 179    for (n = 2; n < series_limit; n++)
 180    {
 181       eval_multiply(x_pow_n_div_n_fact, x);
 182       eval_divide(x_pow_n_div_n_fact, n);
 183       eval_increment(ap);
 184       eval_multiply(pochham_a, ap);
 185       eval_multiply(term, pochham_a, x_pow_n_div_n_fact);
 186       eval_add(H1F0, term);
 187       if (eval_get_sign(term) < 0)
 188          term.negate();
 189       if (lim.compare(term) >= 0)
 190          break;
 191    }
 192    if (n >= series_limit)
 193       BOOST_MP_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge"));
 194 }
 195
 196 template <class T>
 197 void eval_exp(T& result, const T& x)
 198 {
 199    static_assert(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types.");
 200    if (&x == &result)
 201    {
 202       T temp;
 203       eval_exp(temp, x);
 204       result = temp;
 205       return;
 206    }
 207    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
 208    using si_type = typename boost::multiprecision::detail::canonical<int, T>::type     ;
 209    using exp_type = typename T::exponent_type                                           ;
 210    using canonical_exp_type = typename boost::multiprecision::detail::canonical<exp_type, T>::type;
 211
 212    // Handle special arguments.
 213    int  type  = eval_fpclassify(x);
 214    bool isneg = eval_get_sign(x) < 0;
 215    if (type == static_cast<int>(FP_NAN))
 216    {
 217       result = x;
 218       errno  = EDOM;
 219       return;
 220    }
 221    else if (type == static_cast<int>(FP_INFINITE))
 222    {
 223       if (isneg)
 224          result = ui_type(0u);
 225       else
 226          result = x;
 227       return;
 228    }
 229    else if (type == static_cast<int>(FP_ZERO))
 230    {
 231       result = ui_type(1);
 232       return;
 233    }
 234
 235    // Get local copy of argument and force it to be positive.
 236    T xx = x;
 237    T exp_series;
 238    if (isneg)
 239       xx.negate();
 240
 241    // Check the range of the argument.
 242    if (xx.compare(si_type(1)) <= 0)
 243    {
 244       //
 245       // Use series for exp(x) - 1:
 246       //
 247       T lim;
 248       BOOST_IF_CONSTEXPR(std::numeric_limits<number<T, et_on> >::is_specialized)
 249          lim = std::numeric_limits<number<T, et_on> >::epsilon().backend();
 250       else
 251       {
 252          result = ui_type(1);
 253          eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 254       }
 255       unsigned k = 2;
 256       exp_series = xx;
 257       result     = si_type(1);
 258       if (isneg)
 259          eval_subtract(result, exp_series);
 260       else
 261          eval_add(result, exp_series);
 262       eval_multiply(exp_series, xx);
 263       eval_divide(exp_series, ui_type(k));
 264       eval_add(result, exp_series);
 265       while (exp_series.compare(lim) > 0)
 266       {
 267          ++k;
 268          eval_multiply(exp_series, xx);
 269          eval_divide(exp_series, ui_type(k));
 270          if (isneg && (k & 1))
 271             eval_subtract(result, exp_series);
 272          else
 273             eval_add(result, exp_series);
 274       }
 275       return;
 276    }
 277
 278    // Check for pure-integer arguments which can be either signed or unsigned.
 279    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type ll;
 280    eval_trunc(exp_series, x);
 281    eval_convert_to(&ll, exp_series);
 282    if (x.compare(ll) == 0)
 283    {
 284       detail::pow_imp(result, get_constant_e<T>(), ll, std::integral_constant<bool, true>());
 285       return;
 286    }
 287    else if (exp_series.compare(x) == 0)
 288    {
 289       // We have a value that has no fractional part, but is too large to fit
 290       // in a long long, in this situation the code below will fail, so
 291       // we're just going to assume that this will overflow:
 292       if (isneg)
 293          result = ui_type(0);
 294       else
 295          result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
 296       return;
 297    }
 298
 299    // The algorithm for exp has been taken from MPFUN.
 300    // exp(t) = [ (1 + r + r^2/2! + r^3/3! + r^4/4! ...)^p2 ] * 2^n
 301    // where p2 is a power of 2 such as 2048, r = t_prime / p2, and
 302    // t_prime = t - n*ln2, with n chosen to minimize the absolute
 303    // value of t_prime. In the resulting Taylor series, which is
 304    // implemented as a hypergeometric function, |r| is bounded by
 305    // ln2 / p2. For small arguments, no scaling is done.
 306
 307    // Compute the exponential series of the (possibly) scaled argument.
 308
 309    eval_divide(result, xx, get_constant_ln2<T>());
 310    exp_type n;
 311    eval_convert_to(&n, result);
 312
 313    if (n == (std::numeric_limits<exp_type>::max)())
 314    {
 315       // Exponent is too large to fit in our exponent type:
 316       if (isneg)
 317          result = ui_type(0);
 318       else
 319          result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
 320       return;
 321    }
 322
 323    // The scaling is 2^11 = 2048.
 324    const si_type p2 = static_cast<si_type>(si_type(1) << 11);
 325
 326    eval_multiply(exp_series, get_constant_ln2<T>(), static_cast<canonical_exp_type>(n));
 327    eval_subtract(exp_series, xx);
 328    eval_divide(exp_series, p2);
 329    exp_series.negate();
 330    hyp0F0(result, exp_series);
 331
 332    detail::pow_imp(exp_series, result, p2, std::integral_constant<bool, true>());
 333    result = ui_type(1);
 334    eval_ldexp(result, result, n);
 335    eval_multiply(exp_series, result);
 336
 337    if (isneg)
 338       eval_divide(result, ui_type(1), exp_series);
 339    else
 340       result = exp_series;
 341 }
 342
 343 template <class T>
 344 void eval_log(T& result, const T& arg)
 345 {
 346    static_assert(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
 347    //
 348    // We use a variation of http://dlmf.nist.gov/4.45#i
 349    // using frexp to reduce the argument to x * 2^n,
 350    // then let y = x - 1 and compute:
 351    // log(x) = log(2) * n + log1p(1 + y)
 352    //
 353    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
 354    using exp_type = typename T::exponent_type                                           ;
 355    using canonical_exp_type = typename boost::multiprecision::detail::canonical<exp_type, T>::type;
 356    using fp_type = typename std::tuple_element<0, typename T::float_types>::type                  ;
 357    int                                                                          s = eval_signbit(arg);
 358    switch (eval_fpclassify(arg))
 359    {
 360    case FP_NAN:
 361       result = arg;
 362       errno  = EDOM;
 363       return;
 364    case FP_INFINITE:
 365       if (s)
 366          break;
 367       result = arg;
 368       return;
 369    case FP_ZERO:
 370       result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
 371       result.negate();
 372       errno = ERANGE;
 373       return;
 374    }
 375    if (s)
 376    {
 377       result = std::numeric_limits<number<T> >::quiet_NaN().backend();
 378       errno  = EDOM;
 379       return;
 380    }
 381
 382    exp_type e;
 383    T        t;
 384    eval_frexp(t, arg, &e);
 385    bool alternate = false;
 386
 387    if (t.compare(fp_type(2) / fp_type(3)) <= 0)
 388    {
 389       alternate = true;
 390       eval_ldexp(t, t, 1);
 391       --e;
 392    }
 393
 394    eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e));
 395    INSTRUMENT_BACKEND(result);
 396    eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */
 397    if (!alternate)
 398       t.negate(); /* 0 <= t <= 0.33333 */
 399    T pow = t;
 400    T lim;
 401    T t2;
 402
 403    if (alternate)
 404       eval_add(result, t);
 405    else
 406       eval_subtract(result, t);
 407
 408    BOOST_IF_CONSTEXPR(std::numeric_limits<number<T, et_on> >::is_specialized)
 409       eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend());
 410    else
 411       eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 412    if (eval_get_sign(lim) < 0)
 413       lim.negate();
 414    INSTRUMENT_BACKEND(lim);
 415
 416    ui_type k = 1;
 417    do
 418    {
 419       ++k;
 420       eval_multiply(pow, t);
 421       eval_divide(t2, pow, k);
 422       INSTRUMENT_BACKEND(t2);
 423       if (alternate && ((k & 1) != 0))
 424          eval_add(result, t2);
 425       else
 426          eval_subtract(result, t2);
 427       INSTRUMENT_BACKEND(result);
 428    } while (lim.compare(t2) < 0);
 429 }
 430
 431 template <class T>
 432 const T& get_constant_log10()
 433 {
 434    static BOOST_MP_THREAD_LOCAL T             result;
 435    static BOOST_MP_THREAD_LOCAL long digits = 0;
 436    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 437    {
 438       using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
 439       T                                                                            ten;
 440       ten = ui_type(10u);
 441       eval_log(result, ten);
 442       digits = boost::multiprecision::detail::digits2<number<T> >::value();
 443    }
 444
 445    return result;
 446 }
 447
 448 template <class T>
 449 void eval_log10(T& result, const T& arg)
 450 {
 451    static_assert(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types.");
 452    eval_log(result, arg);
 453    eval_divide(result, get_constant_log10<T>());
 454 }
 455
 456 template <class R, class T>
 457 inline void eval_log2(R& result, const T& a)
 458 {
 459    eval_log(result, a);
 460    eval_divide(result, get_constant_ln2<R>());
 461 }
 462
 463 template <typename T>
 464 inline void eval_pow(T& result, const T& x, const T& a)
 465 {
 466    static_assert(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types.");
 467    using si_type = typename boost::multiprecision::detail::canonical<int, T>::type;
 468    using fp_type = typename std::tuple_element<0, typename T::float_types>::type             ;
 469
 470    if ((&result == &x) || (&result == &a))
 471    {
 472       T t;
 473       eval_pow(t, x, a);
 474       result = t;
 475       return;
 476    }
 477
 478    if ((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
 479    {
 480       result = x;
 481       return;
 482    }
 483    if (a.compare(si_type(0)) == 0)
 484    {
 485       result = si_type(1);
 486       return;
 487    }
 488
 489    int type = eval_fpclassify(x);
 490
 491    switch (type)
 492    {
 493    case FP_ZERO:
 494       switch (eval_fpclassify(a))
 495       {
 496       case FP_ZERO:
 497          result = si_type(1);
 498          break;
 499       case FP_NAN:
 500          result = a;
 501          break;
 502       case FP_NORMAL: {
 503          // Need to check for a an odd integer as a special case:
 504          BOOST_MP_TRY
 505          {
 506             typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type i;
 507             eval_convert_to(&i, a);
 508             if (a.compare(i) == 0)
 509             {
 510                if (eval_signbit(a))
 511                {
 512                   if (i & 1)
 513                   {
 514                      result = std::numeric_limits<number<T> >::infinity().backend();
 515                      if (eval_signbit(x))
 516                         result.negate();
 517                      errno = ERANGE;
 518                   }
 519                   else
 520                   {
 521                      result = std::numeric_limits<number<T> >::infinity().backend();
 522                      errno  = ERANGE;
 523                   }
 524                }
 525                else if (i & 1)
 526                {
 527                   result = x;
 528                }
 529                else
 530                   result = si_type(0);
 531                return;
 532             }
 533          }
 534          BOOST_MP_CATCH(const std::exception&)
 535          {
 536             // fallthrough..
 537          }
 538          BOOST_MP_CATCH_END
 539          BOOST_FALLTHROUGH;
 540       }
 541       default:
 542          if (eval_signbit(a))
 543          {
 544             result = std::numeric_limits<number<T> >::infinity().backend();
 545             errno  = ERANGE;
 546          }
 547          else
 548             result = x;
 549          break;
 550       }
 551       return;
 552    case FP_NAN:
 553       result = x;
 554       errno  = ERANGE;
 555       return;
 556    default:;
 557    }
 558
 559    int s = eval_get_sign(a);
 560    if (s == 0)
 561    {
 562       result = si_type(1);
 563       return;
 564    }
 565
 566    if (s < 0)
 567    {
 568       T t, da;
 569       t = a;
 570       t.negate();
 571       eval_pow(da, x, t);
 572       eval_divide(result, si_type(1), da);
 573       return;
 574    }
 575
 576    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type an;
 577    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type max_an =
 578        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);
 579    typename boost::multiprecision::detail::canonical<std::intmax_t, T>::type min_an =
 580        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;
 581
 582    T fa;
 583    BOOST_MP_TRY
 584    {
 585       eval_convert_to(&an, a);
 586       if (a.compare(an) == 0)
 587       {
 588          detail::pow_imp(result, x, an, std::integral_constant<bool, true>());
 589          return;
 590       }
 591    }
 592    BOOST_MP_CATCH(const std::exception&)
 593    {
 594       // conversion failed, just fall through, value is not an integer.
 595       an = (std::numeric_limits<std::intmax_t>::max)();
 596    }
 597    BOOST_MP_CATCH_END
 598    if ((eval_get_sign(x) < 0))
 599    {
 600       typename boost::multiprecision::detail::canonical<std::uintmax_t, T>::type aun;
 601       BOOST_MP_TRY
 602       {
 603          eval_convert_to(&aun, a);
 604          if (a.compare(aun) == 0)
 605          {
 606             fa = x;
 607             fa.negate();
 608             eval_pow(result, fa, a);
 609             if (aun & 1u)
 610                result.negate();
 611             return;
 612          }
 613       }
 614       BOOST_MP_CATCH(const std::exception&)
 615       {
 616          // conversion failed, just fall through, value is not an integer.
 617       }
 618       BOOST_MP_CATCH_END
 619
 620       eval_floor(result, a);
 621       // -1^INF is a special case in C99:
 622       if ((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
 623       {
 624          result = si_type(1);
 625       }
 626       else if (a.compare(result) == 0)
 627       {
 628          // exponent is so large we have no fractional part:
 629          if (x.compare(si_type(-1)) < 0)
 630          {
 631             result = std::numeric_limits<number<T, et_on> >::infinity().backend();
 632          }
 633          else
 634          {
 635             result = si_type(0);
 636          }
 637       }
 638       else if (type == FP_INFINITE)
 639       {
 640          result = std::numeric_limits<number<T, et_on> >::infinity().backend();
 641       }
 642       else BOOST_IF_CONSTEXPR (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
 643       {
 644          result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
 645          errno  = EDOM;
 646       }
 647       else
 648       {
 649          BOOST_MP_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type."));
 650       }
 651       return;
 652    }
 653
 654    T t, da;
 655
 656    eval_subtract(da, a, an);
 657
 658    if ((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
 659    {
 660       if (a.compare(fp_type(1e-5f)) <= 0)
 661       {
 662          // Series expansion for small a.
 663          eval_log(t, x);
 664          eval_multiply(t, a);
 665          hyp0F0(result, t);
 666          return;
 667       }
 668       else
 669       {
 670          // Series expansion for moderately sized x. Note that for large power of a,
 671          // the power of the integer part of a is calculated using the pown function.
 672          if (an)
 673          {
 674             da.negate();
 675             t = si_type(1);
 676             eval_subtract(t, x);
 677             hyp1F0(result, da, t);
 678             detail::pow_imp(t, x, an, std::integral_constant<bool, true>());
 679             eval_multiply(result, t);
 680          }
 681          else
 682          {
 683             da = a;
 684             da.negate();
 685             t = si_type(1);
 686             eval_subtract(t, x);
 687             hyp1F0(result, da, t);
 688          }
 689       }
 690    }
 691    else
 692    {
 693       // Series expansion for pow(x, a). Note that for large power of a, the power
 694       // of the integer part of a is calculated using the pown function.
 695       if (an)
 696       {
 697          eval_log(t, x);
 698          eval_multiply(t, da);
 699          eval_exp(result, t);
 700          detail::pow_imp(t, x, an, std::integral_constant<bool, true>());
 701          eval_multiply(result, t);
 702       }
 703       else
 704       {
 705          eval_log(t, x);
 706          eval_multiply(t, a);
 707          eval_exp(result, t);
 708       }
 709    }
 710 }
 711
 712 template <class T, class A>
 713 #if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
 714 inline typename std::enable_if<!boost::multiprecision::detail::is_integral<A>::value, void>::type
 715 #else
 716 inline typename std::enable_if<is_compatible_arithmetic_type<A, number<T> >::value && !boost::multiprecision::detail::is_integral<A>::value, void>::type
 717 #endif
 718 eval_pow(T& result, const T& x, const A& a)
 719 {
 720    // Note this one is restricted to float arguments since pow.hpp already has a version for
 721    // integer powers....
 722    using canonical_type = typename boost::multiprecision::detail::canonical<A, T>::type         ;
 723    using cast_type = typename std::conditional<std::is_same<A, canonical_type>::value, T, canonical_type>::type;
 724    cast_type                                                                      c;
 725    c = a;
 726    eval_pow(result, x, c);
 727 }
 728
 729 template <class T, class A>
 730 #if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
 731 inline void
 732 #else
 733 inline typename std::enable_if<is_compatible_arithmetic_type<A, number<T> >::value, void>::type
 734 #endif
 735 eval_pow(T& result, const A& x, const T& a)
 736 {
 737    using canonical_type = typename boost::multiprecision::detail::canonical<A, T>::type         ;
 738    using cast_type = typename std::conditional<std::is_same<A, canonical_type>::value, T, canonical_type>::type;
 739    cast_type                                                                      c;
 740    c = x;
 741    eval_pow(result, c, a);
 742 }
 743
 744 template <class T>
 745 void eval_exp2(T& result, const T& arg)
 746 {
 747    static_assert(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
 748
 749    // Check for pure-integer arguments which can be either signed or unsigned.
 750    typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
 751    T                                                                                     temp;
 752    BOOST_MP_TRY
 753    {
 754       eval_trunc(temp, arg);
 755       eval_convert_to(&i, temp);
 756       if (arg.compare(i) == 0)
 757       {
 758          temp = static_cast<typename std::tuple_element<0, typename T::unsigned_types>::type>(1u);
 759          eval_ldexp(result, temp, i);
 760          return;
 761       }
 762    }
 763    #ifdef BOOST_MP_MATH_AVAILABLE
 764    BOOST_MP_CATCH(const boost::math::rounding_error&)
 765    { /* Fallthrough */
 766    }
 767    #endif
 768    BOOST_MP_CATCH(const std::runtime_error&)
 769    { /* Fallthrough */
 770    }
 771    BOOST_MP_CATCH_END
 772
 773    temp = static_cast<typename std::tuple_element<0, typename T::unsigned_types>::type>(2u);
 774    eval_pow(result, temp, arg);
 775 }
 776
 777 namespace detail {
 778
 779 template <class T>
 780 void small_sinh_series(T x, T& result)
 781 {
 782    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
 783    bool                                                                         neg = eval_get_sign(x) < 0;
 784    if (neg)
 785       x.negate();
 786    T p(x);
 787    T mult(x);
 788    eval_multiply(mult, x);
 789    result    = x;
 790    ui_type k = 1;
 791
 792    T lim(x);
 793    eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 794
 795    do
 796    {
 797       eval_multiply(p, mult);
 798       eval_divide(p, ++k);
 799       eval_divide(p, ++k);
 800       eval_add(result, p);
 801    } while (p.compare(lim) >= 0);
 802    if (neg)
 803       result.negate();
 804 }
 805
 806 template <class T>
 807 void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
 808 {
 809    using ui_type = typename boost::multiprecision::detail::canonical<unsigned, T>::type;
 810    using fp_type = typename std::tuple_element<0, typename T::float_types>::type                  ;
 811
 812    switch (eval_fpclassify(x))
 813    {
 814    case FP_NAN:
 815       errno = EDOM;
 816       // fallthrough...
 817    case FP_INFINITE:
 818       if (p_sinh)
 819          *p_sinh = x;
 820       if (p_cosh)
 821       {
 822          *p_cosh = x;
 823          if (eval_get_sign(x) < 0)
 824             p_cosh->negate();
 825       }
 826       return;
 827    case FP_ZERO:
 828       if (p_sinh)
 829          *p_sinh = x;
 830       if (p_cosh)
 831          *p_cosh = ui_type(1);
 832       return;
 833    default:;
 834    }
 835
 836    bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;
 837
 838    if (p_cosh || !small_sinh)
 839    {
 840       T e_px, e_mx;
 841       eval_exp(e_px, x);
 842       eval_divide(e_mx, ui_type(1), e_px);
 843       if (eval_signbit(e_mx) != eval_signbit(e_px))
 844          e_mx.negate(); // Handles lack of signed zero in some types
 845
 846       if (p_sinh)
 847       {
 848          if (small_sinh)
 849          {
 850             small_sinh_series(x, *p_sinh);
 851          }
 852          else
 853          {
 854             eval_subtract(*p_sinh, e_px, e_mx);
 855             eval_ldexp(*p_sinh, *p_sinh, -1);
 856          }
 857       }
 858       if (p_cosh)
 859       {
 860          eval_add(*p_cosh, e_px, e_mx);
 861          eval_ldexp(*p_cosh, *p_cosh, -1);
 862       }
 863    }
 864    else
 865    {
 866       small_sinh_series(x, *p_sinh);
 867    }
 868 }
 869
 870 } // namespace detail
 871
 872 template <class T>
 873 inline void eval_sinh(T& result, const T& x)
 874 {
 875    static_assert(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types.");
 876    detail::sinhcosh(x, &result, static_cast<T*>(0));
 877 }
 878
 879 template <class T>
 880 inline void eval_cosh(T& result, const T& x)
 881 {
 882    static_assert(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types.");
 883    detail::sinhcosh(x, static_cast<T*>(0), &result);
 884 }
 885
 886 template <class T>
 887 inline void eval_tanh(T& result, const T& x)
 888 {
 889    static_assert(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
 890    T c;
 891    detail::sinhcosh(x, &result, &c);
 892    if ((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE))
 893    {
 894       bool s = eval_signbit(result) != eval_signbit(c);
 895       result = static_cast<typename std::tuple_element<0, typename T::unsigned_types>::type>(1u);
 896       if (s)
 897          result.negate();
 898       return;
 899    }
 900    eval_divide(result, c);
 901 }
 902
 903 #ifdef BOOST_MSVC
 904 #pragma warning(pop)
 905 #endif