diff --git a/Makefile b/Makefile
index a67fcca..1210ce3 100644
--- a/Makefile
+++ b/Makefile
@@ -8,8 +8,8 @@ CXX = clang++ -fno-color-diagnostics --stdlib=libc++ --std=c++11 -I/u/cseed/llvm
 
 INCLUDES = -I/opt/local/include -I.
 
-# OPTFLAGS = -g
-OPTFLAGS = -O2 -g
+OPTFLAGS = -g
+# OPTFLAGS = -O2 -g
 # OPTFLAGS = -O2 -DNDEBUG
 
 LDFLAGS = -L/opt/local/lib -L/u/cseed/llvm-3.1/lib
diff --git a/algebra/Z2.h b/algebra/Z2.h
index 4245c7f..d38237b 100644
--- a/algebra/Z2.h
+++ b/algebra/Z2.h
@@ -41,6 +41,12 @@ class Z2
   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)
+  {
+    return operator += (z1 * z2);
+  }
+  
   bool divides (Z2 x) const { return v || !x.v; }
   bool operator | (const Z2 x) const { return divides (x); }
   
diff --git a/algebra/module.h b/algebra/module.h
index 5395c76..69faef8 100644
--- a/algebra/module.h
+++ b/algebra/module.h
@@ -63,7 +63,12 @@ class module : public refcounted
     if (i <= free_rank ())
       return c == 0;
     else
-      return generator_ann (i) | c;
+      {
+	// ???
+	abort ();
+	
+	// return generator_ann (i) | c;
+      }
   }
   
   R annihilator (R c, unsigned i) const
@@ -1853,7 +1858,15 @@ mod_map<R>::kernel () const
 	  linear_combination<R> &to_x = to_xs[j],
 	    &from_x = from_xs[j];
 	  R to_xc = to_x(i);
