# Maple file for the computations in # C.Hillar, S. Armstrong, Solvability of symmetric word equations in positive definite # letters, Journal of the London Mathematical Society, 76 (2007), no. 3, 777-796. with(linalg): kronprod := proc (A, B) 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([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; A1 := matrix(3,3,[1,20,210,20,402,4240,210,4240,44903]); B1 := matrix(3,3,[36501,-3820,190,-3820,401,-20,190,-20,1]); bigjacobian := wordjacobian(B1,A1,[1,1,2,3,2,1,1],3,7); smalljacobian := evalm(constructM(3)&*bigjacobian&*constructN(3)); det(smalljacobian);