diff --git a/Makefile b/Makefile
index 8340db3..18109a1 100644
--- a/Makefile
+++ b/Makefile
@@ -5,14 +5,14 @@ BISON = bison
 FLEX = flex
 
 CXX = g++ -std=c++11
-# CXX = OMPI_CXX=clang++ mpic++ -fno-color-diagnostics --stdlib=libc++ --std=c++11 -I/u/cseed/llvm-3.1/lib/c++/v1
+#CXX = OMPI_CXX=clang++ mpic++ -fno-color-diagnostics --stdlib=libc++ --std=c++11 -I/u/cseed/llvm-3.1/lib/c++/v1
 #CXX = clang++ -fno-color-diagnostics --stdlib=libc++ --std=c++11
 
 INCLUDES = -I. -I/opt/local/include 
 
 # OPTFLAGS = -g
 OPTFLAGS = -O2 -g
-# OPTFLAGS = -O2 -DNDEBUG
+# OPTFLAGS = -O2 -g -DNDEBUG
 
 LDFLAGS = -L/opt/local/lib
 
@@ -44,7 +44,7 @@ KNOTKIT_HEADERS = knotkit.h planar_diagram.h dt_code.h knot_diagram.h \
   smoothing.h cobordism.h cube.h steenrod_square.h \
   spanning_tree_complex.h cube_impl.h sseq.h simplify_chain_complex.h
 
-LIBS = -lgmp -lz
+LIBS = -lgmpxx -lgmp  -lz
 
 all: kk
 
