Loading...
Searching...
No Matches
mp_boost.cpp
1#include "mp_class.h"
2#include <boost/multiprecision/miller_rabin.hpp>
3#include <boost/mpl/int.hpp>
4#include <utility>
5#include <cmath>
6#include <symengine/symengine_assert.h>
7#include <symengine/symengine_exception.h>
8#include <symengine/prime_sieve.h>
9
10using boost::multiprecision::denominator;
11using boost::multiprecision::miller_rabin_test;
12using boost::multiprecision::numerator;
13using boost::multiprecision::detail::find_lsb;
14
15#if SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
16
17namespace SymEngine
18{
19
20integer_class pow(const integer_class &a, unsigned long b)
21{
22 return boost::multiprecision::pow(a, numeric_cast<unsigned>(b));
23}
24
25void mp_fdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
26 const integer_class &b)
27{
28 /*boost::multiprecision doesn't have a built-in fdiv_qr (floored division).
29 Its divide_qr uses truncated division, as does its
30 modulus operator. Thus, using boost::multiprecision::divide_qr we get:
31 divide_qr(-5, 3, quo, rem) //quo == -1, rem == -2
32 divide_qr(5, -3, quo, rem) //quo == -1, rem == 2
33 but we want:
34 mp_fdiv_r(quo, rem, -5, 3) //quo == -2, rem == 1
35 mp_fdiv_r(quo, rem, 5, -3) //rem == -2, rem == -1
36 The problem arises only when the quotient is negative. To convert
37 a truncated result into a floored result in this case, simply subtract
38 one from the truncated quotient and add the divisor to the truncated
39 remainder.
40 */
41
42 // must copy a and b before calling divide_qr because a or b may refer to
43 // the same
44 // object as q or r, causing incorrect results
45 integer_class a_cpy = a, b_cpy = b;
46 bool neg_quotient = ((a < 0 && b > 0) || (a > 0 && b < 0)) ? true : false;
47 boost::multiprecision::divide_qr(a_cpy, b_cpy, q, r);
48 // floor the quotient if necessary
49 if (neg_quotient && r != 0) {
50 q -= 1;
51 }
52 // remainder should have same sign as divisor
53 if ((b_cpy > 0 && r < 0) || (b_cpy < 0 && r > 0)) {
54 r += b_cpy;
55 return;
56 }
57}
58
59void mp_cdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
60 const integer_class &b)
61{
62 integer_class a_cpy = a, b_cpy = b;
63 bool pos_quotient = ((a < 0 && b < 0) || (a > 0 && b > 0)) ? true : false;
64 boost::multiprecision::divide_qr(a_cpy, b_cpy, q, r);
65 // ceil the quotient if necessary
66 if (pos_quotient && r != 0) {
67 q += 1;
68 }
69 // remainder should have opposite sign as divisor
70 if ((b_cpy > 0 && r > 0) || (b_cpy < 0 && r < 0)) {
71 r -= b_cpy;
72 return;
73 }
74}
75
76void mp_gcdext(integer_class &gcd, integer_class &s, integer_class &t,
77 const integer_class &a, const integer_class &b)
78{
79
80 integer_class this_s(1);
81 integer_class this_t(0);
82 integer_class next_s(0);
83 integer_class next_t(1);
84 integer_class this_r(a);
85 integer_class next_r(b);
86 integer_class q;
87 while (next_r != 0) {
88 // should use truncated division, so use
89 // boost::multiprecision::divide_qr
90 // beware of overwriting this_r during internal operations of divide_qr
91 // copy it first
92 integer_class this_r_cpy = this_r;
93 boost::multiprecision::divide_qr(this_r_cpy, next_r, q, this_r);
94 this_s -= q * next_s;
95 this_t -= q * next_t;
96 std::swap(this_s, next_s);
97 std::swap(this_t, next_t);
98 std::swap(this_r, next_r);
99 }
100 // normalize the gcd, s and t
101 if (this_r < 0) {
102 this_r *= -1;
103 this_s *= -1;
104 this_t *= -1;
105 }
106 gcd = std::move(this_r);
107 s = std::move(this_s);
108 t = std::move(this_t);
109}
110
111bool mp_invert(integer_class &res, const integer_class &a,
112 const integer_class &m)
113{
114 integer_class gcd, s, t;
115 mp_gcdext(gcd, s, t, a, m);
116 if (gcd != 1) {
117 res = 0;
118 return false;
119 } else {
120 mp_fdiv_r(s, s, m); // reduce s modulo m. undefined behavior when m ==
121 // 0, so don't need to check
122 if (s < 0) {
123 s += mp_abs(m);
124 } // give the canonical representative of s
125 res = s;
126 return true;
127 }
128}
129
130// floored modulus
131integer_class fmod(const integer_class &a, const integer_class &mod)
132{
133 integer_class res = a % mod;
134 if (res < 0) {
135 res += mod;
136 }
137 return res;
138}
139
140void mp_pow_ui(rational_class &res, const rational_class &i, unsigned long n)
141{
142 integer_class num = numerator(i); // copy
143 integer_class den = denominator(i); // copy
144 num = pow(num, n);
145 den = pow(den, n);
146 res = rational_class(std::move(num), std::move(den));
147}
148
149void mp_powm(integer_class &res, const integer_class &base,
150 const integer_class &exp, const integer_class &m)
151{
152 // if exp is negative, interpret as follows
153 // base**(exp) mod m == (base**(-1))**abs(exp) mod m
154 // == (base**(-1) mod m) ** abs(exp) mod m
155 // where base**(-1) mod m is the modular inverse
156 if (exp < 0) {
157 integer_class base_inverse;
158 if (!mp_invert(base_inverse, base, m)) {
159 throw SymEngine::SymEngineException("negative exponent undefined "
160 "in powm if base is not "
161 "invertible mod m");
162 }
163 res = boost::multiprecision::powm(base_inverse, mp_abs(exp), m);
164 return;
165 } else {
166 res = boost::multiprecision::powm(base, exp, m);
167 // boost's powm calculates base**exp % m, but uses truncated
168 // modulus, e.g. powm(-2,3,5) == -3. We want powm(-2,3,5) == 2
169 if (res < 0) {
170 res += m;
171 }
172 }
173}
174
175integer_class step(const unsigned long &n, const integer_class &i,
176 integer_class &x)
177{
178 SYMENGINE_ASSERT(n > 1);
179 unsigned long m = n - 1;
180 integer_class &&x_m = pow(x, m);
181 return integer_class((integer_class(m * x) + integer_class(i / x_m)) / n);
182}
183
184bool positive_root(integer_class &res, const integer_class &i,
185 const unsigned long n)
186{
187 integer_class x
188 = 1; // TODO: make a better starting guess based on (number of bits)/n
189 integer_class y = step(n, i, x);
190 do {
191 x = y;
192 y = step(n, i, x);
193 } while (y < x);
194 res = x;
195 if (pow(x, n) == i) {
196 return true;
197 }
198 return false;
199}
200
201// return true if i is a perfect nth power, i.e. res**i == n
202bool mp_root(integer_class &res, const integer_class &i, const unsigned long n)
203{
204 if (n == 0) {
205 throw std::runtime_error("0th root is undefined");
206 }
207 if (n == 1) {
208 res = i;
209 return true;
210 }
211 if (i == 0) {
212 res = 0;
213 return true;
214 }
215 if (i > 0) {
216 return positive_root(res, i, n);
217 }
218 if (i < 0 && (n % 2 == 0)) {
219 throw std::runtime_error("even root of a negative is non-real");
220 }
221 bool b = positive_root(res, -i, n);
222 res *= -1;
223 return b;
224}
225
226integer_class mp_sqrt(const integer_class &i)
227{
228 // as of 11/1/2016, boost::multiprecision::sqrt() is buggy:
229 // https://svn.boost.org/trac/boost/ticket/12559
230 // implement with mp_root for now
231 integer_class res;
232 mp_root(res, i, 2);
233 return res;
234}
235
236void mp_rootrem(integer_class &a, integer_class &b, const integer_class &i,
237 unsigned long n)
238{
239 mp_root(a, i, n);
240 integer_class p = pow(a, n);
241 ;
242 b = i - p;
243}
244
245void mp_sqrtrem(integer_class &a, integer_class &b, const integer_class &i)
246{
247 a = mp_sqrt(i);
248 b = i - boost::multiprecision::pow(a, 2);
249}
250
251// return nonzero if i is probably prime.
252int mp_probab_prime_p(const integer_class &i, unsigned retries)
253{
254 if (i % 2 == 0)
255 return (i == 2);
256 return miller_rabin_test(i, retries);
257}
258
259void mp_nextprime(integer_class &res, const integer_class &i)
260{
261 // simple implementation: just check all odds bigger than i for primality
262 if (i < 2) {
263 res = 2;
264 return;
265 }
266 integer_class candidate;
267 candidate = (i % 2 == 0) ? i + 1 : i + 2;
268 // Knuth recommends 25 trials for a pretty strong likelihood that candidate
269 // is prime
270 while (!mp_probab_prime_p(candidate, 25)) {
271 candidate += 2;
272 }
273 res = std::move(candidate);
274}
275
276unsigned long mp_scan1(const integer_class &i)
277{
278 if (i == 0) {
279 return ULONG_MAX;
280 }
281 return find_lsb(i, {});
282}
283
284// define simple 2x2 matrix with exponentiation by repeated squaring
285// to use in logarithmic-time fibonacci calculation
286
288 integer_class data[2][2]; // data[1][0] is row 1, column 0 entry of matrix
289
290 two_by_two_matrix(integer_class a, integer_class b, integer_class c,
291 integer_class d)
292 : data{{a, b}, {c, d}}
293 {
294 }
295 two_by_two_matrix() : data{{0, 0}, {0, 0}} {}
296 two_by_two_matrix(const two_by_two_matrix &other) = default;
297 two_by_two_matrix &operator=(const two_by_two_matrix &other)
298 {
299 this->data[0][0] = other.data[0][0];
300 this->data[0][1] = other.data[0][1];
301 this->data[1][0] = other.data[1][0];
302 this->data[1][1] = other.data[1][1];
303 return *this;
304 }
305 two_by_two_matrix operator*(const two_by_two_matrix &other)
306 {
308 res.data[0][0] = this->data[0][0] * other.data[0][0]
309 + this->data[0][1] * other.data[1][0];
310 res.data[0][1] = this->data[0][0] * other.data[0][1]
311 + this->data[0][1] * other.data[1][1];
312 res.data[1][0] = this->data[1][0] * other.data[0][0]
313 + this->data[1][1] * other.data[1][0];
314 res.data[1][1] = this->data[1][0] * other.data[0][1]
315 + this->data[1][1] * other.data[1][1];
316 return res;
317 }
318 // recursive repeated squaring
319 two_by_two_matrix pow(unsigned long n)
320 {
321 if (n == 0) {
322 return two_by_two_matrix::identity();
323 }
324 if (n == 1) {
325 return *this;
326 }
327 if (n == 2) {
328 return (*this) * (*this);
329 }
330 if (n % 2 == 0) {
331 return (this->pow(n / 2)).pow(2);
332 }
333 return ((this->pow((n - 1) / 2)).pow(2)) * (*this);
334 }
335
336 static two_by_two_matrix identity()
337 {
338 return two_by_two_matrix(1, 0, 0, 1);
339 }
340};
341
342inline two_by_two_matrix fib_matrix(unsigned long n)
343{
344 two_by_two_matrix x(1, 1, 1, 0);
345 return x.pow(n);
346}
347
348void mp_fib_ui(integer_class &res, unsigned long n)
349{
350 // reference: https://www.nayuki.io/page/fast-fibonacci-algorithms
351 res = fib_matrix(n).data[0][1];
352}
353
354// sets a = Fibonacci(n) and b = Fibonacci(n-1)
355void mp_fib2_ui(integer_class &a, integer_class &b, unsigned long n)
356{
357 // reference: https://www.nayuki.io/page/fast-fibonacci-algorithms
358 two_by_two_matrix result_matrix = fib_matrix(n);
359 a = result_matrix.data[0][1];
360 b = result_matrix.data[1][1];
361}
362
363inline two_by_two_matrix luc_matrix(unsigned long n)
364{
365 two_by_two_matrix multiplier(1, 1, 1, 0);
366 two_by_two_matrix start(1, 0, 2, 0);
367 return multiplier.pow(n) * start;
368}
369
370void mp_lucnum_ui(integer_class &res, unsigned long n)
371{
372 // implementation based on the following fact:
373 // [[1,1],[1,0]]^(n-1)*[2,0,1,0] = [[L(n+1),0][L(n),0]]
374 // where L(n) is the nth Lucas number
375 res = luc_matrix(n).data[1][0];
376}
377
378void mp_lucnum2_ui(integer_class &a, integer_class &b, unsigned long n)
379{
380 if (n == 0) {
381 throw std::runtime_error("index of lucas number cannot be negative");
382 }
383 two_by_two_matrix result_matrix = luc_matrix(n - 1);
384 a = result_matrix.data[0][0];
385 b = result_matrix.data[1][0];
386}
387
388void mp_fac_ui(integer_class &res, unsigned long n)
389{
390 // couldn't make boost's template version of factorial work,
391 // so implement slow, naive version for now
392 res = 1;
393 for (unsigned long i = 2; i <= n; ++i) {
394 res *= i;
395 }
396}
397
398void mp_bin_ui(integer_class &res, const integer_class &n, unsigned long r)
399{
400 // slow, naive implementation
401 integer_class x = n - r;
402 res = 1;
403 for (unsigned long i = 1; i <= r; ++i) {
404 res *= x + i;
405 res /= i;
406 }
407}
408
409// this is extremely slow!
410bool mp_perfect_power_p(const integer_class &i)
411{
412 if (i == 0 || i == 1 || i == -1) {
413 return true;
414 }
415 // if i == a**k, with k == pq for some integers p,q
416 // then i == a**(pq) == (a**p)**q == (a**q)**p
417 // hence i is a pth power and a qth power
418 // Hence if i == p**k, then i is a pth power for any p
419 // in the prime factorization of k, and it suffices
420 // to check whether i is a prime power.
421
422 // the largest possible prime p would arise from the
423 // case where i is a prime power of two
424 // so check all prime roots up to log(i) base 2.
425
426 unsigned long max = std::ilogb(i.convert_to<double>());
427
428 // treat case p=2 separately b/c mp_root throws exception
429 // with an even root of a negative
430 if (mp_perfect_square_p(i)) {
431 return true;
432 }
433 integer_class p(2);
434 integer_class root(0);
435 while (true) {
436 mp_nextprime(p, p);
437 if (p > max) {
438 return false;
439 }
440 if (mp_root(root, i, p.convert_to<unsigned long>())) {
441 return true;
442 }
443 }
444}
445
446bool mp_perfect_square_p(const integer_class &i)
447{
448 if (i < 0) {
449 return false;
450 }
451 integer_class root;
452 return mp_root(root, i, 2);
453}
454
455integer_class mp_primorial(unsigned long n)
456{
457 integer_class res = 1;
458 Sieve::iterator pi(static_cast<unsigned>(n));
459 unsigned int p;
460 while ((p = pi.next_prime()) <= n) {
461 res *= p;
462 }
463 return res;
464}
465
466// according to the gmp documentation, the behavior of the
467// corresponding function mpz_legendre is
468// undefined if n is not a positive odd prime.
469// hence we treat n as though it were a positive odd prime
470// but it is up to the caller to very this.
471int mp_legendre(const integer_class &a, const integer_class &n)
472{
473 integer_class res;
474 mp_powm(res, a, integer_class((n - 1) / 2), n);
475 return res <= 1 ? res.convert_to<int>() : -1;
476}
477
478// private function that computes jacobi symbols
479// without checking that arguments satisfy a >= 0
480// and n is odd.
481int unchecked_jacobi(const integer_class &a, const integer_class &n)
482{
483 // https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol
484 if (a == 1) {
485 return 1;
486 }
487 integer_class num = a;
488 integer_class den = n;
489 // (1) Reduce the "numerator" modulo the "denominator"
490 num = fmod(num, den);
491
492 // (2) Extract any factors of 2 from the "numerator"
493 unsigned long factors_of_two = 0;
494 while (num % 2 == 0 && num != 0) {
495 num /= 2; // use a shift instead of division here? faster?
496 ++factors_of_two;
497 }
498 int product_of_twos = 1; // (2 | den)**factors_of_two
499
500 // (2 | den) is -1 iff den % 8 = 3 or den % 8 = 5.
501 // If factors_of_two is odd and (2 | den) is -1,
502 // then (2 | den)**factors_of_two is -1.
503 // Otherwise, (2 | den)**factors_of_two is 1.
504 int den_mod_8 = fmod(den, 8).convert_to<int>();
505 if ((factors_of_two % 2 == 1) && (den_mod_8 == 3 || den_mod_8 == 5)) {
506 product_of_twos = -1;
507 }
508
509 // (3) If the remaining "numerator" (after extraction of twos) is 1,
510 // then the remaining jacobi symbol is 1.
511 // If the "numerator" and "denominator" are not coprime, result is 0.
512 if (num == 1) {
513 return product_of_twos;
514 }
515 if (boost::multiprecision::gcd(num, den) != 1) {
516 return 0;
517 }
518
519 // (4) Otherwise, the "numerator" and "denominator" are now odd
520 // positive coprime integers, so we can flip the symbol
521 // with quadratic reciprocity, then return to step (1)
522
523 // (num | den) == (den | num), unless num % 4 and den %4 are both
524 // 3, in which case (num | den) == (-1) * (den | num)
525 int quadratic_reciprocity_factor = 1;
526 if ((fmod(num, 4) == 3) && (fmod(den, 4) == 3)) {
527 quadratic_reciprocity_factor = -1;
528 }
529 return product_of_twos * quadratic_reciprocity_factor
530 * unchecked_jacobi(den, num);
531}
532
533// public interface for computing jacobi symbols. performs checking.
534int mp_jacobi(const integer_class &a, const integer_class &n)
535{
536 if (n < 0) {
537 throw std::runtime_error("jacobi denominator must be positive");
538 }
539 if (n % 2 == 0) {
540 throw std::runtime_error("jacobi denominator must be odd");
541 }
542 return unchecked_jacobi(a, n);
543}
544
545int mp_kronecker(const integer_class &a, const integer_class &n)
546{
547 /*
548 https://en.wikipedia.org/wiki/Kronecker_symbol
549 We compute the Kronecker symbol in terms of the Jacobi symbol.
550 For an integer n!=0, let n = u * PRODUCT (p_i**e_i), where u = -1 or 1 and
551 the p_i's and e_i's are the primes and corresponding powers
552 in the prime factorization of |n|.
553 Then for an integer a, the Kronecker symbol is given by
554 (a | n) = (a | u) * PRODUCT [(a | p_i)**e_i]
555 where
556 (a | u) is -1 if u is -1 and a < 0. Otherwise it is 1.
557 (a | 2) is 0 if a is even, 1 if (a mod 8) == 1 or 7, and -1 if (a mod 8) ==
558 3 or 5.
559 (a | p) is the Legendre symbol if p is an odd prime.
560
561 Let j be the power of the prime 2 in n's prime factorization, and let
562 m = |n|/(2**j) -- that is, m is |n| with all factors of 2 extracted.
563
564 Notice that, if p_1 is 2, then
565 PRODUCT [(a | p_i)**e_i]
566 == (a | 2)**j * PRODUCT [(a | p_i)**e_i] here the p_i's are the remaining
567 (odd) prime factors
568 == (a | 2)**j * (a | m), here (a | m) is the Jacobi
569 symbol, because m>0 is odd
570
571 Thus,
572 (a | n) == (a | u) * (a | 2)**j * (a | m) if n is even
573 (a | n) == (a | u) * (a | m) if n is odd
574 */
575
576 if (n == 0) {
577 throw std::runtime_error("second arg of Kronecker cannot be zero");
578 }
579
580 // Compute (a | u)
581 int kr_a_u = 1;
582 if (n.sign() == -1 && a < 0) {
583 kr_a_u = -1;
584 }
585
586 // Compute m, j
587 integer_class m = boost::multiprecision::abs(n);
588 unsigned long j = 0;
589 while (m % 2 == 0 && m != 0) { // while m is even
590 m /= 2; // implement with shift? faster?
591 ++j;
592 }
593
594 // if n is even, compute (a | 2)**j
595 int kr_a_2 = 0;
596 int kr_a_2_to_j = 0;
597 int a_mod_8 = fmod(a, 8).convert_to<int>();
598
599 if (a % 2 != 0) {
600 kr_a_2 = (a_mod_8 == 1 || a_mod_8 == 7) ? 1 : -1;
601 kr_a_2_to_j = (kr_a_2 == -1 && (j % 2 != 0))
602 ? -1
603 : 1; //(-1)**odd == -1; (-1)**even=(1)**integer=1
604 }
605
606 if (n % 2 == 0) {
607 return kr_a_u * kr_a_2_to_j * unchecked_jacobi(a, m);
608 } else {
609 return kr_a_u * unchecked_jacobi(a, m);
610 }
611}
612
613} // namespace SymEngine
614
615#endif // SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
T fmod(T... args)
T ilogb(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const Integer > gcd(const Integer &a, const Integer &b)
Greatest Common Divisor.
Definition: ntheory.cpp:31
RCP< const Basic > max(const vec_basic &arg)
Canonicalize Max:
Definition: functions.cpp:3555
RCP< const Integer > mod(const Integer &n, const Integer &d)
modulo round toward zero
Definition: ntheory.cpp:66
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:271
T pow(T... args)
T swap(T... args)