[ceph.git] / ceph / src / boost / libs / multiprecision / include / boost / multiprecision / detail / functions / pow.hpp


// Copyright Christopher Kormanyos 2002 - 2013.
// Copyright 2011 - 2013 John Maddock. Distributed under the Boost
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

// This work is based on an earlier work:
// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
//
// This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
// 

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:6326)  // comparison of two constants
#endif

namespace detail{

template<typename T, typename U> 
inline void pow_imp(T& result, const T& t, const U& p, const mpl::false_&)
{
   // Compute the pure power of typename T t^p.
   // Use the S-and-X binary method, as described in
   // D. E. Knuth, "The Art of Computer Programming", Vol. 2,
   // Section 4.6.3 . The resulting computational complexity
   // is order log2[abs(p)].

   typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type;

   if(&result == &t)
   {
      T temp;
      pow_imp(temp, t, p, mpl::false_());
      result = temp;
      return;
   }

   // This will store the result.
   if(U(p % U(2)) != U(0))
   {
      result = t;
   }
   else
      result = int_type(1);

   U p2(p);

   // The variable x stores the binary powers of t.
   T x(t);

   while(U(p2 /= 2) != U(0))
   {
      // Square x for each binary power.
      eval_multiply(x, x);

      const bool has_binary_power = (U(p2 % U(2)) != U(0));

      if(has_binary_power)
      {
         // Multiply the result with each binary power contained in the exponent.
         eval_multiply(result, x);
      }
   }
}

template<typename T, typename U> 
inline void pow_imp(T& result, const T& t, const U& p, const mpl::true_&)
{
   // Signed integer power, just take care of the sign then call the unsigned version:
   typedef typename boost::multiprecision::detail::canonical<U, T>::type  int_type;
   typedef typename make_unsigned<U>::type                                ui_type;

   if(p < 0)
   {
      T temp;
      temp = static_cast<int_type>(1);
      T denom;
      pow_imp(denom, t, static_cast<ui_type>(-p), mpl::false_());
      eval_divide(result, temp, denom);
      return;
   }
   pow_imp(result, t, static_cast<ui_type>(p), mpl::false_());
}

} // namespace detail

template<typename T, typename U> 
inline typename enable_if<is_integral<U> >::type eval_pow(T& result, const T& t, const U& p)
{
   detail::pow_imp(result, t, p, boost::is_signed<U>());
}

template <class T>
void hyp0F0(T& H0F0, const T& x)
{
   // Compute the series representation of Hypergeometric0F0 taken from
   // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F0/06/01/
   // There are no checks on input range or parameter boundaries.

   typedef typename mpl::front<typename T::unsigned_types>::type ui_type;

   BOOST_ASSERT(&H0F0 != &x);
   long tol = boost::multiprecision::detail::digits2<number<T, et_on> >::value();
   T t;

   T x_pow_n_div_n_fact(x);

   eval_add(H0F0, x_pow_n_div_n_fact, ui_type(1));

   T lim;
   eval_ldexp(lim, H0F0, 1 - tol);
   if(eval_get_sign(lim) < 0)
      lim.negate();

   ui_type n;

   const unsigned series_limit = 
      boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
      ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
   // Series expansion of hyperg_0f0(; ; x).
   for(n = 2; n < series_limit; ++n)
   {
      eval_multiply(x_pow_n_div_n_fact, x);
      eval_divide(x_pow_n_div_n_fact, n);
      eval_add(H0F0, x_pow_n_div_n_fact);
      bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0;
      if(neg)
         x_pow_n_div_n_fact.negate();
      if(lim.compare(x_pow_n_div_n_fact) > 0)
         break;
      if(neg)
         x_pow_n_div_n_fact.negate();
   }
   if(n >= series_limit)
      BOOST_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge"));
}