diff --git a/algebra/Q.h b/algebra/Q.h
index 21a1d06..24c8509 100644
--- a/algebra/Q.h
+++ b/algebra/Q.h
@@ -1,70 +1,154 @@
+#ifndef _KNOTKIT_ALGEBRA_Q_H
+#define _KNOTKIT_ALGEBRA_Q_H
+#include <gmpxx.h>
+#include <iostream>
+#include <memory>
+#include <tuple>
+#include <cassert>
 
 class Q
 {
  public:
-  typedef ::linear_combination<Q> linear_combination;
-  typedef ::linear_combination_const_iter<Q> linear_combination_const_iter;
+  using linear_combination = ::linear_combination<Q>;
+  using linear_combination_const_iter = ::linear_combination_const_iter<Q>;
 
  private:
-  enum steal { STEAL };
+  // enum steal { STEAL };
+  // enum copy { COPY };
+  // class Q_impl : public refcounted
+  // {
+  // public:
+  //   mpq_t x;
+    
+  // public:
+  //   Q_impl () { mpq_init (x); }
+  //   Q_impl (int init)
+  //   {
+  //     mpq_init (x);
+  //     mpq_set_si (x, init, 1);
+  //   }
+    
+  //   Q_impl (copy, mpq_srcptr init)
+  //   {
+  //     mpq_init (x);
+  //     mpq_set (x, init);
+  //   }
+    
+  //   Q_impl (steal, mpq_srcptr init) { x[0] = *init; }
+  //   Q_impl (reader &r)
+  //   {
+  //     mpq_init (x);
+  //     r.read_mpz (mpq_numref (x));
+  //     r.read_mpz (mpq_denref (x));
+  //   }
+    
+  //   ~Q_impl () { mpq_clear (x); }
+    
+  //   void write_self (writer &w) const
+  //   {
+  //     w.write_mpz (mpq_numref (x));
+  //     w.write_mpz (mpq_denref (x));
+  //   }
+  // };
   
-  class Q_impl : public refcounted
-  {
-  public:
-    mpq_t x;
-    
-  public:
-    Q_impl () { mpq_init (x); }
-    Q_impl (int init)
-    {
-      mpq_init (x);
-      mpq_set_si (x, init, 1);
-    }
-    
-    Q_impl (copy, mpq_srcptr init)
-    {
-      mpq_init (x);
-      mpq_set (x, init);
-    }
-    
-    Q_impl (steal, mpq_srcptr init) { x[0] = *init; }
-    Q_impl (reader &r)
-    {
-      mpq_init (x);
-      r.read_mpz (mpq_numref (x));
-      r.read_mpz (mpq_denref (x));
-    }
-    
-    ~Q_impl () { mpq_clear (x); }
-    
-    void write_self (writer &w) const
-    {
-      w.write_mpz (mpq_numref (x));
-      w.write_mpz (mpq_denref (x));
-    }
-  };
+  // ptr<Q_impl> impl;
+  std::shared_ptr<mpq_class> impl;
+    void write_state() const {
+    std::cout << "I store the following value " << *this << "\n";
+    std::cout << "Number of objects pointing to the same value " << impl.use_count() << "\n";
+    /* std::cout << "I point to " << impl.get() << "\n"; */
+    /* std::cout << "My size " << sizeof(*this) << "\n"; */
+    /* std::cout << "Size of std::shared_ptr<mpz_class> " << sizeof(impl) << "\n"; */
+    /* std::cout << "Size of mpz_class " << sizeof(*impl) << "\n"; */
+  }
   
-  ptr<Q_impl> impl;
-  
-  Q (steal, mpq_srcptr init) : impl(new Q_impl (STEAL, init)) { }
-
  public:
-  Q () : impl(new Q_impl) { }
-  Q (int init) : impl(new Q_impl (init)) { }
-  Q (const Q &q) : impl(q.impl) { }
-  Q (copy, const Q &q) : impl(new Q_impl (COPY, q.impl->x)) { }
-  Q (reader &r) : impl(new Q_impl (r)) { }
-  ~Q () { }
+  Q() : impl(new mpq_class()) {
+#ifdef DEBUG_Q
+    std::cout << "Q()" << "\n";
+    write_state();
+#endif
+  }
+  Q(mpq_t q) : impl(new mpq_class(q)) {
+#ifdef DEBUG_Q
+    std::cout << "Q(mpq_t)" << "\n";
+    write_state();
+#endif
+  }
+  Q(mpq_class q) : impl(new mpq_class(q)) {
+#ifdef DEBUG_Q
+    std::cout << "Q(mpq_class)" << "\n";
+    write_state();
+#endif
+  }
+  Q(int init) : impl(new mpq_class (init)) {
+#ifdef DEBUG_Q
+    std::cout << "Q(int)" << "\n";
+    write_state();
+#endif
+  }
+  Q(const Q& q) : impl(q.impl) {
+#ifdef DEBUG_Q
+    std::cout << "Q(const Q&)" << "\n";
+    write_state();
+#endif
+  }
+  Q(Q&& q) : impl(std::move (q.impl)) {
+    q.impl = std::make_shared<mpq_class>(0);
+#ifdef DEBUG_Q
+    std::cout << "Q(Q&& q)" << "\n";
+    write_state();
+#endif
+  }
+  Q(copy, const Q &q) : impl(new mpq_class (*q.impl.get())) {
+#ifdef DEBUG_Q
+    std::cout << "Q(copy, const Q&)" << "\n";
+    write_state();
+#endif
+  }
+  //Q(reader &r) : impl(new Q_impl (r)) { }
+  ~Q() {
+#ifdef DEBUG_Q
+    std::cout << "~Q()" << "\n";
+    write_state();
+    if(get_count() == 1)
+      std::cout << "Destroying..." << "\n";
+#endif
+  }
   
-  Q &operator = (const Q &q) { impl = q.impl; return *this; }
-  Q &operator = (int x) { impl = new Q_impl (x); return *this; }
+  Q &operator = (const Q &q) {
+    impl = q.impl;
+    return *this;
+  }
+  Q& operator = (Q&& q) {
+    impl = std::move(q.impl);
+    q.impl = std::make_shared<mpq_class>(0);
+    return *this;
+  }
+  Q &operator = (int x) {
+    impl = std::make_shared<mpq_class>(0);
+    return *this;
+  }
   
-  bool operator == (const Q &q) const { return mpq_cmp (impl->x,q.impl->x) == 0; }
+  bool operator == (const Q &q) const {
+    return *impl.get() == *q.impl.get();
+  }
+  bool operator == (const int r) const {
+    return *impl.get() == mpq_class(r);
+  }
+  bool operator != (const Q& q) const {
+    return ! operator == (q);
+  }
+  bool operator != (const int r) const {
+    return !operator == (r);
+  }
   
-  bool operator == (int r) const { return mpq_cmp_si (impl->x, r, 1) == 0; }
-  bool operator != (int r) const { return !operator == (r); }
-  
-  bool operator < (const Q &q) const { return mpq_cmp (impl->x, q.impl->x) < 0; }
+  bool operator < (const Q &q) const {
+    return *impl.get() < *q.impl.get();
+  }
+  // bool operator > (const Q& q) const {
+  //   return *impl.get() > *q.impl.get();
+  // }
   
   bool is_unit () const
   {
@@ -73,87 +157,63 @@ class Q
   
   Q operator - () const
   {
-    mpq_t x;
-    mpq_init (x);
-    mpq_neg (x, impl->x);
-    return Q (STEAL, x);
+    return Q(-*impl.get());
   }
 
   Q recip () const
   {
-    mpq_t x;
-    mpq_init (x);
-    mpq_inv (x, impl->x);
-    return Q (STEAL, x);
+    mpq_t q;
+    mpq_init(q);
+    mpq_inv(q, impl.get()->get_mpq_t());
+    return Q(q);
   }
   
-  Q operator + (const Q &q) const
+  Q operator + (const Q& q) const
   {
-    mpq_t x;
-    mpq_init (x);
-    mpq_add (x, impl->x, q.impl->x);
-    return Q (STEAL, x);
+    return Q(*impl.get() + *q.impl.get());
   }
   
-  Q operator - (const Q &q) const
+  Q operator - (const Q& q) const
   {
-    mpq_t x;
-    mpq_init (x);
-    mpq_sub (x, impl->x, q.impl->x);
-    return Q (STEAL, x);
+    return Q(*impl.get() - *q.impl.get());
   }
   
-  Q operator * (const Q &q) const
+  Q operator * (const Q& q) const
   {
-    mpq_t x;
-    mpq_init (x);
-    mpq_mul (x, impl->x, q.impl->x);
-    return Q (STEAL, x);
+    return Q(*impl.get() * (*q.impl.get()));
   }
   
-  Q operator / (const Q &q) const
+  Q operator / (const Q& q) const
   {
-    assert (q != 0);
-    
-    mpq_t x;
-    mpq_init (x);
-    mpq_div (x, impl->x, q.impl->x);
-    return Q (STEAL, x);
+    return Q(*impl.get() / *q.impl.get());
   }
   
-  Q &muladdeq (const Q &q1, const Q &q2)
+  Q &muladdeq (const Q& q1, const Q& q2)
   {
-    // ??? do inline saves refcount overhead
     return operator += (q1 * q2);
   }
   
-  Q &operator += (const Q &q)
+  Q &operator += (const Q& q)
   {
-    mpq_add (impl->x, impl->x, q.impl->x);
-    return *this;
+    return *this = *this + q;
   }
   
-  Q &operator -= (const Q &q)
+  Q &operator -= (const Q& q)
   {
-    mpq_sub (impl->x, impl->x, q.impl->x);
-    return *this;
+    return *this = *this - q;
   }
   
   Q &operator *= (const Q &q)
   {
-    mpq_mul (impl->x, impl->x, q.impl->x);
-    return *this;
+    return *this = *this * q;
   }
   
   Q &operator /= (const Q &q)
   {
     assert (q != 0);
-    
-    mpq_div (impl->x, impl->x, q.impl->x);
-    return *this;
+    return *this = *this / q;
   }
   
-  // d | n, d.divides (n)
   bool divides (const Q &num) const
   {
     return *this != 0 || num == 0;
@@ -163,12 +223,12 @@ class Q
   
   Q div (const Q &d) const { return operator / (d); }
   
-  tuple<Q, Q, Q> extended_gcd (const Q &q) const
+  std::tuple<Q, Q, Q> extended_gcd (const Q &q) const
   {
     if (*this != 0)
-      return tuple<Q, Q, Q> (*this, 1, 0);
+      return std::tuple<Q, Q, Q> (*this, 1, 0);
     else
-      return tuple<Q, Q, Q> (q, 0, 1);
+      return std::tuple<Q, Q, Q> (q, 0, 1);
   }
   
   Q gcd (const Q &q) const
@@ -176,9 +236,16 @@ class Q
     assert (*this != 0 || q != 0);
     return 1;
   }
-  
+  friend std::ostream& operator << (std::ostream& os, const Q& q) {
+    return os << *q.impl.get();
+  }
   static void show_ring () { printf ("Q"); }
-  void show_self () const { mpq_out_str (stdout, 10, impl->x); }
-  void display_self () const { show_self (); newline (); }
-  void write_self (writer &w) const { write (w, *impl); }
+  void show_self () const { std::cout << *this; }
+  void display_self () const { std::cout << *this << "\n"; }
+  // void write_self (writer &w) const { write (w, *impl); }
+  int get_count() const {
+    return impl.use_count();
+  }
 };
+
+#endif // _KNOTKIT_ALGEBRA_Q_H
diff --git a/algebra/Z.h b/algebra/Z.h
index 8478ebf..b15856a 100644
--- a/algebra/Z.h
+++ b/algebra/Z.h
@@ -1,209 +1,239 @@
+#ifndef _KNOTKIT_ALGEBRA_Z_H
+#define _KNOTKIT_ALGEBRA_Z_H
+#include<gmpxx.h>
+#include<tuple>
+#include<memory>
+#include<iostream>
+#include<cassert>
 
 class Z
 {
  public:
-  typedef ::linear_combination<Z> linear_combination;
-  typedef ::linear_combination_const_iter<Z> linear_combination_const_iter;
-  
- private:
-  enum steal { STEAL };
-  
-  class Z_impl : public refcounted
-  {
-  public:
-    mpz_t x;
-    
-  public:
-    Z_impl () { mpz_init (x); }
-    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);
-      r.read_mpz (x);
-    }
-    
-    ~Z_impl () { mpz_clear (x); }
-    
-    void write_self (writer &w) const
-    {
-      w.write_mpz (x);
-    }
-  };
-  
-  ptr<Z_impl> impl;
-  
-  Z (steal, mpz_srcptr init) : impl(new Z_impl (STEAL, init)) { }
+  using linear_combination = ::linear_combination<Z>;
+  using linear_combination_const_iter = ::linear_combination_const_iter<Z>;
+
+ private:  
+  std::shared_ptr<mpz_class> impl;
+  void write_state() const {
+    std::cout << "I store the following value " << *this << "\n";
+    std::cout << "Number of objects pointing to the same value " << impl.use_count() << "\n";
+    /* std::cout << "I point to " << impl.get() << "\n"; */
+    /* std::cout << "My size " << sizeof(*this) << "\n"; */
+    /* std::cout << "Size of std::shared_ptr<mpz_class> " << sizeof(impl) << "\n"; */
+    /* std::cout << "Size of mpz_class " << sizeof(*impl) << "\n"; */
+  }
   
  public:
