def getMat(k):
    vs = [(m, M) for m in range(k) for M in range(m,k)]
    vToInd = {v:i for i,v in enumerate(vs)}
    mat = matrix(QQ, len(vs))
    for i1, (m,M) in enumerate(vs):
        for x1 in range(m, k):
            for x2 in range(m, k):
                if max(x1,x2)<M: continue
                v2 = (min(x1,x2), max(x1,x2))
                mat.add_to_entry(i1, vToInd[v2], 1)
    return mat

def getGraph(k):
    g = DiGraph()
    g.allow_loops(True)
    mat = getMat(k)
    vs = [(m, M) for m in range(k) for M in range(m,k)]
    for i1, v1 in enumerate(vs):
        for i2, v2 in enumerate(vs):
            if mat[i1,i2]>0:
                g.add_edge(v1,v2,mat[i1,i2])
    return g


def f(n, k):
    mat = getMat(k)
    return sum((mat^n)[0])


#g = getGraph(3)
#show(g.plot(figsize=6, edge_labels=True))
#show(getMat(3))






def getFormula(k, verbose=0):
    if verbose: print ("Finding formula for k = ", k)

    #b = 2^n
    R.<n, b> = PolynomialRing(QQ, 2)

    #ell is the eigenvalue, either 1 or 2
    def jord(m, ell):
        return matrix([[binomial(n,j-i)*(1 if ell==1 else b/ell^(j-i)) if j>=i else 0 for j in range(m)] for i in range(m)])


    A = getMat(k)
    matN = A.dimensions()[0]
    D, J = A.jordan_form(transformation=True)
    divs = D.subdivisions()[0]
    #assuming len(divs)>0
    ms = [divs[0]]
    eigs = [D[0,0]]
    for j in range(len(divs)-1):
        ms.append(divs[j+1]-divs[j])
        eigs.append(D[divs[j], divs[j]])
    ms.append(matN-divs[-1])
    eigs.append(D[matN-1,matN-1])
    if verbose: print ("Jordan blocks: " + ", ".join("J_{%d}(%d)" %(m,l) for m,l in zip(ms,eigs)))
    if verbose==2: show(D)
    #show(J*D*J^(-1), A)
    vMat = block_diagonal_matrix(*(jord(m,l) for m,l in zip(ms, eigs)))
    return sum((J*vMat*J^(-1))[0])


#testK = 7
#g = getFormula(testK)
#show(factor(g))
#test
#print ([f(n,testK) for n in range(1, 10)])
#print ([g(n=n, b=2^n) for n in range(1, 10)])

for k in range(2, 10):
    show (k, ": ", factor(getFormula(k)))


#print ([2^(n+2)-n-3 for n in range(1,10)])
#print ([1/2*(n+2)*(2^(n+2)*(n-1)+n+5) for n in range(1,10)])
#print ([sum(j*2^(n+1-j) for j in range(n+2)) for n in range(1, 10)])

#print ([list(r) for r in getMat(3)])

#a = matrix([[f(n,k) for n in range(1, 10)] for k in range(1, 10)])
#show(a)