code:
// Copyright Christopher Kormanyos 2002 - 2011.
// Copyright 2011 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

template <class T>
void hyp0F1(T& result, const T& b, const T& x)
{
   typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;

   // Compute the series representation of Hypergeometric0F1 taken from
   // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/06/01/01/
   // There are no checks on input range or parameter boundaries.

   T x_pow_n_div_n_fact(x);
   T pochham_b         (b);
   T bp                (b);

   eval_divide(result, x_pow_n_div_n_fact, pochham_b);
   eval_add(result, ui_type(1));

   si_type n;

   T tol;
   tol = ui_type(1);
   eval_ldexp(tol, tol, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
   eval_multiply(tol, result);
   if(eval_get_sign(tol) < 0)
      tol.negate();
   T term;

   const int 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_0f1(; b; 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(bp);
      eval_multiply(pochham_b, bp);

      eval_divide(term, x_pow_n_div_n_fact, pochham_b);
      eval_add(result, term);

      bool neg_term = eval_get_sign(term) < 0;
      if(neg_term)
         term.negate();
      if(term.compare(tol) <= 0)
         break;
   }

   if(n >= series_limit)
      BOOST_THROW_EXCEPTION(std::runtime_error("H0F1 Failed to Converge"));
}

template <class T>
void eval_sin(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sin function is only valid for floating point types.");
   if(&result == &x)
   {
      T temp;
      eval_sin(temp, x);
      result = temp;
      return;
   }

   typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
   typedef typename mpl::front<typename T::float_types>::type fp_type;

   switch(eval_fpclassify(x))
   {
   case FP_INFINITE:
   case FP_NAN:
      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 is undefined or complex and there is no NaN for this number type."));
      return;
   case FP_ZERO:
      result = ui_type(0);
      return;
   default: ;
   }

   // Local copy of the argument
   T xx = x;

   // Analyze and prepare the phase of the argument.
   // Make a local, positive copy of the argument, xx.
   // The argument xx will be reduced to 0 <= xx <= pi/2.
   bool b_negate_sin = false;

   if(eval_get_sign(x) < 0)
   {
      xx.negate();
      b_negate_sin = !b_negate_sin;
   }

   T n_pi, t;
   // Remove even multiples of pi.
   if(xx.compare(get_constant_pi<T>()) > 0)
   {
      eval_divide(n_pi, xx, get_constant_pi<T>());
      eval_trunc(n_pi, n_pi);
      t = ui_type(2);
      eval_fmod(t, n_pi, t);
      const bool b_n_pi_is_even = eval_get_sign(t) == 0;
      eval_multiply(n_pi, get_constant_pi<T>());
      eval_subtract(xx, n_pi);

      BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
      BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));

      // Adjust signs if the multiple of pi is not even.
      if(!b_n_pi_is_even)
      {
         b_negate_sin = !b_negate_sin;
      }
   }

   // Reduce the argument to 0 <= xx <= pi/2.
   eval_ldexp(t, get_constant_pi<T>(), -1);
   if(xx.compare(t) > 0)
   {
      eval_subtract(xx, get_constant_pi<T>(), xx);
      BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
   }

   eval_subtract(t, xx);
   const bool b_zero    = eval_get_sign(xx) == 0;
   const bool b_pi_half = eval_get_sign(t) == 0;

   // Check if the reduced argument is very close to 0 or pi/2.
   const bool    b_near_zero    = xx.compare(fp_type(1e-1)) < 0;
   const bool    b_near_pi_half = t.compare(fp_type(1e-1)) < 0;;

   if(b_zero)
   {
      result = ui_type(0);
   }
   else if(b_pi_half)
   {
      result = ui_type(1);
   }
   else if(b_near_zero)
   {
      eval_multiply(t, xx, xx);
      eval_divide(t, si_type(-4));
      T t2;
      t2 = fp_type(1.5);
      hyp0F1(result, t2, t);
      BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
      eval_multiply(result, xx);
   }
   else if(b_near_pi_half)
   {
      eval_multiply(t, t);
      eval_divide(t, si_type(-4));
      T t2;
      t2 = fp_type(0.5);
      hyp0F1(result, t2, t);
      BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
   }
   else
   {
      // Scale to a small argument for an efficient Taylor series,
      // implemented as a hypergeometric function. Use a standard
      // divide by three identity a certain number of times.
      // Here we use division by 3^9 --> (19683 = 3^9).

      static const si_type n_scale = 9;
      static const si_type n_three_pow_scale = static_cast<si_type>(19683L);

      eval_divide(xx, n_three_pow_scale);

      // Now with small arguments, we are ready for a series expansion.
      eval_multiply(t, xx, xx);
      eval_divide(t, si_type(-4));
      T t2;
      t2 = fp_type(1.5);
      hyp0F1(result, t2, t);
      BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
      eval_multiply(result, xx);

      // Convert back using multiple angle identity.
      for(boost::int32_t k = static_cast<boost::int32_t>(0); k < n_scale; k++)
      {
         // Rescale the cosine value using the multiple angle identity.
         eval_multiply(t2, result, ui_type(3));
         eval_multiply(t, result, result);
         eval_multiply(t, result);
         eval_multiply(t, ui_type(4));
         eval_subtract(result, t2, t);
      }
   }

   if(b_negate_sin)
      result.negate();
}

template <class T>
void eval_cos(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cos function is only valid for floating point types.");
   if(&result == &x)
   {
      T temp;
      eval_cos(temp, x);
      result = temp;
      return;
   }

   typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
   typedef typename mpl::front<typename T::float_types>::type fp_type;

   switch(eval_fpclassify(x))
   {
   case FP_INFINITE:
   case FP_NAN:
      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 is undefined or complex and there is no NaN for this number type."));
      return;
   case FP_ZERO:
      result = ui_type(1);
      return;
   default: ;
   }

   // Local copy of the argument
   T xx = x;

   // Analyze and prepare the phase of the argument.
   // Make a local, positive copy of the argument, xx.
   // The argument xx will be reduced to 0 <= xx <= pi/2.
   bool b_negate_cos = false;

   if(eval_get_sign(x) < 0)
   {
      xx.negate();
   }

   T n_pi, t;
   // Remove even multiples of pi.
   if(xx.compare(get_constant_pi<T>()) > 0)
   {
      eval_divide(t, xx, get_constant_pi<T>());
      eval_trunc(n_pi, t);
      BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
      eval_multiply(t, n_pi, get_constant_pi<T>());
      BOOST_MATH_INSTRUMENT_CODE(t.str(0, std::ios_base::scientific));
      eval_subtract(xx, t);
      BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));

      // Adjust signs if the multiple of pi is not even.
      t = ui_type(2);
      eval_fmod(t, n_pi, t);
      const bool b_n_pi_is_even = eval_get_sign(t) == 0;

      if(!b_n_pi_is_even)
      {
         b_negate_cos = !b_negate_cos;
      }
   }

   // Reduce the argument to 0 <= xx <= pi/2.
   eval_ldexp(t, get_constant_pi<T>(), -1);
   int com = xx.compare(t);
   if(com > 0)
   {
      eval_subtract(xx, get_constant_pi<T>(), xx);
      b_negate_cos = !b_negate_cos;
      BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
   }

   const bool b_zero    = eval_get_sign(xx) == 0;
   const bool b_pi_half = com == 0;

   // Check if the reduced argument is very close to 0.
   const bool    b_near_zero    = xx.compare(fp_type(1e-1)) < 0;

   if(b_zero)
   {
      result = si_type(1);
   }
   else if(b_pi_half)
   {
      result = si_type(0);
   }
   else if(b_near_zero)
   {
      eval_multiply(t, xx, xx);
      eval_divide(t, si_type(-4));
      n_pi = fp_type(0.5f);
      hyp0F1(result, n_pi, t);
      BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
   }
   else
   {
      eval_subtract(t, xx);
      eval_sin(result, t);
   }
   if(b_negate_cos)
      result.negate();
}