template <class T>
void hyp1F0(T& H1F0, const T& a, const T& x)
{
   // Compute the series representation of Hypergeometric1F0 taken from
   // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F0/06/01/01/
   // and also see the corresponding section for the power function (i.e. x^a).
   // There are no checks on input range or parameter boundaries.

   typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;

   BOOST_ASSERT(&H1F0 != &x);
   BOOST_ASSERT(&H1F0 != &a);

   T x_pow_n_div_n_fact(x);
   T pochham_a         (a);
   T ap                (a);

   eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact);
   eval_add(H1F0, si_type(1));
   T lim;
   eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
   if(eval_get_sign(lim) < 0)
      lim.negate();

   si_type n;
   T term, part;

   const si_type series_limit =
      boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
      ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
   // Series expansion of hyperg_1f0(a; ; x).
   for(n = 2; n < series_limit; n++)
   {
      eval_multiply(x_pow_n_div_n_fact, x);
      eval_divide(x_pow_n_div_n_fact, n);
      eval_increment(ap);
      eval_multiply(pochham_a, ap);
      eval_multiply(term, pochham_a, x_pow_n_div_n_fact);
      eval_add(H1F0, term);
      if(eval_get_sign(term) < 0)
         term.negate();
      if(lim.compare(term) >= 0)
         break;
   }
   if(n >= series_limit)
      BOOST_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge"));
}

template <class T>
void eval_exp(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types.");
   if(&x == &result)
   {
      T temp;
      eval_exp(temp, x);
      result = temp;
      return;
   }
   typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
   typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
   typedef typename T::exponent_type exp_type;
   typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;

   // Handle special arguments.
   int type = eval_fpclassify(x);
   bool isneg = eval_get_sign(x) < 0;
   if(type == (int)FP_NAN)
   {
      result = x;
      return;
   }
   else if(type == (int)FP_INFINITE)
   {
      result = x;
      if(isneg)
         result = ui_type(0u);
      else 
         result = x;
      return;
   }
   else if(type == (int)FP_ZERO)
   {
      result = ui_type(1);
      return;
   }

   // Get local copy of argument and force it to be positive.
   T xx = x;
   T exp_series;
   if(isneg)
      xx.negate();

   // Check the range of the argument.
   if(xx.compare(si_type(1)) <= 0)
   {
      //
      // Use series for exp(x) - 1:
      //
      T lim;
      if(std::numeric_limits<number<T, et_on> >::is_specialized)
         lim = std::numeric_limits<number<T, et_on> >::epsilon().backend();
      else
      {
         result = ui_type(1);
         eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
      }
      unsigned k = 2;
      exp_series = xx;
      result = si_type(1);
      if(isneg)
         eval_subtract(result, exp_series);
      else
         eval_add(result, exp_series);
      eval_multiply(exp_series, xx);
      eval_divide(exp_series, ui_type(k));
      eval_add(result, exp_series);
      while(exp_series.compare(lim) > 0)
      {
         ++k;
         eval_multiply(exp_series, xx);
         eval_divide(exp_series, ui_type(k));
         if(isneg && (k&1))
            eval_subtract(result, exp_series);
         else
            eval_add(result, exp_series);
      }
      return;
   }

   // Check for pure-integer arguments which can be either signed or unsigned.
   typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type ll;
   eval_trunc(exp_series, x);
   eval_convert_to(&ll, exp_series);
   if(x.compare(ll) == 0)
   {
      detail::pow_imp(result, get_constant_e<T>(), ll, mpl::true_());
      return;
   }

   // The algorithm for exp has been taken from MPFUN.
   // exp(t) = [ (1 + r + r^2/2! + r^3/3! + r^4/4! ...)^p2 ] * 2^n
   // where p2 is a power of 2 such as 2048, r = t_prime / p2, and
   // t_prime = t - n*ln2, with n chosen to minimize the absolute
   // value of t_prime. In the resulting Taylor series, which is
   // implemented as a hypergeometric function, |r| is bounded by
   // ln2 / p2. For small arguments, no scaling is done.

   // Compute the exponential series of the (possibly) scaled argument.

   eval_divide(result, xx, get_constant_ln2<T>());
   exp_type n;
   eval_convert_to(&n, result);

   // The scaling is 2^11 = 2048.
   const si_type p2 = static_cast<si_type>(si_type(1) << 11);

   eval_multiply(exp_series, get_constant_ln2<T>(), static_cast<canonical_exp_type>(n));
   eval_subtract(exp_series, xx);
   eval_divide(exp_series, p2);
   exp_series.negate();
   hyp0F0(result, exp_series);

   detail::pow_imp(exp_series, result, p2, mpl::true_());
   result = ui_type(1);
   eval_ldexp(result, result, n);
   eval_multiply(exp_series, result);

   if(isneg)
      eval_divide(result, ui_type(1), exp_series);
   else
      result = exp_series;
}

