Replaced laurentpoly with multivariate_laurentpoly and updated uses

(sseq, module) appropriately; removed algebra/laurentpoly.h and added
algebra/multivariate_laurentpoly.h.  Added io, muladdeq to Z.  Added
functions open_file, close_file to lib/io.h.  Added operator < to
map_wrapper.
This commit is contained in:
Cotton Seed 2011-12-22 15:35:49 -05:00
parent 5d4acb7a86
commit 4443c8e331
13 changed files with 594 additions and 440 deletions

View File

@ -34,6 +34,7 @@ ALGEBRA_HEADERS = algebra/algebra.h algebra/grading.h algebra/module.h \
algebra/Z2.h algebra/linear_combination.h \
algebra/Z.h algebra/Zp.h algebra/Q.h \
algebra/polynomial.h algebra/multivariate_polynomial.h \
algebra/multivariate_laurentpoly.h \
algebra/laurentpoly.h algebra/fraction_field.h
KNOTKIT_HEADERS = knotkit.h planar_diagram.h dt_code.h knot_diagram.h \
smoothing.h cobordism.h cube.h spanning_tree_complex.h cube_impl.h sseq.h

View File

@ -18,7 +18,18 @@ class Z
Z_impl (int init) { mpz_init_set_si (x, init); }
Z_impl (copy, mpz_srcptr init) { mpz_init_set (x, init); }
Z_impl (steal, mpz_srcptr init) { x[0] = *init; }
Z_impl (reader &r)
{
mpz_init (x);
mpz_inp_raw (x, r.fp);
}
~Z_impl () { mpz_clear (x); }
void write_self (writer &w) const
{
mpz_out_raw (w.fp, x);
}
};
ptr<Z_impl> impl;
@ -30,6 +41,7 @@ class Z
Z (int init) : impl(new Z_impl (init)) { }
Z (const Z &z) : impl(z.impl) { }
Z (copy, const Z &z) : impl(new Z_impl (COPY, z.impl->x)) { }
Z (reader &r) : impl(new Z_impl (r)) { }
~Z () { }
Z &operator = (const Z &z) { impl = z.impl; return *this; }
@ -100,6 +112,12 @@ class Z
return Z (COPY, *this);
}
Z &muladdeq (const Z &z1, const Z &z2)
{
// ??? use muladd primitive
return operator += (z1 * z2);
}
Z &operator += (const Z &z)
{
mpz_add (impl->x, impl->x, z.impl->x);
@ -176,4 +194,5 @@ class Z
static void show_ring () { printf ("Z"); }
void show_self () const { mpz_out_str (stdout, 10, impl->x); }
void display_self () const { show_self (); newline (); }
void write_self (writer &w) const { write (w, *impl); }
};

View File

@ -51,15 +51,16 @@ template<class R> class linear_combination_const_iter;
template<class R> class module;
#include <algebra/multivariate_polynomial.h>
#include <algebra/laurentpoly.h>
/* constructor tag for multivariate_polynomial,
multivariate_laurentpoly. */
enum variable { VARIABLE };
#include <algebra/grading.h>
#include <algebra/module.h>
#include <algebra/Z2.h>
#include <algebra/linear_combination.h>
#include <algebra/multivariate_polynomial.h>
#include <algebra/multivariate_laurentpoly.h>
#include <algebra/Z.h>
#include <algebra/Zp.h>
@ -67,3 +68,6 @@ template<class R> class module;
#include <algebra/polynomial.h>
#include <algebra/fraction_field.h>
#include <algebra/module.h>
#include <algebra/linear_combination.h>

View File

