(* GAMMA: A Mathematica package for performing \[CapitalGamma]-matrix algebra and Fierz transformations in arbitrary dimensions    *)
(* Version 2.0\[Beta]                                                       *)
(* Author: Ulf Gran                                                    *)
(*         ulf.gran@chalmers.se                                        *)
(*         Copyright 2001                                              *)
(*         Check http://www.gran.name for updates                      *)

(* Redefined built-in functions:                                      *)
(*                                                                    *)
(* - NonCommutativeMultiply                                           *)
(* - Power                                                            *)
(* - Subscript                                                        *)
(* - Times                                                            *)
(* - the symbol used as TensorSpinorIndicator                         *)


(*BeginPackage["Gamma`",{"Combinatorica`"}];*)
BeginPackage["GAMMA`"];
Print["GAMMA Version 2.0\[Beta]"];
Print["9 December 2014"];
Print["\[Copyright] Ulf Gran, 2001"];


Clear["GAMMA`*"];Clear["GAMMA`Private`*"];
Dim=11;
SpinorDim=32;
ChiralSpinors=False;
EpsilonConvention=-1;
SpaceTimeSignature=-1;
CustomFontSize=16;
TensorIndicator=Identity;
TensorSpinorIndicator=OverTilde;
HookTensorIndicator=OverHat;
HookTensorSpinorIndicator=OverTilde;
$RecursionLimit=\[Infinity];
$IterationLimit=\[Infinity];
SetDim::"usage"="SetDim[dim] sets the space-time dimension to dim.";
SetSpinorDim::"usage"="SetSpinorDim[dim] sets the spinor dimension to dim.";
SetChiralSpinors::"usage"="SetChiralSpinors[True].";
SetTensorInd::"usage"="SetTensorInd[func] sets the Mathematica function to be applied to the symbol designating a tensor to func.";
SetTensorSpinorInd::"usage"="SetTensorSpinorInd[func] sets the Mathematica function to be applied to the symbol designating a tensor-spinor, in order to discern it from a tensor, to func.";
SetHookTensorInd::"usage"="SetHookTensorInd[func] sets the Mathematica function to be applied to the symbol designating a hook tensor to func.";
SetHookTensorSpinorInd::"usage"="SetHookTensorSpinorInd[func] sets the Mathematica function to be applied to the symbol designating a hook tensor-spinor, in order to discern it from a hook tensor, to func.";
SetFontSize::"usage"="SetFontSize[size] sets the font size to size, the default value is 16.";
Delta::"usage"="Delta[list1,list2] represents a Kronecker delta with the two groups of indices given by list1 and list2.";
DualiseGammas::"usage"=".";
Flip::"usage"="Transposes expressions where two spinor are contracted with a \[CapitalGamma]-matrix in between in order to be able to use GammaExtract to remove \[CapitalGamma]-traces.";
GammaProd::"usage"="GammaProd[list1,...,listn] represents the product of n \[CapitalGamma]-matrices with indices given by list1 to listn.";
GammaContract::"usage"="GammaContract[expr] simplifies the expression expr containing contracted indices on adjacent \[CapitalGamma]-matrices.";
GC::"usage"="Shorthand for GammaContract.";
GammaLongContract::"usage"="GammaContract[expr] simplifies the expression expr containing contracted indices on non-adjacent \[CapitalGamma]-matrices.";
GLC::"usage"="Shorthand for GammaLongContract.";
GammaExpand::"usage"="GammaExpand[expr] expands the product of the first two \[CapitalGamma]-matrices in expr.";
GE::"usage"="Shorthand for GammaExpand.";
GammaTrace::"usage"="GammaTrace[expr,list] computes the \[CapitalGamma]-trace of an expression containing \[CapitalGamma]-matrices, enforcing antisymmetry in the indices contained in list.";
GT::"usage"="Shorthand for GammaTrace.";
Tensor::"usage"="Tensor[sym,list] represents a tensor, denoted by the symbol sym, with indices given by list.";
HookTensor::"usage"="Hook[sym,list1,list2] represents a hook tensor, denoted by the symbol sym, with tensor indices given by list1 and hook index given by list2.";
TensorSpinor::"usage"="TensorSpinor[sym,list] represents a tensor-spinor, denoted by the symbol sym, with (tensor) indices given by list.";
HookTensorSpinor::"usage"="HookTensorSpinor[sym,list1,list2] represents a hook tensor-spinor, denoted by the symbol sym, with tensor indices given by list1 and hook index given by list2.";
GammaExtract::"usage"="GammaExtract[expr,list] extracts \[CapitalGamma]-matrices contracting indices on a tensor-spinor, assuming antisymmetry in the indices contained in list, in order to be able to invoke the \[CapitalGamma]-tracelessness of the tensor-spinor.";
GET::"usage"="Shorthand for GammaExtract.";
GammaExtractOld::"usage"="";
RemoveGammaTrace::"usage"="RemoveGammaTrace[expr] removes the \[CapitalGamma]-trace of expr assuming that expr is a tensor-spinor.";
RGT::"usage"="Shorthand for RemoveGammaTrace.";
RemoveHookGammaTrace::"usage"="RemoveHookGammaTrace[expr,list1,list2] removes the \[CapitalGamma]-trace of expr assuming that expr lies in the hook representation having list1 as tensor indices and list2 as hook index.";
RHGT::"usage"="Shorthand for RemoveHookGammaTrace.";ACanonicalOrder::"usage"="ACanonicalOrder[expr,list] assumes that expr is antisymmetric in the indices given by list and puts them in canonical order, thus enabling additional simplification.";
AC::"usage"="Shorthand for ACanonicalOrder.";
SCanonicalOrder::"usage"="SCanonicalOrder[expr,list] assumes that expr is symmetric in the indices given by list and puts them in canonical order, thus enabling additional simplification.";
SC::"usage"="Shorthand for SCanonicalOrder.";
TensorCanonicalOrder::"usage"="";
TC::"usage"="Shorthand for TensorCanonicalOrder.";
TensorProject::"usage"="";
TP::"usage"="Shorthand for TensorProject.";
HookCanonicalOrder::"usage"="";
HC::"usage"="Shorthand for HookCanonicalOrder.";
RemoveTrace::"usage"="RemoveTrace[expr,list1,list2] removes the trace of expr w.r.t. the indices in list1 and list2.";
RenameDummy::"usage"="RenameDummy[expr] renames dummy (contracted) indices in expr, thus enabling additional simplification.";
RD::"usage"="Shorthand for RenameDummy.";
RDR::"usage"="Shorthand for RenameDummyRight.";
RDL::"usage"="Shorthand for RenameDummyRight.";
RenameDummyLeft::"usage"="RenameDummyLeft[expr] renames dummy (contracted) indices in expr, starting from the left,thus enabling additional simplification.";RenameDummyRight::"usage"="RenameDummyRight[expr] renames dummy (contracted) indices in expr, starting from the right,thus enabling additional simplification.";
ASym::"usage"="ASym[expr, list] explicitly antisymmetrizes expr in the indices contained in list.";
Sym::"usage"="Sym[expr, list] explicitly symmetrizes expr in the indices contained in list.";
Fierz::"usage"="Fierz[basislist,structures,asymlist,sym,function] performs the\[IndentingNewLine]Fierz transformation of the terms specified in structures using the terms\[IndentingNewLine]in the list basislist as basis elements. Antisymmetry is assumed in the\[IndentingNewLine]indices contained in the list asymlist and sym (A or S) gives the symmetry of the spinor indices to be fierzed.Finally, function is applied to generate explicit \[CapitalGamma]-matrix expressions given a\[IndentingNewLine]term in structures.";
FierzTerm::"usage"="FierzTerm[termlist,asymlist,sym,function] as Fierz[] but\[IndentingNewLine]does not perform the complete Fierz transformation. Instead it produces a list\[IndentingNewLine]of expressions obtained by contracting two of the spinor indices, as specified\[IndentingNewLine]by ind, with the basis of \[CapitalGamma]-matrices specified by sym.";
FierzSolve::"usage"="FierzSolve[eqlist] solves the list of equations generated by Fierz[]";
MakeTensor::"usage"="MakeTensor[term,basiselement] makes an explicit \[CapitalGamma]-matrix\[IndentingNewLine]expression from term contracting the \[CapitalGamma]-matrix basiselement with\[IndentingNewLine]two of the spinor indices,as specified by ind.";
MakeHook::"usage"="MakeHook[term,basiselement] makes an explicit \[CapitalGamma]-matrix\[IndentingNewLine]expression from term contracting the \[CapitalGamma]-matrix basiselement with\[IndentingNewLine]two of the spinor indices, as specified by ind.";
LinearPart::"usage"="LinearPart[expr] picks out the part of expr that is linear, or less, in the fields";
Multiply::"usage"="Multiply[term1,..., termn] multiplies the n terms taking into account that some dummy indices may be the same in some terms.";
HookSimplify::"usage"="";
HS::"usage"="Shorthand for HookSimplify";
HookSimplifyFirst::"usage"="";
HookSimplifyLast::"usage"="";
Pick::"usage"="Pick[expr,obj,length] picks out terms containing obj of expr whrere in the cases when obj has length anti-symmetric indices.";
HookProject::"usage"="HookProject[expr,list1,list2] projects expr onto the hook with list1 as tensor indices and hook index given by list2.";
HP::"usage"="Shorthand for HookProject";
Epsilon::"usage"="";
EpsilonConvention::"usage"="";
SetEpsilonConvention::"usage"="";
SetSignature::"usage"="";
PhiSimplify::"usage"="";

CStar::"usage"="";
ComplexConjugate::"usage"="";
CC::"usage"="";

Begin["`Private`"];

SetAttributes[Subscript,HoldRest];
SetAttributes[Power,HoldRest];
NCList={GammaProd,TensorSpinor,HookTensorSpinor,Global`\[CapitalGamma],Global`\[CapitalGamma]11,CStar,CS,cs};
ComplexFluxes={Global`H,Global`L,Global`\[CapitalPhi], Global`\[Xi],Global`G, Global`P,Global`\[Alpha],Global`\[Sigma],Global`\[Tau]};
FieldList={Tensor,TensorSpinor,HookTensor,HookTensorSpinor};
SetDim[x_]:=Module[{},Dim=x;];
SetSpinorDim[x_]:=Module[{},SpinorDim=x;];
SetFontSize[x_]:=Module[{},CustomFontSize=x;];
SetChiralSpinors[x_]:=Module[{},ChiralSpinors=x;];
SetEpsilonConvention[x_]:=Module[{},EpsilonConvention=x;];
SetSignature[x_]:=Module[{},SpaceTimeSignature=x;];

(* Add support for the Spin(7) form \[Phi] *)

(* Add support for the Subscript[G, 2] form \[CurlyPhi] *)

(* Three types of terms: 1) Two \[CurlyPhi]s contracted 
						2) One \[CurlyPhi] and one \[FivePointedStar]\[CurlyPhi] contracted 
							3 Two \[FivePointedStar]\[CurlyPhi] contracted *)

(* Unresolved problem: two \[CurlyPhi]s can combine to a power *)

PhiSimplify[x_]:=Distribute[PSt[FixedPoint[Expand,x]],Plus,PSt,Plus,PS];

PS[v_. Tensor[Global`\[CurlyPhi],x_List]Tensor[Global`\[CurlyPhi],y_List]]:=Module[{dummy,xcompl,ycompl},
	dummy=Last[inter];
	xcompl=Complement[x,{dummy}];
	ycompl=Complement[y,{dummy}];
	PhiSimplify[v Signature[x]Signature[Join[xcompl,{dummy}]]Signature[y]Signature[Join[ycompl,{dummy}]](2 Delta[xcompl,ycompl]-Tensor[Global`\[FivePointedStar]\[CurlyPhi],Join[xcompl,ycompl]])]]/;((inter=Intersection[x,y])=!={});

PS[v_. Tensor[Global`\[CurlyPhi],x_List]Tensor[Global`\[FivePointedStar]\[CurlyPhi],y_List]]:=Module[{dummy,xcompl,ycompl,x1,x2,y1,y2,y3,xlist,ylist},xlist={x1,x2};ylist={y1,y2,y3};dummy=Last[inter];xcompl=Complement[x,{dummy}];ycompl=Complement[y,{dummy}];PhiSimplify[v Signature[x]Signature[Join[xcompl,{dummy}]]Signature[y]Signature[Join[ycompl,{dummy}]](ASym[ASym[6Delta[{x1},{y1}]Tensor[Global`\[CurlyPhi],{x2,y2,y3}],xlist],ylist])/.Thread[Rule[Join[xlist,ylist],Join[xcompl,ycompl]]]]]/;((inter=Intersection[x,y])=!={});

PS[v_. Tensor[Global`\[FivePointedStar]\[CurlyPhi],x_List]Tensor[Global`\[FivePointedStar]\[CurlyPhi],y_List]]:=Module[{dummy,xcompl,ycompl,x1,x2,x3,y1,y2,y3,xlist,ylist},xlist={x1,x2,x3};ylist={y1,y2,y3};dummy=Last[inter];xcompl=Complement[x,{dummy}];ycompl=Complement[y,{dummy}];PhiSimplify[v Signature[x]Signature[Join[xcompl,{dummy}]]Signature[y]Signature[Join[ycompl,{dummy}]](ASym[ASym[3!Delta[xlist,ylist]-9Delta[{x1},{y1}]Tensor[Global`\[FivePointedStar]\[CurlyPhi],{x2,x3,y2,y3}]-Tensor[Global`\[CurlyPhi],xlist]Tensor[Global`\[CurlyPhi],ylist],xlist],ylist])/.Thread[Rule[Join[xlist,ylist],Join[xcompl,ycompl]]]]]/;((inter=Intersection[x,y])=!={});