template <class T>
void eval_tan(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tan function is only valid for floating point types.");
   if(&result == &x)
   {
      T temp;
      eval_tan(temp, x);
      result = temp;
      return;
   }
   T t;
   eval_sin(result, x);
   eval_cos(t, x);
   eval_divide(result, t);
}

template <class T>
void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x)
{
  // Compute the series representation of hyperg_2f1 taken from
  // Abramowitz and Stegun 15.1.1.
  // There are no checks on input range or parameter boundaries.

   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;

   T x_pow_n_div_n_fact(x);
   T pochham_a         (a);
   T pochham_b         (b);
   T pochham_c         (c);
   T ap                (a);
   T bp                (b);
   T cp                (c);

   eval_multiply(result, pochham_a, pochham_b);
   eval_divide(result, pochham_c);
   eval_multiply(result, x_pow_n_div_n_fact);
   eval_add(result, ui_type(1));

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

   if(eval_get_sign(lim) < 0)
      lim.negate();

   ui_type n;
   T term;

   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_2f1(a, b; c; 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_increment(bp);
      eval_multiply(pochham_b, bp);
      eval_increment(cp);
      eval_multiply(pochham_c, cp);

      eval_multiply(term, pochham_a, pochham_b);
      eval_divide(term, pochham_c);
      eval_multiply(term, x_pow_n_div_n_fact);
      eval_add(result, 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("H2F1 failed to converge."));
}

template <class T>
void eval_asin(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The asin function is only valid for floating point types.");
   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
   typedef typename mpl::front<typename T::float_types>::type fp_type;

   if(&result == &x)
   {
      T t(x);
      eval_asin(result, t);
      return;
   }

   switch(eval_fpclassify(x))
   {
   case FP_NAN:
   case FP_INFINITE:
      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 is undefined or complex and there is no NaN for this number type."));
      return;
   case FP_ZERO:
      result = ui_type(0);
      return;
   default: ;
   }

   const bool b_neg = eval_get_sign(x) < 0;

   T xx(x);
   if(b_neg)
      xx.negate();

   int c = xx.compare(ui_type(1));
   if(c > 0)
   {
      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 is undefined or complex and there is no NaN for this number type."));
      return;
   }
   else if(c == 0)
   {
      result = get_constant_pi<T>();
      eval_ldexp(result, result, -1);
      if(b_neg)
         result.negate();
      return;
   }

   if(xx.compare(fp_type(1e-4)) < 0)
   {
      // http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
      eval_multiply(xx, xx);
      T t1, t2;
      t1 = fp_type(0.5f);
      t2 = fp_type(1.5f);
      hyp2F1(result, t1, t1, t2, xx);
      eval_multiply(result, x);
      return;
   }
   else if(xx.compare(fp_type(1 - 1e-4f)) > 0)
   {
      T dx1;
      T t1, t2;
      eval_subtract(dx1, ui_type(1), xx);
      t1 = fp_type(0.5f);
      t2 = fp_type(1.5f);
      eval_ldexp(dx1, dx1, -1);
      hyp2F1(result, t1, t1, t2, dx1);
      eval_ldexp(dx1, dx1, 2);
      eval_sqrt(t1, dx1);
      eval_multiply(result, t1);
      eval_ldexp(t1, get_constant_pi<T>(), -1);
      result.negate();
      eval_add(result, t1);
      if(b_neg)
         result.negate();
      return;
   }
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
   typedef typename boost::multiprecision::detail::canonical<long double, T>::type guess_type;
#else
   typedef fp_type guess_type;
#endif
   // Get initial estimate using standard math function asin.
   guess_type dd;
   eval_convert_to(&dd, xx);

   result = (guess_type)(std::asin(dd));

   // Newton-Raphson iteration, we should double our precision with each iteration,
   // in practice this seems to not quite work in all cases... so terminate when we
   // have at least 2/3 of the digits correct on the assumption that the correction
   // we've just added will finish the job...

   boost::intmax_t current_precision = eval_ilogb(result);
   boost::intmax_t target_precision = current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3;

   // Newton-Raphson iteration
   while(current_precision > target_precision)
   {
      T sine, cosine;
      eval_sin(sine, result);
      eval_cos(cosine, result);
      eval_subtract(sine, xx);
      eval_divide(sine, cosine);
      eval_subtract(result, sine);
      current_precision = eval_ilogb(sine);
#ifdef FP_ILOGB0
      if(current_precision == FP_ILOGB0)
         break;
#endif
   }
   if(b_neg)
      result.negate();
}

template <class T>
inline void eval_acos(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;

   switch(eval_fpclassify(x))
   {
   case FP_NAN:
   case FP_INFINITE:
      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 is undefined or complex and there is no NaN for this number type."));
      return;
   case FP_ZERO:
      result = get_constant_pi<T>();
      eval_ldexp(result, result, -1); // divide by two.
      return;
   }

   eval_abs(result, x);
   int c = result.compare(ui_type(1));

   if(c > 0)
   {
      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 is undefined or complex and there is no NaN for this number type."));
      return;
   }
   else if(c == 0)
   {
      if(eval_get_sign(x) < 0)
         result = get_constant_pi<T>();
      else
         result = ui_type(0);
      return;
   }

   eval_asin(result, x);
   T t;
   eval_ldexp(t, get_constant_pi<T>(), -1);
   eval_subtract(result, t);
   result.negate();
}

template <class T>
void eval_atan(T& result, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
   typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
   typedef typename mpl::front<typename T::float_types>::type fp_type;

   switch(eval_fpclassify(x))
   {
   case FP_NAN:
      result = x;
      return;
   case FP_ZERO:
      result = ui_type(0);
      return;
   case FP_INFINITE:
      if(eval_get_sign(x) < 0)
      {
         eval_ldexp(result, get_constant_pi<T>(), -1);
         result.negate();
      }
      else
         eval_ldexp(result, get_constant_pi<T>(), -1);
      return;
   default: ;
   }

   const bool b_neg = eval_get_sign(x) < 0;

   T xx(x);
   if(b_neg)
      xx.negate();

   if(xx.compare(fp_type(0.1)) < 0)
   {
      T t1, t2, t3;
      t1 = ui_type(1);
      t2 = fp_type(0.5f);
      t3 = fp_type(1.5f);
      eval_multiply(xx, xx);
      xx.negate();
      hyp2F1(result, t1, t2, t3, xx);
      eval_multiply(result, x);
      return;
   }

   if(xx.compare(fp_type(10)) > 0)
   {
      T t1, t2, t3;
      t1 = fp_type(0.5f);
      t2 = ui_type(1u);
      t3 = fp_type(1.5f);
      eval_multiply(xx, xx);
      eval_divide(xx, si_type(-1), xx);
      hyp2F1(result, t1, t2, t3, xx);
      eval_divide(result, x);
      if(!b_neg)
         result.negate();
      eval_ldexp(t1, get_constant_pi<T>(), -1);
      eval_add(result, t1);
      if(b_neg)
         result.negate();
      return;
   }

   // Get initial estimate using standard math function atan.
   fp_type d;
   eval_convert_to(&d, xx);
   result = fp_type(std::atan(d));

   // Newton-Raphson iteration, we should double our precision with each iteration,
   // in practice this seems to not quite work in all cases... so terminate when we
   // have at least 2/3 of the digits correct on the assumption that the correction
   // we've just added will finish the job...

   boost::intmax_t current_precision = eval_ilogb(result);
   boost::intmax_t target_precision = current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3;

   T s, c, t;
   while(current_precision > target_precision)
   {
      eval_sin(s, result);
      eval_cos(c, result);
      eval_multiply(t, xx, c);
      eval_subtract(t, s);
      eval_multiply(s, t, c);
      eval_add(result, s);
      current_precision = eval_ilogb(s);
#ifdef FP_ILOGB0
      if(current_precision == FP_ILOGB0)
         break;
#endif
   }
   if(b_neg)
      result.negate();
}

template <class T>
void eval_atan2(T& result, const T& y, const T& x)
{
   BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
   if(&result == &y)
   {
      T temp(y);
      eval_atan2(result, temp, x);
      return;
   }
   else if(&result == &x)
   {
      T temp(x);
      eval_atan2(result, y, temp);
      return;
   }

   typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;

   switch(eval_fpclassify(y))
   {
   case FP_NAN:
      result = y;
      return;
   case FP_ZERO:
      {
         int c = eval_get_sign(x);
         if(c < 0)
            result = get_constant_pi<T>();
         else if(c >= 0)
            result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
         return;
      }
   case FP_INFINITE:
      {
         if(eval_fpclassify(x) == FP_INFINITE)
         {
            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 is undefined or complex and there is no NaN for this number type."));
         }
         else
         {
            eval_ldexp(result, get_constant_pi<T>(), -1);
            if(eval_get_sign(y) < 0)
               result.negate();
         }
         return;
      }
   }

   switch(eval_fpclassify(x))
   {
   case FP_NAN:
      result = x;
      return;
   case FP_ZERO:
      {
         eval_ldexp(result, get_constant_pi<T>(), -1);
         if(eval_get_sign(y) < 0)
            result.negate();
         return;
      }
   case FP_INFINITE:
      if(eval_get_sign(x) > 0)
         result = ui_type(0);
      else
         result = get_constant_pi<T>();
      if(eval_get_sign(y) < 0)
         result.negate();
      return;
   }

   T xx;
   eval_divide(xx, y, x);
   if(eval_get_sign(xx) < 0)
      xx.negate();

   eval_atan(result, xx);

   // Determine quadrant (sign) based on signs of x, y
   const bool y_neg = eval_get_sign(y) < 0;
   const bool x_neg = eval_get_sign(x) < 0;

   if(y_neg != x_neg)
      result.negate();

   if(x_neg)
   {
      if(y_neg)
         eval_subtract(result, get_constant_pi<T>());
      else
         eval_add(result, get_constant_pi<T>());
   }
}
template<class T, class A>
inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& 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 = a;
   eval_atan2(result, x, c);
}

