with(linalg):unassign('C','Chip','Chia','eta','q','n','k','pmax','R','Z','qv','pv','polyp','polya'):Digits:=15:q:=17: n:=5: k:=3: pmax:=5:ch:=theta->z^2-2*cos(2*Pi*theta)*z+1:qv:=[1,3,4];pv:=convert(({1,2,3,4,5,6,7,8} minus convert(qv,set)),list);C:=array(-pmax..(2*k+pmax),0..pmax):Chip:=array(0..pmax,1..q):Chia:=array(0..k,1..q):eta:=array(-pmax..(k+1),0..pmax):for t from 1 to q do polyp:=simplify(simplify(ch(t*pv[1]/q)*ch(t*pv[2]/q))*simplify(ch(t*pv[3]/q)*ch(t*pv[4]/q)*ch(t*pv[5]/q))); for p from 0 to pmax do Chip[p,t]:=evalf((-1)^p*coeff(polyp,z,p)) od od:for t from 1 to q do polya:=1: for m1 from 1 to k do polya:=simplify(polya*ch(t*qv[m1]/q)) od: for a from 0 to k do Chia[a,t]:=(-1)^a*evalf(eval(coeff(polya,z,a))) od od:for a from 0 to k do for p from 0 to pmax do C[a,p]:=sum(eval(Chia[a,m3]*Chip[p,m3]), m3=1..q) od od:for a from -pmax to -1 do for p from 0 to pmax do C[a,p]:=0 od od:for a from (k+1) to 2*k do for p from 0 to pmax do C[a,p]:=C[2*k-a,p] od od:for a from (2*k+1) to (2*k+pmax) do for p from 0 to pmax do C[a,p]:=0 od od:for p from 0 to pmax do for b from -p to 1 do eta[b,p]:=(-1)^b*sum(C[m4+b,p-m4],m4=0..p) od od:for p from 0 to pmax do for b from 2 to (k+1) do eta[b,p]:=(-1)^b*sum(simplify(C[m5+b,p-m5]-C[b-2-m5,p-m5]),m5=0..p) od od:for p from 0 to pmax do for b from -p to (k+1) do eta[b,p]:=round(evalf(eta[b,p])) od od:for a from -pmax to 2*k+pmax do for t from 0 to pmax do C[a,t]:=round(evalf(C[a,t])) od od:convert(C,matrix); e:=convert(eta,matrix): transpose(e);print(qv,pv);