import numpy as np
import numpy.linalg as alg

A = np.array([[64,0,-28,17],[12,-15,0,38],[-4,1,64,17],[-3,5,0,32]])
A1 = A
L = [1]
for k in range(1,5) :
    a = -np.trace(A1)/k
    L.append(a)
    A1 = np.dot(A,A1+a*np.eye(4))
    
B = np.zeros((4,4))
for k in range(4,-1,-1) :
    B += L[k]*alg.matrix_power(A,4-k)
print(B)
    
print(L,np.poly(A))