(*BeginPackage["liedim`"];

Clear["liedim`*"];*)

sign::usage="sign[x] gives 1 or 0 depending on the odd or even character of
the expression x, which may be a lie-expression or a module-expression.";

weight::usage="weight[x] gives w[..], the weight of
the expression x, which may be a lie-expression or a module-expression.";

w::usage="w[..] is an element of the semigroup of weights, see also plus and deg.";

plus::usage="plus[w[..],w[..]] is the operation in the semigroup of
weights. It may be redefined in input. By default it is defined by
applying Plus to each coordinate (supposing the semigroup is a direct
product of semigroups on which Plus is defined) ";

deg::usage="deg[w[..]] deg is a homomorphism from the semigroup of
weights to the semigroup of positive integers. By default it is the
first coordinate in w[..], but may be redefined in input.";

liedim::usage="liedim[n] gives the dimension of the Lie-algebra in
degree n (if it is computed), liedim[] gives the list of dimensions
computed.";


lie::usage="lie[x,y] is the operator that generates expressions in the free
lie algebra, no rules are defined, see also sq.";

sq::usage="sq[x] is an operator that generates expressions in the free
lie algebra, no rules are defined, see also lie.";

liegen::usage="liegen[n] is the set of generators
having degree n.";

gens::usage="gens[n] is the number of generators of degree n.";

rel::usage="rel[n] is the list of relations of degree n.";

modbas::usage="modbas[x] is the basis element in the module
corresponding to the generator x in the Lie algebra; modbas[n,k] is
the extra basis element number k in degree n; modbas[n], n a positive
integer, is the set of basis elements in degree n.";

msq::usage="msq[modbas[..]] is the squaring operator in the module. If
the characteristic is not 2, it is the same as 1/2
mult[modbas[..],modbas[..]].";

def::usage="def[modbas[..]] def is a homomorphism from the module to
the Lie algebra; def[m] shows how m may be obtained from the
generators of the Lie algebra";

fed::usage="fed[lie[..]] fed is a homomorphism from the Lie algebra to
the module, it is the inverse of def; def[fed[a]] gives a normal form
for the element a in the Lie algebra";


mult::usage="mult[modbas[..],modbas[..]] mult is the induced Lie
operation on the module. It satisfies all axioms for a Lie
multiplication up to the degree computed. See also \"multtable\"";

op::usage="op[lie[..],modbas[..]] this gives the operation of the Lie
algebra on the module. The first argument is an expression generated
by lie, sq and generators, the second argument is any element in the
module (i.e. a linear combination of basis elements modbas[..]).";

newdegree::usage="newdegree[n] this makes all computations in degree
n, if degree up to n-1 is done.";

maxdegree::usage="maxdegree[n] this makes all computations up to (and
including) degree n.";

moredegrees::usage="moredegrees[k,n] this makes all computations from
degree k to degree n, if degree up to k-1 is done.";

input::usage=" you should give definitions of the following functions:
generators, relations and (optionally) genweights (if the generators
are not all of weight 1), gensigns (if the generators are not all
odd), modulus (if not 0), withsquares, cubeszero,
weightpartition. Write ?<function> and ?help for more information.";

help::usage=" write ?input for information of how to give an input. A
list of all functions, which you can use in this program is obtained
by writing \"functions\". You may run the program up to degree n by
executing maxdegree[n]. You can get one more degree by executing
newdegree[n] and the degrees from n to k by executing
moredegrees[n,k]. When the computation is ready, you get the basis
elements in degree k by writing modbas[k]; the binary Lie operation is
obtained by \"mult\" (and the unary square operation in characteristic
2 by \"msq\" ). ";

generators::usage=" generators={a,b,c,...} you give as input a list of
generators named as you prefer (not just integers).";

genweights::usage=" genweights={...} you give as input the weights of
the generators in the same order as \"generators\", the weights may be
positive integers (=degree) or of the form w[n1,n2,...] where n1 is
the degree and n2,... are integers.";

gensigns::usage="gensigns={...} you give as input the signs of the
generators in the same order as \"generators\", the signs are either 0
(for an even generator) or 1 (for an odd generator).";

gendiffl::usage="gendiffl={...} you give as input the differential of
the generators in the same order as \"generators\". The differential
should be of degree -1 and homogeneous of some weight w[-1,..] and
change \"sign\". The weight of the differential is given in input as
difflweight, the default value is w[-1]";

difflweight::usage="difflweight is given in input if not w[-1] which
is the default value; if e.g. difflweight=w[1,2] it means that
weight[diffl[x]]=weight[x]+ w[1,2]";

diffl::usage="diffl[..] this is the extension as a derivation of the
given differential. Its argument could be a lie-expression and then
the output is a linear combination of normalized lie-expressions. It
may also be a modbas-element and then the output is a linear
combination of modbas-elements";

homology::usage="homology[n] gives the dimension of the homology of
the Lie algebra with given differential diffl. The argument may be a
weight, w[..], instead of n.";

homologyweights::usage="homologyweights[w[..]] gives a list of list of
dimensions of homology groups, the first list is homological degree
zero, the next homological degree one etc., the numbers in list number
i gives the dimensions in homological degree i-1 and internal degree
i,i+1,...";

newhomology::usage="newhomology[w[]] gives a list of Lie elements
which are cycles and represent a basis for the homology in the given
weight.";

relations::usage=" relations={...} you give as input the relations,
which are linear combinations of expressions generated by the
generators, the binary operator lie and the unary operator sq (not
just generators however). The relations must be homogeneous with
respect to weight and sign.";

modulus::usage=" this is by default =0, you may change it in input to
any prime number.";

withsquares::usage=" the default value is True. You may change it to
False in input, if you don't want to have special axioms in
characteristic 2 concerning the squaring operator sq.";
		
cubeszero::usage=" the default value is True. You may change it to
False in input, if you don't want to have the special axiom in
characteristic 3 that [x,[x,x]]=0 for odd x.";

weightpartition::usage=" the default value is False. You may change it
to True in input, if you want the program to compute in each weight
separately. This might in some cases be more efficient.";

relmatrix::usage="relmatrix[V,R] is the matrix of coefficients in the
list R of linear expressions in the variables, given by the list V.";

multtable::usage="multtable[n,m] is a 3-dimensional matrix of
coefficients. In place {i,j,k} is the coefficient for modbas[n+m][[k]]
when modbas[n][[i]] is multiplied with modbas[m][[j]]. One may also
use multtable[n,i,m] which gives the 2-dimensional matrix of
coefficients which in place {j,k} has the coefficient for
modbas[n+m][[k]] when modbas[n][[i]] is multiplied with
modbas[m][[j]]. When x is a generator multtable[x,m] gives the matrix
of coefficients for modbas[m+deg[x],j] when modbas[x] is multiplied
with modbas[m][[i]].";

