q:=Indeterminate(Rationals,1);
SetName(q,"q");
loud:=true;
LogTo("spechthom.log");LogTo();

Print("________________________________________________________________________________________\n");
Print("\nThe magic homomorphism calculator 2.2.9                                   16th May 2012\n");
Print("________________________________________________________________________________________\n");
Print("\nuseful functions:\n");
Print("   workouthoms(lambda,mu,e,p)     (computes hom space between Specht modules)\n");
Print("   makesemi(T[,e,p])              (rewrites t as a linear combination of semistandards)\n");
Print("   makesemil(l,[,e,p])            (the same for l a linear combination of homs)\n");
Print("   compsi(T,d,t[,e,p])            (composes the hom with psi_{d,t} and semistandardises)\n");
Print("   compsil(l,d,t[,e,p])           (the same for l a linear combination of homs)\n");
Print("   quiet()                        (suppress printed output)\n");
Print("   noisy()                        (restore printed output)\n");
Print("________________________________________________________________________________________\n");

quiet:=function()
loud:=false;
end;

noisy:=function()
loud:=true;
end;

################################ simple partition functions

conjugate:=function(a)
#assumes a is a partition
if a=[] then
	return [];
else
	return List([1..a[1]],i->Number(a,j->j>=i));
fi;
end;

samecore:=function(la,mu,e)
#tells whether la and mu have the same e-core
local m,j,k;
if e=0 then
	return la=mu;
else
	m:=List([1..e],i->0);
	for j in [1..Maximum(Length(la),Length(mu))] do
		if j<=Length(la) then k:=(la[j]-j) mod e; else k:=(-j) mod e;fi;
		m[k+1]:=m[k+1]+1;
		if j<=Length(mu) then k:=(mu[j]-j) mod e; else k:=(-j) mod e;fi;
		m[k+1]:=m[k+1]-1;
	od;
	return Number(m,i->i<>0)=0;
fi;
end;

dominates:=function(la,mu)
#assumes the two partitions have the same size
if Length(la)>Length(mu) then
	return false;
else
	return Number([1..Length(la)],i->Sum(la{[1..i]})<Sum(mu{[1..i]}))=0;
fi;
end;

################################ quantum integers and binomial coefficients

specialq:=function(e,p)
#e is at least 2, p is 0 or prime, and if p is positive then e=p or is not divisible by p
local n;
if p=0 then
	if e=0 then
		return 1;
	else
		return E(e);
	fi;
elif p=e then
	return Z(p)^0;
elif e=0 then
	return Indeterminate(GF(p),"q");
else
	n:=p;
	while (n-1) mod e<>0 do
		n:=n*p;
	od;
	return Z(n)^((n-1)/e);
fi;
end;

qinteger:=function(m)
return (q^m-q^0)/(q-q^0);
end;

qfactorial:=function(m)
return Product(List([1..m],n->qinteger(n)));
end;

qbc:=function(n,r)
#evalutates the quantum binomial coefficient n choose r with the 1+q+..+q^(m-1) convention
return qfactorial(n)/(qfactorial(r)*qfactorial(n-r));
end;

evallist:=function(li,qq)
local lis;
lis:=Filtered(List(li,m->[m[1],Value(m[2],qq)]),m->m[2]<>0*qq);
return lis;
end;

################################################################# tableaux

usualtab:=function(t)
#turns a tableau in matrix form into one in tableau form
return List(t,m->Concatenation(List([1..Length(m)],i->List([1..m[i]],j->i))));
end;

unusualtab:=function(t)
#turns a tableau in tableau form into one in typical matrix form (so row sums are la, columns sums are mu)
return List(t,m->List([1..Maximum(m)],i->Number(m,j->j=i)));
end;

sstd:=function(la,mu)
#computes all semi-standard la-tableaux of type mu.  Assumes la is a partition and mu a composition
#*********returns a tableau as a matrix, which the i,j-entry is the number of is in row j (so the row sums are mu, and the column sums are la)**********
local ss,i,j,t,l,m,good;
if Sum(la)<>Sum(mu) then
	return [];
elif Sum(la)=0 then
	return [[]];
else
	l:=Length(la);
	m:=Length(mu);
	if l=1 then
		return [List(mu,k->[k])];
	else
		ss:=[];
		t:=List([1..m],x->List([1..l],y->0));
		i:=1;
		j:=1;
		while i>0 do
			if j<i and j<l then
				t[i][j]:=t[i][j]+1;
				good:=true;
				if Sum(t[i])>mu[i] or Sum(TransposedMat(t)[j])>la[j] then
					good:=false;
				elif j>1 then #semistandard condition
					if Sum(TransposedMat(t)[j]{[1..i]})>Sum(TransposedMat(t)[j-1]{[1..i-1]}) then good:=false;fi;
				fi;
				if good=false then
					t[i][j]:=0;
					j:=j-1;
					if j=0 then
						i:=i-1;
						if i>0 then
							t[i][Minimum([i,l])]:=0;
							j:=Minimum([i,l])-1;
							if j=0 then i:=0;fi;
						fi;
					fi;
				else
					j:=j+1;
					t[i][j]:=-1;
				fi;
			else
				if i=1 then t[i][j]:=mu[1]; else t[i][j]:=mu[i]-Sum(t[i]{[1..j-1]});fi;
				good:=true;
				if Sum(TransposedMat(t)[j])>la[j] then
					good:=false;
				elif j>1 then #semistandard condition
					if Sum(TransposedMat(t)[j]{[1..i]})>Sum(TransposedMat(t)[j-1]{[1..i-1]}) then good:=false;fi;
				fi;
				if good=false then
					t[i][j]:=0;
					j:=j-1;
					if j=0 then
						i:=i-1;
						if i>0 then
							t[i][Minimum([i,l])]:=0;
							j:=Minimum(i,l)-1;
							if j=0 then i:=0;fi;
						fi;
					fi;
				else
					i:=i+1;
					if i=m+1 then
						Add(ss,ShallowCopy(List(t,u->ShallowCopy(u))));
						i:=m;
						t[i][Minimum([i,l])]:=0;
						j:=Minimum([i,l])-1;
					else
						j:=1;
						t[i][j]:=-1;
					fi;
				fi;
			fi;
		od;
		return ss;
	fi;