template<class T, class A>
inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(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_atan2(result, c, a);
}

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif


code:
///////////////////////////////////////////////////////////////////////////////
//  Copyright 2011 John Maddock. 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)

#ifndef BOOST_MP_GENERIC_INTERCONVERT_HPP
#define BOOST_MP_GENERIC_INTERCONVERT_HPP

#include <boost/multiprecision/detail/default_ops.hpp>

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4127 6326)
#endif

namespace boost{ namespace multiprecision{ namespace detail{

template <class To, class From>
inline To do_cast(const From & from)
{
   return static_cast<To>(from);
}
template <class To, class B, ::boost::multiprecision::expression_template_option et>
inline To do_cast(const number<B, et>& from)
{
   return from.template convert_to<To>();
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
{
   using default_ops::eval_get_sign;
   using default_ops::eval_bitwise_and;
   using default_ops::eval_convert_to;
   using default_ops::eval_right_shift;
   using default_ops::eval_ldexp;
   using default_ops::eval_add;
   using default_ops::eval_is_zero;
   // smallest unsigned type handled natively by "From" is likely to be it's limb_type:
   typedef typename canonical<unsigned char, From>::type   l_limb_type;
   // get the corresponding type that we can assign to "To":
   typedef typename canonical<l_limb_type, To>::type         to_type;
   From t(from);
   bool is_neg = eval_get_sign(t) < 0;
   if(is_neg)
      t.negate();
   // Pick off the first limb:
   l_limb_type limb;
   l_limb_type mask = static_cast<l_limb_type>(~static_cast<l_limb_type>(0));
   From fl;
   eval_bitwise_and(fl, t, mask);
   eval_convert_to(&limb, fl);
   to = static_cast<to_type>(limb);
   eval_right_shift(t, std::numeric_limits<l_limb_type>::digits);
   //
   // Then keep picking off more limbs until "t" is zero:
   //
   To l;
   unsigned shift = std::numeric_limits<l_limb_type>::digits;
   while(!eval_is_zero(t))
   {
      eval_bitwise_and(fl, t, mask);
      eval_convert_to(&limb, fl);
      l = static_cast<to_type>(limb);
      eval_right_shift(t, std::numeric_limits<l_limb_type>::digits);
      eval_ldexp(l, l, shift);
      eval_add(to, l);
      shift += std::numeric_limits<l_limb_type>::digits;
   }
   //
   // Finish off by setting the sign:
   //
   if(is_neg)
      to.negate();
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
{
   using default_ops::eval_get_sign;
   using default_ops::eval_bitwise_and;
   using default_ops::eval_convert_to;
   using default_ops::eval_right_shift;
   using default_ops::eval_left_shift;
   using default_ops::eval_bitwise_or;
   using default_ops::eval_is_zero;
   // smallest unsigned type handled natively by "From" is likely to be it's limb_type:
   typedef typename canonical<unsigned char, From>::type   limb_type;
   // get the corresponding type that we can assign to "To":
   typedef typename canonical<limb_type, To>::type         to_type;
   From t(from);
   bool is_neg = eval_get_sign(t) < 0;
   if(is_neg)
      t.negate();
   // Pick off the first limb:
   limb_type limb;
   limb_type mask = static_cast<limb_type>(~static_cast<limb_type>(0));
   From fl;
   eval_bitwise_and(fl, t, mask);
   eval_convert_to(&limb, fl);
   to = static_cast<to_type>(limb);
   eval_right_shift(t, std::numeric_limits<limb_type>::digits);
   //
   // Then keep picking off more limbs until "t" is zero:
   //
   To l;
   unsigned shift = std::numeric_limits<limb_type>::digits;
   while(!eval_is_zero(t))
   {
      eval_bitwise_and(fl, t, mask);
      eval_convert_to(&limb, fl);
      l = static_cast<to_type>(limb);
      eval_right_shift(t, std::numeric_limits<limb_type>::digits);
      eval_left_shift(l, shift);
      eval_bitwise_or(to, l);
      shift += std::numeric_limits<limb_type>::digits;
   }
   //
   // Finish off by setting the sign:
   //
   if(is_neg)
      to.negate();
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/)
{
#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4127)
#endif
   //
   // The code here only works when the radix of "From" is 2, we could try shifting by other
   // radixes but it would complicate things.... use a string conversion when the radix is other
   // than 2:
   //
   if(std::numeric_limits<number<From> >::radix != 2)
   {
      to = from.str(0, std::ios_base::fmtflags()).c_str();
      return;
   }

   typedef typename canonical<unsigned char, To>::type ui_type;

   using default_ops::eval_fpclassify;
   using default_ops::eval_add;
   using default_ops::eval_subtract;
   using default_ops::eval_convert_to;
   using default_ops::eval_get_sign;
   using default_ops::eval_is_zero;

   //
   // First classify the input, then handle the special cases:
   //
   int c = eval_fpclassify(from);

   if(c == (int)FP_ZERO)
   {
      to = ui_type(0);
      return;
   }
   else if(c == (int)FP_NAN)
   {
      to = static_cast<const char*>("nan");
      return;
   }
   else if(c == (int)FP_INFINITE)
   {
      to = static_cast<const char*>("inf");
      if(eval_get_sign(from) < 0)
         to.negate();
      return;
   }

   typename From::exponent_type e;
   From f, term;
   to = ui_type(0);

   eval_frexp(f, from, &e);

   static const int shift = std::numeric_limits<boost::intmax_t>::digits - 1;

   while(!eval_is_zero(f))
   {
      // extract int sized bits from f:
      eval_ldexp(f, f, shift);
      eval_floor(term, f);
      e -= shift;
      eval_ldexp(to, to, shift);
      typename boost::multiprecision::detail::canonical<boost::intmax_t, To>::type ll;
      eval_convert_to(&ll, term);
      eval_add(to, ll);
      eval_subtract(f, term);
   }
   typedef typename To::exponent_type to_exponent;
   if((e > (std::numeric_limits<to_exponent>::max)()) || (e < (std::numeric_limits<to_exponent>::min)()))
   {
      to = static_cast<const char*>("inf");
      if(eval_get_sign(from) < 0)
         to.negate();
      return;
   }
   eval_ldexp(to, to, static_cast<to_exponent>(e));
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
{
   typedef typename component_type<number<To> >::type     to_component_type;

   number<From> t(from);
   to_component_type n(numerator(t)), d(denominator(t));
   using default_ops::assign_components;
   assign_components(to, n.backend(), d.backend());
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
{
   typedef typename component_type<number<To> >::type     to_component_type;

   number<From> t(from);
   to_component_type n(t), d(1);
   using default_ops::assign_components;
   assign_components(to, n.backend(), d.backend());
}

template <class R, class LargeInteger>
R safe_convert_to_float(const LargeInteger& i)
{
   using std::ldexp;
   if(!i)
      return R(0);
   if(std::numeric_limits<R>::is_specialized && std::numeric_limits<R>::max_exponent)
   {
      LargeInteger val(i);
      if(val.sign() < 0)
         val = -val;
      unsigned mb = msb(val);
      if(mb >= std::numeric_limits<R>::max_exponent)
      {
         int scale_factor = (int)mb + 1 - std::numeric_limits<R>::max_exponent;
         BOOST_ASSERT(scale_factor >= 1);
         val >>= scale_factor;
         R result = val.template convert_to<R>();
         if(std::numeric_limits<R>::digits == 0 || std::numeric_limits<R>::digits >= std::numeric_limits<R>::max_exponent)
         {
            //
            // Calculate and add on the remainder, only if there are more
            // digits in the mantissa that the size of the exponent, in
            // other words if we are dropping digits in the conversion
            // otherwise:
            //
            LargeInteger remainder(i);
            remainder &= (LargeInteger(1) << scale_factor) - 1;
            result += ldexp(safe_convert_to_float<R>(remainder), -scale_factor);
         }
         return i.sign() < 0 ? static_cast<R>(-result) : result;
      }
   }
   return i.template convert_to<R>();
}

template <class To, class Integer>
inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
   generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
{
   //
   // If we get here, then there's something about one type or the other
   // that prevents an exactly rounded result from being calculated
   // (or at least it's not clear how to implement such a thing).
   //
   using default_ops::eval_divide;
   number<To> fn(safe_convert_to_float<number<To> >(n)), fd(safe_convert_to_float<number<To> >(d));
   eval_divide(result, fn.backend(), fd.backend());
}
template <class To, class Integer>
inline typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
   generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
{
   //
   // If we get here, then there's something about one type or the other
   // that prevents an exactly rounded result from being calculated
   // (or at least it's not clear how to implement such a thing).
   //
   To fd(safe_convert_to_float<To>(d));
   result = safe_convert_to_float<To>(n);
   result /= fd;
}

template <class To, class Integer>
typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
   generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_&)
{
   //
   // If we get here, then the precision of type To is known, and the integer type is unbounded
   // so we can use integer division plus manipulation of the remainder to get an exactly
   // rounded result.
   //
   if(num == 0)
   {
      result = 0;
      return;
   }
   bool s = false;
   if(num < 0)
   {
      s = true;
      num = -num;
   }
   int denom_bits = msb(denom);
   int shift = std::numeric_limits<To>::digits + denom_bits - msb(num);
   if(shift > 0)
      num <<= shift;
   else if(shift < 0)
      denom <<= boost::multiprecision::detail::unsigned_abs(shift);
   Integer q, r;
   divide_qr(num, denom, q, r);
   int q_bits = msb(q);
   if(q_bits == std::numeric_limits<To>::digits - 1)
   {
      //
      // Round up if 2 * r > denom:
      //
      r <<= 1;
      int c = r.compare(denom);
      if(c > 0)
         ++q;
      else if((c == 0) && (q & 1u))
      {
         ++q;
      }
   }
   else
   {
      BOOST_ASSERT(q_bits == std::numeric_limits<To>::digits);
      //
      // We basically already have the rounding info:
      //
      if(q & 1u)
      {
         if(r || (q & 2u))
            ++q;
      }
   }
   using std::ldexp;
   result = do_cast<To>(q);
   result = ldexp(result, -shift);
   if(s)
      result = -result;
}
template <class To, class Integer>
inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
   generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_& tag)
{
   number<To> t;
   generic_convert_rational_to_float_imp(t, num, denom, tag);
   result = t.backend();
}

template <class To, class From>
inline void generic_convert_rational_to_float(To& result, const From& f)
{
   //
   // Type From is always a Backend to number<>, or an
   // instance of number<>, but we allow
   // To to be either a Backend type, or a real number type,
   // that way we can call this from generic conversions, and
   // from specific conversions to built in types.
   //
   typedef typename mpl::if_c<is_number<From>::value, From, number<From> >::type actual_from_type;
   typedef typename mpl::if_c<is_number<To>::value || is_floating_point<To>::value, To, number<To> >::type actual_to_type;
   typedef typename component_type<actual_from_type>::type integer_type;
   typedef mpl::bool_<!std::numeric_limits<integer_type>::is_specialized
                      || std::numeric_limits<integer_type>::is_bounded
                      || !std::numeric_limits<actual_to_type>::is_specialized
                      || !std::numeric_limits<actual_to_type>::is_bounded
                      || (std::numeric_limits<actual_to_type>::radix != 2)> dispatch_tag;

   integer_type n(numerator(static_cast<actual_from_type>(f))), d(denominator(static_cast<actual_from_type>(f)));
   generic_convert_rational_to_float_imp(result, n, d, dispatch_tag());
}

template <class To, class From>
inline void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
{
   generic_convert_rational_to_float(to, from);
}

template <class To, class From>
void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<2>& /*radix*/)
{
   typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
   static const int shift = std::numeric_limits<boost::long_long_type>::digits;
   typename From::exponent_type e;
   typename component_type<number<To> >::type num, denom;
   number<From> val(from);
   val = frexp(val, &e);
   while(val)
   {
      val = ldexp(val, shift);
      e -= shift;
      boost::long_long_type ll = boost::math::lltrunc(val);
      val -= ll;
      num <<= shift;
      num += ll;
   }
   denom = ui_type(1u);
   if(e < 0)
      denom <<= -e;
   else if(e > 0)
      num <<= e;
   assign_components(to, num.backend(), denom.backend());
}

template <class To, class From, int Radix>
void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
{
   //
   // This is almost the same as the binary case above, but we have to use
   // scalbn and ilogb rather than ldexp and frexp, we also only extract
   // one Radix digit at a time which is terribly inefficient!
   //
   typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
   typename From::exponent_type e;
   typename component_type<number<To> >::type num, denom;
   number<From> val(from);
   e = ilogb(val);
   val = scalbn(val, -e);
   while(val)
   {
      boost::long_long_type ll = boost::math::lltrunc(val);
      val -= ll;
      val = scalbn(val, 1);
      num *= Radix;
      num += ll;
      --e;
   }
   ++e;
   denom = ui_type(Radix);
   denom = pow(denom, abs(e));
   if(e > 0)
   {
      num *= denom;
      denom = 1;
   }
   assign_components(to, num.backend(), denom.backend());
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/)
{
   generic_interconvert_float2rational(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
{
   number<From> t(from);
   number<To> result(numerator(t) / denominator(t));
   to = result.backend();
}

template <class To, class From>
void generic_interconvert_float2int(To& to, const From& from, const mpl::int_<2>& /*radix*/)
{
   typedef typename From::exponent_type exponent_type;
   static const exponent_type shift = std::numeric_limits<boost::long_long_type>::digits;
   exponent_type e;
   number<To>   num(0u);
   number<From> val(from);
   val = frexp(val, &e);
   while(e > 0)
   {
      int s = (std::min)(e, shift);
      val = ldexp(val, s);
      e -= s;
      boost::long_long_type ll = boost::math::lltrunc(val);
      val -= ll;
      num <<= s;
      num += ll;
   }
   to = num.backend();
}

template <class To, class From, int Radix>
void generic_interconvert_float2int(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
{
   //
   // This is almost the same as the binary case above, but we have to use
   // scalbn and ilogb rather than ldexp and frexp, we also only extract
   // one Radix digit at a time which is terribly inefficient!
   //
   typename From::exponent_type e;
   number<To> num(0u);
   number<From> val(from);
   e = ilogb(val);
   val = scalbn(val, -e);
   while(e >= 0)
   {
      boost::long_long_type ll = boost::math::lltrunc(val);
      val -= ll;
      val = scalbn(val, 1);
      num *= Radix;
      num += ll;
      --e;
   }
   to = num.backend();
}

template <class To, class From>
void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/)
{
   generic_interconvert_float2int(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
}

}
}
} // namespaces

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

#endif  // BOOST_MP_GENERIC_INTERCONVERT_HPP

license: mit


code:
///////////////////////////////////////////////////////////////
//  Copyright 2012 John Maddock. 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_

#ifndef BOOST_MP_INT_FUNC_HPP
#define BOOST_MP_INT_FUNC_HPP

#include <boost/multiprecision/number.hpp>

namespace boost{ namespace multiprecision{

namespace default_ops
{

template <class Backend>
inline void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r)
{
   eval_divide(q, x, y);
   eval_modulus(r, x, y);
}

template <class Backend, class Integer>
inline Integer eval_integer_modulus(const Backend& x, Integer val)
{
   BOOST_MP_USING_ABS
   using default_ops::eval_modulus;
   using default_ops::eval_convert_to;
   typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type int_type;
   Backend t;
   eval_modulus(t, x, static_cast<int_type>(val));
   Integer result;
   eval_convert_to(&result, t);
   return abs(result);
}

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4127)
#endif

template <class B>
inline void eval_gcd(B& result, const B& a, const B& b)
{
   using default_ops::eval_lsb;
   using default_ops::eval_is_zero;
   using default_ops::eval_get_sign;

   int shift;

   B u(a), v(b);

   int s = eval_get_sign(u);

   /* GCD(0,x) := x */
   if(s < 0)
   {
      u.negate();
   }
   else if(s == 0)
   {
      result = v;
      return;
   }
   s = eval_get_sign(v);
   if(s < 0)
   {
      v.negate();
   }
   else if(s == 0)
   {
      result = u;
      return;
   }

   /* Let shift := lg K, where K is the greatest power of 2
   dividing both u and v. */

   unsigned us = eval_lsb(u);
   unsigned vs = eval_lsb(v);
   shift = (std::min)(us, vs);
   eval_right_shift(u, us);
   eval_right_shift(v, vs);

   do
   {
      /* Now u and v are both odd, so diff(u, v) is even.
      Let u = min(u, v), v = diff(u, v)/2. */
      s = u.compare(v);
      if(s > 0)
         u.swap(v);
      if(s == 0)
         break;
      eval_subtract(v, u);
      vs = eval_lsb(v);
      eval_right_shift(v, vs);
   }
   while(true);

   result = u;
   eval_left_shift(result, shift);
}

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

template <class B>
inline void eval_lcm(B& result, const B& a, const B& b)
{
   typedef typename mpl::front<typename B::unsigned_types>::type ui_type;
   B t;
   eval_gcd(t, a, b);

   if(eval_is_zero(t))
   {
      result = static_cast<ui_type>(0);
   }
   else
   {
      eval_divide(result, a, t);
      eval_multiply(result, b);
   }
   if(eval_get_sign(result) < 0)
      result.negate();
}

}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
   divide_qr(const number<Backend, ExpressionTemplates>& x, const number<Backend, ExpressionTemplates>& y,
   number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
{
   using default_ops::eval_qr;
   eval_qr(x.backend(), y.backend(), q.backend(), r.backend());
}

template <class Backend, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
   divide_qr(const number<Backend, ExpressionTemplates>& x, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& y,
   number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
{
   divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r);
}

template <class tag, class A1, class A2, class A3, class A4, class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
   divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const number<Backend, ExpressionTemplates>& y,
   number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
{
   divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r);
}

template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b, class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
   divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& y,
   number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
{
   divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r);
}

template <class Backend, expression_template_option ExpressionTemplates, class Integer>
inline typename enable_if<mpl::and_<is_integral<Integer>, mpl::bool_<number_category<Backend>::value == number_kind_integer> >, Integer>::type
   integer_modulus(const number<Backend, ExpressionTemplates>& x, Integer val)
{
   using default_ops::eval_integer_modulus;
   return eval_integer_modulus(x.backend(), val);
}

template <class tag, class A1, class A2, class A3, class A4, class Integer>
inline typename enable_if<mpl::and_<is_integral<Integer>, mpl::bool_<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer> >, Integer>::type
   integer_modulus(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, Integer val)
{
   typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type result_type;
   return integer_modulus(result_type(x), val);
}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type
   lsb(const number<Backend, ExpressionTemplates>& x)
{
   using default_ops::eval_lsb;
   return eval_lsb(x.backend());
}

template <class tag, class A1, class A2, class A3, class A4>
inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type
   lsb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
{
   typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
   number_type n(x);
   using default_ops::eval_lsb;
   return eval_lsb(n.backend());
}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type
   msb(const number<Backend, ExpressionTemplates>& x)
{
   using default_ops::eval_msb;
   return eval_msb(x.backend());
}

template <class tag, class A1, class A2, class A3, class A4>
inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type
   msb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
{
   typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
   number_type n(x);
   using default_ops::eval_msb;
   return eval_msb(n.backend());
}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, bool>::type
   bit_test(const number<Backend, ExpressionTemplates>& x, unsigned index)
{
   using default_ops::eval_bit_test;
   return eval_bit_test(x.backend(), index);
}

template <class tag, class A1, class A2, class A3, class A4>
inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, bool>::type
   bit_test(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, unsigned index)
{
   typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
   number_type n(x);
   using default_ops::eval_bit_test;
   return eval_bit_test(n.backend(), index);
}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
   bit_set(number<Backend, ExpressionTemplates>& x, unsigned index)
{
   using default_ops::eval_bit_set;
   eval_bit_set(x.backend(), index);
   return x;
}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
   bit_unset(number<Backend, ExpressionTemplates>& x, unsigned index)
{
   using default_ops::eval_bit_unset;
   eval_bit_unset(x.backend(), index);
   return x;
}

template <class Backend, expression_template_option ExpressionTemplates>
inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
   bit_flip(number<Backend, ExpressionTemplates>& x, unsigned index)
{
   using default_ops::eval_bit_flip;
   eval_bit_flip(x.backend(), index);
   return x;
}

namespace default_ops{

//
// Within powm, we need a type with twice as many digits as the argument type, define
// a traits class to obtain that type:
//
template <class Backend>
struct double_precision_type
{
   typedef Backend type;
};

//
// If the exponent is a signed integer type, then we need to
// check the value is positive:
//
template <class Backend>
inline void check_sign_of_backend(const Backend& v, const mpl::true_)
{
   if(eval_get_sign(v) < 0)
   {
      BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
   }
}
template <class Backend>
inline void check_sign_of_backend(const Backend&, const mpl::false_){}
//
// Calculate (a^p)%c:
//
template <class Backend>
void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c)
{
   using default_ops::eval_bit_test;
   using default_ops::eval_get_sign;
   using default_ops::eval_multiply;
   using default_ops::eval_modulus;
   using default_ops::eval_right_shift;

   typedef typename double_precision_type<Backend>::type double_type;
   typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;

   check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());