@ -1,399 +0,0 @@
template<class C, class E> class laurentpoly;
class exponent1
{
public:
int first;
public:
exponent1 () : first(0) { }
explicit exponent1 (unsigned first_) : first(first_) { }
exponent1 (const exponent1 &e) : first(e.first) { }
~exponent1 () { }
exponent1 &operator = (const exponent1 &e) { first = e.first; return *this; }
inline laurentpoly<int, exponent1> operator + (const exponent1 &e) const;
exponent1 &operator *= (const exponent1 &e) { first += e.first; return *this; }
exponent1 operator * (const exponent1 &e) const { return exponent1 (first + e.first); }
template<class T> laurentpoly<T, exponent1> operator * (const T &c) const
{
laurentpoly<T, exponent1> r;
r.addeq (c, *this);
return r;
}
bool operator == (int x) const { assert (x == 0); return !first; }
bool operator != (int x) const { return !operator == (x); }
bool operator == (const exponent1 &e) const { return first == e.first; }
bool operator != (const exponent1 &e) const { return !operator == (e); }
bool operator < (const exponent1 &e) const { return first < e.first; }
exponent1 operator ^ (int p) const { return exponent1 (first * p); }
void show_self () const { printf ("(%d)", first); }
void display_self () const { show_self (); newline (); }
};
class exponent2
{
public:
int first,
second;
public:
exponent2 () : first(0), second(0) { }
explicit exponent2 (int i) : first(0), second(0) { assert (i == 0); }
exponent2 (int first_, int second_) : first(first_), second(second_) { }
exponent2 (const exponent2 &e) : first(e.first), second(e.second) { }
~exponent2 () { }
exponent2 &operator = (const exponent2 &e) { first = e.first; second = e.second; return *this; }
inline laurentpoly<int, exponent2> operator + (const exponent2 &e) const;
exponent2 &operator *= (const exponent2 &e) { first += e.first; second += e.second; return *this; }
exponent2 operator * (const exponent2 &e) const { return exponent2 (first + e.first, second + e.second); }
template<class T> laurentpoly<T, exponent2> operator * (const T &c) const
{
laurentpoly<T, exponent2> r;
r.addeq (c, *this);
return r;
}
bool operator == (int x) const { assert (x == 0); return !first && !second; }
bool operator != (int x) const { return !operator == (x); }
bool operator == (const exponent2 &e) const { return first == e.first && second == e.second; }
bool operator != (const exponent2 &e) const { return !operator == (e); }
bool operator < (const exponent2 &e) const { return first < e.first
|| (first == e.first && second < e.second); }
exponent2 operator ^ (int p) const { return exponent2 (first * p, second * p); }
void show_self () const { printf ("(%d,%d)", first, second); }
void display_self () const { show_self (); newline (); }
unsigned texwidth () const
{
if (!first && !second)
return 1;
unsigned w = 0;
if (first == 1)
w ++;
else if (first < 0)
w += 3;
else if (first)
w += 2;
if (second == 1)
w ++;
else if (second < 0)
w += 3;
else if (second)
w += 2;
return w;
}
void texshow (FILE *fp) const
{
if (!first && !second)
printf ("1");
else
{
if (first == 1)
printf ("t");
else if (first)
printf ("t^{%d}", first);
if (second == 1)
printf ("q");
else if (second)
printf ("q^{%d}", second);
}
}
};
#if 0
template<class T, class E> laurentpoly<T, E> operator * (const T &c, const E &e)
{
return e*c;
}
#endif
template<class C, class E>
class laurentpoly
{
public:
typedef typename map<E, C>::iter map_iter;
typedef typename map<E, C>::const_iter map_const_iter;
map<E, C> m;
public:
laurentpoly () { }
explicit laurentpoly (const C &c) { m.push (E (0), c); }
explicit laurentpoly (const E &e) { m.push (e, 1); }
laurentpoly (const laurentpoly &p) : m(p.m) { }
laurentpoly (copy, const laurentpoly &p) : m(COPY, p.m) { }
~laurentpoly () { }
laurentpoly &operator = (const laurentpoly &p) { m = p.m; return *this; }
laurentpoly &operator = (const C &c) { m.clear (); m.push (E (0), c); }
bool operator == (const laurentpoly &p) { return m == p.m; }
bool operator != (const laurentpoly &p) { return m != p.m; }
laurentpoly &operator += (const laurentpoly &p);
laurentpoly &operator -= (const laurentpoly &p);
laurentpoly &operator *= (const laurentpoly &p) { operator = (*this * p); }
laurentpoly &operator *= (const C &s);
laurentpoly &subeq (const C &s, const E &e);
laurentpoly &addeq (const C &s, const E &e);
laurentpoly &muleq (const C &s, const E &e);
laurentpoly &muladdeq (const C &s, const E &e, const laurentpoly &p);
laurentpoly operator - () const;
laurentpoly operator + (const E &e) const { return *this + laurentpoly (e); }
laurentpoly operator + (const C &c) const { return *this + laurentpoly (c); }
laurentpoly operator + (const laurentpoly &p) const;
laurentpoly operator - (const laurentpoly &p) const;
laurentpoly operator * (const laurentpoly &p) const;
laurentpoly operator * (const E &e) const;
laurentpoly operator * (const C &s) const;
C aug () const
{
C r = 0;
for (map_const_iter i = m; i; i ++)
r += i.val ();
return r;
}
void show_self () const;
void display_self () const { show_self (); newline (); }
void texshow (FILE *fp);
};
template<class C, class E> inline laurentpoly<C, E> operator + (const E &e, const laurentpoly<C, E> &p) { return p + e; }
template<class C, class E> inline laurentpoly<C, E> operator + (const C &c, const laurentpoly<C, E> &p) { return p + c; }
template<class C, class E> inline laurentpoly<C, E> operator * (const C &s, const laurentpoly<C, E> &p) { return p * s; }
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::operator += (const laurentpoly &p)
{
for (map_const_iter i = p.m; i; i ++)
addeq (i.val (), i.key ());
return *this;
}
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::operator -= (const laurentpoly &p)
{
for (map_const_iter i = p.m; i; i ++)
subeq (i.val (), i.key ());
return *this;
}
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::operator *= (const C &s)
{
for (map_iter i = m; i; i ++)
i.val () *= s;
return *this;
}
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::addeq (const C &s, const E &e)
{
pair<C &, bool> x = m.find (e);
if (!x.second)
x.first = s;
else
x.first += s;
return *this;
}
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::subeq (const C &s, const E &e)
{
pair<C &, bool> x = m.find (e);
if (!x.second)
x.first = -s;
else
x.first -= s;
return *this;
}
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::muleq (const C &s, const E &e)
{
map<C, E> new_m;
for (map_const_iter i = m; i; i ++)
m.push (e * m.key (), s * m.val ());
m = new_m;
return *this;
}
template<class C, class E> laurentpoly<C, E> &
laurentpoly<C, E>::muladdeq (const C &s, const E &e, const laurentpoly &p)
{
for (map_const_iter i = p.m; i; i ++)
addeq (s * i.val (), e * i.key ());
}
template<class C, class E> laurentpoly<C, E>
laurentpoly<C, E>::operator - () const
{
map<C, E> new_m;
for (map_const_iter i = m; i; i ++)
new_m.push (m.key (), -m.val ());
return laurentpoly (new_m);
}
template<class C, class E> laurentpoly<C, E>
laurentpoly<C, E>::operator + (const laurentpoly &p) const
{
laurentpoly r (COPY, *this);
r += p;
return r;
}
template<class C, class E> laurentpoly<C, E>
laurentpoly<C, E>::operator - (const laurentpoly &p) const
{
laurentpoly r (COPY, *this);
r -= p;
return r;
}
template<class C, class E> laurentpoly<C, E>
laurentpoly<C, E>::operator * (const laurentpoly &p) const
{
laurentpoly r (*this);
for (map_const_iter i = m; i; i ++)
r.muladdeq (i.val (), i.key (), p.m);
return r;
}
template<class C, class E> laurentpoly<C, E>
laurentpoly<C, E>::operator * (const C &s) const
{
laurentpoly r (COPY, *this);
for (map_iter i = r.m; i; i ++)
i.val () *= s;
return r;
}
template<class C, class E> laurentpoly<C, E>
laurentpoly<C, E>::operator * (const E &e) const
{
laurentpoly r;
for (map_const_iter i = m; i; i ++)
r.addeq (i.val (), i.key () * e);
return r;
}
template<class C, class E> void
laurentpoly<C, E>::show_self () const
{
if (m.card () == 0)
{
printf ("0");
return;
}
bool first = 1;
for (map_const_iter i = m; i; i ++)
{
if (i.val () != 0)
{
if (!first)
printf (" + ");
if (i.val () == 1)
show (i.key ());
else
{
show (i.val ());
printf ("*");
show (i.key ());
}
first = 0;
}
}
}
template<class C, class E> void
laurentpoly<C, E>::texshow (FILE *fp)
{
bool first = 1;
unsigned w = 0;
for (map_const_iter i = m; i; i ++)
{
if (i.val () == 0)
continue;
if (w > 47)
{
fprintf (fp, " \\\\\n");
fprintf (fp, " & & & ");
w = 0;
}
if (first)
first = 0;
else
{
printf ("+");
w += 3;
}
if (i.val () == 1)
{
i.key ().texshow (fp);
w += i.key ().texwidth ();
}
else
{
fprintf (fp, "%d", i.val ());
w ++;
if (i.key () != 0)
{
i.key ().texshow (fp);
w += i.key ().texwidth ();
}
}
}
}
typedef laurentpoly<int, exponent1> intpoly1;
typedef laurentpoly<int, exponent2> intpoly2;
inline laurentpoly<int, exponent1>
exponent1::operator + (const exponent1 &e) const
{
intpoly1 r;
r.addeq (1, *this);
r.addeq (1, e);
return r;
}
inline laurentpoly<int, exponent2>
exponent2::operator + (const exponent2 &e) const
{
intpoly2 r;
r.addeq (1, *this);
r.addeq (1, e);
return r;
}

