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
multiply(Integer & result,const I2 & a,const I2 & b)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
add(Integer & result,const I2 & a,const I2 & b)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
subtract(Integer & result,const I2 & a,const I2 & b)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>
divide_qr(const Integer & x,const Integer & y,Integer & q,Integer & r)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>
integer_modulus(const I1 & x,I2 val)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
powm(const I1 & a,I2 b,I3 c)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
powm(const I1 & a,I2 b,I3 c)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>
lsb(const Integer & val)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>
msb(Integer val)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>
bit_test(const Integer & val,unsigned index)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>
bit_set(Integer & val,unsigned index)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>
bit_unset(Integer & val,unsigned index)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>
bit_flip(Integer & val,unsigned index)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>
sqrt(const Integer & x,Integer & r)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>
sqrt(const Integer & x)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
247