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_
6 #ifndef BOOST_MP_INT_FUNC_HPP
7 #define BOOST_MP_INT_FUNC_HPP
9 #include <boost/multiprecision/number.hpp>
11 namespace boost{ namespace multiprecision{
16 template <class Backend>
17 inline void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r)
20 eval_modulus(r, x, y);
23 template <class Backend, class Integer>
24 inline Integer eval_integer_modulus(const Backend& x, Integer val)
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;
31 eval_modulus(t, x, static_cast<int_type>(val));
33 eval_convert_to(&result, t);
39 #pragma warning(disable:4127)
43 inline void eval_gcd(B& result, const B& a, const B& b)
45 using default_ops::eval_lsb;
46 using default_ops::eval_is_zero;
47 using default_ops::eval_get_sign;
53 int s = eval_get_sign(u);
76 /* Let shift := lg K, where K is the greatest power of 2
77 dividing both u and v. */
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);
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. */
96 eval_right_shift(v, vs);
101 eval_left_shift(result, shift);
109 inline void eval_lcm(B& result, const B& a, const B& b)
111 typedef typename mpl::front<typename B::unsigned_types>::type ui_type;
117 result = static_cast<ui_type>(0);
121 eval_divide(result, a, t);
122 eval_multiply(result, b);
124 if(eval_get_sign(result) < 0)
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)
135 using default_ops::eval_qr;
136 eval_qr(x.backend(), y.backend(), q.backend(), r.backend());
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)
144 divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r);
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)
152 divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r);
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)
160 divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r);
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)
167 using default_ops::eval_integer_modulus;
168 return eval_integer_modulus(x.backend(), val);
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)
175 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type result_type;
176 return integer_modulus(result_type(x), val);
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)
183 using default_ops::eval_lsb;
184 return eval_lsb(x.backend());
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)
191 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
193 using default_ops::eval_lsb;
194 return eval_lsb(n.backend());
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)
201 using default_ops::eval_msb;
202 return eval_msb(x.backend());
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)
209 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
211 using default_ops::eval_msb;
212 return eval_msb(n.backend());
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)
219 using default_ops::eval_bit_test;
220 return eval_bit_test(x.backend(), index);
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)
227 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
229 using default_ops::eval_bit_test;
230 return eval_bit_test(n.backend(), index);
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)
237 using default_ops::eval_bit_set;
238 eval_bit_set(x.backend(), index);
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)
246 using default_ops::eval_bit_unset;
247 eval_bit_unset(x.backend(), index);
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)
255 using default_ops::eval_bit_flip;
256 eval_bit_flip(x.backend(), index);
260 namespace default_ops{
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:
266 template <class Backend>
267 struct double_precision_type
269 typedef Backend type;
273 // If the exponent is a signed integer type, then we need to
274 // check the value is positive:
276 template <class Backend>
277 inline void check_sign_of_backend(const Backend& v, const mpl::true_)
279 if(eval_get_sign(v) < 0)
281 BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
284 template <class Backend>
285 inline void check_sign_of_backend(const Backend&, const mpl::false_){}
287 // Calculate (a^p)%c:
289 template <class Backend>
290 void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c)
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;
298 typedef typename double_precision_type<Backend>::type double_type;
299 typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
301 check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());
303 double_type x, y(a), b(p), t;
306 while(eval_get_sign(b) > 0)
308 if(eval_bit_test(b, 0))
310 eval_multiply(t, x, y);
311 eval_modulus(x, t, c);
313 eval_multiply(t, y, y);
314 eval_modulus(y, t, c);
315 eval_right_shift(b, ui_type(1));
318 eval_modulus(result, x2, c);
321 template <class Backend, class Integer>
322 void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c)
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;
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;
335 check_sign_of_backend(p, mpl::bool_<std::numeric_limits<number<Backend> >::is_signed>());
337 if(eval_get_sign(p) < 0)
339 BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
342 double_type x, y(a), b(p), t;
345 while(eval_get_sign(b) > 0)
347 if(eval_bit_test(b, 0))
349 eval_multiply(t, x, y);
350 eval_modulus(x, t, static_cast<i1_type>(c));
352 eval_multiply(t, y, y);
353 eval_modulus(y, t, static_cast<i1_type>(c));
354 eval_right_shift(b, ui_type(1));
357 eval_modulus(result, x2, static_cast<i2_type>(c));
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)
363 typedef typename double_precision_type<Backend>::type double_type;
364 typedef typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type ui_type;
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;
372 double_type x, y(a), t;
379 eval_multiply(t, x, y);
380 eval_modulus(x, t, c);
382 eval_multiply(t, y, y);
383 eval_modulus(y, t, c);
387 eval_modulus(result, x2, c);
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)
395 BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
397 eval_powm(result, a, static_cast<typename make_unsigned<Integer>::type>(b), c);
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)
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;
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;
414 double_type x, y(a), t;
421 eval_multiply(t, x, y);
422 eval_modulus(x, t, static_cast<i1_type>(c));
424 eval_multiply(t, y, y);
425 eval_modulus(y, t, static_cast<i1_type>(c));
429 eval_modulus(result, x2, static_cast<i2_type>(c));
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)
437 BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
439 eval_powm(result, a, static_cast<typename make_unsigned<Integer1>::type>(b), c);
444 template <class T, class U, class V>
445 void operator()(T& result, const T& b, const U& p, const V& m)const
447 eval_powm(result, b, p, m);
453 template <class T, class U, class V>
454 inline typename enable_if<
456 mpl::bool_<number_category<T>::value == number_kind_integer>,
459 is_number_expression<T>
463 is_number_expression<U>,
468 is_number_expression<V>,
481 detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type
485 powm(const T& b, const U& p, const V& mod)
487 return detail::expression<detail::function, default_ops::powm_func, T, U, V>(
488 default_ops::powm_func(), b, p, mod);