PS[v_. Tensor[Global`\[Phi],x_List]Tensor[Global`\[Phi],y_List]]:=Module[{dummy,xcompl,ycompl,xtcompl={x1,x2,x3},ytcompl={y1,y2,y3}},
	dummy=Last[inter];
	xcompl=Complement[x,{dummy}];
	ycompl=Complement[y,{dummy}];
	PhiSimplify[v Signature[x]Signature[Join[xcompl,{dummy}]]Signature[y]Signature[Join[ycompl,{dummy}]]
		(6 Delta[xcompl,ycompl]-9 ASym[ASym[Delta[{First[xtcompl]},{First[ytcompl]}]Tensor[Global`\[Phi],Join[Rest[xtcompl],Rest[ytcompl]]],xtcompl],ytcompl])/.MapThread[Rule, {xtcompl, xcompl}]/.MapThread[Rule, {ytcompl, ycompl}]]]/;((inter=Intersection[x,y])=!={});

PS[x_]:=x;


(* Delta *)
Attributes[Delta]={Orderless};

Delta/:Delta[x_List,y_List]^2:=Delta[x,x];
Delta/:(Delta[{a_,b_},{c_,d_}]-Delta[{a_,c_},{b_,d_}]+Delta[{a_,d_},{b_,c_}])=0;
(* Lists *)
Delta/:Delta[x_List,y_List]z_[w___List,v_List,s___List]u_.:=Module[{same=ReplaceList[{x,v},{{___,t_,___},{___,t_,___}}:>t],poslist1, poslist2,xxleft,lensame,yycomb,lenyycomb,actyycomb,actyycompo,i}, xxleft=Delete[x,poslist1=PositionList[x,same]];
poslist2=PositionList[v,same];
yycomb=Subsets[y,{lensame=Length[same]}];Power[-1,Evaluate[((Flatten[poslist1]/.List->Plus)-lensame (lensame+1)/2)]]/(lenyycomb=Length[yycomb])Sum[Signature[Join[actyycomb=yycomb[[i]],actyycompo=Complement[y,actyycomb]]]Delta[xxleft, actyycompo]z[w,ReplacePart[v,actyycomb,poslist2,(Range[Length[poslist2]]/.a_Integer:>List[a])] ,s]u,{i,lenyycomb}]]/;(Intersection[x,v]=!={});
(* Symbol + Lists *)
Delta/:Delta[x_List,y_List]z_[sym_,w___List,v_List,s___List]u_.:=Module[{same=ReplaceList[{x,v},{{___,t_,___},{___,t_,___}}:>t],poslist1, poslist2,xxleft,lensame,yycomb,lenyycomb,actyycomb,actyycompo,i}, xxleft=Delete[x,poslist1=PositionList[x,same]];
poslist2=PositionList[v,same];
yycomb=Subsets[y,{lensame=Length[same]}];Power[-1,Evaluate[((Flatten[poslist1]/.List->Plus)-lensame (lensame+1)/2)]]/(lenyycomb=Length[yycomb])Sum[Signature[Join[actyycomb=yycomb[[i]],actyycompo=Complement[y,actyycomb]]]Delta[xxleft, actyycompo]z[sym,w,ReplacePart[v,actyycomb,poslist2,(Range[Length[poslist2]]/.a_Integer:>List[a])] ,s]u,{i,lenyycomb}]]/;(Intersection[x,v]=!={});

Delta/:Delta[x_List,y_List]Unevaluated[NonCommutativeMultiply[w___,(d_)[b___,v_List,c___],s___]]u_.:=Module[{same=ReplaceList[{x,v},{{___,t_,___},{___,t_,___}}:>t],poslist1, poslist2,xxleft,lensame,yycomb,lenyycomb,actyycomb,actyycompo,i},
xxleft=Delete[x,poslist1=PositionList[x,same]];
poslist2=PositionList[v,same];
yycomb=Subsets[y,{lensame=Length[same]}];Power[-1,Evaluate[((Flatten[poslist1]/.List->Plus)-lensame (lensame+1)/2)]]/(lenyycomb=Length[yycomb])Sum[Signature[Join[actyycomb=yycomb[[i]],actyycompo=Complement[y,actyycomb]]]Delta[xxleft, actyycompo]NonCommutativeMultiply[w,d[b,ReplacePart[v,actyycomb,poslist2,(Range[Length[poslist2]]/.a_Integer:>List[a])],c] ,s]u,{i,lenyycomb}]]/;(Intersection[x,v]=!={});
Delta[{},{}]:=1;
Delta[x_List,y_List]:=0/;(Length[x]>Dim);
Delta[{___,x_,x_,___},_List]:=0;Delta[x_List,y_List]:=Signature[x]Signature[y]Delta[Sort[x],Sort[y]]/;!(OrderedQ[x]&&OrderedQ[y]);
Delta[x_List,y_List]:=Module[{inter,xxleft,yyleft,lensame,n=Length[x]},Delta[xxleft=Complement[x,inter],yyleft=Complement[y,inter]] Signature[Join[inter,xxleft]]Signature[Join[inter,yyleft]] Binomial[Dim-n+(lensame=Length[inter]),lensame]/Binomial[n,lensame]/;((inter=Intersection[x,y])=!={})];

Format[Delta[x_List,y_List],{StandardForm,TeXForm}]:=DisplayForm[SubsuperscriptBox[Global`\[Delta],RowBox[y],RowBox[x]]];


Unprotect[Power];
Power[Subscript[Global`\[Delta],x_Symbol],y_Symbol]:=Delta[{y},{x}];
Power[Subscript[Global`\[Delta],x_],y_]:=Delta[Unevaluated[y]/.Times->List,Unevaluated[x]/.Times->List];
Protect[Power];


(* GammaProd, GammaContract, GammaExpand & GammaTrace *)

GammaProd[]:=1;
GammaProd[___,{___,x_,x_,___},___]=0;
(* sometimes this is better left out *)
GammaProd[___,x_List,___]:=0/;(Length[x]>Dim);
GammaProd[x___,{},y___]:=GammaProd[x,y];
GammaProd[y___,x_List,z___]:=Signature[x]GammaProd[y,Sort[x],z]/;\[Not]OrderedQ[x];

(* DualiseGammas do not seem correct! Wrong sign factor! *)

(* DualiseGammas[expr_]:=Module[{u0,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16,u17,u18,u19,tempind},
	tempind={u0,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16,u17,u18,u19};
	expr/.v_. GammaProd[x_]:>v (-1)^((Dim-Length[x])(Dim-Length[x]-1)/2)/(Dim-Length[x])! EpsilonConvention 
	Epsilon[Join[x,Take[tempind,Dim-Length[x]]]]GammaProd[Take[tempind,Dim-Length[x]]]/;(Length[x]>Dim/2)] 
	
	*)

Unprotect[NonCommutativeMultiply];

NonCommutativeMultiply/:(_)[___,x_List,___]Unevaluated[NonCommutativeMultiply[___,HookTensorSpinor[_Symbol, z_List, a_List] ,___]]:=0/;Intersection[x,Join[z,a]]===Sort[Join[z,a]];

NonCommutativeMultiply[1,1]:=1;
(*NonCommutativeMultiply[x___,y_+z_,w___]:=Distribute[NonCommutativeMultiply[x,y+z,w]];*)

NonCommutativeMultiply[x___,y_+z_,w___]:=NonCommutativeMultiply[x,y,w]+NonCommutativeMultiply[x,z,w];

NonCommutativeMultiply[x_,1]:=x;
NonCommutativeMultiply[1,x_]:=x;
NonCommutativeMultiply[x__,1,y__]:=NonCommutativeMultiply[x,y];
NonCommutativeMultiply[1,x_,y__]:=NonCommutativeMultiply[x,y];
NonCommutativeMultiply[x__,y_,1]:=NonCommutativeMultiply[x,y];

NonCommutativeMultiply[x__,GammaProd[y__List],GammaProd[z__List],w___]:=NonCommutativeMultiply[x,GammaProd[y,z],w];
NonCommutativeMultiply[x___,GammaProd[y__List],GammaProd[z__List],w__]:=NonCommutativeMultiply[x,GammaProd[y,z],w];
NonCommutativeMultiply[GammaProd[y__List],GammaProd[z__List]]:=GammaProd[y,z];

(*NonCommutativeMultiply[x___,(v_Integer|v_Rational|v_Symbol|v_Complex)z_,y___]:=v NonCommutativeMultiply[x,z,y];*)

NonCommutativeMultiply[x_ /;(Map[FreeQ[x,#]&,NCList]/.List:>And),y_]:=x y;
NonCommutativeMultiply[x_,y_/;(Map[FreeQ[y,#]&,NCList]/.List:>And)]:=x y;
NonCommutativeMultiply[x___,z_. y_/;(Map[FreeQ[y,#]&,NCList]/.List:>And)&&y=!=1,w___]:=y NonCommutativeMultiply[x,z,w];

NonCommutativeMultiply[GammaProd[___,y_List],HookTensorSpinor[_Symbol,z_List,w_List]]:=0/;Intersection[y,Join[z,w]]===Sort[Join[z,w]];
NonCommutativeMultiply[x_,TensorSpinorIndicator[sym_Symbol/;(sym=!=\[Delta]&&sym=!=Global`\[CapitalGamma]&&FreeQ[sym,RowBox])]]:=NonCommutativeMultiply[x,TensorSpinor[sym,{}]];
NonCommutativeMultiply[___,HookTensorSpinor[_Symbol, z_List, a_List],___,(_)[___,x_List,___],___]:=0/;Intersection[x,Join[z,a]]===Sort[Join[z,a]];NonCommutativeMultiply[___,(_)[___,x_List,___],___,HookTensorSpinor[_Symbol, z_List, a_List],___]:=0/;Intersection[x,Join[z,a]]===Sort[Join[z,a]];
Protect[NonCommutativeMultiply];

(* GammaCortract takes a list of lists as the second argument *)
GammaContract[x_,l_List:{{}}]:=Distribute[Gct[FixedPoint[Expand,x],l],Plus,Gct,Plus,Gc];
GC[x_,l_List:{{}}]:=Distribute[Gct[FixedPoint[Expand,x],l],Plus,Gct,Plus,Gc];

Gc[x_]:=Gc[x,{{}}];
Gc[Unevaluated[ NonCommutativeMultiply[x__]v_],l_List]:=v Gc[NonCommutativeMultiply[x],Join[l,Union[DeleteCases[AllIndexLists[v ],{_}]]]];
Gc[GammaProd[x__]v_.,_List]:=GammaProd[x]v/;(ReplaceList[{x},{___List,{___, y_,___},{___,y_,___},___List}:>y]==={});

Gc[GammaProd[z___,x_List,y_List,w___]v_.,l_List]:=Module[{inter,xxleft,yyleft,k,n,r,ltotal=Join[l,{z},{w},Union[DeleteCases[AllIndexLists[v ],{_}]]],xxtotal,yytotal,xxcompo,yycompo},(xxleft=Complement[x,inter];
yyleft=Complement[y,inter];
n=Length[xxleft]; r=Length[yyleft];k=Length[inter];
GammaContract[Signature[Join[inter,xxleft]]Signature[Join[inter,yyleft]]Sum[Module[{xxcomb=Subsets[xxleft,{p}],yycomb=Subsets[yyleft,{p}]},xxcompo=Map[Complement[xxleft,#]&,xxcomb];yycompo=Map[Complement[yyleft,#]&,yycomb];xxtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[xxcomb[[i]],xxcompo[[i]]]],{i,Length[xxcomb]}]];yytotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[yycomb[[i]],yycompo[[i]]]],{i,Length[yycomb]}]];
Power[-1,Evaluate[((2n+k-p-1)(p+k)/2)]]p!(Dim+p-n-r)!/(Dim+p-n-r-k)! v Apply[Plus,Flatten[Outer[Count[xxtotal,#1]Count[yytotal,#2]Delta[Take[#1,p],Take[#2,p]]GammaProd[z,Join[Drop[#1,p],Drop[#2,p]],w]Signature[#1]Signature[#2]&,Union[xxtotal],Union[yytotal],1]]]],{p,Max[0,n+r+k-Dim],Min[n,r]}],l])/;(inter=Intersection[x,y])!={}];

Gc[Unevaluated[NonCommutativeMultiply[x___,GammaProd[z__] ,y___]],l_List]:=NonCommutativeMultiply[x,Gc[GammaProd[z],Join[l,Union[DeleteCases[AllIndexLists[x y],{_}]]]],y];

Gc[x_/;FreeQ[x,GammaProd],_List]:=x;


(* GammaLongContract takes a list of lists as the second argument *)
GammaLongContract[x_,l_List:{{}}]:=Distribute[Glct[FixedPoint[Expand,x],l],Plus,Glct,Plus,Glc];
GLC[x_,l_List:{{}}]:=Distribute[Glct[FixedPoint[Expand,x],l],Plus,Glct,Plus,Glc];

Glc[x_]:=Glc[x,{{}}];
Glc[Unevaluated[ NonCommutativeMultiply[x__]v_],l_List]:=v Glc[NonCommutativeMultiply[x],Join[l,Union[DeleteCases[AllIndexLists[v ],{_}]]]];
Glc[GammaProd[x__]v_.,_List]:=GammaProd[x]v/;(ReplaceList[{x},{___List,{___, y_,___},__List,{___,y_,___},___List}:>y]==={});

Glc[Unevaluated[NonCommutativeMultiply[x___,GammaProd[z__],y___]],l_List]:=NonCommutativeMultiply[x,Glc[GammaProd[z],Join[l,Union[DeleteCases[AllIndexLists[x y],{_}]]]],y];
Glc[GammaProd[z___,x_List,m_List,y_List,w___]v_.,l_List]:=Module[{inter,yyleft,k,n,r,p,mmleft,ltotal=Join[l,{z},{w},Union[DeleteCases[AllIndexLists[v ],{_}]]],mmtotal,yytotal},inter=Intersection[x,y];
mmleft=m;
yyleft=Complement[y,inter];(*j1...jl*)
n=Length[m];r=Length[yyleft];k=Length[inter];
Signature[Join[inter,yyleft]]Sum[Module[{mmcomb=Subsets[m,{p}],yycomb=Subsets[yyleft,{p}],actmmcompo,actyycompo,mmcompo,yycompo},
mmcompo=Map[Complement[mmleft,#]&,mmcomb];yycompo=Map[Complement[yyleft,#]&,yycomb];mmtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[mmcomb[[i]],mmcompo[[i]]]],{i,Length[mmcomb]}]];yytotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[yycomb[[i]],yycompo[[i]]]],{i,Length[yycomb]}]];
c[k,n,r,p] v Apply[Plus,Flatten[Outer[Count[mmtotal,#1]Count[yytotal,#2]Delta[Take[#1,p],Take[#2,p]]GammaProd[z,Join[actmmcompo=Drop[#1,p],actyycompo=Drop[#2,p]],w]Signature[#1]Signature[#2]&,Union[mmtotal],Union[yytotal],1]]]],{p,0,Min[n,r]}]]/;Intersection[x,y]===x;
Glc[GammaProd[z___,x_List,m_List,y_List,w___]v_.,l_List]:=Module[{inter,xxleft,k,n,r,p,mmleft,ltotal=Join[l,{z},{w},Union[DeleteCases[AllIndexLists[v ],{_}]]],mmtotal,xxtotal},inter=Intersection[x,y];
xxleft=Complement[x,inter];(*j1...jl*)mmleft=m;
n=Length[m];r=Length[xxleft];k=Length[inter];
Signature[Join[inter,xxleft]]Sum[Module[{mmcomb=Subsets[m,{p}],xxcomb=Subsets[xxleft,{p}],xxcompo,mmcompo},mmcompo=Map[Complement[mmleft,#]&,mmcomb];xxcompo=Map[Complement[xxleft,#]&,xxcomb];mmtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[mmcomb[[i]],mmcompo[[i]]]],{i,Length[mmcomb]}]];xxtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[xxcomb[[i]],xxcompo[[i]]]],{i,Length[xxcomb]}]];d[k,n,r,p] v Apply[Plus,Flatten[Outer[Count[mmtotal,#2]Count[xxtotal,#1]Delta[Take[#1,p],Take[#2,p]]GammaProd[z,Join[Drop[#1,p],Drop[#2,p]],w]Signature[#1]Signature[#2]&,Union[xxtotal],Union[mmtotal],1]]]],{p,0,Min[n,r]}]]/;Intersection[x,y]===y;
Glc[GammaProd[z___,x_List,m_List,y_List,w___]v_.,l_List]:=Module[{inter,xxleft,yyleft,lenxxleft,lenyyleft,k,p,ltotal=Join[l,{z},{w},Union[DeleteCases[AllIndexLists[v ],{_}]]]},inter=Intersection[x,y];k=Length[inter];xxleft=Complement[x,inter];yyleft=Complement[y,inter];
If[(lenxxleft=Length[xxleft])<=(lenyyleft=Length[yyleft]),
Signature[Join[xxleft,inter]]Sum[Module[{xxcomb=Subsets[xxleft,{p}],xxcompo,xxtotal},
xxcompo=Map[Complement[xxleft,#]&,xxcomb];
xxtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[xxcomb[[i]],xxcompo[[i]]]],{i,Length[xxcomb]}]];Coefficient[G[lenxxleft,k],G1[lenxxleft-p]G2[k-p]] v /Length[xxcomb]Apply[Plus,Map[Count[xxtotal,#]Signature[#]Delta[Take[#,p],Take[inter,p]]NonCommutativeMultiply[GammaProd[z,Drop[#,p]],GammaLongContract[GammaProd[Drop[inter,p],m,y],ltotal],GammaProd[w]]&,Union[xxtotal]]]],{p,0,Min[k,lenxxleft]}],
Signature[Join[inter,yyleft]] Sum[Module[{yycomb=Subsets[yyleft,{p}],yycompo,yytotal},yycompo=Map[Complement[yyleft,#]&,yycomb];
yytotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[yycomb[[i]],yycompo[[i]]]],{i,Length[yycomb]}]];Coefficient[G[k,lenyyleft],G1[k-p]G2[lenyyleft-p]] v /Length[yycomb]Apply[Plus,Map[Count[yytotal,#]Signature[#]Delta[Take[#,p],Take[inter,p]]NonCommutativeMultiply[GammaProd[z],GammaLongContract[GammaProd[x,m,Drop[inter,p]],ltotal],GammaProd[Drop[#,p],w]]&,Union[yytotal]]]],{p,0,Min[k,lenyyleft]}]]]/;Intersection[x,y]=!={};
Glc[GammaProd[z___,x_List,m_List,n__List,y_List,w___]v_.,l_List]:=v GammaLongContract[GammaProd[z,x]**GammaExpand[GammaProd[m,n],Join[l,{z},{x},{y},{w},Union[DeleteCases[AllIndexLists[v],{_}]]]]**GammaProd[y,w],Join[l,Union[DeleteCases[AllIndexLists[v],{_}]]]]/;(Intersection[x,y]=!={});
Glc[x_/;FreeQ[x,GammaProd],_List]:=x;


c[n_,k_,l_,q_]:=Power[-1,Evaluate[((q(2k-1)-q^2+n(2k+n-1))/2)]]q!n!Sum[Power[-1,Evaluate[p]] Binomial[Dim+q-k-l,n-p]Binomial[k-q,p],{p,0,Min[n,k]}];
d[n_,k_,l_,q_]:=Power[-1,Evaluate[((q(2k-1)-q^2+n(2k+n-1))/2)]]q!n!Sum[Power[-1,Evaluate[(p +(4+k^2+l^2-3 n+n^2+l (-3+n-3 q)+3 q+3 q^2+k (2 l-3 (1+q))))]]Binomial[Dim+q-k-l,n-p]Binomial[k-q,p],{p,0,Min[n,k]}];


(* GammaExpand takes a list of lists as the second argument *)
GammaExpand[x_,l_List:{{}}]:=Distribute[Gett[FixedPoint[Expand,x],l],Plus,Gett,Plus,Ge];
GE[x_,l_List:{{}}]:=Distribute[Gett[FixedPoint[Expand,x],l],Plus,Gett,Plus,Ge];
Ge[x_]:=Ge[x,{{}}];

Ge[HoldPattern[NonCommutativeMultiply[x___,GammaProd[y__List],z___]v_.],l_List]:=GE[v] NonCommutativeMultiply[x,Ge[GammaProd[y],Join[l,Union[DeleteCases[AllIndexLists[v x z],{_}]]]],z];
Ge[GammaProd[x_List,y_List,z___] v_.,l_List]:=Module[{lenxx=Length[x],lenyy=Length[y],ltotal=Join[l,{z},Union[DeleteCases[AllIndexLists[v ],{_}]]],xxcompo,yycompo,xxtotal,yytotal},
	Sum[Module[{xxcomb=Subsets[x,{p}],yycomb=Subsets[y,{p}]},
		xxcompo=Map[Complement[x,#]&,xxcomb];
		yycompo=Map[Complement[y,#]&,yycomb];
		xxtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[xxcomb[[i]],xxcompo[[i]]]],{i,Length[xxcomb]}]];
		yytotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[yycomb[[i]],yycompo[[i]]]],{i,Length[yycomb]}]];
		Power[-1,Evaluate[(lenxx p-p (p+1)/2)]]p! v Apply[Plus,Flatten[Outer[Count[xxtotal,#1]Count[yytotal,#2]Delta[Take[#1,p],Take[#2,p]]GammaProd[Join[Drop[#1,p],Drop[#2,p]],z]Signature[#1]Signature[#2]&,Union[xxtotal],Union[yytotal],1]]]],
	{p,0,Min[lenxx,lenyy]}]
	];
Ge[GammaProd[x_List] v_.,_List]:=v GammaProd[x];
Ge[x_/;FreeQ[x,GammaProd],_List]:=x;

Ge[x_]:=Module[{},Message[GammaExpand::arg];Print[x];Abort[]];

GammaExpand::arg="Unsupported argument.";


(* GammaTrace takes a list of lists as the second argument *)
GammaTrace[x_,{y__Symbol}]:=GT[x,{{y}}];
GammaTrace[x_,asym_List:{{}}]:=Distribute[Gtt[FixedPoint[Expand,x],asym],Plus,Gtt,Plus,Gt];
GT[x_,asym_List:{{}}]:=Distribute[Gtt[FixedPoint[Expand,x],asym],Plus,Gtt,Plus,Gt];

Gt[GammaProd[x_List,y_List]v_.,asym_List]:=Module[{lenxx,lenyy},If[(lenxx=Length[x])==(lenyy=Length[y]),SpinorDim Power[-1,Evaluate[(lenxx (lenxx-1)/2)]]lenxx! ACanonicalOrder[ v Delta[x,y],asym],0]+If[(OddQ[Dim]||ChiralSpinors)&&lenxx+lenyy>=Dim&&EvenQ[lenxx+lenyy-Dim],
Module[{p=(lenxx+lenyy-Dim)/2,xxcomb,yycomb,xxcompo,yycompo,xxtotal,yytotal},ltotal=Join[asym,Union[DeleteCases[AllIndexLists[v ],{_}]]];xxcomb=Subsets[x,{p}];yycomb=Subsets[y,{p}];xxcompo=Map[Complement[x,#]&,xxcomb];yycompo=Map[Complement[y,#]&,yycomb];xxtotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[xxcomb[[i]],xxcompo[[i]]]],{i,Length[xxcomb]}]];yytotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[yycomb[[i]],yycompo[[i]]]],{i,Length[yycomb]}]];Power[-1,Evaluate[(lenxx p-p (p+1)/2)]]p! v Apply[Plus,Flatten[Outer[Count[xxtotal,#1]Count[yytotal,#2]Delta[Take[#1,p],Take[#2,p]]SpinorDim EpsilonConvention Epsilon[Join[Drop[#1,p],Drop[#2,p]]]Signature[#1]Signature[#2]&,Union[xxtotal],Union[yytotal],1]]]],0]];
Gt[GammaProd[x_List,y_List,z__List] v_.,asym_List]:=Module[{len=Length[Flatten[{x,y,z}]]},If[(EvenQ[len])||(EvenQ[len-Dim]&&len>=Dim),GammaTrace[ACanonicalOrder[GammaContract[ACanonicalOrder[GammaExpand[v GammaProd[x,y,z],asym],asym],asym],asym],asym],0]];
(* In case of tracing a Subscript[\[CapitalGamma], (D)] when D is odd just keep it since there are different conventions (a sign) *)
Gt[GammaProd[x_List]v_.,_List]:=If[Length[x]==Dim&&(OddQ[Dim]||ChiralSpinors),v SpinorDim EpsilonConvention Epsilon[x],0];
Gt[x_/;FreeQ[x,GammaProd],_List]:=SpinorDim x;

Epsilon[{___,x_,x_,___}]=0;
Epsilon[x_List]:=Signature[x]Epsilon[Sort[x]]/;\[Not]OrderedQ[x];
Epsilon/:Epsilon[x_List]Epsilon[y_List]:=SpaceTimeSignature Length[x]!Delta[x,y]

Format[Epsilon[x_List],{StandardForm,TeXForm}]:=DisplayForm[SubscriptBox[Global`\[Epsilon],RowBox[x]]];

Format[GammaProd[x_List],{StandardForm,TeXForm}]:=DisplayForm[SubscriptBox[Global`\[CapitalGamma],RowBox[x/.{Global`$minus->"- ",Global`$plus->"+ "}]]];
Format[GammaProd[x__List],{StandardForm,TeXForm}]:=DisplayForm[NonCommutativeMultiply[Apply[Sequence,Map[SubscriptBox[Global`\[CapitalGamma],RowBox[#/.{Global`$minus->"- ",Global`$plus->"+ "}]]&,{x}]]]];

MakeExpression[SubscriptBox[y_,n_],StandardForm]:=Module[{indices=RowBox/@DeleteCases[Split[(Flatten[{n}/.RowBox->List]/.{"+"->"$plus","-"->"$minus"}),#2=!=","&],",",2]},If[Length[indices]===2,indices=Insert[indices,",",2]];MakeExpression[RowBox[{"Subscript[",y,",",Sequence@@indices,"]"}],StandardForm]];

Unprotect[Subscript];
Subscript[Global`\[Epsilon],x_]:=Epsilon[Unevaluated[x]/.Times->List];
Subscript[Global`\[CapitalGamma],x_Symbol|x_Integer]:=GammaProd[{x}];
Subscript[Global`\[CapitalGamma],$plus]:=GammaProd[{}];
Subscript[Global`\[CapitalGamma],x_]:=GammaProd[Unevaluated[x]/.Times->List];
Protect[Subscript];


(* Tensor *)
(* Try to generalise to more general types of tensors *)

Tensor/:Tensor[sym_,___List,z_List,___List](f_)[{v___,x_,w___}] (f_)[{v___,y_,w___}]:=0/;(MemberQ[z,x]&&MemberQ[z,y]);
Tensor/:Tensor[sym_,x___,y_Symbol,z___]:=Tensor[sym,x,{y},z];

Tensor[_,___List,{___,x_,___,x_,___},___List]:=0;
(*Tensor[y_Symbol/;y=!=Global`\[CapitalOmega],___List,{___,x_,___},___List,{___,x_,___},___List]:=0;*)
Tensor[sym_,y___List,x_List,z___List]:=Signature[x] Apply[Tensor,{sym,y,Sort[x],z}]/; !OrderedQ[x];

SetTensorInd[x_]:=Module[{},TensorIndicator=x;];

Format[Tensor[sym_,x__List],{StandardForm,TeXForm}]:=DisplayForm[Subscript[TensorIndicator[sym],Evaluate[Sequence@@Map[RowBox,{x}]]]];

Unprotect[Subscript];
Subscript[TensorIndicator[sym_/;(sym=!=Global`\[Delta]&&sym=!=Global`\[CapitalGamma])],x__/;FreeQ[{x},RowBox]]:=Tensor[sym,(Unevaluated[x]/.Times->List)];
Protect[Subscript];

TensorCanonicalOrder[x_,ind_List]:=Distribute[TCOt[FixedPoint[Expand,x],ind],Plus,TCOt,Plus,TCO];
TC[x_,ind_List]:=Distribute[TCOt[FixedPoint[Expand,x],ind],Plus,TCOt,Plus,TCO];

(* (Use recursion for simplicity since speed is (hopefully) not important for this command) *)
TCO[v_,_List]:=v/;FreeQ[v,Tensor];
(* The end condition of the recursion *)
TCO[Tensor[s_Symbol,x__List]v_.,ind_List]:=Tensor[s,x]v/;{x}===ind;

(* The index lists must be sorted! *)
TCO[Tensor[s_Symbol,x___List,a_List,b___List,{c1___,y2_,c2___},d___List]v_.,{x___List,{y1___,y2_,y3___},z___List}]:=v Module[{inplace=Intersection[a,{y1,y2,y3}],outofplace,lenoutofplace},outofplace=Complement[a,inplace];lenoutofplace=Length[outofplace];Power[-1,Evaluate[lenoutofplace+1+Length[{c1}]]] lenoutofplace/(Length[inplace]+1) Signature[Flatten[{inplace,outofplace}]]TC[ASym[Tensor[s,x,Flatten[{inplace,y2,Drop[outofplace,-1]}],b,{Last[outofplace],c1,c2},d],outofplace],{x,{y1,y2,y3},z}]];

TensorProject[x_,ind__List]:=Distribute[TPRt[FixedPoint[Expand,x],ind],Plus,TPRt,Plus,TPR];
TP[x_,ind__List]:=Distribute[TPRt[FixedPoint[Expand,x],ind],Plus,TPRt,Plus,TPR];
TP2[x_,ind__List]:=Distribute[TPR2t[FixedPoint[Expand,x],ind],Plus,TPR2t,Plus,TPR2];


(* Recursive implementation *)
(*
TPR[0,__List]:=0;
TPR[x_,{},__List]:=x;
TPR[x_ y_./;FreeQ[y,_List],a_List,b_List]:=y Module[{ajoinb=Join[a,b],minlen=Min[Length[a],Length[b]],size=AllSize[x],ind=AllIndices[x],intlist,d},intlist=(Range[Length[ajoinb]]/.d_Integer:>List[d]);ASym[ASym[Sum[Binomial[Length[b],k]Reconstruct[size,ReplacePart[ind,Join[Drop[a,-k],Take[b,k],Take[a,-k],Drop[b,k]],PositionList[ind,ajoinb],intlist]],{k,0,minlen}],a],b]];
TPR[x_ y_./;FreeQ[y,_List],a_List,b_List,c__List]:=y TP2[TPR[x,b,c],a,b,c];

TPR2[0,__List]:=0;
TPR2[x_,{},__List]:=x;
TPR2[x_ y_./;FreeQ[y,_List],a_List,b_List]:=y TPR[x,a,b];
TPR2[x_ y_./;FreeQ[y,_List],a_List,b_List,c__List]:=y Module[{size=AllSize[x],ajoinb=Join[a,b],ind=AllIndices[x],minlen=Min[Length[a],Length[b]],intlist,d},intlist=(Range[Length[ajoinb]]/.d_Integer:>List[d]);ASym[ASym[Sum[Binomial[Length[b],k]TPR2[Reconstruct[size,ReplacePart[ind,Join[Drop[a,-k],Take[b,k],Take[a,-k],Drop[b,k]],PositionList[ind,ajoinb],intlist]],Take[b,k],c],{k,0,minlen}],a],b]];
*)

(*

(* non-recursive implementation for the n-tensor case *)

TPR[x_ y_./;FreeQ[y,_List],a__List]:=y Module[{nolists=Length[{a}],size=AllSize[x],ind=AllIndices[x]},Apply[Plus,Map[#[[1]]Reconstruct[size,#[[2]]]&
,Fold[Function[{structs,blocks},Flatten[Map[TPreshuffle[#1,Flatten[PositionList[ind,{a}[[blocks[[1]]]]],1],Flatten[PositionList[ind,{a}[[blocks[[2]]]]],1]]&,structs],1]],{{1,ind}},Reverse[Flatten[Table[{l,k},{k,1,Length[{a}]},{l,1,k-1}],1]]]]]];

TPreshuffle[{x_,str_List},a_List,b_List]:=Module[{ajoinb=Join[a,b],ajoinblist,intlist,d},intlist=(Range[Length[ajoinb]]/.d_Integer:>List[d]);ajoinblist=(ajoinb/.d_Integer:>List[d]);Flatten[Map[Module[{k=#,asub=Subsets[a,{#}],bsub=Subsets[b,{#}],acomp,bcomp,weight=Binomial[Length[a],#],d},acomp=Map[Complement[a,#]&,asub];bcomp=Map[Complement[b,#]&,bsub];Map[{Signature[str[[a]]]Signature[str[[b]]]Signature[Join[Drop[#[[1]],k],Take[#[[1]],k]]]Signature[Join[Take[#[[2]],k],Drop[#[[2]],k]]]x/weight,ReplacePart[str,str[[Join[Drop[#[[1]],k],Take[#[[2]],k],Take[#[[1]],k],Drop[#[[2]],k]]]],ajoinblist,intlist]}&,Flatten[Outer[List,MapThread[Join,{asub,acomp}],MapThread[Join,{bsub,bcomp}],1],1]]]&,Range[0,Length[b]]],1]];

*)

(* 3-tensor case non-recursive implementation *)

TPR[x_,a_List,b_List,e_List]:=Module[{p=Length[a],q=Length[b],r=Length[e],size=AllSize[x],indices=AllIndices[x],pos},pos=PositionList[indices,Join[a,b,e]];c[q,0,0,r,p,q,r]=1;Rd[Expand[Sum[If[c[q2,q3,r2,r3,p,q,r]!=0,c[q2,q3,r2,r3,p,q,r]Reconstruct[size,ReplacePart[indices,Join[Take[a,p1[q2,q3,r2,r3,p,q,r]],Take[b,q1[q2,q3,r2,r3,p,q,r]],Take[e,r1[q2,q3,r2,r3,p,q,r]],Take[a,{p1[q2,q3,r2,r3,p,q,r]+1,p1[q2,q3,r2,r3,p,q,r]+p2[q2,q3,r2,r3,p,q,r]}],Take[b,{q1[q2,q3,r2,r3,p,q,r]+1,q1[q2,q3,r2,r3,p,q,r]+q2}],Take[e,{r1[q2,q3,r2,r3,p,q,r]+1,r1[q2,q3,r2,r3,p,q,r]+r2}],Take[a,-p3[q2,q3,r2,r3,p,q,r]],Take[b,-q3],Take[e,-r3]],pos,(Range[p+q+r]/.f_Integer:>List[f])]],0],{r3,0,r},{r2,0,r-r3},{q3,0,q},{q2,0,q-q3}]],(#1[[1]][[1]]<#2[[1]][[1]])&]];

p1[q2_,q3_,r2_,r3_,p_,q_,r_]:=p-q+q2+q3-r+r2+r3;
p2[q2_,q3_,r2_,r3_,p_,q_,r_]:=q-q2-r2;
p3[q2_,q3_,r2_,r3_,p_,q_,r_]:=-q3+r-r3;
q1[q2_,q3_,r2_,r3_,p_,q_,r_]:=q-q2-q3;
r1[q2_,q3_,r2_,r3_,p_,q_,r_]:=r-r2-r3;

(* recursive implementation of c(q2,q3,r2,r3) *)

c[___,x_/;x<0,___]:=0;
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;q1[q2,q3,r2,r3,p,q,r]+q2+q3-q>0
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;r1[q2,q3,r2,r3,p,q,r]+r2+r3-r>0
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;p2[q2,q3,r2,r3,p,q,r]+q2+r2-q>0
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;p3[q2,q3,r2,r3,p,q,r]+q3+r3-r>0
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;p1[q2,q3,r2,r3,p,q,r]+p2[q2,q3,r2,r3,p,q,r]+p3[q2,q3,r2,r3,p,q,r]-p>0
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;p1[q2,q3,r2,r3,p,q,r]+q1[q2,q3,r2,r3,p,q,r]+r1[q2,q3,r2,r3,p,q,r]-p>0
c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=0)/;(p1[q2,q3,r2,r3,p,q,r]<0||p2[q2,q3,r2,r3,p,q,r]<0||p3[q2,q3,r2,r3,p,q,r]<0||q1[q2,q3,r2,r3,p,q,r]<0||r1[q2,q3,r2,r3,p,q,r]<0)

c[q2_,q3_,r2_/;r2>0,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=1/r2 (-Power[-1,Evaluate[q2+r1[q2,q3,r2,r3,p,q,r]-1]](r1[q2,q3,r2,r3,p,q,r]+1)c[q2+1,q3,r2-1,r3,p,q,r]-Power[-1,Evaluate[q3+r2+1]](r3+1) c[q2+1,q3-1,r2-1,r3+1,p,q,r]));

c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=1/q1[q2,q3,r2,r3,p,q,r] (-Power[-1,Evaluate[p2[q2,q3,r2,r3,p,q,r]+p3[q2,q3,r2,r3,p,q,r]+q1[q2,q3,r2,r3,p,q,r]+q2+1]](q3+1) c[q2,q3+1,r2,r3,p,q,r]-Power[-1,Evaluate[p2[q2,q3,r2,r3,p,q,r]+q1[q2,q3,r2,r3,p,q,r]+1]](q2+1)c[q2+1,q3,r2,r3,p,q,r]))/;q1[q2,q3,r2,r3,p,q,r]>0;

c[q2_,q3_,r2_,r3_,p_,q_,r_]:=(c[q2,q3,r2,r3,p,q,r]=1/r1[q2,q3,r2,r3,p,q,r] (-Power[-1,Evaluate[p2[q2,q3,r2,r3,p,q,r]+r1[q2,q3,r2,r3,p,q,r]+1+q1[q2,q3,r2,r3,p,q,r]+q2]](r2+1) c[q2,q3,r2+1,r3,p,q,r]-Power[-1,Evaluate[p2[q2,q3,r2,r3,p,q,r]+p3[q2,q3,r2,r3,p,q,r]+r1[q2,q3,r2,r3,p,q,r]+r2+1+q1[q2,q3,r2,r3,p,q,r]+q3]](r3+1)c[q2,q3,r2,r3+1,p,q,r]))/;r1[q2,q3,r2,r3,p,q,r]>0;

(* 4-tensor case non-recursive implementation *)

TPR[x_,a_List,b_List,e_List,f_List]:=Module[{p=Length[a],q=Length[b],r=Length[e],s=Length[f],size=AllSize[x],indices=AllIndices[x],pos},pos=PositionList[indices,Join[a,b,e,f]];c[q,0,0,0,r,0,0,0,s,p,q,r,s]=1;
Rd[Expand[Sum[If[c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]!=0,c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]Reconstruct[size,ReplacePart[indices,Join[Take[a,p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[b,q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[e,r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[f,s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[a,{p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]}],Take[b,{q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2}],Take[e,{r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2}],Take[f,{s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2}],
Take[a,{p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]}],Take[b,{q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+1,q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+q3}],Take[e,{r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+1,r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+r3}],Take[f,{s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+1,s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+s3}],Take[a,-p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[b,-q4],Take[e,-r4],Take[f,-s4]],pos,(Range[p+q+r+s]/.g_Integer:>List[g])]],0],{s4,0,s},{s3,0,s-s4},{s2,0,s-s3-s4},{r4,0,r},{r3,0,r-r4},{r2,0,r-r3-r4},{q4,0,q},{q3,0,q-q4},{q2,0,q-q3-q4}]],(#1[[1]][[1]]<#2[[1]][[1]])&]];

p1[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=p-q+q2+q3+q4-r+r2+r3+r4-s+s2+s3+s4
p2[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=q-q2-r2-s2
p3[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=-q3+r-r3-s3
p4[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=-q4-r4+s-s4
q1[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=q-q2-q3-q4
r1[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=r-r2-r3-r4
s1[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=s-s2-s3-s4

(* recursive implementation of c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s] *)

c[___,x_/;x<0,___]:=0;

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]-p>0
c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+q3+q4-q>0
c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+r3+r4-r>0
c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+s3+s4-s>0

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]-p>0
c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+r2+s2-q>0
c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q3+r3+s3-r>0
c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q4+r4+s4-s>0

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=0)/;(p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0||p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0||p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0||p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0||q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0||r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0||s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<0)

c[q2_,q3_,q4_,r2_/;r2>0,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=1/r2 (-Power[-1,Evaluate[q2+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]-1]](r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1)c[q2+1,q3,q4,r2-1,r3,r4,s2,s3,s4,p,q,r,s]-Power[-1,Evaluate[q3+r2+1]](r3+1) c[q2+1,q3-1,q4,r2-1,r3+1,r4,s2,s3,s4,p,q,r,s]-Power[-1,q3+q4+r2+r3+1](r4+1) c[q2+1,q3,q4-1,r2-1,r3,r4+1,s2,s3,s4,p,q,r,s]));


c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_/;s3>0,s4_,p_,q_,r_,s_]:=( c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=1/s3 (-Power[-1,Evaluate[r2+r3+s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2-1]](s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1)c[q2,q3,q4,r2,r3+1,r4,s2,s3-1,s4,p,q,r,s]-Power[-1,Evaluate[r3+s2-1]](s2+1) c[q2,q3,q4,r2-1,r3+1,r4,s2+1,s3-1,s4,p,q,r,s]-Power[-1,Evaluate[r4+s3+1]](s4+1) c[q2,q3,q4,r2,r3+1,r4-1,s2,s3-1,s4+1,p,q,r,s]));

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=( c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=1/s2 (-Power[-1,Evaluate[q2+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]-1]](s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1)c[q2+1,q3,q4,r2,r3,r4,s2-1,s3,s4,p,q,r,s]-Power[-1,Evaluate[q3+r2+r3+1+s2]](s3+1)c[q2+1,q3-1,q4,r2,r3,r4,s2-1,s3+1,s4,p,q,r,s]-Power[-1,Evaluate[q3+q4+r2+r4+s2+s3+1]](s4+1)c[q2+1,q3,q4-1,r2,r3,r4,s2-1,s3,s4+1,p,q,r,s]))/;s2>0;

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=1/q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s] (-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+q3+1]](q4+1) c[q2,q3,q4+1,r2,r3,r4,s2,s3,s4,p,q,r,s]-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+1]](q3+1) c[q2,q3+1,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1]](q2+1)c[q2+1,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]))/;q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>0;

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=1/r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s] (-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2]](r2+1) c[q2,q3,q4,r2+1,r3,r4,s2,s3,s4,p,q,r,s]-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+1+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q3]](r3+1)c[q2,q3,q4,r2,r3+1,r4,s2,s3,s4,p,q,r,s]-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+r3+1+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q4]](r4+1)c[q2,q3,q4,r2,r3,r4+1,s2,s3,s4,p,q,r,s]))/;r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>0;

c[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=(c[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]=1/s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s] (-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1]](s2+1) c[q2,q3,q4,r2,r3,r4,s2+1,s3,s4,p,q,r,s]-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+1+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q3+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r3]](s3+1)c[q2,q3,q4,r2,r3,r4,s2,s3+1,s4,p,q,r,s]-Power[-1,Evaluate[p2[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p3[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+p4[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+s3+1+q1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q4+r1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r4]](s4+1)c[q2,q3,q4,r2,r3,r4,s2,s3,s4+1,p,q,r,s]))/;s1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>0;

RemoveTrace[x_,a_List,b_List]:=Module[{lena=Length[a],lenb=Length[b],dind=Intersection[ddList,AllIndices[x]],offset},offset=Max[{0},Flatten[PositionList[ddList,dind]]];x+Expand[Sum[ASym[ASym[Power[-1,Evaluate[i]] (Binomial[lena,i]Binomial[lenb,i])/Binomial[Dim-lena-lenb+i+1,i] Delta[Take[a,i],Take[b,i]]Rd[Delta[Take[a,i],Take[b,i]]x,(#1[[1]][[1]]<#2[[1]][[1]])&,offset],a],b],{i,1,lenb}]]];

(* NB: the output is _not_ antisymmetrised! *)

RemoveTrace[x_,a_List,b_List,c_List]:=Distribute[RTt[FixedPoint[Expand,x],a,b,c],Plus,RTt,Plus,RT];


(*
(* 3-tensor case recursive implementation *)

(* there should probably be another ASym in the first block of reshuffled indices *)
RT[x_,a_List,b_List,c_List]:=Module[{p=Length[a],q=Length[b],r=Length[c],dind=Intersection[ddList,Join[a,b,c]],size=AllSize[x],indices=AllIndices[x],offset,pos,lastind},offset=Max[{0},Flatten[PositionList[ddList,dind]]];pos=PositionList[indices,Join[a,b,c]];Rd[Expand[Sum[Sum[If[(Dim-p-q-r+cac+1)\[GreaterEqual]0,lastind=Join[Take[b,-r2],Drop[c,r2]];ASym[Delta[Take[a,cac],Take[lastind,cac]](RemoveTrace[Reconstruct[size,ReplacePart[indices,Join[Take[ddList,{offset+1,offset+cac}],Drop[a,cac],Drop[b,-r2],Take[c,r2],Take[ddList,{offset+1+cac,offset+2cac}],Drop[lastind,cac]],pos,(Range[p+q+r]/.e_Integer:>List[e])]],Join[Take[ddList,{offset+1,offset+cac}],Drop[a,cac]],Join[Drop[b,-r2],Take[c,r2]]]/.Inner[Rule,Take[ddList,{offset+cac+1,offset+2cac}],Take[ddList,{offset+1,offset+cac}],List])Power[-1,Evaluate[r2^2+cac]]Binomial[q,r2]Binomial[r,r2]Binomial[r,cac]Binomial[p,cac]cac!Max[1,Dim-p-q-r+1]!/(Dim-p-q-r+cac+1)!,lastind],0],{r2,0,r}],{cac,0,Min[p,r]}]],(#1[[1]][[1]]<#2[[1]][[1]])&,offset]];
*)


(* NB: the output is _not_ antisymmetrised! *)

(* 3-tensor case non-recursive implementation *)

RT[x_,a_List,b_List,c_List]:=Module[{p=Length[a],q=Length[b],r=Length[c],dind=Intersection[ddList,Join[a,b,c]],size=AllSize[x],indices=AllIndices[x],offset,pos},offset=Max[{0},Flatten[PositionList[ddList,dind]]];pos=PositionList[indices,Join[a,b,c]];Rd[Expand[Sum[Sum[Sum[Sum[If[iac[iab,p1,q2,r2,p,q,r]>=0&&iac[iab,p1,q2,r2,p,q,r]<=Min[p,r]&&cab[iab,p1,q2,r2,p,q,r]>=0&&cab[iab,p1,q2,r2,p,q,r]<=Min[p,q]&&cac[iab,p1,q2,r2,p,q,r]>=0&&cac[iab,p1,q2,r2,p,q,r]<=Min[p,r]&&q3[iab,p1,q2,r2,p,q,r]>=0&&q3[iab,p1,q2,r2,p,q,r]<=q&&r3[iab,p1,q2,r2,p,q,r]>=0&&r3[iab,p1,q2,r2,p,q,r]<=r&&((Dim-p-q-r+2)>=1||(Dim-p-q-r+iab+iac[iab,p1,q2,r2,p,q,r]+1)>=0),Delta[Take[a,iab+iac[iab,p1,q2,r2,p,q,r]],Join[Take[b,iab],Take[c,iac[iab,p1,q2,r2,p,q,r]]]]Reconstruct[size,ReplacePart[indices,Join[Take[ddList,{offset+1,offset+cab[iab,p1,q2,r2,p,q,r]}],Take[ddList,{offset+cab[iab,p1,q2,r2,p,q,r]+1,offset+cab[iab,p1,q2,r2,p,q,r]+cac[iab,p1,q2,r2,p,q,r]}],Drop[a,iab+iac[iab,p1,q2,r2,p,q,r]],Take[ddList,{offset+1,offset+cab[iab,p1,q2,r2,p,q,r]}],Take[b,{iab+1,iab+q2}],Take[c,{iac[iab,p1,q2,r2,p,q,r]+1,iac[iab,p1,q2,r2,p,q,r]+r2}],Take[ddList,{offset+cab[iab,p1,q2,r2,p,q,r]+1,offset+cab[iab,p1,q2,r2,p,q,r]+cac[iab,p1,q2,r2,p,q,r]}],Take[b,-q3[iab,p1,q2,r2,p,q,r]],Take[c,-r3[iab,p1,q2,r2,p,q,r]]],pos,(Range[p+q+r]/.e_Integer:>List[e])]]Power[-1,Evaluate[iab(1+p+p1+q+r2)+(q+q2+r2+1) p1+q2(p+q+r2+1)+(p+1) r2+p+q+p q+q^2]]p! q! r!Binomial[iab+iac[iab,p1,q2,r2,p,q,r],iab]/(p1!q2!q3[iab,p1,q2,r2,p,q,r]!r2!r3[iab,p1,q2,r2,p,q,r]!cab[iab,p1,q2,r2,p,q,r]!cac[iab,p1,q2,r2,p,q,r]!(Range[Max[1,(Dim-p-q-r+2)],(Dim-p-q-r+iab+iac[iab,p1,q2,r2,p,q,r]+1)]/.List->Times)),0],{r2,0,r}],{q2,0,q-iab}],{p1,0,p-iab}],{iab,0,Min[p,q]}]],(#1[[1]][[1]]<#2[[1]][[1]])&,offset]];

iac[iab_,p1_,q2_,r2_,p_,q_,r_]:=-iab+p-p1;
cab[iab_,p1_,q2_,r2_,p_,q_,r_]:=q-q2-r2;
cac[iab_,p1_,q2_,r2_,p_,q_,r_]:=p-p1-q+q2+r2;
q3[iab_,p1_,q2_,r2_,p_,q_,r_]:=-iab+q-q2;
r3[iab_,p1_,q2_,r2_,p_,q_,r_]:=iab-p+p1+r-r2; 

(* 4-tensor case, non-recursive implementation *)

RemoveTrace[x_,a_List,b_List,c_List,d_List]:=Distribute[RTt[FixedPoint[Expand,x],a,b,c,d],Plus,RTt,Plus,RT];

RT[x_,a_List,b_List,c_List,d_List]:=Module[{p=Length[a],q=Length[b],r=Length[c],s=Length[d],dind=Intersection[ddList,Join[a,b,c,d]],size=AllSize[x],indices=AllIndices[x],offset,pos},offset=Max[{0},Flatten[PositionList[ddList,dind]]];pos=PositionList[indices,Join[a,b,c,d]];Rd[Expand[Sum[Sum[Sum[Sum[Sum[Sum[Sum[Sum[Sum[If[iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=Min[p,q]&&iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=Min[p,r]&&iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=Min[p,s]&&cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=Min[p,q]&&cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=Min[p,r]&&cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=Min[p,s]&&p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]>=0&&p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]<=p&&((Dim-p-q-r-s+2)>=1||(Dim-p-q-r-s+iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1)>=0),Delta[Take[a,iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Join[Take[b,iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[c,iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[d,iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]]]]Reconstruct[size,ReplacePart[indices,Join[Take[ddList,{offset+1,offset+cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]}],Drop[a,iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]],Take[ddList,{offset+1,offset+cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]}],Take[b,{iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2}],Take[c,{iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2}],Take[d,{iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2}],Take[ddList,{offset+cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,offset+cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]}],Take[b,{iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+1,iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+q2+q3}],Take[c,{iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+1,iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r2+r3}],Take[d,{iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+1,iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+s2+s3}],Take[ddList,{offset+cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1,offset+cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]}],Take[b,-q4],Take[c,-r4],Take[d,-s4]],pos,(Range[p+q+r+s]/.e_Integer:>List[e])]]Power[-1,Evaluate[(iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s])(q2+q3+q4)+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s](r2+r3+r4)+(q2+r2+s2)(cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s])+(q3+r3+s3) cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+r3 q4+r2(q3+q4)+s3(q4+r4)+s2(q3+q4+r3+r4)+iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]]]p! q! r!s!Multinomial[iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s],iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s],iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]]/(p1[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]!q2!q3!q4!r2!r3!r4!s2!s3!s4!cab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]!cac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]!cad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]!(Range[Max[1,(Dim-p-q-r-s+2)],(Dim-p-q-r-s+iab[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iac[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+iad[q2,q3,q4,r2,r3,r4,s2,s3,s4,p,q,r,s]+1)]/.List->Times)),0],{s4,0,s-s2-s3}],{s3,0,s-s2}],{s2,0,s}],{r4,0,r-r2-r3}],{r3,0,r-r2}],{r2,0,r}],{q4,0,q-q2-q3}],{q3,0,q-q2}],{q2,0,q}]],(#1[[1]][[1]]<#2[[1]][[1]])&,offset]];

iab[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=q-q2-q3-q4;
iac[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=r-r2-r3-r4;
iad[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=s-s2-s3-s4;
cab[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=q-q2-r2-s2;
cac[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=r-q3-r3-s3;
cad[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=s-q4-r4-s4;
p1[q2_,q3_,q4_,r2_,r3_,r4_,s2_,s3_,s4_,p_,q_,r_,s_]:=p-q+q2+q3+q4-r+r2+r3+r4-s+s2+s3+s4;



(* Hook   *)

HookTensor/:HookTensor[_Symbol,z_List,_List](f_)[{v___,x_,w___}] (f_)[{v___,y_,w___}]:=0/;(MemberQ[z,x]&&MemberQ[z,y]);
HookTensor/:HookTensor[_Symbol, z_List, a_List](_)[___,x_List,___]:=0/;Intersection[x,Join[z,a]]===Sort[Join[z,a]];
HookTensor/:HookTensor[_Symbol, z_List, a_List] Unevaluated[NonCommutativeMultiply[___,(_)[___,x_List,___],___]]:=0/;Intersection[x,Join[z,a]]===Sort[Join[z,a]];

HookTensor[sym_Symbol,x_List,y_List]:=Signature[x] Apply[HookTensor,{sym,Sort[x],y}]/; !OrderedQ[x];
HookTensor[sym_Symbol,{x_},{y_}]:=HookTensor[sym,{y},{x}]/;!OrderedQ[{x,y}];
HookTensor[_Symbol,{___,x_,x_,___},y_List]:=0;
HookTensor[_Symbol,{___,x_,___},{x_}]:=0;

HS[x_,l_List:{{}}]:=HookSimplify[x,l];

HookSimplify[x_,l_List:{{}}]:=RD[Distribute[HookSimplifyLastt[FixedPoint[Expand,RD[Distribute[HookSimplifyFirstt[FixedPoint[Expand,x],l],Plus,HookSimplifyFirstt,Plus,HookSimplifyFirst]]],l],Plus,HookSimplifyLastt,Plus,HookSimplifyLast]];
HookSimplifyFirst[object_Symbol[sym_Symbol,x_List,y_List]v_.,l_List]:=Module[{indexlists=AllIndexLists[v]},If[Length[indexlists]!=0&&MemberQ[First[indexlists],First[y]],Module[{ltotal=Join[l,Union[DeleteCases[indexlists,{_}]]],hookcomb=Subsets[Join[x,y],{1}],hookcompo,hooktotal,lrewritten},hookcompo=Map[Complement[Join[x,y],#]&,hookcomb];hooktotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[hookcompo[[i]],hookcomb[[i]]]],{i,Length[hookcomb]}]];hookleft=Select[hooktotal,!MemberQ[First[indexlists],Last[#]]&];lrewritten=Length[hooktotal]-Length[hookleft];-Signature[Join[x,y]]Apply[Plus,Map[Signature[#]Count[hooktotal,#]HookTensor[sym,Drop[#,-1],Take[#,-1]]v/lrewritten&,Union[hookleft]]]],HookTensor[sym,x,y] v]]/;MemberQ[{HookTensor,HookTensorSpinor},object];
HookSimplifyLast[object_Symbol[sym_Symbol,x_List,y_List]v_.,l_List]:=Module[{indexlists=AllIndexLists[v]},If[Length[indexlists]!=0&&MemberQ[Last[indexlists],First[y]],Module[{ltotal=Join[l,Union[DeleteCases[indexlists,{_}]]],hookcomb=Subsets[Join[x,y],{1}],hookcompo,hooktotal,lrewritten},hookcompo=Map[Complement[Join[x,y],#]&,hookcomb];hooktotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[hookcompo[[i]],hookcomb[[i]]]],{i,Length[hookcomb]}]];hookleft=Select[hooktotal,!MemberQ[Last[indexlists],Last[#]]&];lrewritten=Length[hooktotal]-Length[hookleft];-Signature[Join[x,y]]Apply[Plus,Map[Signature[#]Count[hooktotal,#]HookTensor[sym,Drop[#,-1],Take[#,-1]]v/lrewritten&,Union[hookleft]]]],HookTensor[sym,x,y] v]]/;MemberQ[{HookTensor,HookTensorSpinor},object];
HookSimplifyLast[Unevaluated[NonCommutativeMultiply[z___,HookTensorSpinor[sym_Symbol,x_List,y_List]]v_.],l_List]:=Module[{indexlists=AllIndexLists[v NonCommutativeMultiply[z]]},If[MemberQ[Last[indexlists],First[y]],Module[{ltotal=Join[l,Union[DeleteCases[indexlists,{_}]]],hookcomb=Subsets[Join[x,y],{1}],hookcompo,hooktotal,lrewritten},hookcompo=Map[Complement[Join[x,y],#]&,hookcomb];hooktotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[hookcompo[[i]],hookcomb[[i]]]],{i,Length[hookcomb]}]];hookleft=Select[hooktotal,!MemberQ[Last[indexlists],Last[#]]&];lrewritten=Length[hooktotal]-Length[hookleft];-Signature[Join[x,y]]Apply[Plus,Map[Signature[#]Count[hooktotal,#]NonCommutativeMultiply[z,HookTensorSpinor[sym,Drop[#,-1],Take[#,-1]]]v/lrewritten&,Union[hookleft]]]],NonCommutativeMultiply[z,HookTensorSpinor[sym,x,y]] v]];
HookSimplifyFirst[Unevaluated[NonCommutativeMultiply[z___,HookTensorSpinor[sym_Symbol,x_List,y_List]]v_.],l_List]:=Module[{indexlists=AllIndexLists[v NonCommutativeMultiply[z]]},If[MemberQ[First[indexlists],First[y]],Module[{ltotal=Join[l,Union[DeleteCases[indexlists,{_}]]],hookcomb=Subsets[Join[x,y],{1}],hookcompo,hooktotal,lrewritten},hookcompo=Map[Complement[Join[x,y],#]&,hookcomb];hooktotal=Map[ACOList[#,ltotal]&,Table[Flatten[Append[hookcompo[[i]],hookcomb[[i]]]],{i,Length[hookcomb]}]];hookleft=Select[hooktotal,!MemberQ[First[indexlists],Last[#]]&];lrewritten=Length[hooktotal]-Length[hookleft];-Signature[Join[x,y]]Apply[Plus,Map[Signature[#]Count[hooktotal,#]NonCommutativeMultiply[z,HookTensorSpinor[sym,Drop[#,-1],Take[#,-1]]]v/lrewritten&,Union[hookleft]]]],NonCommutativeMultiply[z,HookTensorSpinor[sym,x,y]] v]];
HookSimplifyFirst[Unevaluated[NonCommutativeMultiply[HookTensorSpinor[s_Symbol,x_List,y_List],z_]v_.],_List]:=v NonCommutativeMultiply[HookTensorSpinor[s,x,y],z];
HookSimplifyLast[Unevaluated[NonCommutativeMultiply[HookTensorSpinor[s_Symbol,x_List,y_List],z_]v_.],_List]:=v NonCommutativeMultiply[HookTensorSpinor[s,x,y],z];
HookSimplifyFirst[x_/;FreeQ[x,HookTensor]&&FreeQ[x,HookTensorSpinor],_List]:=x;
HookSimplifyLast[x_/;FreeQ[x,HookTensor]&&FreeQ[x,HookTensorSpinor],_List]:=x;

HookProject[expr_,x_List,a_List]:=Distribute[Hpt[FixedPoint[Expand,expr],x,a],Plus,Hpt,Plus,Hp];
HP[expr_,x_List,a_List]:=Distribute[Hpt[FixedPoint[Expand,expr],x,a],Plus,Hpt,Plus,Hp];

Hp[(v_Integer|v_Rational|v_Symbol|v_Complex) x_,tensor_List,hook_List]:=v Hp[x,tensor,hook];
Hp[(obj_Symbol)[sym_Symbol,tensor_List,hook_List],tensor_List,hook_List]:=obj[sym,tensor,hook]/;obj===HookTensor||obj===HookTensorSpinor;

Hp[expr_,x_List,a_List]:=Module[{d=Length[x],temp1=(ASym[expr,x]-ASym[expr,Join[x,a]])},RD[temp1-d/(Dim-d+1) ASym[Delta[{First[x]},a](temp1/.{First[x]:>dtemp,First[a]:>dtemp}),x]]];

RemoveHookGammaTrace[expr_,x_List,a_List]:=Distribute[Rhgtt[FixedPoint[Expand,expr],x,a],Plus,Rhgtt,Plus,Rhgt];
RHGT[expr_,x_List,a_List]:=Distribute[Rhgtt[FixedPoint[Expand,expr],x,a],Plus,Rhgtt,Plus,Rhgt];

Rhgt[(v_Integer|v_Rational|v_Symbol|v_Complex) x_,tensor_List,hook_List]:=v Rhgt[x,tensor,hook];
Rhgt[x_/;(FreeQ[x,TensorSpinor]&&FreeQ[x,HookTensorSpinor]),_List,_List]=x;

Rhgt[expr_,x_List,a_List]:=Module[{allind=AllIndices[expr], lentensorind=Length[x],tensorpos,hookpos,p,term1,term2,term3},tensorpos=PositionList[allind,x];hookpos=PositionList[allind,a];
term1=Sum[Module[{temp1=Take[tensorpos,p],temp3=Take[ppList,p]},k[lentensorind,p]GammaProd[Take[x,p],temp3]**Reconstruct[AllSize[expr],ReplacePart[allind,temp3,temp1,(Range[p]/.y_Integer:>List[y])]]],{p,1,lentensorind}];term2=Sum[Module[{temp1=Take[tensorpos,p],temp2=Drop[tensorpos,p],temp3=Take[ppList,p],temp4=Drop[x,p-1]},l[lentensorind,p]GammaProd[Join[a,Take[x,p-1]],temp3]**Reconstruct[AllSize[expr],ReplacePart[ReplacePart[allind,temp3,temp1,(Range[p]/.y_Integer:>List[y])],temp4,Join[temp2,hookpos],(Range[lentensorind-p+1]/.z_Integer:>List[z])]]],{p,1,lentensorind}];term3=Sum[Module[{temp1=Take[tensorpos,p],temp2=Drop[tensorpos,p],temp3=Take[ppList,p],temp4=Drop[x,p-1]},m[lentensorind,p] Delta[a,{First[x]}]GammaProd[Take[x,{2,p-1}],temp3]**Reconstruct[AllSize[expr],ReplacePart[ReplacePart[allind,temp3,temp1,(Range[p]/.y_Integer:>List[y])],temp4,Join[temp2,hookpos],(Range[lentensorind-p+1]/.z_Integer:>List[z])]]],{p,2,lentensorind}];expr+ASym[term1+term2+term3,x]-ASym[term1+term2+term3,Join[x,a]]
];

k[n_,p_]:=Power[-1,Evaluate[(p (p+1)/2)]] Binomial[n,p]/(p!Binomial[Dim-2n+p+1,p]);
l[n_,p_]:=Power[-1,Evaluate[n+1]] ((Dim-n+1)(Dim-2n+1)-(n+1)p)/((Dim+2)(Dim-n+2)) k[n,p];
m[n_,p_]:=Power[-1,Evaluate[n]] ((Dim-2n+p+1)(n+1)(p-1))/((Dim+2)(Dim-n+2)) k[n,p];



SetHookTensorInd[x_]:=Module[{},HookTensorIndicator=x;];

Format[HookTensor[sym_Symbol,x_List,y_List],{StandardForm,TeXForm}]:=DisplayForm[Subscript[HookTensorIndicator[sym],RowBox[x],RowBox[y]]];

Unprotect[Subscript];
Subscript[HookTensorIndicator[sym_Symbol/;(sym=!=Global`\[Delta]&&sym=!=Global`\[CapitalGamma])],x_Symbol/;FreeQ[x,RowBox],y_Symbol]:=HookTensor[sym,{x},{y}];
Subscript[HookTensorIndicator[sym_Symbol/;(sym=!=Global`\[Delta]&&sym=!=Global`\[CapitalGamma])],x_/;FreeQ[x,RowBox],y_Symbol]:=HookTensor[sym,Unevaluated[x]/.Times->List,{y}];
Protect[Subscript];


(* HookTensorSpinor *)

HookTensorSpinor/:HookTensorSpinor[_Symbol,z_List,_List](f_)[{v___,x_,w___}] (f_)[{v___,y_,w___}]:=0/;(MemberQ[z,x]&&MemberQ[z,y]);

HookTensorSpinor/:(y_+z_)**HookTensorSpinor[s_Symbol,x_List,a_List]:=y**HookTensorSpinor[s,x,a]+z**HookTensorSpinor[s,x,a];
HookTensorSpinor/:HookTensorSpinor[s_Symbol,x_List,a_List]**(y_+z_):=HookTensorSpinor[s,x,a]**y+HookTensorSpinor[s,x,a]**z;
HookTensorSpinor/:(y_/;(Map[FreeQ[y,#]&,NCList]/.List:>And)w_.)**HookTensorSpinor[s_Symbol,x_List,a_List]:=y (w**HookTensorSpinor[s,x,a]);
HookTensorSpinor/:HookTensorSpinor[s_Symbol,x_List,a_List]**(y_/;(Map[FreeQ[y,#]&,NCList]/.List:>And)w_.):=y (HookTensorSpinor[s,x,a]**w);
HookTensorSpinor/: GammaProd[___,x_List]**HookTensorSpinor[s_Symbol,y_List,a_List]:=0/;Intersection[x,y]===x;

HookTensorSpinor[sym_Symbol,x_List,y_List]:=Signature[x] Apply[HookTensorSpinor,{sym,Sort[x],y}]/; !OrderedQ[x];
HookTensorSpinor[sym_Symbol,{x_},{y_}]:=HookTensorSpinor[sym,{y},{x}]/;!OrderedQ[{x,y}];
HookTensorSpinor[_Symbol,{___,x_,___,x_,___},y_List]:=0;
HookTensorSpinor[_Symbol,{___,x_,___},{x_}]:=0;

SetHookTensorSpinorInd[x_]:=Module[{},HookTensorSpinorIndicator=x;];

Format[HookTensorSpinor[sym_Symbol,x_List,y_List],{StandardForm,TeXForm}]:=DisplayForm[Subscript[HookTensorSpinorIndicator[sym],RowBox[x],RowBox[y]]];

Unprotect[Subscript];
Subscript[HookTensorSpinorIndicator[sym_Symbol/;(sym=!=Global`\[Delta]&&sym=!=Global`\[CapitalGamma])],x_Symbol/;FreeQ[x,RowBox],y_Symbol]:=HookTensorSpinor[sym,{x},{y}];
Subscript[HookTensorSpinorIndicator[sym_Symbol/;(sym=!=Global`\[Delta]&&sym=!=Global`\[CapitalGamma])],x_/;FreeQ[x,RowBox],y_Symbol]:=HookTensorSpinor[sym,Unevaluated[x]/.Times->List,{y}];

Subscript[TensorSpinor[sym_Symbol,{}],x_Symbol/;FreeQ[x,RowBox],y_Symbol]:=HookTensorSpinor[sym,{x},{y}];
Subscript[TensorSpinor[sym_Symbol,{}],x_/;FreeQ[x,RowBox],y_Symbol]:=HookTensorSpinor[sym,Unevaluated[x]/.Times->List,{y}];
Protect[Subscript];


(* TensorSpinor & GammaExtract *)

GammaExtract[x_,y_List:{}]:=ACanonicalOrder[Distribute[GExtt[FixedPoint[Expand,x]],Plus,GExtt,Plus,GExt],y];
GET[x_,y_List:{}]:=ACanonicalOrder[Distribute[GExtt[FixedPoint[Expand,x]],Plus,GExtt,Plus,GExt],y];

GExt[x_/;(FreeQ[x,GammaProd]||(FreeQ[x,TensorSpinor]&&FreeQ[x,HookTensorSpinor]))]:=x;

GExt[Unevaluated[NonCommutativeMultiply[x___,GammaProd[y___,z_List],(obj_)[s_Symbol,t__List]]v_.]]:=v NonCommutativeMultiply[x,GammaProd[y,z],obj[s,t]]/;Intersection[z,Join[t]]==={};
GExt[Unevaluated[NonCommutativeMultiply[___,GammaProd[___List,{w_}],(_)[_Symbol,x__List]]v_.]]:=0/;MemberQ[Join[x],w];
GExt[Unevaluated[NonCommutativeMultiply[z___,GammaProd[y___List,x_List],(obj_)[s_Symbol,t__List]]v_.]]:=Module[
{inter,xxleft,n,r},inter=Intersection[x,Join[t]];
xxleft=Complement[x,inter];
n=Length[xxleft]; r=Length[inter];
If[n<r,0,Signature[Join[xxleft,inter]]Module[{xxcomb=Subsets[xxleft,{r}],actxxcompo}, v (g[n,r]/Length[xxcomb])Apply[Plus,Map[Signature[Join[#[[1]],actxxcompo=Complement[xxleft,#[[1]]]]]Delta[#[[1]],#[[2]]]NonCommutativeMultiply[z,GammaProd[y,actxxcompo],obj[s,t]]&,Flatten[Outer[List,xxcomb,{inter},1],1]]]]]];

G[n_,r_]:=G[n,r]=G1[n]G2[r]-Sum[Power[-1,Evaluate[(n p-p (p+1)/2)]]p! Binomial[n,p] Binomial[r,p]G[n-p,r-p],{p,1,Min[n,r]}];
g[n_,r_]:=g[n,r]=Coefficient[G[n,r],G1[n-r]G2[0]];

RemoveGammaTrace[x_,free_List:{}]:=Distribute[Rgtt[FixedPoint[Expand,x],free],Plus,Rgtt,Plus,Rgt];
RGT[x_,free_List:{}]:=Distribute[Rgtt[FixedPoint[Expand,x],free],Plus,Rgtt,Plus,Rgt];

Rgt[(v_Integer|v_Rational|v_Symbol|v_Complex) x_,free_List:{}]:=v Rgt[x,free];
Rgt[x_/;(FreeQ[x,TensorSpinor]&&FreeQ[x,HookTensorSpinor]),_List:{}]=x;

Rgt[x_,free_List:{}]:=Module[{allind=AllIndices[x],freeind, lenfreeind,pos,sortpos,p},If[free==={},freeind=Complement[allind,ReplaceList[allind,{___,t_,___,t_,___}:>t]];,freeind=free];lenfreeind=Length[freeind];pos=PositionList[allind,freeind];sortpos=Sort[pos];
	(AC[x,freeind]+Signature[pos]Sum[Module[{temp1=Take[sortpos,p],temp2=Drop[sortpos,p],temp3=Take[ppList,p],temp4=Drop[freeind,p],temp5},permtotal=Map[ACOList[#,{temp3,temp4}]&,Permutations[temp5=Join[temp3,temp4]]];Power[-1,Evaluate[(p (p+1)/2)]] Binomial[lenfreeind,p]/(p!lenfreeind!Binomial[Dim-2lenfreeind+p+1,p]) Signature[temp5]GammaProd[Take[freeind,p],temp3]**Apply[Plus,Map[Count[permtotal,#]Signature[#]Reconstruct[AllSize[x],If[p<lenfreeind,ReplacePart[ReplacePart[allind,Take[#,p],temp1,(Range[Length[temp1]]/.y_Integer:>List[y])],Drop[#,p],temp2,(Range[Length[temp2]]/.z_Integer:>List[z])],temp4={};ReplacePart[allind,ppList,sortpos,(Range[Length[sortpos]]/.y_Integer:>List[y])]]]&,Union[permtotal]]]],{p,1,lenfreeind}])];
(* Sufficient up to Dim=20 *)
ppList={Global`p1,Global`p2,Global`p3,Global`p4,Global`p5,Global`p6,Global`p7,Global`p8,Global`p9,Global`p10,Global`p11,Global`p12,Global`p13,Global`p14,Global`p15,Global`p16,Global`p17,Global`p18,Global`p19,Global`p20};

TensorSpinor/:Unevaluated[(y_+z_)**TensorSpinor[s_Symbol,x_List]]:=y**TensorSpinor[s,x]+z**TensorSpinor[s,x];
TensorSpinor/:Unevaluated[TensorSpinor[s_Symbol,x_List]**(y_+z_)]:=TensorSpinor[s,x]**y+TensorSpinor[s,x]**z;
TensorSpinor/:Unevaluated[(y_/;(Map[FreeQ[y,#]&,NCList]/.List:>And)w_.)**TensorSpinor[s_Symbol,x_List]]:=y (w**TensorSpinor[s,x]);
TensorSpinor/:Unevaluated[TensorSpinor[s_Symbol,x_List]**(y_/;(Map[FreeQ[y,#]&,NCList]/.List:>And)w_.)]:=y (TensorSpinor[s,x]**w);
(*TensorSpinor/: GammaProd[___,x_List]**TensorSpinor[s_Symbol,y_List]:=0/;Intersection[x,y]===x;*)

TensorSpinor[s_Symbol,x_List]:=Signature[x] Apply[TensorSpinor,{s,Sort[x]}]/; !OrderedQ[x];
TensorSpinor[_Symbol,{___,x_,___,x_,___}]:=0;

SetTensorSpinorInd[x_]:=Module[{},TensorSpinorIndicator=x;];

Format[TensorSpinor[s_Symbol,{}],{StandardForm,TeXForm}]:=DisplayForm[TensorSpinorIndicator[RowBox[{s}]]];
Format[TensorSpinor[s_Symbol,x_List],{StandardForm,TeXForm}]:=DisplayForm[Subscript[TensorSpinorIndicator[s],RowBox[x]]];

Unprotect[Times];
Times[x__,TensorSpinorIndicator[sym_Symbol/;(sym=!=\[Delta]&&sym=!=Global`\[CapitalGamma]&&FreeQ[sym,RowBox])]]:=Times[x, TensorSpinor[sym,{}]];
Protect[Times];

Unprotect[Subscript];
Subscript[TensorSpinor[sym_Symbol,{}],x_Symbol]:=TensorSpinor[sym,{x}];
Subscript[TensorSpinor[sym_Symbol,{}],x_/;FreeQ[x,RowBox]]:=TensorSpinor[sym,Unevaluated[x]/.Times->List];
Protect[Subscript];

Evaluate[TensorSpinorIndicator]/:Subscript[TensorSpinorIndicator[sym_Symbol/;(sym=!=\[Delta]&&sym=!=Global`\[CapitalGamma])],x_Symbol/;FreeQ[x,RowBox]]:=TensorSpinor[sym,{x}];Evaluate[TensorSpinorIndicator]/:Subscript[TensorSpinorIndicator[sym_Symbol/;(sym=!=\[Delta]&&sym=!=Global`\[CapitalGamma])],x_/;FreeQ[x,RowBox]]:=TensorSpinor[sym,Unevaluated[x]/.Times->List];
Evaluate[TensorSpinorIndicator][sym_Symbol/;(sym=!=\[Delta]&&sym=!=Global`\[CapitalGamma])]:=TensorSpinor[sym,{}];


(* Auxiliary functions *)

AllIndices[Power[(x_)[w___Symbol,v__List],n_]y_.]:=Join[Flatten[Table[v,{n}]],AllIndices[y]];
AllIndices[(x_)[System`OverBar[_],v__List]y_.]:=Join[v,AllIndices[y]];
AllIndices[(x_)[w___Symbol,v__List]y_.]:=Join[v,AllIndices[y]];
AllIndices[u_. HoldPattern[NonCommutativeMultiply[(x_)[w___Symbol,v__List],y__]]]:=Join[v,AllIndices[u NonCommutativeMultiply[y]]];
AllIndices[u_. HoldPattern[NonCommutativeMultiply[(x_)[w___Symbol,v__List]]]]:=Join[v,AllIndices[u]];
AllIndices[(_Integer|_Rational|_Symbol|_Complex)]:={};
(*AllIndices[Power[___]]:={};*)
(*AllIndices[x_]:={};*)

AllIndices[x_]:=Module[{},Message[AllIndices::arg];Print[x];Abort[]];

AllIndices::arg="Unsupported argument.";

AllIndexLists[x:(_Integer|_Rational|_Symbol|_Complex)v_]:=AllIndexLists[v];
AllIndexLists[Power[(x_)[w___Symbol,v__List],n_]y_.]:=Join[Table[v,{n}],AllIndexLists[y]];
AllIndexLists[(x_)[System`OverBar[_],v__List]y_.]:=Join[{v},AllIndexLists[y]];
AllIndexLists[(x_)[w___,v__List]y_.]:=Join[{v},AllIndexLists[y]];
AllIndexLists[u_. HoldPattern[NonCommutativeMultiply[(x_)[w___Symbol,v__List],y__]]]:=Join[{v},AllIndexLists[u NonCommutativeMultiply[y]]];
AllIndexLists[u_. HoldPattern[NonCommutativeMultiply[(x_)[w___Symbol,v__List]]]]:=Join[{v},AllIndexLists[u]];
AllIndexLists[(_Integer|_Rational|_Symbol|_Complex)]:={};

AllIndexLists[x_]:=Module[{},Message[AllIndexLists::arg];Print[x];Abort[]];

AllIndexLists::arg="Unsupported argument.";

AllSize[Power[(x_)[w___Symbol,v__List],n_]y_.]:=Join[Flatten[Table[{{x,{w},Map[Length,{v}]}},{n}],1],AllSize[y]];
AllSize[(x_)[System`OverBar[w_],v__List]y_.]:=Join[{{x,{System`OverBar[w]},Map[Length,{v}]}},AllSize[y]];
AllSize[(x_)[w___,v__List]y_.]:=Join[{{x,{w},Map[Length,{v}]}},AllSize[y]];
AllSize[u_. HoldPattern[NonCommutativeMultiply[(x_)[w___Symbol,v__List],y__]]]:=Join[{{x,{w},Map[Length,{v}]}},AllSize[u NonCommutativeMultiply[y]]];
AllSize[u_. HoldPattern[NonCommutativeMultiply[(x_)[w___Symbol,v__List]]]]:=Join[{{x,{w},Map[Length,{v}]}},AllSize[u]];
AllSize[(_Integer|_Rational|_Symbol|_Complex)]:={};
(*AllSize[Power[___]]:={};*)

AllSize[x_]:=Module[{},Message[AllSize::arg];Print[x];Abort[]];

AllSize::arg="Unsupported argument.";

Reconstruct[{struct_, rest___},ind_]:=If[struct[[2]]==={},Apply[struct[[1]],MyPartition[struct[[3]],ind]] **Reconstruct[{rest},Drop[ind,struct[[3]]/.List->Plus]],Apply[struct[[1]],Join[struct[[2]],MyPartition[struct[[3]],ind]]] **Reconstruct[{rest},Drop[ind,struct[[3]]/.List->Plus]]];
Reconstruct[{},{}]:=1;

MyPartition[{x_, y__},ind_]:=Prepend[MyPartition[{y},Drop[ind,x]],Take[ind,x]];
MyPartition[{x_},ind_]:={Take[ind,x]};

PositionList[lis_,find_]:=If[Intersection[lis,find]===Sort[find],Map[Position[lis,#][[1]]&,find],Message[PositionList::arg];Hold[PositionList[lis,find]];Print[{lis,find}]Abort[]];
PositionList::arg="An index sought for is not in the expression.";

HookCanonicalOrder[x_,u_List,v_List]:=Distribute[HCOt[FixedPoint[Expand,x],u,v],Plus,HCOt,Plus,HCO];
HC[x_,u_List,v_List]:=Distribute[HCOt[FixedPoint[Expand,x],u,v],Plus,HCOt,Plus,HCO];

HCO[HookTensor[s_Symbol,x_List,y_List]v_.,asymlist_List,hookindex_List]:=If[(Sort[Join[x,y]]===Sort[Join[asymlist,hookindex]])&&y=!=hookindex,-Signature[Join[x,y]]Signature[Join[asymlist,hookindex]]v HookTensor[s,asymlist,hookindex]/Length[asymlist],v HookTensor[s,x,y]];

HCO[HookTensorSpinor[s_Symbol,x_List,y_List]v_.,asymlist_List,hookindex_List]:=If[(Sort[Join[x,y]]===Sort[Join[asymlist,hookindex]])&&y=!=hookindex,-Signature[Join[x,y]]Signature[Join[asymlist,hookindex]]v HookTensor[s,asymlist,hookindex]/Length[asymlist],v HookTensor[s,x,y]];

ACanonicalOrder[x_,{y__Symbol}]:=AC[x,{{y}}];
ACanonicalOrder[x_,{{}}]:=x;
ACanonicalOrder[x_,v_List]:=Distribute[ACOt[FixedPoint[Expand,x],v],Plus,ACOt,Plus,FoldACO];
AC[x_,{y__Symbol}]:=AC[x,{{y}}];
AC[x_,{{}}]:=x;
AC[x_,v_List]:=Distribute[ACOt[FixedPoint[Expand,x],v],Plus,ACOt,Plus,FoldACO];

FoldACO[x_,v_List]:=Fold[ACO,x,v];

ACO[(x_Integer|x_Rational|x_Symbol|x_Complex),_List]:=x;
ACO[x_ y_./;FreeQ[x,_List],v_List]:=x ACO[y,v];

ACO[Delta[a1_List,a2_List] Delta[b1_List,b2_List] x_.,v_List]:=ACO[Delta[Join[a1,b1],Join[a2,b2]]x,v]/;(Complement[a1,v]==={}&&Complement[b1,v]==={});

ACO[Delta[x_List,y_List] z_.,v_List]:=Module[{lx,xxint,lenxxint,yyint,lenyyint},Which[lenxxint+lenyyint>lx,0,lenxxint+lenyyint==lx,Power[-1,Evaluate[((PositionList[x,xxint]/.List->Plus)-lenxxint (lenxxint+1)/2-(lx+1)lenyyint+(PositionList[y,yyint]/.List->Plus)+lenyyint (lenyyint+1)/2)]]lenxxint!lenyyint!/lx! ACO[Delta[Join[xxint,yyint],Join[Complement[y,yyint],Complement[x,xxint]]]z,v]]/;((lenxxint=Length[xxint=Intersection[x,v]])+(lenyyint=Length[yyint=Intersection[y,v]])>=(lx=Length[x])&&lenxxint!=0&&lenyyint!=0)];

ACO[(_)[x_Symbol/;x=!=Global`\[CapitalOmega],l_List,h_List]v_.,a_List]:=0/;Apply[And,Map[MemberQ[a,#]&,Join[l,h]]];

ACO[HoldPattern[NonCommutativeMultiply[___,(_)[_Symbol,l_List,h_List],___]v_.],a_List]:=0/;Apply[And,Map[MemberQ[a,#]&,Join[l,h]]];

ACO[expr_ ,sel_]:=Module[{ind=AllIndices[expr],poslist},
	Signature[poslist=PositionList[ind,sel]]Reconstruct[AllSize[expr],ReplacePart[ind,sel,Sort[poslist],(Range[Length[poslist]]/.x_Integer:>List[x])]]];

ACOList[l_List,order_List]:=Module[{retlist=l},
	Do[retlist=ACOSingleList[retlist,order[[i]]],{i,Length[order]}];retlist];
	
ACOSingleList[l_List,order_List]:=Module[{inter=Intersection[l,order],poslist},
	If[inter==={},
		l,
		poslist=PositionList[l,inter];
		ReplacePart[l,inter,Sort[poslist],(Range[Length[poslist]]/.x_Integer:>List[x])]
		]
	];

SCanonicalOrder[x_,{}]:=x;
SCanonicalOrder[x_,v_List]:=Distribute[SCOt[FixedPoint[Expand,x],v],Plus,SCOt,Plus,SCO];
SC[x_,{}]:=x;
SC[x_,v_List]:=Distribute[SCOt[FixedPoint[Expand,x],v],Plus,SCOt,Plus,SCO];

SCO[(x_Integer|x_Rational|x_Symbol|x_Complex),_List]:=x;
SCO[x_ y_./;FreeQ[x,_List],v_List]:=x SCO[y,v];

SCO[expr_,sel_List]:=0/;Max[Map[Length[Intersection[#,sel]]&,MyPartition[Flatten[Part[AllSize[expr],All,3]],AllIndices[expr]]]]>1;

SCO[expr_ ,sel_List]:=Module[{ind=AllIndices[expr],poslist},
poslist=PositionList[ind,sel];
Reconstruct[AllSize[expr],ReplacePart[ind,sel,Sort[poslist],(Range[Length[poslist]]/.x_Integer:>List[x])]]];

RenameDummy[x_,off_:0]:=RenameDummyLeft[RenameDummyRight[x,off],off];
RD[x_,off_:0]:=RenameDummyLeft[RenameDummyRight[x,off],off];
RDR[x_,off_:0]:=RenameDummyRight[x,off];
RDL[x_,off_:0]:=RenameDummyLeft[x,off];
RenameDummyLeft[x_,off_:0]:=Distribute[Rdt[FixedPoint[Expand,x],(#1[[1]][[1]]<#2[[1]][[1]])&,off],Plus,Rdt,Plus,Rd];
RenameDummyRight[x_,off_:0]:=Distribute[Rdt[FixedPoint[Expand,x],(#1[[2]][[1]]>#2[[2]][[1]])&,off],Plus,Rdt,Plus,Rd];

Rd[1,_,_:0]:=1;
Rd[x_  y_./;FreeQ[x,_List],test_,offset_:0]:=x Rd[y,test,offset];

Rd[expr_ ,test_,offset_:0]:=Module[{ind=AllIndices[expr],poslist},If[(poslist=Flatten[Sort[Map[Function[x,Position[ind,x]],ReplaceList[ind,{___,x_,___,x_,___}:>x]],test]])!={},Reconstruct[AllSize[expr],ReplacePart[ind,Drop[ddList,offset],(poslist/.x_Integer:>List[x]),(IntegerPart[(1+Range[Length[poslist]])/2]/.x_Integer:>List[x])]],expr ]];

ddList={Global`d0,Global`d1,Global`d2,Global`d3,Global`d4,Global`d5,Global`d6,Global`d7,Global`d8,Global`d9,Global`e0,Global`e1,Global`e2,Global`e3,Global`e4,Global`e5,Global`e6,Global`e7,Global`e8,Global`e9,Global`f0,Global`f1,Global`f2,Global`f3,Global`f4,Global`f5,Global`f6,Global`f7,Global`f8,Global`f9,Global`g0,Global`g1,Global`g2,Global`g3,Global`g4,Global`g5,Global`g6,Global`g7,Global`g8,Global`g9};

DummyCount[x_]:=Distribute[DCt[FixedPoint[Expand,x]],Plus,DCt,List,DC]

DC[1]:={0};
DC[x_  y_./;FreeQ[x,_List]]:=DummyCount[y];

DC[expr_ ]:=Module[{ind=AllIndices[expr]},{Length[Flatten[ReplaceList[ind,{___,x_,___,x_,___}:>x]]]}];

Multiply[x__]:=Distribute[Mut[Apply[Sequence,FixedPoint[Expand,{x}]]],Plus,Mut,Plus,Mu];
Mu[x_,y_,z__]:=Mu[Mu[x,y],z];
Mu[x_,y_]:=Module[{offset=Max[Flatten[DummyCount[x]]]},Rd[x,(#1[[1]][[1]]<#2[[1]][[1]])&]**Rd[y,(#1[[1]][[1]]<#2[[1]][[1]])&,offset]];
Mu[x_]:=x;


ASym[x_,{}]:=x;
ASym[x_,v_List]:=Distribute[ASt[FixedPoint[Expand,x],v],Plus,ASt,Plus,AS];

AS[0,_List]:=0;
AS[(x_Integer|x_Rational|x_Symbol|x_Complex)y_,v_List]:=x AS[y,v];
AS[Power[x1_,x2_]y_,v_List]:=Power[x1,x2] AS[y,v];

AS[Delta[x_List,y_List]z_.,w_List]:=0/;(Length[Intersection[x,w]]+Length[Intersection[y,w]]>Length[x]);
AS[Delta[a1_List,a2_List] Delta[b1_List,b2_List] x_.,v_List]:=AS[Delta[Join[a1,b1],Join[a2,b2]]x,v]/;(Complement[a1,v]==={}&&Complement[b1,v]==={});
AS[HookTensor[_Symbol,x_List,y_List] _.,z_List]:=0/;Intersection[Union[x,y],z]===Union[x,y];
AS[HookTensorSpinor[_Symbol,x_List,y_List] _.,z_List]:=0/;Intersection[Union[x,y],z]===Union[x,y];

AS[expr_, sel_]:=Module[{ind=AllIndices[expr],size=AllSize[expr],poslist,perms,len,i,posrange},
perms=Permutations[sel]; poslist=PositionList[ind,sel];
posrange=(Range[Length[poslist]]/.u_Integer:>List[u]);
1/ (len=Length[perms])Signature[sel]Sum[Signature[perms[[i]]]Reconstruct[size,ReplacePart[ind,perms[[i]],poslist,posrange]],{i,len}]];

Sym[x_,{}]:=x;
Sym[x_,v_List]:=Distribute[St[FixedPoint[Expand,x],v],Plus,St,Plus,S];

S[0,_List]:=0;
S[(x_Integer|x_Rational|x_Symbol|x_Complex)y_,v_List]:=x S[y,v];
S[Power[x1_,x2_]y_,v_List]:=Power[x1,x2] S[y,v];

S[expr_ , sel_]:=Module[{ind=AllIndices[expr],size=AllSize[expr],poslist,perms,len,i,posrange},
perms=Permutations[sel]; poslist=PositionList[ind,sel];
posrange=(Range[Length[poslist]]/.u_Integer:>List[u]);
1/ (len=Length[perms])Sum[Reconstruct[size,ReplacePart[ind,perms[[i]],poslist,posrange]],{i,len}]];

LinearPart[x_]:=Distribute[LPTt[FixedPoint[Expand,x]],Plus,LPTt,Plus,LPT];
LPT[x_]:=If[Count[Map[!FreeQ[x,#]&,FieldList],True]>1,0,x];

(*
Pick[expr_,obj_,len_]:=Apply[Plus,Flatten[Map [ReplaceList[#,x_/;\[Not]FreeQ[x,obj[_Symbol,l_/;Length[l]===len,___List]]\[Rule]x]&,Apply[List,expr]]]]
*)

Flip[x_]:=Distribute[Ft[FixedPoint[Expand,x]],Plus,Ft,Plus,F];

F[v_. NonCommutativeMultiply[TensorSpinor[s1_Symbol,x_List],GammaProd[y_List],TensorSpinor[s2_Symbol,z_List]]]:=-v Power[-1,(Length[y]+3)(Length[y]+2)/2] NonCommutativeMultiply[TensorSpinor[s2,z],GammaProd[y],TensorSpinor[s1,x]]/;x=!={};
F[v_. NonCommutativeMultiply[TensorSpinor[s1_Symbol,x_List],GammaProd[y_List],HookTensorSpinor[s2_Symbol,z_List,w_List]]]:=-v Power[-1,(Length[y]+3)(Length[y]+2)/2] NonCommutativeMultiply[HookTensorSpinor[s2,z,w],GammaProd[y],TensorSpinor[s1,x]]/;x=!={};
F[v_. NonCommutativeMultiply[HookTensorSpinor[s1_Symbol,x_List,w_List],GammaProd[y_List],TensorSpinor[s2_Symbol,z_List]]]:=-v Power[-1,(Length[y]+3)(Length[y]+2)/2] NonCommutativeMultiply[TensorSpinor[s2,z],GammaProd[y],HookTensorSpinor[s1,x,w]];
F[v_. NonCommutativeMultiply[HookTensorSpinor[s1_Symbol,x_List,u_List],GammaProd[y_List],HookTensorSpinor[s2_Symbol,z_List,w_List]]]:=-v Power[-1,(Length[y]+3)(Length[y]+2)/2] NonCommutativeMultiply[HookTensorSpinor[s2,z,w],GammaProd[y],HookTensorSpinor[s1,x,u]];
F[x_]:=x;

(* Complex conjugation *)

ComplexConjugate[expr_]:=CC[expr];
CC[expr_]:=Distribute[cct[FixedPoint[Expand,expr]],Plus,cct,Plus,cc];

cc[0]:=0;
cc[1]:=1;
cc[u_. Delta[x_List,y_List]]:=cc[u]Delta[x,y];

(* The Gamma matrices have real indices *)
cc[u_. GammaProd[x__]]:=cc[u]GammaProd[x]

cc[u_. Tensor[s_Symbol,x__List]]:=cc[u]Tensor[System`OverBar[s],x]/;MemberQ[ComplexFluxes,s];
cc[u_. Tensor[System`OverBar[s_Symbol],x__List]]:=cc[u]Tensor[s,x]/;MemberQ[ComplexFluxes,s];

cc[u_. Tensor[s_Symbol[t__],x__List]]:=cc[u]Tensor[System`OverBar[s][t],x]/;MemberQ[ComplexFluxes,s];
cc[u_. Tensor[System`OverBar[s_Symbol][t__],x__List]]:=cc[u]Tensor[s[t],x]/;MemberQ[ComplexFluxes,s];

cc[u_. Tensor[s_,x__List]]:=cc[u]Tensor[s,x];

cc[x:(_Integer|_Rational|_Complex)]:=(x/.Complex[z_,y_]->Complex[z,-y]);

cc[u_. Global`conj[s_Symbol]]:=cc[u] s
cc[u_. s_Symbol]:=cc[u] Global`conj[s]
cc[u_. System`OverBar[s_Symbol]]:=cc[u] s

cc[u_./s_Symbol]:=cc[u]/ Global`conj[s]
cc[u_. /Global`conj[s_Symbol]]:=cc[u]/ s

cc[u_. s_Symbol[x_]]:=cc[u]CC[s][x];

cc[u_. Power[x_,y_]]:=cc[u] Power[CC[x],y]

cc[u_. Complex[z_,y_]]:=cc[u]Complex[z,-y]

cc[x_]:=Module[{},Message[ComplexConjugate::arg];Print[x];Abort[]];

ComplexConjugate::arg="Unsupported argument.";

(* Formulae implementing the C* operator *)
(* (Action on HookTensorSpinor not implemented since never used.) *)

CStar[expr_]:=CS[expr];
CS[expr_]:=Distribute[cst[FixedPoint[Expand,expr]],Plus,cst,Plus,cs];

cs[0]:=0;
cs[1]:=1;

cs[u_. HoldPattern[NonCommutativeMultiply[x__]]]:=cc[u]NonCommutativeMultiply[Sequence@@Map[CS,{x}]]

cs[u_. GammaProd[x__List]]:=cc[u] GammaProd[x]

cs[u_. TensorSpinor[s_ /; ! StringMatchQ[ToString[s], "C" ~~ __], x_List]] := 
	cs[u] TensorSpinor[ToExpression["C" ~~ ToString[s]], x]
 
cs[u_. TensorSpinor[s_ /; StringMatchQ[ToString[s], "C" ~~ __], x_List]] := 
	cs[u] TensorSpinor[ToExpression[StringDrop[ToString[s], 1]], x]
	
cs[u_. Tensor[s_,x__List]]:=cs[u]cc[Tensor[s,x]];

cs[x:(_Integer|_Rational|_Complex)]:=cc[x];

cs[x_]:=Module[{},Message[CStar::arg];Print[x];Abort[]];

CStar::arg="Unsupported argument.";

(* Fierz *)

SetAttributes[Fierz, HoldFirst];

Fierz[indepvars_List,struct_List,asymlist_List,sym_,type_]:=Module[{indepfierz,deplist,eqlist},
indepfierz=Transpose[FierzTerm[indepvars,asymlist,sym,type]];
deplist=Delete[struct,PositionList[struct,indepvars]];{Length[deplist],Length[indepfierz],eqlist=MapIndexed[Function[{argument,pos},FierzTerm[{argument},asymlist,sym,type]-Sum[Global`a[pos[[1]],k]indepfierz[[k]],{k,Length[indepfierz]}]],deplist]}
(*{MatrixForm[Table[a[i,j],{i,Length[deplist]},{j,Length[indepfierz]}]]/.Flatten[SolveAlways[Map[#\[Equal]0&,eqlist/.Delta[__]\[RuleDelayed]1/;simp],vars=VarList[eqlist/.Delta[__]\[RuleDelayed]1/;simp]]],Length[vars]}*)];
FierzTerm[struct_List,asymlist_List,sym_,type_]:=Module[{expr},expr=Map[Function[{bas},Map[RenameDummy[GammaExtract[GammaExtract[Expand[ACanonicalOrder[ACanonicalOrder[GammaExpand[Expand[ACanonicalOrder[ACanonicalOrder[GammaExpand[ACanonicalOrder[ACanonicalOrder[GammaContract[RenameDummy[ACanonicalOrder[ACanonicalOrder[GammaExpand[type[#,GammaProd[bas]]/Length[Intersection[#[[1]],#[[2]]]]!],bas],asymlist]]],asymlist],bas]],bas],asymlist]]],bas],asymlist]],bas],asymlist]]&,struct]],Basis[sym]]
/.TensorSpinor[__]->1]/;(sym===Global`A||sym===Global`S);

FierzSolve[eqlist_List]:=MatrixForm[Table[Global`a[i,j],{i,eqlist[[1]]},{j,eqlist[[2]]}]/.Flatten[Map[SolveAlways[{#[[1]]==0,#[[2]]==0,#[[3]]==0}/.Delta[___]->0,VarList[#]]&,eqlist[[3]]/.Delta[___]->0]]];

MakeTensor[s_,basis_]:=ReleaseHold[(Delta[s[[4]][[1]],s[[4]][[2]]]s[[3]]/.{Global`P:>GammaProd[s[[1]]],Global`Q:>GammaProd[s[[2]]],Global`X:>basis,Global`tr:>GammaTrace})];

MakeHook[s_,basis_]:=Module[{n=Length[s[[3]]]},ACanonicalOrder[ReleaseHold[(s[[6]] Delta[s[[7]][[1]],s[[7]][[2]]]/.{Global`P->GammaProd[Join[s[[1]],Take[s[[3]],s[[5]]]]],Global`Q->GammaProd[Join[s[[2]],Drop[s[[3]],s[[5]]],s[[4]]]],Global`X:>basis,Global`tr:>GammaTrace})-Power[-1,Evaluate[n ]] (s[[6]] Delta[s[[7]][[1]],s[[7]][[2]]]/.{Global`P->GammaProd[Join[s[[1]],s[[4]],Take[s[[3]],s[[5]]-1]]],Global`Q->GammaProd[Join[s[[2]],Drop[s[[3]],s[[5]]-1]]],Global`X:>basis,Global`tr:>GammaTrace})+(Power[-1,Evaluate[(n+1+s[[5]])]](n+1)/(Dim-n+1)) (s[[6]] Delta[s[[7]][[1]],s[[7]][[2]]]Delta[s[[4]],{Last[s[[3]]]}]/.{Global`P->GammaProd[Join[s[[1]],Take[s[[3]],s[[5]]-1],{Global`d}]],Global`Q->GammaProd[Join[s[[2]],Drop[Drop[s[[3]],s[[5]]-1],-1],{Global`d}]],Global`X:>basis,Global`tr:>GammaTrace})],s[[3]]]];

Basis[Global`A]={{},{Global`h1,Global`h2,Global`h3},{Global`h1,Global`h2,Global`h3,Global`h4}};
Basis[Global`S]={{Global`h1},{Global`h1,Global`h2},{Global`h1,Global`h2,Global`h3,Global`h4,Global`h5}};
VarList[x_ ]:=Union[Flatten[Map[DeleteCases[Variables[#],Global`a[_Integer,_Integer]]&,x]]];


End[];
EndPackage[];



