{

integer_class g;

mp_gcd(g, a.as_integer_class(), b.as_integer_class());

return integer(std::move(g));

const Ptr<RCP<const Integer>> &t, const Integer &a,

const Integer &b)

{

integer_class g_, s_, t_;

mp_gcdext(g_, s_, t_, a.as_integer_class(), b.as_integer_class());

*g = integer(std::move(g_));

*s = integer(std::move(s_));

*t = integer(std::move(t_));

{

integer_class c;

mp_lcm(c, a.as_integer_class(), b.as_integer_class());

return integer(std::move(c));

const Integer &m)

{

int ret_val;

integer_class inv_t;

ret_val = mp_invert(inv_t, a.as_integer_class(), m.as_integer_class());

*b = integer(std::move(inv_t));

return ret_val;

const Ptr<RCP<const Integer>> &r, const Integer &n,

const Integer &d)

{

integer_class _q, _r;

mp_tdiv_qr(_q, _r, n.as_integer_class(), d.as_integer_class());

*q = integer(std::move(_q));

*r = integer(std::move(_r));

{

integer_class q;

mp_fdiv_r(q, n.as_integer_class(), d.as_integer_class());

return integer(std::move(q));

{

integer_class q;

mp_fdiv_q(q, n.as_integer_class(), d.as_integer_class());

return integer(std::move(q));

const Ptr<RCP<const Integer>> &r, const Integer &n,

const Integer &d)

{

integer_class _q, _r;

mp_fdiv_qr(_q, _r, n.as_integer_class(), d.as_integer_class());

*q = integer(std::move(_q));

*r = integer(std::move(_r));

{

integer_class f;

mp_fib_ui(f, n);

return integer(std::move(f));

const Ptr<RCP<const Integer>> &s, unsigned long n)

{

integer_class g_t;

integer_class s_t;

mp_fib2_ui(g_t, s_t, n);

*g = integer(std::move(g_t));

*s = integer(std::move(s_t));

{

integer_class f;

mp_lucnum_ui(f, n);

return integer(std::move(f));

unsigned long n)

{

integer_class g_t;

integer_class s_t;

mp_lucnum2_ui(g_t, s_t, n);

*g = integer(std::move(g_t));

*s = integer(std::move(s_t));

{

integer_class f;

mp_bin_ui(f, n.as_integer_class(), k);

return integer(std::move(f));

{

integer_class f;

mp_fac_ui(f, n);

return integer(std::move(f));

{

integer_class c;

mp_nextprime(c, a.as_integer_class());

return integer(std::move(c));

{

integer_class sqrtN = mp_sqrt(N);

unsigned long limit = mp_get_ui(sqrtN);

if (limit > std::numeric_limits<unsigned>::max())

throw SymEngineException("N too large to factor");

Sieve::iterator pi(numeric_cast<unsigned>(limit));

unsigned p;

while ((p = pi.next_prime()) <= limit) {

if (N % p == 0) {

factor = p;

return 1;

}

}

return 0;

{

if (n < 21)

throw SymEngineException("Require n >= 21 to use lehman method");

int ret_val = 0;

integer_class u_bound;

mp_root(u_bound, n, 3);

u_bound = u_bound + 1;

Sieve::iterator pi(numeric_cast<unsigned>(mp_get_ui(u_bound)));

unsigned p;

while ((p = pi.next_prime()) <= mp_get_ui(u_bound)) {

if (n % p == 0) {

rop = n / p;

ret_val = 1;

break;

}

}

if (not ret_val) {

integer_class k, a, b, l;

k = 1;

while (k <= u_bound) {

a = mp_sqrt(4 * k * n);

mp_root(b, n, 6);

mp_root(l, k, 2);

b = b / (4 * l);

b = b + a;

while (a <= b) {

l = a * a - 4 * k * n;

if (mp_perfect_square_p(l)) {

b = a + mp_sqrt(l);

mp_gcd(rop, n, b);

ret_val = 1;

break;

}

a = a + 1;

}

if (ret_val)

break;

k = k + 1;

}

}

return ret_val;

{

int ret_val;

integer_class rop;

ret_val = _factor_lehman_method(rop, n.as_integer_class());

*f = integer(std::move(rop));

return ret_val;

const integer_class &c, unsigned B)

{

if (n < 4 or B < 3)

throw SymEngineException(

"Require n > 3 and B > 2 to use Pollard's p-1 method");

integer_class m, _c;

_c = c;

Sieve::iterator pi(B);

unsigned p;

while ((p = pi.next_prime()) <= B) {

m = 1;

// calculate log(p, B), this can be improved

while (m <= B / p) {

m = m * p;

}

mp_powm(_c, _c, m, n);

}

_c = _c - 1;

mp_gcd(rop, _c, n);

if (rop == 1 or rop == n)

return 0;

else

return 1;

const Integer &n, unsigned B, unsigned retries)

{

int ret_val = 0;

integer_class rop, nm4, c;

mp_randstate state;

nm4 = n.as_integer_class() - 4;

for (unsigned i = 0; i < retries and ret_val == 0; ++i) {

state.urandomint(c, nm4);

c += 2;

ret_val = _factor_pollard_pm1_method(rop, n.as_integer_class(), c, B);

}

if (ret_val != 0)

*f = integer(std::move(rop));

return ret_val;

const integer_class &a, const integer_class &s,

unsigned steps = 10000)

{

if (n < 5)

throw SymEngineException("Require n > 4 to use pollard's-rho method");

integer_class u, v, g, m;

u = s;

v = s;

for (unsigned i = 0; i < steps; ++i) {

u = (u * u + a) % n;

v = (v * v + a) % n;

v = (v * v + a) % n;

m = u - v;

mp_gcd(g, m, n);

if (g == n)

return 0;

if (g == 1)

continue;

rop = g;

return 1;

}

return 0;

const Integer &n, unsigned retries)

{

int ret_val = 0;

integer_class rop, nm1, nm4, a, s;

mp_randstate state;

nm1 = n.as_integer_class() - 1;

nm4 = n.as_integer_class() - 4;

for (unsigned i = 0; i < retries and ret_val == 0; ++i) {

state.urandomint(a, nm1);

state.urandomint(s, nm4);

s += 1;

ret_val = _factor_pollard_rho_method(rop, n.as_integer_class(), a, s);

}

if (ret_val != 0)

*f = integer(std::move(rop));

return ret_val;

{

int ret_val = 0;

integer_class _n, _f;

_n = n.as_integer_class();

#ifdef HAVE_SYMENGINE_ECM

if (mp_perfect_power_p(_n)) {

unsigned long int i = 1;

integer_class m, rem;

rem = 1; // Any non zero number

m = 2; // set `m` to 2**i, i = 1 at the begining

// calculate log2n, this can be improved

for (; m < _n; ++i)

m = m * 2;

// eventually `rem` = 0 zero as `n` is a perfect power. `f_t` will

// be set to a factor of `n` when that happens

while (i > 1 and rem != 0) {

mp_rootrem(_f, rem, _n, i);

--i;

}

ret_val = 1;

} else {

if (mp_probab_prime_p(_n, 25) > 0) { // most probably, n is a prime

ret_val = 0;

_f = _n;

} else {

for (int i = 0; i < 10 and not ret_val; ++i)

ret_val = ecm_factor(get_mpz_t(_f), get_mpz_t(_n), B1, nullptr);

mp_demote(_f);

if (not ret_val)

throw SymEngineException(

"ECM failed to factor the given number");

}

}

#else

// B1 is discarded if gmp-ecm is not installed

ret_val = _factor_trial_division_sieve(_f, _n);

#endif // HAVE_SYMENGINE_ECM

*f = integer(std::move(_f));

return ret_val;

{

int ret_val;

integer_class factor;

ret_val = _factor_trial_division_sieve(factor, n.as_integer_class());

if (ret_val == 1)

*f = integer(std::move(factor));

return ret_val;

const Integer &n)

{

integer_class sqrtN;

integer_class _n = n.as_integer_class();

if (_n == 0)

return;

if (_n < 0)

_n *= -1;

sqrtN = mp_sqrt(_n);

auto limit = mp_get_ui(sqrtN);

if (not mp_fits_ulong_p(sqrtN)

or limit > std::numeric_limits<unsigned>::max())

throw SymEngineException("N too large to factor");

Sieve::iterator pi(numeric_cast<unsigned>(limit));

unsigned p;

while ((p = pi.next_prime()) <= limit) {

while (_n % p == 0) {

prime_list.push_back(integer(p));

_n = _n / p;

}

if (_n == 1)

break;

}

if (not(_n == 1))

prime_list.push_back(integer(std::move(_n)));

{

integer_class sqrtN;

integer_class _n = n.as_integer_class();

unsigned count;

if (_n == 0)

return;

if (_n < 0)

_n *= -1;

sqrtN = mp_sqrt(_n);

auto limit = mp_get_ui(sqrtN);

if (not mp_fits_ulong_p(sqrtN)

or limit > std::numeric_limits<unsigned>::max())

throw SymEngineException("N too large to factor");

Sieve::iterator pi(numeric_cast<unsigned>(limit));

unsigned p;

while ((p = pi.next_prime()) <= limit) {

count = 0;

while (_n % p == 0) { // when a prime factor is found, we divide

++count; // _n by that prime as much as we can

_n = _n / p;

}

if (count > 0) {

insert(primes_mul, integer(p), count);

if (_n == 1)

break;

}

}

if (not(_n == 1))

insert(primes_mul, integer(std::move(_n)), 1);

{

#ifdef HAVE_SYMENGINE_ARB

fmpq_t res;

fmpq_init(res);

bernoulli_fmpq_ui(res, n);

mpq_t a;

mpq_init(a);

fmpq_get_mpq(a, res);

rational_class b(a);

fmpq_clear(res);

mpq_clear(a);

return Rational::from_mpq(std::move(b));

#else

// TODO: implement a faster algorithm

std::vector<rational_class> v(n + 1);

for (unsigned m = 0; m <= n; ++m) {

v[m] = rational_class(1u, m + 1);

for (unsigned j = m; j >= 1; --j) {

v[j - 1] = j * (v[j - 1] - v[j]);

}

}

return Rational::from_mpq(v[0]);

#endif

{

rational_class res(0);

if (m == 1) {

for (unsigned i = 1; i <= n; ++i) {

res += rational_class(1u, i);

}

return Rational::from_mpq(res);

} else {

for (unsigned i = 1; i <= n; ++i) {

if (m > 0) {

rational_class t(1u, i);

#if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP

mp_pow_ui(get_den(t), get_den(t), m);

#else

mp_pow_ui(t, t, m);

#endif

res += t;

} else {

integer_class t(i);

mp_pow_ui(t, t, static_cast<unsigned long>(-m));

res += t;

}

}

return Rational::from_mpq(res);

}

const std::vector<RCP<const Integer>> &rem,

const std::vector<RCP<const Integer>> &mod)

{

if (mod.size() > rem.size())

throw SymEngineException("Too few remainders");

if (mod.size() == 0)

throw SymEngineException("Moduli vector cannot be empty");

integer_class m, r, g, s, t;

m = mod[0]->as_integer_class();

r = rem[0]->as_integer_class();

for (unsigned i = 1; i < mod.size(); ++i) {

mp_gcdext(g, s, t, m, mod[i]->as_integer_class());

// g = s * m + t * mod[i]

t = rem[i]->as_integer_class() - r;

if (not mp_divisible_p(t, g))

return false;

r += m * s * (t / g); // r += m * (m**-1 mod[i]/g)* (rem[i] - r) / g

m *= mod[i]->as_integer_class() / g;

mp_fdiv_r(r, r, m);

}

*R = integer(std::move(r));

return true;

const std::vector<std::vector<RCP<const Integer>>> &rem,

const std::vector<RCP<const Integer>> &mod)

{

if (mod.size() > rem.size())

throw SymEngineException("Too few remainders");

if (mod.size() == 0)

throw SymEngineException("Moduli vector cannot be empty");

integer_class m, _m, r, s, t;

m = mod[0]->as_integer_class();

R = rem[0];

for (unsigned i = 1; i < mod.size(); ++i) {

std::vector<RCP<const Integer>> rem2;

mp_invert(s, m, mod[i]->as_integer_class());

_m = m;

m *= mod[i]->as_integer_class();

for (auto &elem : R) {

for (auto &_k : rem[i]) {

r = elem->as_integer_class();

r += _m * s * (_k->as_integer_class() - r);

mp_fdiv_r(r, r, m);

rem2.push_back(integer(r));

}

}

R = rem2;

}

{

if (n < 2)

return false;

integer_class _n = n, temp;

e = 1;

unsigned i = 2;

while (mp_perfect_power_p(_n) and _n >= 2) {

if (mp_root(temp, _n, i)) {

e *= i;

_n = temp;

} else {

++i;

}

}

if (mp_probab_prime_p(_n, 25)) {

p = _n;

return true;

}

return false;

const integer_class &e, bool even = false)

{

std::vector<RCP<const Integer>> primes;

prime_factors(primes, *integer(p - 1));

integer_class t;

g = 2;

while (g < p) {

bool root = true;

for (const auto &it : primes) {

t = it->as_integer_class();

t = (p - 1) / t;

mp_powm(t, g, t, p);

if (t == 1) { // If g**(p-1)/q is 1 then g is not a primitive root.

root = false;

break;

}

}

if (root)

break;

++g;

}

if (e > 1) {

t = p * p;

integer_class pm1 = p - 1;

mp_powm(t, g, pm1, t);

if (t == 1) { // If g**(p-1) mod (p**2) == 1 then g + p is a primitive

// root.

g += p;

}

}

if (even and g % 2 == 0) {

mp_pow_ui(t, p, mp_get_ui(e));

g += t; // If g is even then root of 2*p**e is g + p**e.

}

{

integer_class _n = n.as_integer_class();

if (_n < 0)

_n = -_n;

if (_n <= 1)

return false;

if (_n < 5) {

*g = integer(_n - 1);

return true;

}

bool even = false;

if (_n % 2 == 0) {

if (_n % 4 == 0) {

return false; // If n mod 4 == 0 and n > 4, then no primitive roots.

}

_n /= 2;

even = true;

}

integer_class p, e;

if (not _prime_power(p, e, _n))

return false;

_primitive_root(_n, p, e, even);

*g = integer(std::move(_n));

return true;

const integer_class &p, const integer_class &e,

bool even = false)

{

integer_class g, h, d, t, pe2, n, pm1;

_primitive_root(g, p, integer_class(1),

false); // Find one primitive root for p.

h = 1;

pm1 = p - 1;

// Generate other primitive roots for p. h = g**i and gcd(i, p-1) = 1.

// Ref[2]

mp_pow_ui(n, p, mp_get_ui(e));

for (unsigned long i = 1; i < p; ++i) {

h *= g;

h %= p;

mp_gcd(d, pm1, integer_class(i));

if (d == 1) {

if (e == 1) {

if (even and h % 2 == 0)

roots.push_back(integer(h + n));

else

roots.push_back(integer(h));

} else {

integer_class pp = p * p;

// Find d such that (h + d*p)**(p-1) mod (p**2) == 1. Ref[1]

// h**(p-1) - 1 = d*p*h**(p-2)

// d = (h - h**(2-p)) / p

t = 2 - p;

mp_powm(d, h, t, pp);

d = ((h - d) / p + p) % p;

t = h;

// t = h + i * p + j * p * p and i != d

mp_pow_ui(pe2, p, mp_get_ui(e) - 2);

for (unsigned long j = 0; j < pe2; ++j) {

for (unsigned long i = 0; i < p; ++i) {

if (i != d) {

if (even and t % 2 == 0)

roots.push_back(integer(t + n));

else

roots.push_back(integer(t));

}

t += p;

}

}

}

}

}

const Integer &n)

{

integer_class _n = n.as_integer_class();

if (_n < 0)

_n = -_n;

if (_n <= 1)

return;

if (_n < 5) {

roots.push_back(integer(_n - 1));

return;

}

bool even = false;

if (_n % 2 == 0) {

if (_n % 4 == 0) {

return; // If n%4 == 0 and n > 4, then no primitive roots.

}

_n /= 2;

even = true;

}

integer_class p, e;

if (not _prime_power(p, e, _n))

return;

_primitive_root_list(roots, p, e, even);

std::sort(roots.begin(), roots.end(), SymEngine::RCPIntegerKeyLess());

return;

{

if (n->is_zero())

return integer(1);

integer_class phi = n->as_integer_class(), p;

if (phi < 0)

phi = -phi;

map_integer_uint prime_mul;

prime_factor_multiplicities(prime_mul, *n);

for (const auto &it : prime_mul) {

p = it.first->as_integer_class();

mp_divexact(phi, phi, p);

// phi is exactly divisible by p.

phi *= p - 1;

}

return integer(std::move(phi));

{

if (n->is_zero())

return integer(1);

map_integer_uint prime_mul;

integer_class lambda, t, p;

unsigned multiplicity;

prime_factor_multiplicities(prime_mul, *n);

lambda = 1;

for (const auto &it : prime_mul) {

p = it.first->as_integer_class();

multiplicity = it.second;

if (p == 2

and multiplicity

> 2) { // For powers of 2 greater than 4 divide by 2.

multiplicity--;

}

t = p - 1;

mp_lcm(lambda, lambda, t);

mp_pow_ui(t, p, multiplicity - 1);

// lambda and p are relatively prime.

lambda = lambda * t;

}

return integer(std::move(lambda));

const RCP<const Integer> &a,

const RCP<const Integer> &n)

{

integer_class order, p, t;

integer_class _a = a->as_integer_class(),

_n = mp_abs(n->as_integer_class());

mp_gcd(t, _a, _n);

if (t != 1)

return false;

RCP<const Integer> lambda = carmichael(n);

map_integer_uint prime_mul;

prime_factor_multiplicities(prime_mul, *lambda);

_a %= _n;

order = lambda->as_integer_class();

for (const auto &it : prime_mul) {

p = it.first->as_integer_class();

mp_pow_ui(t, p, it.second);

mp_divexact(order, order, t);

mp_powm(t, _a, order, _n);

while (t != 1) {

mp_powm(t, t, p, _n);

order *= p;

}

}

*o = integer(std::move(order));

return true;

const integer_class &p)

{

mp_randstate state;

integer_class n, y, b, q, pm1, t(1);

pm1 = p - 1;

unsigned e, m;

e = numeric_cast<unsigned>(mp_scan1(pm1));

q = pm1 >> e; // p - 1 = 2**e*q

while (t != -1) {

state.urandomint(n, p);

t = mp_legendre(n, p);

}

mp_powm(y, n, q, p); // y = n**q mod p

mp_powm(b, a, q, p); // b = a**q mod p

t = (q + 1) / 2;

mp_powm(rop, a, t, p); // rop = a**((q + 1) / 2) mod p

while (b != 1) {

m = 0;

t = b;

while (t != 1) {

mp_powm(t, t, integer_class(2), p);

++m; // t = t**2 = b**2**(m)

}

if (m == e)

return false;

mp_pow_ui(q, integer_class(2), e - m - 1); // q = 2**(e - m - 1)

mp_powm(t, y, q, p); // t = y**(2**(e - m - 1))

mp_powm(y, t, integer_class(2), p); // y = t**2

e = m;

rop = (rop * t) % p;

b = (b * y) % p;

}

return true;

const integer_class &p)

{

if (p == 2) {

rop = a % p;

return true;

}

int l = mp_legendre(a, p);

integer_class t;

if (l == -1) {

return false;

} else if (l == 0) {

rop = 0;

} else if (p % 4 == 3) {

t = (p + 1) / 4;

mp_powm(rop, a, t, p);

} else if (p % 8 == 5) {

t = (p - 1) / 4;

mp_powm(t, a, t, p);

if (t == 1) {

t = (p + 3) / 8;

mp_powm(rop, a, t, p);

} else {

t = (p - 5) / 8;

integer_class t1 = 4 * a;

mp_powm(t, t1, t, p);

rop = (2 * a * t) % p;

}

} else {

if (p < 10000) { // If p < 10000, brute force is faster.

integer_class sq = integer_class(1), _a;

mp_fdiv_r(_a, a, p);

for (unsigned i = 1; i < p; ++i) {

if (sq == _a) {

rop = i;

return true;

}

sq += 2 * i + 1;

mp_fdiv_r(sq, sq, p);

}

return false;

} else {

return _sqrt_mod_tonelli_shanks(rop, a, p);

}

}

return true;

const integer_class &g, const integer_class &n,

const integer_class &q, const unsigned &k,

const integer_class &p)

{

log = 0;

integer_class gamma = a, alpha, _n, t, beta, qj(1), m, l;

_n = n / q;

mp_powm(alpha, g, _n, p);

mp_sqrtrem(m, t, q);

if (t != 0)

++m; // m = ceiling(sqrt(q)).

map_integer_uint

table; // Table for lookup in baby-step giant-step algorithm

integer_class alpha_j(1), d, s;

s = -m;

mp_powm(s, alpha, s, p);

for (unsigned j = 0; j < m; ++j) {

insert(table, integer(alpha_j), j);

alpha_j = (alpha_j * alpha) % p;

}