# Maple file for the computations in "A degree theoretic approach to solvability of
# symmetric word equations in two positive definite letters"with(linalg):kronprod := proc (A, B)
options `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 1998 by Robert B. Israel. All rights reserved`;
local Ap, Bp, i,j;
if nargs > 2 then RETURN(kronprod(kronprod(A,B),args[3..nargs])) fi;
if type(A,{vector,list(algebraic)}) and
type(B,{vector,list(algebraic)}) then
# vector x vector = vector
vector([seq(seq(A[i]*B[j], j=1..linalg[vectdim](B)),
i=1..linalg[vectdim](A))])
else # otherwise result is matrix
if type(A,matrix) then Ap:= A
elif type(A,listlist) then Ap:= convert(A,matrix)
elif type(A,list) then Ap:= matrix(map(t->[t],A))
elif type(A,specfunc(list,transpose)) then Ap:= matrix([op(A)])
else Ap:= convert(A,matrix)
fi;
if type(B,matrix) then Bp:= B
elif type(B,listlist) then Bp:= convert(B,matrix)
elif type(B,list) then Bp:= matrix(map(t->[t],B))
elif type(B,specfunc(list,transpose)) then Bp:= matrix([op(B)])
else Bp:= convert(B,matrix)
fi;
linalg[stackmatrix](seq(linalg[augment](
seq(linalg[scalarmul](Bp,Ap[i,j]), j = 1 .. linalg[coldim](Ap))),
i = 1 .. linalg[rowdim](Ap)));
fi
end;isequal := proc (a,b) # are a and b equal?
if (a = b) then return(1); end;
return(0);
end;matrixpower := proc( M, n, p ) # finds the n-by-n matrix power M^p
local temp, temp2;
if (p = 1) then return(M); end;
if (p = 0) then return(matrix(n,n,(i,j)->isequal(i,j))); end;
if (modp(p,2) = 0) then
temp := matrixpower(M,n,p/2);
return(evalm(temp&*temp));
else
temp := M;
temp2 := matrixpower(M,n,(p-1)/2);
return(evalm(temp2&*temp2&*temp));
end;
end;evalword := proc (X,B,W,n,k,start) # evalautes a word with k products of powers at two n-by-n X and B
# word W = X^p_1 B^q1 X^p2... is given as an integer vector [p_1,q_1,...]
# e.g. word XB^2XB^4 is [1,2,1,4]
# start designates where in the vector to start multiplying left to right,
local i,j, totalProd;
totalProd := evalm(matrix(n,n,(i,j)->isequal(i,j))); # identity matrix
for i from start to k do # construct product
if (modp(i,2) = 1) then
totalProd := evalm(totalProd&*matrixpower(X,n,W[i]));
else
totalProd := evalm(totalProd&*matrixpower(B,n,W[i]));
end;
od;
return evalm(totalProd);
end; wordjacobian := proc (X,B,W,n,k) # finds jacobian of word W with k products of powers evaluated at n-by-n X and B
# word W = X^p_1 B^q1 X^p2... is given as an integer vector [p_1,q_1,...]
# e.g. word XB^2XB^4X is [1,2,1,4,1]
local i,j, subword, totalSum, leftSubword, leftProd, rightSubword, rightProd;
totalSum := matrix(n^2,n^2,0);
for i from 1 to k do # construct tensor product sum
if (modp(i,2) = 1) then
for j from 0 to W[i]-1 do
leftSubword := W; # form left subword vector
leftSubword[i] := j;
leftProd := evalword(X,B,leftSubword,n,i,1);
rightSubword := W; # form right subword vector
rightSubword[i] := W[i]-j-1;
if ( W[i]-j-1 = 0 ) then
rightProd := transpose(evalword(X,B,W,n,k,i+1));
else
rightProd := transpose(evalword(X,B,rightSubword,n,k,i));
end;
totalSum := evalm(totalSum) + evalm(kronprod(rightProd,leftProd));
od;
end;
od;
return(evalm(totalSum));
end; constructM := proc(n) # Finds the matrix M
local i,j,l,k,m, M,N,alpha,beta;
m := n*(n+1)/2;
M := matrix(m,n*n,0); N := matrix(m*n*n,m*n*n,0);
for i from 1 to n do
for j from 1 to n do
for l from 1 to n do
for k from l to n do
alpha := n*(j-1)+i;
beta := (2*n-l)*(l-1)/2+k;
if (i = k and j = l) then M[beta,alpha] := 1; end;
od;
od;
od;
od;
return M;
end;constructN := proc(n) # Finds the matrix N
local i,j,l,k,m, N,alpha,beta;
m := n*(n+1)/2;
N := matrix(n*n,m,0);
for i from 1 to n do
for j from 1 to n do
for l from 1 to n do
for k from l to n do
alpha := n*(j-1)+i;
beta := (2*n-l)*(l-1)/2+k;
if ((i = k and j = l) or (i = l and j = k)) then N[alpha,beta] := 1; end;
od;
od;
od;
od;
return N;
end;B := diag(a,b);X := matrix(2,2,[x,y,y,(1+y^2)/x]);bigjacobian := wordjacobian(B,X,[1,1,2,3,2,1,1],2,7):smalljacobian := evalm(constructM(2)&*bigjacobian&*constructN(2)):factor(det(%));# all postive coefficients