   double_type x, y(a), b(p), t;
   x = ui_type(1u);

   while(eval_get_sign(b) > 0)
   {
      if(eval_bit_test(b, 0))
      {
         eval_multiply(t, x, y);
         eval_modulus(x, t, c);
      }
      eval_multiply(t, y, y);
      eval_modulus(y, t, c);
      eval_right_shift(b, ui_type(1));
   }
   Backend x2(x);
   eval_modulus(result, x2, c);
}

template <class Backend, class Integer>
void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c)
{
   typedef typename double_precision_type<Backend>::type double_type;
   typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
   typedef typename boost::multiprecision::detail::canonical<Integer, double_type>::type i1_type;
   typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type i2_type;

   using default_ops::eval_bit_test;
   using default_ops::eval_get_sign;
   using default_ops::eval_multiply;
   using default_ops::eval_modulus;
   using default_ops::eval_right_shift;

   check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());

   if(eval_get_sign(p) < 0)
   {
      BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
   }

   double_type x, y(a), b(p), t;
   x = ui_type(1u);

   while(eval_get_sign(b) > 0)
   {
      if(eval_bit_test(b, 0))
      {
         eval_multiply(t, x, y);
         eval_modulus(x, t, static_cast<i1_type>(c));
      }
      eval_multiply(t, y, y);
      eval_modulus(y, t, static_cast<i1_type>(c));
      eval_right_shift(b, ui_type(1));
   }
   Backend x2(x);
   eval_modulus(result, x2, static_cast<i2_type>(c));
}

template <class Backend, class Integer>
typename enable_if<is_unsigned<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
{
   typedef typename double_precision_type<Backend>::type double_type;
   typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;

   using default_ops::eval_bit_test;
   using default_ops::eval_get_sign;
   using default_ops::eval_multiply;
   using default_ops::eval_modulus;
   using default_ops::eval_right_shift;

   double_type x, y(a), t;
   x = ui_type(1u);

   while(b > 0)
   {
      if(b & 1)
      {
         eval_multiply(t, x, y);
         eval_modulus(x, t, c);
      }
      eval_multiply(t, y, y);
      eval_modulus(y, t, c);
      b >>= 1;
   }
   Backend x2(x);
   eval_modulus(result, x2, c);
}

template <class Backend, class Integer>
typename enable_if<is_signed<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
{
   if(b < 0)
   {
      BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
   }
   eval_powm(result, a, static_cast<typename make_unsigned<Integer>::type>(b), c);
}

template <class Backend, class Integer1, class Integer2>
typename enable_if<is_unsigned<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
{
   typedef typename double_precision_type<Backend>::type double_type;
   typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
   typedef typename boost::multiprecision::detail::canonical<Integer1, double_type>::type i1_type;
   typedef typename boost::multiprecision::detail::canonical<Integer2, Backend>::type i2_type;

   using default_ops::eval_bit_test;
   using default_ops::eval_get_sign;
   using default_ops::eval_multiply;
   using default_ops::eval_modulus;
   using default_ops::eval_right_shift;

   double_type x, y(a), t;
   x = ui_type(1u);

   while(b > 0)
   {
      if(b & 1)
      {
         eval_multiply(t, x, y);
         eval_modulus(x, t, static_cast<i1_type>(c));
      }
      eval_multiply(t, y, y);
      eval_modulus(y, t, static_cast<i1_type>(c));
      b >>= 1;
   }
   Backend x2(x);
   eval_modulus(result, x2, static_cast<i2_type>(c));
}

