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;
g := proc(d,q,x)
local M,a,b,i,j,k:
k := simplify(q);
M := matrix(d,d):
a := matrix(1,d,1):
b := matrix(d,1,1):
for i from 0 to 2^(d^2)-1 do
for j from 0 to d^2-1 do M[iquo(j,d)+1,(j mod d)+1] := iquo(i,2^j) mod 2: end do:
if f(a,M,b,x) = k then return M: end if:
end do:
return matrix(d,d,0);
end proc;
g(3,(2+3*x-x^2)/(x^3-3*x^2-3*x+1),x);