-  Z () : impl(new Z_impl) { }
-  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; }
-  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); }
-  
-  bool operator < (const Z &z) const { return mpz_cmp (impl->x, z.impl->x) < 0; }
+  Z() : impl(new mpz_class) {
+#ifdef DEBUG_Z
+    std::cout << "Z()" << "\n";
+    write_state();
+#endif
+  }
+  Z(mpz_t z) : impl(new mpz_class(z)) {
+#ifdef DEBUG_Z
+    std::cout << "Z(mpz_t)" << "\n";
+    write_state();
+#endif
+  }
+  Z(mpz_class z) : impl(new mpz_class(z)) {
+#ifdef DEBUG_Z
+    std::cout << "Z(mpz_class)" << "\n";
+    write_state();
+#endif
+  }
+  Z(int init) : impl(new mpz_class(init)) {
+#ifdef DEBUG_Z
+    std::cout << "Z(int)" << "\n";
+    write_state();
+#endif
+  }
+  Z(const Z& z) : impl(z.impl) {
+#ifdef DEBUG_Z
+    std::cout << "Z(const Z& z)" << "\n";
+    write_state();
+#endif
+  }
+  Z(Z&& z) : impl(std::move(z.impl)) {
+#ifdef DEBUG_Z
+    std::cout << "Z(Z&& z)" << "\n";
+    write_state();
+#endif
+    z.impl = std::make_shared<mpz_class>(mpz_class(0));
+  }
+  ~Z() {
+#ifdef DEBUG_Z
+    std::cout << "~Z()" << "\n";
+    write_state();
+    if(impl.use_count() == 1)
+      std::cout << "Destroying..." << "\n";
+#endif
+  }
+
+  Z& operator = (const Z& z) {
+    impl = z.impl;
+#ifdef DEBUG_Z
+    std::cout << "Z& operator = (const Z&)" << "\n";
+    write_state();
+#endif
+    return *this;
+  }
   