template <class Backend, class Integer1, class Integer2>
typename enable_if<is_signed<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
{
   if(b < 0)
   {
      BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
   }
   eval_powm(result, a, static_cast<typename make_unsigned<Integer1>::type>(b), c);
}

struct powm_func
{
   template <class T, class U, class V>
   void operator()(T& result, const T& b, const U& p, const V& m)const
   {
      eval_powm(result, b, p, m);
   }
};

}

template <class T, class U, class V>
inline typename enable_if<
   mpl::and_<
      mpl::bool_<number_category<T>::value == number_kind_integer>,
      mpl::or_<
         is_number<T>,
         is_number_expression<T>
      >,
      mpl::or_<
         is_number<U>,
         is_number_expression<U>,
         is_integral<U>
      >,
      mpl::or_<
         is_number<V>,
         is_number_expression<V>,
         is_integral<V>
      >
   >,
   typename mpl::if_<
      is_no_et_number<T>,
      T,
      typename mpl::if_<
         is_no_et_number<U>,
         U,
         typename mpl::if_<
            is_no_et_number<V>,
            V,
            detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type
         >::type
      >::type
   >::type
   powm(const T& b, const U& p, const V& mod)
{
   return detail::expression<detail::function, default_ops::powm_func, T, U, V>(
      default_ops::powm_func(), b, p, mod);
}

}} //namespaces

#endif


code:
///////////////////////////////////////////////////////////////////////////////
//  Copyright 2016 John Maddock. 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)

#ifndef BOOST_MP_MIN_MAX_HPP
#define BOOST_MP_MIN_MAX_HPP

#include <boost/multiprecision/traits/is_backend.hpp>