fi;
end;

issemi:=function(t)
#tells us whether the tableau t is semistandard (in usual matrix form (i.e. row i is row i), may have trailing 0s etc)
local iss,i,j;
iss:=true;
for i in [2..Length(t)] do
	for j in [1..Length(t[i])] do
		if j<i then
			if t[i][j]>0 then iss:=false;fi;
		else
			if Sum(t[i]{[1..j]})>Sum(t[i-1]{[1..Minimum([j-1,Length(t[i-1])])]}) then iss:=false;fi;
		fi;
	od;
od;
return iss;		
end;

isreversesemi:=function(t)
#tells us whether the tableau t is reverse semistandard (in usual matrix form, may have trailing 0s etc)
local iss,i,j;
iss:=true;
for i in [2..Length(t)] do
	for j in [1..Length(t[i])] do
		if Sum(t[i]{[j..Length(t[i])]})>Sum(t[i-1]{[j+1..Length(t[i-1])]}) then iss:=false;fi;
	od;
od;
return iss;		
end;

allconjsemi:=function(la,mu)
return List(sstd(la,mu),t->List(TransposedMat(t),k->k{[1..Maximum(Filtered([1..Length(k)],j->k[j]>0))]}));
end;

################################################################### semistandardising

subcomps:=function(mu,c)
#mu is a composition and this returns all subcompositions of size c
local s,sc,i;
sc:=[];
s:=List(mu,k->0);
s[1]:=-1;
i:=1;
while i>0 do
	s[i]:=s[i]+1;
	if s[i]>mu[i] then
		s[i]:=0;
		i:=i-1;
	elif Sum(s)=c then
		Add(sc,ShallowCopy(s));
		s[i]:=0;
		i:=i-1;
	else
		if Sum(mu{[i+1..Length(mu)]})>=c-Sum(s) then
			i:=i+1;
			s[i]:=-1;
		fi;
	fi;
od;
return sc;
end;

coee:=function(c,d)
#c and d are contents, and this gives the coefficient needed for the garnir relation
local co,i,j;
co:=q^0;
for i in [1..Minimum([Length(c),Length(d)])] do
	co:=co*qbc(c[i]+d[i],c[i]);
od;
for i in [2..Length(c)] do
	for j in [1..Minimum([i-1,Length(d)])] do
		co:=co*q^(c[i]*d[j]);
	od;
od;
return co;
end;

abbrev:=function(c)#strips off trailing zeroes
local f;
f:=First([Length(c),Length(c)-1..1],i->c[i]>0);
if f=fail then
	return [];
else
	return c{[1..f]};
fi;
end;

homlist:=[];
#homlist is a global variable which will have entries [tableau in matrix form, coefficient], and the coefficient can have q specialised

condense:=function()
#makes sure homlist is collected and sorted
local h;
h:=Set(List(homlist,k->k[1]));
homlist:=List(h,k->[k,Sum(Filtered(homlist,j->j[1]=k),j->j[2])]);
homlist:=Filtered(homlist,k->k[2]<>k[2]+k[2]);
end;

semiststep:=function(arg)
#carries out one step in the algorithm for semistandardising homlist
#first argument is the element of homlist which is going to be addressed - we assume this really is an element of the list, and that its first element is non-semistandard
# optional second and third arguments are e,p
#assumes cleaner has been carried out
#carries out condense at the end
local qq,f,r,u,i,topc,botc,midc,c,mcoe,sc,newu,newc,pp;
f:=arg[1];
if IsBound(arg[3]) then qq:=specialq(arg[2],arg[3]);else qq:=q;fi;
r:=First([1..Length(f[1])-1],r->not issemi(f[1]{[r,r+1]}));
u:=[ShallowCopy(f[1][r]),ShallowCopy(f[1][r+1])];#u is a copy of the two-row subtableau
i:=First([1..Length(u[2])],i->Sum(u[2]{[1..i]})>Sum(u[1]{[1..i-1]}));#i is the *value* in the bottom row (n.b. not necessarily the top!) of the first column where semistandardness goes wrong (the value on the top row is at least i, and is the first thing in the top row which is at least i)
topc:=u[1]{[1..i-1]};#content which is being fixed at the top
if Length(u[2])>i then botc:=Concatenation(List([1..i],j->0),u[2]{[i+1..Length(u[2])]});else botc:=[];fi;
c:=Sum(u[2])-Sum(botc);
midc:=u[1]+u[2]-topc-botc;#remaining content
mcoe:=Value(coee(topc,u[1]-topc)*coee(u[2]-botc,botc),qq);
for sc in subcomps(midc,c) do
	if abbrev(sc)<>abbrev(u[2]-botc) then
		newu:=Concatenation(f[1]{[1..r-1]},[abbrev(topc+midc-sc),abbrev(botc+sc)],f[1]{[r+2..Length(f[1])]});
		newc:=f[2]*Value(coee(topc,midc-sc)*coee(sc,botc),qq);
		if newc<>qq-qq then
			pp:=Position(List(homlist,k->k[1]),newu);
			if pp=fail then
				Add(homlist,[ShallowCopy(newu),ShallowCopy(-newc/mcoe)]);
			else
				homlist[pp][2]:=homlist[pp][2]-newc/mcoe;
			fi;
		fi;
	fi;