+  Z& operator = (const int x) {
+    impl = std::make_shared<mpz_class>(mpz_class(x));
+#ifdef DEBUG_Z
+    std::cout << "Z& operator = (int)" << "\n";
+    write_state();
+#endif
+    return *this;
+  }
+  Z& operator = (Z&& z) {
+    impl = std::move(z.impl);
+    z.impl = std::make_shared<mpz_class>(mpz_class(0));
+#ifdef DEBUG_Z
+    std::cout << "Z& operator = (Z&&)" << "\n";
+    write_state();
+#endif
+    return *this;
+  }
+
+  bool operator == (const Z& z) const {
+    return *impl.get() == *z.impl.get();
+  }
+  bool operator != (const Z& z) const {
+    return !operator == (z);
+  }
+
+  bool operator == (const int y) const {
+    return *impl.get() == mpz_class(y);
+  }
+  bool operator != (const int y) const {
+    return ! operator == (y);
+  }
+
+  bool operator < (const Z& z) const {
+    return *impl.get() < *z.impl.get();
+  }
+
   bool is_unit () const
   {
     return *this == 1 || *this == -1;
   }
-  
-  Z operator + (const Z &z) const
-  {
-    mpz_t x;
-    mpz_init (x);
-    mpz_add (x, impl->x, z.impl->x);
-    return Z (STEAL, x);
-  }
-  
-  Z operator - () const
-  {
-    mpz_t x;
-    mpz_init (x);
-    mpz_neg (x, impl->x);
-    return Z (STEAL, x);
-  }
-  
-  Z operator - (const Z &z) const
-  {
-    mpz_t x;
-    mpz_init (x);
-    mpz_sub (x, impl->x, z.impl->x);
-    return Z (STEAL, x);
-  }
-  
-  Z operator * (const Z &z) const
-  {
-    mpz_t x;
-    mpz_init (x);
-    mpz_mul (x, impl->x, z.impl->x);
-    return Z (STEAL, x);
-  }
-  
-  Z operator / (const Z &denom) const
-  {
-    if (mpz_cmp_si (denom.impl->x, 1) == 0)
-      return *this;
-    else
-      {
-	assert (mpz_cmp_si (denom.impl->x, -1) == 0);
-	
-	mpz_t x;
-	mpz_init (x);
-	mpz_neg (x, impl->x);
-	return Z (STEAL, x);
-      }
-  }
-  
-  Z recip () const
-  {
-    assert (is_unit ());
-    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);
-    return *this;
-  }
-  
-  Z &operator += (const Z &z)
-  {
-    mpz_add (impl->x, impl->x, z.impl->x);
-    return *this;
-  }
-  
-  Z &operator -= (const Z &z)
-  {
-    mpz_sub (impl->x, impl->x, z.impl->x);
-    return *this;
-  }
-  
-  Z &operator *= (const Z &z)
-  {
-    mpz_mul (impl->x, impl->x, z.impl->x);
-    return *this;
-  }
-  
-  Z &operator /= (const Z &z)
-  {
-    if (mpz_cmp_si (z.impl->x, 1) == 0)
-      ;
-    else
-      {
-	assert (mpz_cmp_si (z.impl->x, -1) == 0);
-	mpz_neg (impl->x, impl->x);
-      }
-    return *this;
-  }
-  
-  bool divides (const Z &num) const
-  {
-    return mpz_divisible_p (num.impl->x, impl->x);
-  }
-  bool operator | (const Z &num) const { return divides (num); }
 