log::usage="log[p,z,y,n] is the \"logarithm\" of the series  p(z,y)=p_0(z)+yp_1(z) up to
degree n, where p_0(z) is the series for the even (super)degrees and
p_1(z) is the  series for the odd degrees. p(z,y) is the infinite
product of  terms (1+z^i)^e_i with the odd output series as
exponents and terms (1-z^i)^(-e_i) with the even output series as
exponents.  E.g., if you want to get the series of the free Lie
algebra on  two even superdegree generators of degree one and three
odd  superdegree generators of degree two, you let
p=1/(1-2z-3yz^2). If you have a series q(z) and want the z-degree determine the
superdegree, you let p=q(yz)";

hcfree::usage="hcfree[p,z,y,n] gives the series for the reduced cyclic
homology of a free algebra on generators having series p(z,y)=p_0(z)+yp_1(z) up to
degree n, where p_0(z) is the series for the even (super)degrees and
p_1(z) is the  series for the odd (super)degrees. The reduced cyclic homology
is T(V)^+/[T(V),T(V)], where V is the vector space with series p(z,y),
this is the same as T(V)^+ modulo the equivalence relation defined on
words as cyclic permutaion. E.g., if you want to get this series for the free 
algebra on  two even superdegree generators of degree one and three
odd  superdegree generators of degree two, you let
p=2z+3yz^2. If you have a series q(z) and want the z-degree determine the
superdegree, you let p=q(yz)";

centre::usage="centre[n] gives generators for the centre of the Lie algebra in degree n";

ideal::usage="ideal[n,x] gives generators in degree n for the ideal generated by a list x consisting of homogeneous elements (linear combinations of modbas[k,.])";

intersection::usage="intersection[x,y] gives the intersection of the subspaces generated by the lists x and y consisting of homogeneous elements of the same degree";

divisor::usage="divisor[a,b,t,s] gives a basis for the elements in \
degree s which multiply the modbas-elements, which should be of degree t, 
in the list a 
to the subspace generated by the modbas-elements in b, which are  of degree s+t";

ann::usage="ann[a,t,s] gives a basis for the elements in degree s \
which multiply the modbas-elements, which are of degree t, in the list a to zero"; 

invimage::usage="invimage[A,B] gives a basis for the vectors x such that Ax is in the rowspace of B, A and B are matrices such that BA is defined";

subalg::usage="subalg[n,x] gives a basis in degree n for the Lie subalgebra generated by the list x of modbas-elements";
 
(*Begin["`Private`"];*)

Clear["Global`*"];

(*Clear["liedim`*"];*)

$RecursionLimit=Infinity;

$IterationLimit=Infinity;

(*functions=Names["liedim`*"];*)

functions={deg,lie,sq,liedim,gens,liegen,gensigns,rel,modbas,generators,relations,mult,msq,def,fed,op,genweights,maxdegree,moredegrees,newdegree,cubeszero,withsquares,weightpartition,help,input,diffl,difflweight,gendiffl,homology,newhomology,homologyweights,modulus,relmatrix,ideal,ann,subalg,log,divisor,invimage,intersection,centre,multtable,sign,weight,plus,w};

