diff --git a/Makefile b/Makefile
index f207be9..0558a56 100644
--- a/Makefile
+++ b/Makefile
@@ -2,12 +2,13 @@
 BISON = /opt/local/bin/bison
 FLEX = /opt/local/bin/flex
 
-CXX = g++
+# CXX = g++
+CXX = clang++ -fno-color-diagnostics --stdlib=libc++ --std=c++11
 
 INCLUDES = -I/opt/local/include -I.
 
-# OPTFLAGS = -g
-OPTFLAGS = -O2 -g
+OPTFLAGS = -g
+# OPTFLAGS = -O2 -g
 # OPTFLAGS = -O2 -DNDEBUG
 
 LDFLAGS = -L/opt/local/lib
@@ -38,7 +39,7 @@ ALGEBRA_HEADERS = algebra/algebra.h algebra/grading.h algebra/module.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
 
-LIBS = -lgmp
+LIBS = -lgmp -lz
 
 all: gss
 
diff --git a/algebra/Q.h b/algebra/Q.h
index b6ee350..21a1d06 100644
--- a/algebra/Q.h
+++ b/algebra/Q.h
@@ -31,16 +31,16 @@ class Q
     Q_impl (reader &r)
     {
       mpq_init (x);
-      mpz_inp_raw (mpq_numref (x), r.fp);
-      mpz_inp_raw (mpq_denref (x), r.fp);
+      r.read_mpz (mpq_numref (x));
+      r.read_mpz (mpq_denref (x));
     }
     
     ~Q_impl () { mpq_clear (x); }
     
     void write_self (writer &w) const
     {
-      mpz_out_raw (w.fp, mpq_numref (x));
-      mpz_out_raw (w.fp, mpq_denref (x));
+      w.write_mpz (mpq_numref (x));
+      w.write_mpz (mpq_denref (x));
     }
   };
   
@@ -163,12 +163,12 @@ class Q
   
   Q div (const Q &d) const { return operator / (d); }
   
-  triple<Q, Q, Q> extended_gcd (const Q &q) const
+  tuple<Q, Q, Q> extended_gcd (const Q &q) const
   {
     if (*this != 0)
-      return triple<Q, Q, Q> (*this, 1, 0);
+      return tuple<Q, Q, Q> (*this, 1, 0);
     else
-      return triple<Q, Q, Q> (q, 0, 1);
+      return tuple<Q, Q, Q> (q, 0, 1);
   }
   
   Q gcd (const Q &q) const
diff --git a/algebra/Z.h b/algebra/Z.h
index 3575992..d23c0fc 100644
--- a/algebra/Z.h
+++ b/algebra/Z.h
@@ -21,14 +21,14 @@ class Z
     Z_impl (reader &r)
     {
       mpz_init (x);
-      mpz_inp_raw (x, r.fp);
+      r.read_mpz (x);
     }
     
     ~Z_impl () { mpz_clear (x); }
     
     void write_self (writer &w) const
     {
-      mpz_out_raw (w.fp, x);
+      w.write_mpz (x);
     }
   };
   
@@ -48,6 +48,7 @@ class Z
   Z &operator = (int x) { impl = new Z_impl (x); return *this; }
   
   bool operator == (const Z &z) const { return mpz_cmp (impl->x, z.impl->x) == 0; }
+  bool operator != (const Z &z) const { return !operator == (z); }
   
   bool operator == (int y) const { return mpz_cmp_si (impl->x, y) == 0; }
   bool operator != (int y) const { return !operator == (y); }
@@ -91,13 +92,13 @@ class Z
     return Z (STEAL, x);
   }
   
-  Z operator / (const Z &z) const
+  Z operator / (const Z &denom) const
   {
-    if (mpz_cmp_si (z.impl->x, 1) == 0)
+    if (mpz_cmp_si (denom.impl->x, 1) == 0)
       return *this;
     else
       {
-	assert (mpz_cmp_si (z.impl->x, -1) == 0);
+	assert (mpz_cmp_si (denom.impl->x, -1) == 0);
 	
 	mpz_t x;
 	mpz_init (x);
@@ -112,6 +113,7 @@ class Z
     return Z (COPY, *this);
   }
   
+  // *this += z1*z2
   Z &muladdeq (const Z &z1, const Z &z2)
   {
     mpz_addmul (impl->x, z1.impl->x, z2.impl->x);
@@ -148,33 +150,30 @@ class Z
     return *this;
   }
   
-  bool divides (const Z &n) const
+  bool divides (const Z &num) const
   {
-    // d = *this
-    return mpz_divisible_p (n.impl->x, impl->x);
+    return mpz_divisible_p (num.impl->x, impl->x);
   }
+  bool operator | (const Z &num) const { return divides (num); }
   
-  bool operator | (const Z &z) const { return divides (z); }
-  
-  Z div (const Z &d) const
+  Z divide_exact (const Z &denom) const
   {
-    // n = *this
+    // num = *this
     mpz_t q;
     mpz_init (q);
-    mpz_divexact (q, impl->x, d.impl->x);
+    mpz_divexact (q, impl->x, denom.impl->x);
     return Z (STEAL, q);
   }
   
-  triple<Z, Z, Z> extended_gcd (const Z &z) const
+  tuple<Z, Z> divide_with_remainder (const Z &denom) const
   {
-    mpz_t d, s, t;
-    mpz_init (d);
-    mpz_init (s);
-    mpz_init (t);
-    mpz_gcdext (d, s, t, impl->x, z.impl->x);
-    return triple<Z, Z, Z> (Z (STEAL, d),
-			    Z (STEAL, s),
-			    Z (STEAL, t));
+    // *this = num
+    mpz_t q, r;
+    mpz_init (q);
+    mpz_init (r);
+    mpz_tdiv_qr (q, r, impl->x, denom.impl->x);
+    return make_tuple (Z (STEAL, q),
+		       Z (STEAL, r));
   }
   
   Z gcd (const Z &z) const
@@ -191,6 +190,18 @@ class Z
     return Z (STEAL, m);
   }
   
+  tuple<Z, Z, Z> extended_gcd (const Z &z) const
+  {
+    mpz_t d, s, t;
+    mpz_init (d);
+    mpz_init (s);
+    mpz_init (t);
+    mpz_gcdext (d, s, t, impl->x, z.impl->x);
+    return make_tuple (Z (STEAL, d),
+		       Z (STEAL, s),
+		       Z (STEAL, t));
+  }
+  
   static void show_ring () { printf ("Z"); }
   void show_self () const { mpz_out_str (stdout, 10, impl->x); }
   void display_self () const { show_self (); newline (); }
diff --git a/algebra/Z2.h b/algebra/Z2.h
index f13f192..0cb8799 100644
--- a/algebra/Z2.h
+++ b/algebra/Z2.h
@@ -54,12 +54,12 @@ class Z2
     return Z2 (1);
   }
   
-  triple<Z2, Z2, Z2> extended_gcd (Z2 x) const
+  tuple<Z2, Z2, Z2> extended_gcd (Z2 x) const
   {
     if (v)
-      return triple<Z2, Z2, Z2> (Z2 (1), Z2 (1), Z2 (0));
+      return make_tuple (Z2 (1), Z2 (1), Z2 (0));
     else
-      return triple<Z2, Z2, Z2> (Z2 (1), Z2 (0), Z2 (1));
+      return make_tuple (Z2 (1), Z2 (0), Z2 (1));
   }
   
   static void show_ring () { printf ("Z2"); }
diff --git a/algebra/Zp.h b/algebra/Zp.h
index dc53609..bd3fe0a 100644
--- a/algebra/Zp.h
+++ b/algebra/Zp.h
@@ -32,6 +32,7 @@ class Zp
   }
   
   bool operator == (const Zp &x) const { return v == x.v; }
+  bool operator != (const Zp &x) const { return !operator == (x); }
   
   bool operator == (int x) const { return operator == (Zp (x)); }
   bool operator != (int x) const { return !operator == (x); }
@@ -56,11 +57,11 @@ class Zp
   
   Zp recip () const
   {
-    triple<unsigned, int, int> t = unsigned_extended_gcd (v, p);
-    assert (t.first == 1);
-    assert ((int)t.first == t.second*(int)v  + t.third*(int)p);
+    tuple<unsigned, int, int> t = unsigned_extended_gcd (v, p);
+    assert (get<0> (t) == 1);
+    assert ((int)get<0> (t) == get<1> (t)*(int)v  + get<2> (t)*(int)p);
     
-    return Zp (t.second);
+    return Zp (get<1> (t));
   }
   
   Zp &operator += (const Zp &x)
@@ -92,12 +93,12 @@ class Zp
   
   Zp div (const Zp &d) const { return operator / (d); }
   
-  triple<Zp, Zp, Zp> extended_gcd (const Zp &x) const
+  tuple<Zp, Zp, Zp> extended_gcd (const Zp &x) const
   {
     if (v)
-      return triple<Zp, Zp, Zp> (v, 1, 0);
+      return make_tuple (v, Zp (1), Zp (0));
     else
-      return triple<Zp, Zp, Zp> (x, 0, 1);
+      return make_tuple (x, Zp (0), Zp (1));
   }
   
   Zp gcd (const Zp &x) const
diff --git a/algebra/algebra.cpp b/algebra/algebra.cpp
index e467674..451f4cb 100644
--- a/algebra/algebra.cpp
+++ b/algebra/algebra.cpp
@@ -39,43 +39,43 @@ int64_gcd (int64 a, int64 b)
 		     (uint64)std::abs (b));
 }
 
-static triple<unsigned, int, int>
+static tuple<unsigned, int, int>
 extended_gcd_1 (int a, int b)
 {
   if (b == 0)
-    return triple<unsigned, int, int> (a, 1, 0);
+    return tuple<unsigned, int, int> (a, 1, 0);
   
   unsigned t = a % b;
   if (!t)
-    return triple<unsigned, int, int> (b, 0, 1);
+    return make_tuple (b, 0, 1);
   else
     {
-      triple<unsigned, int, int> s = unsigned_extended_gcd (b, t);
+      tuple<unsigned, int, int> s = unsigned_extended_gcd (b, t);
       
-      unsigned d = s.first;
-      int x = s.third,
-	y = s.second - s.third * (a / b);
+      unsigned d = get<0> (s);
+      int x = get<2> (s),
+	y = get<1> (s) - get<2> (s) * (a / b);
       
       assert ((int)d == a*x + b*y);
       
-      return triple<unsigned, int, int> (d, x, y);
+      return make_tuple (d, x, y);
     }
 }
 
-triple<unsigned, int, int>
+tuple<unsigned, int, int>
 unsigned_extended_gcd (unsigned a, unsigned b)
 {
   return extended_gcd_1 ((int)a, (int)b);
 }
 
-triple<unsigned, int, int>
+tuple<unsigned, int, int>
 extended_gcd (int a, int b)
 {
-  triple<unsigned, int, int> t = extended_gcd_1 (std::abs (a),
-						 std::abs (b));
-  unsigned d = t.first;
-  int x = t.second,
-    y = t.third;
+  tuple<unsigned, int, int> t = extended_gcd_1 (std::abs (a),
+						std::abs (b));
+  unsigned d = get<0> (t);
+  int x = get<1> (t),
+    y = get<2> (t);
   if (a < 0)
     x *= -1;
   if (b < 0)
@@ -83,7 +83,7 @@ extended_gcd (int a, int b)
   
   assert ((int)d == a*x + b*y);
   
-  return triple<unsigned, int, int> (d, x, y);
+  return make_tuple (d, x, y);
 }
 
 unsigned
