ceph/src/boost/libs/multiprecision/include/boost/multiprecision/detail/integer_ops.hpp

   1 ///////////////////////////////////////////////////////////////
   2 //  Copyright 2012 John Maddock. Distributed under the Boost
   3 //  Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
   5
   6 #ifndef BOOST_MP_INT_FUNC_HPP
   7 #define BOOST_MP_INT_FUNC_HPP
   8
   9 #include <boost/multiprecision/number.hpp>
  10
  11 namespace boost{ namespace multiprecision{
  12
  13 namespace default_ops
  14 {
  15
  16 template <class Backend>
  17 inline void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r)
  18 {
  19    eval_divide(q, x, y);
  20    eval_modulus(r, x, y);
  21 }
  22
  23 template <class Backend, class Integer>
  24 inline Integer eval_integer_modulus(const Backend& x, Integer val)
  25 {
  26    BOOST_MP_USING_ABS
  27    using default_ops::eval_modulus;
  28    using default_ops::eval_convert_to;
  29    typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type int_type;
  30    Backend t;
  31    eval_modulus(t, x, static_cast<int_type>(val));
  32    Integer result;
  33    eval_convert_to(&result, t);
  34    return abs(result);
  35 }
  36
  37 #ifdef BOOST_MSVC
  38 #pragma warning(push)
  39 #pragma warning(disable:4127)
  40 #endif
  41
  42 template <class B>
  43 inline void eval_gcd(B& result, const B& a, const B& b)
  44 {
  45    using default_ops::eval_lsb;
  46    using default_ops::eval_is_zero;
  47    using default_ops::eval_get_sign;
  48
  49    int shift;
  50
  51    B u(a), v(b);
  52
  53    int s = eval_get_sign(u);
  54
  55    /* GCD(0,x) := x */
  56    if(s < 0)
  57    {
  58       u.negate();
  59    }
  60    else if(s == 0)
  61    {
  62       result = v;
  63       return;
  64    }
  65    s = eval_get_sign(v);
  66    if(s < 0)
  67    {
  68       v.negate();
  69    }
  70    else if(s == 0)
  71    {
  72       result = u;
  73       return;
  74    }
  75
  76    /* Let shift := lg K, where K is the greatest power of 2
  77    dividing both u and v. */
  78
  79    unsigned us = eval_lsb(u);
  80    unsigned vs = eval_lsb(v);
  81    shift = (std::min)(us, vs);
  82    eval_right_shift(u, us);
  83    eval_right_shift(v, vs);
  84
  85    do
  86    {
  87       /* Now u and v are both odd, so diff(u, v) is even.
  88       Let u = min(u, v), v = diff(u, v)/2. */
  89       s = u.compare(v);
  90       if(s > 0)
  91          u.swap(v);
  92       if(s == 0)
  93          break;
  94       eval_subtract(v, u);
  95       vs = eval_lsb(v);
  96       eval_right_shift(v, vs);
  97    }
  98    while(true);
  99
 100    result = u;
 101    eval_left_shift(result, shift);
 102 }
 103
 104 #ifdef BOOST_MSVC
 105 #pragma warning(pop)
 106 #endif
 107
 108 template <class B>
 109 inline void eval_lcm(B& result, const B& a, const B& b)
 110 {
 111    typedef typename mpl::front<typename B::unsigned_types>::type ui_type;
 112    B t;
 113    eval_gcd(t, a, b);
 114
 115    if(eval_is_zero(t))
 116    {
 117       result = static_cast<ui_type>(0);
 118    }
 119    else
 120    {
 121       eval_divide(result, a, t);
 122       eval_multiply(result, b);
 123    }
 124    if(eval_get_sign(result) < 0)
 125       result.negate();
 126 }
 127
 128 }
 129
 130 template <class Backend, expression_template_option ExpressionTemplates>
 131 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
 132    divide_qr(const number<Backend, ExpressionTemplates>& x, const number<Backend, ExpressionTemplates>& y,
 133    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
 134 {
 135    using default_ops::eval_qr;
 136    eval_qr(x.backend(), y.backend(), q.backend(), r.backend());
 137 }
 138
 139 template <class Backend, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
 140 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
 141    divide_qr(const number<Backend, ExpressionTemplates>& x, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& y,
 142    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
 143 {
 144    divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r);
 145 }
 146
 147 template <class tag, class A1, class A2, class A3, class A4, class Backend, expression_template_option ExpressionTemplates>
 148 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
 149    divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const number<Backend, ExpressionTemplates>& y,
 150    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
 151 {
 152    divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r);
 153 }
 154
 155 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>
 156 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type
 157    divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& y,
 158    number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
 159 {
 160    divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r);
 161 }
 162
 163 template <class Backend, expression_template_option ExpressionTemplates, class Integer>
 164 inline typename enable_if<mpl::and_<is_integral<Integer>, mpl::bool_<number_category<Backend>::value == number_kind_integer> >, Integer>::type
 165    integer_modulus(const number<Backend, ExpressionTemplates>& x, Integer val)
 166 {
 167    using default_ops::eval_integer_modulus;
 168    return eval_integer_modulus(x.backend(), val);
 169 }
 170
 171 template <class tag, class A1, class A2, class A3, class A4, class Integer>
 172 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
 173    integer_modulus(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, Integer val)
 174 {
 175    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type result_type;
 176    return integer_modulus(result_type(x), val);
 177 }
 178
 179 template <class Backend, expression_template_option ExpressionTemplates>
 180 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type
 181    lsb(const number<Backend, ExpressionTemplates>& x)
 182 {
 183    using default_ops::eval_lsb;
 184    return eval_lsb(x.backend());
 185 }
 186
 187 template <class tag, class A1, class A2, class A3, class A4>
 188 inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type
 189    lsb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
 190 {
 191    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
 192    number_type n(x);
 193    using default_ops::eval_lsb;
 194    return eval_lsb(n.backend());
 195 }
 196
 197 template <class Backend, expression_template_option ExpressionTemplates>
 198 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type
 199    msb(const number<Backend, ExpressionTemplates>& x)
 200 {
 201    using default_ops::eval_msb;
 202    return eval_msb(x.backend());
 203 }
 204
 205 template <class tag, class A1, class A2, class A3, class A4>
 206 inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type
 207    msb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
 208 {
 209    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
 210    number_type n(x);
 211    using default_ops::eval_msb;
 212    return eval_msb(n.backend());
 213 }
 214
 215 template <class Backend, expression_template_option ExpressionTemplates>
 216 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, bool>::type
 217    bit_test(const number<Backend, ExpressionTemplates>& x, unsigned index)
 218 {
 219    using default_ops::eval_bit_test;
 220    return eval_bit_test(x.backend(), index);
 221 }
 222
 223 template <class tag, class A1, class A2, class A3, class A4>
 224 inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, bool>::type
 225    bit_test(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, unsigned index)
 226 {
 227    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
 228    number_type n(x);
 229    using default_ops::eval_bit_test;
 230    return eval_bit_test(n.backend(), index);
 231 }
 232
 233 template <class Backend, expression_template_option ExpressionTemplates>
 234 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
 235    bit_set(number<Backend, ExpressionTemplates>& x, unsigned index)
 236 {
 237    using default_ops::eval_bit_set;
 238    eval_bit_set(x.backend(), index);
 239    return x;
 240 }
 241
 242 template <class Backend, expression_template_option ExpressionTemplates>
 243 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
 244    bit_unset(number<Backend, ExpressionTemplates>& x, unsigned index)
 245 {
 246    using default_ops::eval_bit_unset;
 247    eval_bit_unset(x.backend(), index);
 248    return x;
 249 }
 250
 251 template <class Backend, expression_template_option ExpressionTemplates>
 252 inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
 253    bit_flip(number<Backend, ExpressionTemplates>& x, unsigned index)
 254 {
 255    using default_ops::eval_bit_flip;
 256    eval_bit_flip(x.backend(), index);
 257    return x;
 258 }
 259
 260 namespace default_ops{
 261
 262 //
 263 // Within powm, we need a type with twice as many digits as the argument type, define
 264 // a traits class to obtain that type:
 265 //
 266 template <class Backend>
 267 struct double_precision_type
 268 {
 269    typedef Backend type;
 270 };
 271
 272 //
 273 // If the exponent is a signed integer type, then we need to
 274 // check the value is positive:
 275 //
 276 template <class Backend>
 277 inline void check_sign_of_backend(const Backend& v, const mpl::true_)
 278 {
 279    if(eval_get_sign(v) < 0)
 280    {
 281       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
 282    }
 283 }
 284 template <class Backend>
 285 inline void check_sign_of_backend(const Backend&, const mpl::false_){}
 286 //
 287 // Calculate (a^p)%c:
 288 //
 289 template <class Backend>
 290 void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c)
 291 {
 292    using default_ops::eval_bit_test;
 293    using default_ops::eval_get_sign;
 294    using default_ops::eval_multiply;
 295    using default_ops::eval_modulus;
 296    using default_ops::eval_right_shift;
 297
 298    typedef typename double_precision_type<Backend>::type double_type;
 299    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
 300
 301    check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());
 302
 303    double_type x, y(a), b(p), t;
 304    x = ui_type(1u);
 305
 306    while(eval_get_sign(b) > 0)
 307    {
 308       if(eval_bit_test(b, 0))
 309       {
 310          eval_multiply(t, x, y);
 311          eval_modulus(x, t, c);
 312       }
 313       eval_multiply(t, y, y);
 314       eval_modulus(y, t, c);
 315       eval_right_shift(b, ui_type(1));
 316    }
 317    Backend x2(x);
 318    eval_modulus(result, x2, c);
 319 }
 320
 321 template <class Backend, class Integer>
 322 void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c)
 323 {
 324    typedef typename double_precision_type<Backend>::type double_type;
 325    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
 326    typedef typename boost::multiprecision::detail::canonical<Integer, double_type>::type i1_type;
 327    typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type i2_type;
 328
 329    using default_ops::eval_bit_test;
 330    using default_ops::eval_get_sign;
 331    using default_ops::eval_multiply;
 332    using default_ops::eval_modulus;
 333    using default_ops::eval_right_shift;
 334
 335    check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());
 336
 337    if(eval_get_sign(p) < 0)
 338    {
 339       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
 340    }
 341
 342    double_type x, y(a), b(p), t;
 343    x = ui_type(1u);
 344
 345    while(eval_get_sign(b) > 0)
 346    {
 347       if(eval_bit_test(b, 0))
 348       {
 349          eval_multiply(t, x, y);
 350          eval_modulus(x, t, static_cast<i1_type>(c));
 351       }
 352       eval_multiply(t, y, y);
 353       eval_modulus(y, t, static_cast<i1_type>(c));
 354       eval_right_shift(b, ui_type(1));
 355    }
 356    Backend x2(x);
 357    eval_modulus(result, x2, static_cast<i2_type>(c));
 358 }
 359
 360 template <class Backend, class Integer>
 361 typename enable_if<is_unsigned<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
 362 {
 363    typedef typename double_precision_type<Backend>::type double_type;
 364    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
 365
 366    using default_ops::eval_bit_test;
 367    using default_ops::eval_get_sign;
 368    using default_ops::eval_multiply;
 369    using default_ops::eval_modulus;
 370    using default_ops::eval_right_shift;
 371
 372    double_type x, y(a), t;
 373    x = ui_type(1u);
 374
 375    while(b > 0)
 376    {
 377       if(b & 1)
 378       {
 379          eval_multiply(t, x, y);
 380          eval_modulus(x, t, c);
 381       }
 382       eval_multiply(t, y, y);
 383       eval_modulus(y, t, c);
 384       b >>= 1;
 385    }
 386    Backend x2(x);
 387    eval_modulus(result, x2, c);
 388 }
 389
 390 template <class Backend, class Integer>
 391 typename enable_if<is_signed<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
 392 {
 393    if(b < 0)
 394    {
 395       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
 396    }
 397    eval_powm(result, a, static_cast<typename make_unsigned<Integer>::type>(b), c);
 398 }
 399
 400 template <class Backend, class Integer1, class Integer2>
 401 typename enable_if<is_unsigned<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
 402 {
 403    typedef typename double_precision_type<Backend>::type double_type;
 404    typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
 405    typedef typename boost::multiprecision::detail::canonical<Integer1, double_type>::type i1_type;
 406    typedef typename boost::multiprecision::detail::canonical<Integer2, Backend>::type i2_type;
 407
 408    using default_ops::eval_bit_test;
 409    using default_ops::eval_get_sign;
 410    using default_ops::eval_multiply;
 411    using default_ops::eval_modulus;
 412    using default_ops::eval_right_shift;
 413
 414    double_type x, y(a), t;
 415    x = ui_type(1u);
 416
 417    while(b > 0)
 418    {
 419       if(b & 1)
 420       {
 421          eval_multiply(t, x, y);
 422          eval_modulus(x, t, static_cast<i1_type>(c));
 423       }
 424       eval_multiply(t, y, y);
 425       eval_modulus(y, t, static_cast<i1_type>(c));
 426       b >>= 1;
 427    }
 428    Backend x2(x);
 429    eval_modulus(result, x2, static_cast<i2_type>(c));
 430 }
 431
 432 template <class Backend, class Integer1, class Integer2>
 433 typename enable_if<is_signed<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
 434 {
 435    if(b < 0)
 436    {
 437       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
 438    }
 439    eval_powm(result, a, static_cast<typename make_unsigned<Integer1>::type>(b), c);
 440 }
 441
 442 struct powm_func
 443 {
 444    template <class T, class U, class V>
 445    void operator()(T& result, const T& b, const U& p, const V& m)const
 446    {
 447       eval_powm(result, b, p, m);
 448    }
 449 };
 450
 451 }
 452
 453 template <class T, class U, class V>
 454 inline typename enable_if<
 455    mpl::and_<
 456       mpl::bool_<number_category<T>::value == number_kind_integer>,
 457       mpl::or_<
 458          is_number<T>,
 459          is_number_expression<T>
 460       >,
 461       mpl::or_<
 462          is_number<U>,
 463          is_number_expression<U>,
 464          is_integral<U>
 465       >,
 466       mpl::or_<
 467          is_number<V>,
 468          is_number_expression<V>,
 469          is_integral<V>
 470       >
 471    >,
 472    typename mpl::if_<
 473       is_no_et_number<T>,
 474       T,
 475       typename mpl::if_<
 476          is_no_et_number<U>,
 477          U,
 478          typename mpl::if_<
 479             is_no_et_number<V>,
 480             V,
 481             detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type
 482          >::type
 483       >::type
 484    >::type
 485    powm(const T& b, const U& p, const V& mod)
 486 {
 487    return detail::expression<detail::function, default_ops::powm_func, T, U, V>(
 488       default_ops::powm_func(), b, p, mod);
 489 }
 490
 491 }} //namespaces
 492
 493 #endif
 494
 495