GraphSequenceBase.stabilize() return also a formula

graph_tools/base.py

@@ -224,7 +224,7 @@ def _init_():
       return self._graph_cache.get(n, lambda: self._make_graph(self.gadget(n)))
 
     def stabilize(self, parameter, start_at=None):
-      from sage.all import matrix, QQ
+      from sage.all import matrix, QQ, binomial, SR
 
       if start_at is None: start_at = self.sequence_start
       start_gadget = self.gadget(start_at)
@@ -267,11 +267,44 @@ def _init_():
         for col in sorted_boundaries
       })
 
+      J, P = m.jordan_form(transformation=True)
+      assert P * J * P**(-1) == m, "Error %s != %s" % (P * J * P**(-1), m)
+      l = fin * P
+      r = P**(-1) * iv
+
+      k = SR.symbol('k') - start_at
+      f = 0
+
+      subdivisions = J.subdivisions()
+      assert subdivisions[0] == subdivisions[1]
+      subdivisions = [ 0 ] + subdivisions[0] + [ J.nrows() ]
+      to_delete = []
+
+      for i in range(len(subdivisions) - 1):
+        if J.subdivision(i, i).is_zero():
+          to_delete.extend(range(subdivisions[i], subdivisions[i+1]))
+        else:
+          s = subdivisions[i+1] - subdivisions[i]
+          for j in range(s):
+            ind = subdivisions[i] + j
+            f += l[0, ind] * sum(
+              binomial(k, j) * J[ind, ind]**(k - j) * r[ind + j, 0]  for j in range(s)
+            )
+
+      l = l.delete_columns(to_delete)
+      J = J.delete_rows(to_delete)
+      J = J.delete_columns(to_delete)
+      r = r.delete_rows(to_delete)
+
+      f = f.simplify_full()
+
       return {
         'variables': sorted_boundaries,
         'step_matrix': m,
         'initial_vector': iv,
-        'finalize': fin
+        'finalize': fin,
+        'simplified': (l, J, r),
+        'formula': f
       }
 
 

graph_tools/tests.py

@@ -20,7 +20,11 @@ True
 True
 >>> Necklace.stabilize(CircuitDoubleCover)
 Loop 1 done - 1 boundaries (0 new)
-{'variables': [((1, 2), (1, 2), None)], 'step_matrix': [4], 'initial_vector': [2], 'finalize': [1]}
+{'variables': [((1, 2), (1, 2), None)], 'step_matrix': [4], 'initial_vector': [2], \
+'finalize': [1], 'simplified': ([1], [4], [2]), 'formula': 1/2*4^k}
+>>> Flower.stabilize(CircuitDoubleCover, 2)['formula']
+Loop 1 done - 3 boundaries (0 new)
+1/3*2^(k - 1) + 1/3*(-1)^k + 1
 >>> sun_cdc = lambda n: sum( b.value for b in sun(n).eval_gadget(CircuitDoubleCover) )
 >>> all( sun_cdc(n) == 2**n - 2*n for n in range(3, 10) )
 True