od;
Unbind(homlist[Position(homlist,f)]);
condense();
end;

rsemiststep:=function(arg)
#carries out one step in the algorithm for semistandardising homlist
#first argument is the element of homlist which is going to be addressed - we assume this really is an element of the list, and that its first element is non-reverse-semistandard
# optional second and third arguments are e,p
#assumes cleaner has been carried out
#carries out condense at the end
local qq,f,r,u,i,topc,botc,midc,c,mcoe,sc,newu,newc,pp;
f:=arg[1];
if IsBound(arg[3]) then qq:=specialq(arg[2],arg[3]);else qq:=q;fi;
r:=First([1..Length(f[1])-1],r->not isreversesemi(f[1]{[r,r+1]}));
u:=[ShallowCopy(f[1][r]),ShallowCopy(f[1][r+1])];#u is a copy of the two-row subtableau
i:=First([Length(u[2]),Length(u[2])-1..1],i->Sum(u[2]{[i..Length(u[2])]})>Sum(u[1]{[i+1..Length(u[1])]}));#i is the *value* in the bottom row (n.b. not necessarily the top!) of the first column where reverse-semistandardness goes wrong (the value on the top row is at most i, and is the first thing in the top row which is at most i)
botc:=u[2]{[1..i-1]};#content which is being fixed at the bottom
if Length(u[1])>i then topc:=Concatenation(List([1..i],j->0),u[1]{[i+1..Length(u[1])]});else topc:=[];fi;
c:=Sum(u[2])-Sum(botc);
midc:=u[1]+u[2]-topc-botc;#remaining content
mcoe:=Value(coee(topc,u[1]-topc)*coee(u[2]-botc,botc),qq);
for sc in subcomps(midc,c) do
	if abbrev(sc)<>abbrev(u[2]-botc) then
		newu:=Concatenation(f[1]{[1..r-1]},[abbrev(topc+midc-sc),abbrev(botc+sc)],f[1]{[r+2..Length(f[1])]});
		newc:=f[2]*Value(coee(topc,midc-sc)*coee(sc,botc),qq);
		if newc<>qq-qq then
			pp:=Position(List(homlist,k->k[1]),newu);
			if pp=fail then
				Add(homlist,[ShallowCopy(newu),ShallowCopy(-newc/mcoe)]);
			else
				homlist[pp][2]:=homlist[pp][2]-newc/mcoe;
			fi;
		fi;
	fi;
od;
Unbind(homlist[Position(homlist,f)]);
condense();
end;

cleaner:=function()
#removes trailing zeroes from rows of tableaux appearing in homlist
homlist:=List(homlist,k->[unusualtab(usualtab(k[1])),k[2]]);
condense();
end;

semist:=function(arg)
#semistandardises homlist - two optional arguments are e,p
cleaner();
while First(homlist,k->not issemi(k[1]))<>fail do
	if IsBound(arg[2]) then semiststep(First(homlist,k->not issemi(k[1])),arg[1],arg[2]);else semiststep(First(homlist,k->not issemi(k[1])));fi;
od;
end;

rsemist:=function(arg)
#semistandardises homlist - two optional arguments are e,p
cleaner();
while First(homlist,k->not isreversesemi(k[1]))<>fail do
	if IsBound(arg[2]) then rsemiststep(First(homlist,k->not isreversesemi(k[1])),arg[1],arg[2]);else rsemiststep(First(homlist,k->not isreversesemi(k[1])));fi;
od;
end;

domcomp:=function(mu,nu)
#tells us whether the composition mu dominates the composition nu in the usual sense
#assumes they have the same sum
#can have zeroes anywhere
return First([1..Length(mu)-1],i->Sum(mu{[1..i]})<Sum(nu{[1..Minimum([i,Length(nu)])]}))=fail;
end;

domtab:=function(t,u)
#t and u are row-standard tableaux in typical matrix form
#tells us whether t dominates u in the sense of having smaller entries higher up
#assumes they have the same shape and content
return First([1..Length(t)-1],i->not domcomp(Sum(t{[1..i]}),Sum(u{[1..i]})))=fail;
end;

isnonzero:=function(arg)
#tells us whether the hom represented by homlist is non-zero, wuthout necessarily semistandardising it all the way
#optional arguments are e,p
local done,i,j,inz,f;
cleaner();
done:=false;
while done=false do
	i:=First([1..Length(homlist)],i->issemi(homlist[i][1]));
	if i<>fail then
		j:=First([1..i-1],j->domtab(homlist[i][1],homlist[j][1]));
		if j<>fail then
			f:=StructuralCopy(homlist[j]);
			if IsBound(arg[2]) then semiststep(f,arg[1],arg[2]);else semiststep(f);fi;
		else
			inz:=true;
			done:=true;
		fi;
	elif homlist=[] then
		inz:=false;
		done:=true;
	else
		f:=StructuralCopy(homlist[1]);
		if IsBound(arg[2]) then semiststep(f,arg[1],arg[2]);else semiststep(f);fi;
	fi;
od;
return inz;
end;