-  Z divide_exact (const Z &denom) const
-  {
-    // num = *this
+  Z operator + (const Z& z) const {
+    return Z(*impl.get() + *z.impl.get());
+  }
+
+  Z operator - () const {
+    return Z(-*impl.get());
+  }
+
+  Z operator - (const Z& z) const {
+    return Z(*impl.get() - *z.impl.get());
+  }
+
+  Z operator * (const Z& z) const {
+    return Z(*impl.get() * *z.impl.get());
+  }
+
+  Z operator / (const Z& z) const {
+    if(z == 0)
+      return *this;
+    else {
+      assert(z != 0);
+      return Z(*impl.get() / *z.impl.get());
+    }
+  }
+
+  Z recip () const {
+    assert(is_unit());
+    return *this;
+  }
+
+  Z& muladdeq(const Z& z1, const Z& z2) {
+    return *this += z1 * z2;
+  }
+
+  Z& operator += (const Z& z) {
+    *this = *this + z;
+    return *this;
+  }
+
+  Z& operator -= (const Z& z) {
+    *this = *this - z;
+    return *this;
+  }
+
+  Z& operator *= (const Z& z) {
+    *this = *this * z;
+    return *this;
+  }
+
+  Z& operator /= (const Z& z) {
+    *this = *this / z;
+    return *this;
+  }
+
+  bool divides(const Z& num) const {
+    return mpz_divisible_p(num.impl.get()->get_mpz_t(),impl.get()->get_mpz_t());
+  }
+
+  bool operator | (const Z& num) const {
+    return divides(num);
+  }
+
+  Z divide_exact(const Z& denom) const {
     mpz_t q;
-    mpz_init (q);
-    mpz_divexact (q, impl->x, denom.impl->x);
-    return Z (STEAL, q);
+    mpz_init(q);
+    mpz_divexact(q,impl.get()->get_mpz_t(), denom.impl.get()->get_mpz_t());
+    return Z(q);
   }
-  
-  tuple<Z, Z> divide_with_remainder (const Z &denom) const
-  {
-    // *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));
+
+  std::tuple<Z,Z> divide_with_remainder(const Z& denom) const {
+    mpz_t q,r;
+    mpz_init(q);
+    mpz_init(r);
+    mpz_tdiv_qr(q, r, impl.get()->get_mpz_t(), denom.impl.get()->get_mpz_t());
+    return std::make_tuple(Z(q), Z(r));
   }
-  
-  Z gcd (const Z &z) const
-  {
+
+  Z gcd (const Z& z) const {
     mpz_t d;
-    mpz_gcd (d, impl->x, z.impl->x);
-    return Z (STEAL, d);
+    mpz_init(d);
+    mpz_gcd(d, impl.get()->get_mpz_t(), z.impl.get()->get_mpz_t());
+    return Z(d);
   }
-  
-  Z lcm (const Z &z) const
-  {
+
+  Z lcm (const Z& z) const {
     mpz_t m;
-    mpz_lcm (m, impl->x, z.impl->x);
-    return Z (STEAL, m);
+    mpz_init(m);
+    mpz_lcm(m, impl.get()->get_mpz_t(), z.impl.get()->get_mpz_t());
+    return Z(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));
+
+  std::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.get()->get_mpz_t(), z.impl.get()->get_mpz_t());
+    return std::make_tuple<Z, Z, Z>(Z(d), Z(s), Z(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 (); }
-  void write_self (writer &w) const { write (w, *impl); }
+  void show_self () const {
+    std::cout << *this;
+  }
+  void display_self () const {
+    std::cout << *this << "\n";
+  }
+
+  friend std::ostream& operator << (std::ostream& os, const Z& z) {
+    return os << *z.impl;
+  }
+  int get_count() const {
+    return impl.use_count();
+  }
 };
+
+#endif // _KNOTKIT_ALGEBRA_Z_H
diff --git a/algebra/Z2.h b/algebra/Z2.h
index 09be75a..8cb5df5 100644
--- a/algebra/Z2.h
+++ b/algebra/Z2.h
@@ -1,9 +1,15 @@
+#ifndef KNOTKIT_ALGEBRA_Z2_H
+#define KNOTKIT_ALGEBRA_Z2_H
+#include <iostream>
+#include <tuple>
+#include <cassert>
+
 class Z2
 {
  public:
-  typedef ::linear_combination<Z2> linear_combination;
+  using linear_combination = ::linear_combination<Z2>;
   // typedef linear_combination_iter<Z2> linear_combination_iter;
-  typedef ::linear_combination_const_iter<Z2> linear_combination_const_iter;
+  using linear_combination_const_iter = ::linear_combination_const_iter<Z2>;
   
  private:
   bool v;
@@ -13,37 +19,40 @@ class Z2
   Z2 (int x) : v((bool)(x & 1)) { }
   Z2 (unsigned x) : v((bool)(x & 1)) { }
   Z2 (bool v_) : v(v_) { }
-  Z2 (const Z2 &x) : v(x.v) { }
-  Z2 (copy, const Z2 &x) : v(x.v) { }
+  Z2 (const Z2& x) : v(x.v) { }
+  Z2 (Z2&& x) : v(std::move(x.v)) {
+    x.v = 0;
+  }
   Z2 (reader &r) { v = r.read_bool (); }
   ~Z2 () { }
   
-  Z2 &operator = (const Z2 &x) { v = x.v; return *this; }
-  Z2 &operator = (int x) { v = (bool)(x & 1); return *this; }
+  Z2& operator = (const Z2& x) { v = x.v; return *this; }
+  Z2& operator = (Z2&& x) { v = x.v; x.v = 0; return *this; }
+  Z2& operator = (int x) { v = (bool)(x & 1); return *this; }
   
-  bool operator == (const Z2 &x) const { return v == x.v; }
-  bool operator != (const Z2 &x) const { return v != x.v; }
+  bool operator == (const Z2& x) const { return v == x.v; }
+  bool operator != (const Z2& x) const { return v != x.v; }
   
   bool operator == (int x) const { return v == (bool)(x & 1); }
   bool operator != (int x) const { return !operator == (x); }
   
   bool is_unit () const { return v; }
   
-  Z2 operator + (const Z2 &x) const { return Z2 (v ^ x.v); }
-  Z2 operator - (const Z2 &x) const { return Z2 (v ^ x.v); }
+  Z2 operator + (const Z2& x) const { return Z2 (v ^ x.v); }
+  Z2 operator - (const Z2& x) const { return Z2 (v ^ x.v); }
   Z2 operator - () const { return *this; }
-  Z2 operator * (const Z2 &x) const { return Z2 (v & x.v); }
-  Z2 operator / (const Z2 &x) const { assert (x.v); return *this; }
+  Z2 operator * (const Z2& x) const { return Z2 (v & x.v); }
+  Z2 operator / (const Z2& x) const { assert (x.v); return *this; }
   
   Z2 recip () const { assert (v); return *this; }
   
-  Z2 &operator += (const Z2 &x) { v ^= x.v; return *this; }
-  Z2 &operator -= (const Z2 &x) { v ^= x.v; return *this; }
-  Z2 &operator *= (const Z2 &x) { v &= x.v; return *this; }
-  Z2 &operator /= (const Z2 &x) { assert (x.v); return *this; }
+  Z2& operator += (const Z2& x) { v ^= x.v; return *this; }
+  Z2& operator -= (const Z2& x) { v ^= x.v; return *this; }
+  Z2& operator *= (const Z2& x) { v &= x.v; return *this; }
+  Z2& operator /= (const Z2& x) { assert (x.v); return *this; }
   
   // *this += z1*z2
-  Z2 &muladdeq (const Z2 &z1, const Z2 &z2)
+  Z2& muladdeq (const Z2& z1, const Z2& z2)
   {
     return operator += (z1 * z2);
   }
@@ -63,16 +72,21 @@ class Z2
     return Z2 (1);
   }
   
-  tuple<Z2, Z2, Z2> extended_gcd (Z2 x) const
+  std::tuple<Z2, Z2, Z2> extended_gcd (Z2 x) const
   {
     if (v)
-      return make_tuple (Z2 (1), Z2 (1), Z2 (0));
+      return std::make_tuple (Z2 (1), Z2 (1), Z2 (0));
     else
-      return make_tuple (Z2 (1), Z2 (0), Z2 (1));
+      return std::make_tuple (Z2 (1), Z2 (0), Z2 (1));
+  }
+
+  friend std::ostream& operator << (std::ostream& os, const Z2& x) {
+    return os << (int)x.v;
   }
-  
   static void show_ring () { printf ("Z2"); }
   void write_self (writer &w) const { w.write_bool (v); }
-  void show_self () const { printf ("%d", (int)v); }
-  void display_self () const { printf ("%d\n", (int)v); }
+  void show_self () const { std::cout << *this; }
+  void display_self () const { std::cout << *this << "\n"; }
 };
+
+#endif //KNOTKIT_ALGEBRA_Z2_H
diff --git a/algebra/Zp.h b/algebra/Zp.h
index bd3fe0a..80414a0 100644
--- a/algebra/Zp.h
+++ b/algebra/Zp.h
@@ -1,10 +1,16 @@
+#ifndef KNOTKIT_ALGEBRA_ZP_H
+#define KNOTKIT_ALGEBRA_ZP_H
+
+#include <iostream>
+#include <tuple>
+#include <cassert>
 
 template<unsigned p>
 class Zp
 {
  public:
-  typedef ::linear_combination<Zp<p> > linear_combination;
-  typedef ::linear_combination_const_iter<Zp<p> > linear_combination_const_iter;
+  using linear_combination = ::linear_combination<Zp<p>>;
+  using linear_combination_const_iter = ::linear_combination_const_iter<Zp<p>>;
   
  private:
   unsigned v;
@@ -21,84 +27,88 @@ class Zp
     v = i;
   }
   
-  Zp (const Zp &x) : v(x.v) { }
-  Zp (copy, const Zp &x) : v(x.v) { }
+  Zp (const Zp& x) : v(x.v) { }
+  Zp (copy, const Zp& x) : v(x.v) {}
+  Zp (Zp&& x) : v(std::move(x.v)) {
+    x.v = 0;
+  }
   ~Zp () { }
   
-  Zp &operator = (const Zp &x) { v = x.v; return *this; }
-  Zp &operator = (int x)
+  Zp& operator = (const Zp& x) { v = x.v; return *this; }
+  Zp& operator = (Zp&& x) { v = x.v; x.v = 0; return *this; }
+  Zp& operator = (int x)
   {
     return operator = (Zp (x));
   }
   
-  bool operator == (const Zp &x) const { return v == x.v; }
-  bool operator != (const Zp &x) const { return !operator == (x); }
+  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); }
   
