# 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) 
  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;


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 positive coefficients