multiprecision/detail/integer_ops.hpp

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2012 John Maddock. Distributed under the Boost
0003 //  Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
0005
0006 #ifndef BOOST_MP_DETAIL_INTEGER_OPS_HPP
0007 #define BOOST_MP_DETAIL_INTEGER_OPS_HPP
0008
0009 #include <boost/multiprecision/number.hpp>
0010 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0011
0012 namespace boost { namespace multiprecision {
0013
0014 namespace default_ops {
0015
0016 template <class Backend>
0017 inline BOOST_MP_CXX14_CONSTEXPR void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r)
0018 {
0019    eval_divide(q, x, y);
0020    eval_modulus(r, x, y);
0021 }
0022
0023 template <class Backend, class Integer>
0024 inline BOOST_MP_CXX14_CONSTEXPR Integer eval_integer_modulus(const Backend& x, Integer val)
0025 {
0026    BOOST_MP_USING_ABS
0027    using default_ops::eval_convert_to;
0028    using default_ops::eval_modulus;
0029    using int_type = typename boost::multiprecision::detail::canonical<Integer, Backend>::type;
0030    Backend                                                                           t;
0031    eval_modulus(t, x, static_cast<int_type>(val));
0032    Integer result(0);
0033    eval_convert_to(&result, t);
0034    return abs(result);
0035 }
0036
0037 template <class B>
0038 inline BOOST_MP_CXX14_CONSTEXPR void eval_gcd(B& result, const B& a, const B& b)
0039 {
0040    using default_ops::eval_get_sign;
0041    using default_ops::eval_is_zero;
0042    using default_ops::eval_lsb;
0043
0044    std::ptrdiff_t shift(0);
0045
0046    B u(a), v(b);
0047
0048    int s = eval_get_sign(u);
0049
0050    /* GCD(0,x) := x */
0051    if (s < 0)
0052    {
0053       u.negate();
0054    }
0055    else if (s == 0)
0056    {
0057       result = v;
0058       return;
0059    }
0060    s = eval_get_sign(v);
0061    if (s < 0)
0062    {
0063       v.negate();
0064    }
0065    else if (s == 0)
0066    {
0067       result = u;
0068       return;
0069    }
0070
0071    /* Let shift := lg K, where K is the greatest power of 2
0072    dividing both u and v. */
0073
0074    std::size_t us = eval_lsb(u);
0075    std::size_t vs = eval_lsb(v);
0076    shift       = static_cast<std::ptrdiff_t>((std::min)(us, vs));
0077    eval_right_shift(u, us);
0078    eval_right_shift(v, vs);
0079
0080    do
0081    {
0082       /* Now u and v are both odd, so diff(u, v) is even.
0083       Let u = min(u, v), v = diff(u, v)/2. */
0084       s = u.compare(v);
0085       if (s > 0)
0086          u.swap(v);
0087       if (s == 0)
0088          break;
0089       eval_subtract(v, u);
0090       vs = eval_lsb(v);
0091       eval_right_shift(v, vs);
0092    } while (true);
0093
0094    result = u;
0095    eval_left_shift(result, shift);
0096 }
0097
0098 template <class B>
0099 inline BOOST_MP_CXX14_CONSTEXPR void eval_lcm(B& result, const B& a, const B& b)
0100 {
0101    using ui_type = typename std::tuple_element<0, typename B::unsigned_types>::type;
0102    B                                                             t;
0103    eval_gcd(t, a, b);
0104
0105    if (eval_is_zero(t))
0106    {
0107       result = static_cast<ui_type>(0);
0108    }
0109    else
0110    {
0111       eval_divide(result, a, t);
0112       eval_multiply(result, b);
0113    }
0114    if (eval_get_sign(result) < 0)
0115       result.negate();
0116 }
0117
0118 } // namespace default_ops
0119
0120 template <class Backend, expression_template_option ExpressionTemplates>
0121 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
0122 divide_qr(const number<Backend, ExpressionTemplates>& x, const number<Backend, ExpressionTemplates>& y,
0123           number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
0124 {
0125    using default_ops::eval_qr;
0126    eval_qr(x.backend(), y.backend(), q.backend(), r.backend());
0127 }
0128
0129 template <class Backend, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
0130 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
0131 divide_qr(const number<Backend, ExpressionTemplates>& x, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& y,
0132           number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
0133 {
0134    divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r);
0135 }
0136
0137 template <class tag, class A1, class A2, class A3, class A4, class Backend, expression_template_option ExpressionTemplates>
0138 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
0139 divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const number<Backend, ExpressionTemplates>& y,
0140           number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
0141 {
0142    divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r);
0143 }
0144
0145 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>
0146 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
0147 divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& y,
0148           number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
0149 {
0150    divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r);
0151 }
0152
0153 template <class Backend, expression_template_option ExpressionTemplates, class Integer>
0154 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && (number_category<Backend>::value == number_kind_integer), Integer>::type
0155 integer_modulus(const number<Backend, ExpressionTemplates>& x, Integer val)
0156 {
0157    using default_ops::eval_integer_modulus;
0158    return eval_integer_modulus(x.backend(), val);
0159 }
0160
0161 template <class tag, class A1, class A2, class A3, class A4, class Integer>
0162 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && (number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer), Integer>::type
0163 integer_modulus(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, Integer val)
0164 {
0165    using result_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
0166    return integer_modulus(result_type(x), val);
0167 }
0168
0169 template <class Backend, expression_template_option ExpressionTemplates>
0170 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, std::size_t>::type
0171 lsb(const number<Backend, ExpressionTemplates>& x)
0172 {
0173    using default_ops::eval_lsb;
0174    return eval_lsb(x.backend());
0175 }
0176
0177 template <class tag, class A1, class A2, class A3, class A4>
0178 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, std::size_t>::type
0179 lsb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
0180 {
0181    using number_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
0182    number_type                                                                           n(x);
0183    using default_ops::eval_lsb;
0184    return eval_lsb(n.backend());
0185 }
0186
0187 template <class Backend, expression_template_option ExpressionTemplates>
0188 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, std::size_t>::type
0189 msb(const number<Backend, ExpressionTemplates>& x)
0190 {
0191    using default_ops::eval_msb;
0192    return eval_msb(x.backend());
0193 }
0194
0195 template <class tag, class A1, class A2, class A3, class A4>
0196 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, std::size_t>::type
0197 msb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
0198 {
0199    using number_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
0200    number_type                                                                           n(x);
0201    using default_ops::eval_msb;
0202    return eval_msb(n.backend());
0203 }
0204
0205 template <class Backend, expression_template_option ExpressionTemplates>
0206 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, bool>::type
0207 bit_test(const number<Backend, ExpressionTemplates>& x, std::size_t index)
0208 {
0209    using default_ops::eval_bit_test;
0210    return eval_bit_test(x.backend(), index);
0211 }
0212
0213 template <class tag, class A1, class A2, class A3, class A4>
0214 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, bool>::type
0215 bit_test(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, std::size_t index)
0216 {
0217    using number_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
0218    number_type                                                                           n(x);
0219    using default_ops::eval_bit_test;
0220    return eval_bit_test(n.backend(), index);
0221 }
0222
0223 template <class Backend, expression_template_option ExpressionTemplates>
0224 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
0225 bit_set(number<Backend, ExpressionTemplates>& x, std::size_t index)
0226 {
0227    using default_ops::eval_bit_set;
0228    eval_bit_set(x.backend(), index);
0229    return x;
0230 }
0231
0232 template <class Backend, expression_template_option ExpressionTemplates>
0233 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
0234 bit_unset(number<Backend, ExpressionTemplates>& x, std::size_t index)
0235 {
0236    using default_ops::eval_bit_unset;
0237    eval_bit_unset(x.backend(), index);
0238    return x;
0239 }
0240
0241 template <class Backend, expression_template_option ExpressionTemplates>
0242 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
0243 bit_flip(number<Backend, ExpressionTemplates>& x, std::size_t index)
0244 {
0245    using default_ops::eval_bit_flip;
0246    eval_bit_flip(x.backend(), index);
0247    return x;
0248 }
0249
0250 namespace default_ops {
0251
0252 //
0253 // Within powm, we need a type with twice as many digits as the argument type, define
0254 // a traits class to obtain that type:
0255 //
0256 template <class Backend>
0257 struct double_precision_type
0258 {
0259    using type = Backend;
0260 };
0261
0262 //
0263 // If the exponent is a signed integer type, then we need to
0264 // check the value is positive:
0265 //
0266 template <class Backend>
0267 inline BOOST_MP_CXX14_CONSTEXPR void check_sign_of_backend(const Backend& v, const std::integral_constant<bool, true>)
0268 {
0269    if (eval_get_sign(v) < 0)
0270    {
0271       BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
0272    }
0273 }
0274 template <class Backend>
0275 inline BOOST_MP_CXX14_CONSTEXPR void check_sign_of_backend(const Backend&, const std::integral_constant<bool, false>) {}
0276 //
0277 // Calculate (a^p)%c:
0278 //
0279 template <class Backend>
0280 BOOST_MP_CXX14_CONSTEXPR void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c)
0281 {
0282    using default_ops::eval_bit_test;
0283    using default_ops::eval_get_sign;
0284    using default_ops::eval_modulus;
0285    using default_ops::eval_multiply;
0286    using default_ops::eval_right_shift;
0287
0288    using double_type = typename double_precision_type<Backend>::type                                      ;
0289    using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
0290
0291    check_sign_of_backend(p, std::integral_constant<bool, std::numeric_limits<number<Backend> >::is_signed>());
0292
0293    double_type x, y(a), b(p), t;
0294    x = ui_type(1u);
0295
0296    while (eval_get_sign(b) > 0)
0297    {
0298       if (eval_bit_test(b, 0))
0299       {
0300          eval_multiply(t, x, y);
0301          eval_modulus(x, t, c);
0302       }
0303       eval_multiply(t, y, y);
0304       eval_modulus(y, t, c);
0305       eval_right_shift(b, ui_type(1));
0306    }
0307    Backend x2(x);
0308    eval_modulus(result, x2, c);
0309 }
0310
0311 template <class Backend, class Integer>
0312 BOOST_MP_CXX14_CONSTEXPR void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c)
0313 {
0314    using double_type = typename double_precision_type<Backend>::type                                      ;
0315    using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
0316    using i1_type = typename boost::multiprecision::detail::canonical<Integer, double_type>::type      ;
0317    using i2_type = typename boost::multiprecision::detail::canonical<Integer, Backend>::type          ;
0318
0319    using default_ops::eval_bit_test;
0320    using default_ops::eval_get_sign;
0321    using default_ops::eval_modulus;
0322    using default_ops::eval_multiply;
0323    using default_ops::eval_right_shift;
0324
0325    check_sign_of_backend(p, std::integral_constant<bool, std::numeric_limits<number<Backend> >::is_signed>());
0326
0327    if (eval_get_sign(p) < 0)
0328    {
0329       BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
0330    }
0331
0332    double_type x, y(a), b(p), t;
0333    x = ui_type(1u);
0334
0335    while (eval_get_sign(b) > 0)
0336    {
0337       if (eval_bit_test(b, 0))
0338       {
0339          eval_multiply(t, x, y);
0340          eval_modulus(x, t, static_cast<i1_type>(c));
0341       }
0342       eval_multiply(t, y, y);
0343       eval_modulus(y, t, static_cast<i1_type>(c));
0344       eval_right_shift(b, ui_type(1));
0345    }
0346    Backend x2(x);
0347    eval_modulus(result, x2, static_cast<i2_type>(c));
0348 }
0349
0350 template <class Backend, class Integer>
0351 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
0352 {
0353    using double_type = typename double_precision_type<Backend>::type                                      ;
0354    using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
0355
0356    using default_ops::eval_bit_test;
0357    using default_ops::eval_get_sign;
0358    using default_ops::eval_modulus;
0359    using default_ops::eval_multiply;
0360    using default_ops::eval_right_shift;
0361
0362    double_type x, y(a), t;
0363    x = ui_type(1u);
0364
0365    while (b > 0)
0366    {
0367       if (b & 1)
0368       {
0369          eval_multiply(t, x, y);
0370          eval_modulus(x, t, c);
0371       }
0372       eval_multiply(t, y, y);
0373       eval_modulus(y, t, c);
0374       b >>= 1;
0375    }
0376    Backend x2(x);
0377    eval_modulus(result, x2, c);
0378 }
0379
0380 template <class Backend, class Integer>
0381 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value>::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
0382 {
0383    if (b < 0)
0384    {
0385       BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
0386    }
0387    eval_powm(result, a, static_cast<typename boost::multiprecision::detail::make_unsigned<Integer>::type>(b), c);
0388 }
0389
0390 template <class Backend, class Integer1, class Integer2>
0391 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer1>::value >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
0392 {
0393    using double_type = typename double_precision_type<Backend>::type                                      ;
0394    using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
0395    using i1_type = typename boost::multiprecision::detail::canonical<Integer1, double_type>::type     ;
0396    using i2_type = typename boost::multiprecision::detail::canonical<Integer2, Backend>::type         ;
0397
0398    using default_ops::eval_bit_test;
0399    using default_ops::eval_get_sign;
0400    using default_ops::eval_modulus;
0401    using default_ops::eval_multiply;
0402    using default_ops::eval_right_shift;
0403
0404    double_type x, y(a), t;
0405    x = ui_type(1u);
0406
0407    while (b > 0)
0408    {
0409       if (b & 1)
0410       {
0411          eval_multiply(t, x, y);
0412          eval_modulus(x, t, static_cast<i1_type>(c));
0413       }
0414       eval_multiply(t, y, y);
0415       eval_modulus(y, t, static_cast<i1_type>(c));
0416       b >>= 1;
0417    }
0418    Backend x2(x);
0419    eval_modulus(result, x2, static_cast<i2_type>(c));
0420 }
0421
0422 template <class Backend, class Integer1, class Integer2>
0423 BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer1>::value && boost::multiprecision::detail::is_integral<Integer1>::value>::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
0424 {
0425    if (b < 0)
0426    {
0427       BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
0428    }
0429    eval_powm(result, a, static_cast<typename boost::multiprecision::detail::make_unsigned<Integer1>::type>(b), c);
0430 }
0431
0432 struct powm_func
0433 {
0434    template <class T, class U, class V>
0435    BOOST_MP_CXX14_CONSTEXPR void operator()(T& result, const T& b, const U& p, const V& m) const
0436    {
0437       eval_powm(result, b, p, m);
0438    }
0439    template <class R, class T, class U, class V>
0440    BOOST_MP_CXX14_CONSTEXPR void operator()(R& result, const T& b, const U& p, const V& m) const
0441    {
0442       T temp;
0443       eval_powm(temp, b, p, m);
0444       result = std::move(temp);
0445    }
0446 };
0447
0448 } // namespace default_ops
0449
0450 template <class T, class U, class V>
0451 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
0452         (number_category<T>::value == number_kind_integer) &&
0453         (is_number<T>::value || is_number_expression<T>::value) &&
0454         (is_number<U>::value || is_number_expression<U>::value || boost::multiprecision::detail::is_integral<U>::value) &&
0455         (is_number<V>::value || is_number_expression<V>::value || boost::multiprecision::detail::is_integral<V>::value),
0456     typename std::conditional<
0457         is_no_et_number<T>::value,
0458         T,
0459         typename std::conditional<
0460             is_no_et_number<U>::value,
0461             U,
0462             typename std::conditional<
0463                 is_no_et_number<V>::value,
0464                 V,
0465                 detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type>::type>::type>::type
0466 powm(const T& b, const U& p, const V& mod)
0467 {
0468    return detail::expression<detail::function, default_ops::powm_func, T, U, V>(
0469        default_ops::powm_func(), b, p, mod);
0470 }
0471
0472 }} // namespace boost::multiprecision
0473
0474 #endif