namespace boost{ namespace multiprecision{

//
// Expression template overloads for (min) and (max):
//
// Introduced in response to https://svn.boost.org/trac/boost/ticket/11149
// note that these can not legally be injected into namespace std, and that doing so
// may break future enhancements to the standard.  None the less adding
// namespace std{ using boost::multiprecision::(min); using boost::multiprecision::(max); }
// to your code may get some generic code working that wouldn't work otherwise.
//
// The use of enable_if on the return type is to avoid poisoning std::min/max,
// otherwise attempting to make an explicit call to min<long>(a, b) when these and std
// versions are in scope, will cause the compiler to try to instantiate the signatures
// for our versions as well as the std ones, which in turn instantiates number<long>
// which fails to compile as "long" is not a valid backend type.
//
template <class Backend>
inline typename boost::enable_if_c < boost::multiprecision::detail::is_backend<Backend>::value, const number<Backend, et_on>&>::type
   (min)(const number<Backend, et_on>& a, const number<Backend, et_on>& b)
{
   return a < b ? a : b;
}
template <class Backend, class tag, class A1, class A2, class A3, class A4>
inline typename boost::enable_if_c < boost::multiprecision::detail::is_backend<Backend>::value, const number<Backend, et_on> >::type
   (min)(const number<Backend, et_on>& a, const detail::expression<tag, A1, A2, A3, A4>& b)
{
   number<Backend, et_on> t(b);
   if(a < t)
      return a;
   return BOOST_MP_MOVE(t);
}
template <class tag, class A1, class A2, class A3, class A4, class Backend>
inline typename boost::enable_if_c < boost::multiprecision::detail::is_backend<Backend>::value, const number<Backend, et_on> >::type
   (min)(const detail::expression<tag, A1, A2, A3, A4>& a, const number<Backend, et_on>& b)
{
   number<Backend, et_on> t(a);
   if(t < b)
      return BOOST_MP_MOVE(t);
   return b;
}
template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b>
inline typename detail::expression<tag, A1, A2, A3, A4>::result_type
      (min)(const detail::expression<tag, A1, A2, A3, A4>& a, const detail::expression<tagb, A1b, A2b, A3b, A4b>& b)
{
   typename detail::expression<tag, A1, A2, A3, A4>::result_type t1(a), t2(b);
   if(t1 < t2)
      return BOOST_MP_MOVE(t1);
   return BOOST_MP_MOVE(t2);
}
template <class tag, class A1, class A2, class A3, class A4>
inline typename detail::expression<tag, A1, A2, A3, A4>::result_type (min)(const detail::expression<tag, A1, A2, A3, A4>& a, const detail::expression<tag, A1, A2, A3, A4>& b)
{
   typename detail::expression<tag, A1, A2, A3, A4>::result_type t1(a), t2(b);
   if(t1 < t2)
      return BOOST_MP_MOVE(t1);
   return BOOST_MP_MOVE(t2);
}

template <class Backend>
inline typename boost::enable_if_c < boost::multiprecision::detail::is_backend<Backend>::value, const number<Backend, et_on>&>::type
   (max)(const number<Backend, et_on>& a, const number<Backend, et_on>& b)
{
   return a > b ? a : b;
}
template <class Backend, class tag, class A1, class A2, class A3, class A4>
inline typename boost::enable_if_c < boost::multiprecision::detail::is_backend<Backend>::value, const number<Backend, et_on> >::type
   (max)(const number<Backend, et_on>& a, const detail::expression<tag, A1, A2, A3, A4>& b)
{
   number<Backend, et_on> t(b);
   if(a > t)
      return a;
   return BOOST_MP_MOVE(t);
}
template <class tag, class A1, class A2, class A3, class A4, class Backend>
inline typename boost::enable_if_c < boost::multiprecision::detail::is_backend<Backend>::value, const number<Backend, et_on> >::type
   (max)(const detail::expression<tag, A1, A2, A3, A4>& a, const number<Backend, et_on>& b)
{
   number<Backend, et_on> t(a);
   if(t > b)
      return BOOST_MP_MOVE(t);
   return b;
}
template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b>
inline typename detail::expression<tag, A1, A2, A3, A4>::result_type
      (max)(const detail::expression<tag, A1, A2, A3, A4>& a, const detail::expression<tagb, A1b, A2b, A3b, A4b>& b)
{
   typename detail::expression<tag, A1, A2, A3, A4>::result_type t1(a), t2(b);
   if(t1 > t2)
      return BOOST_MP_MOVE(t1);
   return BOOST_MP_MOVE(t2);
}
template <class tag, class A1, class A2, class A3, class A4>
inline typename detail::expression<tag, A1, A2, A3, A4>::result_type (max)(const detail::expression<tag, A1, A2, A3, A4>& a, const detail::expression<tag, A1, A2, A3, A4>& b)
{
   typename detail::expression<tag, A1, A2, A3, A4>::result_type t1(a), t2(b);
   if(t1 > t2)
      return BOOST_MP_MOVE(t1);
   return BOOST_MP_MOVE(t2);
}

}} // namespaces

#endif


code:
///////////////////////////////////////////////////////////////////////////////
//  Copyright 2012 John Maddock. 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)

#ifndef BOOST_MP_NO_ET_OPS_HPP
#define BOOST_MP_NO_ET_OPS_HPP

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable: 4714)
#endif