-  bool operator < (const Zp &x) const { return v < x.v; }
+  bool operator < (const Zp& x) const { return v < x.v; }
   
   bool is_unit () const
   {
     return v != 0;
   }
   
-  Zp operator + (const Zp &x) const { return Zp (v + x.v); }
-  Zp operator - (const Zp &x) const { return Zp ((int)v - (int)x.v); }
+  Zp operator + (const Zp& x) const { return Zp (v + x.v); }
+  Zp operator - (const Zp& x) const { return Zp ((int)v - (int)x.v); }
   Zp operator - () const { return Zp (- (int)v);  }
   
-  Zp operator * (const Zp &x) const { return Zp (v * x.v); }
+  Zp operator * (const Zp& x) const { return Zp (v * x.v); }
   
-  Zp operator / (const Zp &x) const
+  Zp operator / (const Zp& x) const
   {
     return operator * (x.recip ());
   }
   
   Zp recip () const
   {
-    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);
+    std::tuple<unsigned, int, int> t = unsigned_extended_gcd (v, p);
+    assert (std::get<0> (t) == 1);
+    assert ((int)std::get<0> (t) == std::get<1> (t)*(int)v  + std::get<2> (t)*(int)p);
     
-    return Zp (get<1> (t));
+    return Zp (std::get<1> (t));
   }
   