View File

@ -56,8 +56,8 @@ class module : public refcounted
ptr<const free_submodule<R> > submodule (const Rmod_span &span) const;
intpoly2 free_poincare_polynomial () const;
intpoly1 free_delta_poincare_polynomial () const;
multivariate_laurentpoly<Z> free_poincare_polynomial () const;
multivariate_laurentpoly<Z> free_delta_poincare_polynomial () const;
void show_self () const;
void display_self () const;
@ -779,26 +779,31 @@ module<R>::submodule (const Rmod_span &span) const
span.pivots);
}
template<class R> intpoly2
template<class R> multivariate_laurentpoly<Z>
module<R>::free_poincare_polynomial () const
{
intpoly2 r;
multivariate_laurentpoly<Z> r;
for (unsigned i = 1; i <= free_rank (); i ++)
{
grading hq = generator_grading (i);
r.addeq (1, exponent2 (hq.h, hq.q));
multivariate_laurent_monomial m;
m.push_exponent (1, hq.h);
m.push_exponent (2, hq.q);
r.muladdeq (1, m);
}
return r;
}
template<class R> intpoly1
template<class R> multivariate_laurentpoly<Z>
module<R>::free_delta_poincare_polynomial () const
{
intpoly1 r;
multivariate_laurentpoly<Z> r;
for (unsigned i = 1; i <= free_rank (); i ++)
{
grading hq = generator_grading (i);
r.addeq (1, exponent1 (hq.q - 2*hq.h));
multivariate_laurent_monomial m;
m.push_exponent (1, hq.q - 2*hq.h);
r.muladdeq (1, m);
}
return r;
}

