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_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 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 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 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 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 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
66 >::type,
67 typename mpl::if_c<
68 2 * sizeof(I) <= sizeof(double_limb_type),
69 typename mpl::if_c<
70 is_signed<I>::value,
71 signed_double_limb_type,
72 double_limb_type
73 >::type,
74 number<cpp_int_backend<sizeof(I)*CHAR_BIT*2, sizeof(I)*CHAR_BIT*2, (is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> >
75 >::type
76 >::type type;
77 };
78
79 }
80
81 template <class I1, class I2, class I3>
82 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)83 powm(const I1& a, I2 b, I3 c)
84 {
85 typedef typename detail::double_integer<I1>::type double_type;
86
87 I1 x(1), y(a);
88 double_type result;
89
90 while(b > 0)
91 {
92 if(b & 1)
93 {
94 multiply(result, x, y);
95 x = integer_modulus(result, c);
96 }
97 multiply(result, y, y);
98 y = integer_modulus(result, c);
99 b >>= 1;
100 }
101 return x % c;
102 }
103
104 template <class I1, class I2, class I3>
105 inline 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)106 powm(const I1& a, I2 b, I3 c)
107 {
108 if(b < 0)
109 {
110 BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
111 }
112 return powm(a, static_cast<typename make_unsigned<I2>::type>(b), c);
113 }
114
115 template <class Integer>
lsb(const Integer & val)116 typename enable_if_c<is_integral<Integer>::value, unsigned>::type lsb(const Integer& val)
117 {
118 if(val <= 0)
119 {
120 if(val == 0)
121 {
122 BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
123 }
124 else
125 {
126 BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
127 }
128 }
129 return detail::find_lsb(val);
130 }
131
132 template <class Integer>
msb(Integer val)133 typename enable_if_c<is_integral<Integer>::value, unsigned>::type msb(Integer val)
134 {
135 if(val <= 0)
136 {
137 if(val == 0)
138 {
139 BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
140 }
141 else
142 {
143 BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
144 }
145 }
146 return detail::find_msb(val);
147 }
148
149 template <class Integer>
bit_test(const Integer & val,unsigned index)150 typename enable_if_c<is_integral<Integer>::value, bool>::type bit_test(const Integer& val, unsigned index)
151 {
152 Integer mask = 1;
153 if(index >= sizeof(Integer) * CHAR_BIT)
154 return 0;
155 if(index)
156 mask <<= index;
157 return val & mask ? true : false;
158 }
159
160 template <class Integer>
bit_set(Integer & val,unsigned index)161 typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, unsigned index)
162 {
163 Integer mask = 1;
164 if(index >= sizeof(Integer) * CHAR_BIT)
165 return val;
166 if(index)
167 mask <<= index;
168 val |= mask;
169 return val;
170 }
171
172 template <class Integer>
bit_unset(Integer & val,unsigned index)173 typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, unsigned index)
174 {
175 Integer mask = 1;
176 if(index >= sizeof(Integer) * CHAR_BIT)
177 return val;
178 if(index)
179 mask <<= index;
180 val &= ~mask;
181 return val;
182 }
183
184 template <class Integer>
bit_flip(Integer & val,unsigned index)185 typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, unsigned index)
186 {
187 Integer mask = 1;
188 if(index >= sizeof(Integer) * CHAR_BIT)
189 return val;
190 if(index)
191 mask <<= index;
192 val ^= mask;
193 return val;
194 }
195
196 template <class Integer>
sqrt(const Integer & x,Integer & r)197 typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r)
198 {
199 //
200 // This is slow bit-by-bit integer square root, see for example
201 // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
202 // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
203 // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
204 // at some point.
205 //
206 Integer s = 0;
207 if(x == 0)
208 {
209 r = 0;
210 return s;
211 }
212 int g = msb(x);
213 if(g == 0)
214 {
215 r = 1;
216 return s;
217 }
218
219 Integer t = 0;
220 r = x;
221 g /= 2;
222 bit_set(s, g);
223 bit_set(t, 2 * g);
224 r = x - t;
225 --g;
226 do
227 {
228 t = s;
229 t <<= g + 1;
230 bit_set(t, 2 * g);
231 if(t <= r)
232 {
233 bit_set(s, g);
234 r -= t;
235 }
236 --g;
237 }
238 while(g >= 0);
239 return s;
240 }
241
242 template <class Integer>
sqrt(const Integer & x)243 typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x)
244 {
245 Integer r;
246 return sqrt(x, r);
247 }
248
249 }} // namespaces
250
251 #endif
252