-	  if (! (to_vc | to_xc))
+	  if (to_xc == 0)
+	    continue;
+	  
+	  if (to_vc == 0 && to_xc != 0)
+	    {
+	      to_v += to_x;
+	      from_v += from_x;
+	    }
+	  else if (! (to_vc | 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);
@@ -1959,8 +1972,14 @@ mod_span<R>::mod_span (ptr<const module<R> > mod,
 	  
 	  linear_combination<R> &x = xs[j];
 	  R xc = x(i);
+	  if (xc == 0)
+	    continue;
 	  
-	  if (! (vc | xc))
+	  if (vc == 0)
+	    {
+	      v += x;
+	    }
+	  else if (! (vc | xc))
 	    {
 	      tuple<R, R, R> t = vc.extended_gcd (xc);
 	      assert (get<0> (t) == vc*get<1> (t) + get<2> (t)*xc);
diff --git a/algebra/multivariate_laurentpoly.h b/algebra/multivariate_laurentpoly.h
index 8a08828..8e9f3a1 100644
--- a/algebra/multivariate_laurentpoly.h
+++ b/algebra/multivariate_laurentpoly.h
@@ -219,7 +219,8 @@ class multivariate_laurentpoly
   
   multivariate_laurentpoly (const multivariate_laurentpoly &p) : coeffs(p.coeffs) { }
   multivariate_laurentpoly (copy, const multivariate_laurentpoly &p)
-    : coeffs(COPY2, p.coeffs)
+    // ??? COPY2?
+    : coeffs(COPY, p.coeffs)
   {
   }
 
diff --git a/algebra/polynomial.h b/algebra/polynomial.h
index 5c3abcd..4d1771f 100644
--- a/algebra/polynomial.h
+++ b/algebra/polynomial.h
@@ -67,7 +67,7 @@ public:
   
   polynomial recip () const
   {
-    assert (coeffs.card () != 1);
+    assert (coeffs.card () == 1);
     
     pair<unsigned, T> p = coeffs.head ();
     assert (p.first == 0);
@@ -134,13 +134,16 @@ public:
   polynomial operator * (const polynomial &p) const;
   
   bool divides (const polynomial &num) const;
+  
+  // *this | num
   bool operator | (const polynomial &num) const { return divides (num); }
   
   tuple<polynomial, polynomial> divide_with_remainder (const polynomial &denom) const;
   
   polynomial mod (const polynomial &denom) const;
-  
+
   polynomial divide_exact (const polynomial &denom) const;
+  polynomial div (const polynomial &denom) const { return divide_exact (denom); }
   
   polynomial gcd (const polynomial &b) const;
   polynomial lcm (const polynomial &b) const;
@@ -392,14 +395,14 @@ class polynomial<Z2>
   
  public:
   polynomial () { }
-  explicit polynomial (int x)
+  polynomial (int x)
   {
     Z2 c (x);
     if (c != 0)
       coeffs.push (0);
   }
   
-  explicit polynomial (Z2 c)
+  polynomial (Z2 c)
   {
     if (c != 0)
       coeffs.push (0);
@@ -450,6 +453,17 @@ class polynomial<Z2>
   bool operator != (int x) const { return !operator == (x); }
   
   unsigned degree () const { return coeffs.tail (); }
+  bool is_unit () const
+  {
+    return coeffs.card () == 1
+      && coeffs.head () == 0;
+  }
+  
+  polynomial recip () const
+  {
+    assert (is_unit ());
+    return 1;
+  }
   
   polynomial &operator += (polynomial p) { coeffs ^= p.coeffs; return *this; }
   polynomial &operator -= (polynomial p) { coeffs ^= p.coeffs; return *this; }
diff --git a/cube_impl.h b/cube_impl.h
index 24fd900..56d7662 100644
--- a/cube_impl.h
+++ b/cube_impl.h
@@ -166,7 +166,7 @@ cube<R>::compute_map (unsigned dh, unsigned max_n,
 			continue;
 		      assert (!unsigned_bittest (v_to, p));
 		    }
-
+		  
 		  b[generator (fromstate, v_from)].muladd
 		    (sign, generator (tostate, v_to));
 		}
@@ -302,7 +302,7 @@ cube<R>::compute_X (unsigned p) const
 template<class R> mod_map<R>
 cube<R>::H_i (unsigned c)
 {
-  mod_map<R> b (khC, 0);
+  map_builder<R> b (khC, 0);
   for (unsigned i = 0; i < n_resolutions; i ++)
     {
       if (unsigned_bittest (i, c))
@@ -319,6 +319,14 @@ cube<R>::H_i (unsigned c)
       
       unsigned from = from_s.crossing_from_circle (kd, c),
 	to = from_s.crossing_to_circle (kd, c);
+      
+      int sign = 1;
+      for (unsigned j = 1; j < c; j ++)
+        {
+          if (unsigned_bittest (i, j))
+            sign *= -1;
+        }
+      
       for (unsigned j = 0; j < from_s.num_monomials (); j ++)
 	{
 	  unsigned j2 = 0;
@@ -336,7 +344,7 @@ cube<R>::H_i (unsigned c)
 				    to_s.crossing_from_circle (kd, c));
 	      j2 = unsigned_bitset (j2,
 				    to_s.crossing_to_circle (kd, c));
-	      v.muladd (1, generator (i2, j2));
+	      v.muladd (-sign, generator (i2, j2));
 	    }
 	  else if (from != to
 		   && !unsigned_bittest (j, from)
@@ -344,7 +352,7 @@ cube<R>::H_i (unsigned c)
 	    {
 	      assert (! unsigned_bittest (j2,
 					  to_s.crossing_from_circle (kd, c)));
-	      v.muladd (1, generator (i2, j2));
+	      v.muladd (sign, generator (i2, j2));
 	    }
 	}
     }
diff --git a/main.cpp b/main.cpp
index f43ee8b..d098570 100644
--- a/main.cpp
+++ b/main.cpp
@@ -850,18 +850,104 @@ test_forgetful_ss ()
 }
 
 void
-test_forgetful_signs ()
+compute_forgetful_torsion ()
 {
-  typedef Zp<97> R;
+  typedef polynomial<Zp<2> > R;
   
-#if 0
-  for (unsigned i = 1; i <= 14; i ++)
+  for (unsigned i = 1; i <= 12; i ++)
     for (unsigned j = 1; j <= mt_links (i, 0); j ++)
       {
 	knot_diagram kd (mt_link (i, 0, j));
-      }
+	unsigned n = kd.num_components ();
+	if (n != 2)
+	  continue;
+	
+	show (kd); newline ();
+	
+	unionfind<1> u (kd.num_edges ());
+	for (unsigned i = 1; i <= kd.n_crossings; i ++)
+	  {
+	    u.join (kd.ept_edge (kd.crossings[i][1]),
+		    kd.ept_edge (kd.crossings[i][3]));
+	    u.join (kd.ept_edge (kd.crossings[i][2]),
+		    kd.ept_edge (kd.crossings[i][4]));
+	  }
+	
+	map<unsigned, unsigned> root_comp;
+	unsigned t = 0;
+	for (unsigned i = 1; i <= kd.num_edges (); i ++)
+	  {
+	    if (u.find (i) == i)
+	      {
+		++ t;
+		root_comp.push (i, t);
+	      }
+	  }
+	assert (t == n);
+	
+	cube<R> c (kd);
+#if 0
+	for (unsigned i = 0; i < c.n_resolutions; i ++)
+	  {
+	    smallbitset state (kd.n_crossings, i);
+	    smoothing s (kd, state);
+	    s.show_self (kd, state);
+	    newline ();
+	  }
 #endif
+	
+	mod_map<R> untwisted_d = c.compute_d (1, 0, 0, 0, 0);
+	assert (untwisted_d.compose (untwisted_d) == 0);
+	
+	mod_map<R> d = untwisted_d;
+	for (unsigned x = 1; x <= kd.n_crossings; x ++)
+	  {
+	    unsigned r1 = u.find (kd.ept_edge (kd.crossings[x][1])),
+	      r2 = u.find (kd.ept_edge (kd.crossings[x][2]));
+	    
+	    if (r1 != r2)
+	      d = d + c.compute_dinv (x)*(R (1, 1));
+	  }
+	
+#if 0
+	mod_map<R> h = d.compose (d);
+	display ("h:\n", h);
+#endif
+	
+	assert (d.compose (d) == 0);
+	
+	// display ("d:\n", d);
+	
+	chain_complex_simplifier<R> s (c.khC, d, 0);
+	printf ("|s.new_C| = %d\n", s.new_C->dim ());
+	display ("s.new_d:\n", s.new_d);
+	
+#if 0
+	ptr<const free_submodule<R> > ker = s.new_d.kernel ();
+	printf ("|ker| = %d\n", ker->dim ());
+	
+	ptr<const free_submodule<R> > im = s.new_d.image ();
+	printf ("|im| = %d\n", im->dim ());
+#endif
+	
+	ptr<const module<R> > H = s.new_d.homology ();
+	display ("H:\n", *H);
+	// printf ("dim H = %d, free rank H = %d\n", H->dim (), H->free_rank ());
+      }
+}
+
+void
+test_forgetful_signs ()
+{
+  typedef Q R;
+  
 #if 1
+  for (unsigned i = 1; i <= 14; i ++)
+    for (unsigned j = 1; j <= mt_links (i, 1); j ++)
+      {
+	knot_diagram kd (mt_link (i, 1, j));
+#endif
+#if 0
   for (unsigned i = 2; i <= 16; i ++)
     for (unsigned j = i; j <= 16; j ++)
       {
@@ -870,6 +956,7 @@ test_forgetful_signs ()
 		&& kd.n_crossings <= i * (j - 1));
 	if (kd.n_crossings > 13)
 	  continue;
+      }
 #endif
 	
 	unsigned n = kd.num_components ();
@@ -1026,6 +1113,110 @@ test_forgetful_signs ()
       }
 }
 
+void
+compute_lee_bound ()
+{
+  typedef Z2 R;
+  
+#if 0
+  for (unsigned i = 2; i <= 16; i ++)
+    for (unsigned j = i; j <= 16; j ++)
+      {
+	knot_diagram kd (torus_knot (i, j));
+	assert (kd.n_crossings == (i - 1) * j
+		&& kd.n_crossings <= i * (j - 1));
+	if (kd.n_crossings > 13)
+	  continue;
+      }
+#endif
+#if 1
+  for (unsigned i = 1; i <= 10; i ++)
+    for (unsigned j = 1; j <= rolfsen_crossing_knots (i); j ++)
+      {
+	knot_diagram kd (rolfsen_knot (i, j));
+#endif
+	
+	show (kd); newline ();
+	
+	cube<R> c (kd);
+#if 0
+	for (unsigned i = 0; i < c.n_resolutions; i ++)
+	  {
+	    smallbitset state (kd.n_crossings, i);
+	    smoothing s (kd, state);
+	    s.show_self (kd, state);
+	    newline ();
+	  }
+#endif
+	
+	mod_map<R> d = c.compute_d (1, 0, 0, 0, 0);
+	assert (d.compose (d) == 0);
+	
+	for (unsigned x = 1; x <= kd.n_crossings; x ++)
+	  d = d + c.H_i (x);
+	
+	assert (d.compose (d) == 0);
+	
+	ptr<const module<R> > Ek = c.khC;
+	mod_map<R> dk = d;
+	
+	unsigned kinf;
+	for (int dq = 0;; dq += 2)
+	  {
+	    chain_complex_simplifier<R> s (Ek, dk, dq);
+	    Ek = s.new_C;
+	    dk = s.new_d;
+	    
+	    printf ("|E%d| = %d\n", (dq / 2) + 1, Ek->dim ());
+	    if (dk == 0)
+	      {
+		kinf = (dq / 2);
+		break;
+	      }
+	  }
+	
+	printf ("kinf = %d\n", kinf);
+      }
+}
+
+void
+test_crossing_ss ()
+{
+  typedef multivariate_laurentpoly<Z2> R;
+  
+  for (unsigned i = 1; i <= 14; i ++)
+    for (unsigned j = 1; j <= htw_knots (i); j ++)
+      {
+	knot_diagram kd (htw_knot (i, j));
+	show (kd); newline ();
+	
+	cube<R> c (kd);
+#if 0
+	for (unsigned i = 0; i < c.n_resolutions; i ++)
+	  {
+	    smallbitset state (kd.n_crossings, i);
+	    smoothing s (kd, state);
+	    s.show_self (kd, state);
+	    newline ();
+	  }
+#endif
+	
+	mod_map<R> untwisted_d = c.compute_d (1, 0, 0, 0, 0);
+	assert (untwisted_d.compose (untwisted_d) == 0);
+	
+	mod_map<R> d = untwisted_d;
+	for (unsigned x = 1; x <= kd.n_crossings; x ++)
+	  {
+	    d = d + c.compute_dinv (x)*(R (1, VARIABLE, x));
+	  }
+	
+	display ("d:\n", d);
+	
+	mod_map<R> h = d.compose (d);
+	display ("h:\n", h);
+      }
+}
+
 /* Kh homotopy type:
    - check for CP^2,
    - mutation invariance
@@ -1752,12 +1943,138 @@ convert_15 ()
   write (w, knot15_sq);
 }
 
+void
+compare_gss_splitting ()
+{
+  typedef Z2 R;
+  
+#if 1
+  for (unsigned i = 1; i <= 14; i ++)
+    for (unsigned j = 1; j <= mt_links (i, 1); j ++)
+      {
+	knot_diagram kd (mt_link (i, 1, j));
+#endif
+	
+	unsigned n = kd.num_components ();
+	if (n < 2)
+	  continue;
+	
+	show (kd); newline ();
+	
+	unionfind<1> u (kd.num_edges ());
+	for (unsigned i = 1; i <= kd.n_crossings; i ++)
+	  {
+	    u.join (kd.ept_edge (kd.crossings[i][1]),
+		    kd.ept_edge (kd.crossings[i][3]));
+	    u.join (kd.ept_edge (kd.crossings[i][2]),
+		    kd.ept_edge (kd.crossings[i][4]));
+	  }
+	
+	map<unsigned, unsigned> root_comp;
+	unsigned t = 0;
+	for (unsigned i = 1; i <= kd.num_edges (); i ++)
+	  {
+	    if (u.find (i) == i)
+	      {
+		++ t;
+		root_comp.push (i, t);
+	      }
+	  }
+	assert (t == n);
+	
+	cube<R> c (kd);
+#if 1
+	for (unsigned i = 0; i < c.n_resolutions; i ++)
+	  {
+	    smallbitset state (kd.n_crossings, i);
+	    smoothing s (kd, state);
+	    s.show_self (kd, state);
+	    newline ();
+	  }
+#endif
+
+	printf ("inter-component x:\n");
+	mod_map<R> ds (c.khC);
+	for (unsigned x = 1; x <= kd.n_crossings; x ++)
+	  {
+	    unsigned p1 = kd.crossings[x][1],
+	      p2 = kd.crossings[x][2];
+	    unsigned r1 = u.find (kd.ept_edge (p1)),
+	      r2 = u.find (kd.ept_edge (p2));
+	    unsigned c1 = root_comp(r1),
+	      c2 = root_comp(r2);
+
+	    if (c1 != c2)
+	      {
+		printf ("  %d\n", x);
+		
+		ds = ds + c.compute_dinv (x);
+	      }
+	  }
+	
+	mod_map<R> d1 = c.compute_d (1, 0, 0, 0, 0);
+	assert (d1.compose (ds) == ds.compose (d1));
+
+	mod_map<R> d2 = c.compute_d (2, 0, 0, 0, 0);
+	mod_map<R> h2 = d2.compose (ds) + ds.compose (d2);
+	if (h2 != 0)
+	  display ("h2:\n", h2);
+	
+	mod_map<R> d3 = c.compute_d (3, 0, 0, 0, 0);
+	mod_map<R> h3 = d3.compose (ds) + ds.compose (d3);
+	if (h3 != 0)
+	  display ("h3:\n", h3);
+	
+	mod_map<R> d = c.compute_d (0, 0, 0, 0, 0);
+	mod_map<R> h = d.compose (ds) + ds.compose (d);
+	if (h != 0)
+	  display ("h:\n", h);
+      }
+}
+
 int
 main ()
 {
-  test_forgetful_signs ();
+  compute_lee_bound ();
   return 0;
   
+  compare_gss_splitting ();
+  
+  {
+    knot_diagram kd (rolfsen_knot (3, 1));
+    cube<Q> c (kd);
+    mod_map<Q> d = c.compute_d (1, 0, 0, 0, 0);
+    chain_complex_simplifier<Q> s (c.khC, d, 0);
+    assert (s.new_d == 0);
+    display ("s.new_C:\n", *s.new_C);
+  }
+  return 0;
+  
+  test_forgetful_signs ();
+  
+  {
+    knot_diagram kd (torus_knot (2, 2));
+    printf ("n+ = %d, n- = %d\n", kd.nplus, kd.nminus);
+    cube<Z2> c (kd);
+    sseq ss = compute_szabo_sseq (c);
+    ss.texshow (stdout, "Hopf");
+  }
+  {
+    knot_diagram kd (MIRROR, knot_diagram (torus_knot (2, 2)));
+    printf ("n+ = %d, n- = %d\n", kd.nplus, kd.nminus);
+    cube<Z2> c (kd);
+    sseq ss = compute_szabo_sseq (c);
+    ss.texshow (stdout, "Hopf");
+  }
+  
+  return 0;
+  
+  compute_forgetful_torsion ();
+  return 0;
+  
+  test_forgetful_signs ();
+  test_crossing_ss ();
+  
   find_width4_in_h0 ();
   // convert_15 ();
   // test_knot_sq ();