-  Zp &operator += (const Zp &x)
+  Zp& operator += (const Zp& x)
   {
     v = (v + x.v) % p;
     return *this;
   }
   
-  Zp &operator -= (const Zp &x) { return operator = (Zp ((int)v - (int)x.v)); }
+  Zp& operator -= (const Zp& x) { return operator = (Zp ((int)v - (int)x.v)); }
   
-  Zp &operator *= (const Zp &x)
+  Zp& operator *= (const Zp& x)
   {
     v = (v * x.v) % p;
     return *this;
   }
   
-  Zp &operator /= (const Zp &x)
+  Zp& operator /= (const Zp& x)
   {
     return operator *= (x.recip ());
   }
   
   // d | n, d.divides (n)
-  bool divides (const Zp &n) const
+  bool divides (const Zp& n) const
   {
     return v || !n.v;
   }
   
-  bool operator | (const Zp &n) const { return divides (n); }
+  bool operator | (const Zp& n) const { return divides (n); }
   
-  Zp div (const Zp &d) const { return operator / (d); }
+  Zp div (const Zp& d) const { return operator / (d); }
   
-  tuple<Zp, Zp, Zp> extended_gcd (const Zp &x) const
+  std::tuple<Zp, Zp, Zp> extended_gcd (const Zp &x) const
   {
     if (v)
-      return make_tuple (v, Zp (1), Zp (0));
+      return std::make_tuple (v, Zp (1), Zp (0));
     else
-      return make_tuple (x, Zp (0), Zp (1));
+      return std::make_tuple (x, Zp (0), Zp (1));
   }
   
   Zp gcd (const Zp &x) const
