diff --git a/Makefile b/Makefile
index d64d919..8f10a71 100644
--- a/Makefile
+++ b/Makefile
@@ -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
diff --git a/algebra/Z.h b/algebra/Z.h
index 19b8a9a..7f5c058 100644
--- a/algebra/Z.h
+++ b/algebra/Z.h
@@ -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); }
 };
diff --git a/algebra/algebra.h b/algebra/algebra.h
index 7634b6e..42319f5 100644
--- a/algebra/algebra.h
+++ b/algebra/algebra.h
@@ -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>
diff --git a/algebra/laurentpoly.h b/algebra/laurentpoly.h
deleted file mode 100644
index ad2e1d9..0000000
--- a/algebra/laurentpoly.h
+++ /dev/null
@@ -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;
-}
diff --git a/algebra/module.h b/algebra/module.h
index edbff66..8cb71e7 100644
--- a/algebra/module.h
+++ b/algebra/module.h
@@ -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;
 }
diff --git a/algebra/multivariate_laurentpoly.h b/algebra/multivariate_laurentpoly.h
new file mode 100644
index 0000000..cf4df00
--- /dev/null
+++ b/algebra/multivariate_laurentpoly.h
@@ -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");
+}
diff --git a/algebra/multivariate_polynomial.h b/algebra/multivariate_polynomial.h
index 84c7d8f..75aac48 100644
--- a/algebra/multivariate_polynomial.h
+++ b/algebra/multivariate_polynomial.h
@@ -2,8 +2,6 @@
 /* multivariate polynomial in a (vector) variable x with coefficients
    in T. */
 
-enum variable { VARIABLE };
-
 template<unsigned n>
 class multivariate_monomial
 {
diff --git a/lib/io.cpp b/lib/io.cpp
index 273de3a..c900a73 100644
--- a/lib/io.cpp
+++ b/lib/io.cpp
@@ -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 ()
diff --git a/lib/io.h b/lib/io.h
index dfadb45..5ce6f9c 100644
--- a/lib/io.h
+++ b/lib/io.h
@@ -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 ();
   
diff --git a/lib/map_wrapper.h b/lib/map_wrapper.h
index 386f996..8aae708 100644
--- a/lib/map_wrapper.h
+++ b/lib/map_wrapper.h
@@ -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
 {
diff --git a/main.cpp b/main.cpp
index 5aaad92..dacb973 100644
--- a/main.cpp
+++ b/main.cpp
@@ -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  
 }
diff --git a/sseq.cpp b/sseq.cpp
index 0b76612..662c8ab 100644
--- a/sseq.cpp
+++ b/sseq.cpp
@@ -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;
diff --git a/sseq.h b/sseq.h
index 2bdb124..822ddbf 100644
--- a/sseq.h
+++ b/sseq.h
@@ -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);