const integer_class &b)

{

/*boost::multiprecision doesn't have a built-in fdiv_qr (floored division).

// must copy a and b before calling divide_qr because a or b may refer to

// the same

// object as q or r, causing incorrect results

integer_class a_cpy = a, b_cpy = b;

bool neg_quotient = ((a < 0 && b > 0) || (a > 0 && b < 0)) ? true : false;

boost::multiprecision::divide_qr(a_cpy, b_cpy, q, r);

// floor the quotient if necessary

if (neg_quotient && r != 0) {

q -= 1;

}

// remainder should have same sign as divisor

if ((b_cpy > 0 && r < 0) || (b_cpy < 0 && r > 0)) {

r += b_cpy;

return;

}

const integer_class &b)

{

integer_class a_cpy = a, b_cpy = b;

bool pos_quotient = ((a < 0 && b < 0) || (a > 0 && b > 0)) ? true : false;

boost::multiprecision::divide_qr(a_cpy, b_cpy, q, r);

// ceil the quotient if necessary

if (pos_quotient && r != 0) {

q += 1;

}

// remainder should have opposite sign as divisor

if ((b_cpy > 0 && r > 0) || (b_cpy < 0 && r < 0)) {

r -= b_cpy;

return;

}

const integer_class &a, const integer_class &b)

{

integer_class this_s(1);

integer_class this_t(0);

integer_class next_s(0);

integer_class next_t(1);

integer_class this_r(a);

integer_class next_r(b);

integer_class q;

while (next_r != 0) {

// should use truncated division, so use

// boost::multiprecision::divide_qr

// beware of overwriting this_r during internal operations of divide_qr

// copy it first

integer_class this_r_cpy = this_r;

boost::multiprecision::divide_qr(this_r_cpy, next_r, q, this_r);

this_s -= q * next_s;

this_t -= q * next_t;

std::swap(this_s, next_s);

std::swap(this_t, next_t);

std::swap(this_r, next_r);

}

// normalize the gcd, s and t

if (this_r < 0) {

this_r *= -1;

this_s *= -1;

this_t *= -1;

}

gcd = std::move(this_r);

s = std::move(this_s);

t = std::move(this_t);

const integer_class &m)

{

integer_class gcd, s, t;

mp_gcdext(gcd, s, t, a, m);

if (gcd != 1) {

res = 0;

return false;

} else {

mp_fdiv_r(s, s, m); // reduce s modulo m. undefined behavior when m ==

// 0, so don't need to check

if (s < 0) {

s += mp_abs(m);

} // give the canonical representative of s

res = s;

return true;

}

{

integer_class res = a % mod;

if (res < 0) {

res += mod;

}

return res;

{

integer_class num = numerator(i); // copy

integer_class den = denominator(i); // copy

num = pow(num, n);

den = pow(den, n);

res = rational_class(std::move(num), std::move(den));

const integer_class &exp, const integer_class &m)

{

// if exp is negative, interpret as follows

// base**(exp) mod m == (base**(-1))**abs(exp) mod m

// == (base**(-1) mod m) ** abs(exp) mod m

// where base**(-1) mod m is the modular inverse

if (exp < 0) {

integer_class base_inverse;

if (!mp_invert(base_inverse, base, m)) {

throw SymEngine::SymEngineException("negative exponent undefined "

"in powm if base is not "

"invertible mod m");

}

res = boost::multiprecision::powm(base_inverse, mp_abs(exp), m);

return;

} else {

res = boost::multiprecision::powm(base, exp, m);

// boost's powm calculates base**exp % m, but uses truncated

// modulus, e.g. powm(-2,3,5) == -3. We want powm(-2,3,5) == 2

if (res < 0) {

res += m;

}

}

integer_class &x)

{

SYMENGINE_ASSERT(n > 1);

unsigned long m = n - 1;

integer_class &&x_m = pow(x, m);

return integer_class((integer_class(m * x) + integer_class(i / x_m)) / n);

const unsigned long n)

{

integer_class x

= 1; // TODO: make a better starting guess based on (number of bits)/n

integer_class y = step(n, i, x);

do {

x = y;

y = step(n, i, x);

} while (y < x);

res = x;

if (pow(x, n) == i) {

return true;

}

return false;

{

if (n == 0) {

throw std::runtime_error("0th root is undefined");

}

if (n == 1) {

res = i;

return true;

}

if (i == 0) {

res = 0;

return true;

}

if (i > 0) {

return positive_root(res, i, n);

}

if (i < 0 && (n % 2 == 0)) {

throw std::runtime_error("even root of a negative is non-real");

}

bool b = positive_root(res, -i, n);

res *= -1;

return b;

{

// as of 11/1/2016, boost::multiprecision::sqrt() is buggy:

// https://svn.boost.org/trac/boost/ticket/12559

// implement with mp_root for now

integer_class res;

mp_root(res, i, 2);

return res;

unsigned long n)

{

mp_root(a, i, n);

integer_class p = pow(a, n);

;

b = i - p;

{

a = mp_sqrt(i);

b = i - boost::multiprecision::pow(a, 2);

{

if (i % 2 == 0)

return (i == 2);

return miller_rabin_test(i, retries);

{

// simple implementation: just check all odds bigger than i for primality

if (i < 2) {

res = 2;

return;

}

integer_class candidate;

candidate = (i % 2 == 0) ? i + 1 : i + 2;

// Knuth recommends 25 trials for a pretty strong likelihood that candidate

// is prime

while (!mp_probab_prime_p(candidate, 25)) {

candidate += 2;

}

res = std::move(candidate);

{

if (i == 0) {

return ULONG_MAX;

}

return find_lsb(i, {});

integer_class data[2][2]; // data[1][0] is row 1, column 0 entry of matrix

two_by_two_matrix(integer_class a, integer_class b, integer_class c,

integer_class d)

: data{{a, b}, {c, d}}

{

}

two_by_two_matrix() : data{{0, 0}, {0, 0}} {}

two_by_two_matrix(const two_by_two_matrix &other) = default;

two_by_two_matrix &operator=(const two_by_two_matrix &other)

{

this->data[0][0] = other.data[0][0];

this->data[0][1] = other.data[0][1];

this->data[1][0] = other.data[1][0];

this->data[1][1] = other.data[1][1];

return *this;

}

two_by_two_matrix operator*(const two_by_two_matrix &other)

{

two_by_two_matrix res;

res.data[0][0] = this->data[0][0] * other.data[0][0]

+ this->data[0][1] * other.data[1][0];

res.data[0][1] = this->data[0][0] * other.data[0][1]

+ this->data[0][1] * other.data[1][1];

res.data[1][0] = this->data[1][0] * other.data[0][0]

+ this->data[1][1] * other.data[1][0];

res.data[1][1] = this->data[1][0] * other.data[0][1]

+ this->data[1][1] * other.data[1][1];

return res;

}

// recursive repeated squaring

two_by_two_matrix pow(unsigned long n)

{

if (n == 0) {

return two_by_two_matrix::identity();

}

if (n == 1) {

return *this;

}

if (n == 2) {

return (*this) * (*this);

}

if (n % 2 == 0) {

return (this->pow(n / 2)).pow(2);

}

return ((this->pow((n - 1) / 2)).pow(2)) * (*this);

}

static two_by_two_matrix identity()

{

return two_by_two_matrix(1, 0, 0, 1);

}

{

two_by_two_matrix x(1, 1, 1, 0);

return x.pow(n);

{

// reference: https://www.nayuki.io/page/fast-fibonacci-algorithms

res = fib_matrix(n).data[0][1];

{

// reference: https://www.nayuki.io/page/fast-fibonacci-algorithms

two_by_two_matrix result_matrix = fib_matrix(n);

a = result_matrix.data[0][1];

b = result_matrix.data[1][1];

{

two_by_two_matrix multiplier(1, 1, 1, 0);

two_by_two_matrix start(1, 0, 2, 0);

return multiplier.pow(n) * start;

{

// implementation based on the following fact:

// [[1,1],[1,0]]^(n-1)*[2,0,1,0] = [[L(n+1),0][L(n),0]]

// where L(n) is the nth Lucas number

res = luc_matrix(n).data[1][0];

{

if (n == 0) {

throw std::runtime_error("index of lucas number cannot be negative");

}

two_by_two_matrix result_matrix = luc_matrix(n - 1);

a = result_matrix.data[0][0];

b = result_matrix.data[1][0];

{

// couldn't make boost's template version of factorial work,

// so implement slow, naive version for now

res = 1;

for (unsigned long i = 2; i <= n; ++i) {

res *= i;

}

{

// slow, naive implementation

integer_class x = n - r;

res = 1;

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

res *= x + i;

res /= i;

}

{

if (i == 0 || i == 1 || i == -1) {

return true;

}

// if i == a**k, with k == pq for some integers p,q

// then i == a**(pq) == (a**p)**q == (a**q)**p

// hence i is a pth power and a qth power

// Hence if i == p**k, then i is a pth power for any p

// in the prime factorization of k, and it suffices

// to check whether i is a prime power.

// the largest possible prime p would arise from the

// case where i is a prime power of two

// so check all prime roots up to log(i) base 2.

unsigned long max = std::ilogb(i.convert_to<double>());

// treat case p=2 separately b/c mp_root throws exception

// with an even root of a negative

if (mp_perfect_square_p(i)) {

return true;

}

integer_class p(2);

integer_class root(0);

while (true) {

mp_nextprime(p, p);

if (p > max) {

return false;

}

if (mp_root(root, i, p.convert_to<unsigned long>())) {

return true;

}

}

{

if (i < 0) {

return false;

}

integer_class root;

return mp_root(root, i, 2);

{

integer_class res = 1;

Sieve::iterator pi(static_cast<unsigned>(n));

unsigned int p;

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

res *= p;

}

return res;

{

integer_class res;

mp_powm(res, a, integer_class((n - 1) / 2), n);

return res <= 1 ? res.convert_to<int>() : -1;

{

// https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol

if (a == 1) {

return 1;

}

integer_class num = a;

integer_class den = n;

// (1) Reduce the "numerator" modulo the "denominator"

num = fmod(num, den);

// (2) Extract any factors of 2 from the "numerator"

unsigned long factors_of_two = 0;

while (num % 2 == 0 && num != 0) {

num /= 2; // use a shift instead of division here? faster?

++factors_of_two;

}

int product_of_twos = 1; // (2 | den)**factors_of_two

// (2 | den) is -1 iff den % 8 = 3 or den % 8 = 5.

// If factors_of_two is odd and (2 | den) is -1,

// then (2 | den)**factors_of_two is -1.

// Otherwise, (2 | den)**factors_of_two is 1.

int den_mod_8 = fmod(den, 8).convert_to<int>();

if ((factors_of_two % 2 == 1) && (den_mod_8 == 3 || den_mod_8 == 5)) {

product_of_twos = -1;

}

// (3) If the remaining "numerator" (after extraction of twos) is 1,

// then the remaining jacobi symbol is 1.

// If the "numerator" and "denominator" are not coprime, result is 0.

if (num == 1) {

return product_of_twos;

}

if (boost::multiprecision::gcd(num, den) != 1) {

return 0;

}

// (4) Otherwise, the "numerator" and "denominator" are now odd

// positive coprime integers, so we can flip the symbol

// with quadratic reciprocity, then return to step (1)

// (num | den) == (den | num), unless num % 4 and den %4 are both

// 3, in which case (num | den) == (-1) * (den | num)

int quadratic_reciprocity_factor = 1;

if ((fmod(num, 4) == 3) && (fmod(den, 4) == 3)) {

quadratic_reciprocity_factor = -1;

}

return product_of_twos * quadratic_reciprocity_factor

* unchecked_jacobi(den, num);

{

if (n < 0) {

throw std::runtime_error("jacobi denominator must be positive");

}

if (n % 2 == 0) {

throw std::runtime_error("jacobi denominator must be odd");

}

return unchecked_jacobi(a, n);

{

/*

if (n == 0) {

throw std::runtime_error("second arg of Kronecker cannot be zero");

}

// Compute (a | u)

int kr_a_u = 1;

if (n.sign() == -1 && a < 0) {

kr_a_u = -1;

}

// Compute m, j

integer_class m = boost::multiprecision::abs(n);

unsigned long j = 0;

while (m % 2 == 0 && m != 0) { // while m is even

m /= 2; // implement with shift? faster?

++j;

}

// if n is even, compute (a | 2)**j

int kr_a_2 = 0;

int kr_a_2_to_j = 0;

int a_mod_8 = fmod(a, 8).convert_to<int>();

if (a % 2 != 0) {

kr_a_2 = (a_mod_8 == 1 || a_mod_8 == 7) ? 1 : -1;

kr_a_2_to_j = (kr_a_2 == -1 && (j % 2 != 0))

? -1

: 1; //(-1)**odd == -1; (-1)**even=(1)**integer=1

}

if (n % 2 == 0) {

return kr_a_u * kr_a_2_to_j * unchecked_jacobi(a, m);

} else {

return kr_a_u * unchecked_jacobi(a, m);

}