namespace boost{
namespace multiprecision{

//
// Operators for non-expression template enabled number.
// NOTE: this is not a complete header - really just a suffix to default_ops.hpp.
// NOTE: these operators have to be defined after the methods in default_ops.hpp.
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator - (const number<B, et_off>& v)
{
   BOOST_STATIC_ASSERT_MSG(is_signed_number<B>::value, "Negating an unsigned type results in ill-defined behavior.");
   number<B, et_off> result(v);
   result.backend().negate();
   return BOOST_MP_MOVE(result);
}
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator ~ (const number<B, et_off>& v)
{
   number<B, et_off> result;
   eval_complement(result.backend(), v.backend());
   return BOOST_MP_MOVE(result);
}
//
// Addition:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator + (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_add;
   eval_add(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator + (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_add;
   eval_add(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator + (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_add;
   eval_add(result.backend(), b.backend(), number<B, et_off>::canonical_value(a));
   return BOOST_MP_MOVE(result);
}
//
// Subtraction:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator - (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_subtract;
   eval_subtract(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator - (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_subtract;
   eval_subtract(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator - (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_subtract;
   eval_subtract(result.backend(), number<B, et_off>::canonical_value(a), b.backend());
   return BOOST_MP_MOVE(result);
}
//
// Multiply:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator * (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_multiply;
   eval_multiply(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator * (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_multiply;
   eval_multiply(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator * (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_multiply;
   eval_multiply(result.backend(), b.backend(), number<B, et_off>::canonical_value(a));
   return BOOST_MP_MOVE(result);
}
//
// divide:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator / (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_divide;
   eval_divide(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator / (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_divide;
   eval_divide(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator / (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_divide;
   eval_divide(result.backend(), number<B, et_off>::canonical_value(a), b.backend());
   return BOOST_MP_MOVE(result);
}
//
// modulus:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator % (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_modulus;
   eval_modulus(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator % (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_modulus;
   eval_modulus(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator % (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_modulus;
   eval_modulus(result.backend(), number<B, et_off>::canonical_value(a), b.backend());
   return BOOST_MP_MOVE(result);
}
//
// Bitwise or:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator | (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator | (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator | (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(result.backend(), b.backend(), number<B, et_off>::canonical_value(a));
   return BOOST_MP_MOVE(result);
}
//
// Bitwise xor:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator ^ (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator ^ (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator ^ (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(result.backend(), b.backend(), number<B, et_off>::canonical_value(a));
   return BOOST_MP_MOVE(result);
}
//
// Bitwise and:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator & (const number<B, et_off>& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(result.backend(), a.backend(), b.backend());
   return BOOST_MP_MOVE(result);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator & (const number<B, et_off>& a, const V& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(result.backend(), a.backend(), number<B, et_off>::canonical_value(b));
   return BOOST_MP_MOVE(result);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator & (const V& a, const number<B, et_off>& b)
{
   number<B, et_off> result;
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(result.backend(), b.backend(), number<B, et_off>::canonical_value(a));
   return BOOST_MP_MOVE(result);
}
//
// shifts:
//
template <class B, class I>
BOOST_MP_FORCEINLINE typename enable_if_c<is_integral<I>::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator << (const number<B, et_off>& a, const I& b)
{
   number<B, et_off> result(a);
   using default_ops::eval_left_shift;
   detail::check_shift_range(b, mpl::bool_<(sizeof(I) > sizeof(std::size_t))>(), is_signed<I>());
   eval_left_shift(result.backend(), b);
   return BOOST_MP_MOVE(result);
}
template <class B, class I>
BOOST_MP_FORCEINLINE typename enable_if_c<is_integral<I>::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator >> (const number<B, et_off>& a, const I& b)
{
   number<B, et_off> result(a);
   using default_ops::eval_right_shift;
   detail::check_shift_range(b, mpl::bool_<(sizeof(I) > sizeof(std::size_t))>(), is_signed<I>());
   eval_right_shift(result.backend(), b);
   return BOOST_MP_MOVE(result);
}

#if !defined(BOOST_NO_CXX11_RVALUE_REFERENCES) && !(defined(__GNUC__) && ((__GNUC__ == 4) && (__GNUC_MINOR__ < 5)))
//
// If we have rvalue references go all over again with rvalue ref overloads and move semantics.
// Note that while it would be tempting to implement these so they return an rvalue reference
// (and indeed this would be optimally efficient), this is unsafe due to users propensity to
// write:
//
// const T& t = a * b;
//
// which would lead to a dangling reference if we didn't return by value.  Of course move
// semantics help a great deal in return by value, so performance is still pretty good...
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator - (number<B, et_off>&& v)
{
   BOOST_STATIC_ASSERT_MSG(is_signed_number<B>::value, "Negating an unsigned type results in ill-defined behavior.");
   v.backend().negate();
   return static_cast<number<B, et_off>&&>(v);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator ~ (number<B, et_off>&& v)
{
   eval_complement(v.backend(), v.backend());
   return static_cast<number<B, et_off>&&>(v);
}
//
// Addition:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator + (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_add;
   eval_add(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator + (const number<B, et_off>& a, number<B, et_off>&& b)
{
   using default_ops::eval_add;
   eval_add(b.backend(), a.backend());
   return static_cast<number<B, et_off>&&>(b);
}
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator + (number<B, et_off>&& a, number<B, et_off>&& b)
{
   using default_ops::eval_add;
   eval_add(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator + (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_add;
   eval_add(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator + (const V& a, number<B, et_off>&& b)
{
   using default_ops::eval_add;
   eval_add(b.backend(), number<B, et_off>::canonical_value(a));
   return static_cast<number<B, et_off>&&>(b);
}
//
// Subtraction:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator - (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_subtract;
   eval_subtract(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if<is_signed_number<B>, number<B, et_off> >::type operator - (const number<B, et_off>& a, number<B, et_off>&& b)
{
   using default_ops::eval_subtract;
   eval_subtract(b.backend(), a.backend());
   b.backend().negate();
   return static_cast<number<B, et_off>&&>(b);
}
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator - (number<B, et_off>&& a, number<B, et_off>&& b)
{
   using default_ops::eval_subtract;
   eval_subtract(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator - (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_subtract;
   eval_subtract(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<(is_compatible_arithmetic_type<V, number<B, et_off> >::value && is_signed_number<B>::value), number<B, et_off> >::type
   operator - (const V& a, number<B, et_off>&& b)
{
   using default_ops::eval_subtract;
   eval_subtract(b.backend(), number<B, et_off>::canonical_value(a));
   b.backend().negate();
   return static_cast<number<B, et_off>&&>(b);
}
//
// Multiply:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator * (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_multiply;
   eval_multiply(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator * (const number<B, et_off>& a, number<B, et_off>&& b)
{
   using default_ops::eval_multiply;
   eval_multiply(b.backend(), a.backend());
   return static_cast<number<B, et_off>&&>(b);
}
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator * (number<B, et_off>&& a, number<B, et_off>&& b)
{
   using default_ops::eval_multiply;
   eval_multiply(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator * (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_multiply;
   eval_multiply(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator * (const V& a, number<B, et_off>&& b)
{
   using default_ops::eval_multiply;
   eval_multiply(b.backend(), number<B, et_off>::canonical_value(a));
   return static_cast<number<B, et_off>&&>(b);
}
//
// divide:
//
template <class B>
BOOST_MP_FORCEINLINE number<B, et_off> operator / (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_divide;
   eval_divide(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if<is_compatible_arithmetic_type<V, number<B, et_off> >, number<B, et_off> >::type
   operator / (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_divide;
   eval_divide(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
//
// modulus:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator % (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_modulus;
   eval_modulus(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator % (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_modulus;
   eval_modulus(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
//
// Bitwise or:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator | (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator | (const number<B, et_off>& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(b.backend(), a.backend());
   return static_cast<number<B, et_off>&&>(b);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator | (number<B, et_off>&& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator | (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator | (const V& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_or;
   eval_bitwise_or(b.backend(), number<B, et_off>::canonical_value(a));
   return static_cast<number<B, et_off>&&>(b);
}
//
// Bitwise xor:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator ^ (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator ^ (const number<B, et_off>& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(b.backend(), a.backend());
   return static_cast<number<B, et_off>&&>(b);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator ^ (number<B, et_off>&& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator ^ (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator ^ (const V& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_xor;
   eval_bitwise_xor(b.backend(), number<B, et_off>::canonical_value(a));
   return static_cast<number<B, et_off>&&>(b);
}
//
// Bitwise and:
//
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator & (number<B, et_off>&& a, const number<B, et_off>& b)
{
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator & (const number<B, et_off>& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(b.backend(), a.backend());
   return static_cast<number<B, et_off>&&>(b);
}
template <class B>
BOOST_MP_FORCEINLINE typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, et_off> >::type operator & (number<B, et_off>&& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(a.backend(), b.backend());
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class V>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator & (number<B, et_off>&& a, const V& b)
{
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(a.backend(), number<B, et_off>::canonical_value(b));
   return static_cast<number<B, et_off>&&>(a);
}
template <class V, class B>
BOOST_MP_FORCEINLINE typename enable_if_c<is_compatible_arithmetic_type<V, number<B, et_off> >::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator & (const V& a, number<B, et_off>&& b)
{
   using default_ops::eval_bitwise_and;
   eval_bitwise_and(b.backend(), number<B, et_off>::canonical_value(a));
   return static_cast<number<B, et_off>&&>(b);
}
//
// shifts:
//
template <class B, class I>
BOOST_MP_FORCEINLINE typename enable_if_c<is_integral<I>::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator << (number<B, et_off>&& a, const I& b)
{
   using default_ops::eval_left_shift;
   eval_left_shift(a.backend(), b);
   return static_cast<number<B, et_off>&&>(a);
}
template <class B, class I>
BOOST_MP_FORCEINLINE typename enable_if_c<is_integral<I>::value && (number_category<B>::value == number_kind_integer), number<B, et_off> >::type
   operator >> (number<B, et_off>&& a, const I& b)
{
   using default_ops::eval_right_shift;
   eval_right_shift(a.backend(), b);
   return static_cast<number<B, et_off>&&>(a);
}

#endif

}} // namespaces

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

#endif // BOOST_MP_NO_ET_OPS_HPP