diff --git a/algebra/algebra.h b/algebra/algebra.h
index 42319f5..4acd0ca 100644
--- a/algebra/algebra.h
+++ b/algebra/algebra.h
@@ -1,6 +1,4 @@
 
-#include <gmp.h>
-
 #include <lib/lib.h>
 
 inline int recip (int x)
@@ -42,9 +40,9 @@ unsigned int_lcm (int a, int b);
 uint64 int64_lcm (int64 a, int64 b);
 
 // (d, x, y) = gcd (a, b) where x*a + y*b = d
-triple<unsigned, int, int> extended_gcd (int a, int b);
+tuple<unsigned, int, int> extended_gcd (int a, int b);
 
-triple<unsigned, int, int> unsigned_extended_gcd (unsigned a, unsigned b);
+tuple<unsigned, int, int> unsigned_extended_gcd (unsigned a, unsigned b);
 
 template<class R> class linear_combination;
 template<class R> class linear_combination_const_iter;
diff --git a/algebra/fraction_field.h b/algebra/fraction_field.h
index d2ea6a4..b96ac84 100644
--- a/algebra/fraction_field.h
+++ b/algebra/fraction_field.h
@@ -36,6 +36,7 @@ template<class T> class fraction_field
   fraction_field &operator = (int x) { num = x; denom = 1; return *this; }
   
   bool operator == (const fraction_field &q) const { return num * q.denom == q.num * denom; }
+  bool operator != (const fraction_field &q) const { return !operator == (q); }
   
   bool operator == (int x) const { return num == denom * T (x); }
   bool operator != (int x) const { return !operator == (x); }
@@ -187,12 +188,12 @@ template<class T> class fraction_field
     return fraction_field (1);
   }
   
-  triple<fraction_field, fraction_field, fraction_field> extended_gcd (const fraction_field &x) const
+  tuple<fraction_field, fraction_field, fraction_field> extended_gcd (const fraction_field &x) const
   {
     if (*this != 0)
-      return triple<fraction_field, fraction_field, fraction_field> (*this, 1, 0);
+      return make_tuple (*this, fraction_field (1), fraction_field (0));
     else
-      return triple<fraction_field, fraction_field, fraction_field> (x, 0, 1);
+      return make_tuple (x, fraction_field (0), fraction_field (1));
   }
   
   static void show_ring () { printf ("fraction_field("); T::show_ring (); printf (")"); }
diff --git a/algebra/grading.h b/algebra/grading.h
index 936aeb3..4fec9b0 100644
--- a/algebra/grading.h
+++ b/algebra/grading.h
@@ -36,6 +36,8 @@ class grading
   }
   
   bool operator == (const grading &gr) const { return h == gr.h && q == gr.q; }
+  bool operator != (const grading &gr) const { return !operator == (gr); }
+  
   bool operator < (const grading &gr) const
   {
     return h < gr.h
diff --git a/algebra/linear_combination.h b/algebra/linear_combination.h
index dda27ce..c3d6bba 100644
--- a/algebra/linear_combination.h
+++ b/algebra/linear_combination.h
@@ -55,6 +55,7 @@ class linear_combination
     lc3 -= lc2;
     return lc3 == 0;
   }
+  bool operator != (const linear_combination &lc) const { return !operator == (lc); }
   
   bool operator == (int x) const
   {
diff --git a/algebra/module.h b/algebra/module.h
index 610802b..71f3eff 100644
--- a/algebra/module.h
+++ b/algebra/module.h
@@ -31,10 +31,10 @@ class module : public refcounted
     id = id_counter;
     id_counter ++;
   }
-  module (const module &); // doesn't exist
+  module (const module &) = delete;
   virtual ~module () { }
   
-  module &operator = (const module &); // doesn't exist
+  module &operator = (const module &) = delete;
   
  public:
   // the number of generators; n
@@ -466,11 +466,11 @@ class base_module : public module<R>
   G g;
   
  public:
-  base_module (); // doesn't exist
+  base_module () = delete;
   base_module (const G &g_) : g(g_) { }
   ~base_module () { }
   
-  base_module &operator = (const base_module &); // doesn't exist
+  base_module &operator = (const base_module &) = delete;
   
   unsigned dim () const { return g.dim (); }
   unsigned free_rank () const { return g.free_rank (); }
@@ -487,7 +487,7 @@ class explicit_module : public module<R>
   basedvector<grading, 1> hq;
   
  public:
-  explicit_module (); // doesn't exist
+  explicit_module () = delete;
   explicit_module (unsigned r_,
 		   basedvector<R, 1> ann_,
 		   basedvector<grading, 1> hq_)
@@ -497,7 +497,7 @@ class explicit_module : public module<R>
   explicit explicit_module (unsigned r_, basedvector<grading, 1> hq_) : r(r_), hq(hq_) { }
   ~explicit_module () { }
   
-  explicit_module &operator = (const explicit_module &); // doesn't exist
+  explicit_module &operator = (const explicit_module &) = delete;
   
   unsigned dim () const { return r + ann.size (); }
   unsigned free_rank () const { return r; }
@@ -525,7 +525,7 @@ class free_submodule : public module<R>
   { }
   ~free_submodule () { }
   
-  free_submodule &operator = (const free_submodule &); // doesn't exist
+  free_submodule &operator = (const free_submodule &) = delete;
   
   ptr<const module<R> > parent_module () const { return parent; }
   
@@ -562,13 +562,13 @@ class quotient_module : public module<R>
   basedvector<map<unsigned, R>, 1> pi;
   
  public:
-  quotient_module (const quotient_module &); // doesn't exist
+  quotient_module (const quotient_module &) = delete;
   quotient_module (ptr<const module<R> > parent_)
     : parent(parent_)
   { }
   ~quotient_module () { }
   
-  quotient_module &operator = (const quotient_module &); // doesn't exsit
+  quotient_module &operator = (const quotient_module &) = delete;
   
   ptr<const module<R> > parent_module () const { return parent; }
   
@@ -617,7 +617,7 @@ class map_impl : public refcounted
   ptr<const module<R> > to;
   
  public:
-  map_impl (const map_impl &); // doesn't exist
+  map_impl (const map_impl &) = delete;
   map_impl (ptr<const module<R> > fromto)
     : from(fromto), to(fromto)
   { }
@@ -626,7 +626,7 @@ class map_impl : public refcounted
   { }
   virtual ~map_impl () { }
   
-  map_impl &operator = (const map_impl &); // doesn't exist
+  map_impl &operator = (const map_impl &) = delete;
   
   virtual linear_combination<R> column (unsigned i) const = 0;
   
@@ -867,6 +867,7 @@ class mod_map
       }
     return 1;
   }
+  bool operator != (const mod_map &m) const { return !operator == (m); }
   
   bool operator == (int x) const
   {
@@ -1095,16 +1096,16 @@ quotient_helper<R>::improve_pivot_row (unsigned i, unsigned j, unsigned i2)
     }
 #endif
   
-  triple<R, R, R> t = rc.extended_gcd (r2c);
-  assert (t.first == rc*t.second + t.third*r2c);
+  tuple<R, R, R> t = rc.extended_gcd (r2c);
+  assert (get<0> (t) == rc*get<1> (t) + get<2> (t)*r2c);
   
-  rows[i] = r*t.second + r2*t.third;
-  rows[i2] = (rc.div (t.first))*r2 - (r2c.div (t.first))*r;
+  rows[i] = r*get<1> (t) + r2*get<2> (t);
+  rows[i2] = (rc.div (get<0> (t)))*r2 - (r2c.div (get<0> (t)))*r;
   
-  assert ((rc | r2c) == rc.divides (t.first));
-  assert (!rc.divides (t.first) || rows[i2](j) == 0);
+  assert ((rc | r2c) == rc.divides (get<0> (t)));
+  assert (!rc.divides (get<0> (t)) || rows[i2](j) == 0);
   
-  return !rc.divides (t.first);
+  return !rc.divides (get<0> (t));
 }
 
 template<class R> bool
@@ -1149,8 +1150,8 @@ quotient_helper<R>::improve_pivot_column (unsigned i, unsigned j, unsigned j2)
     }
 #endif
   
-  triple<R, R, R> t = rc.extended_gcd (rc2);
-  assert (t.first == rc*t.second + t.third*rc2);
+  tuple<R, R, R> t = rc.extended_gcd (rc2);
+  assert (get<0> (t) == rc*get<1> (t) + get<2> (t)*rc2);
   
   for (unsigned k = 1; k <= rows.size (); k ++)
     {
@@ -1158,9 +1159,9 @@ quotient_helper<R>::improve_pivot_column (unsigned i, unsigned j, unsigned j2)
       R rkc = rk(j),
 	rkc2 = rk(j2);
       
-      rk.set_coeff (rkc*t.second + rkc2*t.third,
+      rk.set_coeff (rkc*get<1> (t) + rkc2*get<2> (t),
 		    j);
-      rk.set_coeff (rkc2*(rc.div (t.first)) - rkc*(rc2.div (t.first)),
+      rk.set_coeff (rkc2*(rc.div (get<0> (t))) - rkc*(rc2.div (get<0> (t))),
 		    j2);
     }
   
@@ -1169,8 +1170,8 @@ quotient_helper<R>::improve_pivot_column (unsigned i, unsigned j, unsigned j2)
   
   assert (g.hq () == g2.hq ());
   
-  generators[j] = (rc.div (t.first)) * g + (rc2.div (t.first)) * g2;
-  generators[j2] = t.second * g2 - t.third * g;
+  generators[j] = (rc.div (get<0> (t))) * g + (rc2.div (get<0> (t))) * g2;
+  generators[j2] = get<1> (t) * g2 - get<2> (t) * g;
   
  #if 0
   for (unsigned k = 1; k <= rows.size (); k ++)
@@ -1189,8 +1190,8 @@ quotient_helper<R>::improve_pivot_column (unsigned i, unsigned j, unsigned j2)
       R d = ginv(j),
 	d2 = ginv(j2);
       
-      ginv.set_coeff (t.second*d + t.third*d2, j);
-      ginv.set_coeff (rc.div (t.first) * d2 - rc2.div (t.first) * d, j2);
+      ginv.set_coeff (get<1> (t)*d + get<2> (t)*d2, j);
+      ginv.set_coeff (rc.div (get<0> (t)) * d2 - rc2.div (get<0> (t)) * d, j2);
     }
   
 #if 0
@@ -1207,10 +1208,10 @@ quotient_helper<R>::improve_pivot_column (unsigned i, unsigned j, unsigned j2)
     }
 #endif
   
-  assert ((rc | rc2) == rc.divides (t.first));
-  assert (!rc.divides (t.first) || r(j2) == 0);
+  assert ((rc | rc2) == rc.divides (get<0> (t)));
+  assert (!rc.divides (get<0> (t)) || r(j2) == 0);
   
-  return !rc.divides (t.first);
+  return !rc.divides (get<0> (t));
 }
 
 template<class R> void
@@ -1595,7 +1596,7 @@ free_submodule<R>::restrict_submodule (ptr<const free_submodule<R> > m) const
     span[i] = restrict (m->inject_generator (i));
   
   mod_span<R> span2 (this, span);
-  return submodule (span2);
+  return this->submodule (span2);
 }
 
 template<class R> bool
