4#include <symengine/prime_sieve.h>
10#ifdef HAVE_SYMENGINE_ECM
13#ifdef HAVE_SYMENGINE_PRIMESIEVE
14#include <primesieve.hpp>
16#ifdef 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>
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,
42 integer_class g_, s_, t_;
77 const Ptr<RCP<const Integer>> &r,
const Integer &n,
101 const Ptr<RCP<const Integer>> &r,
const Integer &n,
104 integer_class _q, _r;
118 const Ptr<RCP<const Integer>> &s,
unsigned long n)
122 mp_fib2_ui(g_t, s_t, n);
127RCP<const Integer>
lucas(
unsigned long n)
134void lucas2(
const Ptr<RCP<const Integer>> &g,
const Ptr<RCP<const Integer>> &s,
139 mp_lucnum2_ui(g_t, s_t, n);
182int _factor_trial_division_sieve(integer_class &
factor,
const integer_class &N)
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));
190 while ((p = pi.next_prime()) <= limit) {
199int _factor_lehman_method(integer_class &rop,
const integer_class &n)
202 throw SymEngineException(
"Require n >= 21 to use lehman method");
205 integer_class u_bound;
207 mp_root(u_bound, n, 3);
208 u_bound = u_bound + 1;
210 Sieve::iterator pi(numeric_cast<unsigned>(mp_get_ui(u_bound)));
212 while ((p = pi.next_prime()) <= mp_get_ui(u_bound)) {
222 integer_class k, a, b, l;
226 while (k <= u_bound) {
227 a = mp_sqrt(4 * k * n);
234 l = a * a - 4 * k * n;
235 if (mp_perfect_square_p(l)) {
266int _factor_pollard_pm1_method(integer_class &rop,
const integer_class &n,
267 const integer_class &c,
unsigned B)
270 throw SymEngineException(
271 "Require n > 3 and B > 2 to use Pollard's p-1 method");
276 Sieve::iterator pi(B);
278 while ((p = pi.next_prime()) <= B) {
284 mp_powm(_c, _c, m, n);
289 if (rop == 1 or rop == n)
297 const Integer &n,
unsigned B,
unsigned retries)
300 integer_class rop, nm4, c;
305 for (
unsigned i = 0; i < retries and ret_val == 0; ++i) {
306 state.urandomint(c, nm4);
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)
324 throw SymEngineException(
"Require n > 4 to use pollard's-rho method");
326 integer_class u, v, g, m;
330 for (
unsigned i = 0; i < steps; ++i) {
349 const Integer &n,
unsigned retries)
352 integer_class rop, nm1, nm4, a, s;
357 for (
unsigned i = 0; i < retries and ret_val == 0; ++i) {
358 state.urandomint(a, nm1);
359 state.urandomint(s, nm4);
373 integer_class _n, _f;
377#ifdef HAVE_SYMENGINE_ECM
378 if (mp_perfect_power_p(_n)) {
380 unsigned long int i = 1;
381 integer_class m, rem;
391 while (i > 1 and rem != 0) {
392 mp_rootrem(_f, rem, _n, i);
399 if (mp_probab_prime_p(_n, 25) > 0) {
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);
408 throw SymEngineException(
409 "ECM failed to factor the given number");
414 ret_val = _factor_trial_division_sieve(_f, _n);
442 auto limit = mp_get_ui(sqrtN);
443 if (not mp_fits_ulong_p(sqrtN)
445 throw SymEngineException(
"N too large to factor");
449 while ((p = pi.next_prime()) <= limit) {
450 while (_n % p == 0) {
451 prime_list.push_back(
integer(p));
472 auto limit = mp_get_ui(sqrtN);
473 if (not mp_fits_ulong_p(sqrtN)
475 throw SymEngineException(
"N too large to factor");
479 while ((p = pi.next_prime()) <= limit) {
481 while (_n % p == 0) {
497#ifdef HAVE_SYMENGINE_ARB
500 bernoulli_fmpq_ui(res, n);
503 fmpq_get_mpq(a, res);
511 for (
unsigned m = 0; m <= n; ++m) {
512 v[m] = rational_class(1u, m + 1);
514 for (
unsigned j = m; j >= 1; --j) {
515 v[j - 1] = j * (v[j - 1] - v[j]);
524 rational_class res(0);
526 for (
unsigned i = 1; i <= n; ++i) {
527 res += rational_class(1u, i);
531 for (
unsigned i = 1; i <= n; ++i) {
533 rational_class t(1u, i);
534#if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP
535 mp_pow_ui(get_den(t), get_den(t), m);
542 mp_pow_ui(t, t,
static_cast<unsigned long>(-m));
552bool crt(
const Ptr<RCP<const Integer>> &R,
556 if (
mod.size() > rem.size())
557 throw SymEngineException(
"Too few remainders");
559 throw SymEngineException(
"Moduli vector cannot be empty");
561 integer_class m, r, g, s, t;
562 m =
mod[0]->as_integer_class();
563 r = rem[0]->as_integer_class();
565 for (
unsigned i = 1; i <
mod.size(); ++i) {
566 mp_gcdext(g, s, t, m,
mod[i]->as_integer_class());
568 t = rem[i]->as_integer_class() - r;
569 if (not mp_divisible_p(t, g))
571 r += m * s * (t / g);
572 m *=
mod[i]->as_integer_class() / g;
583void _crt_cartesian(
std::vector<RCP<const Integer>> &R,
587 if (
mod.size() > rem.size())
588 throw SymEngineException(
"Too few remainders");
590 throw SymEngineException(
"Moduli vector cannot be empty");
591 integer_class m, _m, r, s, t;
592 m =
mod[0]->as_integer_class();
595 for (
unsigned i = 1; i <
mod.size(); ++i) {
597 mp_invert(s, m,
mod[i]->as_integer_class());
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);
614bool _prime_power(integer_class &p, integer_class &e,
const integer_class &n)
618 integer_class _n = n, temp;
621 while (mp_perfect_power_p(_n) and _n >= 2) {
622 if (mp_root(temp, _n, i)) {
629 if (mp_probab_prime_p(_n, 25)) {
639void _primitive_root(integer_class &g,
const integer_class &p,
640 const integer_class &e,
bool even =
false)
649 for (
const auto &it : primes) {
650 t = it->as_integer_class();
665 integer_class pm1 = p - 1;
666 mp_powm(t, g, pm1, t);
672 if (even and g % 2 == 0) {
673 mp_pow_ui(t, p, mp_get_ui(e));
700 if (not _prime_power(p, e, _n))
702 _primitive_root(_n, p, e, even);
714void _primitive_root_list(
std::vector<RCP<const Integer>> &roots,
715 const integer_class &p,
const integer_class &e,
718 integer_class g, h, d, t, pe2, n, pm1;
719 _primitive_root(g, p, integer_class(1),
725 mp_pow_ui(n, p, mp_get_ui(e));
726 for (
unsigned long i = 1; i < p; ++i) {
729 mp_gcd(d, pm1, integer_class(i));
732 if (even and h % 2 == 0)
733 roots.push_back(
integer(h + n));
737 integer_class pp = p * p;
742 mp_powm(d, h, t, pp);
743 d = ((h - d) / p + p) % p;
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) {
750 if (even and t % 2 == 0)
751 roots.push_back(
integer(t + n));
773 roots.push_back(
integer(_n - 1));
785 if (not _prime_power(p, e, _n))
787 _primitive_root_list(roots, p, e, even);
792RCP<const Integer>
totient(
const RCP<const Integer> &n)
797 integer_class phi = n->as_integer_class(), p;
803 for (
const auto &it : prime_mul) {
804 p = it.first->as_integer_class();
805 mp_divexact(phi, phi, p);
818 integer_class lambda, t, p;
819 unsigned multiplicity;
823 for (
const auto &it : prime_mul) {
824 p = it.first->as_integer_class();
825 multiplicity = it.second;
832 mp_lcm(lambda, lambda, t);
833 mp_pow_ui(t, p, multiplicity - 1);
843 const RCP<const Integer> &a,
844 const RCP<const Integer> &n)
846 integer_class order, p, t;
847 integer_class _a = a->as_integer_class(),
848 _n = mp_abs(n->as_integer_class());
857 order = lambda->as_integer_class();
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);
865 mp_powm(t, t, p, _n);
889bool _sqrt_mod_tonelli_shanks(integer_class &rop,
const integer_class &a,
890 const integer_class &p)
893 integer_class n, y, b, q, pm1, t(1);
896 e = numeric_cast<unsigned>(mp_scan1(pm1));
900 state.urandomint(n, p);
901 t = mp_legendre(n, p);
906 mp_powm(rop, a, t, p);
912 mp_powm(t, t, integer_class(2), p);
917 mp_pow_ui(q, integer_class(2), e - m - 1);
919 mp_powm(y, t, integer_class(2), p);
927bool _sqrt_mod_prime(integer_class &rop,
const integer_class &a,
928 const integer_class &p)
934 int l = mp_legendre(a, p);
940 }
else if (p % 4 == 3) {
942 mp_powm(rop, a, t, p);
943 }
else if (p % 8 == 5) {
948 mp_powm(rop, a, t, p);
951 integer_class t1 = 4 * a;
952 mp_powm(t, t1, t, p);
953 rop = (2 * a * t) % p;
957 integer_class sq = integer_class(1), _a;
959 for (
unsigned i = 1; i < p; ++i) {
965 mp_fdiv_r(sq, sq, p);
969 return _sqrt_mod_tonelli_shanks(rop, a, p);
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)
984 integer_class
gamma = a, alpha, _n, t,
beta, qj(1), m, l;
986 mp_powm(alpha, g, _n, p);
992 integer_class alpha_j(1), d, s;
994 mp_powm(s, alpha, s, p);
996 for (
unsigned j = 0; j < m; ++j) {
998 alpha_j = (alpha_j * alpha) % p;
1001 for (
unsigned long j = 0; j < k; ++j) {
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)];
1018 mp_powm(t, g, t, p);
1026bool _nthroot_mod1(
std::vector<RCP<const Integer>> &roots,
1027 const integer_class &a,
const integer_class &n,
1028 const integer_class &p,
const unsigned k,
1029 bool all_roots =
false)
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;
1036 mp_powm(t, a, t, pk);
1043 mp_gcdext(_n, r, s, n, t);
1045 mp_fdiv_r(r, r, t / _n);
1047 mp_powm(s, a, r, p);
1052 }
else if (_n == 2) {
1053 _sqrt_mod_prime(root, s, p);
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));
1060 for (
const auto &it : prime_mul) {
1061 q = it.first->as_integer_class();
1062 mp_pow_ui(qt, q, it.second);
1065 while (h % q == 0) {
1069 mp_invert(t, h, qt);
1072 mp_powm(v, s1, x, p);
1074 if (c == it.second) {
1077 mp_powm(x, s1, h, p);
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);
1083 mp_powm(r, g, t, p);
1093 while (r % p == 0) {
1094 mp_divexact(r, r, p);
1099 integer_class pc = n / r, pd = pc * p;
1101 mp_powm(s, root, r, p);
1104 for (
unsigned d = c + 2; d <= k; ++d) {
1107 mp_powm(t, s, t, pd);
1108 t = (a * t - s) / pc;
1118 for (
unsigned d = 2; d < 2 * k; d *= 2) {
1123 mp_powm(u, root, t, pd);
1125 mp_invert(t, t, pd);
1126 root += (s - u * root) * t;
1127 mp_fdiv_r(root, root, pd);
1129 if (m != 1 and all_roots) {
1135 _primitive_root(g, p, integer_class(2));
1137 mp_powm(t, g, t, pk);
1139 for (
unsigned j = 0; j < m; ++j) {
1140 roots.push_back(
integer(root));
1142 mp_fdiv_r(root, root, pk);
1145 roots.push_back(
integer(root));
1152bool _is_nthroot_mod1(
const integer_class &a,
const integer_class &n,
1153 const integer_class &p,
const unsigned k)
1155 integer_class t, pk, m, phi;
1156 mp_pow_ui(pk, p, k);
1157 phi = pk * (p - 1) / p;
1160 mp_powm(t, a, t, pk);
1169bool _nthroot_mod_prime_power(
std::vector<RCP<const Integer>> &roots,
1170 const integer_class &a,
const integer_class &n,
1171 const integer_class &p,
const unsigned k,
1172 bool all_roots =
false)
1174 integer_class pk, root;
1178 integer_class r = n, t, s, pc, pj;
1179 pk = integer_class(1) << k;
1180 unsigned c = numeric_cast<unsigned>(mp_scan1(n));
1189 if (c > 0 and a % 4 == 3) {
1192 roots.push_back(
integer(a % 4));
1193 if (all_roots and c > 0)
1201 t = integer_class(1) << (k - 2);
1202 pc = integer_class(1) << c;
1208 mp_powm(root, a, s, pk);
1209 roots.push_back(
integer(root));
1214 t = integer_class(1) << (c + 2);
1222 for (
unsigned j = c + 2; j < k; ++j) {
1224 mp_powm(t, root, pc, pj);
1228 root += integer_class(1) << (j - c);
1231 mp_powm(root, root, s, pk);
1236 for (
unsigned i = 0; i < 2; ++i) {
1237 for (
unsigned long j = 0; j < pc; ++j) {
1238 roots.push_back(
integer(root));
1244 roots.push_back(
integer(root));
1248 return _nthroot_mod1(roots, a, n, p, k, all_roots);
1252 mp_pow_ui(pk, p, k);
1257 if (not all_roots) {
1265 m = numeric_cast<unsigned>(k - 1 - (k - 1) / mp_get_ui(n));
1266 mp_pow_ui(pm, p, m);
1269 mp_divexact(_a, _a, p);
1270 while (_a % p == 0) {
1271 mp_divexact(_a, _a, p);
1274 if (r < n or r % n != 0
1275 or not _nthroot_mod_prime_power(_roots, _a, n, p, k - r,
1279 m = numeric_cast<unsigned>(r / mp_get_ui(n));
1280 mp_pow_ui(pm, p, m);
1281 if (not all_roots) {
1286 for (
auto &it : _roots) {
1287 it =
integer(it->as_integer_class() * pm);
1289 m = numeric_cast<unsigned>(r - r / mp_get_ui(n));
1290 mp_pow_ui(pm, p, m);
1293 mp_pow_ui(pkm, p, k - m);
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));
1308bool _is_nthroot_mod_prime_power(
const integer_class &a,
const integer_class &n,
1309 const integer_class &p,
const unsigned k)
1315 unsigned c = numeric_cast<unsigned>(mp_scan1(n));
1322 if (c > 0 and a % 4 == 3) {
1338 t = integer_class(1) << (c + 2);
1345 return _is_nthroot_mod1(a, n, p, k);
1349 mp_pow_ui(pk, p, k);
1356 mp_divexact(_a, _a, p);
1357 while (_a % p == 0) {
1358 mp_divexact(_a, _a, p);
1361 if (r < n or r % n != 0
1362 or not _is_nthroot_mod_prime_power(_a, n, p, k - r)) {
1372 const RCP<const Integer> &a,
const RCP<const Integer> &n,
1373 const RCP<const Integer> &
mod)
1375 if (
mod->as_integer_class() <= 0) {
1377 }
else if (
mod->as_integer_class() == 1) {
1387 for (
const auto &it : prime_mul) {
1389 mp_pow_ui(_mod, it.first->as_integer_class(), it.second);
1391 ret_val = _nthroot_mod_prime_power(
1392 rem, a->as_integer_class(), n->as_integer_class(),
1393 it.first->as_integer_class(), it.second,
false);
1397 crt(root, rem, moduli);
1402 const RCP<const Integer> &a,
const RCP<const Integer> &n,
1403 const RCP<const Integer> &m)
1405 if (m->as_integer_class() <= 0) {
1407 }
else if (m->as_integer_class() == 1) {
1417 for (
const auto &it : prime_mul) {
1419 mp_pow_ui(_mod, it.first->as_integer_class(), it.second);
1422 ret_val = _nthroot_mod_prime_power(
1423 rem1, a->as_integer_class(), n->as_integer_class(),
1424 it.first->as_integer_class(), it.second,
true);
1429 _crt_cartesian(roots, rem, moduli);
1433bool powermod(
const Ptr<RCP<const Integer>> &powm,
const RCP<const Integer> &a,
1434 const RCP<const Number> &b,
const RCP<const Integer> &m)
1436 if (is_a<Integer>(*b)) {
1437 integer_class t = down_cast<const Integer &>(*b).as_integer_class();
1438 if (b->is_negative())
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());
1448 }
else if (is_a<Rational>(*b)) {
1449 RCP<const Integer> num, den, r;
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);
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());
1469 const RCP<const Integer> &a,
const RCP<const Number> &b,
1470 const RCP<const Integer> &m)
1472 if (is_a<Integer>(*b)) {
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());
1482 }
else if (is_a<Rational>(*b)) {
1483 RCP<const Integer> num, den, r;
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));
1489 integer_class t = num->as_integer_class();
1490 if (num->is_negative())
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());
1514 throw SymEngineException(
"quadratic_residues: Input must be > 0");
1518 for (integer_class i = integer_class(0); i <= a.
as_int() / 2; i++) {
1522 sort(residue.
begin(), residue.
end());
1538 throw SymEngineException(
1539 "is_quad_residue: Second parameter must be non-zero");
1552 const RCP<const Integer> a1 =
integer(a_final);
1553 const RCP<const Integer> p1 =
integer(p2);
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);
1569 return mp_legendre(a_final, p2) == 1;
1578 integer_class _mod =
mod.as_integer_class();
1582 }
else if (_mod == 1) {
1589 RCP<const Integer> mod2 =
integer(_mod);
1594 for (
const auto &it : prime_mul) {
1595 ret_val = _is_nthroot_mod_prime_power(
1597 it.first->as_integer_class(), it.second);
1607 throw SymEngineException(
"mobius: Integer <= 0");
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;
1616 is_square_free =
false;
1620 if (!is_square_free) {
1622 }
else if (num_prime_factors % 2 == 0) {
1629long mertens(
const unsigned long a)
1632 for (
unsigned long i = 1; i <= a; ++i) {
1649 const integer_class &n)
1651 auto res = ((s - 2) * n * n - (s - 4) * n) / 2;
1666 const integer_class &x)
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));
1680 unsigned long p = 2;
1681 integer_class intone, i, j, m, res;
1686 while ((intone << p) <= n) {
1691 mp_pow_ui(res, m, p);
1697 mp_pow_ui(res, i, p);
1700 if (lowest_exponent) {
Classes and functions relating to the binary operation of addition.
const integer_class & as_integer_class() const
Convert to integer_class.
signed long int as_int() const
Convert to int, raise an exception if it does not fit.
static RCP< const Number > from_mpq(const rational_class &i)
Main namespace for SymEngine package.
integer_class mp_polygonal_number(const integer_class &s, const integer_class &n)
Numeric calculation of the n:th s-gonal number.
bool divides(const Integer &a, const Integer &b)
RCP< const Integer > nextprime(const Integer &a)
int kronecker(const Integer &a, const Integer &n)
Kronecker Function.
void quotient_mod_f(const Ptr< RCP< const Integer > > &q, const Ptr< RCP< const Integer > > &r, const Integer &n, const Integer &d)
void prime_factor_multiplicities(map_integer_uint &primes_mul, const Integer &n)
Find multiplicities of prime factors of n
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.
RCP< const Integer > mod(const Integer &n, const Integer &d)
modulo round toward zero
int mobius(const Integer &a)
Mobius Function.
void lucas2(const Ptr< RCP< const Integer > > &g, const Ptr< RCP< const Integer > > &s, unsigned long n)
Lucas number n and n-1.
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.
RCP< const Integer > lucas(unsigned long n)
Lucas number.
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.
RCP< const Integer > quotient_f(const Integer &n, const Integer &d)
void prime_factors(std::vector< RCP< const Integer > > &prime_list, const Integer &n)
Find prime factors of n
RCP< const Integer > fibonacci(unsigned long n)
Fibonacci number.
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.
bool is_quad_residue(const Integer &a, const Integer &p)
Returns true if 'a' is a quadratic residue of 'p'.
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.
RCP< const Integer > mod_f(const Integer &n, const Integer &d)
modulo round toward -inf
RCP< const Integer > quotient(const Integer &n, const Integer &d)
vec_integer_class quadratic_residues(const Integer &a)
Finds all Quadratic Residues of a Positive Integer.
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>
int jacobi(const Integer &a, const Integer &n)
Jacobi Function.
bool powermod(const Ptr< RCP< const Integer > > &powm, const RCP< const Integer > &a, const RCP< const Number > &b, const RCP< const Integer > &m)
RCP< const Integer > lcm(const Integer &a, const Integer &b)
Least Common Multiple.
RCP< const Integer > binomial(const Integer &n, unsigned long k)
Binomial Coefficient.
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.
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.
int factor_lehman_method(const Ptr< RCP< const Integer > > &f, const Integer &n)
Factor using lehman's methods.
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.
void insert(T1 &m, const T2 &first, const T3 &second)
std::pair< integer_class, integer_class > mp_perfect_power_decomposition(const integer_class &n, bool lowest_exponent)
Decompose a positive integer into perfect powers.
int factor(const Ptr< RCP< const Integer > > &f, const Integer &n, double B1)
int factor_trial_division(const Ptr< RCP< const Integer > > &f, const Integer &n)
void fibonacci2(const Ptr< RCP< const Integer > > &g, const Ptr< RCP< const Integer > > &s, unsigned long n)
Fibonacci n and n-1.
int probab_prime_p(const Integer &a, unsigned reps)
Probabilistic Prime.
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.
bool is_nth_residue(const Integer &a, const Integer &n, const Integer &mod)
Returns true if 'a' is a nth power residue of 'mod'.
int factor_pollard_rho_method(const Ptr< RCP< const Integer > > &f, const Integer &n, unsigned retries)
Factor using Pollard's rho methods.
RCP< const Integer > factorial(unsigned long n)
Factorial.
RCP< const Integer > carmichael(const RCP< const Integer > &n)
Carmichael function.
integer_class mp_principal_polygonal_root(const integer_class &s, const integer_class &x)
Numeric calculation of the principal s-gonal root of x.
void quotient_mod(const Ptr< RCP< const Integer > > &q, const Ptr< RCP< const Integer > > &r, const Integer &n, const Integer &d)
void powermod_list(std::vector< RCP< const Integer > > &pows, const RCP< const Integer > &a, const RCP< const Number > &b, const RCP< const Integer > &m)
std::enable_if< std::is_integral< T >::value, RCP< constInteger > >::type integer(T i)
RCP< const Number > bernoulli(unsigned long n)
int legendre(const Integer &a, const Integer &n)
Legendre Function.
RCP< const Number > harmonic(unsigned long n, long m)
Computes the sum of the inverses of the first perfect mth powers.
int mod_inverse(const Ptr< RCP< const Integer > > &b, const Integer &a, const Integer &m)
inverse modulo
void primitive_root_list(std::vector< RCP< const Integer > > &roots, const Integer &n)
less operator (<) for Integers