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