@@ -1672,11 +1673,11 @@ mod_map<R>::kernel () const
 	  R to_xc = to_x(i);
 	  if (! (to_vc | to_xc))
 	    {
-	      triple<R, R, R> t = to_vc.extended_gcd (to_xc);
-	      assert (t.first == to_vc*t.second + t.third*to_xc);
+	      tuple<R, R, R> t = to_vc.extended_gcd (to_xc);
+	      assert (get<0> (t) == to_vc*get<1> (t) + get<2> (t)*to_xc);
 	      
-	      to_v = t.second*to_v + t.third*to_x;
-	      from_v = t.second*from_v + t.third*from_x;
+	      to_v = get<1> (t)*to_v + get<2> (t)*to_x;
+	      from_v = get<1> (t)*from_v + get<2> (t)*from_x;
 	      
 	      assert (to_v(i) != 0);
 	    }
@@ -1770,10 +1771,10 @@ mod_span<R>::mod_span (ptr<const module<R> > mod,
 	  
 	  if (! (vc | xc))
 	    {
-	      triple<R, R, R> t = vc.extended_gcd (xc);
-	      assert (t.first == vc*t.second + t.third*xc);
+	      tuple<R, R, R> t = vc.extended_gcd (xc);
+	      assert (get<0> (t) == vc*get<1> (t) + get<2> (t)*xc);
 	      
-	      v = t.second*v + t.third*x;
+	      v = get<1> (t)*v + get<2> (t)*x;
 	      
 	      assert (v(i) != 0);
 	    }
diff --git a/algebra/multivariate_laurentpoly.h b/algebra/multivariate_laurentpoly.h
index 49c2137..42747c3 100644
--- a/algebra/multivariate_laurentpoly.h
+++ b/algebra/multivariate_laurentpoly.h
@@ -55,6 +55,7 @@ class multivariate_laurent_monomial
 #endif
     return m == e.m;
   }
+  bool operator != (const multivariate_laurent_monomial &e) const { return !operator == (e); }
   
   bool operator == (int x) const
   {
@@ -243,6 +244,7 @@ class multivariate_laurentpoly
 #endif
     return coeffs == p.coeffs;
   }
+  bool operator != (const multivariate_laurentpoly &p) const { return !operator == (p); }
   
   bool operator == (int x) const
   {
diff --git a/algebra/multivariate_polynomial.h b/algebra/multivariate_polynomial.h
index 75aac48..309a03e 100644
--- a/algebra/multivariate_polynomial.h
+++ b/algebra/multivariate_polynomial.h
@@ -56,6 +56,7 @@ class multivariate_monomial
       }
     return 1;
   }