inversimage[A_List,f_,y_] := Select[A,f[#]===y&];

relations={};

genweights=1;

gensigns=1;

difflweight=w[-1];

cubeszero=True;

weightpartition=False;

withsquares=True;

modulus=0;




liedim[]={};

liedim[0]=1;

liedimop[1]=0;

modbassq[1]={};

modbas[x_w] := {}/;deg[x]<=0;

modbas[n_Integer] := {}/;n<0;

(*liedim[x_w] :=0;deg[x]<=0;*)

liegen[n_] := liegen[n]=inversimage[generators,deg,n];

liegen[x_w] := liegen[x]=inversimage[generators,weight,x];

gens[n_] := Length[liegen[n]];

rel[k_] := rel[k]=inversimage[relations,deg,k];

rel[k_w] := rel[k]=inversimage[relations,weight,k];


sign[sq[s_]] := 0;

sign[sqmod[s_]] := 0;

sign[msq[s_]] := 0;

sign[t_lie] := Mod[Plus@@Map[sign,List@@t],2];

sign[t_liemod] := Mod[Plus@@Map[sign,List@@t],2];

sign[t_op] := Mod[Plus@@Map[sign,List@@t],2];

sign[t_modbas] := sign[t]=sign[deflazy[t]];

sign[0]=0;

sign[t_Plus] := sign[t[[1]]];

sign[t_Times] := sign[t[[2]]];

sign[x_] := sign[x]=
	gensigns[[Position[generators,x][[1,1]]]];



msq[0] := 0;

msq[t_Plus] := Block[{t1=t[[1]],t2=Drop[t,1]},
	msq[t1]+msq[t2]+imult[t1,t2]];

msq[t_Times] := t[[1]]^2 msq[t[[2]]]//Expand;


weight[t_lie] := plus@@Map[weight,List@@t];

weight[t_liemod] := plus@@Map[weight,List@@t];

weight[sq[t_]] := plus[weight[t], weight[t]];

weight[msq[t_]] := plus[weight[t], weight[t]];

weight[sqmod[t_]] := plus[weight[t], weight[t]];

weight[t_op] := plus@@Map[weight,List@@t];

weight[t_modbas] := weight[t]=weight[deflazy[t]];

weight[t_Plus] := weight[t[[1]]];

weight[t_Times] := weight[t[[2]]];

weight[x_] := weight[x]=
	If[MemberQ[generators,x],genweights[[Position[generators,x][[1,1]]]],w[0]];

w/:plus[x__w] := w@@Plus@@(Apply[List,#]&/@List[x]);

w/:minus[x_w,y_w] := w@@(Apply[List,x]-Apply[List,y]);



(*deg[x_w] := Plus@@x;*)

deg[x_w] := x[[1]];

deg[x_List] := deg/@x;

deg[x_] := deg[weight[x]]/;!Head[x]===weight;


diffl[lie[x_,y_]]:=def[fed[
	lie[diffl[x],y]+(-1)^sign[x] lie[x,diffl[y]]]];

diffl[x_modbas]:=fed[diffl[def[x]]];

diffl[x_List] := diffl/@x;

diffl[0] := 0;

diffl[t_Plus] := Expand[diffl /@ t];

diffl[t_Times] := Expand[t[[1]]*diffl[t[[2]]]];


def[0] := 0;

def[t_Plus] := Map[def,t]//Expand;

def[t_Times] := t[[1]] def[t[[2]]]//Expand;

def[op[x_,y_]] := lie[x,def[y]];

def[msq[x_]] := sq[def[x]];

def[modbas[x_]] := x;


deflazy[op[x_,y_]] := liemod[x,y];

deflazy[msq[x_]] := sqmod[x];

deflazy[modbas[x_]] := x;


imult[0,t_] := 0;

imult[s_,0] := 0;

imult[s_,t_Plus] := Map[imult[s,#]&,t]//Expand;

imult[s_,t_Times] := t[[1]] imult[s,t[[2]]]//Expand;

imult[s_Times,t_] := s[[1]] imult[s[[2]],t]//Expand;

imult[s_Plus,t_] := Map[imult[#,t]&,s]//Expand;

imult[a_List,t_] := imult[#,t]&/@a;

imult[a_,t_List] := imult[a,#]&/@t;

imult[b_,c_] := imult[b,c]=op[deflazy[b],c]; 

premult[b_,c_] := preop[deflazy[b],c]; 

preop[b_sqmod,c_] := premult[b[[1]],imult[b[[1]],c]];

preop[b_liemod,c_] := Block[{sfir=b[[1]],sseco=b[[2]]},
	Expand[op[sfir,imult[sseco,c]]-
	(-1)^(sign[sfir] sign[sseco]) premult[sseco,op[sfir,c]]]];

preop[b_,c_] := op[b,c];



mult[x_,y_] := If[modulus==0,imult[x,y],mod[imult[x,y]]];


fed[lie[x_,b_]] :=  op[x,fed[b]];

fed[sq[x_]] := msq[fed[x]];

fed[t_Times] := t[[1]] fed[t[[2]]]//Expand;

fed[a_Plus] := Map[fed,a]//Expand;

fed[0]=0;

fed[x_] := modbas[x];
	 


op[sq[s_],t_modbas] := 
	Expand[op[s,op[s,t]]];

op[s_lie,t_modbas] := 
	Block[{sfir=s[[1]],sseco=s[[2]]},
	Expand[op[sfir,op[sseco,t]]-
	(-1)^(sign[sfir] sign[sseco]) op[sseco,op[sfir,t]]]];

op[sqmod[s_],t_modbas] := 
	imult[s,imult[s,t]];

op[s_liemod,t_modbas] :=
	Block[{sfir=s[[1]],sseco=s[[2]]},
	Expand[op[sfir,imult[sseco,t]]-
	(-1)^(sign[sfir] sign[sseco]) imult[sseco,op[sfir,t]]]];

op[s_,0] := 0;

op[0,t_] := 0;

op[s_,t_Plus] := Map[op[s,#]&,t]//Expand;

op[s_,t_Times] := t[[1]] op[s,t[[2]]]//Expand;

op[s_Times,t_] := s[[1]] op[s[[2]],t]//Expand;

op[s_Plus,t_] := Map[op[#,t]&,s]//Expand;

op[a_List,t_] := op[#,t]&/@a;

op[a_,t_List] := op[a,#]&/@t;

op[x_,mb0] := fed[x];

op[x_,msq[y_]] := (-1)^(sign[x]) premult[y,op[x,y]];


mod[0] := 0;

mod[t_Plus] := Map[mod,t];

mod[t_Times] := Mod[t[[1]],modulus] t[[2]];

mod[n_Integer] := Mod[n,modulus];

mod[x_List] := Map[mod,x];

mod[x_] := x;

checkrel[{},n_Integer] := True;

checkrel[x_List,n_Integer] := 
	Block[{xfir=x[[1]]},If[checksq[xfir],
	If[!Head[xfir]===Plus,checkrel[Drop[x,1],n+1],
	If[!Length[Union[weight/@List@@xfir]]==1,
	Print["The relation number ",n," 
	is not weight-homogeneous"],
	If[!Length[Union[sign/@List@@xfir]]==1,
	Print["The relation number ",n," 
	is not sign-homogeneous"],
	checkrel[Drop[x,1],n+1]]]]]];

checkrel[x_List] := Block[{badx},
	badx=Select[Complement[Level[x,{-1}],generators],
	Head[#]===Symbol&];If[!badx==={},
	Print["You have used ",badx," in the relations without 
	declaring them as generators"],checkrel[x,1]]];

checksq[{}] := True;

checksq[x_Plus] := 
	Block[{xfir=x[[1]]},
	If[checksq[xfir],checksq[Drop[x,1]]]];

checksq[x_Times] := checksq[x[[2]]];

checksq[sq[x_]] := 
	If[sign[x]==1,True,
	Print[x," is even and cannot be squared"];False];

checksq[x_lie] := True;

	
newdegree[1] := Block[
	{weights1,
	weights1num,weights1union,weightdim},
	If[Head[generators]===Symbol||Head[relations]===Symbol,
	Print["You must give an input, 
	write ?input for help"],
	gennumber=Length[generators];
	If[genweights==1,
	genweights=Table[w[1],{gennumber}]];
	If[!Length[genweights]==gennumber,
	Print["The number of weights should
	be equal to the number of generators"],
	If[Union[Head/@genweights]==={Integer},
	genweights=w/@genweights];
	If[!Select[genweights,!Head[#]===w&]==={},
	Print["All weights should
	be numbers or of the form w[...]"],
	gendeg=deg/@genweights;
	gendegmax=Max[gendeg];
	If[gensigns==0||gensigns==1,
	gensigns=Table[gensigns,{gennumber}]];
	If[!liedim[]==={},
	Print["You have already computed degree 1"],
	If[!Length[gensigns]==gennumber,
	Print["The number of signs should
	be equal to the number of generators"],
	If[!Select[gensigns,!(#==0||#==1)&]==={},
	Print["The signs should be either 0 or 1"],
	If[checkrel[relations],
	reldegmax=Max[Max@@(deg/@relations),0];
	weights1=inversimage[genweights,deg,1];
	liedim[1]=Length[weights1];
	weights1union=Union[weights1];
	weights1num=Length[weights1union];
	weightdim=Count[weights1,#]&/@weights1union;
	(*ReleaseHold/@Thread[Hold[(liedim[#1]=#2)&][
	weights1union,weightdim]];*)
	Evaluate[liedim/@weights1union]=weightdim;
	Evaluate[modbas/@weights1union]=modbas/@#&/@(liegen/@weights1union);
	liedim[x_w]:=0/;deg[x]==1;
	modbas[x_w]:={}/;deg[x]==1;
	Print["weight=",weights1union[[#]],
	": dim=",weightdim[[#]]]&/@Range[weights1num];
	If[modulus>0,badgen=
	Select[generators,mod[deg[#]]==0&]];
	If[modulus==0,msq[x_modbas] := 1/2 imult[x,x];
	(*sq[x_] := 1/2 lie[x,x]*)];
	If[modulus>2,msq[x_modbas] := 
	PowerMod[2,-1,modulus] imult[x,x];
	(*sq[x_] := 
	PowerMod[2,-1,modulus] lie[x,x]*)];
	liedim[]={liedim[1]}];
	If[Head[gendiffl]===List,
	tilld[x_,y_] := diffl[x]=y;
	ReleaseHold/@Thread[Hold[tilld][generators,gendiffl]];
	]]]]]]]];
	

modbas[0]={mb0};
	
modbas[n_Integer] := If[modulus==0,modbasop[n],
	Join[modbasop[n],modbassq[n]]];
		
modbasop[n_Integer] := Join[modbas/@liegen[n],
		modbas[n,#]&/@Range[liedimop[n]]];
	
relcomm[x_Plus] := relcomm/@List@@x;

relcomm[x_Times] := relcomm[x[[2]]];

relcomm[x_lie] :=Expand[op[x[[1]],fed[x[[2]]]]+
		(-1)^
		(sign[x[[1]]] sign[x[[2]]]) op[x[[2]],
		fed[x[[1]]]]];
		 
relcomm[x_sq] := 0;

badrel[x_,n_]:=If[deg[x]>=n,{},
	mod[Expand[op[x,#]+(-1)^(sign[x] sign[#]) premult[#,
	modbas[x]]]]&/@modbas[n-deg[x]]];
	
coefflist[x_,y_List] :=
	Coefficient[x,#]&/@y;
	
relmatrix[{},{}] := {{0}};
	
relmatrix[v_List,{}] := {Table[0,{Length[v]}]};

relmatrix[var_List,rel_List] :=
	coefflist[#,var]&/@rel;
	
skipzeroes[m_List] := m/;Length[m]==1;

skipzeroes[m_List] :=
	If[Union[Last[m]]==={0},
	skipzeroes[Drop[m,-1]],m]/;Length[m]>1;
	
decomp[{}] := {};

decomp[{m1_List}] := 
	If[Union[m1]==={0},{},
	Position[m1,1][[1]]];
	
decomp[m_List] := 
	Join[Position[m[[1]],1][[1]],
	decomp[Drop[m,1]]];
	
newop[0] = 0;

newop[x_] := If[Head[x]===op||Head[x]===msq,x=0,
	If[modulus==0,
	Evaluate[Last[x]]=Expand[-Drop[x,-1]],
	Evaluate[Last[x]]=mod[Expand[-Drop[x,-1]]]]];
	
newop[x_List] := newop/@x;


		
weightorder[A_List,w_List] := 
	inversimage[A,weight,#]&/@w;	
				
newdegree[n_] := Block[
	{rellist,basiselts,relm,basis,newbasiselts,
	newrel,lied,modb,genwei,wbas,
	diffweights,dims,modbasweights,basiseltssq,
	basiseltsall,basisop,basissq,basiseltsop,
	liedop,weightset},
	
	If[Length[liedim[]]<n-1,
		Print["You must compute degree ",n-1," first"],
	If[Length[liedim[]]>n-1,
		Print["You have already computed degree ",n],
	
If[!weightpartition,
		
If[modulus==0,
		
	rellist=Flatten[
		Outer[op,rel[#],modbas[n-#]]&/@
		Range[2,Min[reldegmax,n]]];
			
	rellist=Join[rellist,Flatten[relcomm/@rel[n]]];
	
	basiselts=Flatten[
		Outer[op,liegen[#],
		modbas[n-#]]&/@
		Range[Min[gendegmax,n-1]]];

If[basiselts==={},genwei=inversimage[genweights,deg,n];
		diffweights=Union[genwei];
		dims=Count[genwei,#]&/@diffweights;
		liedimop[n]=0;liedim[n]=Length[genwei];
		liedim[]=Append[liedim[],liedim[n]];
		Evaluate[liedim/@diffweights]=dims;
		Evaluate[modbas/@diffweights]=
			(modbas/@liegen[#])&/@diffweights;
		liedim[x_w]:=0/;deg[x]==n;
		modbas[x_w]:={}/;deg[x]==n;
		Print["weight=",diffweights[[#]],
		": dim=",dims[[#]]]&/@Range[Length[dims]],
	
	relm=skipzeroes[RowReduce[relmatrix[
		basiselts,rellist]]];

	basis=Complement[
		Range[Length[basiselts]],decomp[relm]];
	
	newbasiselts=basiselts[[basis]];

	newrel=Expand[premult[#[[2]],modbas[#[[1]]]]+
		(-1)^(sign[#[[1]]] sign[#[[2]]]) #]&/@
		newbasiselts;
		
	relm=skipzeroes[RowReduce[Join[
		relm,relmatrix[basiselts,newrel]]]];
		
	basis=Complement[basis,decomp[relm]];
	
	newbasiselts=basiselts[[basis]];
		
	liedimop[n]=lied=Length[newbasiselts];
	
	liedim[n]=lied+gens[n];
	
	liedim[]=Append[liedim[],liedim[n]];
	
	modb=modbas[n,#]&/@Range[lied];

	genwei=inversimage[genweights,deg,n];
	
	wbas=Join[weight/@newbasiselts,genwei];
	
	diffweights=Union[wbas];
	
	dims=Count[wbas,#]&/@diffweights;
	
	modbasweights=Join[modbas/@liegen[#],
		inversimage[newbasiselts,weight,#]]&/@diffweights;

	Print["weight=",diffweights[[#]],
		": dim=",dims[[#]]]&/@Range[Length[dims]];
		
	Evaluate[deflazy/@modb]=deflazy/@newbasiselts;
	
	Evaluate[def/@modb]=def/@newbasiselts;
	
	Evaluate[newbasiselts]=modb;
				
	Evaluate[liedim/@diffweights]=dims;

	Evaluate[modbas/@diffweights]=modbasweights;

	liedim[x_w]:=0/;deg[x]==n;
	
	modbas[x_w]:={}/;deg[x]==n;
	
	newop[relm.basiselts];liedim[]],

	
(* else, if weightpartition=False and modulus>0 *)
	

	rellist=mod/@Flatten[
		Outer[op,rel[#],modbas[n-#]]&/@
		Range[2,Min[reldegmax,n]]];
		
	rellist=Join[rellist,
		Flatten[mod[relcomm[#]]&/@rel[n]],
		Flatten[badrel[#,n]&/@badgen]];

	If[modulus==2&&Mod[n,2]==0,
		If[withsquares,
		rellist=Join[rellist,
		mod[premult[#,#]]&/@modbasop[n/2]],
		rellist=Join[rellist,
		mod[premult[#,#]]&/@
		inversimage[modbasop[n/2],sign,0]]]];
 
	If[modulus==3&&cubeszero&&Mod[n,3]==0,
		rellist=Join[rellist,
		mod[premult[imult[#,#],#]]&/@
		inversimage[modbasop[n/3],sign,1]]];
		
	basiselts=
		Flatten[Outer[op,liegen[#],
		modbasop[n-#]]&/@
		Range[Min[gendegmax,n-1]]];
		
	basiseltssq=
		If[modulus==2&&EvenQ[n]&&withsquares,
		msq/@inversimage[modbas[n/2],sign,1],{}];
		
	basiseltsall=Join[basiseltssq,basiselts];

If[basiseltsall==={},genwei=inversimage[genweights,deg,n];
		diffweights=Union[genwei];
		dims=Count[genwei,#]&/@diffweights;
		liedimop[n]=0;liedim[n]=Length[genwei];
		modbassq[n]={};
		liedim[]=Append[liedim[],liedim[n]];
		Evaluate[liedim/@diffweights]=dims;
		Evaluate[modbas/@diffweights]=
			(modbas/@liegen[#])&/@diffweights;
		liedim[x_w]:=0/;deg[x]==n;
		modbas[x_w]:={}/;deg[x]==n;
		Print["weight=",diffweights[[#]],
		": dim=",dims[[#]]]&/@Range[Length[dims]],
			
	relm=skipzeroes[RowReduce[relmatrix[
		basiseltsall,rellist],Modulus->modulus]];

	basis=Complement[
		Range[Length[basiseltsall]],decomp[relm]];
		
	basissq=Intersection[basis,Range[Length[basiseltssq]]];
	
	basisop=Complement[basis,basissq];
	
	newbasiselts=basiseltsall[[basisop]];

	newrel=mod[Expand[premult[#[[2]],modbas[#[[1]]]]+
		(-1)^(sign[#[[1]]] sign[#[[2]]]) #]]&/@
		newbasiselts;
		
	relm=skipzeroes[RowReduce[Join[relm,
		relmatrix[basiseltsall,newrel]],Modulus->modulus]];

	basisop=Complement[basisop,decomp[relm]];
	
	basissq=Complement[basissq,decomp[relm]];
	
	basiseltsop=basiseltsall[[basisop]];
	
	basiseltssq=modbassq[n]=basiseltsall[[basissq]];

	newbasiselts=Join[basiseltsop,basiseltssq];
	
	liedimop[n]=liedop=Length[basiseltsop];
	
	liedim[n]=liedop+Length[basiseltssq]+gens[n];
	
	liedim[]=Append[liedim[],liedim[n]];
	
	modb=modbas[n,#]&/@Range[liedop];

	genwei=inversimage[genweights,deg,n];
	
	wbas=Join[weight/@newbasiselts,genwei];
	
	diffweights=Union[wbas];
	
	dims=Count[wbas,#]&/@diffweights;

	modbasweights=Join[modbas/@liegen[#],
		inversimage[newbasiselts,weight,#]]&/@diffweights;

	Print["weight=",diffweights[[#]],
		": dim=",dims[[#]]]&/@Range[Length[dims]];

	Evaluate[deflazy/@modb]=deflazy/@basiseltsop;
	
	Evaluate[def/@modb]=def/@basiseltsop;
	
	Evaluate[basiseltsop]=modb;
			
	Evaluate[liedim/@diffweights]=dims;
	
	Evaluate[modbas/@diffweights]=modbasweights;

	liedim[x_w]:=0/;deg[x]==n;
	
	modbas[x_w]:={}/;deg[x]==n;
	
	newop[relm.basiseltsall];liedim[]]],

	
(*else, if weightpartition=True*)


If[modulus==0,
	
	rellist=Flatten[
		Outer[op,rel[#],modbas[n-#]]&/@
		Range[2,Min[reldegmax,n]]];
			
	rellist=Join[rellist,Flatten[relcomm/@rel[n]]];
	
	basiselts=Flatten[
		Outer[op,liegen[#],
		modbas[n-#]]&/@
		Range[Min[gendegmax,n-1]]];

	genwei=inversimage[genweights,deg,n];
	
	weightset=Union[Join[weight/@basiselts,genwei]];

	basiselts=weightorder[basiselts,weightset];

	rellist=weightorder[rellist,weightset];

	liedimop[n]=0;

	ReleaseHold/@Thread[Hold[newweight]
		[n,basiselts,rellist,weightset]];

	liedim[n]=liedimop[n]+gens[n]; 

	liedim[x_w] := 0/;deg[x]==n;

	modbas[x_w] := {}/;deg[x]==n;
	
	liedim[]=Append[liedim[],liedim[n]],


(* else, if weightpartition=True and modulus>0 *)


	rellist=mod/@Flatten[
		Outer[op,rel[#],modbas[n-#]]&/@
		Range[2,Min[reldegmax,n]]];
		
	rellist=Join[rellist,
		Flatten[mod[relcomm[#]]&/@rel[n]],
		Flatten[badrel[#,n]&/@badgen]];

	If[modulus==2&&Mod[n,2]==0,
		If[withsquares,
		rellist=Join[rellist,
		mod[premult[#,#]]&/@modbasop[n/2]],
		rellist=Join[rellist,
		mod[premult[#,#]]&/@
		inversimage[modbasop[n/2],sign,0]]]];
 
	If[modulus==3&&cubeszero&&Mod[n,3]==0,
		rellist=Join[rellist,
		mod[premult[imult[#,#],#]]&/@
		inversimage[modbasop[n/3],sign,1]]];

		
	basiselts=
		Flatten[Outer[op,liegen[#],
		modbasop[n-#]]&/@
		Range[Min[gendegmax,n-1]]];
		
	basiseltssq=
		If[modulus==2&&EvenQ[n]&&withsquares,
		msq/@inversimage[modbas[n/2],sign,1],{}];

	genwei=inversimage[genweights,deg,n];
	
	weightset=Union[Join[weight/@basiselts,genwei]];

	weightset=Union[Join[weight/@basiselts,
		weight/@basiseltssq,genwei]];

	basiselts=weightorder[basiselts,weightset];

	basiseltssq=weightorder[basiseltssq,weightset];

	rellist=weightorder[rellist,weightset];

	liedimop[n]=0;

	modbassq[n]={};

	ReleaseHold/@Thread[Hold[newweight]
		[n,basiselts,basiseltssq,rellist,weightset]];

	liedim[n]=liedimop[n]+Length[modbassq[n]]+gens[n]; 
	
	liedim[x_w] := 0/;deg[x]==n;

	modbas[x_w] := {}/;deg[x]==n;
	
	liedim[]=Append[liedim[],liedim[n]]]]]]];


(* if modulus=0 *)
	
newweight[n_,ba_,re_,we_] := Block[
	{relm,basis,newbasiselts,newrel,lied,gendim,modb},

	gendim=Count[genweights,we];

	If[ba==={},
		liedim[we]=gendim;
		modbas[we]=modbas/@liegen[we];
		If[!gendim==0,Print["weight=",we,": dim=",gendim]],
	
	relm=skipzeroes[RowReduce[relmatrix[
		ba,re]]];

	basis=Complement[
		Range[Length[ba]],decomp[relm]];
	
	newbasiselts=ba[[basis]];

	newrel=Expand[premult[#[[2]],modbas[#[[1]]]]+
		(-1)^(sign[#[[1]]] sign[#[[2]]]) #]&/@
		newbasiselts;
		
	relm=skipzeroes[RowReduce[Join[
		relm,relmatrix[ba,newrel]]]];
		
	basis=Complement[basis,decomp[relm]];
	
	newbasiselts=ba[[basis]];

	lied=Length[newbasiselts];

	liedim[we]=lied+gendim;

	modb=modbas[n,#]&/@Range[liedimop[n]+1,liedimop[n]+lied];

	Evaluate[deflazy/@modb]=deflazy/@newbasiselts;
	
	Evaluate[def/@modb]=def/@newbasiselts;
	
	Evaluate[newbasiselts]=modb;

	modbas[we]=Join[modbas/@liegen[we],newbasiselts];

	liedimop[n]=lied+liedimop[n];

	If[!liedim[we]==0,Print["weight=",we,": dim=",liedim[we]]];
	
	newop[relm.ba]]];

	
newweight[n_,ba_,basq_,re_,we_] := Block[
	{basiseltsall,relm,basis,basissq,basisop,newbasiselts,newrel,
	basiseltssq,lied,gendim,modb},

	gendim=Count[genweights,we];
	
	basiseltsall=Join[basq,ba];

	If[basiseltsall==={},
		liedim[we]=gendim;
		modbas[we]=modbas/@liegen[we];
		If[!gendim==0,Print["weight=",we,": dim=",gendim]],

	relm=skipzeroes[RowReduce[relmatrix[
			basiseltsall,re],Modulus->modulus]];

	basis=Complement[
		Range[Length[basiseltsall]],decomp[relm]];
		
	basissq=Intersection[basis,Range[Length[basq]]];
	
	basisop=Complement[basis,basissq];
	
	newbasiselts=basiseltsall[[basisop]];

	newrel=mod[Expand[premult[#[[2]],modbas[#[[1]]]]+
		(-1)^(sign[#[[1]]] sign[#[[2]]]) #]]&/@
		newbasiselts;
		
	relm=skipzeroes[RowReduce[Join[relm,
		relmatrix[basiseltsall,newrel]],Modulus->modulus]];

	basisop=Complement[basisop,decomp[relm]];
	
	basissq=Complement[basissq,decomp[relm]];
	
	newbasiselts=basiseltsall[[basisop]];
	
	basiseltssq=basiseltsall[[basissq]];

	modbassq[n]=Join[modbassq[n],basiseltssq];
	
	lied=Length[newbasiselts];

	liedim[we]=lied+Length[basiseltssq]+gendim;
	
	modb=modbas[n,#]&/@Range[liedimop[n]+1,liedimop[n]+lied];

	Evaluate[deflazy/@modb]=deflazy/@newbasiselts;
	
	Evaluate[def/@modb]=def/@newbasiselts;
	
	Evaluate[newbasiselts]=modb;

	modbas[we]=Join[modbas/@liegen[we],newbasiselts,basiseltssq];

	liedimop[n]=lied+liedimop[n];

	If[!liedim[we]==0,Print["weight=",we,": dim=",liedim[we]]];

	newop[relm.basiseltsall]]];


	
moredegrees[i_,j_] := 
	((newdegree[#];)&/@Range[i,j];liedim[]);

maxdegree[n_] := moredegrees[1,n];

multtable[n_Integer,m_Integer] := 
	relmatrix[modbas[n+m],#]&/@mult[modbas[n],modbas[m]];

multtable[n_Integer,k_Integer,m_Integer] := 
	relmatrix[modbas[n+m],mult[modbas[n][[k]],modbas[m]]];

multtable[x_,m_Integer] := 
	relmatrix[modbas[m+deg[x]],mult[modbas[x],modbas[m]]];
	
liedim[j_,k_]:= 
	Length[skipzeroes[RowReduce[relmatrix[
	modbas[j+k],Join@@mult[modbas[j],modbas[k]]]]]];



(*Format[lie[x_lie,y_lie]]:=
  Block[{prex=ToString[Format[x]],prey=ToString[Format[y]]},
    StringJoin["[",prex,",",prey,"]"]];

Format[lie[x_lie,y_]]:=
  Block[{prex=ToString[Format[x]],prey=ToString[y]},
    StringJoin["[",prex,",",prey,"]"]];

Format[lie[x_,y_lie]]:=
  Block[{prex=ToString[x],prey=ToString[Format[y]]},
    StringJoin["[",prex,",",prey,"]"]];

Format[lie[x_,y_]]:=
  StringJoin["[",ToString[x],",",ToString[y],"]"];*)



cycles[1]=liedim[1];

cycles[n_]:=
  Block[{modut=modbas[n-1],modin=modbas[n]},
    If[modut==={},liedim[n],
      If[modin==={},0,
        NullSpace[Transpose[relmatrix[modut,diffl[modin]]]]//Length]]];

boundaries[n_]:=liedim[n+1]-cycles[n+1];

homology[n_]:= cycles[n]-liedim[n+1]+cycles[n+1];

allhomology[n_]:=homology/@Range[n-1];

cycles[x_w]:=
  Block[{modut=modbas[plus[x,difflweight]],modin=modbas[x]},
    If[modut==={},liedim[x],
      If[modin==={},0,If[modulus==0,
        NullSpace[Transpose[relmatrix[modut,diffl[modin]]]]//Length,
	NullSpace[Transpose[relmatrix[modut,diffl[modin]]],
	Modulus->modulus]//Length
	]]]];

realcycles[x_w]:=
  Block[{modut=modbas[plus[x,difflweight]],modin=modbas[x]},
    If[modut==={},IdentityMatrix[liedim[x]],
      If[modin==={},{0},If[modulus==0,
        NullSpace[Transpose[relmatrix[modut,diffl[modin]]]],
	NullSpace[Transpose[relmatrix[modut,diffl[modin]]],
	Modulus->modulus]
	]]]];

realboundaries[x_w]:=
  Block[{modin=modbas[minus[x,difflweight]],modut=modbas[x]},
    If[modin==={}||modut==={},{},If[modulus==0,
        skipzeroes[RowReduce[relmatrix[modut,diffl[modin]]]],
	skipzeroes[RowReduce[relmatrix[modut,diffl[modin]],Modulus->modulus]]
	]]];

realhomology[x_w]:=Block[{modi=modbas[x],realb=realboundaries[x],conum,
	coba,modout=modbas[plus[x,difflweight]]},
	conum=Complement[Range[Length[modi]],decomp[realb]];
	coba=modi[[conum]];
	If[coba==={},{0},If[modout==={},coba,If[modulus==0,
        NullSpace[Transpose[relmatrix[modout,diffl[coba]]]],
	NullSpace[Transpose[relmatrix[modout,diffl[coba]]],
	Modulus->modulus]]]]];

newhomology[x_w]:=Block[{modi=modbas[x],realb=realboundaries[x],conum,
	coba,modout=modbas[plus[x,difflweight]]},
	conum=Complement[Range[Length[modi]],decomp[realb]];
	coba=modi[[conum]];
	If[coba==={},{0},If[modout==={},def/@coba,If[modulus==0,
        def[Dot[coba, #]] & /@ 
	NullSpace[Transpose[relmatrix[modout,diffl[coba]]]],
	def[Dot[coba, #]] & /@ 
	NullSpace[Transpose[relmatrix[modout,diffl[coba]]],
	Modulus->modulus]]]]];

modbashomology[x_w]:=Block[{modi=modbas[x],realb=realboundaries[x],conum,
	coba,modout=modbas[plus[x,difflweight]]},
	conum=Complement[Range[Length[modi]],decomp[realb]];
	coba=modi[[conum]];
	If[coba==={},{0},If[modout==={},coba,If[modulus==0,
        Dot[coba, #] & /@ 
	NullSpace[Transpose[relmatrix[modout,diffl[coba]]]],
	Dot[coba, #] & /@ 
	NullSpace[Transpose[relmatrix[modout,diffl[coba]]],
	Modulus->modulus]]]]];



boundaries[x_w]:=liedim[minus[x,difflweight]]-cycles[minus[x,difflweight]];

homology[x_w]:= cycles[x]-boundaries[x];

homdeghomology[n_,max_]:= homology[w[#,n]]&/@Range[n+1,max];

homologyweights[max_]:=homdeghomology[#,max]&/@Range[0,max-1];

newcycles[x_w] := def[Dot[modbas[x], #]] & /@ realcycles[x];

newboundaries[x_w] := def[Dot[modbas[x], #]] & /@ realboundaries[x];

logeven[f_, z_, k_] := 
  Sum[MoebiusMu[i]/i*
        Series[Log[(Normal[Series[f /. z -> z^i, {z, 0, k}]])], {z, 0, k}] // 
      Normal, {i, k}];

log[f_, z_,y_, k_] := Block[{pluslog,minuslog},
  	pluslog=Sum[MoebiusMu[i]/i*
        Series[Log[(Normal[Series[f /. {y->(-1)^(i+1),z -> z^i}, {z, 0, k}]])],
	{z, 0, k}] //Normal, {i, k}];
	minuslog=Sum[MoebiusMu[i]/i*
        Series[Log[(Normal[Series[f /. {y->(-1),z -> z^i},
	{z, 0, k}]])],
	{z, 0, k}] //Normal, {i, k}];
	Expand[(pluslog+minuslog)/2]+y (Expand[(pluslog-minuslog)/2])];

hclog[f_, z_,y_, k_] := Block[{pluslog,minuslog},
  	pluslog=Sum[EulerPhi[i]/i*
        Series[Log[(Normal[Series[f /. {y->(-1)^(i+1),z -> z^i}, {z, 0, k}]])],
	{z, 0, k}] //Normal, {i, k}];
	minuslog=Sum[EulerPhi[i]/i*
        Series[Log[(Normal[Series[f /. {y->(-1),z -> z^i},
	{z, 0, k}]])],
	{z, 0, k}] //Normal, {i, k}];
	Expand[(pluslog+minuslog)/2]+y
	(Expand[(pluslog-minuslog)/2])];

hcfree[v_,z_,y_,k_] := hclog[1/(1-v),z,y,k]/.y->1;

liefree[v_,z_,y_,k_] := log[1/(1-v),z,y,k];

logodd[p_,z_,n_]:= Block[{prel=log[p/.z->y z,z,y,n]},prel/.y->1];


(*log[p_,z_,n_]:=loglocal[Normal[Series[p,{z,0,n}]]-1,z,n,0];

loglocal[0,z_,n_,R_]:=R;
loglocal[p_,z_,n_,R_]:=Block[{exp,coef,term},If[!Head[p]===Plus,term=p,
	term=First[p]]; 
	exp=Exponent[term,z];
	coef=Coefficient[term,z^exp];
	If[EvenQ[exp],loglocal[Normal[Series[(1+p) (1-z^exp)^coef,{z,0,n}]]-1,
	z,n,R+coef z^exp],
	loglocal[Normal[Series[(1+p)/(1+z^exp)^coef,{z,0,n}]]-1,
	z,n,R+coef z^exp]]];*)


linext[1,f_] := {f[0,z___] := 0;
    f[t_Plus,z___] := Map[f[#,z]&, t] // Expand;
    f[t_Times,z___] := t[[1]]*f[t[[2]],z] // Expand;}


linext[2, f_] := {f[0, y_, z___] := 0; f[x_, 0, z___] := 0;
    f[s_, t_Plus, z___] := Map[f[s, #, z] &, t] // Expand;
    f[s_, Times[t1_,t2_], z___] := t1 f[s, t2, z] // Expand;
    f[Times[s1_,s2_], t_, z___] := s1 f[s2, t, z] // Expand;
    f[s_, Times[t1_,t2_,t3_], z___] := t1 t2 f[s, t3, z] // Expand;
    f[Times[s1_,s2_,s3_], t_, z___] := s1 s2 f[s2, t, z] // Expand;
    f[s_Plus, t_, z___] := Map[f[#, t,z] &, s] // Expand;}

(*linext[2, f_] := {f[0, y_] := 0; f[x_, 0] := 0;
    f[s_, t_Plus] := Map[f[s, #] &, t] // Expand;
    f[s_, t_Times] := t[[1]] f[s, t[[2]]] // Expand;
    f[s_Times, t_] := s[[1]] f[s[[2]], t] // Expand;
    f[s_Plus, t_] := Map[f[#, t] &, s] // Expand;}*)  


(*linext[2, f_, 1] := {f[0, y_, z_] := 0; f[x_, 0, z_] := 0;
    f[s_, t_Plus, z_] := Map[f[s, #, z] &, t] // Expand;
    f[s_, t_Times, z_] := t[[1]] f[s, t[[2]], z] // Expand;
    f[s_Times, t_, z_] := s[[1]] f[s[[2]], t, z] // Expand;
    f[s_Plus, t_, z_] := Map[f[#, t,z] &, s] // Expand;}*)


 (* the following is multiplication in the exterior algebra (extmult)
(the elements are of the form  e.g. e[1,3,5]) together
with a basis for the exterior algebra (extbasis[n]) and also a basis in
each degree (extbasis[d,n]) *)
extmult[x_e, y_e] :=
	Signature[Join[x, y]]*Union[x, y];
 
linext[2, extmult];
 
extmult[x_e,y_List] := extmult[x, #] & /@ y;
 
extmult[x_List, y_e] := extmult[#, y]& /@ x;
 
extmult[x_List, y_List] := Flatten[extmult[#, y] & /@ x];

extbasis[0, n_] := {e[]};
 
extbasis[1, n_] := Table[e[i], {i, n}];

extbasis[k_, n_] := Union[abs[extmult[extbasis[k - 1, n],
	extbasis[1,n]]]];
 
extbasis[n_] := Flatten[extbasis[#, n] & /@ Range[0, n]];

abs[{}] = {};
 
abs[x_e] := x;
 
abs[x_Times] := x[[2]]; 

abs[x_List] :=
    If[First[x] === 0, abs[Drop[x, 1]], 
      Prepend[abs[Drop[x, 1]], abs[First[x]]]];

genform[d_,m_]:=Plus@@(Times[Random[Integer,100000],#]&/@extbasis[d,m]);

genextforms[d_,n_, t_] := Table[genform[d, n], {t}];

cubeszerotest[n_,t_,p_] :=
	NullSpace[relmatrix[extbasis[3, n],extmult[
	genextforms[2,n,t],extbasis[1,n]]], Modulus -> p];

(* Kac-Moody; Cartanmatrisen ges som ca={{..},{..},...}; operationen av "cf" på den 
positiva sidan ges av cf[i,modbas[..]]. Generatorerna för Lie-algebran skall heta g[1],...*)

cf[i_, modbas[g[j_]]] := If[i === j, -h[i], 0];
cf[i_, x_modbas] := 
    Block[{xx = deflazy[x]}, 
        If[xx[[1, 1]] === i, -cmult[h[i], xx[[2]]] + 
            cmult[modbas[g[i]], cf[i, xx[[2]]]], 
          cmult[modbas[xx[[1]]], cf[i, xx[[2]]]]]] /; ! x[[1, 0]] === Symbol;

cmult[h[i_], x_modbas] := Block[{ww = weight[x], le}, le = Length[ww];
      Dot[List @@ Take[ww, -le + 1], ca[[i]]] x // Expand];
cmult[x_modbas, h[i_]] := -cmult[h[i], x];
linext[2, cmult];
cmult[x_modbas, y_modbas] := mult[x, y];

(*centre[n_] := 
  modbas[n].# & /@ (Join @@ ((relmatrix[modbas[n + 1], op[#, modbas[n]]] // 
                    Transpose) & /@ generators) // NullSpace);*)

(*intersection[x,y] is just linear algebra; x and y are lists of
vectors in the same degree n expressed as linear combinations of the
basis elements in modbas[n]; the result is a basis for the
intersection of the subspaces generated by x and y*)

intersection[x_List, y_List] := Block[{nn},If[x==={}||y==={},{},nn=deg[x[[1]]];
   modbas[nn].#&/@ NullSpace[
      Join[NullSpace[relmatrix[modbas[nn], x]],
        NullSpace[relmatrix[modbas[nn], y]]]]]];

intersection[x_List] := intersection[Join[{intersection[x[[1]],x[[2]]]},Drop[x,2]]]/;Length[x]>2;

intersection[{x_List,y_List}] := intersection[x,y];

intersection[{x_List}] := x;

ideal[n_,x_List] := ideal[n,Min@@deg/@x,x]; 


ideal[n_, p_, x_List] := Block[{xp=Select[x,deg[#]===p&],xplus},xplus=Complement[x,xp];
    ideal[n, p + 1,Join[xplus,
        modbas[p + 1].# & /@
          skipzeroes[
            RowReduce[relmatrix[modbas[p + 1], Join @@ op[generators,xp]]]]]]] /;p < n;
ideal[n_,p_,x_List]:= {0}/; p>n;

ideal[n_, n_, x_List] := modbas[n].# & /@
          skipzeroes[
            RowReduce[relmatrix[modbas[n],Select[x,deg[#]===n&]]]];

(* centre[n,k] consist of those elements in degree 
n-k+1 such that multiplication by any generator 
(which should be of degree 1) give elements in 
centre[n,k-1] and centre[n,0]=0. Moreover
 "invimage[A,B]" is defined for an 
nxk-matrix A and a pxn-matrix B and gives a 
basis for the set of elements x such that 
Ax belongs to the rowspace of B. Also nullspace 
is a variant of NullSpace where if the output 
is {}, it is changed to a zero matrix of the 
right dimension. centre[n,k] gives vectors as 
answer, while "centremodbas[n,k]" gives the 
answer as linear combination of modbas-elements. *)

centre[n_, k_] := 
   ( centre[n, k] = 
      Block[{ce = nullspace[centre[n, k - 1]]}, 
          nullspace[
            Join @@ Map[
                ce.# &, (Transpose[
                        relmatrix[modbas[n - k + 2], 
                          op[#1, modbas[n - k + 1]]]] &) /@ generators]]]) /; 
        k > 0;

centre[n_, 0] := {Table[0, {liedim[n + 1]}]};
centre[n_] := centre[n, 1];
nullspace[A_] := 
    Block[{nu = NullSpace[A,Modulus->modulus]}, 
      If[nu === {}, {Table[0, {Length[A[[1]]]}]}, nu]];
invimage[A_, B_] := nullspace[nullspace[B].A];
centremodbas[n_, k_] := (modbas[n + 1 - k].#1 &) /@ centre[n, k];
centremodbas[n_] := (modbas[n].#1 &) /@ centre[n]; 

 (* The program "divisor" below computes a basis for the elements in \
degree s which multiply the modbas-elements in the list a of degree t \
to the subspace generated by the modbas-elements in b of degree s+t. \
The program "ann" below computes a basis for the elements in degree s \
which multiply the modbas-elements in the list a of degree t to zero. \
*)

divisor[a_List, b_List, t_, s_] := 
  Module[{bmat = relmatrix[modbas[s + t], b]}, 
   modbas[s].# & /@ (nullspace[
      Join @@ (nullspace[
            bmat].# & /@ (Transpose[
             relmatrix[modbas[s + t], mult[modbas[s], #]]] & /@ 
           a))])];
ann[a_List, t_, s_] := divisor[a, {}, t, s]; 


skipz[x_List] := If[Last[x] === 0, skipz[Drop[x, -1]], x];
skipz[{}] := {};

subalglist[n_,x_List]:= Select[x,deg[#]===n&];

subalg[n_,x_List] := (subalg[n,x] = Join[subalglist[n,x],modbas[n].# & /@
          skipzeroes[
            RowReduce[relmatrix[modbas[n],Join@@Join@@(mult[subalglist[#,x],subalg[n-#,x]]&/@Range[1,n-1])]]]])/; n>1; 

subalg[1, x_List] := modbas[1].# & /@
          skipzeroes[
            RowReduce[relmatrix[modbas[1],subalglist[1,x]]]];


(*End[];*) 

 (*EndPackage[]*)