ssiseu:=function(arg)
#given a tableau in usual tableau form, makes it semistandard
local qq;
if IsBound(arg[3]) then qq:=specialq(arg[2],arg[3]);else qq:=q;fi;
homlist:=[[unusualtab(arg[1]),qq^0]];
if IsBound(arg[3]) then semist(arg[2],arg[3]);else semist();fi;
return ShallowCopy(homlist);
end;

rssiseu:=function(arg)
#given a tableau in usual tableau form, makes it semistandard
local qq;
if IsBound(arg[3]) then qq:=specialq(arg[2],arg[3]);else qq:=q;fi;
homlist:=[[unusualtab(arg[1]),qq^0]];
if IsBound(arg[3]) then rsemist(arg[2],arg[3]);else semist();fi;
return ShallowCopy(homlist);
end;

makesemi:=function(arg)
#############user function##################
#given a tableau in the usual tableau form, expresses it as a linear combination of semistandards
############################################
if IsBound(arg[3]) then
	return List(ssiseu(arg[1],arg[2],arg[3]),m->[usualtab(m[1]),m[2]]);
else
	return List(ssiseu(arg[1]),m->[usualtab(m[1]),m[2]]);
fi;
end;

makesemil:=function(arg)
##############user function#################
#given a linear combination of tableaux in the usual tableau form, expresses it as a linear combination of semistandards
############################################
local qq,k;
if IsBound(arg[3]) then
	qq:=specialq(arg[2],arg[3]);
	homlist:=[];
	for k in arg[1] do
		if IsLaurentPolynomial(k[2]) then
			Add(homlist,[unusualtab(k[1]),Value(k[2],qq)]);
		else
			Add(homlist,[unusualtab(k[1]),k[2]]);
		fi;
	od;
	semist(arg[2],arg[3]);
else
	homlist:=List(arg[1],k->[unusualtab(k[1]),k[2]]);
	semist();
fi;
return List(homlist,m->[usualtab(m[1]),m[2]]);
end;

rev:=function(t)
return List(t,j->Reversed(j));
end;

makereversesemi:=function(arg)
#############user function##################
#given a tableau in the usual tableau form, expresses it as a linear combination of semistandards
############################################
if IsBound(arg[3]) then
	return List(rssiseu(arg[1],arg[2],arg[3]),m->[rev(usualtab(m[1])),m[2]]);
else
	return List(rssiseu(arg[1]),m->[rev(usualtab(m[1])),m[2]]);
fi;
end;

makereversesemil:=function(arg)
##############user function#################
#given a linear combination of tableaux in the usual tableau form, expresses it as a linear combination of semistandards
############################################
local qq,k;
if IsBound(arg[3]) then
	qq:=specialq(arg[2],arg[3]);
	homlist:=[];
	for k in arg[1] do
		if IsLaurentPolynomial(k[2]) then
			Add(homlist,[unusualtab(k[1]),Value(k[2],qq)]);
		else
			Add(homlist,[unusualtab(k[1]),k[2]]);
		fi;
	od;
	rsemist(arg[2],arg[3]);
else
	homlist:=List(arg[1],k->[unusualtab(k[1]),k[2]]);
	rsemist();
fi;
return List(homlist,m->[rev(usualtab(m[1])),m[2]]);
end;

nonzero:=function(arg)
#arg consists of a tableau in usual form, possibly followed by e,p, and this tells us whether the corresponding homomorphism is non-zero
local qq;
if IsBound(arg[3]) then qq:=specialq(arg[2],arg[3]);else qq:=q;fi;
homlist:=[[unusualtab(arg[1]),qq^0]];
if IsBound(arg[3]) then return isnonzero(arg[2],arg[3]);else return isnonzero();fi;
end;

nonzerol:=function(arg)
homlist:=List(arg[1],k->[unusualtab(k[1]),k[2]]);
#Print(homlist,"\n\n\n",isnonzero(2,2));
if IsBound(arg[3]) then return isnonzero(arg[2],arg[3]);else return isnonzero();fi;
end;

################################################################# applying psi-homs

psidt:=function(s,d,t)
#takes a tableau s and composes with the map psi(d,t) (q indeterminate).  returns a list of pairs [tableau,coefficient].
#n.b. s is in matrix form
local ss,l,i,ps,g,u,j,coe;
ss:=ShallowCopy(List(s,m->ShallowCopy(m)));
l:=Length(ss);
for i in [1..l] do
	while Length(ss[i])<d+1 do
		Add(ss[i],0);
	od;
od;
ps:=[];
g:=[-1];
i:=1;
while i>0 do
	g[i]:=g[i]+1;
	if g[i]>ss[i][d+1] then
		Unbind(g[i]);
		i:=i-1;
	elif Sum(g)=t then
		u:=ShallowCopy(List(ss,m->ShallowCopy(m)));
		for j in [1..Length(g)] do
			u[j][d]:=u[j][d]+g[j];
			u[j][d+1]:=u[j][d+1]-g[j];
		od;
		coe:=Product(List([1..Length(g)],m->qbc(u[m][d],g[m])))*q^(Sum(List([1..Minimum([l-1,Length(g)])],m->g[m]*Sum(List([m+1..l],n->ss[n][d])))));
		for j in [1..l] do
			while u[j][Length(u[j])]=0 do Unbind(u[j][Length(u[j])]);od;
		od;
		Add(ps,[ShallowCopy(List(u,m->ShallowCopy(m))),coe]);
		Unbind(g[i]);
		i:=i-1;
	else
		i:=i+1;
		if i>l then i:=i-1;else Add(g,-1);fi;
	fi;
od;
return ps;
end;

psidtu:=function(s,d,t)
#takes a tableau s and composes with the map psi(d,t) (q indeterminate).  returns a list of pairs [tableau,coefficient].
#n.b. s is in the 'usual' form
local ss,l,i,ps,g,u,j,coe;
ss:=unusualtab(ShallowCopy(List(s,m->ShallowCopy(m))));
l:=Length(ss);
for i in [1..l] do
	while Length(ss[i])<d+1 do
		Add(ss[i],0);
	od;
od;
ps:=[];
g:=[-1];
i:=1;
while i>0 do
	g[i]:=g[i]+1;
	if g[i]>ss[i][d+1] then
		Unbind(g[i]);
		i:=i-1;
	elif Sum(g)=t then
		u:=ShallowCopy(List(ss,m->ShallowCopy(m)));
		for j in [1..Length(g)] do
			u[j][d]:=u[j][d]+g[j];
			u[j][d+1]:=u[j][d+1]-g[j];
		od;
		coe:=Product(List([1..Length(g)],m->qbc(u[m][d],g[m])))*q^(Sum(List([1..Minimum([l-1,Length(g)])],m->g[m]*Sum(List([m+1..l],n->ss[n][d])))));
		for j in [1..l] do
			while u[j][Length(u[j])]=0 do Unbind(u[j][Length(u[j])]);od;
		od;
		Add(ps,[ShallowCopy(List(u,m->ShallowCopy(m))),coe]);
		Unbind(g[i]);
		i:=i-1;
	else
		i:=i+1;
		if i>l then i:=i-1;else Add(g,-1);fi;
	fi;
od;
return List(ps,m->[usualtab(m[1]),m[2]]);
end;

psilist:=function(arg)
# optional fourth and fifth arguments are e,p
# tableaux in matrix form, coefficients unspecialised
# this semistandardises as well
################################################################
local qq,mm,m;
if IsBound(arg[5]) then
	qq:=specialq(arg[4],arg[5]);
	homlist:=[];
	for m in arg[1] do
		if IsLaurentPolynomial(m[2]) then mm:=Value(m[2],qq);else mm:=m[2];fi;
		Append(homlist,List(psidt(m[1],arg[2],arg[3]),k->[k[1],Value(k[2],qq)*mm]));
	od;
	semist(arg[4],arg[5]);
	return List(homlist,j->[j[1],j[2]]);
else
	homlist:=Concatenation(List(arg[1],m->List(psidt(m[1],arg[2],arg[3]),k->[k[1],k[2]*m[2]])));
	semist();
	return List(homlist,j->[j[1],j[2]]);
fi;
end;

psilistnz:=function(arg)
# optional fourth and fifth arguments are e,p
# tableaux in matrix form, coefficients unspecialised
# this semistandardises as well
################################################################
local qq,mm,m;
if IsBound(arg[5]) then
	qq:=specialq(arg[4],arg[5]);
	homlist:=[];
	for m in arg[1] do
		if IsLaurentPolynomial(m[2]) then mm:=Value(m[2],qq);else mm:=m[2];fi;
		Append(homlist,List(psidt(m[1],arg[2],arg[3]),k->[k[1],Value(k[2],qq)*mm]));
	od;
	return isnonzero(arg[4],arg[5]);
else
	homlist:=Concatenation(List(arg[1],m->List(psidt(m[1],arg[2],arg[3]),k->[k[1],k[2]*m[2]])));
	return isnonzero();
fi;
end;

psiu:=function(arg)
#input and output tabs are in usual form
local l;
l:=arg[1];
if IsBound(arg[4]) then
	return List(psilist(List(l,m->[unusualtab(m[1]),m[2]]),arg[2],arg[3],arg[4],arg[5]),k->[usualtab(k[1]),k[2]]);
else
	return List(psilist(List(l,m->[unusualtab(m[1]),m[2]]),arg[2],arg[3]),k->[usualtab(k[1]),k[2]]);
fi;
end;

psiunz:=function(arg)
#input and output tabs are in usual form
local l;
l:=arg[1];
if IsBound(arg[4]) then
	return psilistnz(List(l,m->[unusualtab(m[1]),m[2]]),arg[2],arg[3],arg[4],arg[5]);
else
	return psilistnz(List(l,m->[unusualtab(m[1]),m[2]]),arg[2],arg[3]);
fi;
end;

compsil:=function(arg)
############################### user function ##################
if IsBound(arg[5]) then
	return psiu(arg[1],arg[2],arg[3],arg[4],arg[5]);
else
	return psiu(arg[1],arg[2],arg[3]);
fi;	
end;

compsi:=function(arg)
############################### user function ##################
if IsBound(arg[5]) then
	return compsil([[arg[1],1]],arg[2],arg[3],arg[4],arg[5]);
else
	return compsil([[arg[1],1]],arg[2],arg[3]);
fi;	
end;

compsilnz:=function(arg)
#like compsil, but just tells us whether the result is non-zero
if IsBound(arg[5]) then
	return psiunz(arg[1],arg[2],arg[3],arg[4],arg[5]);
else
	return psiunz(arg[1],arg[2],arg[3]);
fi;	
end;

compsinz:=function(arg)
if IsBound(arg[5]) then
	return compsilnz([[arg[1],1]],arg[2],arg[3],arg[4],arg[5]);
else
	return compsilnz([[arg[1],1]],arg[2],arg[3]);
fi;	
end;


relations:=[];
#relations will be a global variable, so that we can add to it in any routine
#we will store the relations in a memory-sensible way: the non-zero entries will be labelled by [coordinate,coefficient]
acs:=[];
#this is also a global variable, which stores the semistandard tableaux (and gradually has these knocked out if they can't contribute to a hom)

