Chris@16: ///////////////////////////////////////////////////////////////
Chris@16: //  Copyright 2012 John Maddock. Distributed under the Boost
Chris@16: //  Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
Chris@16: 
Chris@16: #ifndef BOOST_MP_INT_FUNC_HPP
Chris@16: #define BOOST_MP_INT_FUNC_HPP
Chris@16: 
Chris@16: #include <boost/multiprecision/number.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace multiprecision{
Chris@16: 
Chris@16: namespace default_ops
Chris@16: {
Chris@16: 
Chris@16: template <class Backend>
Chris@16: inline void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r)
Chris@16: {
Chris@16:    eval_divide(q, x, y);
Chris@16:    eval_modulus(r, x, y);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, class Integer>
Chris@16: inline Integer eval_integer_modulus(const Backend& x, Integer val)
Chris@16: {
Chris@16:    BOOST_MP_USING_ABS
Chris@16:    using default_ops::eval_modulus;
Chris@16:    using default_ops::eval_convert_to;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type int_type;
Chris@16:    Backend t;
Chris@16:    eval_modulus(t, x, static_cast<int_type>(val));
Chris@16:    Integer result;
Chris@16:    eval_convert_to(&result, t);
Chris@16:    return abs(result);
Chris@16: }
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(push)
Chris@16: #pragma warning(disable:4127)
Chris@16: #endif
Chris@16: 
Chris@16: template <class B>
Chris@16: inline void eval_gcd(B& result, const B& a, const B& b)
Chris@16: {
Chris@16:    using default_ops::eval_lsb;
Chris@16:    using default_ops::eval_is_zero;
Chris@16:    using default_ops::eval_get_sign;
Chris@16: 
Chris@16:    int shift;
Chris@16: 
Chris@16:    B u(a), v(b);
Chris@16: 
Chris@16:    int s = eval_get_sign(u);
Chris@16: 
Chris@16:    /* GCD(0,x) := x */
Chris@16:    if(s < 0)
Chris@16:    {
Chris@16:       u.negate();
Chris@16:    }
Chris@16:    else if(s == 0)
Chris@16:    {
Chris@16:       result = v;
Chris@16:       return;
Chris@16:    }
Chris@16:    s = eval_get_sign(v);
Chris@16:    if(s < 0)
Chris@16:    {
Chris@16:       v.negate();
Chris@16:    }
Chris@16:    else if(s == 0)
Chris@16:    {
Chris@16:       result = u;
Chris@16:       return;
Chris@16:    }
Chris@16: 
Chris@16:    /* Let shift := lg K, where K is the greatest power of 2
Chris@16:    dividing both u and v. */
Chris@16: 
Chris@16:    unsigned us = eval_lsb(u);
Chris@16:    unsigned vs = eval_lsb(v);
Chris@16:    shift = (std::min)(us, vs);
Chris@16:    eval_right_shift(u, us);
Chris@16:    eval_right_shift(v, vs);
Chris@16: 
Chris@16:    do 
Chris@16:    {
Chris@16:       /* Now u and v are both odd, so diff(u, v) is even.
Chris@16:       Let u = min(u, v), v = diff(u, v)/2. */
Chris@16:       s = u.compare(v);
Chris@16:       if(s > 0)
Chris@16:          u.swap(v);
Chris@16:       if(s == 0)
Chris@16:          break;
Chris@16:       eval_subtract(v, u);
Chris@16:       vs = eval_lsb(v);
Chris@16:       eval_right_shift(v, vs);
Chris@16:    } 
Chris@16:    while(true);
Chris@16: 
Chris@16:    result = u;
Chris@16:    eval_left_shift(result, shift);
Chris@16: }
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(pop)
Chris@16: #endif
Chris@16: 
Chris@16: template <class B>
Chris@16: inline void eval_lcm(B& result, const B& a, const B& b)
Chris@16: {
Chris@16:    typedef typename mpl::front<typename B::unsigned_types>::type ui_type;
Chris@16:    B t;
Chris@16:    eval_gcd(t, a, b);
Chris@16: 
Chris@16:    if(eval_is_zero(t))
Chris@16:    {
Chris@16:       result = static_cast<ui_type>(0);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       eval_divide(result, a, t);
Chris@16:       eval_multiply(result, b);
Chris@16:    }
Chris@16:    if(eval_get_sign(result) < 0)
Chris@16:       result.negate();
Chris@16: }
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type 
Chris@16:    divide_qr(const number<Backend, ExpressionTemplates>& x, const number<Backend, ExpressionTemplates>& y,
Chris@16:    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
Chris@16: {
Chris@16:    using default_ops::eval_qr;
Chris@16:    eval_qr(x.backend(), y.backend(), q.backend(), r.backend());
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type 
Chris@16:    divide_qr(const number<Backend, ExpressionTemplates>& x, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& y,
Chris@16:    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
Chris@16: {
Chris@16:    divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r);
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class A1, class A2, class A3, class A4, class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type 
Chris@16:    divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const number<Backend, ExpressionTemplates>& y,
Chris@16:    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
Chris@16: {
Chris@16:    divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r);
Chris@16: }
Chris@16: 
Chris@16: 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>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type 
Chris@16:    divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& y,
Chris@16:    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
Chris@16: {
Chris@16:    divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates, class Integer>
Chris@16: inline typename enable_if<mpl::and_<is_integral<Integer>, mpl::bool_<number_category<Backend>::value == number_kind_integer> >, Integer>::type 
Chris@16:    integer_modulus(const number<Backend, ExpressionTemplates>& x, Integer val)
Chris@16: {
Chris@16:    using default_ops::eval_integer_modulus;
Chris@16:    return eval_integer_modulus(x.backend(), val);
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class A1, class A2, class A3, class A4, class Integer>
Chris@16: 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 
Chris@16:    integer_modulus(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, Integer val)
Chris@16: {
Chris@16:    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type result_type;
Chris@16:    return integer_modulus(result_type(x), val);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type 
Chris@16:    lsb(const number<Backend, ExpressionTemplates>& x)
Chris@16: {
Chris@16:    using default_ops::eval_lsb;
Chris@16:    return eval_lsb(x.backend());
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class A1, class A2, class A3, class A4>
Chris@16: inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type 
Chris@16:    lsb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
Chris@16: {
Chris@16:    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
Chris@16:    number_type n(x);
Chris@16:    using default_ops::eval_lsb;
Chris@16:    return eval_lsb(n.backend());
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type 
Chris@16:    msb(const number<Backend, ExpressionTemplates>& x)
Chris@16: {
Chris@16:    using default_ops::eval_msb;
Chris@16:    return eval_msb(x.backend());
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class A1, class A2, class A3, class A4>
Chris@16: inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type 
Chris@16:    msb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
Chris@16: {
Chris@16:    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
Chris@16:    number_type n(x);
Chris@16:    using default_ops::eval_msb;
Chris@16:    return eval_msb(n.backend());
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, bool>::type 
Chris@16:    bit_test(const number<Backend, ExpressionTemplates>& x, unsigned index)
Chris@16: {
Chris@16:    using default_ops::eval_bit_test;
Chris@16:    return eval_bit_test(x.backend(), index);
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class A1, class A2, class A3, class A4>
Chris@16: inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, bool>::type 
Chris@16:    bit_test(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, unsigned index)
Chris@16: {
Chris@16:    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
Chris@16:    number_type n(x);
Chris@16:    using default_ops::eval_bit_test;
Chris@16:    return eval_bit_test(n.backend(), index);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type 
Chris@16:    bit_set(number<Backend, ExpressionTemplates>& x, unsigned index)
Chris@16: {
Chris@16:    using default_ops::eval_bit_set;
Chris@16:    eval_bit_set(x.backend(), index);
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type 
Chris@16:    bit_unset(number<Backend, ExpressionTemplates>& x, unsigned index)
Chris@16: {
Chris@16:    using default_ops::eval_bit_unset;
Chris@16:    eval_bit_unset(x.backend(), index);
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, expression_template_option ExpressionTemplates>
Chris@16: inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type 
Chris@16:    bit_flip(number<Backend, ExpressionTemplates>& x, unsigned index)
Chris@16: {
Chris@16:    using default_ops::eval_bit_flip;
Chris@16:    eval_bit_flip(x.backend(), index);
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: namespace default_ops{
Chris@16: 
Chris@16: //
Chris@16: // Within powm, we need a type with twice as many digits as the argument type, define
Chris@16: // a traits class to obtain that type:
Chris@16: //
Chris@16: template <class Backend>
Chris@16: struct double_precision_type
Chris@16: {
Chris@16:    typedef Backend type;
Chris@16: };
Chris@16: 
Chris@16: //
Chris@101: // If the exponent is a signed integer type, then we need to
Chris@101: // check the value is positive:
Chris@101: //
Chris@101: template <class Backend>
Chris@101: inline void check_sign_of_backend(const Backend& v, const mpl::true_)
Chris@101: {
Chris@101:    if(eval_get_sign(v) < 0)
Chris@101:    {
Chris@101:       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
Chris@101:    }
Chris@101: }
Chris@101: template <class Backend>
Chris@101: inline void check_sign_of_backend(const Backend&, const mpl::false_){}
Chris@101: //
Chris@16: // Calculate (a^p)%c:
Chris@16: //
Chris@16: template <class Backend>
Chris@16: void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c)
Chris@16: {
Chris@16:    using default_ops::eval_bit_test;
Chris@16:    using default_ops::eval_get_sign;
Chris@16:    using default_ops::eval_multiply;
Chris@16:    using default_ops::eval_modulus;
Chris@16:    using default_ops::eval_right_shift;
Chris@16: 
Chris@16:    typedef typename double_precision_type<Backend>::type double_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
Chris@101: 
Chris@101:    check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());
Chris@16:    
Chris@16:    double_type x, y(a), b(p), t;
Chris@16:    x = ui_type(1u);
Chris@16: 
Chris@16:    while(eval_get_sign(b) > 0)
Chris@16:    {
Chris@16:       if(eval_bit_test(b, 0))
Chris@16:       {
Chris@16:          eval_multiply(t, x, y);
Chris@16:          eval_modulus(x, t, c);
Chris@16:       }
Chris@16:       eval_multiply(t, y, y);
Chris@16:       eval_modulus(y, t, c);
Chris@16:       eval_right_shift(b, ui_type(1));
Chris@16:    }
Chris@16:    Backend x2(x);
Chris@16:    eval_modulus(result, x2, c);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, class Integer>
Chris@16: void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c)
Chris@16: {
Chris@16:    typedef typename double_precision_type<Backend>::type double_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<Integer, double_type>::type i1_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type i2_type;
Chris@16: 
Chris@16:    using default_ops::eval_bit_test;
Chris@16:    using default_ops::eval_get_sign;
Chris@16:    using default_ops::eval_multiply;
Chris@16:    using default_ops::eval_modulus;
Chris@16:    using default_ops::eval_right_shift;
Chris@16: 
Chris@101:    check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());
Chris@101: 
Chris@16:    if(eval_get_sign(p) < 0)
Chris@16:    {
Chris@16:       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
Chris@16:    }
Chris@16: 
Chris@16:    double_type x, y(a), b(p), t;
Chris@16:    x = ui_type(1u);
Chris@16: 
Chris@16:    while(eval_get_sign(b) > 0)
Chris@16:    {
Chris@16:       if(eval_bit_test(b, 0))
Chris@16:       {
Chris@16:          eval_multiply(t, x, y);
Chris@16:          eval_modulus(x, t, static_cast<i1_type>(c));
Chris@16:       }
Chris@16:       eval_multiply(t, y, y);
Chris@16:       eval_modulus(y, t, static_cast<i1_type>(c));
Chris@16:       eval_right_shift(b, ui_type(1));
Chris@16:    }
Chris@16:    Backend x2(x);
Chris@16:    eval_modulus(result, x2, static_cast<i2_type>(c));
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, class Integer>
Chris@16: typename enable_if<is_unsigned<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
Chris@16: {
Chris@16:    typedef typename double_precision_type<Backend>::type double_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
Chris@16: 
Chris@16:    using default_ops::eval_bit_test;
Chris@16:    using default_ops::eval_get_sign;
Chris@16:    using default_ops::eval_multiply;
Chris@16:    using default_ops::eval_modulus;
Chris@16:    using default_ops::eval_right_shift;
Chris@16: 
Chris@16:    double_type x, y(a), t;
Chris@16:    x = ui_type(1u);
Chris@16: 
Chris@16:    while(b > 0)
Chris@16:    {
Chris@16:       if(b & 1)
Chris@16:       {
Chris@16:          eval_multiply(t, x, y);
Chris@16:          eval_modulus(x, t, c);
Chris@16:       }
Chris@16:       eval_multiply(t, y, y);
Chris@16:       eval_modulus(y, t, c);
Chris@16:       b >>= 1;
Chris@16:    }
Chris@16:    Backend x2(x);
Chris@16:    eval_modulus(result, x2, c);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, class Integer>
Chris@16: typename enable_if<is_signed<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
Chris@16: {
Chris@16:    if(b < 0)
Chris@16:    {
Chris@16:       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
Chris@16:    }
Chris@16:    eval_powm(result, a, static_cast<typename make_unsigned<Integer>::type>(b), c);
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, class Integer1, class Integer2>
Chris@16: typename enable_if<is_unsigned<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
Chris@16: {
Chris@16:    typedef typename double_precision_type<Backend>::type double_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<Integer1, double_type>::type i1_type;
Chris@16:    typedef typename boost::multiprecision::detail::canonical<Integer2, Backend>::type i2_type;
Chris@16: 
Chris@16:    using default_ops::eval_bit_test;
Chris@16:    using default_ops::eval_get_sign;
Chris@16:    using default_ops::eval_multiply;
Chris@16:    using default_ops::eval_modulus;
Chris@16:    using default_ops::eval_right_shift;
Chris@16: 
Chris@16:    double_type x, y(a), t;
Chris@16:    x = ui_type(1u);
Chris@16: 
Chris@16:    while(b > 0)
Chris@16:    {
Chris@16:       if(b & 1)
Chris@16:       {
Chris@16:          eval_multiply(t, x, y);
Chris@16:          eval_modulus(x, t, static_cast<i1_type>(c));
Chris@16:       }
Chris@16:       eval_multiply(t, y, y);
Chris@16:       eval_modulus(y, t, static_cast<i1_type>(c));
Chris@16:       b >>= 1;
Chris@16:    }
Chris@16:    Backend x2(x);
Chris@16:    eval_modulus(result, x2, static_cast<i2_type>(c));
Chris@16: }
Chris@16: 
Chris@16: template <class Backend, class Integer1, class Integer2>
Chris@16: typename enable_if<is_signed<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
Chris@16: {
Chris@16:    if(b < 0)
Chris@16:    {
Chris@16:       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
Chris@16:    }
Chris@16:    eval_powm(result, a, static_cast<typename make_unsigned<Integer1>::type>(b), c);
Chris@16: }
Chris@16: 
Chris@16: struct powm_func
Chris@16: {
Chris@16:    template <class T, class U, class V>
Chris@16:    void operator()(T& result, const T& b, const U& p, const V& m)const
Chris@16:    {
Chris@16:       eval_powm(result, b, p, m);
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: template <class T, class U, class V>
Chris@16: inline typename enable_if<
Chris@16:    mpl::and_<
Chris@16:       mpl::bool_<number_category<T>::value == number_kind_integer>, 
Chris@16:       mpl::or_<
Chris@16:          is_number<T>,
Chris@16:          is_number_expression<T>
Chris@16:       >,
Chris@16:       mpl::or_<
Chris@16:          is_number<U>,
Chris@16:          is_number_expression<U>,
Chris@16:          is_integral<U>
Chris@16:       >,
Chris@16:       mpl::or_<
Chris@16:          is_number<V>,
Chris@16:          is_number_expression<V>,
Chris@16:          is_integral<V>
Chris@16:       >
Chris@16:    >,
Chris@16:    detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type 
Chris@16:    powm(const T& b, const U& p, const V& mod)
Chris@16: {
Chris@16:    return detail::expression<detail::function, default_ops::powm_func, T, U, V>(
Chris@16:       default_ops::powm_func(), b, p, mod);
Chris@16: }
Chris@16: 
Chris@16: }} //namespaces
Chris@16: 
Chris@16: #endif
Chris@16: 
Chris@16: