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>
104 integer_class
_q,
_r;
139 mp_lucnum2_ui(
g_t,
s_t, n);
184 integer_class
sqrtN = mp_sqrt(
N);
187 throw SymEngineException(
"N too large to factor");
190 while ((p = pi.next_prime()) <=
limit) {
202 throw SymEngineException(
"Require n >= 21 to use lehman method");
212 while ((p = pi.next_prime()) <= mp_get_ui(
u_bound)) {
222 integer_class k, a, b,
l;
227 a = mp_sqrt(4 * k * n);
234 l = a * a - 4 * k * n;
235 if (mp_perfect_square_p(
l)) {
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);
306 state.urandomint(
c,
nm4);
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) {
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) {
399 if (mp_probab_prime_p(
_n, 25) > 0) {
408 throw SymEngineException(
409 "ECM failed to factor the given number");
445 throw SymEngineException(
"N too large to factor");
449 while ((p = pi.next_prime()) <=
limit) {
450 while (
_n % p == 0) {
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
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));
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;
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)) {
640 const integer_class &
e,
bool even =
false)
650 t =
it->as_integer_class();
665 integer_class
pm1 = p - 1;
666 mp_powm(t,
g,
pm1, t);
673 mp_pow_ui(t, p, mp_get_ui(
e));
715 const integer_class &p,
const integer_class &
e,
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));
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) {
797 integer_class
phi = n->as_integer_class(), p;
804 p =
it.first->as_integer_class();
818 integer_class lambda, t, p;
824 p =
it.first->as_integer_class();
832 mp_lcm(lambda, lambda, t);
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();
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);
890 const integer_class &p)
893 integer_class n,
y, b,
q,
pm1, t(1);
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);
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);
979 const integer_class &
g,
const integer_class &n,
980 const integer_class &
q,
const unsigned &k,
981 const integer_class &p)
986 mp_powm(alpha,
g,
_n, p);
994 mp_powm(s, alpha, s, p);
996 for (
unsigned j = 0;
j < m; ++
j) {
1001 for (
unsigned long j = 0;
j < k; ++
j) {
1006 for (
unsigned i = 0;
not found && i < m; ++i) {
1018 mp_powm(t,
g, t, p);
1027 const integer_class &a,
const integer_class &n,
1028 const integer_class &p,
const unsigned k,
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) {
1057 integer_class
h,
q,
qt,
z, v, x,
s1 = s;
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);
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);
1137 mp_powm(t,
g, t,
pk);
1139 for (
unsigned j = 0;
j < m; ++
j) {
1142 mp_fdiv_r(root, root,
pk);
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);
1170 const integer_class &a,
const integer_class &n,
1171 const integer_class &p,
const unsigned k,
1174 integer_class
pk, root;
1178 integer_class
r = n, t, s,
pc,
pj;
1179 pk = integer_class(1) << k;
1189 if (
c > 0
and a % 4 == 3) {
1201 t = integer_class(1) << (k - 2);
1202 pc = integer_class(1) <<
c;
1208 mp_powm(root, a, s,
pk);
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) {
1252 mp_pow_ui(
pk, p, k);
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
1280 mp_pow_ui(
pm, p, m);
1290 mp_pow_ui(
pm, p, m);
1293 mp_pow_ui(
pkm, p, k - m);
1296 root =
it->as_integer_class();
1297 for (
unsigned long i = 0; i <
pm; ++i) {
1308bool _is_nthroot_mod_prime_power(
const integer_class &a,
const integer_class &n,
1309 const integer_class &p,
const unsigned k)
1322 if (
c > 0
and a % 4 == 3) {
1338 t = integer_class(1) << (
c + 2);
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)) {
1375 if (
mod->as_integer_class() <= 0) {
1377 }
else if (
mod->as_integer_class() == 1) {
1389 mp_pow_ui(
_mod,
it.first->as_integer_class(),
it.second);
1392 rem, a->as_integer_class(), n->as_integer_class(),
1393 it.first->as_integer_class(),
it.second,
false);
1405 if (m->as_integer_class() <= 0) {
1407 }
else if (m->as_integer_class() == 1) {
1419 mp_pow_ui(
_mod,
it.first->as_integer_class(),
it.second);
1423 rem1, a->as_integer_class(), n->as_integer_class(),
1424 it.first->as_integer_class(),
it.second,
true);
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());
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());
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());
1485 if (
den->is_negative()) {
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++) {
1538 throw SymEngineException(
1539 "is_quad_residue: Second parameter must be non-zero");
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);
1578 integer_class
_mod =
mod.as_integer_class();
1582 }
else if (
_mod == 1) {
1595 ret_val = _is_nthroot_mod_prime_power(
1597 it.first->as_integer_class(),
it.second);
1607 throw SymEngineException(
"mobius: Integer <= 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;
1686 while ((
intone << p) <= n) {
1691 mp_pow_ui(
res, m, p);
1697 mp_pow_ui(
res, i, p);
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
void hash_combine(hash_t &seed, const T &v)
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