template <class T>
void eval_log(T& result, const T& arg)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
   //
   // We use a variation of http://dlmf.nist.gov/4.45#i
   // using frexp to reduce the argument to x * 2^n,
   // then let y = x - 1 and compute:
   // log(x) = log(2) * n + log1p(1 + y)
   //
   typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
   typedef typename T::exponent_type exp_type;
   typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
   typedef typename mpl::front<typename T::float_types>::type fp_type;

   exp_type e;
   T t;
   eval_frexp(t, arg, &e);
   bool alternate = false;

   if(t.compare(fp_type(2) / fp_type(3)) <= 0)
   {
      alternate = true;
      eval_ldexp(t, t, 1);
      --e;
   }
   
   eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e));
   INSTRUMENT_BACKEND(result);
   eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */
   if(!alternate)
      t.negate(); /* 0 <= t <= 0.33333 */
   T pow = t;
   T lim;
   T t2;

   if(alternate)
      eval_add(result, t);
   else
      eval_subtract(result, t);

   if(std::numeric_limits<number<T, et_on> >::is_specialized)
      eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend());
   else
      eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
   if(eval_get_sign(lim) < 0)
      lim.negate();
   INSTRUMENT_BACKEND(lim);

   ui_type k = 1;
   do
   {
      ++k;
      eval_multiply(pow, t);
      eval_divide(t2, pow, k);
      INSTRUMENT_BACKEND(t2);
      if(alternate && ((k & 1) != 0))
         eval_add(result, t2);
      else
         eval_subtract(result, t2);
      INSTRUMENT_BACKEND(result);
   }while(lim.compare(t2) < 0);
}

template <class T>
const T& get_constant_log10()
{
   static BOOST_MP_THREAD_LOCAL T result;
   static BOOST_MP_THREAD_LOCAL bool b = false;
   static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value();
   if(!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
   {
      typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
      T ten;
      ten = ui_type(10u);
      eval_log(result, ten);
      b = true;
      digits = boost::multiprecision::detail::digits2<number<T> >::value();
   }

   constant_initializer<T, &get_constant_log10<T> >::do_nothing();

   return result;
}

template <class T>
void eval_log10(T& result, const T& arg)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types.");
   eval_log(result, arg);
   eval_divide(result, get_constant_log10<T>());
}

template <class R, class T>
inline void eval_log2(R& result, const T& a)
{
   eval_log(result, a);
   eval_divide(result, get_constant_ln2<R>());
}

template<typename T> 
inline void eval_pow(T& result, const T& x, const T& a)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types.");
   typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
   typedef typename mpl::front<typename T::float_types>::type fp_type;

   if((&result == &x) || (&result == &a))
   {
      T t;
      eval_pow(t, x, a);
      result = t;
      return;
   }

   if(a.compare(si_type(1)) == 0)
   {
      result = x;
      return;
   }

   int type = eval_fpclassify(x);

   switch(type)
   {
   case FP_INFINITE:
      result = x;
      return;
   case FP_ZERO:
      switch(eval_fpclassify(a))
      {
      case FP_ZERO:
         result = si_type(1);
         break;
      case FP_NAN:
         result = a;
         break;
      default:
         result = x;
         break;
      }
      return;
   case FP_NAN:
      result = x;
      return;
   default: ;
   }

   int s = eval_get_sign(a);
   if(s == 0)
   {
      result = si_type(1);
      return;
   }

   if(s < 0)
   {
      T t, da;
      t = a;
      t.negate();
      eval_pow(da, x, t);
      eval_divide(result, si_type(1), da);
      return;
   }
   
   typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an;
   T fa;
#ifndef BOOST_NO_EXCEPTIONS
   try
   {
#endif
      eval_convert_to(&an, a);
      if(a.compare(an) == 0)
      {
         detail::pow_imp(result, x, an, mpl::true_());
         return;
      }
#ifndef BOOST_NO_EXCEPTIONS
   }
   catch(const std::exception&)
   {
      // conversion failed, just fall through, value is not an integer.
      an = (std::numeric_limits<boost::intmax_t>::max)();
   }
