politarczyk
/
knotkit
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
knotkit/algebra/multivariate_laurentpoly.h

575 lines
13 KiB
C++
Raw Permalink Blame History

#ifndef _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
#define _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
#include <sstream>
#include <string>
/* multivariate polynomial in a (vector) variable x with coefficients
in T. */
class multivariate_laurent_monomial
{
public:
map<unsigned, int> m;
explicit multivariate_laurent_monomial (const map<unsigned, int> &m_) : m(m_) { }
public:
multivariate_laurent_monomial () { }
multivariate_laurent_monomial (variable, unsigned j)
{
m.push (j, 1);
}
multivariate_laurent_monomial (variable, unsigned j, int e)
{
if (e != 0)
m.push (j, e);
}
multivariate_laurent_monomial (const multivariate_laurent_monomial &e) : m(e.m) { }
multivariate_laurent_monomial (copy, const multivariate_laurent_monomial &e)
: m(COPY, e.m)
{
}
multivariate_laurent_monomial (reader &r) : m(r) { }
~multivariate_laurent_monomial () { }
multivariate_laurent_monomial &operator = (const multivariate_laurent_monomial &e)
{
m = e.m;
return *this;
}
unsigned degree () const
{
int d = 0;
for (map<unsigned, int>::const_iter i = m; i; i ++)
d += i.val ();
return d;
}
bool operator == (const multivariate_laurent_monomial &e) const
{
#ifndef NDEBUG
check ();
e.check ();
#endif
return m == e.m;
}
bool operator != (const multivariate_laurent_monomial &e) const { return !operator == (e); }
bool operator == (int x) const
{
assert (x == 1);
#ifndef NDEBUG
check ();
#endif
return m.is_empty ();
}
bool operator != (int x) const { return !operator == (x); }
bool operator < (const multivariate_laurent_monomial &e) const
{
#ifndef NDEBUG
check ();
e.check ();
#endif
return m < e.m;
}
bool operator <= (const multivariate_laurent_monomial &e) const
{
#ifndef NDEBUG
check ();
e.check ();
#endif
return m <= e.m;
}
multivariate_laurent_monomial &operator *= (const multivariate_laurent_monomial &e)
{
for (map<unsigned, int>::const_iter i = e.m; i; i ++)
{
pair<int &, bool> p = m.find (i.key ());
if (p.second)
{
p.first += i.val ();
if (p.first == 0)
m -= i.key ();
}
else
p.first = i.val ();
}
#ifndef NDEBUG
check ();
#endif
return *this;
}
multivariate_laurent_monomial operator * (const multivariate_laurent_monomial &e) const
{
multivariate_laurent_monomial r;
r *= *this;
r *= e;
return r;
}
int operator () (unsigned j) const
{
return m(j, 0);
}
void set_exponent (unsigned j, int e)
{
if (e)
m[j] = e;
else
m -= j;
}
void push_exponent (unsigned j, int e)
{
assert (! (m % j));
if (e != 0)
m.push (j, e);
}
std::string to_string() const {
std::ostringstream res;
for (map<unsigned, int>::const_iter i = m; i; i ++) {
assert (i.val () != 0);
if (i.val () == 1) {
res << "x" << std::to_string(i.key());
}
else {
res << "x"
<< std::to_string(i.key())
<< "^"
<< std::to_string(i.val());
}
}
return res.str();
}
friend std::ostream& operator << (std::ostream& os, const multivariate_laurent_monomial& m) {
return os << m.to_string();
}
void show_self () const
{
std::cout << *this;
}
void write_self (writer &w) const { write (w, m); }
#ifndef NDEBUG
void check () const
{
for (map<unsigned, int>::const_iter i = m; i; i ++)
assert (i.val () != 0);
}
#endif
};
template<class T>
class multivariate_laurentpoly
{
public:
using linear_combination = ::linear_combination<multivariate_laurentpoly<T>>;
using linear_combination_const_iter = ::linear_combination_const_iter<multivariate_laurentpoly<T> >;
public:
using monomial = multivariate_laurent_monomial;
map<monomial, T> coeffs;
explicit multivariate_laurentpoly (const map<monomial, T> &coeffs_) : coeffs(coeffs_) { }
void normalize ()
{
for (typename map<monomial, T>::iter i = coeffs; i; i ++)
{
if (i.val () == 0)
i.del ();
}
#ifndef NDEBUG
check ();
#endif
}
public:
multivariate_laurentpoly () { }
multivariate_laurentpoly (int x)
{
T c (x);
if (c != 0)
coeffs.push (monomial (), c);
}
multivariate_laurentpoly (T c)
{
if (c != 0)
coeffs.push (monomial (), c);
}
multivariate_laurentpoly (T c, variable, unsigned i)
{
if (c != 0)
coeffs.push (monomial (VARIABLE, i), c);
}
multivariate_laurentpoly (T c, variable, unsigned i, int e)
{
if (c != 0)
coeffs.push (monomial (VARIABLE, i, e), c);
}
multivariate_laurentpoly (T c, const monomial &m)
{
if (c != 0)
coeffs.push (m, c);
}
template<class U>
multivariate_laurentpoly (const multivariate_laurentpoly<U>& pol);
multivariate_laurentpoly (const multivariate_laurentpoly &p) : coeffs(p.coeffs) { }
multivariate_laurentpoly (copy, const multivariate_laurentpoly &p)
// ??? COPY2?
: coeffs(COPY, p.coeffs)
{}
multivariate_laurentpoly (reader &r) : coeffs(r) { }
~multivariate_laurentpoly () { }
multivariate_laurentpoly &operator = (const multivariate_laurentpoly &p)
{
coeffs = p.coeffs;
return *this;
}
multivariate_laurentpoly &operator = (int x) { return operator = (T (x)); }
multivariate_laurentpoly &operator = (T c)
{
coeffs = map<monomial, T> ();
if (c != 0)
coeffs.push (monomial (), c);
return *this;
}
bool operator == (const multivariate_laurentpoly &p) const
{
#ifndef NDEBUG
check ();
p.check ();
#endif
return coeffs == p.coeffs;
}
bool operator != (const multivariate_laurentpoly &p) const { return !operator == (p); }
bool operator == (int x) const
{
#ifndef NDEBUG
check ();
#endif
T c (x);
if (c == 0)
return coeffs.is_empty ();
else
{
if (coeffs.card () != 1)
return 0;
pair<monomial, T> p = coeffs.head ();
return p.first == 1
&& p.second == c;
}
}
bool operator != (int x) const { return !operator == (x); }
bool operator < (const multivariate_laurentpoly &p) const
{
#ifndef NDEBUG
check ();
p.check ();
#endif
return coeffs < p.coeffs;
}
bool operator <= (const multivariate_laurentpoly &p) const
{
#ifndef NDEBUG
check ();
p.check ();
#endif
return coeffs <= p.coeffs;
}
unsigned card () const { return coeffs.card (); }
pair<monomial, T> head () const { return coeffs.head (); }
pair<monomial, T> tail() const { return coeffs.tail(); }
multivariate_laurentpoly &operator += (const multivariate_laurentpoly &p);
multivariate_laurentpoly &operator -= (const multivariate_laurentpoly &p);
multivariate_laurentpoly &operator *= (const multivariate_laurentpoly &p)
{
return operator = (*this * p);
}
multivariate_laurentpoly &operator *= (T s);
multivariate_laurentpoly &muladdeq (T c, variable, unsigned i)
{
monomial m (VARIABLE, i);
T &c2 = coeffs[m];
c2 += c;
if (c2 == 0)
coeffs -= m;
return *this;
}
multivariate_laurentpoly &muladdeq (T c, const monomial &m)
{
T &c2 = coeffs[m];
c2 += c;
if (c2 == 0)
coeffs -= m;
return *this;
}
multivariate_laurentpoly &muladdeq (multivariate_laurentpoly a, multivariate_laurentpoly b);
multivariate_laurentpoly operator - () const { return multivariate_laurentpoly () - *this; }
multivariate_laurentpoly operator + (const multivariate_laurentpoly &p) const
{
multivariate_laurentpoly r (COPY, *this);
r += p;
return r;
}
multivariate_laurentpoly operator - (const multivariate_laurentpoly &p) const
{
multivariate_laurentpoly r (COPY, *this);
r -= p;
return r;
}
multivariate_laurentpoly operator * (const multivariate_laurentpoly &p) const;
#ifndef NDEBUG
void check () const;
#endif
multivariate_laurentpoly& evaluate(T val, unsigned index);
static void show_ring ()
{
T::show_ring ();
printf ("[x_1^+/-1, ..., x_n^+/-1]");
}
std::string to_string() const;
void display_self () const {
std::cout << *this << "\n";
}
void show_self () const {
std::cout << *this;
}
void write_self (writer &w) const { write (w, coeffs); }
};
template<class T> template<class U>
multivariate_laurentpoly<T>::multivariate_laurentpoly(const multivariate_laurentpoly<U>& pol) {
for(typename map<monomial, U>::const_iter i = pol.coeffs; i; i++) {
if(T(i.val()) != 0)
coeffs.push(i.key(), T(i.val()));
}
}
template<class T> multivariate_laurentpoly<T>
operator * (const T &s, const multivariate_laurentpoly<T> &p)
{
multivariate_laurentpoly<T> r (COPY, p);
r *= s;
return r;
}
template<class T> multivariate_laurentpoly<T> &
multivariate_laurentpoly<T>::operator += (const multivariate_laurentpoly &p)
{
for (typename map<monomial, T>::const_iter i = p.coeffs; i; i ++)
{
monomial m = i.key ();
T &c = coeffs[m];
c += i.val ();
if (c == 0)
coeffs -= m;
}
return *this;
}
template<class T> multivariate_laurentpoly<T> &
multivariate_laurentpoly<T>::operator -= (const multivariate_laurentpoly &p)
{
for (typename map<monomial, T>::const_iter i = p.coeffs; i; i ++)
{
monomial m = i.key ();
T &c = coeffs[m];
c -= i.val ();
if (c == 0)
coeffs -= m;
}
return *this;
}
template<class T> multivariate_laurentpoly<T> &
multivariate_laurentpoly<T>::operator *= (T s)
{
if (s == 0)
coeffs.clear ();
else
{
for (typename map<monomial, T>::iter i = coeffs; i; i ++)
i.val () *= s;
normalize ();
}
return *this;
}
template<class T> multivariate_laurentpoly<T>
multivariate_laurentpoly<T>::operator * (const multivariate_laurentpoly &p) const
{
multivariate_laurentpoly r;
for (typename map<monomial, T>::const_iter i = coeffs; i; i ++)
for (typename map<monomial, T>::const_iter j = p.coeffs; j; j ++)
{
monomial m = i.key () * j.key ();
r.coeffs[m].muladdeq (i.val (), j.val ());
}
r.normalize ();
return r;
}
template<class T> multivariate_laurentpoly<T> &
multivariate_laurentpoly<T>::muladdeq (multivariate_laurentpoly a, multivariate_laurentpoly b)
{
for (typename map<monomial, T>::const_iter i = a.coeffs; i; i ++)
for (typename map<monomial, T>::const_iter j = b.coeffs; j; j ++)
{
monomial m = i.key () * j.key ();
coeffs[m].muladdeq (i.val (), j.val ());
}
normalize ();
return *this;
}
#ifndef NDEBUG
template<class T> void
multivariate_laurentpoly<T>::check () const
{
for (typename map<monomial, T>::const_iter i = coeffs; i; i ++)
assert (i.val () != 0);
}
#endif
template<class T>
std::string multivariate_laurentpoly<T>::to_string() const {
std::ostringstream res;
unsigned first = 1;
for (typename map<monomial, T>::const_iter i = coeffs; i; i ++) {
monomial m = i.key ();
T c = i.val ();
assert (c != 0);
if (first)
first = 0;
else
res << " + ";
if (m == 1) {
if (c == 1)
res << "1";
else
res << c;
}
else {
if (c != 1) {
res << c << "*";
}
res << m;
}
}
if (first)
res << "0";
return res.str();
}
template<class T>
std::ostream& operator << (std::ostream& os, const multivariate_laurentpoly<T>& pol) {
return os << pol.to_string();
}
// functions below were added to verify several periodicity criteria
// function below inverts every occurence of a variable x_index
// i.e. this means that x_index is replaced by x_index^{-1}
template<class T>
multivariate_laurentpoly<T> invert_variable(const multivariate_laurentpoly<T>& p, unsigned index = 1) {
multivariate_laurentpoly<T> result;
for (typename map<multivariate_laurent_monomial, T>::const_iter i = p.coeffs; i; i ++) {
multivariate_laurent_monomial mon = i.key();
T c = i.val();
multivariate_laurent_monomial mon2;
for(map<unsigned, int>::const_iter i = mon.m; i; i++) {
if(i.key() == index) {
mon2 *= multivariate_laurent_monomial(VARIABLE, i.key(), -i.val());
}
else {
mon2 *= multivariate_laurent_monomial(VARIABLE, i.key(), i.val());
}
}
result += multivariate_laurentpoly<T>(c, mon2);
}
return result;
}
template<class T>
multivariate_laurentpoly<T>
evaluate_with_copy(multivariate_laurentpoly<T> pol,
T val, unsigned index) {
using polynomial = multivariate_laurentpoly<T>;
using monomial = multivariate_laurent_monomial;
polynomial res;
for(typename map<monomial, T>::const_iter i = pol.coeffs; i; ++i) {
if(i.key().m % index) {
int exp = i.key().m[index];
monomial mon;
for(map<unsigned, int>::const_iter j = i.key().m; j; ++j) {
if(j.key() != index)
mon *= monomial(VARIABLE, j.key(), j.val());
}
res += polynomial(i.val() * pow(val, exp), mon);
}
else {
monomial mon(COPY, i.key());
res += polynomial(i.val(), mon);
}
}
return res;
}
template<class T>
multivariate_laurentpoly<T>&
multivariate_laurentpoly<T>::evaluate(T val, unsigned index) {
*this = evaluate_with_copy<T>(*this, val, index);
return *this;
}
#endif // _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
Reference in New Issue View Git Blame Copy Permalink