improve matrix_to_formula()

graph_tools/utils.py

@@ -162,6 +162,30 @@ def _init_():
 
 
   def matrix_to_formula(l, A, r, k=None, base_ring=None):
+    """Tranform matrix formula l A^k r to simple formula.
+
+    First rewrite A as P J P^-1 using Jordan normal form.
+    Then convert to formula using sage symolic expressions.
+
+    Sage (backed by Maxima) simplifies 0^x to 0 which is
+    not correct because we need 0^0 == 1. So we check correctness
+    of resulting formula for values of k less than size of
+    the largest block of J.
+
+    EXAMPLE:
+    We derive formula for Fibonacci numbers:
+    >>> from sage.all import *
+    >>> m = matrix(QQ, 2, 2, [1, 1, 1, 0])
+    >>> r = matrix(QQ, 2, 1, [1, 0])
+    >>> l = matrix(QQ, 1, 2, [0, 1])
+    >>> [ l * m**i * r for i in range(10) ]
+    [[0], [1], [1], [2], [3], [5], [8], [13], [21], [34]]
+    >>> f = matrix_to_formula(l, m, r, base_ring=QQbar)['formula']; f
+    0.4472135954999580?*1.618033988749895?^k - 0.4472135954999580?*(-0.618033988749895?)^k
+    >>> [ f(k=i) for i in range(10) ]
+    [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]
+    """
+
     from sage.all import binomial, SR
 
     J, P = A.jordan_form(transformation=True, sparse=True, base_ring=base_ring)
@@ -169,6 +193,7 @@ def _init_():
     l = l * P
     r = P**(-1) * r
 
+    if k is None: k = SR.symbol('k')
     out_k = k
     k = SR.symbol('k')
     f = 0
@@ -188,7 +213,9 @@ def _init_():
           binomial(k, ll) * J[ind, ind]**(k - ll) * r[ind + ll, 0] for ll in range(s-j)
         )
 
+    try:
       f = f.simplify_full()
+    except: pass
 
     while max_s > 0:
       val_r = (l * J**(max_s - 1) * r)[0,0]
@@ -196,7 +223,9 @@ def _init_():
       if val_r != val_f: break
       max_s -= 1
 
+    try:
       f = f.subs(k=out_k).simplify_full()
+    except: pass
 
     return {
       'simplified': (l, J, r),