boost/multiprecision/detail/functions/pow.hpp

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