addechelon:=function(r)
#assumes relations is in reduced echelon form, and adds the new relation r
#n.b. r will be destroyed
local le,lastech,topp,midd,bott,s,i,k,l,p,ra;
if acs<>[] then
	#first make sure r is in the right form
	r:=List(Set(List(r,k->k[1])),i->[i,Sum(List(Filtered(r,k->k[1]=i),k->k[2]))]);
	r:=Filtered(r,k->k[2]<>k[2]-k[2] and IsBound(acs[k[1]]));
	le:=Length(relations);
	if le=0 then
		if r<>[] then
			relations:=[List(r,l->[l[1],l[2]/r[1][2]])];
		fi;
	else
		lastech:=relations[le][1][1];
		#start by stripping off existing relations
		i:=1;#we're going to look for a relation which takes off r[i];
		topp:=1;
		while i<=Length(r) and r[i][1]<=lastech do
			bott:=le;
			while bott-topp>1 and relations[topp][1][1]<r[i][1] and relations[bott][1][1]>r[i][1] do
				midd:=Int((topp+bott)/2);
				if relations[midd][1][1]<r[i][1] then topp:=midd;else bott:=midd;fi;
			od;
			if relations[topp][1][1]=r[i][1] then
				s:=relations[topp];
			elif relations[bott][1][1]=r[i][1] then
				s:=relations[bott];
			else
				s:=[];
			fi;
			if s=[] then
				i:=i+1;
			else
				ra:=r[i][2];
				for k in Reversed(s) do
					l:=Position(List(r,m->m[1]),k[1]);
					if l=fail then
						Add(r,[k[1],-k[2]*ra]);
					else
						r[l][2]:=r[l][2]-k[2]*ra;
					fi;
				od;
				r:=Set(Filtered(r,k->k[2]<>k[2]-k[2]));
			fi;
		od;
		if r<>[] then
			r:=List(r,l->[l[1],l[2]/r[1][2]]);
			#now strip the leader of r from the existing relations
			i:=1;
			while i<=Length(relations) and relations[i][1][1]<r[1][1] do
				p:=Position(List(relations[i],k->k[1]),r[1][1]);
				if p<>fail then
					s:=ShallowCopy(relations[i]);
					RemoveSet(relations,s);
					ra:=s[p][2];
					for k in Reversed(r) do
						l:=Position(List(s,m->m[1]),k[1]);
						if l=fail then
							Add(s,[k[1],-k[2]*ra]);
						else
							s[l][2]:=s[l][2]-k[2]*ra;
						fi;
					od;
					s:=Set(Filtered(s,k->k[2]<>k[2]-k[2]));
					if Length(s)=1 then
						Unbind(acs[s[1][1]]);
					else
						AddSet(relations,s);
						i:=i+1;
					fi;
				else
					i:=i+1;
				fi;
			od;
			if Length(r)=1 then
				Unbind(acs[r[1][1]]);
			else
				AddSet(relations,r);
			fi;
		fi;
	fi;
fi;
end;

workingouthoms:=function(la,mu,e,p)
#computes the space of homs between Specht modules S^la and S^mu at q an eth root of unity in a field of characteristic p
#p can be zero, but e is assumed non-zero
#when e=2, this is of course only those homs which are linear combinations of semistandards
local tabb,mx,maxd,alll,i,j,pivots,c,r,kv,kern,m,qq,o,rr,rels,d,t,s,y,pp,maxt,k,nu,notab;
if Sum(la)<>Sum(mu) then
	if loud then Print("Partitions should have the same size\n\n");fi;
	return [];
elif not dominates(la,mu) then
	if loud then Print(la," does not dominate ",\mu,"\n\n");fi;
	return [];
elif not samecore(la,mu,e) then
	if loud then Print(la," and ",\mu," have different ",e,"-cores\n\n");fi;
	return [];