View File

@ -0,0 +1,447 @@
/* 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
{
#ifndef NDEBUG
check ();
e.check ();
#endif
return m < e.m;
}
bool operator == (int x) const
{
assert (x == 1);
#ifndef NDEBUG
check ();
#endif
return m.is_empty ();
}
bool operator != (int x) const { return !operator == (x); }
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);
}
void show_self () const
{
for (map<unsigned, int>::const_iter i = m; i; i ++)
{
assert (i.val () != 0);
if (i.val () == 1)
printf ("x%d", i.key ());
else
printf ("x%d^%d", i.key (), i.val ());
}
}
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:
typedef ::linear_combination<multivariate_laurentpoly<T> > linear_combination;
typedef ::linear_combination_const_iter<multivariate_laurentpoly<T> >
linear_combination_const_iter;
public:
typedef multivariate_laurent_monomial 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, const monomial &m)
{
if (c != 0)
coeffs.push (m, c);
}
multivariate_laurentpoly (const multivariate_laurentpoly &p) : coeffs(p.coeffs) { }
multivariate_laurentpoly (copy, const multivariate_laurentpoly &p)
: 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 == (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 (); }
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); }
multivariate_laurentpoly &operator += (multivariate_laurentpoly p);
multivariate_laurentpoly &operator -= (multivariate_laurentpoly p);
multivariate_laurentpoly &operator *= (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 + (multivariate_laurentpoly p) const
{
multivariate_laurentpoly r (COPY, *this);
r += p;
return r;
}
multivariate_laurentpoly operator - (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
static void show_ring ()
{
T::show_ring ();
printf ("[x_1^+/-1, ..., x_n^+/-1]");
}
void display_self () const { show_self (); newline (); }
void show_self () const;
void write_self (writer &w) const { write (w, coeffs); }
};
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 += (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 -= (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> void
multivariate_laurentpoly<T>::show_self () const
{
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
printf (" + ");
if (m == 1)
{
if (c == 1)
printf ("1");
else
show (c);
}
else
{
if (c != 1)
{
show (c);
printf ("*");
}
show (m);
}
}
if (first)
printf ("0");
}

View File

@ -2,8 +2,6 @@
/* multivariate polynomial in a (vector) variable x with coefficients
in T. */
enum variable { VARIABLE };
template<unsigned n>
class multivariate_monomial
{

View File

@ -1,15 +1,27 @@
#include <lib/lib.h>
writer::writer (const std::string &filename)
: fp(0)
FILE *open_file (const std::string &file, const char *mode)
{
fp = fopen (filename.c_str (), "w");
FILE *fp = fopen (file.c_str (), mode);
if (fp == 0)
{
stderror ("fopen: %s", filename.c_str ());
stderror ("fopen: %s", file.c_str ());
exit (EXIT_FAILURE);
}
return fp;
}
void close_file (FILE *fp)
{
fclose (fp);
}
writer::writer (const std::string &file)
: fp(0)
{
fp = open_file (file, "w");
}
writer::~writer ()
@ -21,15 +33,10 @@ writer::~writer ()
}
}
reader::reader (const std::string &filename)
reader::reader (const std::string &file)
: fp(0)
{
fp = fopen (filename.c_str (), "r");
if (fp == 0)
{
stderror ("fopen: %s", filename.c_str ());
exit (EXIT_FAILURE);
}
fp = open_file (file, "r");
}
reader::~reader ()

