E := matrix(16,16,[[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0],[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0],[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1],[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0],[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0],[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0],[0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1],[0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0]]);
a := matrix(1,16,[1,1,1,0,1,0,1,0,1,0,1,0,0,0,0,0]);
b := matrix(16,1,[[1],[1],[0],[1],[0],[1],[0],[1],[0],[1],[0],[1],[0],[0],[0],[0]]);
evalm(a &* E^3 &* b);
f := proc(a,M,c,x)
local m,q,V,i,g:
m := linalg[charpoly](M,x):
q := 0:
V := array(0..degree(m,x)):
for i from 0 to degree(m,x) do
V[i] := evalm(a &* M^(i+1) &* c)[1,1]*x^i:
end do:
for i from 0 to degree(m,x) do
q := q + coeff(m,x,i)*(g-(sum(V[j],j=0..i)))/x^i:
end do:
simplify(solve(q = 0,g));
end proc;
f(a,E,b,x);
g := %;
taylor(g,x,10);
k := (2+3*x-x^2)/(x^3-3*x^2-3*x+1);
taylor(k,x,10);
g := g*x+2;
g := simplify(g);
k := simplify(k);