else
	qq:=specialq(e,p);
	acs:=allconjsemi(la,mu);
	notab:=ShallowCopy(Length(String(Length(acs))));
	if loud then
		Print("\nlambda=(");
		for j in la do Print(j,",");od;
		Print("\b)                         semistandard tableaux: ",Length(acs),"\nmu=(");
		for j in mu do Print(j,",");od;
		Print("\b)\n");
	fi;
	AppendTo("spechthom.log","\nlambda=(");
	for j in [1..Length(la)-1] do AppendTo("spechthom.log",la[j],",");od;
	AppendTo("spechthom.log",la[Length(la)],")\nmu=(");
	for j in [1..Length(mu)-1] do AppendTo("spechthom.log",mu[j],",");od;
	AppendTo("spechthom.log",mu[Length(mu)],")\nsemistandard tableaux: ",Length(acs),"\n");
	#we're going to try pairs d,t one at a time, stripping out semistandards as we go to save time
	relations:=[];
	rels:=[];
	tabb:=[];
	for t in [1..mu[2]] do#[mu[2],mu[2]-1..1] do
		for d in Reversed(Filtered([1..Length(mu)-1],y->mu[y+1]>=t)) do
			nu:=ShallowCopy(mu);
			nu[d+1]:=nu[d+1]-t;
			nu[d]:=nu[d]+t;
			Sort(nu);
			nu:=Reversed(nu);
			if dominates(la,nu) then
				if acs<>[] then
					if loud then
						Print("d=",d," t=",t,": surviving tableaux: ");
						for s in [1..notab] do Print(" ");od;
						for s in String(Number(acs)) do Print("\b");od;
						Print(Number(acs),"    relations: ");
						for s in [1..notab] do Print(" ");od;
						for s in String(Length(relations)) do Print("\b");od;
						Print(Length(relations),"    hom-space dimension is at most ");
						for s in [1..notab] do Print(" ");od;
						for s in String(Number(acs)-Length(relations)) do Print("\b");od;
						Print(Number(acs)-Length(relations),"    computing tableau number \c");
					fi;
					AppendTo("spechthom.log","\nd=",d," t=",t,": surviving tableaux: ",Number(acs)," relations: ",Length(relations)," hom-space dimension is at most ",Number(acs)-Length(relations));
				fi;
				for s in Filtered([1..Length(acs)],s->IsBound(acs[s])) do
					if loud then
						Print(s,"\c");
						for i in String(s) do Print("\b \b");od;
					fi;
					homlist:=evallist(psidt(acs[s],d,t),qq);
					semist(e,p);
					for m in homlist do
						o:=Position(tabb,m[1]);
						if o=fail then
							Add(tabb,ShallowCopy(List(m[1],z->ShallowCopy(z))));
							Add(rels,[]);
							o:=Length(tabb);
						fi;
						AddSet(rels[o],[s,m[2]]);
					od;
				od;
				if loud then for s in [1..106+3*notab+Length(String(d))+Length(String(t))] do Print("\b \b\c");od;Print("\c");fi;
				#now we have our relations, and we don't need to care what tabb is any more
				tabb:=[];
				if loud then Print("d=",d," t=",t,": adding new relations: \c");fi;
				for rr in Set(rels) do
					s:=Length(Set(rels))+1-Position(Set(rels),rr);
					if loud then Print(s,"\c");fi;
					addechelon(rr);
					if loud then for i in String(s) do Print("\b \b");od;fi;
				od;
				if loud then for s in [1..29+Length(String(d))+Length(String(t))] do Print("\b \b\c");od;Print("\c");fi;
				rels:=[];
			fi;
		od;
	od;
	pivots:=List(relations,k->k[1][1]);
	kern:=[];
	for i in Difference(Filtered([1..Length(acs)],i->IsBound(acs[i])),pivots) do
		kv:=[];
		for j in [Length(acs),Length(acs)-1..1] do
			if not IsBound(acs[i]) then
				kv[j]:=qq-qq;
			elif j=i then
				kv[j]:=qq^0;
			elif not j in pivots then
				kv[j]:=0*qq;
			else
				r:=Position(pivots,j);
				kv[j]:=qq-qq;
				for k in relations[r] do
					if k[1]>j then
						kv[j]:=kv[j]-k[2]*kv[k[1]];
					fi;
				od;
			fi;
		od;
		Add(kern,List(Filtered([1..Length(acs)],j->kv[j]<>0*qq),j->[usualtab(acs[j]),kv[j]]));
	od;
	return kern;
fi;
end;

workouthoms:=function(la,mu,e,p)
#takes into account (generalised) row removal
local qq,rero,coro,larc,murc,adder,cladd,spl,i,woh,s,j,k,l,lla,lala,mumu,lm,sp,ho,co,ad,nz,wk;
qq:=specialq(e,p);
if not Sum(la)=Sum(mu) then
	if loud then Print(la," and ",mu," have different sizes!\n\n");fi;
	return [];
elif not samecore(la,mu,e) then
	if loud then Print(la," and ",mu," have different ",e,"-cores!\n\n");fi;
	return [];
elif not dominates(la,mu) then
	if loud then Print("\nhom-space is 0-dimensional\n\n");fi;
	return [];
elif la=mu then
	if loud then Print("\nhom-space is 1-dimensional\n\n");fi;
	return [[[List([1..Length(la)],i->List([1..la[i]],j->i)),qq^0]]];
else
	#apply generalised row removal
	lla:=ShallowCopy(la);
	while Length(lla)<Length(mu) do Add(lla,0);od;
	l:=Length(lla);
	sp:=Filtered([0..l],i->Sum(lla{[1..i]})=Sum(mu{[1..i]}));
	lm:=[];
	for i in [2..Length(sp)] do
		if sp[i]=sp[i-1]+1 then
			lm[i]:=[[],[],la[sp[i]]];
		else
			lala:=List([sp[i-1]+1..sp[i]],j->lla[j]-lla[sp[i]]);
			mumu:=List([sp[i-1]+1..sp[i]],j->mu[j]-lla[sp[i]]);
			while lala[Length(lala)]=0 do Unbind(lala[Length(lala)]);od;
			lm[i]:=[ShallowCopy(lala),ShallowCopy(mumu),lla[sp[i]]];
		fi;
	od;
	if First(lm,l->not samecore(l[1],l[2],e))<>fail then
		if loud then Print("\nhom-space is 0-dimensional\n\n");fi;
		return [];
	else
		nz:=true;#checks whether there are still plausibly some non-zero homs at each stage
		ho:=[];
		for i in [2..Length(sp)] do
			if nz then
				if lm[i][1]=[] then
					ho[i]:=[[[[],qq^0]]];
				else
					wk:=workingouthoms(lm[i][1],lm[i][2],e,p);
					if wk<>[] then
						ho[i]:=wk;
					else
						nz:=false;
					fi;
				fi;
			fi;
		od;
		if nz then
			#modify the tableaux
			for i in [2..Length(sp)] do
				if lm[i][1]=[] then
					ho[i][1][1][1]:=[List([1..lm[i][3]],j->sp[i])];
				else
					ad:=sp[i-1];
					co:=lm[i][3];
					for j in [1..Length(ho[i])] do#j is a basis vector
						for k in [1..Length(ho[i][j])] do#k is a pair [tab,coeff]
							for l in [1..Length(ho[i][j][k][1])] do#l is a row of a tab
								ho[i][j][k][1][l]:=Concatenation(List([1..co],z->l+ad),ho[i][j][k][1][l]+ad);
							od;
							if co>0 then
								for l in [Length(ho[i][j][k][1])+1..sp[i]-sp[i-1]] do
									Add(ho[i][j][k][1],List([1..co],z->l+ad));
								od;
							fi;
						od;
					od;
				fi;
			od;
			if loud then Print("\n\nhom-space is ",Length(Cartesian(ho{[2..Length(ho)]})),"-dimensional\n\n");fi;
			AppendTo("spechthom.log","\n\nhom-space is ",Length(Cartesian(ho{[2..Length(ho)]})),"-dimensional\n");
			return List(Cartesian(ho{[2..Length(ho)]}),i->List(Cartesian(i),j->[Concatenation(List(j,k->k[1])),Product(List(j,k->k[2]))]));			
		else
			if loud then Print("\nhom-space is 0-dimensional\n\n");fi;
			return [];
		fi;
	fi;