#endif
   if((eval_get_sign(x) < 0))
   {
      typename boost::multiprecision::detail::canonical<boost::uintmax_t, T>::type aun;
#ifndef BOOST_NO_EXCEPTIONS
      try
      {
#endif
         eval_convert_to(&aun, a);
         if(a.compare(aun) == 0)
         {
            fa = x;
            fa.negate();
            eval_pow(result, fa, a);
            if(aun & 1u)
               result.negate();
            return;
         }
#ifndef BOOST_NO_EXCEPTIONS
      }
      catch(const std::exception&)
      {
         // conversion failed, just fall through, value is not an integer.
      }
#endif
      if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
         result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
      else
      {
         BOOST_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type."));
      }
      return;
   }

   T t, da;

   eval_subtract(da, a, an);

   if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0))
   {
      if(a.compare(fp_type(1e-5f)) <= 0)
      {
         // Series expansion for small a.
         eval_log(t, x);
         eval_multiply(t, a);
         hyp0F0(result, t);
         return;
      }
      else
      {
         // Series expansion for moderately sized x. Note that for large power of a,
         // the power of the integer part of a is calculated using the pown function.
         if(an)
         {
            da.negate();
            t = si_type(1);
            eval_subtract(t, x);
            hyp1F0(result, da, t);
            detail::pow_imp(t, x, an, mpl::true_());
            eval_multiply(result, t);
         }
         else
         {
            da = a;
            da.negate();
            t = si_type(1);
            eval_subtract(t, x);
            hyp1F0(result, da, t);
         }
      }
   }
   else
   {
      // Series expansion for pow(x, a). Note that for large power of a, the power
      // of the integer part of a is calculated using the pown function.
      if(an)
      {
         eval_log(t, x);
         eval_multiply(t, da);
         eval_exp(result, t);
         detail::pow_imp(t, x, an, mpl::true_());
         eval_multiply(result, t);
      }
      else
      {
         eval_log(t, x);
         eval_multiply(t, a);
         eval_exp(result, t);
      }
   }
}

template<class T, class A> 
inline typename enable_if<is_floating_point<A>, void>::type eval_pow(T& result, const T& x, const A& a)
{
   // Note this one is restricted to float arguments since pow.hpp already has a version for
   // integer powers....
   typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
   typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
   cast_type c;
   c = a;
   eval_pow(result, x, c);
}

template<class T, class A> 
inline typename enable_if<is_arithmetic<A>, void>::type eval_pow(T& result, const A& x, const T& a)
{
   typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
   typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
   cast_type c;
   c = x;
   eval_pow(result, c, a);
}

template <class T>
void eval_exp2(T& result, const T& arg)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");

   // Check for pure-integer arguments which can be either signed or unsigned.
   typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
   T temp;
   eval_trunc(temp, arg);
   eval_convert_to(&i, temp);
   if(arg.compare(i) == 0)
   {
      temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
      eval_ldexp(result, temp, i);
      return;
   }

   temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u);
   eval_pow(result, temp, arg);
}

namespace detail{

   template <class T>
   void small_sinh_series(T x, T& result)
   {
      typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
      bool neg = eval_get_sign(x) < 0;
      if(neg)
         x.negate();
      T p(x);
      T mult(x);
      eval_multiply(mult, x);
      result = x;
      ui_type k = 1;

      T lim(x);
      eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());

      do
      {
         eval_multiply(p, mult);
         eval_divide(p, ++k);
         eval_divide(p, ++k);
         eval_add(result, p);
      }while(p.compare(lim) >= 0);
      if(neg)
         result.negate();
   }

   template <class T>
   void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
   {
      typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
      typedef typename mpl::front<typename T::float_types>::type fp_type;

      switch(eval_fpclassify(x))
      {
      case FP_NAN:
      case FP_INFINITE:
         if(p_sinh)
            *p_sinh = x;
         if(p_cosh)
         {
            *p_cosh = x;
            if(eval_get_sign(x) < 0)
               p_cosh->negate();
         }
         return;
      case FP_ZERO:
         if(p_sinh)
            *p_sinh = x;
         if(p_cosh)
            *p_cosh = ui_type(1);
         return;
      default: ;
      }

      bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;

      if(p_cosh || !small_sinh)
      {
         T e_px, e_mx;
         eval_exp(e_px, x);
         eval_divide(e_mx, ui_type(1), e_px);

         if(p_sinh) 
         { 
            if(small_sinh)
            {
               small_sinh_series(x, *p_sinh);
            }
            else
            {
               eval_subtract(*p_sinh, e_px, e_mx);
               eval_ldexp(*p_sinh, *p_sinh, -1);
            }
         }
         if(p_cosh) 
         { 
            eval_add(*p_cosh, e_px, e_mx);
            eval_ldexp(*p_cosh, *p_cosh, -1); 
         }
      }
      else
      {
         small_sinh_series(x, *p_sinh);
      }
   }

} // namespace detail