View File

@ -1,11 +1,14 @@
extern FILE *open_file (const std::string &file, const char *mode);
extern void close_file (FILE *fp);
class writer
{
public:
FILE *fp;
public:
writer (const std::string &filename);
writer (const std::string &file);
writer (const writer &); // doesn't exist
~writer ();
@ -24,7 +27,7 @@ class reader
FILE *fp;
public:
reader (const std::string &filename);
reader (const std::string &file);
reader (const reader &); // doesn't exist
~reader ();

View File

@ -129,6 +129,7 @@ class map_wrapper
}
bool operator == (const map_wrapper &m) const;
bool operator < (const map_wrapper &m) const;
void write_self (writer &w) const;
};
@ -151,20 +152,51 @@ map_wrapper<M, K, V>::map_wrapper (reader &r)
template<class M, class K, class V> bool
map_wrapper<M, K, V>::operator == (const map_wrapper &m) const
{
// event bettr: use (ordered) dual iterators
if (card () != m.card ())
return 0;
for (const_iter i = *this; i; i ++)
const_iter i = *this,
j = m;
while (i && j)
{
// ??? inefficient: use operator ^
if (! (m % i.key ())
|| m(i.key ()) != i.val ())
if (i.key () != j.key ()
|| i.val () != j.val ())
return 0;
i ++;
j ++;
}
if (i || j)
return 0;
return 1;
}
template<class M, class K, class V> bool
map_wrapper<M, K, V>::operator < (const map_wrapper &m) const
{
const_iter i = *this,
j = m;
while (i && j)
{
if (i.key () < j.key ())
return 1;
else if (i.key () > j.key ())
return 0;
else
{
assert (i.key () == j.key ());
if (i.val () < j.val ())
return 1;
else if (i.val () > j.val ())
return 0;
else
assert (i.val () == j.val ());
}
i ++;
j ++;
}
return !i && j;
}
template<class M, class K, class V> void
map_wrapper<M, K, V>::write_self (writer &w) const
{

View File

@ -111,6 +111,37 @@ test_field ()
int
main ()
{
knot_diagram kd (rolfsen_knot (8, 19));
cube<Z2> c (kd);
sseq ss = compute_szabo_sseq (c);
multivariate_laurentpoly<Z> ssp = ss.pages[1].poincare_polynomial (ss.bounds);
display ("ssp: ", ssp);
multivariate_laurentpoly<Z> p;
p.muladdeq (5, multivariate_laurent_monomial (VARIABLE, 1, -2));
p.muladdeq (-6, multivariate_laurent_monomial (VARIABLE, 2, 13));
p.muladdeq (14, (multivariate_laurent_monomial (VARIABLE, 1, 5)
* multivariate_laurent_monomial (VARIABLE, 2, -6)));
display ("p: ", p);
display ("p*p: ", p*p);
{
writer w ("dump");
write (w, p*p);
}
{
reader r ("dump");
multivariate_laurentpoly<Z> q (r);
display ("q: ", q);
assert (q == p*p);
}
#if 0
test_ring<Z2> (2);
test_ring<Z> (0);
test_ring<Zp<2> > (2);
@ -123,4 +154,5 @@ main ()
test_field<Zp<3> > ();
test_field<Z2> ();
test_field<Zp<2> > ();
#endif
}

View File

@ -103,33 +103,38 @@ sseq_page::homological_width (const sseq_bounds &b) const
}
}
intpoly2
multivariate_laurentpoly<Z>
sseq_page::poincare_polynomial (const sseq_bounds &b) const
{
intpoly2 r;
multivariate_laurentpoly<Z> r;
for (int i = b.minh; i <= b.maxh; i ++)
{
for (int j = b.minq; j <= b.maxq; j ++)
{
unsigned c = rank[i - b.minh][j - b.minq];
if (c)
r.addeq (c, exponent2 (i, j));
{
multivariate_laurent_monomial m;
m.push_exponent (1, i);
m.push_exponent (2, j);
r.muladdeq (c, m);
}
}
}
return r;
}
intpoly1
multivariate_laurentpoly<Z>
sseq_page::delta_poincare_polynomial (const sseq_bounds &b) const
{
intpoly1 r;
multivariate_laurentpoly<Z> r;
for (int i = b.minh; i <= b.maxh; i ++)
{
for (int j = b.minq; j <= b.maxq; j ++)
{
unsigned c = rank[i - b.minh][j - b.minq];
if (c)
r.addeq (c, exponent1 (j - 2*i));
r.muladdeq (c, multivariate_laurent_monomial (VARIABLE, 1, j - 2*i));
}
}
return r;

4
sseq.h
View File

@ -55,8 +55,8 @@ class sseq_page
bool equal_as_spaces (const sseq_page &pg) const { return rank == pg.rank; }
unsigned total_rank () const;
intpoly2 poincare_polynomial (const sseq_bounds &b) const;
intpoly1 delta_poincare_polynomial (const sseq_bounds &b) const;
multivariate_laurentpoly<Z> poincare_polynomial (const sseq_bounds &b) const;
multivariate_laurentpoly<Z> delta_poincare_polynomial (const sseq_bounds &b) const;
unsigned homological_width (const sseq_bounds &b) const;
void addeq (const sseq_bounds &b, const sseq_bounds &b2, const sseq_page &pg2);