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
12namespace boost{
13namespace multiprecision{
14
15template <class Integer, class I2>
16typename 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}
21template <class Integer, class I2>
22typename 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}
27template <class Integer, class I2>
28typename 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
34template <class Integer>
35typename 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
41template <class I1, class I2>
42typename 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
47namespace 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//
54template <class I>
55struct 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
81template <class I1, class I2, class I3>
82typename 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
104template <class I1, class I2, class I3>
105inline 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
115template <class Integer>
116typename 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
132template <class Integer>
133typename 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
149template <class Integer>
150typename 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
160template <class Integer>
161typename 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
172template <class Integer>
173typename 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
184template <class Integer>
185typename 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
196template <class Integer>
197typename 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
242template <class Integer>
243typename 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

source code of boost/boost/multiprecision/integer.hpp