template <class T>
inline void eval_sinh(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types.");
   detail::sinhcosh(x, &result, static_cast<T*>(0));
}

template <class T>
inline void eval_cosh(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types.");
   detail::sinhcosh(x, static_cast<T*>(0), &result);
}

template <class T>
inline void eval_tanh(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
  T c;
  detail::sinhcosh(x, &result, &c);
  eval_divide(result, c);
}

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
Commit	Line	Data
7c673cae FG	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<is_integral<U> >::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 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
123	// Series expansion of hyperg_0f0(; ; x).
124	for(n = 2; n < series_limit; ++n)
125	{
126	eval_multiply(x_pow_n_div_n_fact, x);
127	eval_divide(x_pow_n_div_n_fact, n);
128	eval_add(H0F0, x_pow_n_div_n_fact);
129	bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0;
130	if(neg)
131	x_pow_n_div_n_fact.negate();
132	if(lim.compare(x_pow_n_div_n_fact) > 0)
133	break;
134	if(neg)
135	x_pow_n_div_n_fact.negate();
136	}
137	if(n >= series_limit)
138	BOOST_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge"));
139	}
140
141	template <class T>
142	void hyp1F0(T& H1F0, const T& a, const T& x)
143	{
144	// Compute the series representation of Hypergeometric1F0 taken from
145	// http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F0/06/01/01/
146	// and also see the corresponding section for the power function (i.e. x^a).
147	// There are no checks on input range or parameter boundaries.
148
149	typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
150
151	BOOST_ASSERT(&H1F0 != &x);
152	BOOST_ASSERT(&H1F0 != &a);
153
154	T x_pow_n_div_n_fact(x);
155	T pochham_a (a);
156	T ap (a);
157
158	eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact);
159	eval_add(H1F0, si_type(1));
160	T lim;
161	eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
162	if(eval_get_sign(lim) < 0)
163	lim.negate();
164
165	si_type n;
166	T term, part;
167
168	const si_type series_limit =
169	boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
170	? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
171	// Series expansion of hyperg_1f0(a; ; x).
172	for(n = 2; n < series_limit; n++)
173	{
174	eval_multiply(x_pow_n_div_n_fact, x);
175	eval_divide(x_pow_n_div_n_fact, n);
176	eval_increment(ap);
177	eval_multiply(pochham_a, ap);
178	eval_multiply(term, pochham_a, x_pow_n_div_n_fact);
179	eval_add(H1F0, term);
180	if(eval_get_sign(term) < 0)
181	term.negate();
182	if(lim.compare(term) >= 0)
183	break;
184	}
185	if(n >= series_limit)
186	BOOST_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge"));
187	}
188
189	template <class T>
190	void eval_exp(T& result, const T& x)
191	{
192	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types.");
193	if(&x == &result)
194	{
195	T temp;
196	eval_exp(temp, x);
197	result = temp;
198	return;
199	}
200	typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
201	typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
202	typedef typename T::exponent_type exp_type;
203	typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
204
205	// Handle special arguments.
206	int type = eval_fpclassify(x);
207	bool isneg = eval_get_sign(x) < 0;
208	if(type == (int)FP_NAN)
209	{
210	result = x;
211	return;
212	}
213	else if(type == (int)FP_INFINITE)
214	{
215	result = x;
216	if(isneg)
217	result = ui_type(0u);
218	else
219	result = x;
220	return;
221	}
222	else if(type == (int)FP_ZERO)
223	{
224	result = ui_type(1);
225	return;
226	}
227
228	// Get local copy of argument and force it to be positive.
229	T xx = x;
230	T exp_series;
231	if(isneg)
232	xx.negate();
233
234	// Check the range of the argument.
235	if(xx.compare(si_type(1)) <= 0)
236	{
237	//
238	// Use series for exp(x) - 1:
239	//
240	T lim;
241	if(std::numeric_limits<number<T, et_on> >::is_specialized)
242	lim = std::numeric_limits<number<T, et_on> >::epsilon().backend();
243	else
244	{
245	result = ui_type(1);
246	eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
247	}
248	unsigned k = 2;
249	exp_series = xx;
250	result = si_type(1);
251	if(isneg)
252	eval_subtract(result, exp_series);
253	else
254	eval_add(result, exp_series);
255	eval_multiply(exp_series, xx);
256	eval_divide(exp_series, ui_type(k));
257	eval_add(result, exp_series);
258	while(exp_series.compare(lim) > 0)
259	{
260	++k;
261	eval_multiply(exp_series, xx);
262	eval_divide(exp_series, ui_type(k));
263	if(isneg && (k&1))
264	eval_subtract(result, exp_series);
265	else
266	eval_add(result, exp_series);
267	}
268	return;
269	}
270
271	// Check for pure-integer arguments which can be either signed or unsigned.
272	typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type ll;
273	eval_trunc(exp_series, x);
274	eval_convert_to(&ll, exp_series);
275	if(x.compare(ll) == 0)
276	{
277	detail::pow_imp(result, get_constant_e<T>(), ll, mpl::true_());
278	return;
279	}
280
281	// The algorithm for exp has been taken from MPFUN.
282	// exp(t) = [ (1 + r + r^2/2! + r^3/3! + r^4/4! ...)^p2 ] * 2^n
283	// where p2 is a power of 2 such as 2048, r = t_prime / p2, and
284	// t_prime = t - n*ln2, with n chosen to minimize the absolute
285	// value of t_prime. In the resulting Taylor series, which is
286	// implemented as a hypergeometric function, \|r\| is bounded by
287	// ln2 / p2. For small arguments, no scaling is done.
288
289	// Compute the exponential series of the (possibly) scaled argument.
290
291	eval_divide(result, xx, get_constant_ln2<T>());
292	exp_type n;
293	eval_convert_to(&n, result);
294
295	// The scaling is 2^11 = 2048.
296	const si_type p2 = static_cast<si_type>(si_type(1) << 11);
297
298	eval_multiply(exp_series, get_constant_ln2<T>(), static_cast<canonical_exp_type>(n));
299	eval_subtract(exp_series, xx);
300	eval_divide(exp_series, p2);
301	exp_series.negate();
302	hyp0F0(result, exp_series);
303
304	detail::pow_imp(exp_series, result, p2, mpl::true_());
305	result = ui_type(1);
306	eval_ldexp(result, result, n);
307	eval_multiply(exp_series, result);
308
309	if(isneg)
310	eval_divide(result, ui_type(1), exp_series);
311	else
312	result = exp_series;
313	}
314
315	template <class T>
316	void eval_log(T& result, const T& arg)
317	{
318	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
319	//
320	// We use a variation of http://dlmf.nist.gov/4.45#i
321	// using frexp to reduce the argument to x * 2^n,
322	// then let y = x - 1 and compute:
323	// log(x) = log(2) * n + log1p(1 + y)
324	//
325	typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
326	typedef typename T::exponent_type exp_type;
327	typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
328	typedef typename mpl::front<typename T::float_types>::type fp_type;
329
330	exp_type e;
331	T t;
332	eval_frexp(t, arg, &e);
333	bool alternate = false;
334
335	if(t.compare(fp_type(2) / fp_type(3)) <= 0)
336	{
337	alternate = true;
338	eval_ldexp(t, t, 1);
339	--e;
340	}
341
342	eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e));
343	INSTRUMENT_BACKEND(result);
344	eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */
345	if(!alternate)
346	t.negate(); /* 0 <= t <= 0.33333 */
347	T pow = t;
348	T lim;
349	T t2;
350
351	if(alternate)
352	eval_add(result, t);
353	else
354	eval_subtract(result, t);
355
356	if(std::numeric_limits<number<T, et_on> >::is_specialized)
357	eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend());
358	else
359	eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
360	if(eval_get_sign(lim) < 0)
361	lim.negate();
362	INSTRUMENT_BACKEND(lim);
363
364	ui_type k = 1;
365	do
366	{
367	++k;
368	eval_multiply(pow, t);
369	eval_divide(t2, pow, k);
370	INSTRUMENT_BACKEND(t2);
371	if(alternate && ((k & 1) != 0))
372	eval_add(result, t2);
373	else
374	eval_subtract(result, t2);
375	INSTRUMENT_BACKEND(result);
376	}while(lim.compare(t2) < 0);
377	}
378
379	template <class T>
380	const T& get_constant_log10()
381	{
382	static BOOST_MP_THREAD_LOCAL T result;
383	static BOOST_MP_THREAD_LOCAL bool b = false;
384	static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value();
385	if(!b \|\| (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
386	{
387	typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
388	T ten;
389	ten = ui_type(10u);
390	eval_log(result, ten);
391	b = true;
392	digits = boost::multiprecision::detail::digits2<number<T> >::value();
393	}
394
395	constant_initializer<T, &get_constant_log10<T> >::do_nothing();
396
397	return result;
398	}
399
400	template <class T>
401	void eval_log10(T& result, const T& arg)
402	{
403	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types.");
404	eval_log(result, arg);
405	eval_divide(result, get_constant_log10<T>());
406	}
407
408	template <class R, class T>
409	inline void eval_log2(R& result, const T& a)
410	{
411	eval_log(result, a);
412	eval_divide(result, get_constant_ln2<R>());
413	}
414
415	template<typename T>
416	inline void eval_pow(T& result, const T& x, const T& a)
417	{
418	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types.");
419	typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
420	typedef typename mpl::front<typename T::float_types>::type fp_type;
421
422	if((&result == &x) \|\| (&result == &a))
423	{
424	T t;
425	eval_pow(t, x, a);
426	result = t;
427	return;
428	}
429
430	if(a.compare(si_type(1)) == 0)
431	{
432	result = x;
433	return;
434	}
435
436	int type = eval_fpclassify(x);
437
438	switch(type)
439	{
440	case FP_INFINITE:
441	result = x;
442	return;
443	case FP_ZERO:
444	switch(eval_fpclassify(a))
445	{
446	case FP_ZERO:
447	result = si_type(1);
448	break;
449	case FP_NAN:
450	result = a;
451	break;
452	default:
453	result = x;
454	break;
455	}
456	return;
457	case FP_NAN:
458	result = x;
459	return;
460	default: ;
461	}
462
463	int s = eval_get_sign(a);
464	if(s == 0)
465	{
466	result = si_type(1);
467	return;
468	}
469
470	if(s < 0)
471	{
472	T t, da;
473	t = a;
474	t.negate();
475	eval_pow(da, x, t);
476	eval_divide(result, si_type(1), da);
477	return;
478	}
479
480	typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an;
481	T fa;
482	#ifndef BOOST_NO_EXCEPTIONS
483	try
484	{
485	#endif
486	eval_convert_to(&an, a);
487	if(a.compare(an) == 0)
488	{
489	detail::pow_imp(result, x, an, mpl::true_());
490	return;
491	}
492	#ifndef BOOST_NO_EXCEPTIONS
493	}
494	catch(const std::exception&)
495	{
496	// conversion failed, just fall through, value is not an integer.
497	an = (std::numeric_limits<boost::intmax_t>::max)();
498	}
499	#endif
500	if((eval_get_sign(x) < 0))
501	{
502	typename boost::multiprecision::detail::canonical<boost::uintmax_t, T>::type aun;
503	#ifndef BOOST_NO_EXCEPTIONS
504	try
505	{
506	#endif
507	eval_convert_to(&aun, a);
508	if(a.compare(aun) == 0)
509	{
510	fa = x;
511	fa.negate();
512	eval_pow(result, fa, a);
513	if(aun & 1u)
514	result.negate();
515	return;
516	}
517	#ifndef BOOST_NO_EXCEPTIONS
518	}
519	catch(const std::exception&)
520	{
521	// conversion failed, just fall through, value is not an integer.
522	}
523	#endif
524	if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
525	result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
526	else
527	{
528	BOOST_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type."));
529	}
530	return;
531	}
532
533	T t, da;
534
535	eval_subtract(da, a, an);
536
537	if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0))
538	{
539	if(a.compare(fp_type(1e-5f)) <= 0)
540	{
541	// Series expansion for small a.
542	eval_log(t, x);
543	eval_multiply(t, a);
544	hyp0F0(result, t);
545	return;
546	}
547	else
548	{
549	// Series expansion for moderately sized x. Note that for large power of a,
550	// the power of the integer part of a is calculated using the pown function.
551	if(an)
552	{
553	da.negate();
554	t = si_type(1);
555	eval_subtract(t, x);
556	hyp1F0(result, da, t);
557	detail::pow_imp(t, x, an, mpl::true_());
558	eval_multiply(result, t);
559	}
560	else
561	{
562	da = a;
563	da.negate();
564	t = si_type(1);
565	eval_subtract(t, x);
566	hyp1F0(result, da, t);
567	}
568	}
569	}
570	else
571	{
572	// Series expansion for pow(x, a). Note that for large power of a, the power
573	// of the integer part of a is calculated using the pown function.
574	if(an)
575	{
576	eval_log(t, x);
577	eval_multiply(t, da);
578	eval_exp(result, t);
579	detail::pow_imp(t, x, an, mpl::true_());
580	eval_multiply(result, t);
581	}
582	else
583	{
584	eval_log(t, x);
585	eval_multiply(t, a);
586	eval_exp(result, t);
587	}
588	}
589	}
590
591	template<class T, class A>
592	inline typename enable_if<is_floating_point<A>, void>::type eval_pow(T& result, const T& x, const A& a)
593	{
594	// Note this one is restricted to float arguments since pow.hpp already has a version for
595	// integer powers....
596	typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
597	typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
598	cast_type c;
599	c = a;
600	eval_pow(result, x, c);
601	}
602
603	template<class T, class A>
604	inline typename enable_if<is_arithmetic<A>, void>::type eval_pow(T& result, const A& x, const T& a)
605	{
606	typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
607	typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
608	cast_type c;
609	c = x;
610	eval_pow(result, c, a);
611	}
612
613	template <class T>
614	void eval_exp2(T& result, const T& arg)
615	{
616	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
617
618	// Check for pure-integer arguments which can be either signed or unsigned.
619	typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
620	T temp;
621	eval_trunc(temp, arg);
622	eval_convert_to(&i, temp);
623	if(arg.compare(i) == 0)
624	{
625	temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
626	eval_ldexp(result, temp, i);
627	return;
628	}
629
630	temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u);
631	eval_pow(result, temp, arg);
632	}
633
634	namespace detail{
635
636	template <class T>
637	void small_sinh_series(T x, T& result)
638	{
639	typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
640	bool neg = eval_get_sign(x) < 0;
641	if(neg)
642	x.negate();
643	T p(x);
644	T mult(x);
645	eval_multiply(mult, x);
646	result = x;
647	ui_type k = 1;
648
649	T lim(x);
650	eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
651
652	do
653	{
654	eval_multiply(p, mult);
655	eval_divide(p, ++k);
656	eval_divide(p, ++k);
657	eval_add(result, p);
658	}while(p.compare(lim) >= 0);
659	if(neg)
660	result.negate();
661	}
662
663	template <class T>
664	void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
665	{
666	typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
667	typedef typename mpl::front<typename T::float_types>::type fp_type;
668
669	switch(eval_fpclassify(x))
670	{
671	case FP_NAN:
672	case FP_INFINITE:
673	if(p_sinh)
674	*p_sinh = x;
675	if(p_cosh)
676	{
677	*p_cosh = x;
678	if(eval_get_sign(x) < 0)
679	p_cosh->negate();
680	}
681	return;
682	case FP_ZERO:
683	if(p_sinh)
684	*p_sinh = x;
685	if(p_cosh)
686	*p_cosh = ui_type(1);
687	return;
688	default: ;
689	}
690
691	bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;
692
693	if(p_cosh \|\| !small_sinh)
694	{
695	T e_px, e_mx;
696	eval_exp(e_px, x);
697	eval_divide(e_mx, ui_type(1), e_px);
698
699	if(p_sinh)
700	{
701	if(small_sinh)
702	{
703	small_sinh_series(x, *p_sinh);
704	}
705	else
706	{
707	eval_subtract(*p_sinh, e_px, e_mx);
708	eval_ldexp(p_sinh, p_sinh, -1);
709	}
710	}
711	if(p_cosh)
712	{
713	eval_add(*p_cosh, e_px, e_mx);
714	eval_ldexp(p_cosh, p_cosh, -1);
715	}
716	}
717	else
718	{
719	small_sinh_series(x, *p_sinh);
720	}
721	}
722
723	} // namespace detail
724
725	template <class T>
726	inline void eval_sinh(T& result, const T& x)
727	{
728	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types.");
729	detail::sinhcosh(x, &result, static_cast<T*>(0));
730	}
731
732	template <class T>
733	inline void eval_cosh(T& result, const T& x)
734	{
735	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types.");
736	detail::sinhcosh(x, static_cast<T*>(0), &result);
737	}
738
739	template <class T>
740	inline void eval_tanh(T& result, const T& x)
741	{
742	BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
743	T c;
744	detail::sinhcosh(x, &result, &c);
745	eval_divide(result, c);
746	}
747
748	#ifdef BOOST_MSVC
749	#pragma warning(pop)
750	#endif