(*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,n] is the logarithm of the series p(z) up to
degree n, i.e. the infinite product with the output series as
exponents (in the numerator for odd degrees and in the denominator for
even degrees) is equal to p";

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";  
 
(*Begin["`Private`"];*)

Clear["Global`*"];

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

$RecursionLimit=Infinity;

$IterationLimit=Infinity;

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

functions=Names["Global`*"];

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_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]]]]]];

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_, n_, x_List] :=Select[x,deg[#]===n&];  



(*End[];*)

 (*EndPackage[]*)

