]>
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_INT_FUNC_HPP | |
7 | #define BOOST_MP_INT_FUNC_HPP | |
8 | ||
9 | #include <boost/multiprecision/number.hpp> | |
10 | ||
11 | namespace boost{ namespace multiprecision{ | |
12 | ||
13 | namespace default_ops | |
14 | { | |
15 | ||
16 | template <class Backend> | |
17 | inline void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r) | |
18 | { | |
19 | eval_divide(q, x, y); | |
20 | eval_modulus(r, x, y); | |
21 | } | |
22 | ||
23 | template <class Backend, class Integer> | |
24 | inline Integer eval_integer_modulus(const Backend& x, Integer val) | |
25 | { | |
26 | BOOST_MP_USING_ABS | |
27 | using default_ops::eval_modulus; | |
28 | using default_ops::eval_convert_to; | |
29 | typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type int_type; | |
30 | Backend t; | |
31 | eval_modulus(t, x, static_cast<int_type>(val)); | |
32 | Integer result; | |
33 | eval_convert_to(&result, t); | |
34 | return abs(result); | |
35 | } | |
36 | ||
37 | #ifdef BOOST_MSVC | |
38 | #pragma warning(push) | |
39 | #pragma warning(disable:4127) | |
40 | #endif | |
41 | ||
42 | template <class B> | |
43 | inline void eval_gcd(B& result, const B& a, const B& b) | |
44 | { | |
45 | using default_ops::eval_lsb; | |
46 | using default_ops::eval_is_zero; | |
47 | using default_ops::eval_get_sign; | |
48 | ||
49 | int shift; | |
50 | ||
51 | B u(a), v(b); | |
52 | ||
53 | int s = eval_get_sign(u); | |
54 | ||
55 | /* GCD(0,x) := x */ | |
56 | if(s < 0) | |
57 | { | |
58 | u.negate(); | |
59 | } | |
60 | else if(s == 0) | |
61 | { | |
62 | result = v; | |
63 | return; | |
64 | } | |
65 | s = eval_get_sign(v); | |
66 | if(s < 0) | |
67 | { | |
68 | v.negate(); | |
69 | } | |
70 | else if(s == 0) | |
71 | { | |
72 | result = u; | |
73 | return; | |
74 | } | |
75 | ||
76 | /* Let shift := lg K, where K is the greatest power of 2 | |
77 | dividing both u and v. */ | |
78 | ||
79 | unsigned us = eval_lsb(u); | |
80 | unsigned vs = eval_lsb(v); | |
81 | shift = (std::min)(us, vs); | |
82 | eval_right_shift(u, us); | |
83 | eval_right_shift(v, vs); | |
84 | ||
85 | do | |
86 | { | |
87 | /* Now u and v are both odd, so diff(u, v) is even. | |
88 | Let u = min(u, v), v = diff(u, v)/2. */ | |
89 | s = u.compare(v); | |
90 | if(s > 0) | |
91 | u.swap(v); | |
92 | if(s == 0) | |
93 | break; | |
94 | eval_subtract(v, u); | |
95 | vs = eval_lsb(v); | |
96 | eval_right_shift(v, vs); | |
97 | } | |
98 | while(true); | |
99 | ||
100 | result = u; | |
101 | eval_left_shift(result, shift); | |
102 | } | |
103 | ||
104 | #ifdef BOOST_MSVC | |
105 | #pragma warning(pop) | |
106 | #endif | |
107 | ||
108 | template <class B> | |
109 | inline void eval_lcm(B& result, const B& a, const B& b) | |
110 | { | |
111 | typedef typename mpl::front<typename B::unsigned_types>::type ui_type; | |
112 | B t; | |
113 | eval_gcd(t, a, b); | |
114 | ||
115 | if(eval_is_zero(t)) | |
116 | { | |
117 | result = static_cast<ui_type>(0); | |
118 | } | |
119 | else | |
120 | { | |
121 | eval_divide(result, a, t); | |
122 | eval_multiply(result, b); | |
123 | } | |
124 | if(eval_get_sign(result) < 0) | |
125 | result.negate(); | |
126 | } | |
127 | ||
128 | } | |
129 | ||
130 | template <class Backend, expression_template_option ExpressionTemplates> | |
131 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type | |
132 | divide_qr(const number<Backend, ExpressionTemplates>& x, const number<Backend, ExpressionTemplates>& y, | |
133 | number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r) | |
134 | { | |
135 | using default_ops::eval_qr; | |
136 | eval_qr(x.backend(), y.backend(), q.backend(), r.backend()); | |
137 | } | |
138 | ||
139 | template <class Backend, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4> | |
140 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type | |
141 | divide_qr(const number<Backend, ExpressionTemplates>& x, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& y, | |
142 | number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r) | |
143 | { | |
144 | divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r); | |
145 | } | |
146 | ||
147 | template <class tag, class A1, class A2, class A3, class A4, class Backend, expression_template_option ExpressionTemplates> | |
148 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type | |
149 | divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const number<Backend, ExpressionTemplates>& y, | |
150 | number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r) | |
151 | { | |
152 | divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r); | |
153 | } | |
154 | ||
155 | template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b, class Backend, expression_template_option ExpressionTemplates> | |
156 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer>::type | |
157 | divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& y, | |
158 | number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r) | |
159 | { | |
160 | divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r); | |
161 | } | |
162 | ||
163 | template <class Backend, expression_template_option ExpressionTemplates, class Integer> | |
164 | inline typename enable_if<mpl::and_<is_integral<Integer>, mpl::bool_<number_category<Backend>::value == number_kind_integer> >, Integer>::type | |
165 | integer_modulus(const number<Backend, ExpressionTemplates>& x, Integer val) | |
166 | { | |
167 | using default_ops::eval_integer_modulus; | |
168 | return eval_integer_modulus(x.backend(), val); | |
169 | } | |
170 | ||
171 | template <class tag, class A1, class A2, class A3, class A4, class Integer> | |
172 | inline typename enable_if<mpl::and_<is_integral<Integer>, mpl::bool_<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer> >, Integer>::type | |
173 | integer_modulus(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, Integer val) | |
174 | { | |
175 | typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type result_type; | |
176 | return integer_modulus(result_type(x), val); | |
177 | } | |
178 | ||
179 | template <class Backend, expression_template_option ExpressionTemplates> | |
180 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type | |
181 | lsb(const number<Backend, ExpressionTemplates>& x) | |
182 | { | |
183 | using default_ops::eval_lsb; | |
184 | return eval_lsb(x.backend()); | |
185 | } | |
186 | ||
187 | template <class tag, class A1, class A2, class A3, class A4> | |
188 | inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type | |
189 | lsb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x) | |
190 | { | |
191 | typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type; | |
192 | number_type n(x); | |
193 | using default_ops::eval_lsb; | |
194 | return eval_lsb(n.backend()); | |
195 | } | |
196 | ||
197 | template <class Backend, expression_template_option ExpressionTemplates> | |
198 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, unsigned>::type | |
199 | msb(const number<Backend, ExpressionTemplates>& x) | |
200 | { | |
201 | using default_ops::eval_msb; | |
202 | return eval_msb(x.backend()); | |
203 | } | |
204 | ||
205 | template <class tag, class A1, class A2, class A3, class A4> | |
206 | inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, unsigned>::type | |
207 | msb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x) | |
208 | { | |
209 | typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type; | |
210 | number_type n(x); | |
211 | using default_ops::eval_msb; | |
212 | return eval_msb(n.backend()); | |
213 | } | |
214 | ||
215 | template <class Backend, expression_template_option ExpressionTemplates> | |
216 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, bool>::type | |
217 | bit_test(const number<Backend, ExpressionTemplates>& x, unsigned index) | |
218 | { | |
219 | using default_ops::eval_bit_test; | |
220 | return eval_bit_test(x.backend(), index); | |
221 | } | |
222 | ||
223 | template <class tag, class A1, class A2, class A3, class A4> | |
224 | inline typename enable_if_c<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, bool>::type | |
225 | bit_test(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, unsigned index) | |
226 | { | |
227 | typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type; | |
228 | number_type n(x); | |
229 | using default_ops::eval_bit_test; | |
230 | return eval_bit_test(n.backend(), index); | |
231 | } | |
232 | ||
233 | template <class Backend, expression_template_option ExpressionTemplates> | |
234 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type | |
235 | bit_set(number<Backend, ExpressionTemplates>& x, unsigned index) | |
236 | { | |
237 | using default_ops::eval_bit_set; | |
238 | eval_bit_set(x.backend(), index); | |
239 | return x; | |
240 | } | |
241 | ||
242 | template <class Backend, expression_template_option ExpressionTemplates> | |
243 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type | |
244 | bit_unset(number<Backend, ExpressionTemplates>& x, unsigned index) | |
245 | { | |
246 | using default_ops::eval_bit_unset; | |
247 | eval_bit_unset(x.backend(), index); | |
248 | return x; | |
249 | } | |
250 | ||
251 | template <class Backend, expression_template_option ExpressionTemplates> | |
252 | inline typename enable_if_c<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type | |
253 | bit_flip(number<Backend, ExpressionTemplates>& x, unsigned index) | |
254 | { | |
255 | using default_ops::eval_bit_flip; | |
256 | eval_bit_flip(x.backend(), index); | |
257 | return x; | |
258 | } | |
259 | ||
260 | namespace default_ops{ | |
261 | ||
262 | // | |
263 | // Within powm, we need a type with twice as many digits as the argument type, define | |
264 | // a traits class to obtain that type: | |
265 | // | |
266 | template <class Backend> | |
267 | struct double_precision_type | |
268 | { | |
269 | typedef Backend type; | |
270 | }; | |
271 | ||
272 | // | |
273 | // If the exponent is a signed integer type, then we need to | |
274 | // check the value is positive: | |
275 | // | |
276 | template <class Backend> | |
277 | inline void check_sign_of_backend(const Backend& v, const mpl::true_) | |
278 | { | |
279 | if(eval_get_sign(v) < 0) | |
280 | { | |
281 | BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); | |
282 | } | |
283 | } | |
284 | template <class Backend> | |
285 | inline void check_sign_of_backend(const Backend&, const mpl::false_){} | |
286 | // | |
287 | // Calculate (a^p)%c: | |
288 | // | |
289 | template <class Backend> | |
290 | void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c) | |
291 | { | |
292 | using default_ops::eval_bit_test; | |
293 | using default_ops::eval_get_sign; | |
294 | using default_ops::eval_multiply; | |
295 | using default_ops::eval_modulus; | |
296 | using default_ops::eval_right_shift; | |
297 | ||
298 | typedef typename double_precision_type<Backend>::type double_type; | |
299 | typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type; | |
300 | ||
301 | check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>()); | |
302 | ||
303 | double_type x, y(a), b(p), t; | |
304 | x = ui_type(1u); | |
305 | ||
306 | while(eval_get_sign(b) > 0) | |
307 | { | |
308 | if(eval_bit_test(b, 0)) | |
309 | { | |
310 | eval_multiply(t, x, y); | |
311 | eval_modulus(x, t, c); | |
312 | } | |
313 | eval_multiply(t, y, y); | |
314 | eval_modulus(y, t, c); | |
315 | eval_right_shift(b, ui_type(1)); | |
316 | } | |
317 | Backend x2(x); | |
318 | eval_modulus(result, x2, c); | |
319 | } | |
320 | ||
321 | template <class Backend, class Integer> | |
322 | void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c) | |
323 | { | |
324 | typedef typename double_precision_type<Backend>::type double_type; | |
325 | typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type; | |
326 | typedef typename boost::multiprecision::detail::canonical<Integer, double_type>::type i1_type; | |
327 | typedef typename boost::multiprecision::detail::canonical<Integer, Backend>::type i2_type; | |
328 | ||
329 | using default_ops::eval_bit_test; | |
330 | using default_ops::eval_get_sign; | |
331 | using default_ops::eval_multiply; | |
332 | using default_ops::eval_modulus; | |
333 | using default_ops::eval_right_shift; | |
334 | ||
335 | check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>()); | |
336 | ||
337 | if(eval_get_sign(p) < 0) | |
338 | { | |
339 | BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); | |
340 | } | |
341 | ||
342 | double_type x, y(a), b(p), t; | |
343 | x = ui_type(1u); | |
344 | ||
345 | while(eval_get_sign(b) > 0) | |
346 | { | |
347 | if(eval_bit_test(b, 0)) | |
348 | { | |
349 | eval_multiply(t, x, y); | |
350 | eval_modulus(x, t, static_cast<i1_type>(c)); | |
351 | } | |
352 | eval_multiply(t, y, y); | |
353 | eval_modulus(y, t, static_cast<i1_type>(c)); | |
354 | eval_right_shift(b, ui_type(1)); | |
355 | } | |
356 | Backend x2(x); | |
357 | eval_modulus(result, x2, static_cast<i2_type>(c)); | |
358 | } | |
359 | ||
360 | template <class Backend, class Integer> | |
361 | typename enable_if<is_unsigned<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c) | |
362 | { | |
363 | typedef typename double_precision_type<Backend>::type double_type; | |
364 | typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type; | |
365 | ||
366 | using default_ops::eval_bit_test; | |
367 | using default_ops::eval_get_sign; | |
368 | using default_ops::eval_multiply; | |
369 | using default_ops::eval_modulus; | |
370 | using default_ops::eval_right_shift; | |
371 | ||
372 | double_type x, y(a), t; | |
373 | x = ui_type(1u); | |
374 | ||
375 | while(b > 0) | |
376 | { | |
377 | if(b & 1) | |
378 | { | |
379 | eval_multiply(t, x, y); | |
380 | eval_modulus(x, t, c); | |
381 | } | |
382 | eval_multiply(t, y, y); | |
383 | eval_modulus(y, t, c); | |
384 | b >>= 1; | |
385 | } | |
386 | Backend x2(x); | |
387 | eval_modulus(result, x2, c); | |
388 | } | |
389 | ||
390 | template <class Backend, class Integer> | |
391 | typename enable_if<is_signed<Integer> >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c) | |
392 | { | |
393 | if(b < 0) | |
394 | { | |
395 | BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); | |
396 | } | |
397 | eval_powm(result, a, static_cast<typename make_unsigned<Integer>::type>(b), c); | |
398 | } | |
399 | ||
400 | template <class Backend, class Integer1, class Integer2> | |
401 | typename enable_if<is_unsigned<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c) | |
402 | { | |
403 | typedef typename double_precision_type<Backend>::type double_type; | |
404 | typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type; | |
405 | typedef typename boost::multiprecision::detail::canonical<Integer1, double_type>::type i1_type; | |
406 | typedef typename boost::multiprecision::detail::canonical<Integer2, Backend>::type i2_type; | |
407 | ||
408 | using default_ops::eval_bit_test; | |
409 | using default_ops::eval_get_sign; | |
410 | using default_ops::eval_multiply; | |
411 | using default_ops::eval_modulus; | |
412 | using default_ops::eval_right_shift; | |
413 | ||
414 | double_type x, y(a), t; | |
415 | x = ui_type(1u); | |
416 | ||
417 | while(b > 0) | |
418 | { | |
419 | if(b & 1) | |
420 | { | |
421 | eval_multiply(t, x, y); | |
422 | eval_modulus(x, t, static_cast<i1_type>(c)); | |
423 | } | |
424 | eval_multiply(t, y, y); | |
425 | eval_modulus(y, t, static_cast<i1_type>(c)); | |
426 | b >>= 1; | |
427 | } | |
428 | Backend x2(x); | |
429 | eval_modulus(result, x2, static_cast<i2_type>(c)); | |
430 | } | |
431 | ||
432 | template <class Backend, class Integer1, class Integer2> | |
433 | typename enable_if<is_signed<Integer1> >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c) | |
434 | { | |
435 | if(b < 0) | |
436 | { | |
437 | BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); | |
438 | } | |
439 | eval_powm(result, a, static_cast<typename make_unsigned<Integer1>::type>(b), c); | |
440 | } | |
441 | ||
442 | struct powm_func | |
443 | { | |
444 | template <class T, class U, class V> | |
445 | void operator()(T& result, const T& b, const U& p, const V& m)const | |
446 | { | |
447 | eval_powm(result, b, p, m); | |
448 | } | |
449 | }; | |
450 | ||
451 | } | |
452 | ||
453 | template <class T, class U, class V> | |
454 | inline typename enable_if< | |
455 | mpl::and_< | |
456 | mpl::bool_<number_category<T>::value == number_kind_integer>, | |
457 | mpl::or_< | |
458 | is_number<T>, | |
459 | is_number_expression<T> | |
460 | >, | |
461 | mpl::or_< | |
462 | is_number<U>, | |
463 | is_number_expression<U>, | |
464 | is_integral<U> | |
465 | >, | |
466 | mpl::or_< | |
467 | is_number<V>, | |
468 | is_number_expression<V>, | |
469 | is_integral<V> | |
470 | > | |
471 | >, | |
472 | typename mpl::if_< | |
473 | is_no_et_number<T>, | |
474 | T, | |
475 | typename mpl::if_< | |
476 | is_no_et_number<U>, | |
477 | U, | |
478 | typename mpl::if_< | |
479 | is_no_et_number<V>, | |
480 | V, | |
481 | detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type | |
482 | >::type | |
483 | >::type | |
484 | >::type | |
485 | powm(const T& b, const U& p, const V& mod) | |
486 | { | |
487 | return detail::expression<detail::function, default_ops::powm_func, T, U, V>( | |
488 | default_ops::powm_func(), b, p, mod); | |
489 | } | |
490 | ||
491 | }} //namespaces | |
492 | ||
493 | #endif | |
494 | ||
495 |