Loading...
Searching...
No Matches
ntheory.cpp
1#include <valarray>
2#include <iterator>
3
4#include <symengine/prime_sieve.h>
5#include <symengine/ntheory.h>
7#include <symengine/add.h>
8#include <symengine/mul.h>
9#include <symengine/pow.h>
10#ifdef HAVE_SYMENGINE_ECM
11#include <ecm.h>
12#endif // HAVE_SYMENGINE_ECM
13#ifdef HAVE_SYMENGINE_PRIMESIEVE
14#include <primesieve.hpp>
15#endif // HAVE_SYMENGINE_PRIMESIEVE
16#ifdef HAVE_SYMENGINE_ARB
17#include "arb.h"
18#include "bernoulli.h"
19#include "rational.h"
20#endif // HAVE_SYMENGINE_ARB
21#ifndef HAVE_SYMENGINE_GMP
22#include <boost/random/uniform_int.hpp>
23#include <boost/random/mersenne_twister.hpp>
24#include <boost/random.hpp>
25#endif // !HAVE_SYMENGINE_GMP
26
27namespace SymEngine
28{
29
30// Basic number theoretic functions
32{
33 integer_class g;
34 mp_gcd(g, a.as_integer_class(), b.as_integer_class());
35 return integer(std::move(g));
36}
37
38void gcd_ext(const Ptr<RCP<const Integer>> &g, const Ptr<RCP<const Integer>> &s,
39 const Ptr<RCP<const Integer>> &t, const Integer &a,
40 const Integer &b)
41{
42 integer_class g_, s_, t_;
43 mp_gcdext(g_, s_, t_, a.as_integer_class(), b.as_integer_class());
44 *g = integer(std::move(g_));
45 *s = integer(std::move(s_));
46 *t = integer(std::move(t_));
47}
48
50{
51 integer_class c;
52 mp_lcm(c, a.as_integer_class(), b.as_integer_class());
53 return integer(std::move(c));
54}
55
56int mod_inverse(const Ptr<RCP<const Integer>> &b, const Integer &a,
57 const Integer &m)
58{
59 int ret_val;
60 integer_class inv_t;
61 ret_val = mp_invert(inv_t, a.as_integer_class(), m.as_integer_class());
62 *b = integer(std::move(inv_t));
63 return ret_val;
64}
65
67{
68 return integer(n.as_integer_class() % d.as_integer_class());
69}
70
72{
73 return integer(n.as_integer_class() / d.as_integer_class());
74}
75
77 const Ptr<RCP<const Integer>> &r, const Integer &n,
78 const Integer &d)
79{
80 integer_class _q, _r;
81 mp_tdiv_qr(_q, _r, n.as_integer_class(), d.as_integer_class());
82 *q = integer(std::move(_q));
83 *r = integer(std::move(_r));
84}
85
87{
88 integer_class q;
89 mp_fdiv_r(q, n.as_integer_class(), d.as_integer_class());
90 return integer(std::move(q));
91}
92
94{
95 integer_class q;
96 mp_fdiv_q(q, n.as_integer_class(), d.as_integer_class());
97 return integer(std::move(q));
98}
99
101 const Ptr<RCP<const Integer>> &r, const Integer &n,
102 const Integer &d)
103{
104 integer_class _q, _r;
105 mp_fdiv_qr(_q, _r, n.as_integer_class(), d.as_integer_class());
106 *q = integer(std::move(_q));
107 *r = integer(std::move(_r));
108}
109
111{
112 integer_class f;
113 mp_fib_ui(f, n);
114 return integer(std::move(f));
115}
116
118 const Ptr<RCP<const Integer>> &s, unsigned long n)
119{
120 integer_class g_t;
121 integer_class s_t;
122 mp_fib2_ui(g_t, s_t, n);
123 *g = integer(std::move(g_t));
124 *s = integer(std::move(s_t));
125}
126
127RCP<const Integer> lucas(unsigned long n)
128{
129 integer_class f;
130 mp_lucnum_ui(f, n);
131 return integer(std::move(f));
132}
133
134void lucas2(const Ptr<RCP<const Integer>> &g, const Ptr<RCP<const Integer>> &s,
135 unsigned long n)
136{
137 integer_class g_t;
138 integer_class s_t;
139 mp_lucnum2_ui(g_t, s_t, n);
140 *g = integer(std::move(g_t));
141 *s = integer(std::move(s_t));
142}
143
144// Binomial Coefficient
145RCP<const Integer> binomial(const Integer &n, unsigned long k)
146{
147 integer_class f;
148 mp_bin_ui(f, n.as_integer_class(), k);
149 return integer(std::move(f));
150}
151
152// Factorial
154{
155 integer_class f;
156 mp_fac_ui(f, n);
157 return integer(std::move(f));
158}
159
160// Returns true if `b` divides `a` without reminder
161bool divides(const Integer &a, const Integer &b)
162{
163 return mp_divisible_p(a.as_integer_class(), b.as_integer_class()) != 0;
164}
165
166// Prime functions
167int probab_prime_p(const Integer &a, unsigned reps)
168{
169 return mp_probab_prime_p(a.as_integer_class(), reps);
170}
171
173{
174 integer_class c;
175 mp_nextprime(c, a.as_integer_class());
176 return integer(std::move(c));
177}
178
179namespace
180{
181// Factoring by Trial division using primes only
182int _factor_trial_division_sieve(integer_class &factor, const integer_class &N)
183{
184 integer_class sqrtN = mp_sqrt(N);
185 unsigned long limit = mp_get_ui(sqrtN);
187 throw SymEngineException("N too large to factor");
188 Sieve::iterator pi(numeric_cast<unsigned>(limit));
189 unsigned p;
190 while ((p = pi.next_prime()) <= limit) {
191 if (N % p == 0) {
192 factor = p;
193 return 1;
194 }
195 }
196 return 0;
197}
198// Factor using lehman method.
199int _factor_lehman_method(integer_class &rop, const integer_class &n)
200{
201 if (n < 21)
202 throw SymEngineException("Require n >= 21 to use lehman method");
203
204 int ret_val = 0;
205 integer_class u_bound;
206
207 mp_root(u_bound, n, 3);
208 u_bound = u_bound + 1;
209
210 Sieve::iterator pi(numeric_cast<unsigned>(mp_get_ui(u_bound)));
211 unsigned p;
212 while ((p = pi.next_prime()) <= mp_get_ui(u_bound)) {
213 if (n % p == 0) {
214 rop = n / p;
215 ret_val = 1;
216 break;
217 }
218 }
219
220 if (not ret_val) {
221
222 integer_class k, a, b, l;
223
224 k = 1;
225
226 while (k <= u_bound) {
227 a = mp_sqrt(4 * k * n);
228 mp_root(b, n, 6);
229 mp_root(l, k, 2);
230 b = b / (4 * l);
231 b = b + a;
232
233 while (a <= b) {
234 l = a * a - 4 * k * n;
235 if (mp_perfect_square_p(l)) {
236 b = a + mp_sqrt(l);
237 mp_gcd(rop, n, b);
238 ret_val = 1;
239 break;
240 }
241 a = a + 1;
242 }
243 if (ret_val)
244 break;
245 k = k + 1;
246 }
247 }
248
249 return ret_val;
250}
251} // anonymous namespace
252
254{
255 int ret_val;
256 integer_class rop;
257
259 *f = integer(std::move(rop));
260 return ret_val;
261}
262
263namespace
264{
265// Factor using Pollard's p-1 method
266int _factor_pollard_pm1_method(integer_class &rop, const integer_class &n,
267 const integer_class &c, unsigned B)
268{
269 if (n < 4 or B < 3)
270 throw SymEngineException(
271 "Require n > 3 and B > 2 to use Pollard's p-1 method");
272
273 integer_class m, _c;
274 _c = c;
275
276 Sieve::iterator pi(B);
277 unsigned p;
278 while ((p = pi.next_prime()) <= B) {
279 m = 1;
280 // calculate log(p, B), this can be improved
281 while (m <= B / p) {
282 m = m * p;
283 }
284 mp_powm(_c, _c, m, n);
285 }
286 _c = _c - 1;
287 mp_gcd(rop, _c, n);
288
289 if (rop == 1 or rop == n)
290 return 0;
291 else
292 return 1;
293}
294} // anonymous namespace
295
297 const Integer &n, unsigned B, unsigned retries)
298{
299 int ret_val = 0;
300 integer_class rop, nm4, c;
301
302 mp_randstate state;
303 nm4 = n.as_integer_class() - 4;
304
305 for (unsigned i = 0; i < retries and ret_val == 0; ++i) {
306 state.urandomint(c, nm4);
307 c += 2;
309 }
310
311 if (ret_val != 0)
312 *f = integer(std::move(rop));
313 return ret_val;
314}
315
316namespace
317{
318// Factor using Pollard's rho method
319int _factor_pollard_rho_method(integer_class &rop, const integer_class &n,
320 const integer_class &a, const integer_class &s,
321 unsigned steps = 10000)
322{
323 if (n < 5)
324 throw SymEngineException("Require n > 4 to use pollard's-rho method");
325
326 integer_class u, v, g, m;
327 u = s;
328 v = s;
329
330 for (unsigned i = 0; i < steps; ++i) {
331 u = (u * u + a) % n;
332 v = (v * v + a) % n;
333 v = (v * v + a) % n;
334 m = u - v;
335 mp_gcd(g, m, n);
336
337 if (g == n)
338 return 0;
339 if (g == 1)
340 continue;
341 rop = g;
342 return 1;
343 }
344 return 0;
345}
346} // namespace
347
349 const Integer &n, unsigned retries)
350{
351 int ret_val = 0;
352 integer_class rop, nm1, nm4, a, s;
353 mp_randstate state;
354 nm1 = n.as_integer_class() - 1;
355 nm4 = n.as_integer_class() - 4;
356
357 for (unsigned i = 0; i < retries and ret_val == 0; ++i) {
358 state.urandomint(a, nm1);
359 state.urandomint(s, nm4);
360 s += 1;
362 }
363
364 if (ret_val != 0)
365 *f = integer(std::move(rop));
366 return ret_val;
367}
368
369// Factorization
370int factor(const Ptr<RCP<const Integer>> &f, const Integer &n, double B1)
371{
372 int ret_val = 0;
373 integer_class _n, _f;
374
375 _n = n.as_integer_class();
376
377#ifdef HAVE_SYMENGINE_ECM
378 if (mp_perfect_power_p(_n)) {
379
380 unsigned long int i = 1;
381 integer_class m, rem;
382 rem = 1; // Any non zero number
383 m = 2; // set `m` to 2**i, i = 1 at the begining
384
385 // calculate log2n, this can be improved
386 for (; m < _n; ++i)
387 m = m * 2;
388
389 // eventually `rem` = 0 zero as `n` is a perfect power. `f_t` will
390 // be set to a factor of `n` when that happens
391 while (i > 1 and rem != 0) {
392 mp_rootrem(_f, rem, _n, i);
393 --i;
394 }
395
396 ret_val = 1;
397 } else {
398
399 if (mp_probab_prime_p(_n, 25) > 0) { // most probably, n is a prime
400 ret_val = 0;
401 _f = _n;
402 } else {
403
404 for (int i = 0; i < 10 and not ret_val; ++i)
405 ret_val = ecm_factor(get_mpz_t(_f), get_mpz_t(_n), B1, nullptr);
406 mp_demote(_f);
407 if (not ret_val)
408 throw SymEngineException(
409 "ECM failed to factor the given number");
410 }
411 }
412#else
413 // B1 is discarded if gmp-ecm is not installed
415#endif // HAVE_SYMENGINE_ECM
416 *f = integer(std::move(_f));
417
418 return ret_val;
419}
420
422{
423 int ret_val;
424 integer_class factor;
426 if (ret_val == 1)
428 return ret_val;
429}
430
432 const Integer &n)
433{
434 integer_class sqrtN;
435 integer_class _n = n.as_integer_class();
436 if (_n == 0)
437 return;
438 if (_n < 0)
439 _n *= -1;
440
441 sqrtN = mp_sqrt(_n);
442 auto limit = mp_get_ui(sqrtN);
443 if (not mp_fits_ulong_p(sqrtN)
445 throw SymEngineException("N too large to factor");
447 unsigned p;
448
449 while ((p = pi.next_prime()) <= limit) {
450 while (_n % p == 0) {
451 prime_list.push_back(integer(p));
452 _n = _n / p;
453 }
454 if (_n == 1)
455 break;
456 }
457 if (not(_n == 1))
458 prime_list.push_back(integer(std::move(_n)));
459}
460
462{
463 integer_class sqrtN;
464 integer_class _n = n.as_integer_class();
465 unsigned count;
466 if (_n == 0)
467 return;
468 if (_n < 0)
469 _n *= -1;
470
471 sqrtN = mp_sqrt(_n);
472 auto limit = mp_get_ui(sqrtN);
473 if (not mp_fits_ulong_p(sqrtN)
475 throw SymEngineException("N too large to factor");
477
478 unsigned p;
479 while ((p = pi.next_prime()) <= limit) {
480 count = 0;
481 while (_n % p == 0) { // when a prime factor is found, we divide
482 ++count; // _n by that prime as much as we can
483 _n = _n / p;
484 }
485 if (count > 0) {
486 insert(primes_mul, integer(p), count);
487 if (_n == 1)
488 break;
489 }
490 }
491 if (not(_n == 1))
493}
494
495RCP<const Number> bernoulli(unsigned long n)
496{
497#ifdef HAVE_SYMENGINE_ARB
498 fmpq_t res;
499 fmpq_init(res);
501 mpq_t a;
502 mpq_init(a);
503 fmpq_get_mpq(a, res);
504 rational_class b(a);
506 mpq_clear(a);
507 return Rational::from_mpq(std::move(b));
508#else
509 // TODO: implement a faster algorithm
511 for (unsigned m = 0; m <= n; ++m) {
512 v[m] = rational_class(1u, m + 1);
513
514 for (unsigned j = m; j >= 1; --j) {
515 v[j - 1] = j * (v[j - 1] - v[j]);
516 }
517 }
518 return Rational::from_mpq(v[0]);
519#endif
520}
521
522RCP<const Number> harmonic(unsigned long n, long m)
523{
524 rational_class res(0);
525 if (m == 1) {
526 for (unsigned i = 1; i <= n; ++i) {
527 res += rational_class(1u, i);
528 }
529 return Rational::from_mpq(res);
530 } else {
531 for (unsigned i = 1; i <= n; ++i) {
532 if (m > 0) {
533 rational_class t(1u, i);
534#if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP
535 mp_pow_ui(get_den(t), get_den(t), m);
536#else
537 mp_pow_ui(t, t, m);
538#endif
539 res += t;
540 } else {
541 integer_class t(i);
542 mp_pow_ui(t, t, static_cast<unsigned long>(-m));
543 res += t;
544 }
545 }
546 return Rational::from_mpq(res);
547 }
548}
549
550// References : Cohen H., A course in computational algebraic number theory
551// (1996), page 21.
552bool crt(const Ptr<RCP<const Integer>> &R,
555{
556 if (mod.size() > rem.size())
557 throw SymEngineException("Too few remainders");
558 if (mod.size() == 0)
559 throw SymEngineException("Moduli vector cannot be empty");
560
561 integer_class m, r, g, s, t;
562 m = mod[0]->as_integer_class();
563 r = rem[0]->as_integer_class();
564
565 for (unsigned i = 1; i < mod.size(); ++i) {
566 mp_gcdext(g, s, t, m, mod[i]->as_integer_class());
567 // g = s * m + t * mod[i]
568 t = rem[i]->as_integer_class() - r;
569 if (not mp_divisible_p(t, g))
570 return false;
571 r += m * s * (t / g); // r += m * (m**-1 mod[i]/g)* (rem[i] - r) / g
572 m *= mod[i]->as_integer_class() / g;
573 mp_fdiv_r(r, r, m);
574 }
575 *R = integer(std::move(r));
576 return true;
577}
578
579namespace
580{
581// Crt over a cartesian product of vectors (Assuming that moduli are pairwise
582// relatively prime).
586{
587 if (mod.size() > rem.size())
588 throw SymEngineException("Too few remainders");
589 if (mod.size() == 0)
590 throw SymEngineException("Moduli vector cannot be empty");
591 integer_class m, _m, r, s, t;
592 m = mod[0]->as_integer_class();
593 R = rem[0];
594
595 for (unsigned i = 1; i < mod.size(); ++i) {
597 mp_invert(s, m, mod[i]->as_integer_class());
598 _m = m;
599 m *= mod[i]->as_integer_class();
600 for (auto &elem : R) {
601 for (auto &_k : rem[i]) {
602 r = elem->as_integer_class();
603 r += _m * s * (_k->as_integer_class() - r);
604 mp_fdiv_r(r, r, m);
605 rem2.push_back(integer(r));
606 }
607 }
608 R = rem2;
609 }
610}
611
612// Tests whether n is a prime power and finds a prime p and e such that n =
613// p**e.
614bool _prime_power(integer_class &p, integer_class &e, const integer_class &n)
615{
616 if (n < 2)
617 return false;
618 integer_class _n = n, temp;
619 e = 1;
620 unsigned i = 2;
621 while (mp_perfect_power_p(_n) and _n >= 2) {
622 if (mp_root(temp, _n, i)) {
623 e *= i;
624 _n = temp;
625 } else {
626 ++i;
627 }
628 }
629 if (mp_probab_prime_p(_n, 25)) {
630 p = _n;
631 return true;
632 }
633 return false;
634}
635
636// Computes a primitive root modulo p**e or 2*p**e where p is an odd prime.
637// References : Cohen H., A course in computational algebraic number theory
638// (2009), pages 25-27.
639void _primitive_root(integer_class &g, const integer_class &p,
640 const integer_class &e, bool even = false)
641{
643 prime_factors(primes, *integer(p - 1));
644
645 integer_class t;
646 g = 2;
647 while (g < p) {
648 bool root = true;
649 for (const auto &it : primes) {
650 t = it->as_integer_class();
651 t = (p - 1) / t;
652 mp_powm(t, g, t, p);
653 if (t == 1) { // If g**(p-1)/q is 1 then g is not a primitive root.
654 root = false;
655 break;
656 }
657 }
658 if (root)
659 break;
660 ++g;
661 }
662
663 if (e > 1) {
664 t = p * p;
665 integer_class pm1 = p - 1;
666 mp_powm(t, g, pm1, t);
667 if (t == 1) { // If g**(p-1) mod (p**2) == 1 then g + p is a primitive
668 // root.
669 g += p;
670 }
671 }
672 if (even and g % 2 == 0) {
673 mp_pow_ui(t, p, mp_get_ui(e));
674 g += t; // If g is even then root of 2*p**e is g + p**e.
675 }
676}
677
678} // anonymous namespace
679
680bool primitive_root(const Ptr<RCP<const Integer>> &g, const Integer &n)
681{
682 integer_class _n = n.as_integer_class();
683 if (_n < 0)
684 _n = -_n;
685 if (_n <= 1)
686 return false;
687 if (_n < 5) {
688 *g = integer(_n - 1);
689 return true;
690 }
691 bool even = false;
692 if (_n % 2 == 0) {
693 if (_n % 4 == 0) {
694 return false; // If n mod 4 == 0 and n > 4, then no primitive roots.
695 }
696 _n /= 2;
697 even = true;
698 }
699 integer_class p, e;
700 if (not _prime_power(p, e, _n))
701 return false;
702 _primitive_root(_n, p, e, even);
703 *g = integer(std::move(_n));
704 return true;
705}
706
707namespace
708{
709// Computes primitive roots modulo p**e or 2*p**e where p is an odd prime.
710// References :
711// [1] Cohen H., A course in computational algebraic number theory (1996), pages
712// 25-27.
713// [2] Hackman P., Elementary number theory (2009), page 28.
715 const integer_class &p, const integer_class &e,
716 bool even = false)
717{
718 integer_class g, h, d, t, pe2, n, pm1;
719 _primitive_root(g, p, integer_class(1),
720 false); // Find one primitive root for p.
721 h = 1;
722 pm1 = p - 1;
723 // Generate other primitive roots for p. h = g**i and gcd(i, p-1) = 1.
724 // Ref[2]
725 mp_pow_ui(n, p, mp_get_ui(e));
726 for (unsigned long i = 1; i < p; ++i) {
727 h *= g;
728 h %= p;
729 mp_gcd(d, pm1, integer_class(i));
730 if (d == 1) {
731 if (e == 1) {
732 if (even and h % 2 == 0)
733 roots.push_back(integer(h + n));
734 else
735 roots.push_back(integer(h));
736 } else {
737 integer_class pp = p * p;
738 // Find d such that (h + d*p)**(p-1) mod (p**2) == 1. Ref[1]
739 // h**(p-1) - 1 = d*p*h**(p-2)
740 // d = (h - h**(2-p)) / p
741 t = 2 - p;
742 mp_powm(d, h, t, pp);
743 d = ((h - d) / p + p) % p;
744 t = h;
745 // t = h + i * p + j * p * p and i != d
746 mp_pow_ui(pe2, p, mp_get_ui(e) - 2);
747 for (unsigned long j = 0; j < pe2; ++j) {
748 for (unsigned long i = 0; i < p; ++i) {
749 if (i != d) {
750 if (even and t % 2 == 0)
751 roots.push_back(integer(t + n));
752 else
753 roots.push_back(integer(t));
754 }
755 t += p;
756 }
757 }
758 }
759 }
760 }
761} //_primitive_root_list
762} // anonymous namespace
763
765 const Integer &n)
766{
767 integer_class _n = n.as_integer_class();
768 if (_n < 0)
769 _n = -_n;
770 if (_n <= 1)
771 return;
772 if (_n < 5) {
773 roots.push_back(integer(_n - 1));
774 return;
775 }
776 bool even = false;
777 if (_n % 2 == 0) {
778 if (_n % 4 == 0) {
779 return; // If n%4 == 0 and n > 4, then no primitive roots.
780 }
781 _n /= 2;
782 even = true;
783 }
784 integer_class p, e;
785 if (not _prime_power(p, e, _n))
786 return;
789 return;
790}
791
793{
794 if (n->is_zero())
795 return integer(1);
796
797 integer_class phi = n->as_integer_class(), p;
798 if (phi < 0)
799 phi = -phi;
802
803 for (const auto &it : prime_mul) {
804 p = it.first->as_integer_class();
805 mp_divexact(phi, phi, p);
806 // phi is exactly divisible by p.
807 phi *= p - 1;
808 }
809 return integer(std::move(phi));
810}
811
813{
814 if (n->is_zero())
815 return integer(1);
816
818 integer_class lambda, t, p;
819 unsigned multiplicity;
820
822 lambda = 1;
823 for (const auto &it : prime_mul) {
824 p = it.first->as_integer_class();
825 multiplicity = it.second;
826 if (p == 2
828 > 2) { // For powers of 2 greater than 4 divide by 2.
829 multiplicity--;
830 }
831 t = p - 1;
832 mp_lcm(lambda, lambda, t);
833 mp_pow_ui(t, p, multiplicity - 1);
834 // lambda and p are relatively prime.
835 lambda = lambda * t;
836 }
837 return integer(std::move(lambda));
838}
839
840// References : Cohen H., A course in computational algebraic number theory
841// (1996), page 25.
843 const RCP<const Integer> &a,
844 const RCP<const Integer> &n)
845{
846 integer_class order, p, t;
847 integer_class _a = a->as_integer_class(),
848 _n = mp_abs(n->as_integer_class());
849 mp_gcd(t, _a, _n);
850 if (t != 1)
851 return false;
852
853 RCP<const Integer> lambda = carmichael(n);
856 _a %= _n;
857 order = lambda->as_integer_class();
858
859 for (const auto &it : prime_mul) {
860 p = it.first->as_integer_class();
861 mp_pow_ui(t, p, it.second);
862 mp_divexact(order, order, t);
863 mp_powm(t, _a, order, _n);
864 while (t != 1) {
865 mp_powm(t, t, p, _n);
866 order *= p;
867 }
868 }
869 *o = integer(std::move(order));
870 return true;
871}
872int legendre(const Integer &a, const Integer &n)
873{
874 return mp_legendre(a.as_integer_class(), n.as_integer_class());
875}
876
877int jacobi(const Integer &a, const Integer &n)
878{
879 return mp_jacobi(a.as_integer_class(), n.as_integer_class());
880}
881
882int kronecker(const Integer &a, const Integer &n)
883{
884 return mp_kronecker(a.as_integer_class(), n.as_integer_class());
885}
886
887namespace
888{
889bool _sqrt_mod_tonelli_shanks(integer_class &rop, const integer_class &a,
890 const integer_class &p)
891{
892 mp_randstate state;
893 integer_class n, y, b, q, pm1, t(1);
894 pm1 = p - 1;
895 unsigned e, m;
896 e = numeric_cast<unsigned>(mp_scan1(pm1));
897 q = pm1 >> e; // p - 1 = 2**e*q
898
899 while (t != -1) {
900 state.urandomint(n, p);
901 t = mp_legendre(n, p);
902 }
903 mp_powm(y, n, q, p); // y = n**q mod p
904 mp_powm(b, a, q, p); // b = a**q mod p
905 t = (q + 1) / 2;
906 mp_powm(rop, a, t, p); // rop = a**((q + 1) / 2) mod p
907
908 while (b != 1) {
909 m = 0;
910 t = b;
911 while (t != 1) {
912 mp_powm(t, t, integer_class(2), p);
913 ++m; // t = t**2 = b**2**(m)
914 }
915 if (m == e)
916 return false;
917 mp_pow_ui(q, integer_class(2), e - m - 1); // q = 2**(e - m - 1)
918 mp_powm(t, y, q, p); // t = y**(2**(e - m - 1))
919 mp_powm(y, t, integer_class(2), p); // y = t**2
920 e = m;
921 rop = (rop * t) % p;
922 b = (b * y) % p;
923 }
924 return true;
925}
926
927bool _sqrt_mod_prime(integer_class &rop, const integer_class &a,
928 const integer_class &p)
929{
930 if (p == 2) {
931 rop = a % p;
932 return true;
933 }
934 int l = mp_legendre(a, p);
935 integer_class t;
936 if (l == -1) {
937 return false;
938 } else if (l == 0) {
939 rop = 0;
940 } else if (p % 4 == 3) {
941 t = (p + 1) / 4;
942 mp_powm(rop, a, t, p);
943 } else if (p % 8 == 5) {
944 t = (p - 1) / 4;
945 mp_powm(t, a, t, p);
946 if (t == 1) {
947 t = (p + 3) / 8;
948 mp_powm(rop, a, t, p);
949 } else {
950 t = (p - 5) / 8;
951 integer_class t1 = 4 * a;
952 mp_powm(t, t1, t, p);
953 rop = (2 * a * t) % p;
954 }
955 } else {
956 if (p < 10000) { // If p < 10000, brute force is faster.
957 integer_class sq = integer_class(1), _a;
958 mp_fdiv_r(_a, a, p);
959 for (unsigned i = 1; i < p; ++i) {
960 if (sq == _a) {
961 rop = i;
962 return true;
963 }
964 sq += 2 * i + 1;
965 mp_fdiv_r(sq, sq, p);
966 }
967 return false;
968 } else {
969 return _sqrt_mod_tonelli_shanks(rop, a, p);
970 }
971 }
972 return true;
973}
974
975// References : Menezes, Alfred J., Paul C. Van Oorschot, and Scott A. Vanstone.
976// Handbook of applied cryptography. CRC press, 2010. pages 104 - 108
977// Calculates log = x mod q**k where g**x == a mod p and order(g, p) = n.
978void _discrete_log(integer_class &log, const integer_class &a,
979 const integer_class &g, const integer_class &n,
980 const integer_class &q, const unsigned &k,
981 const integer_class &p)
982{
983 log = 0;
984 integer_class gamma = a, alpha, _n, t, beta, qj(1), m, l;
985 _n = n / q;
986 mp_powm(alpha, g, _n, p);
987 mp_sqrtrem(m, t, q);
988 if (t != 0)
989 ++m; // m = ceiling(sqrt(q)).
990 map_integer_uint
991 table; // Table for lookup in baby-step giant-step algorithm
992 integer_class alpha_j(1), d, s;
993 s = -m;
994 mp_powm(s, alpha, s, p);
995
996 for (unsigned j = 0; j < m; ++j) {
998 alpha_j = (alpha_j * alpha) % p;
999 }
1000
1001 for (unsigned long j = 0; j < k; ++j) { // Pohlig-Hellman
1002 mp_powm(beta, gamma, _n, p);
1003 // Baby-step giant-step algorithm for l = log_alpha(beta)
1004 d = beta;
1005 bool found = false;
1006 for (unsigned i = 0; not found && i < m; ++i) {
1007 if (table.find(integer(d)) != table.end()) {
1008 l = i * m + table[integer(d)];
1009 found = true;
1010 break;
1011 }
1012 d = (d * s) % p;
1013 }
1014 _n /= q;
1015 t = -l * qj;
1016
1017 log -= t;
1018 mp_powm(t, g, t, p);
1019 gamma *= t; // gamma *= g ** (-l * (q ** j))
1020 qj *= q;
1021 }
1022}
1023
1024// References : Johnston A., A generalised qth root algorithm.
1025// Solution for x**n == a mod p**k where a != 0 mod p and p is an odd prime.
1027 const integer_class &a, const integer_class &n,
1028 const integer_class &p, const unsigned k,
1029 bool all_roots = false)
1030{
1031 integer_class _n, r, root, s, t, g(0), pk, m, phi;
1032 mp_pow_ui(pk, p, k);
1033 phi = pk * (p - 1) / p;
1034 mp_gcd(m, phi, n);
1035 t = phi / m;
1036 mp_powm(t, a, t, pk);
1037 // Check whether a**(phi / gcd(phi, n)) == 1 mod p**k.
1038 if (t != 1) {
1039 return false;
1040 }
1041 // Solve x**n == a mod p first.
1042 t = p - 1;
1043 mp_gcdext(_n, r, s, n, t);
1044 if (r < 0) {
1045 mp_fdiv_r(r, r, t / _n);
1046 }
1047 mp_powm(s, a, r, p);
1048
1049 // Solve x**(_n) == s mod p where _n | p - 1.
1050 if (_n == 1) {
1051 root = s;
1052 } else if (_n == 2) {
1053 _sqrt_mod_prime(root, s, p);
1054 } else { // Ref[1]
1055 map_integer_uint prime_mul;
1057 integer_class h, q, qt, z, v, x, s1 = s;
1058 _primitive_root(g, p, integer_class(2));
1059 unsigned c;
1060 for (const auto &it : prime_mul) {
1061 q = it.first->as_integer_class();
1062 mp_pow_ui(qt, q, it.second);
1063 h = (p - 1) / q;
1064 c = 1;
1065 while (h % q == 0) {
1066 ++c;
1067 h /= q;
1068 }
1069 mp_invert(t, h, qt);
1070 z = t * -h;
1071 x = (1 + z) / qt;
1072 mp_powm(v, s1, x, p);
1073
1074 if (c == it.second) {
1075 s1 = v;
1076 } else {
1077 mp_powm(x, s1, h, p);
1078 t = h * qt;
1079 mp_powm(r, g, t, p);
1080 mp_pow_ui(qt, q, c - it.second);
1081 _discrete_log(t, x, r, qt, q, c - it.second, p);
1082 t = -z * t;
1083 mp_powm(r, g, t, p);
1084 v *= r;
1085 mp_fdiv_r(v, v, p);
1086 s1 = v;
1087 }
1088 }
1089 root = s1;
1090 }
1091 r = n;
1092 unsigned c = 0;
1093 while (r % p == 0) {
1094 mp_divexact(r, r, p);
1095 ++c;
1096 }
1097
1098 // Solve s == x**r mod p**k where (x**r)**(p**c)) == a mod p**k
1099 integer_class pc = n / r, pd = pc * p;
1100 if (c >= 1) {
1101 mp_powm(s, root, r, p);
1102 // s == root**r mod p. Since s**(p**c) == 1 == a mod p**(c + 1), lift
1103 // until p**k.
1104 for (unsigned d = c + 2; d <= k; ++d) {
1105 t = 1 - pc;
1106 pd *= p;
1107 mp_powm(t, s, t, pd);
1108 t = (a * t - s) / pc;
1109 s += t;
1110 }
1111 } else {
1112 s = a;
1113 }
1114
1115 // Solve x**r == s mod p**k given that root**r == s mod p and r % p != 0.
1116 integer_class u;
1117 pd = p;
1118 for (unsigned d = 2; d < 2 * k; d *= 2) { // Hensel lifting
1119 t = r - 1;
1120 pd *= pd;
1121 if (d > k)
1122 pd = pk;
1123 mp_powm(u, root, t, pd);
1124 t = r * u;
1125 mp_invert(t, t, pd);
1126 root += (s - u * root) * t;
1127 mp_fdiv_r(root, root, pd);
1128 }
1129 if (m != 1 and all_roots) {
1130 // All roots are generated by root*(g**(phi / gcd(phi , n)))**j
1131 if (n == 2) {
1132 t = -1;
1133 } else {
1134 if (g == 0)
1135 _primitive_root(g, p, integer_class(2));
1136 t = phi / m;
1137 mp_powm(t, g, t, pk);
1138 }
1139 for (unsigned j = 0; j < m; ++j) {
1140 roots.push_back(integer(root));
1141 root *= t;
1142 mp_fdiv_r(root, root, pk);
1143 }
1144 } else {
1145 roots.push_back(integer(root));
1146 }
1147 return true;
1148}
1149
1150// Checks if Solution for x**n == a mod p**k exists where a != 0 mod p and p is
1151// an odd prime.
1152bool _is_nthroot_mod1(const integer_class &a, const integer_class &n,
1153 const integer_class &p, const unsigned k)
1154{
1155 integer_class t, pk, m, phi;
1156 mp_pow_ui(pk, p, k);
1157 phi = pk * (p - 1) / p;
1158 mp_gcd(m, phi, n);
1159 t = phi / m;
1160 mp_powm(t, a, t, pk);
1161 // Check whether a**(phi / gcd(phi, n)) == 1 mod p**k.
1162 if (t != 1) {
1163 return false;
1164 }
1165 return true;
1166}
1167
1168// Solution for x**n == a mod p**k.
1170 const integer_class &a, const integer_class &n,
1171 const integer_class &p, const unsigned k,
1172 bool all_roots = false)
1173{
1174 integer_class pk, root;
1176 if (a % p != 0) {
1177 if (p == 2) {
1178 integer_class r = n, t, s, pc, pj;
1179 pk = integer_class(1) << k;
1180 unsigned c = numeric_cast<unsigned>(mp_scan1(n));
1181 r = n >> c; // n = 2**c * r where r is odd.
1182
1183 // Handle special cases of k = 1 and k = 2.
1184 if (k == 1) {
1185 roots.push_back(integer(1));
1186 return true;
1187 }
1188 if (k == 2) {
1189 if (c > 0 and a % 4 == 3) {
1190 return false;
1191 }
1192 roots.push_back(integer(a % 4));
1193 if (all_roots and c > 0)
1194 roots.push_back(integer(3));
1195 return true;
1196 }
1197 if (c >= k - 2) {
1198 c = k - 2; // Since x**(2**c) == x**(2**(k - 2)) mod 2**k, let c
1199 // = k - 2.
1200 }
1201 t = integer_class(1) << (k - 2);
1202 pc = integer_class(1) << c;
1203
1204 mp_invert(s, r, t);
1205 if (c == 0) {
1206 // x**r == a mod 2**k and x**2**(k - 2) == 1 mod 2**k, implies
1207 // x**(r * s) == x == a**s mod 2**k.
1208 mp_powm(root, a, s, pk);
1209 roots.push_back(integer(root));
1210 return true;
1211 }
1212
1213 // First, solve for y**2**c == a mod 2**k where y == x**r
1214 t = integer_class(1) << (c + 2);
1215 mp_fdiv_r(t, a, t);
1216 // Check for a == y**2**c == 1 mod 2**(c + 2).
1217 if (t != 1)
1218 return false;
1219 root = 1;
1220 pj = pc * 4;
1221 // 1 is a root of x**2**c == 1 mod 2**(c + 2). Lift till 2**k.
1222 for (unsigned j = c + 2; j < k; ++j) {
1223 pj *= 2;
1224 mp_powm(t, root, pc, pj);
1225 t -= a;
1226 if (t % pj != 0)
1227 // Add 2**(j - c).
1228 root += integer_class(1) << (j - c);
1229 }
1230 // Solve x**r == root mod 2**k.
1231 mp_powm(root, root, s, pk);
1232
1233 if (all_roots) {
1234 // All roots are generated by, root * (j * (2**(k - c) +/- 1)).
1235 t = pk / pc * root;
1236 for (unsigned i = 0; i < 2; ++i) {
1237 for (unsigned long j = 0; j < pc; ++j) {
1238 roots.push_back(integer(root));
1239 root += t;
1240 }
1241 root = t - root;
1242 }
1243 } else {
1244 roots.push_back(integer(root));
1245 }
1246 return true;
1247 } else {
1248 return _nthroot_mod1(roots, a, n, p, k, all_roots);
1249 }
1250 } else {
1251 integer_class _a;
1252 mp_pow_ui(pk, p, k);
1253 _a = a % pk;
1254 unsigned m;
1255 integer_class pm;
1256 if (_a == 0) {
1257 if (not all_roots) {
1258 roots.push_back(integer(0));
1259 return true;
1260 }
1261 _roots.push_back(integer(0));
1262 if (n >= k)
1263 m = k - 1;
1264 else
1265 m = numeric_cast<unsigned>(k - 1 - (k - 1) / mp_get_ui(n));
1266 mp_pow_ui(pm, p, m);
1267 } else {
1268 unsigned r = 1;
1269 mp_divexact(_a, _a, p);
1270 while (_a % p == 0) {
1271 mp_divexact(_a, _a, p);
1272 ++r;
1273 }
1274 if (r < n or r % n != 0
1276 all_roots)) {
1277 return false;
1278 }
1279 m = numeric_cast<unsigned>(r / mp_get_ui(n));
1280 mp_pow_ui(pm, p, m);
1281 if (not all_roots) {
1282 roots.push_back(
1283 integer(_roots.back()->as_integer_class() * pm));
1284 return true;
1285 }
1286 for (auto &it : _roots) {
1287 it = integer(it->as_integer_class() * pm);
1288 }
1289 m = numeric_cast<unsigned>(r - r / mp_get_ui(n));
1290 mp_pow_ui(pm, p, m);
1291 }
1292 integer_class pkm;
1293 mp_pow_ui(pkm, p, k - m);
1294
1295 for (const auto &it : _roots) {
1296 root = it->as_integer_class();
1297 for (unsigned long i = 0; i < pm; ++i) {
1298 roots.push_back(integer(root));
1299 root += pkm;
1300 }
1301 }
1302 }
1303 return true;
1304}
1305} // anonymous namespace
1306
1307// Returns whether Solution for x**n == a mod p**k exists or not
1308bool _is_nthroot_mod_prime_power(const integer_class &a, const integer_class &n,
1309 const integer_class &p, const unsigned k)
1310{
1311 integer_class pk;
1312 if (a % p != 0) {
1313 if (p == 2) {
1314 integer_class t;
1315 unsigned c = numeric_cast<unsigned>(mp_scan1(n));
1316
1317 // Handle special cases of k = 1 and k = 2.
1318 if (k == 1) {
1319 return true;
1320 }
1321 if (k == 2) {
1322 if (c > 0 and a % 4 == 3) {
1323 return false;
1324 }
1325 return true;
1326 }
1327 if (c >= k - 2) {
1328 c = k - 2; // Since x**(2**c) == x**(2**(k - 2)) mod 2**k, let c
1329 // = k - 2.
1330 }
1331 if (c == 0) {
1332 // x**r == a mod 2**k and x**2**(k - 2) == 1 mod 2**k, implies
1333 // x**(r * s) == x == a**s mod 2**k.
1334 return true;
1335 }
1336
1337 // First, solve for y**2**c == a mod 2**k where y == x**r
1338 t = integer_class(1) << (c + 2);
1339 mp_fdiv_r(t, a, t);
1340 // Check for a == y**2**c == 1 mod 2**(c + 2).
1341 if (t != 1)
1342 return false;
1343 return true;
1344 } else {
1345 return _is_nthroot_mod1(a, n, p, k);
1346 }
1347 } else {
1348 integer_class _a;
1349 mp_pow_ui(pk, p, k);
1350 _a = a % pk;
1351 integer_class pm;
1352 if (_a == 0) {
1353 return true;
1354 } else {
1355 unsigned r = 1;
1356 mp_divexact(_a, _a, p);
1357 while (_a % p == 0) {
1358 mp_divexact(_a, _a, p);
1359 ++r;
1360 }
1361 if (r < n or r % n != 0
1362 or not _is_nthroot_mod_prime_power(_a, n, p, k - r)) {
1363 return false;
1364 }
1365 return true;
1366 }
1367 }
1368 return true;
1369}
1370
1371bool nthroot_mod(const Ptr<RCP<const Integer>> &root,
1372 const RCP<const Integer> &a, const RCP<const Integer> &n,
1373 const RCP<const Integer> &mod)
1374{
1375 if (mod->as_integer_class() <= 0) {
1376 return false;
1377 } else if (mod->as_integer_class() == 1) {
1378 *root = integer(0);
1379 return true;
1380 }
1384 bool ret_val;
1385
1387 for (const auto &it : prime_mul) {
1388 integer_class _mod;
1389 mp_pow_ui(_mod, it.first->as_integer_class(), it.second);
1390 moduli.push_back(integer(std::move(_mod)));
1392 rem, a->as_integer_class(), n->as_integer_class(),
1393 it.first->as_integer_class(), it.second, false);
1394 if (not ret_val)
1395 return false;
1396 }
1397 crt(root, rem, moduli);
1398 return true;
1399}
1400
1402 const RCP<const Integer> &a, const RCP<const Integer> &n,
1403 const RCP<const Integer> &m)
1404{
1405 if (m->as_integer_class() <= 0) {
1406 return;
1407 } else if (m->as_integer_class() == 1) {
1408 roots.push_back(integer(0));
1409 return;
1410 }
1414 bool ret_val;
1415
1417 for (const auto &it : prime_mul) {
1418 integer_class _mod;
1419 mp_pow_ui(_mod, it.first->as_integer_class(), it.second);
1420 moduli.push_back(integer(std::move(_mod)));
1423 rem1, a->as_integer_class(), n->as_integer_class(),
1424 it.first->as_integer_class(), it.second, true);
1425 if (not ret_val)
1426 return;
1427 rem.push_back(rem1);
1428 }
1431}
1432
1434 const RCP<const Number> &b, const RCP<const Integer> &m)
1435{
1436 if (is_a<Integer>(*b)) {
1437 integer_class t = down_cast<const Integer &>(*b).as_integer_class();
1438 if (b->is_negative())
1439 t *= -1;
1440 mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1441 if (b->is_negative()) {
1442 bool ret_val = mp_invert(t, t, m->as_integer_class());
1443 if (not ret_val)
1444 return false;
1445 }
1446 *powm = integer(std::move(t));
1447 return true;
1448 } else if (is_a<Rational>(*b)) {
1450 get_num_den(down_cast<const Rational &>(*b), outArg(num), outArg(den));
1451 if (den->is_negative()) {
1452 den = den->mulint(*minus_one);
1453 num = num->mulint(*minus_one);
1454 }
1455 integer_class t = mp_abs(num->as_integer_class());
1456 mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1457 if (num->is_negative()) {
1458 bool ret_val = mp_invert(t, t, m->as_integer_class());
1459 if (not ret_val)
1460 return false;
1461 }
1462 r = integer(std::move(t));
1463 return nthroot_mod(powm, r, den, m);
1464 }
1465 return false;
1466}
1467
1469 const RCP<const Integer> &a, const RCP<const Number> &b,
1470 const RCP<const Integer> &m)
1471{
1472 if (is_a<Integer>(*b)) {
1473 integer_class t
1474 = mp_abs(down_cast<const Integer &>(*b).as_integer_class());
1475 mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1476 if (b->is_negative()) {
1477 bool ret_val = mp_invert(t, t, m->as_integer_class());
1478 if (not ret_val)
1479 return;
1480 }
1481 pows.push_back(integer(std::move(t)));
1482 } else if (is_a<Rational>(*b)) {
1484 get_num_den(down_cast<const Rational &>(*b), outArg(num), outArg(den));
1485 if (den->is_negative()) {
1486 den = den->mulint(*integer(-1));
1487 num = num->mulint(*integer(-1));
1488 }
1489 integer_class t = num->as_integer_class();
1490 if (num->is_negative())
1491 t *= -1;
1492 mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1493 if (num->is_negative()) {
1494 bool ret_val = mp_invert(t, t, m->as_integer_class());
1495 if (not ret_val)
1496 return;
1497 }
1498 r = integer(t);
1499 nthroot_mod_list(pows, r, den, m);
1500 }
1501}
1502
1504{
1505 /*
1506 Returns the list of quadratic residues.
1507 Example
1508 ========
1509 >>> quadratic_residues(7)
1510 [0, 1, 2, 4]
1511 */
1512
1513 if (a.as_integer_class() < 1) {
1514 throw SymEngineException("quadratic_residues: Input must be > 0");
1515 }
1516
1518 for (integer_class i = integer_class(0); i <= a.as_int() / 2; i++) {
1519 residue.push_back((i * i) % a.as_int());
1520 }
1521
1522 sort(residue.begin(), residue.end());
1523 residue.erase(unique(residue.begin(), residue.end()), residue.end());
1524
1525 return residue;
1526}
1527
1528bool is_quad_residue(const Integer &a, const Integer &p)
1529{
1530 /*
1531 Returns true if ``a`` (mod ``p``) is in the set of squares mod ``p``,
1532 i.e a % p in set([i**2 % p for i in range(p)]). If ``p`` is an odd but
1533 not prime, an iterative method is used to make the determination.
1534 */
1535
1536 integer_class p2 = p.as_integer_class();
1537 if (p2 == 0)
1538 throw SymEngineException(
1539 "is_quad_residue: Second parameter must be non-zero");
1540 if (p2 < 0)
1541 p2 = -p2;
1542 integer_class a_final = a.as_integer_class();
1543 if (a.as_integer_class() >= p2 || a.as_integer_class() < 0)
1544 mp_fdiv_r(a_final, a.as_integer_class(), p2);
1545 if (a_final < 2)
1546 return true;
1547
1548 if (!probab_prime_p(*integer(p2))) {
1549 if ((p2 % 2 == 1) && jacobi(*integer(a_final), p) == -1)
1550 return false;
1551
1554
1557 bool ret_val;
1558
1559 for (const auto &it : prime_mul) {
1560 ret_val = _is_nthroot_mod_prime_power(
1561 a1->as_integer_class(), integer(2)->as_integer_class(),
1562 it.first->as_integer_class(), it.second);
1563 if (not ret_val)
1564 return false;
1565 }
1566 return true;
1567 }
1568
1569 return mp_legendre(a_final, p2) == 1;
1570}
1571
1572bool is_nth_residue(const Integer &a, const Integer &n, const Integer &mod)
1573/*
1574Returns true if ``a`` (mod ``mod``) is in the set of nth powers mod ``mod``,
1575i.e a % mod in set([i**n % mod for i in range(mod)]).
1576*/
1577{
1578 integer_class _mod = mod.as_integer_class();
1579
1580 if (_mod == 0) {
1581 return false;
1582 } else if (_mod == 1) {
1583 return true;
1584 }
1585
1586 if (_mod < 0)
1587 _mod = -(_mod);
1588
1592 bool ret_val;
1593
1594 for (const auto &it : prime_mul) {
1595 ret_val = _is_nthroot_mod_prime_power(
1597 it.first->as_integer_class(), it.second);
1598 if (not ret_val)
1599 return false;
1600 }
1601 return true;
1602}
1603
1604int mobius(const Integer &a)
1605{
1606 if (a.as_int() <= 0) {
1607 throw SymEngineException("mobius: Integer <= 0");
1608 }
1610 bool is_square_free = true;
1612 auto num_prime_factors = prime_mul.size();
1613 for (const auto &it : prime_mul) {
1614 int p_freq = it.second;
1615 if (p_freq > 1) {
1616 is_square_free = false;
1617 break;
1618 }
1619 }
1620 if (!is_square_free) {
1621 return 0;
1622 } else if (num_prime_factors % 2 == 0) {
1623 return 1;
1624 } else {
1625 return -1;
1626 }
1627}
1628
1629long mertens(const unsigned long a)
1630{
1631 long mertens = 0;
1632 for (unsigned long i = 1; i <= a; ++i) {
1633 mertens += mobius(*(integer(i)));
1634 }
1635 return mertens;
1636}
1637
1648integer_class mp_polygonal_number(const integer_class &s,
1649 const integer_class &n)
1650{
1651 auto res = ((s - 2) * n * n - (s - 4) * n) / 2;
1652 return res;
1653}
1654
1665integer_class mp_principal_polygonal_root(const integer_class &s,
1666 const integer_class &x)
1667{
1668 integer_class tmp;
1669 mp_pow_ui(tmp, s - 4, 2);
1670 integer_class root = mp_sqrt(8 * x * (s - 2) + tmp);
1671 integer_class n = (root + s - 4) / (2 * (s - 2));
1672 return n;
1673}
1674
1677{
1678 // From
1679 // https://codegolf.stackexchange.com/questions/1935/fastest-algorithm-for-decomposing-a-perfect-power
1680 unsigned long p = 2;
1681 integer_class intone, i, j, m, res;
1682 intone = 1;
1685
1686 while ((intone << p) <= n) {
1687 i = 2;
1688 j = n;
1689 while (j > i + 1) {
1690 m = (i + j) / 2;
1691 mp_pow_ui(res, m, p);
1692 if (res > n)
1693 j = m;
1694 else
1695 i = m;
1696 }
1697 mp_pow_ui(res, i, p);
1698 if (res == n) {
1699 respair = std::make_pair(i, p);
1700 if (lowest_exponent) {
1701 return respair;
1702 }
1703 }
1704 p++;
1705 }
1706 return respair;
1707}
1708
1709} // namespace SymEngine
Classes and functions relating to the binary operation of addition.
Integer Class.
Definition integer.h:19
const integer_class & as_integer_class() const
Convert to integer_class.
Definition integer.h:45
signed long int as_int() const
Convert to int, raise an exception if it does not fit.
Definition integer.cpp:34
static RCP< const Number > from_mpq(const rational_class &i)
Definition rational.cpp:23
T make_pair(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition add.cpp:19
integer_class mp_polygonal_number(const integer_class &s, const integer_class &n)
Numeric calculation of the n:th s-gonal number.
Definition ntheory.cpp:1648
bool divides(const Integer &a, const Integer &b)
Definition ntheory.cpp:161
RCP< const Integer > nextprime(const Integer &a)
Definition ntheory.cpp:172
int kronecker(const Integer &a, const Integer &n)
Kronecker Function.
Definition ntheory.cpp:882
void quotient_mod_f(const Ptr< RCP< const Integer > > &q, const Ptr< RCP< const Integer > > &r, const Integer &n, const Integer &d)
Definition ntheory.cpp:100
void prime_factor_multiplicities(map_integer_uint &primes_mul, const Integer &n)
Find multiplicities of prime factors of n
Definition ntheory.cpp:461
RCP< const Basic > beta(const RCP< const Basic > &x, const RCP< const Basic > &y)
Canonicalize Beta:
RCP< const Integer > gcd(const Integer &a, const Integer &b)
Greatest Common Divisor.
Definition ntheory.cpp:31
RCP< const Integer > mod(const Integer &n, const Integer &d)
modulo round toward zero
Definition ntheory.cpp:66
int mobius(const Integer &a)
Mobius Function.
Definition ntheory.cpp:1604
void lucas2(const Ptr< RCP< const Integer > > &g, const Ptr< RCP< const Integer > > &s, unsigned long n)
Lucas number n and n-1.
Definition ntheory.cpp:134
bool multiplicative_order(const Ptr< RCP< const Integer > > &o, const RCP< const Integer > &a, const RCP< const Integer > &n)
Multiplicative order. Return false if order does not exist.
Definition ntheory.cpp:842
RCP< const Integer > lucas(unsigned long n)
Lucas number.
Definition ntheory.cpp:127
bool nthroot_mod(const Ptr< RCP< const Integer > > &root, const RCP< const Integer > &a, const RCP< const Integer > &n, const RCP< const Integer > &mod)
A solution to x**n == a mod m. Return false if none exists.
Definition ntheory.cpp:1371
RCP< const Integer > quotient_f(const Integer &n, const Integer &d)
Definition ntheory.cpp:93
void prime_factors(std::vector< RCP< const Integer > > &prime_list, const Integer &n)
Find prime factors of n
Definition ntheory.cpp:431
void hash_combine(hash_t &seed, const T &v)
Definition basic-inl.h:95
RCP< const Integer > fibonacci(unsigned long n)
Fibonacci number.
Definition ntheory.cpp:110
RCP< const Basic > gamma(const RCP< const Basic > &arg)
Canonicalize Gamma:
bool primitive_root(const Ptr< RCP< const Integer > > &g, const Integer &n)
Computes a primitive root. Returns false if no primitive root exists.
Definition ntheory.cpp:680
bool is_quad_residue(const Integer &a, const Integer &p)
Returns true if 'a' is a quadratic residue of 'p'.
Definition ntheory.cpp:1528
void gcd_ext(const Ptr< RCP< const Integer > > &g, const Ptr< RCP< const Integer > > &s, const Ptr< RCP< const Integer > > &t, const Integer &a, const Integer &b)
Extended GCD.
Definition ntheory.cpp:38
RCP< const Integer > mod_f(const Integer &n, const Integer &d)
modulo round toward -inf
Definition ntheory.cpp:86
RCP< const Integer > quotient(const Integer &n, const Integer &d)
Definition ntheory.cpp:71
vec_integer_class quadratic_residues(const Integer &a)
Finds all Quadratic Residues of a Positive Integer.
Definition ntheory.cpp:1503
void get_num_den(const Rational &rat, const Ptr< RCP< const Integer > > &num, const Ptr< RCP< const Integer > > &den)
returns the num and den of rational rat as RCP<const Integer>
Definition rational.cpp:130
int jacobi(const Integer &a, const Integer &n)
Jacobi Function.
Definition ntheory.cpp:877
bool powermod(const Ptr< RCP< const Integer > > &powm, const RCP< const Integer > &a, const RCP< const Number > &b, const RCP< const Integer > &m)
Definition ntheory.cpp:1433
RCP< const Integer > lcm(const Integer &a, const Integer &b)
Least Common Multiple.
Definition ntheory.cpp:49
RCP< const Integer > binomial(const Integer &n, unsigned long k)
Binomial Coefficient.
Definition ntheory.cpp:145
void nthroot_mod_list(std::vector< RCP< const Integer > > &roots, const RCP< const Integer > &a, const RCP< const Integer > &n, const RCP< const Integer > &m)
All Solutions to x**n == a mod m. Return false if none exists.
Definition ntheory.cpp:1401
bool crt(const Ptr< RCP< const Integer > > &R, const std::vector< RCP< const Integer > > &rem, const std::vector< RCP< const Integer > > &mod)
Chinese remainder function. Return true when a solution exists.
Definition ntheory.cpp:552
int factor_lehman_method(const Ptr< RCP< const Integer > > &f, const Integer &n)
Factor using lehman's methods.
Definition ntheory.cpp:253
int factor_pollard_pm1_method(const Ptr< RCP< const Integer > > &f, const Integer &n, unsigned B, unsigned retries)
Factor using Pollard's p-1 method.
Definition ntheory.cpp:296
void insert(T1 &m, const T2 &first, const T3 &second)
Definition dict.h:83
std::pair< integer_class, integer_class > mp_perfect_power_decomposition(const integer_class &n, bool lowest_exponent)
Decompose a positive integer into perfect powers.
Definition ntheory.cpp:1676
int factor(const Ptr< RCP< const Integer > > &f, const Integer &n, double B1)
Definition ntheory.cpp:370
int factor_trial_division(const Ptr< RCP< const Integer > > &f, const Integer &n)
Definition ntheory.cpp:421
void fibonacci2(const Ptr< RCP< const Integer > > &g, const Ptr< RCP< const Integer > > &s, unsigned long n)
Fibonacci n and n-1.
Definition ntheory.cpp:117
int probab_prime_p(const Integer &a, unsigned reps)
Probabilistic Prime.
Definition ntheory.cpp:167
RCP< const Basic > log(const RCP< const Basic > &arg)
Returns the Natural Logarithm from argument arg
RCP< const Integer > totient(const RCP< const Integer > &n)
Euler's totient function.
Definition ntheory.cpp:792
bool is_nth_residue(const Integer &a, const Integer &n, const Integer &mod)
Returns true if 'a' is a nth power residue of 'mod'.
Definition ntheory.cpp:1572
int factor_pollard_rho_method(const Ptr< RCP< const Integer > > &f, const Integer &n, unsigned retries)
Factor using Pollard's rho methods.
Definition ntheory.cpp:348
RCP< const Integer > factorial(unsigned long n)
Factorial.
Definition ntheory.cpp:153
RCP< const Integer > carmichael(const RCP< const Integer > &n)
Carmichael function.
Definition ntheory.cpp:812
integer_class mp_principal_polygonal_root(const integer_class &s, const integer_class &x)
Numeric calculation of the principal s-gonal root of x.
Definition ntheory.cpp:1665
void quotient_mod(const Ptr< RCP< const Integer > > &q, const Ptr< RCP< const Integer > > &r, const Integer &n, const Integer &d)
Definition ntheory.cpp:76
void powermod_list(std::vector< RCP< const Integer > > &pows, const RCP< const Integer > &a, const RCP< const Number > &b, const RCP< const Integer > &m)
Definition ntheory.cpp:1468
std::enable_if< std::is_integral< T >::value, RCP< constInteger > >::type integer(T i)
Definition integer.h:197
RCP< const Number > bernoulli(unsigned long n)
Definition ntheory.cpp:495
int legendre(const Integer &a, const Integer &n)
Legendre Function.
Definition ntheory.cpp:872
RCP< const Number > harmonic(unsigned long n, long m)
Computes the sum of the inverses of the first perfect mth powers.
Definition ntheory.cpp:522
int mod_inverse(const Ptr< RCP< const Integer > > &b, const Integer &a, const Integer &m)
inverse modulo
Definition ntheory.cpp:56
void primitive_root_list(std::vector< RCP< const Integer > > &roots, const Integer &n)
Definition ntheory.cpp:764
T push_back(T... args)
T sort(T... args)
less operator (<) for Integers
Definition integer.h:185