@@ -106,8 +116,12 @@ class Zp
     assert (v || x.v);
     return 1;
   }
-  
-  static void show_ring () { printf ("Z%d", p); }
-  void show_self () const { printf ("%d (%d)", v, p); }
-  void display_self () const { show_self (); newline (); }
+  friend std::ostream& operator << (std::ostream& os, const Zp& x) {
+    return os << x.v;
+  }
+  static void show_ring () { std::cout << "Z" << p; }
+  void show_self () const { std::cout << v << "(" << p << ")"; }
+  void display_self () const { std::cout << *this << "\n"; }
 };
+
+#endif //KNOTKIT_ALGEBRA_ZP_H
diff --git a/algebra/algebra.h b/algebra/algebra.h
index 4acd0ca..b23bd79 100644
--- a/algebra/algebra.h
+++ b/algebra/algebra.h
@@ -1,6 +1,10 @@
-
 #include <lib/lib.h>
 
+#ifdef DEBUG_ALGEBRA
+#define DEBUG_Z
+#define DEBUG_Q
+#endif
+
 inline int recip (int x)
 {
   assert (x == 1 || x == -1);
diff --git a/kk.cpp b/kk.cpp
index b3803ef..2ddb2a7 100644
--- a/kk.cpp
+++ b/kk.cpp
@@ -1,4 +1,3 @@
-
 #include <knotkit.h>
 
 const char *program_name;
diff --git a/lib/lib.h b/lib/lib.h
index 1bb300c..4537ade 100644
--- a/lib/lib.h
+++ b/lib/lib.h
@@ -8,7 +8,8 @@
 #include <sys/stat.h>
 #include <sys/errno.h>
 
-#include <gmp.h>
+#include <gmpxx.h>
+//#include <gmp.h>
 #include <zlib.h>
 
 #include <initializer_list>