Program Listing for File mp_boost.cpp¶
↰ Return to documentation for file (symengine/symengine/mp_boost.cpp)
#include "mp_class.h"
#include <boost/multiprecision/miller_rabin.hpp>
#include <boost/mpl/int.hpp>
#include <utility>
#include <cmath>
#include <symengine/symengine_assert.h>
#include <symengine/symengine_exception.h>
using boost::multiprecision::numerator;
using boost::multiprecision::denominator;
using boost::multiprecision::miller_rabin_test;
using boost::multiprecision::detail::find_lsb;
using boost::mpl::int_;
#if SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
namespace SymEngine
{
integer_class pow(const integer_class &a, unsigned long b)
{
return boost::multiprecision::pow(a, numeric_cast<unsigned>(b));
}
void mp_fdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
const integer_class &b)
{
/*boost::multiprecision doesn't have a built-in fdiv_qr (floored division).
Its divide_qr uses truncated division, as does its
modulus operator. Thus, using boost::multiprecision::divide_qr we get:
divide_qr(-5, 3, quo, rem) //quo == -1, rem == -2
divide_qr(5, -3, quo, rem) //quo == -1, rem == 2
but we want:
mp_fdiv_r(quo, rem, -5, 3) //quo == -2, rem == 1
mp_fdiv_r(quo, rem, 5, -3) //rem == -2, rem == -1
The problem arises only when the quotient is negative. To convert
a truncated result into a floored result in this case, simply subtract
one from the truncated quotient and add the divisor to the truncated
remainder.
*/
// 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;
}
}
void mp_cdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
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;
}
}
void mp_gcdext(integer_class &gcd, integer_class &s, integer_class &t,
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);
}
bool mp_invert(integer_class &res, const integer_class &a,
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;
}
}
// floored modulus
integer_class fmod(const integer_class &a, const integer_class &mod)
{
integer_class res = a % mod;
if (res < 0) {
res += mod;
}
return res;
}
void mp_pow_ui(rational_class &res, const rational_class &i, unsigned long n)
{
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));
}
void mp_powm(integer_class &res, const integer_class &base,
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 step(const unsigned long &n, const integer_class &i,
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);
}
bool positive_root(integer_class &res, const integer_class &i,
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;
}
// return true if i is a perfect nth power, i.e. res**i == n
bool mp_root(integer_class &res, const integer_class &i, const unsigned long n)
{
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;
}
integer_class mp_sqrt(const integer_class &i)
{
// 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;
}
void mp_rootrem(integer_class &a, integer_class &b, const integer_class &i,
unsigned long n)
{
mp_root(a, i, n);
integer_class p = pow(a, n);
;
b = i - p;
}
void mp_sqrtrem(integer_class &a, integer_class &b, const integer_class &i)
{
a = mp_sqrt(i);
b = i - boost::multiprecision::pow(a, 2);
}
// return nonzero if i is probably prime.
int mp_probab_prime_p(const integer_class &i, unsigned retries)
{
if (i % 2 == 0)
return (i == 2);
return miller_rabin_test(i, retries);
}
void mp_nextprime(integer_class &res, const integer_class &i)
{
// 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);
}
unsigned long mp_scan1(const integer_class &i)
{
if (i == 0) {
return ULONG_MAX;
}
return find_lsb(i, int_<0>());
}
// define simple 2x2 matrix with exponentiation by repeated squaring
// to use in logarithmic-time fibonacci calculation
struct two_by_two_matrix {
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 &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);
}
};
inline two_by_two_matrix fib_matrix(unsigned long n)
{
two_by_two_matrix x(1, 1, 1, 0);
return x.pow(n);
}
void mp_fib_ui(integer_class &res, unsigned long n)
{
// reference: https://www.nayuki.io/page/fast-fibonacci-algorithms
res = fib_matrix(n).data[0][1];
}
// sets a = Fibonacci(n) and b = Fibonacci(n-1)
void mp_fib2_ui(integer_class &a, integer_class &b, unsigned long n)
{
// 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];
}
inline two_by_two_matrix luc_matrix(unsigned long n)
{
two_by_two_matrix multiplier(1, 1, 1, 0);
two_by_two_matrix start(1, 0, 2, 0);
return multiplier.pow(n) * start;
}
void mp_lucnum_ui(integer_class &res, unsigned long n)
{
// 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];
}
void mp_lucnum2_ui(integer_class &a, integer_class &b, unsigned long n)
{
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];
}
void mp_fac_ui(integer_class &res, unsigned long n)
{
// 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;
}
}
void mp_bin_ui(integer_class &res, const integer_class &n, unsigned long r)
{
// slow, naive implementation
integer_class x = n - r;
res = 1;
for (unsigned long i = 1; i <= r; ++i) {
res *= x + i;
res /= i;
}
}
// this is extremely slow!
bool mp_perfect_power_p(const integer_class &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;
}
}
}
bool mp_perfect_square_p(const integer_class &i)
{
if (i < 0) {
return false;
}
integer_class root;
return mp_root(root, i, 2);
}
// according to the gmp documentation, the behavior of the
// corresponding function mpz_legendre is
// undefined if n is not a positive odd prime.
// hence we treat n as though it were a positive odd prime
// but it is up to the caller to very this.
int mp_legendre(const integer_class &a, const integer_class &n)
{
integer_class res;
mp_powm(res, a, integer_class((n - 1) / 2), n);
return res <= 1 ? res.convert_to<int>() : -1;
}
// private function that computes jacobi symbols
// without checking that arguments satisfy a >= 0
// and n is odd.
int unchecked_jacobi(const integer_class &a, const integer_class &n)
{
// 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);
}
// public interface for computing jacobi symbols. performs checking.
int mp_jacobi(const integer_class &a, const integer_class &n)
{
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);
}
int mp_kronecker(const integer_class &a, const integer_class &n)
{
/*
https://en.wikipedia.org/wiki/Kronecker_symbol
We compute the Kronecker symbol in terms of the Jacobi symbol.
For an integer n!=0, let n = u * PRODUCT (p_i**e_i), where u = -1 or 1 and
the p_i's and e_i's are the primes and corresponding powers
in the prime factorization of |n|.
Then for an integer a, the Kronecker symbol is given by
(a | n) = (a | u) * PRODUCT [(a | p_i)**e_i]
where
(a | u) is -1 if u is -1 and a < 0. Otherwise it is 1.
(a | 2) is 0 if a is even, 1 if (a mod 8) == 1 or 7, and -1 if (a mod 8) ==
3 or 5.
(a | p) is the Legendre symbol if p is an odd prime.
Let j be the power of the prime 2 in n's prime factorization, and let
m = |n|/(2**j) -- that is, m is |n| with all factors of 2 extracted.
Notice that, if p_1 is 2, then
PRODUCT [(a | p_i)**e_i]
== (a | 2)**j * PRODUCT [(a | p_i)**e_i] here the p_i's are the remaining
(odd) prime factors
== (a | 2)**j * (a | m), here (a | m) is the Jacobi
symbol, because m>0 is odd
Thus,
(a | n) == (a | u) * (a | 2)**j * (a | m) if n is even
(a | n) == (a | u) * (a | m) if n is odd
*/
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);
}
}
} // SymEngine
#endif // SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP