with(LinearAlgebra): with(combinat):

### some necessary procedures ###
reverse := proc(L::list) return [seq(L[nops(L)-i+1], i = 1 .. nops(L))] end proc:

mypartition := proc(d ::integer) 
	local X;
	X := reverse(partition(d)):
	return [seq(reverse(X[i]),i=1..numbpart(d))];
end proc:

myconjpart := proc (L::list) local preK, K; 
	description "I needed to have a conjugate partition function that is in decreasing order."; 
	return reverse(conjpart(reverse(L)));
end proc:

### procedures specific to this task ###

makeDets:=proc(a,LL::list,mu::list)
	local L:
	description "this procedure makes a product of determinants of sizes determined by a partition LL. The first index of each column is twisted by a permutation mu.";
	L:=myconjpart(LL):
	return `*`(seq(Determinant(Matrix([seq([seq(a[op(j+( `+`(seq(op(p,L),p=1..k-1))),mu),i],j=1 .. op(k,L))], i = 1 .. op(k,L))])) ,k=1..nops(L)))
end proc:


makeUnsymmetric:=proc(J::list,K::list)
description "this procedure takes in a list of partitions J and a list of permutations K and produces the unsymmetrized (and factored!) tensor";
local alpha;
	if(nops(J)<> nops(K)) then return "uneven";
	else
		alpha:= [a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z]:
		return `*`(seq(makeDets(alpha[i],J[i],K[i]),i=1..nops(J)))
	fi:
end proc:


unfactor:= proc(X,degree,L::list)
description "X is the tensor, d is the degree, L is the list of dimensions of the vector spaces":
local temp,temp2,p,i;

 if nops(L) >8 then return "too many factors"; fi:

if nops(L) = 2 then
temp2 := X; 

for p to degree do 
 temp := 0; 
  for i[1] from 0 to op(1,L)-1 do for i[2] from 0 to op(2,L)-1 do
   temp := coeff(coeff(temp2, a[p, i[1]+1]), b[p, i[2]+1])*Z[[ seq(i[p],p=1..nops(L) )]]+temp 
  end do end do; 
 temp2 := temp; #print(nops(temp2)) 
end do; 
return temp;
fi:

if nops(L) = 3 then
temp2 := X; 
for p to degree do 
 temp := 0; 
  for i[1] from 0 to op(1,L) do for i[2] from 0 to op(2,L) do for i[3] from 0 to op(3,L) do 
   temp := coeff(coeff(coeff(temp2, a[p, i[1]+1]), b[p, i[2]+1]), c[p, i[3]+1])*Z[[ seq(i[p],p=1..nops(L) )]]+temp 
  end do end do end do; 

 temp2 := temp; #print(nops(temp2)) 
end do; 
return temp;
fi:

if nops(L) = 4 then
temp2 := X; 
for p to degree do 
 temp := 0; 
  for i[1] from 0 to op(1,L) do for i[2] from 0 to op(2,L) do for i[3] from 0 to op(3,L) do for i[4] from 0 to op(4,L) do 
   temp := coeff(coeff(coeff(coeff(temp2, a[p, i[1]+1]), b[p, i[2]+1]), c[p, i[3]+1]),d[p, i[4]+1])*Z[[ seq(i[p],p=1..nops(L) )]]+temp ;
  end do end do end do end do; 
 temp2 := temp; #print(nops(temp2)) 
end do; 
return temp;
fi:

if nops(L) = 5 then
temp2 := X; 
for p to degree do 
 temp := 0; 
  for i[1] from 0 to op(1,L) do for i[2] from 0 to op(2,L) do for i[3] from 0 to op(3,L) do for i[4] from 0 to op(4,L) do for i[5] from 0 to op(5,L) do 

   temp := coeff(coeff(coeff(coeff(coeff(temp2, a[p, i[1]+1]), b[p, i[2]+1]), c[p, i[3]+1]),d[p, i[4]+1]),e[p,i[5]+1])*Z[[ seq(i[p],p=1..nops(L) )]]+temp 

  end do end do end do end do end do; 
 temp2 := temp; #print(nops(temp2)) 
end do; 
return temp;
fi:

if nops(L) = 6 then
temp2 := X; 
for p to degree do 
 temp := 0; 
  for i[1] from 0 to op(1,L) do for i[2] from 0 to op(2,L) do for i[3] from 0 to op(3,L) do for i[4] from 0 to op(4,L) do for i[5] from 0 to op(5,L) do for i[6] from 0 to op(6,L) do 

   temp := coeff(coeff(coeff(coeff(coeff(coeff(temp2, a[p, i[1]+1]), b[p, i[2]+1]), c[p, i[3]+1]),d[p, i[4]+1]),e[p,i[5]+1]),f[i[6]+1])*Z[[ seq(i[p],p=1..nops(L) )]]+temp 

  end do end do end do end do end do end do; 
 temp2 := temp; #print(nops(temp2)) 
end do; 
return temp;
fi:

if nops(L) = 7 then
temp2 := X; 
for p to d do 
 temp := 0; 
  for i[1] from 0 to op(1,L) do for i[2] from 0 to op(2,L) do for i[3] from 0 to op(3,L) do for i[4] from 0 to op(4,L) do for i[5] from 0 to op(5,L) do for i[6] from 0 to op(6,L) do for i[7] from 0 to op(7,L) do

   temp := coeff(coeff(coeff(coeff(coeff(coeff(coeff(temp2, a[p, i[1]+1]), b[p, i[2]+1]), c[p, i[3]+1]),d[p, i[4]+1]),e[p,i[5]+1]),f[i[6]+1]),g[i[7]+1])*Z[[ seq(i[p],p=1..nops(L) )]]+temp 

  end do end do end do end do end do end do end do; 
 temp2 := temp; #print(nops(temp2)) 
end do; 
return temp;
fi:

end proc:



## examples ##

#T:=makeUnsymmetric([[2,1,1,1],[3,1,1],[2,1,1,1]],[[1,2,3,4,5],[1,5,3,4,2],[1,4,5,2,3]]):
#unfactor(T,5,[3,2,3]): nops(expand(%));


#T:=makeUnsymmetric([[2,2,2],[2,2,2],[3,1,1,1]],[[1,2,3,4,5,6],[1,5,3,4,2,6],[1,4,5,2,3,6]]):
#unfactor(T,6,[2,2,3]): nops(expand(%));


#T:=makeUnsymmetric([[3,1,1],[3,1,1],[2,2,1]],[[1,2,3,4,5],[1,4,5,2,3],[1,2,3,4,5]]):
#unfactor(T,5,[2,2,2]): nops(expand(%));


#T:=makeUnsymmetric([[3,1,1],[3,1,1],[2,2,1]],[[1,2,3,4,5],[1,2,4,3,5],[1,3,5,2,4]]):
#unfactor(T,5,[2,2,2]): nops(expand(%));


## Here is an example where there is multiplicity greater than 1 and we already have guessed the correct permutations to give linearly independent elements of the highest weight space ##

sigma:=[1,2,3,5,6,4]: 
tau := [3, 4, 5, 1, 2, 6]; mu := [1, 4, 5, 6, 2, 3]; 
T1:=makeUnsymmetric([[3,2,1],[3,2,1],[3,1,1,1]],[sigma,tau,mu]):
P1:=unfactor(T1,6,[2,2,3]): nops(expand(%));

tau := [3, 4, 5, 1, 2, 6]; mu := [2, 3, 5, 6, 1, 4]; 
T2:=makeUnsymmetric([[3,2,1],[3,2,1],[3,1,1,1]],[sigma,tau,mu]):
P2:=unfactor(T2,6,[2,2,3]): nops(expand(%));

tau := [3, 4, 5, 1, 2, 6]; mu := [2, 3, 4, 5, 1, 6]; 
T3:=makeUnsymmetric([[3,2,1],[3,2,1],[3,1,1,1]],[sigma,tau,mu]):
P3:=unfactor(T3,6,[2,2,3]): nops(expand(%));

tau := [3, 4, 6, 1, 2, 5]; mu := [2, 3, 4, 5, 1, 6]; 
T4:=makeUnsymmetric([[3,2,1],[3,2,1],[3,1,1,1]],[sigma,tau,mu]):
P4:=unfactor(T4,6,[2,2,3]): nops(expand(%));

PP:= ss*P1+tt*P2+uu*P3+vv*P4:

#Here we check if any linear combination of these polynomials vanish on the 4th secant variety of the Segre product of P^2 x P^2 x P^3

for count from 1 to 4 do
 for j to 6 do 


mysegrepoint[j] := expand(unfactor(`+`(seq(`+`(seq(a[1,k]*U[[k-1,i]],k=1..3))*`+`(seq(b[1,k]*V[[k-1,i]],k=1..3))*`+`(seq(c[1,k]*W[[k-1,i]],k=1..4)),i=1..j)), 1, [2,2,3] ));



end do:
 for j to 6 do mysegrerandomizer := {}:
  for i to j do 
  x:= RandomVector(3); y := RandomVector(3); z := RandomVector(4); 
  mysegrerandomizer := `union`(mysegrerandomizer, {seq(U[[k-1,i]] = x[k], k = 1 .. 3), seq(V[[k-1,i]] = y[k], k = 1 .. 3), seq(W[[k-1,i]] = z[k], k = 1 .. 4)}) 
end do: 
x := 'x'; y := 'y'; z := 'z':
end do:
for j to 6 do 
 myrandomsegrepoint[j] := subs(mysegrerandomizer, mysegrepoint[j])
end do:
for p to 6 do 
mysegreevaluator || p := {seq(seq(seq(Z[[i, j, k]] = coeff(mysegrepoint[p], Z[[i, j, k]]), i = 0 .. 2), j = 0 .. 2), k = 0 .. 3)}:

myrandomsegreevaluator || p := {seq(seq(seq(Z[[i, j, k]] = coeff(myrandomsegrepoint[p], Z[[ i,j,k]]), i = 0 .. 2), j = 0 .. 2), k = 0 .. 3)}: 

end do:

val||count:=subs(myrandomsegreevaluator || 4, expand(PP));
od:


solve({val1,val2,val3,val4});