fi;
end;

ishom:=function(l,e,p)
#tells us whether l gives a hom, by checking all pairs d,t
#assumes the tableaux all have the same shape and type
local d,t,la,mu,ish,nu;
if l=[] then
	return true;
else
	ish:=true;
	la:=List(l[1][1],i->Length(i));
	mu:=List([1..Maximum(Concatenation(l[1][1]))],i->Number(Concatenation(l[1][1]),j->j=i));
	if Length(mu)>1 then
		for t in [mu[2],mu[2]-1..1] do
			for d in Reversed(Filtered([1..Length(mu)-1],y->mu[y+1]>=t)) do
				if ish then
					nu:=ShallowCopy(mu);
					nu[d+1]:=nu[d+1]-t;
					nu[d]:=nu[d]+t;
					Sort(nu);
					nu:=Reversed(nu);
					if dominates(la,nu) then
						ish:=not compsilnz(l,d,t,e,p);
					fi;
				fi;
			od;
		od;
	fi;
	return ish;
fi;
end;

ishomt:=function(t,e,p)
return ishom([[t,specialq(e,p)^0]],e,p);
end;

onetab:=function(la,mu,e,p)
#finds all row-standard tableaux which individually give non-zero homomorphisms (at the moment, doesn't make any attempt to describe linear dependence of these homomorphisms)
#matrices with non-negative integer entries, row sums la and column sums mu.  Assumes they both have positive length
local rs,l,m,ma,i,j,t,cou;
rs:=[];
cou:=0;
if Sum(la)=Sum(mu) then
	l:=Length(la);
	m:=Length(mu);
	if l=1 then
		t:=usualtab([mu]);
		if ishomt(t,e,p) and makesemi(t,e,p)<>[] then Add(rs,t);fi;
	elif m=1 then
		t:=usualtab([List(la,t->[t])]);
		if ishomt(t,e,p) and makesemi(t,e,p)<>[] then Add(rs,t);fi;
	else
		ma:=List([1..l],i->List([1..m],j->0));
		ma[1][1]:=-1;
		i:=1;
		j:=1;
		while i>0 do
			ma[i][j]:=ma[i][j]+1;
			if Sum(ma[i])>la[i] or Sum(List(ma,k->k[j]))>mu[j] then
				ma[i][j]:=0;
				j:=j-1;
				if j=0 then
					i:=i-1;
					if i>0 then
						ma[i][m]:=0;
						j:=m-1;
					fi;
				fi;
			else
				j:=j+1;
				if j=m then
					ma[i][j]:=la[i]-Sum(ma[i]);
					if Sum(TransposedMat(ma)[m])<=mu[m] then
						i:=i+1;
						j:=1;
						if i=l then
							ma[i]:=List([1..m],k->mu[k]-Sum(TransposedMat(ma)[k]));
							cou:=cou+1;
							t:=usualtab(ma);
							if cou mod 100=1 and loud then
								Print("checking ",t," ...\n");
							fi;
							if ishomt(t,e,p) and nonzero(t,e,p) then
								if loud then Print("                ",t,"\n");fi;
								Add(rs,StructuralCopy(t));
							fi;
#							Add(rs,ShallowCopy(List(ma,t->ShallowCopy(t))));
							ma[i]:=List([1..m],k->0);
							i:=l-1;
							ma[i][m]:=0;
							j:=m-1;
						else
							ma[i][j]:=-1;
						fi;
					else
						ma[i][j]:=0;
						j:=j-1;
					fi;
				else
					ma[i][j]:=-1;
				fi;
			fi;
		od;
	fi;
	return rs;
else
	Print("Partitions have different sizes!");
fi;
end;

rsstd:=function(t)
#tells us whether the tableau obtained by arranging each row in decreasing order is reverse semistandard
#assumes the shape is a partition
local tt,i,rs,j;
tt:=StructuralCopy(t);
for i in [1..Length(tt)] do Sort(tt[i],function(x,y) return x>y;end);od;
rs:=true;
for i in [1..Length(tt)-1] do
	for j in [1..Length(tt[i+1])] do
		if tt[i][j]<=tt[i+1][j] then rs:=false;fi;
	od;
od;
return rs;
end;

allsstd:=function(la,mu)
return List(sstd(la,mu),t->usualtab(TransposedMat(t)));
end;

allrstd:=function(la,mu)
return List(sstd(la,Reversed(mu)),t->usualtab(TransposedMat(Reversed(t))));
end;