]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | |
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 | |
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 | |
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 | 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 | 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 | |
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 | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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 |