boost/multiprecision/integer.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 https://www.boost.org/LICENSE_1_0.txt
   5
   6 #ifndef BOOST_MP_INTEGER_HPP
   7 #define BOOST_MP_INTEGER_HPP
   8
   9 #include <boost/multiprecision/cpp_int.hpp>
  10 #include <boost/multiprecision/detail/bitscan.hpp>
  11
  12 namespace boost {
  13 namespace multiprecision {
  14
  15 template <class Integer, class I2>
  16 inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value && is_integral<I2>::value, Integer&>::type
  17 multiply(Integer& result, const I2& a, const I2& b)
  18 {
  19    return result = static_cast<Integer>(a) * static_cast<Integer>(b);
  20 }
  21 template <class Integer, class I2>
  22 inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value && is_integral<I2>::value, Integer&>::type
  23 add(Integer& result, const I2& a, const I2& b)
  24 {
  25    return result = static_cast<Integer>(a) + static_cast<Integer>(b);
  26 }
  27 template <class Integer, class I2>
  28 inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value && is_integral<I2>::value, Integer&>::type
  29 subtract(Integer& result, const I2& a, const I2& b)
  30 {
  31    return result = static_cast<Integer>(a) - static_cast<Integer>(b);
  32 }
  33
  34 template <class Integer>
  35 inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value>::type divide_qr(const Integer& x, const Integer& y, Integer& q, Integer& r)
  36 {
  37    q = x / y;
  38    r = x % y;
  39 }
  40
  41 template <class I1, class I2>
  42 inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<I1>::value && is_integral<I2>::value, I2>::type integer_modulus(const I1& x, I2 val)
  43 {
  44    return static_cast<I2>(x % val);
  45 }
  46
  47 namespace detail {
  48 //
  49 // Figure out the kind of integer that has twice as many bits as some builtin
  50 // integer type I.  Use a native type if we can (including types which may not
  51 // be recognised by boost::int_t because they're larger than boost::long_long_type),
  52 // otherwise synthesize a cpp_int to do the job.
  53 //
  54 template <class I>
  55 struct double_integer
  56 {
  57    static const unsigned int_t_digits =
  58        2 * sizeof(I) <= sizeof(boost::long_long_type) ? std::numeric_limits<I>::digits * 2 : 1;
  59
  60    typedef typename mpl::if_c<
  61        2 * sizeof(I) <= sizeof(boost::long_long_type),
  62        typename mpl::if_c<
  63            is_signed<I>::value,
  64            typename boost::int_t<int_t_digits>::least,
  65            typename boost::uint_t<int_t_digits>::least>::type,
  66        typename mpl::if_c<
  67            2 * sizeof(I) <= sizeof(double_limb_type),
  68            typename mpl::if_c<
  69                is_signed<I>::value,
  70                signed_double_limb_type,
  71                double_limb_type>::type,
  72            number<cpp_int_backend<sizeof(I) * CHAR_BIT * 2, sizeof(I) * CHAR_BIT * 2, (is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> > >::type>::type type;
  73 };
  74
  75 } // namespace detail
  76
  77 template <class I1, class I2, class I3>
  78 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<I1>::value && is_unsigned<I2>::value && is_integral<I3>::value, I1>::type
  79 powm(const I1& a, I2 b, I3 c)
  80 {
  81    typedef typename detail::double_integer<I1>::type double_type;
  82
  83    I1          x(1), y(a);
  84    double_type result(0);
  85
  86    while (b > 0)
  87    {
  88       if (b & 1)
  89       {
  90          multiply(result, x, y);
  91          x = integer_modulus(result, c);
  92       }
  93       multiply(result, y, y);
  94       y = integer_modulus(result, c);
  95       b >>= 1;
  96    }
  97    return x % c;
  98 }
  99
 100 template <class I1, class I2, class I3>
 101 inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<I1>::value && is_signed<I2>::value && is_integral<I3>::value, I1>::type
 102 powm(const I1& a, I2 b, I3 c)
 103 {
 104    if (b < 0)
 105    {
 106       BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
 107    }
 108    return powm(a, static_cast<typename make_unsigned<I2>::type>(b), c);
 109 }
 110
 111 template <class Integer>
 112 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, unsigned>::type lsb(const Integer& val)
 113 {
 114    if (val <= 0)
 115    {
 116       if (val == 0)
 117       {
 118          BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
 119       }
 120       else
 121       {
 122          BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
 123       }
 124    }
 125    return detail::find_lsb(val);
 126 }
 127
 128 template <class Integer>
 129 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, unsigned>::type msb(Integer val)
 130 {
 131    if (val <= 0)
 132    {
 133       if (val == 0)
 134       {
 135          BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
 136       }
 137       else
 138       {
 139          BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
 140       }
 141    }
 142    return detail::find_msb(val);
 143 }
 144
 145 template <class Integer>
 146 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, bool>::type bit_test(const Integer& val, unsigned index)
 147 {
 148    Integer mask = 1;
 149    if (index >= sizeof(Integer) * CHAR_BIT)
 150       return 0;
 151    if (index)
 152       mask <<= index;
 153    return val & mask ? true : false;
 154 }
 155
 156 template <class Integer>
 157 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, unsigned index)
 158 {
 159    Integer mask = 1;
 160    if (index >= sizeof(Integer) * CHAR_BIT)
 161       return val;
 162    if (index)
 163       mask <<= index;
 164    val |= mask;
 165    return val;
 166 }
 167
 168 template <class Integer>
 169 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, unsigned index)
 170 {
 171    Integer mask = 1;
 172    if (index >= sizeof(Integer) * CHAR_BIT)
 173       return val;
 174    if (index)
 175       mask <<= index;
 176    val &= ~mask;
 177    return val;
 178 }
 179
 180 template <class Integer>
 181 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, unsigned index)
 182 {
 183    Integer mask = 1;
 184    if (index >= sizeof(Integer) * CHAR_BIT)
 185       return val;
 186    if (index)
 187       mask <<= index;
 188    val ^= mask;
 189    return val;
 190 }
 191
 192 template <class Integer>
 193 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r)
 194 {
 195    //
 196    // This is slow bit-by-bit integer square root, see for example
 197    // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
 198    // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
 199    // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
 200    // at some point.
 201    //
 202    Integer s = 0;
 203    if (x == 0)
 204    {
 205       r = 0;
 206       return s;
 207    }
 208    int g = msb(x);
 209    if (g == 0)
 210    {
 211       r = 1;
 212       return s;
 213    }
 214
 215    Integer t = 0;
 216    r         = x;
 217    g /= 2;
 218    bit_set(s, g);
 219    bit_set(t, 2 * g);
 220    r = x - t;
 221    --g;
 222    do
 223    {
 224       t = s;
 225       t <<= g + 1;
 226       bit_set(t, 2 * g);
 227       if (t <= r)
 228       {
 229          bit_set(s, g);
 230          r -= t;
 231       }
 232       --g;
 233    } while (g >= 0);
 234    return s;
 235 }
 236
 237 template <class Integer>
 238 BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x)
 239 {
 240    Integer r(0);
 241    return sqrt(x, r);
 242 }
 243
 244 }} // namespace boost::multiprecision
 245
 246 #endif