+  bool operator != (const multivariate_monomial &e) const { return !operator == (e); }
   
   bool operator < (const multivariate_monomial &e) const
   {
@@ -199,7 +200,7 @@ public:
   
   bool is_unit () const;
   
-  bool operator == (multivariate_polynomial p) const
+  bool operator == (const multivariate_polynomial &p) const
   {
 #ifndef NDEBUG
     check ();
@@ -207,6 +208,7 @@ public:
 #endif
     return coeffs == p.coeffs;
   }
+  bool operator != (const multivariate_polynomial &p) const { return !operator == (p); }
   
   bool operator == (int x) const
   {
diff --git a/algebra/polynomial.h b/algebra/polynomial.h
index 4aa83c5..3952544 100644
--- a/algebra/polynomial.h
+++ b/algebra/polynomial.h
@@ -41,6 +41,7 @@ public:
   
   polynomial (const polynomial &p) : coeffs(p.coeffs) { }
   polynomial (copy, const polynomial &p) : coeffs(COPY, p.coeffs) { }
+  polynomial (reader &r) : coeffs(r) { }
   ~polynomial () { }
   
   polynomial &operator = (const polynomial &p) { coeffs = p.coeffs; return *this; }
@@ -64,7 +65,16 @@ public:
   unsigned degree () const;
   bool is_unit () const;
   
-  bool operator == (polynomial p) const
+  polynomial recip () const
+  {
+    assert (coeffs.card () != 1);
+    
+    pair<unsigned, T> p = coeffs.head ();
+    assert (p.first == 0);
+    return polynomial (p.second.recip ());
+  }
+  
+  bool operator == (const polynomial &p) const
   {
 #ifndef NDEBUG
     check ();
@@ -72,6 +82,7 @@ public:
 #endif
     return coeffs == p.coeffs;
   }
+  bool operator != (const polynomial &p) const { return !operator == (p); }
   
   bool operator == (int x) const
   {
@@ -112,27 +123,33 @@ public:
     return r;
   }
   
-  polynomial operator - (polynomial p) const
+  polynomial operator - (const polynomial &p) const
   {
     polynomial r (COPY, *this);
     r -= p;
     return r;
   }
   
-  polynomial operator * (polynomial p) const;
+  polynomial operator * (const polynomial &p) const;
   
-  pair<polynomial, polynomial> divide_with_remainder (polynomial d) const;
+  bool divides (const polynomial &num) const;
+  bool operator | (const polynomial &num) const { return divides (num); }
   
-  polynomial mod (polynomial d) const;
-  bool divides (polynomial d) const;
-  polynomial divide_exact (polynomial d) const;
+  tuple<polynomial, polynomial> divide_with_remainder (const polynomial &denom) const;
   
-  polynomial gcd (polynomial b) const;
+  polynomial mod (const polynomial &denom) const;
+  
+  polynomial divide_exact (const polynomial &denom) const;
+  
+  polynomial gcd (const polynomial &b) const;
+  polynomial lcm (const polynomial &b) const;
+  tuple<polynomial, polynomial, polynomial> extended_gcd (const polynomial &b) const;
   
 #ifndef NDEBUG
   void check () const;
 #endif
   
+  void write_self (writer &w) const { write (w, coeffs); }
   static void show_ring () { T::show_ring (); printf ("[x]"); }
   void display_self () const { show_self (); newline (); }
   void show_self () const;
@@ -197,7 +214,7 @@ polynomial<T>::operator *= (T s)
 }
 
 template<class T> polynomial<T>
-polynomial<T>::operator * (polynomial p) const
+polynomial<T>::operator * (const polynomial &p) const
 {
   polynomial r;
   
@@ -213,14 +230,15 @@ polynomial<T>::operator * (polynomial p) const
   return r;
 }
 
-template<class T> pair<polynomial<T>, polynomial<T> > 
-polynomial<T>::divide_with_remainder (polynomial d) const
+template<class T> tuple<polynomial<T>, polynomial<T> > 
+polynomial<T>::divide_with_remainder (const polynomial &d) const
 {
+  // num = *this
   assert (d != 0);
   
   polynomial r (COPY, *this);
   polynomial q;
-
+  
   pair<unsigned, T> d_leading_term = d.coeffs.tail ();
   for (;;)
     {
@@ -240,45 +258,80 @@ polynomial<T>::divide_with_remainder (polynomial d) const
   assert (r == 0 || r.degree () < d.degree ());
   // assert (*this == q*d + r);
   
-  return pair<polynomial, polynomial> (q, r);
+  return make_tuple (q, r);
 }
 
 template<class T> polynomial<T>
-polynomial<T>::mod (polynomial d) const
+polynomial<T>::mod (const polynomial &denom) const
 {
-  pair<polynomial<T>, polynomial<T> > qr = divide_with_remainder (d);
-  return qr.second;
+  polynomial q, r;
+  tie (q, r) = divide_with_remainder (denom);
+  return r;
 }
 
 template<class T> bool
-polynomial<T>::divides (polynomial d) const
+polynomial<T>::divides (const polynomial &num) const
 {
-  pair<polynomial<T>, polynomial<T> > qr = divide_with_remainder (d);
-  return qr.second == 0;
+  // denom = *this
+  polynomial q, r;
+  tie (q, r) = num.divide_with_remainder (*this);
+  return r == 0;
 }
 
 template<class T> polynomial<T>
-polynomial<T>::divide_exact (polynomial d) const
+polynomial<T>::divide_exact (const polynomial &denom) const
 {
-  pair<polynomial<T>, polynomial<T> > qr = divide_with_remainder (d);
-  assert (qr.second == 0);
-  return qr.first;
+  polynomial q, r;
+  tie (q, r) = divide_with_remainder (denom);
+  assert (r == 0);
+  return q;
 }
 
 template<class T> polynomial<T>
-polynomial<T>::gcd (polynomial b) const
+polynomial<T>::gcd (const polynomial &b) const
 {
   polynomial a = *this;
-  
   while (b != 0)
     {
-      pair<polynomial<T>, polynomial<T> > a_qr = a.divide_with_remainder (b);
+      polynomial r;
+      tie (ignore, r) = a.divide_with_remainder (b);
       a = b;
-      b = a_qr.second;
+      b = r;
     }
   return a;
 }
 
+template<class T> polynomial<T>
+polynomial<T>::lcm (const polynomial &b) const
+{
+  // a = *this
+  return divide_exact (gcd (b)) * b;
+}
+
+template<class T> tuple<polynomial<T>, polynomial<T>, polynomial<T> >
+polynomial<T>::extended_gcd (const polynomial &b) const
+{
+  // a = *this
+  if (b == 0)
+    return make_tuple (*this, polynomial (1), polynomial (0));
+  
+  polynomial q, r;
+  tie (q, r) = divide_with_remainder (b);
+  if (r == 0)
+    return make_tuple (b, polynomial (0), polynomial (1));
+  else
+    {
+      polynomial d, s, t;
+      
+      tie (d, s, t) = b.extended_gcd (r);
+      
+      polynomial s2 = t,
+	t2 = s - t*q;
+      
+      return make_tuple (d, s2, t2);
+    }
+}
+
 #ifndef NDEBUG
 template<class T> void
 polynomial<T>::check () const
@@ -380,7 +433,8 @@ class polynomial<Z2>
     return *this;
   }
   
-  bool operator == (polynomial p) const { return coeffs == p.coeffs; }
+  bool operator == (const polynomial &p) const { return coeffs == p.coeffs; }
+  bool operator != (const polynomial &p) const { return !operator == (p); }
   
   bool operator == (int x) const
   {
@@ -443,4 +497,3 @@ class polynomial<Z2>
   void display_self () const { show_self (); newline (); }
   void show_self () const;
 };
-
diff --git a/cobordism.cpp b/cobordism.cpp
index 55acf28..72e8200 100644
--- a/cobordism.cpp
+++ b/cobordism.cpp
@@ -19,10 +19,10 @@ public:
   
 public:
   rd_unify ();
-  rd_unify (const rd_unify &); // doesn't exist
+  rd_unify (const rd_unify &) = delete;
   ~rd_unify () { }
   
-  rd_unify &operator = (const rd_unify &); // doesn't exist
+  rd_unify &operator = (const rd_unify &) = delete;
   
   bool unify (const resolution_diagram &rd1, const resolution_diagram &rd2);
   unsigned map_starting_monomial (unsigned m);
diff --git a/cobordism.h b/cobordism.h
index e8dfc47..0365343 100644
--- a/cobordism.h
+++ b/cobordism.h
@@ -160,6 +160,7 @@ class resolution_diagram
   resolution_diagram &operator = (const resolution_diagram &rd);
   
   bool operator == (const resolution_diagram &rd) const;
+  bool operator != (const resolution_diagram &rd) const { return !operator == (rd); }
   
   /* reverse, reverse_orientation leave starting/ending circle numbers
      unchanged. */
@@ -190,10 +191,10 @@ class resolution_diagram_builder
   
  public:
   resolution_diagram_builder ();
-  resolution_diagram_builder (const resolution_diagram_builder &); // doesn't exist
+  resolution_diagram_builder (const resolution_diagram_builder &) = delete;
   ~resolution_diagram_builder () { }
   
-  resolution_diagram_builder &operator = (const resolution_diagram_builder &); // doesn't exist
+  resolution_diagram_builder &operator = (const resolution_diagram_builder &) = delete;
   
   void init (const knot_diagram &d,
 	     smallbitset from_state, const smoothing &from_s,
diff --git a/cube.h b/cube.h
index b3883ac..49bd26e 100644
--- a/cube.h
+++ b/cube.h
@@ -3,10 +3,10 @@ class map_rules
 {
  public:
   map_rules () { }
-  map_rules (const map_rules &); // doesn't exist
+  map_rules (const map_rules &) = delete;
   virtual ~map_rules () { }
   
-  map_rules &operator = (const map_rules &);  // doesn't exist
+  map_rules &operator = (const map_rules &) = delete;
   
   virtual void map (basedvector<pair<unsigned, unsigned>, 1> &out,
 		    resolution_diagram_builder &rdb) const = 0;
@@ -77,10 +77,10 @@ class twisted_map_rules
 {
  public:
   twisted_map_rules () { }
-  twisted_map_rules (const twisted_map_rules &); // doesn't exist
+  twisted_map_rules (const twisted_map_rules &) = delete;
   virtual ~twisted_map_rules () { }
   
-  map_rules &operator = (const twisted_map_rules &);  // doesn't exist
+  map_rules &operator = (const twisted_map_rules &) = delete;
   
   virtual void map (basedvector<triple<unsigned, unsigned, set<unsigned> >, 1> &out,
 		    resolution_diagram_builder &rdb) const = 0;
@@ -100,10 +100,10 @@ class twisted_cube
     : c(c_)
   { }
   
-  twisted_cube (const twisted_cube &); // doesn't exist
+  twisted_cube (const twisted_cube &) = delete;
   ~twisted_cube () { }
   
-  twisted_cube &operator = (const twisted_cube &); // doesn't exist
+  twisted_cube &operator = (const twisted_cube &) = delete;
   
   mod_map<R> compute_twisted_map (basedvector<int, 1> edge_weight,
 				  unsigned dh,
diff --git a/cube_impl.h b/cube_impl.h
index 87ae5f8..b930d68 100644
--- a/cube_impl.h
+++ b/cube_impl.h
@@ -9,7 +9,7 @@ public:
   khC_generators (const cube<R> &c_) : c(c_) { }
   ~khC_generators () { }
   
-  khC_generators &operator = (const khC_generators &); // doesn't exist
+  khC_generators &operator = (const khC_generators &) = delete;
   
   unsigned dim () const { return c.n_generators; }
   unsigned free_rank () const { return c.n_generators; }
@@ -187,10 +187,10 @@ class d_rules : public map_rules
 {
 public:
   d_rules () { }
-  d_rules (const d_rules &); // doesn't exist
+  d_rules (const d_rules &) = delete;
   ~d_rules () { }
   
-  d_rules &operator = (const d_rules &); // doesn't exist
+  d_rules &operator = (const d_rules &) = delete;
   
   void map (basedvector<pair<unsigned, unsigned>, 1> &out,
 	    resolution_diagram_builder &rdb) const
@@ -214,10 +214,10 @@ class twin_arrows_P_rules : public map_rules
 {
 public:
   twin_arrows_P_rules () { }
-  twin_arrows_P_rules (const twin_arrows_P_rules &); // doesn't exist
+  twin_arrows_P_rules (const twin_arrows_P_rules &) = delete;
   ~twin_arrows_P_rules () { }
   
-  twin_arrows_P_rules &operator = (const twin_arrows_P_rules &); // doesn't exist
+  twin_arrows_P_rules &operator = (const twin_arrows_P_rules &) = delete;
   
   void map (basedvector<pair<unsigned, unsigned>, 1> &out,
 	    resolution_diagram_builder &rdb) const
@@ -676,10 +676,10 @@ class twisted_barE_rules : public twisted_map_rules
 {
 public:
   twisted_barE_rules () { }
-  twisted_barE_rules (const twisted_barE_rules &); // doesn't exist
+  twisted_barE_rules (const twisted_barE_rules &) = delete;
   ~twisted_barE_rules () { }
   
-  twisted_barE_rules &operator = (const twisted_barE_rules &); // doesn't exist
+  twisted_barE_rules &operator = (const twisted_barE_rules &) = delete;
   
   void map (basedvector<triple<unsigned, unsigned, set<unsigned> >, 1> &out,
 	    resolution_diagram_builder &rdb) const
diff --git a/knot_diagram.cpp b/knot_diagram.cpp
index 5b5e803..97f4ba8 100644
--- a/knot_diagram.cpp
+++ b/knot_diagram.cpp
@@ -82,10 +82,10 @@ public:
   
 public:
   dt_layout (const dt_code &dt, knot_diagram &kd_);
-  dt_layout (const dt_layout &); //  doesn't exist
+  dt_layout (const dt_layout &) = delete;
   ~dt_layout () { }
 
-  dt_layout &operator = (const dt_layout &); // doesn't exist
+  dt_layout &operator = (const dt_layout &) = delete;
   
   unsigned find_crossing (unsigned prevc, unsigned previ, unsigned target,
 			  bool under, unsigned dir);
diff --git a/lib/hashmap.h b/lib/hashmap.h
index 1616840..845da86 100644
--- a/lib/hashmap.h
+++ b/lib/hashmap.h
@@ -1,16 +1,10 @@
 
 /* wrapper for std::unordered_map */
 
-template<class K>
-struct hasher : public std::unary_function<K, hash_t>
-{
-  hash_t operator () (const K &k) const { return hash (k); }
-};
-
 template<class K, class V>
-class hashmap : public map_wrapper<std::tr1::unordered_map<K, V, hasher<K> >, K, V>
+class hashmap : public map_wrapper<std::unordered_map<K, V, hasher<K> >, K, V>
 {
-  typedef map_wrapper<std::tr1::unordered_map<K, V, hasher<K> >, K, V> base;
+  typedef map_wrapper<std::unordered_map<K, V, hasher<K> >, K, V> base;
   
  public:
   hashmap () { }
@@ -21,3 +15,9 @@ class hashmap : public map_wrapper<std::tr1::unordered_map<K, V, hasher<K> >, K,
   
   hashmap &operator = (const hashmap &m) { base::operator = (m); return *this; }
 };
+
+template<class K, class V>
+using hashmap_iter = map_wrapper_iter<std::unordered_map<K, V>, K, V>;
+
+template<class K, class V>
+using hashmap_const_iter = map_wrapper_const_iter<std::unordered_map<K, V>, K, V>;
diff --git a/lib/hashset.h b/lib/hashset.h
new file mode 100644
index 0000000..2c0bc97
--- /dev/null
+++ b/lib/hashset.h
@@ -0,0 +1,23 @@
+
+/* wrapper for std::unordered_set */
+
+template<class T>
+class hashset : public set_wrapper<std::unordered_set<T, hasher<T> >, T>
+{
+  typedef set_wrapper<std::unordered_set<T, hasher<T> >, T> base;
+  
+ public:
+  hashset () { }
+  hashset (const hashset &m) : base(m) { }
+  hashset (copy, const hashset &m) : base(COPY, m) { }
+  hashset (reader &r) : base(r) { }
+  ~hashset () { }
+  
+  hashset &operator = (const hashset &m) { base::operator = (m); return *this; }
+};
+
+template<class T>
+using hashset_iter = set_wrapper_iter<std::unordered_set<T, hasher<T> >, T>;
+
+template<class T>
+using hashset_const_iter = set_wrapper_const_iter<std::unordered_set<T, hasher<T> >, T>;
diff --git a/lib/io.cpp b/lib/io.cpp
index a02ccd0..3ea5068 100644
--- a/lib/io.cpp
+++ b/lib/io.cpp
@@ -1,10 +1,173 @@
 
 #include <lib/lib.h>
 
+void
+writer::write_int (int x)
+{
+  uint8 buf[5];
+  unsigned n = 0;
+  
+  bool more = 1;
+  while (more)
+    {
+      uint8 b = (uint8)(x & 0x7f);
+      x >>= 7;
+      
+      if ((x == 0
+	   && ! (b & 0x40))
+	  || (x == -1
+	      && (b & 0x40) == 0x40))
+	more = 0;
+      else
+	b |= 0x80;
+      
+      assert (n < 5);
+      buf[n] = b;
+      n ++;
+    }
+  
+  write_raw (buf, sizeof buf[0], n);
+}
+
+void
+writer::write_unsigned (unsigned x)
+{
+  uint8 buf[5];
+  unsigned n = 0;
+  
+  bool more = 1;
+  while (more)
+    {
+      uint8 b = (uint8)(x & 0x7f);
+      x >>= 7;
+      
+      if ((x == 0
+	   && ! (b & 0x40)))
+	more = 0;
+      else
+	b |= 0x80;
+      
+      assert (n < 5);
+      buf[n] = b;
+      n ++;
+    }
+  
+  write_raw (buf, sizeof buf[0], n);
+}
+
+void
+writer::write_uint64 (uint64 x)
+{
+  uint8 buf[10];
+  unsigned n = 0;
+  
+  bool more = 1;
+  while (more)
+    {
+      uint8 b = (uint8)(x & 0x7f);
+      x >>= 7;
+      
+      if ((x == 0
+	   && ! (b & 0x40)))
+	more = 0;
+      else
+	b |= 0x80;
+
+      assert (n < 10);
+      buf[n] = b;
+      n ++;
+    }
+  
+  write_raw (buf, sizeof buf[0], n);
+}
+
+void
+writer::write_mpz (const mpz_t x)
+{
+  size_t count;
+  void *buf = mpz_export (nullptr, &count, -1, 1, -1, 0, x);
+  
+  write_unsigned ((unsigned)count);
+  write_raw (buf, 1, count);
+  
+  free (buf);
+}
+
+int
+reader::read_int ()
+{
+  int x = 0;
+  int shift = 0;
+  for (;;)
+    {
+      uint8 b = read_uint8 ();
+      x |= ((int)(b & 0x7f) << shift);
+      shift += 7;
+      if (! (b & 0x80))
+	{
+	  if (shift < int_bits
+	      && (b & 0x40))
+	    x = (x << (int_bits - shift)) >> (int_bits - shift);
+	  break;
+	}
+    }
+  return x;
+}
+
+unsigned
+reader::read_unsigned ()
+{
+  unsigned x = 0;
+  unsigned shift = 0;
+  for (;;)
+    {
+      uint8 b = read_uint8 ();
+      x |= ((unsigned)(b & 0x7f) << shift);
+      shift += 7;
+      if (! (b & 0x80))
+	break;
+    }
+  return x;
+}
+
+uint64
+reader::read_uint64 ()
+{
+  uint64 x = 0;
+  uint64 shift = 0;
+  for (;;)
+    {
+      uint8 b = read_uint8 ();
+      x |= ((uint64)(b & 0x7f) << shift);
+      shift += 7;
+      if (! (b & 0x80))
+	break;
+    }
+  return x;
+}
+
+void
+reader::read_mpz (mpz_t x)
+{
+  unsigned count = read_unsigned ();
+  void *p = malloc (count);
+  if (!p)
+    {
+      stderror ("malloc");
+      exit (EXIT_FAILURE);
+    }
+  
+  read_raw (p, 1, count);
+  
+  mpz_import (x, count, -1, 1, -1, 0, p);
+  
+  free (p);
+}
+
 FILE *open_file (const std::string &file, const char *mode)
 {
   FILE *fp = fopen (file.c_str (), mode);
-  if (fp == 0)
+  if (!fp)
     {
       stderror ("fopen: %s", file.c_str ());
       exit (EXIT_FAILURE);
@@ -17,41 +180,32 @@ void close_file (FILE *fp)
   fclose (fp);
 }
 
-
-writer::writer (const std::string &file)
-  : fp(0)
+gzFile
+open_gzfile (const std::string &file, const char *mode)
 {
-  fp = open_file (file, "w");
-}
-
-writer::~writer ()
-{
-  if (fp)
+  gzFile gzfp = gzopen (file.c_str (), mode);
+  if (!gzfp)
     {
-      fclose (fp);
-      fp = 0;
+      stderror ("gzopen: %s", file.c_str ());
+      exit (EXIT_FAILURE);
     }
+  return gzfp;
 }
 
-reader::reader (const std::string &file)
-  : fp(0)
+extern void close_gzfile (gzFile gzfp)
 {
-  fp = open_file (file, "r");
-}
-
-reader::~reader ()
-{
-  if (fp)
+  int r = gzclose (gzfp);
+  if (r != Z_OK)
     {
-      fclose (fp);
-      fp = 0;
+      stderror ("gzclose");
+      exit (EXIT_FAILURE);
     }
 }
 
 void
-writer::write_char (char x)
+file_writer::write_raw (const void *p, size_t itemsize, size_t nitems)
 {
-  if (fwrite (&x, sizeof x, 1, fp) != 1)
+  if (fwrite (p, itemsize, nitems, fp) != nitems)
     {
       stderror ("fwrite");
       exit (EXIT_FAILURE);
@@ -59,103 +213,35 @@ writer::write_char (char x)
 }
 
 void
-writer::write_bool (bool x)
+file_reader::read_raw (void *p, size_t itemsize, size_t nitems)
 {
-  if (fwrite (&x, sizeof x, 1, fp) != 1)
+  if (fread (p, itemsize, nitems, fp) != nitems)
     {
-      stderror ("fwrite");
+      stderror ("fread");
       exit (EXIT_FAILURE);
     }
 }
 
 void
-writer::write_int (int x)
+gzfile_writer::write_raw (const void *p, size_t itemsize, size_t nitems)
 {
-  if (fwrite (&x, sizeof x, 1, fp) != 1)
+  unsigned nbytes = itemsize*nitems;
+  if (gzwrite (gzfp, p, nbytes) != nbytes)
     {
-      stderror ("fwrite");
+      stderror ("gzwrite");
       exit (EXIT_FAILURE);
     }
 }
 
 void
-writer::write_unsigned (unsigned x)
+gzfile_reader::read_raw (void *p, size_t itemsize, size_t nitems)
 {
-  if (fwrite (&x, sizeof x, 1, fp) != 1)
-    {
-      stderror ("fwrite");
-      exit (EXIT_FAILURE);
-    }
-}
-
-void
-writer::write_uint64 (uint64 x)
-{
-  if (fwrite (&x, sizeof x, 1, fp) != 1)
-    {
-      stderror ("fwrite");
-      exit (EXIT_FAILURE);
-    }
-}
-
-bool
-reader::read_bool ()
-{
-  bool x;
-  if (fread (&x, sizeof x, 1, fp) != 1)
+  unsigned nbytes = itemsize*nitems;
+  if (gzread (gzfp, p, nbytes) != nitems)
     {
       stderror ("fread");
       exit (EXIT_FAILURE);
     }
-  return x;
-}
-
-char
-reader::read_char ()
-{
-  char x;
-  if (fread (&x, sizeof x, 1, fp) != 1)
-    {
-      stderror ("fread");
-      exit (EXIT_FAILURE);
-    }
-  return x;
-}
-
-int
-reader::read_int ()
-{
-  int x;
-  if (fread (&x, sizeof x, 1, fp) != 1)
-    {
-      stderror ("fread");
-      exit (EXIT_FAILURE);
-    }
-  return x;
-}
-
-unsigned
-reader::read_unsigned ()
-{
-  unsigned x;
-  if (fread (&x, sizeof x, 1, fp) != 1)
-    {
-      stderror ("fread");
-      exit (EXIT_FAILURE);
-    }
-  return x;
-}
-
-uint64
-reader::read_uint64 ()
-{
-  uint64 x;
-  if (fread (&x, sizeof x, 1, fp) != 1)
-    {
-      stderror ("fread");
-      exit (EXIT_FAILURE);
-    }
-  return x;
 }
 
 void
@@ -164,12 +250,7 @@ read (reader &r, std::string &s)
   unsigned n = r.read_unsigned ();
   
   char buf[n + 1];
-  unsigned k = fread (buf, sizeof (char), n + 1, r.fp);
-  if (k != n + 1)
-    {
-      stderror ("fread");
-      exit (EXIT_FAILURE);
-    }
+  r.read_raw (buf, sizeof (char), n + 1);
   assert (buf[n] == 0);
   
   s = std::string (buf);
@@ -181,9 +262,5 @@ write (writer &w, const std::string &s)
   unsigned n = s.length ();
   w.write_unsigned (n);
   
-  if (fwrite (s.c_str (), sizeof (char), n + 1, w.fp) != n + 1)
-    {
-      stderror ("fwrite");
-      exit (EXIT_FAILURE);
-    }
+  w.write_raw (s.c_str (), sizeof (char), n + 1);
 }
diff --git a/lib/io.h b/lib/io.h
index 9f499e9..e08a243 100644
--- a/lib/io.h
+++ b/lib/io.h
@@ -1,43 +1,138 @@
 
-extern FILE *open_file (const std::string &file, const char *mode);
-extern void close_file (FILE *fp);
-
 class writer
 {
  public:
-  FILE *fp;
+  writer (const writer &) = delete;
+  writer () = default;
+  virtual ~writer () = default;
   
- public:
-  writer (const std::string &file);
-  writer (const writer &); // doesn't exist
-  ~writer ();
+  writer &operator = (const writer &) = delete;
   
-  writer &operator = (const writer &); // doesn't exist
+  virtual void write_raw (const void *p, size_t itemsize, size_t nitems) = 0;
   
-  void write_bool (bool x);
-  void write_char (char x);
+  void write_bool (bool x) { write_raw (&x, sizeof x, 1); }
+  void write_char (char x) { write_raw (&x, sizeof x, 1); }
+  void write_uint8 (uint8 x) { write_raw (&x, sizeof x, 1); }
   void write_int (int x);
   void write_unsigned (unsigned x);
   void write_uint64 (uint64 x);
+  void write_mpz (const mpz_t x);
 };
 
 class reader
 {
  public:
-  FILE *fp;
-
- public:
-  reader (const std::string &file);
-  reader (const reader &); // doesn't exist
-  ~reader ();
+  reader (const reader &) = delete;
+  reader () = default;
+  virtual ~reader () = default;
   
-  reader &operator = (const reader &); // doesn't exist
+  reader &operator = (const reader &) = delete;
+  
+  virtual void read_raw (void *p, size_t itemsize, size_t nitems) = 0;
+  
+  bool read_bool ()
+  {
+    bool x;
+    read_raw (&x, sizeof x, 1);
+    return x;
+  }
+  
+  char read_char ()
+  {
+    char x;
+    read_raw (&x, sizeof x, 1);
+    return x;
+  }
+  
+  uint8 read_uint8 ()
+  {
+    uint8 x;
+    read_raw (&x, sizeof x, 1);
+    return x;
+  }
   
-  bool read_bool ();
-  char read_char ();
   int read_int ();
   unsigned read_unsigned ();
   uint64 read_uint64 ();
+  void read_mpz (mpz_t x);
+};
+
+extern FILE *open_file (const std::string &file, const char *mode);
+extern void close_file (FILE *fp);
+
+extern gzFile open_gzfile (const std::string &file, const char *mode);
+extern void close_gzfile (gzFile fp);
+
+class file_writer : public writer
+{
+ public:
+  FILE *fp;
+  
+ public:
+  file_writer (const std::string &file)
+    : fp(open_file (file, "w"))
+  {
+  }
+  file_writer (const file_writer &) = delete;
+  ~file_writer () { close_file (fp); }
+  
+  file_writer &operator = (const file_writer &) = delete;
+  
+  void write_raw (const void *p, size_t itemsize, size_t nitems);
+};
+
+class file_reader : public reader
+{
+ public:
+  FILE *fp;
+  
+ public:
+  file_reader (const std::string &file)
+    : fp(open_file (file, "r"))
+  {
+  }
+  file_reader (const file_reader &) = delete;
+  ~file_reader () { close_file (fp); }
+  
+  file_reader &operator = (const file_reader &) = delete;
+  
+  void read_raw (void *p, size_t itemsize, size_t nitems);
+};
+
+class gzfile_writer : public writer
+{
+ public:
+  gzFile gzfp;
+  
+ public:
+  gzfile_writer (const std::string &file)
+    : gzfp(open_gzfile (file, "w9"))
+  {
+  }
+  gzfile_writer (const gzfile_writer &) = delete;
+  ~gzfile_writer () { close_gzfile (gzfp); }
+  
+  gzfile_writer &operator = (const gzfile_writer &) = delete;
+  
+  void write_raw (const void *p, size_t itemsize, size_t nitems);
+};
+
+class gzfile_reader : public reader
+{
+ public:
+  gzFile gzfp;
+  
+ public:
+  gzfile_reader (const std::string &file)
+    : gzfp(open_gzfile (file, "r"))
+  {
+  }
+  gzfile_reader (const gzfile_reader &) = delete;
+  ~gzfile_reader () { close_gzfile (gzfp); }
+  
+  gzfile_reader &operator = (const gzfile_reader &) = delete;
+  
+  void read_raw (void *p, size_t itemsize, size_t nitems);
 };
 
 inline void read (reader &r, bool &x) { x = r.read_bool (); }
@@ -45,7 +140,7 @@ inline void read (reader &r, char &x) { x = r.read_char (); }
 inline void read (reader &r, int &x) { x = r.read_int (); }
 inline void read (reader &r, unsigned &x) { x = r.read_unsigned (); }
 inline void read (reader &r, uint64 &x) { x = r.read_uint64 (); }
-void read (reader &r, std::string &s);
+extern void read (reader &r, std::string &s);
 template<class T> inline void read (reader &r, T &x) { x = T(r); }
 
 inline void ctor_read (reader &r, bool *p) { *p = r.read_bool (); }
@@ -60,5 +155,5 @@ inline void write (writer &w, char x) { w.write_char (x); }
 inline void write (writer &w, int x) { w.write_int (x); }
 inline void write (writer &w, unsigned x) { w.write_unsigned (x); }
 inline void write (writer &w, uint64 x) { w.write_uint64 (x); }
-void write (writer &w, const std::string &s);
+extern void write (writer &w, const std::string &s);
 template<class T> inline void write (writer &w, const T &x) { x.write_self (w); }
diff --git a/lib/lib.h b/lib/lib.h
index e72b113..070b5e5 100644
--- a/lib/lib.h
+++ b/lib/lib.h
@@ -7,17 +7,20 @@
 #include <sys/stat.h>
 #include <sys/errno.h>
 
+#include <gmp.h>
+#include <zlib.h>
+
 #include <set>
 #include <map>
 #include <string>
 #include <queue>
 #include <algorithm>
 
-#include <tr1/functional>
-#include <tr1/unordered_map>
+#include <functional>
+#include <unordered_set>
+#include <unordered_map>
 
 /* just need to implement ==, < */
-template<class T> bool operator != (const T &a, const T &b) { return ! (a == b); }
 template<class T> bool operator <= (const T &a, const T &b) { return (a < b) || (a == b); }
 template<class T> bool operator > (const T &a, const T &b) { return ! (a <= b); }
 template<class T> bool operator >= (const T &a, const T &b) { return ! (a < b); }
@@ -25,6 +28,8 @@ template<class T> bool operator >= (const T &a, const T &b) { return ! (a < b);
 typedef enum copy { COPY, COPY1 } copy;
 typedef enum copy2 { COPY2 } copy2;
 
+typedef unsigned char uint8;
+
 typedef unsigned char uchar;
 typedef unsigned short ushort;
 typedef unsigned int uint;
@@ -35,6 +40,11 @@ typedef unsigned long long uint64;
 
 typedef size_t hash_t;
 
+static const unsigned uint8_bits = 8;
+static const unsigned char_bits = 8;
+static const unsigned uchar_bits = 8;
+static const unsigned int_bits = 32;
+
 inline hash_t hash (bool x) { return (hash_t)x; }
 inline hash_t hash (char x) { return (hash_t)x; }
 inline hash_t hash (unsigned char x) { return (hash_t)x; }
@@ -151,9 +161,19 @@ int alpha_to_int (char c);
 
 #include <lib/show.h>
 #include <lib/refcount.h>
+
+using std::tuple;
+using std::get;
+using std::make_tuple;
+using std::tie;
+using std::ignore;
+
 #include <lib/pair.h>
+
 #include <lib/maybe.h>
 #include <lib/vector.h>
+
+#include <lib/set_wrapper.h>
 #include <lib/set.h>
 
 class bitset;
@@ -175,7 +195,16 @@ inline int random_int (int from, int to)
 
 #include <lib/map_wrapper.h>
 #include <lib/map.h>
+
+template<class K>
+struct hasher : public std::unary_function<K, hash_t>
+{
+  hash_t operator () (const K &k) const { return hash (k); }
+};
+
 #include <lib/hashmap.h>
+#include <lib/hashset.h>
+
 #include <lib/ullmanmap.h>
 #include <lib/mapcommon.h>
 #include <lib/unionfind.h>
diff --git a/lib/map.h b/lib/map.h
index 3be8881..2521af3 100644
--- a/lib/map.h
+++ b/lib/map.h
@@ -1,5 +1,5 @@
 
-/* wrapper for stl map */
+/* wrapper for std::map */
 
 template<class K, class V>
 class map : public map_wrapper<std::map<K, V>, K, V>
@@ -17,3 +17,9 @@ class map : public map_wrapper<std::map<K, V>, K, V>
   
   map &operator = (const map &m) { base::operator = (m); return *this; }
 };
+
+template<class K, class V>
+using map_iter = map_wrapper_iter<std::map<K, V>, K, V>;
+
+template<class K, class V>
+using map_const_iter = map_wrapper_const_iter<std::map<K, V>, K, V>;
diff --git a/lib/map_wrapper.h b/lib/map_wrapper.h
index 762bfc7..48cfacf 100644
--- a/lib/map_wrapper.h
+++ b/lib/map_wrapper.h
@@ -1,5 +1,5 @@
 
-/* wrapper for stl maps */
+/* wrapper for STL maps */
 
 template<class M, class K, class V> class map_wrapper_iter;
 template<class M, class K, class V> class map_wrapper_const_iter;
@@ -7,10 +7,10 @@ template<class M, class K, class V> class map_wrapper_const_iter;
 template<class M, class K, class V>
 class map_wrapper
 {
- private:
   friend class map_wrapper_iter<M, K, V>;
   friend class map_wrapper_const_iter<M, K, V>;
-  
+
+ protected:
   class map_impl : public refcounted
   {
   public:
@@ -130,6 +130,7 @@ class map_wrapper
   }
   
   bool operator == (const map_wrapper &m) const;
+  bool operator != (const map_wrapper &m) const { return !operator == (m); }
   bool operator < (const map_wrapper &m) const;
   
   void write_self (writer &w) const;
@@ -222,19 +223,19 @@ template<class M, class K, class V>
 class map_wrapper_iter
 {
  private:
-  map_wrapper<M, K, V> &m;
+  ptr<typename map_wrapper<M, K, V>::map_impl> impl;
   typename M::iterator i, end;
   bool deleted;
   
  public:
-  map_wrapper_iter (map_wrapper<M, K, V> &m_) : m(m_), i(m_.impl->t.begin ()), end(m_.impl->t.end ()), deleted(0) { }
+  map_wrapper_iter (map_wrapper<M, K, V> &m) : impl(m.impl), i(m.impl->t.begin ()), end(m.impl->t.end ()), deleted(0) { }
   ~map_wrapper_iter () { }
   
   void del ()
   {
     assert (!deleted);
     typename M::iterator iprev = i ++;
-    m.impl->t.erase (iprev);
+    impl->t.erase (iprev);
     deleted = 1;
   }
   
@@ -249,10 +250,11 @@ template<class M, class K, class V>
 class map_wrapper_const_iter
 {
  private:
+  ptr<const typename map_wrapper<M, K, V>::map_impl> impl;
   typename M::const_iterator i, end;
   
  public:
-  map_wrapper_const_iter (const map_wrapper<M, K, V> &m) : i(m.impl->t.begin ()), end(m.impl->t.end ()) { }
+  map_wrapper_const_iter (const map_wrapper<M, K, V> &m) : impl(m.impl), i(m.impl->t.begin ()), end(m.impl->t.end ()) { }
   ~map_wrapper_const_iter () { }
   
   const K &key () const { return i->first; }
diff --git a/lib/pair.h b/lib/pair.h
index fb83a2f..ce37dbb 100644
--- a/lib/pair.h
+++ b/lib/pair.h
@@ -17,6 +17,7 @@ class pair
   pair &operator = (const pair &p) { first = p.first; second = p.second; return *this; }
   
   bool operator == (const pair &p) const { return first == p.first && second == p.second; }
+  bool operator != (const pair &p) const { return !operator == (p); }
   bool operator < (const pair &p) const { return (first < p.first
 						  || (first == p.first && second < p.second)); }
   
@@ -67,6 +68,8 @@ public:
   {
     return first == t.first && second == t.second && third == t.third;
   }
+  bool operator != (const triple &t) const { return !operator == (t); }
+  
   bool operator < (const triple &t) const
   {
     if (first < t.first)
diff --git a/lib/priority_queue.h b/lib/priority_queue.h
index 6b008f7..c332977 100644
--- a/lib/priority_queue.h
+++ b/lib/priority_queue.h
@@ -18,6 +18,8 @@ class priority_queue
     weighted &operator = (const weighted &w) { v = w.v; pri = w.pri; }
     
     bool operator == (const weighted &w) const { return v == w.v && pri == w.pri; }
+    bool operator != (const weighted &w) const { return !operator == (w); }
+    
     bool operator < (const weighted &w) const { return pri > w.pri || (pri == w.pri && v < w.v); }
   };
   
diff --git a/lib/refcount.h b/lib/refcount.h
index f638e63..b51da0d 100644
--- a/lib/refcount.h
+++ b/lib/refcount.h
@@ -42,10 +42,10 @@ class refcounted
   
  public:
   refcounted () : refcount(0) { }
-  refcounted (const refcounted &rc); // doesn't exist
+  refcounted (const refcounted &rc) = delete;
   ~refcounted () { }
   
-  refcounted &operator = (const refcounted &rc); // doesn't exist
+  refcounted &operator = (const refcounted &rc) = delete;
 };
 
 template<class T>
diff --git a/lib/set.h b/lib/set.h
index 2137655..a2aa03f 100644
--- a/lib/set.h
+++ b/lib/set.h
@@ -1,149 +1,45 @@
 
-/* wrap for std::set */
-
-template<class T> class set_iter;
-template<class T> class set_const_iter;
+/* wrapper for std::set */
 
 template<class T>
-class set
+class set : public set_wrapper<std::set<T>, T>
 {
- private:
-  friend class set_iter<T>;
-  friend class set_const_iter<T>;
-  
-  class set_impl : public refcounted
-  {
-  public:
-    std::set<T> t;
-  };
-  
-  ptr<set_impl> impl;
+  typedef set_wrapper<std::set<T>, T> base;
   
  public:
-  typedef set_iter<T> iter;
-  typedef set_const_iter<T> const_iter;
-  
- public:
-  set () : impl(new set_impl) { }
-  set (unsigned dummy_size) : impl(new set_impl) { }
-  set (const set &s) : impl(s.impl) { }
-  set (copy, const set &s) : impl(new set_impl) { impl->t = std::set<T> (s.impl->t); }
-  set (reader &r);
+  set () { }
+  set (unsigned dummy_size) : base(dummy_size) { }
+  set (const set &m) : base(m) { }
+  set (copy, const set &m) : base(COPY, m) { }
+  set (reader &r) : base(r) { }
   ~set () { }
   
-  set &operator = (const set &s) { impl = s.impl; return *this; }
-  
-  T pop ()
-  {
-    typename std::set<T>::const_iterator i = impl->t.begin ();
-    assert (i != impl->t.end ());
-    T tmp = *i;
-    impl->t.erase (i);
-    return tmp;
-  }
-  
-  bool is_empty () const { return impl->t.empty (); }
-  unsigned card () const { return impl->t.size (); }
-  const T &head () const;
-  const T &tail () const;
+  set &operator = (const set &m) { base::operator = (m); return *this; }
   
   bool operator == (const set &s) const;
+  bool operator != (const set &s) const { return !operator == (s); }
   bool operator < (const set &s) const;
   bool operator <= (const set &s) const;
   
-  void clear () { impl->t.clear (); }
-  void push (const T &v) { assert (impl->t.find (v) == impl->t.end ()); impl->t.insert (v); }
   bool toggle (const T &v);
-  void operator += (const T &v) { impl->t.insert (v); }
-  void operator -= (const T &v) { impl->t.erase (v); }
-  void yank (const T &v) { assert (operator % (v)); impl->t.erase (v); }
   
-  set &operator |= (const set<T> &s);
-  set &operator &= (const set<T> &s);
-  set &operator ^= (const set<T> &s);
-  
-  bool operator % (const T &v) const { return impl->t.find (v) != impl->t.end (); }
-  bool operator () (const T &v) const { return operator % (v); }
-  
-  void write_self (writer &w) const;
-  hash_t hash_self () const;
+  set &operator |= (const set &s);
+  set &operator &= (const set &s);
+  set &operator ^= (const set &s);
 };
 
 template<class T>
-class set_iter
-{
-  set<T> &s;
-  typename std::set<T>::const_iterator i, end;
-  bool deleted;
-  
- public:
-  set_iter (set<T> &s_) : s(s_), i(s_.impl->t.begin ()), end(s_.impl->t.end ()), deleted(0) { }
-  ~set_iter () { }
-  
-  void del ()
-  {
-    assert (!deleted);
-    typename std::set<T>::const_iterator iprev = i ++;
-    s.impl->t.erase (iprev);
-    deleted = 1;
-  }
-  
-  const T &val () const { assert (!deleted); return *i; }
-  operator bool () const { assert (!deleted); return i != end; }
-  set_iter &operator ++ () { if (deleted) deleted = 0; else i ++; return *this; }
-  void operator ++ (int) { operator ++ (); }
-};
+  using set_iter = set_wrapper_iter<std::set<T>, T>;
 
 template<class T>
-class set_const_iter
-{
-  typename std::set<T>::const_iterator i, end;
-  
- public:
-  set_const_iter (const set<T> &s) : i(s.impl->t.begin ()), end(s.impl->t.end ()) { }
-  ~set_const_iter () { }
-  
-  const T &val () const { return *i; }
-  operator bool () const { return i != end; }
-  set_const_iter &operator ++ () { i ++; return *this; }
-  void operator ++ (int) { i ++; }
-};
-
-template<class T>
-set<T>::set (reader &r)
-{
-  unsigned n;
-  read (r, n);
-  T x;
-  for (unsigned i = 0; i < n; i ++)
-    {
-      read (r, x);
-      push (x);
-    }
-}
-
-template<class T> const T &
-set<T>::head () const
-{
-  set_const_iter<T> i = *this;
-  assert (i);
-  return i.val ();
-}
-
-template<class T> const T &
-set<T>::tail () const
-{
-  typename std::set<T>::const_reverse_iterator i = impl->t.rbegin ();
-  assert (i != impl->t.rend ());
-  return *i;
-}  
+  using set_const_iter = set_wrapper_const_iter<std::set<T>, T>;
 
 template<class T> bool
 set<T>::operator == (const set &s) const
 {
-  typename std::set<T>::const_iterator i = impl->t.begin (),
+  typename std::set<T>::const_iterator i = this->impl->t.begin (),
     j = s.impl->t.begin (),
-    iend = impl->t.end (),
+    iend = this->impl->t.end (),
     jend = s.impl->t.end ();
   
   while (i != iend && j != jend && *i == *j)
@@ -155,9 +51,9 @@ set<T>::operator == (const set &s) const
 template<class T> bool
 set<T>::operator < (const set &s) const
 {
-  typename std::set<T>::const_iterator i = impl->t.begin (),
+  typename std::set<T>::const_iterator i = this->impl->t.begin (),
     j = s.impl->t.begin (),
-    iend = impl->t.end (),
+    iend = this->impl->t.end (),
     jend = s.impl->t.end ();
   
   while (i != iend && j != jend && *i == *j)
@@ -169,9 +65,9 @@ set<T>::operator < (const set &s) const
 template<class T> bool
 set<T>::operator <= (const set &s) const
 {
-  typename std::set<T>::const_iterator i = impl->t.begin (),
+  typename std::set<T>::const_iterator i = this->impl->t.begin (),
     j = s.impl->t.begin (),
-    iend = impl->t.end (),
+    iend = this->impl->t.end (),
     jend = s.impl->t.end ();
   
   while (i != iend && j != jend && *i == *j)
@@ -184,56 +80,56 @@ set<T>::operator <= (const set &s) const
 template<class T> bool
 set<T>::toggle (const T &v)
 {
-  typename std::set<T>::const_iterator i = impl->t.lower_bound (v),
-    end = impl->t.end ();
+  typename std::set<T>::const_iterator i = this->impl->t.lower_bound (v),
+    end = this->impl->t.end ();
   if (i != end && *i == v)
     {
-      impl->t.erase (i);
+      this->impl->t.erase (i);
       return 1;
     }
   else
     {
-      if (i == impl->t.begin ())
-	impl->t.insert (v);
+      if (i == this->impl->t.begin ())
+	this->impl->t.insert (v);
       else
 	{
 	  typename std::set<T>::const_iterator j = --i;
 	  assert (*j < v);
-	  impl->t.insert (j, v);
+	  this->impl->t.insert (j, v);
 	}
       return 0;
     }
 }
 
 template<class T> set<T> &
-set<T>::operator |= (const set<T> &s)
+set<T>::operator |= (const set &s)
 {
   std::set<T> news;
-  std::set_union (impl->t.begin (), impl->t.end (),
+  std::set_union (this->impl->t.begin (), this->impl->t.end (),
 		  s.impl->t.begin (), s.impl->t.end (),
 		  inserter (news, news.begin ()));
-  impl->t = news;
+  this->impl->t = news;
   return *this;
 }
 
 template<class T> set<T> &
-set<T>::operator &= (const set<T> &s)
+set<T>::operator &= (const set &s)
 {
   std::set<T> news;
-  std::set_intersection (impl->t.begin (), impl->t.end (),
+  std::set_intersection (this->impl->t.begin (), this->impl->t.end (),
 			 s.impl->t.begin (), s.impl->t.end (),
 			 inserter (news, news.begin ()));
-  impl->t = news;
+  this->impl->t = news;
   return *this;
 }
 
 template<class T> set<T> &
-set<T>::operator ^= (const set<T> &s)
+set<T>::operator ^= (const set &s)
 {
 #if 0
   printf ("before:\n");
   printf ("this:");
-  for (typename std::set<T>::const_iterator i = impl->t.begin (); i != impl->t.end (); i ++)
+  for (typename std::set<T>::const_iterator i = this->impl->t.begin (); i != this->impl->t.end (); i ++)
     printf (" %d", *i);
   printf ("\n");
   printf ("s:");
@@ -243,25 +139,25 @@ set<T>::operator ^= (const set<T> &s)
 #endif
   
 #if 1
-  typename std::set<T>::const_iterator i = impl->t.begin (),
-    iend = impl->t.end (),
+  typename std::set<T>::const_iterator i = this->impl->t.begin (),
+    iend = this->impl->t.end (),
     j = s.impl->t.begin (),
     jend = s.impl->t.end ();
   
   while (i != iend && j != jend && *i == *j)
     {
-      impl->t.erase (i);
-      i = impl->t.begin ();
+      this->impl->t.erase (i);
+      i = this->impl->t.begin ();
       j ++;
     }
   
   if (i == iend)
-    impl->t.insert (j, jend);
+    this->impl->t.insert (j, jend);
   else
     {
       if (j != jend && *j < *i)
 	{
-	  i = impl->t.insert (*j).first;
+	  i = this->impl->t.insert (*j).first;
 	  j ++;
 	}
       
@@ -274,14 +170,14 @@ set<T>::operator ^= (const set<T> &s)
 	    {
 	      while (j != jend)
 		{
-		  iprev = impl->t.insert (iprev, *j);
+		  iprev = this->impl->t.insert (iprev, *j);
 		  j ++;
 		}
 	      break;
 	    }
 	  else if (*i == *j)
 	    {
-	      impl->t.erase (i);
+	      this->impl->t.erase (i);
 	      i = iprev;
 	      j ++;
 	    }
@@ -289,12 +185,12 @@ set<T>::operator ^= (const set<T> &s)
 	    {
 	      while (j != jend && *j < *i)
 		{
-		  iprev = impl->t.insert (iprev, *j);
+		  iprev = this->impl->t.insert (iprev, *j);
 		  j ++;
 		}
 	      if (j != jend && *j == *i)
 		{
-		  impl->t.erase (i);
+		  this->impl->t.erase (i);
 		  i = iprev;
 		  j ++;
 		}
@@ -305,37 +201,19 @@ set<T>::operator ^= (const set<T> &s)
   
 #if 0
   std::set<T> news;
-  std::set_symmetric_difference (impl->t.begin (), impl->t.end (),
+  std::set_symmetric_difference (this->impl->t.begin (), this->impl->t.end (),
 				 s.impl->t.begin (), s.impl->t.end (),
 				 inserter (news, news.begin ()));
-  impl->t = news;
+  this->impl->t = news;
 #endif
   
 #if 0
   printf ("after:\n");
   printf ("this:");
-  for (typename std::set<T>::const_iterator i = impl->t.begin (); i != impl->t.end (); i ++)
+  for (typename std::set<T>::const_iterator i = this->impl->t.begin (); i != this->impl->t.end (); i ++)
     printf (" %d", *i);
   printf ("\n");
 #endif
   
   return *this;
 }
-
-template<class T> void
-set<T>::write_self (writer &w) const
-{
-  write (w, card ());
-  for (const_iter i = *this; i; i ++)
-    write (w, i.val ());
-}
-
-template<class T> hash_t
-set<T>::hash_self () const
-{
-  hash_t h = hash (card ());
-  unsigned j = 0;
-  for (const_iter i = *this; i && j < 10; i ++, j ++)
-    h = hash_combine (h, hash (i.val ()));
-  return h;
-}
diff --git a/lib/set_wrapper.h b/lib/set_wrapper.h
new file mode 100644
index 0000000..105c6bc
--- /dev/null
+++ b/lib/set_wrapper.h
@@ -0,0 +1,149 @@
+
+/* wrapper for STL maps */
+
+template<class S, class T> class set_wrapper_iter;
+template<class S, class T> class set_wrapper_const_iter;
+
+template<class S, class T>
+class set_wrapper
+{
+  friend class set_wrapper_iter<S, T>;
+  friend class set_wrapper_const_iter<S, T>;
+  
+ protected:
+  class set_wrapper_impl : public refcounted
+  {
+  public:
+    S t;
+  };
+  
+  ptr<set_wrapper_impl> impl;
+  
+ public:
+  typedef set_wrapper_iter<S, T> iter;
+  typedef set_wrapper_const_iter<S, T> const_iter;
+  
+ public:
+  set_wrapper () : impl(new set_wrapper_impl) { }
+  set_wrapper (unsigned dummy_size) : impl(new set_wrapper_impl) { }
+  set_wrapper (const set_wrapper &s) : impl(s.impl) { }
+  set_wrapper (copy, const set_wrapper &s) : impl(new set_wrapper_impl) { impl->t = S (s.impl->t); }
+  set_wrapper (reader &r);
+  ~set_wrapper () { }
+  
+  set_wrapper &operator = (const set_wrapper &s) { impl = s.impl; return *this; }
+  
+  T pop ()
+  {
+    typename S::const_iterator i = impl->t.begin ();
+    assert (i != impl->t.end ());
+    T tmp = *i;
+    impl->t.erase (i);
+    return tmp;
+  }
+  
+  bool is_empty () const { return impl->t.empty (); }
+  unsigned card () const { return impl->t.size (); }
+  const T &head () const;
+  const T &tail () const;
+  
+  void clear () { impl->t.clear (); }
+  void push (const T &v) { assert (impl->t.find (v) == impl->t.end ()); impl->t.insert (v); }
+  void operator += (const T &v) { impl->t.insert (v); }
+  void operator -= (const T &v) { impl->t.erase (v); }
+  void yank (const T &v) { assert (operator % (v)); impl->t.erase (v); }
+  
+  bool operator % (const T &v) const { return impl->t.find (v) != impl->t.end (); }
+  bool operator () (const T &v) const { return operator % (v); }
+  
+  void write_self (writer &w) const;
+  hash_t hash_self () const;
+};
+
+template<class S, class T>
+class set_wrapper_iter
+{
+  ptr<typename set_wrapper<S, T>::set_wrapper_impl> impl;
+  typename S::const_iterator i, end;
+  bool deleted;
+  
+ public:
+  set_wrapper_iter (set_wrapper<S, T> &s) : impl(s.impl), i(s.impl->t.begin ()), end(s.impl->t.end ()), deleted(0) { }
+  ~set_wrapper_iter () { }
+  
+  void del ()
+  {
+    assert (!deleted);
+    typename S::const_iterator iprev = i ++;
+    impl->t.erase (iprev);
+    deleted = 1;
+  }
+  
+  const T &val () const { assert (!deleted); return *i; }
+  operator bool () const { assert (!deleted); return i != end; }
+  set_wrapper_iter &operator ++ () { if (deleted) deleted = 0; else i ++; return *this; }
+  void operator ++ (int) { operator ++ (); }
+};
+
+template<class S, class T>
+class set_wrapper_const_iter
+{
+  ptr<const typename set_wrapper<S, T>::set_wrapper_impl> impl;
+  typename S::const_iterator i, end;
+  
+ public:
+  set_wrapper_const_iter (const set_wrapper<S, T> &s) : impl(s.impl), i(s.impl->t.begin ()), end(s.impl->t.end ()) { }
+  ~set_wrapper_const_iter () { }
+  
+  const T &val () const { return *i; }
+  operator bool () const { return i != end; }
+  set_wrapper_const_iter &operator ++ () { i ++; return *this; }
+  void operator ++ (int) { i ++; }
+};
+
+template<class S, class T>
+set_wrapper<S, T>::set_wrapper (reader &r)
+{
+  unsigned n;
+  read (r, n);
+  T x;
+  for (unsigned i = 0; i < n; i ++)
+    {
+      read (r, x);
+      push (x);
+    }
+}
+
+template<class S, class T> const T &
+set_wrapper<S, T>::head () const
+{
+  set_wrapper_const_iter<S, T> i = *this;
+  assert (i);
+  return i.val ();
+}
+
+template<class S, class T> const T &
+set_wrapper<S, T>::tail () const
+{
+  typename S::const_reverse_iterator i = impl->t.rbegin ();
+  assert (i != impl->t.rend ());
+  return *i;
+}  
+
+template<class S, class T> void
+set_wrapper<S, T>::write_self (writer &w) const
+{
+  write (w, card ());
+  for (const_iter i = *this; i; i ++)
+    write (w, i.val ());
+}
+
+template<class S, class T> hash_t
+set_wrapper<S, T>::hash_self () const
+{
+  hash_t h = hash (card ());
+  unsigned j = 0;
+  for (const_iter i = *this; i && j < 10; i ++, j ++)
+    h = hash_combine (h, hash (i.val ()));
+  return h;
+}
diff --git a/lib/vector.h b/lib/vector.h
index 52b1f15..06f4e52 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -36,8 +36,10 @@ public:
   
   vector &operator = (const vector &v) { unref (); ref (v.d); return *this; }
   
-  bool operator < (const vector &v) const;
   bool operator == (const vector &v) const;
+  bool operator != (const vector &v) const { return !operator == (v); }
+  
+  bool operator < (const vector &v) const;
   
   bool operator % (unsigned i) { return d && i < d->n; }
   
diff --git a/main.cpp b/main.cpp
index 4ed6dc1..2ae5f6e 100644
--- a/main.cpp
+++ b/main.cpp
@@ -60,14 +60,14 @@ test_ring (int p)
 	  
 	  if (x != 0 && x | y)
 	    {
-	      R q = y.div (x);
+	      R q = y.divide_exact (x);
 	      assert (y == q * x);
 	    }
 	  
 	  if (x != 0 || y != 0)
 	    {
-	      triple<R, R, R> t = x.extended_gcd (y);
-	      assert (t.first == t.second*x + t.third*y);
+	      tuple<R, R, R> t = x.extended_gcd (y);
+	      assert (get<0> (t) == get<1> (t)*x + get<2> (t)*y);
 	    }
 	  
 	  for (int k = -n; k <= n; k ++)
@@ -127,6 +127,97 @@ rank_lte (multivariate_laurentpoly<Z> p,
 int
 main ()
 {
+  basedvector<int, 1> v (4);
+  v[1] = 2;
+  v[2] = -1;
+  v[3] = 0xbf32813;
+  v[4] = -352182;
+
+  vector<int> sleb (512);
+  for (unsigned i = 0; i < 512; i ++)
+    sleb[i] = ((int)i << 23);
+
+  printf ("sleb:\n");
+  for (unsigned i = 0; i < 512; i ++)
+    printf ("  %d\n", sleb[i]);
+  
+  vector<unsigned> uleb (512);
+  for (unsigned i = 0; i < 512; i ++)
+    uleb[i] = ((unsigned)i << 23);
+  
+  printf ("uleb:\n");
+  for (unsigned i = 0; i < 512; i ++)
+    printf ("  %u\n", uleb[i]);
+  
+  polynomial<Z> one (1);
+  polynomial<Z> x (1, 1);
+  polynomial<Z> bigp (Z (367521));
+  
+  polynomial<Z> t = x + bigp*one;
+  polynomial<Z> p = (t*t*t + t*t + t + one)*(bigp*t*t + bigp*bigp*t + bigp*bigp*bigp)*(bigp*bigp*t*t*t + bigp*t*t + t + one);
+  display ("p: ", p);
+  
+  {
+    file_writer w ("test.dat");
+    write (w, v);
+    write (w, p);
+    write (w, sleb);
+    write (w, uleb);
+  }
+  
+  {
+    gzfile_writer w ("test.2.dat.gz");
+    write (w, v);
+    write (w, p);
+    write (w, sleb);
+    write (w, uleb);
+  }
+  
+  system ("gzip -cd test.2.dat.gz > test.2.dat");
+  
+  {
+    file_reader r ("test.dat");
+    basedvector<int, 1> v2 (r);
+    polynomial<Z> p2 (r);
+    
+    vector<int> sleb2 (r);
+    printf ("sleb2:\n");
+    for (unsigned i = 0; i < 512; i ++)
+      printf ("  %d\n", sleb2[i]);
+    
+    vector<unsigned> uleb2 (r);
+    assert (v == v2);
+    assert (p == p2);
+    assert (sleb == sleb2);
+    assert (uleb == uleb2);
+  }
+  
+  system ("gzip -c9 test.dat > test.dat.gz");
+  
+  {
+    gzfile_reader r ("test.dat.gz");
+    basedvector<int, 1> v2 (r);
+    polynomial<Z> p2 (r);
+    vector<int> sleb2 (r);
+    vector<unsigned> uleb2 (r);
+    assert (v == v2);
+    assert (p == p2);
+    assert (sleb == sleb2);
+    assert (uleb == uleb2);
+  }
+  
+#if 0
+  
+#if 0
+  test_ring<Z> (0);
+  test_ring<polynomial<Z> > (0);
+  
+  test_field<Z2> ();
+  test_field<Zp<7> > ();
+  test_field<Q> ();
+  // test_field<fraction_field<Zp<7> > > ();
+#endif
+  
 #if 0
   for (unsigned i = 1; i <= 10; i ++)
     for (unsigned j = 1; j <= rolfsen_crossing_knots (i); j ++)
@@ -178,4 +269,5 @@ main ()
 	if (! rank_lte (E3_p, tt_p))
 	  printf (" > rank E2 > rank tt!!\n");
       }
+#endif
 }
diff --git a/spanning_tree_complex.h b/spanning_tree_complex.h
index 5a7a091..82abe9d 100644
--- a/spanning_tree_complex.h
+++ b/spanning_tree_complex.h
@@ -39,7 +39,7 @@ public:
   tree_generators (const spanning_tree_complex<F> &c_) : c(c_) { }
   ~tree_generators () { }
   
-  tree_generators &operator = (const tree_generators &); // doesn't exist
+  tree_generators &operator = (const tree_generators &) = delete;
   
   unsigned dim () const { return c.trees.size (); }
   unsigned free_rank () const { return c.trees.size (); }
diff --git a/sseq.h b/sseq.h
index 29fa5a9..d2f5420 100644
--- a/sseq.h
+++ b/sseq.h
@@ -34,6 +34,7 @@ class sseq_bounds
     && minq == b.minq
     && maxq == b.maxq;
   }
+  bool operator != (const sseq_bounds &b) const { return !operator == (b); }
 };
 
 class sseq_page
@@ -52,6 +53,8 @@ class sseq_page
   ~sseq_page () { }
   
   bool operator == (const sseq_page &pg) const { return rank == pg.rank && im_rank == pg.im_rank; }
+  bool operator != (const sseq_page &pg) const { return !operator == (pg); }
+  
   bool equal_as_spaces (const sseq_page &pg) const { return rank == pg.rank; }
   
   unsigned total_rank () const;
@@ -95,6 +98,9 @@ class sseq
   unsigned homological_width () const { return pages[1].homological_width (bounds); }
   
   bool operator == (const sseq &ss) const { return bounds == ss.bounds && pages == ss.pages; }
+  bool operator != (const sseq &ss) const { return !operator == (ss); }
+
+  
   bool equal_as_spaces (const sseq &ss) const;
   
   void texshow (FILE *fp, std::string name);
@@ -121,7 +127,7 @@ public:
   { }
   ~simplified_complex_generators () { }
   
-  simplified_complex_generators &operator = (const simplified_complex_generators &); // doesn't exist
+  simplified_complex_generators &operator = (const simplified_complex_generators &) = delete;
   
   unsigned dim () const { return new_n; }
   unsigned free_rank () const { return new_n; }
diff --git a/todo.txt b/todo.txt
index ad38695..9b28950 100644
--- a/todo.txt
+++ b/todo.txt
@@ -13,10 +13,8 @@ in knotkit/
    deepcopy
 
 in lib/
- - add hashset
  - between interoperation between sets that require bounds (bitset,
    ullmanset) and those that don't (set, hashset)
- - add make_pair, etc.
 
 in algebra/
  - linear_combnation can be more efficient (eg no searching when
@@ -25,3 +23,63 @@ in algebra/
  - cleaner interface for gcd, etc. in fraction_field
  - support kernel, cokernel of torsion modules
  - change gcd et al so gcd is always positive
+ - split linear_combination into internal and external versions
+
+git:
+ - fix bug in Josh's version
+ - make this master, merge into Josh's version
+ - merge in square to master => then write paper!
+ - merge in other improvements (sped up linear algebra?  multivariate (laurent) polynomial?)
+ - should exponents on polynomials be Z or unsigned?
+
+general:
+ - unify simplify_chain_complex, sseq, kernel/image/homology (as much as possible)
+ - sseq over Z
+ - quadratic linear solver (enumerator?)
+ - classes should either COPY or fully SHARE, but not some of each
+   (e.g. knot_diagram, resolution_diagram, etc.)  -- this is the right
+   semantics.
+ - clean up/delete setcommon, mapcommon
+
+c++11:
+ - make vector wrap std::vector (probably faster)
+ 
+ - use tuple instead of pair, triple; and use tie, a great idiom
+ - standardize ring, field interface so test_ring works for everything
+ - get clear on default semantics: prefer default when possible
+ 
+ - standardize use of {} for initializers
+ - initializer_list for container types
+ - support for iterator
+
+notes:
+
+all rings (R): multivariate_laurentpoly<R>
+
+fields (F): Z2 Zp Q fraction_field<E> (for just pids?)
+
+euclidean domains (E): Z polynomial<F>
+
+euclidean domain interface:
+
+num/demon : denom must be a unit
+ = num.operator / (denom)
+
+denom.divides (num)
+denom | num
+
+  num.divide_exact (denom)
+
+(q, r) = num.divide_with_remainder (denom)
+  (q, r) = num / denom
+  (q, r) = num.operator / (denom)
+satisfies r = 0 or f(r) < f(denom) and num = q*denom + r
+
+ gcd (a, b)
+ lcm (a, b)
+ r = num.mod (denom)
+satisfies r = 0 or f(r) < f(denom)
+
+
+(d, s, t) = a.extended_gcd (b)
+ satifies d is the gcd and d = a*s + b*t