Setup for (DDA)2DDDAeP2

leG=10;(* rank of gamma matrix *)
rank={1,4,1,4,1,1,4,2,leG};(* rank of (((DF))^2)(Subscript[DDFϵψ, 2])γ *)
pm=Permutations[{1,2}];(* symmetry *)
pm=Table[Flatten[Table[{2pm[[u,i]]-1,2pm[[u,i]]},{i,1,2}],1],{u,1,Length[pm]}];
lfn=(18+leG)/2;
fn=Take[{f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14},lfn];(* index list *)
ind={{0},{0,0,0,0},{0},{0,0,0,0},{0},{0},{0,0,0,0},{0,0},Table[0,{i,1,leG}]};
Tns={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11};(* (((DF))^2)(Subscript[DDFeψ, 2])γ *)
$Assumptions=Table[Tns[[2u+1]]∈Arrays[Table[11,{i,1,18+2u}],Reals,{Symmetric[{1,2}],Antisymmetric[{3,4,5}],Symmetric[{6,7}],Antisymmetric[{8,9,10}],Symmetric[{11,12,13}],Antisymmetric[{14,15,16}],Antisymmetric[{17,18}],Antisymmetric[Table[18+i,{i,1,2u}]],{Join[{6,7,8,9,10,1,2,3,4,5},Table[i,{i,11,18+2u}]],1}}],{u,0,5}];(* 2u = rank of gamma matrix *)
tns={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}[[1+leG]];(* choose tensor *)
cmp1=NestList[RotateLeft,{2,3,4,5},3];(* F=dA in 2nd position *)
cmp2=NestList[RotateLeft,{7,8,9,10},3];(* F=dA in 4th position *)
cmp3=NestList[RotateLeft,{13,14,15,16},3];(* F=dA in 7th position *)
rotl=Flatten[Table[{(((-1))^(u1+u2+u3+1)),{2->cmp1[[u1,1]],3->cmp1[[u1,2]],4->cmp1[[u1,3]],5->cmp1[[u1,4]],7->cmp2[[u2,1]],8->cmp2[[u2,2]],9->cmp2[[u2,3]],10->cmp2[[u2,4]],13->cmp3[[u3,1]],14->cmp3[[u3,2]],15->cmp3[[u3,3]],16->cmp3[[u3,4]]}},{u1,1,4},{u2,1,4},{u3,1,4}],2];(* use to construct F from dA *)

Setup for (DF)2DDFeP2

leG=6;(* rank of gamma matrix *)
rank={1,4,1,4,1,1,4,2,leG};(* rank of (((DF))^2)(Subscript[DDFϵψ, 2])γ *)
pm=Permutations[{1,2}];(* symmetry *)
pm=Table[Flatten[Table[{2pm[[u,i]]-1,2pm[[u,i]]},{i,1,2}],1],{u,1,Length[pm]}];
lfn=(18+leG)/2;
fn=Take[{f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14},lfn];
ind={{0},{0,0,0,0},{0},{0,0,0,0},{0},{0},{0,0,0,0},{0,0},Table[0,{i,1,leG}]};
Tns={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11};(* (((DF))^2)(Subscript[DDFϵγψ, 2]) *)
$Assumptions=Table[Tns[[2u+1]]∈Arrays[Table[11,{i,1,18+2u}],Reals,{Antisymmetric[{2,3,4,5}],Antisymmetric[{7,8,9,10}],Antisymmetric[{13,14,15,16}],Antisymmetric[{17,18}],Antisymmetric[Table[18+i,{i,1,2u}]],{Join[{6,7,8,9,10,1,2,3,4,5},Table[i,{i,11,18+2u}]],1}}],{u,0,5}];(* 2u = rank of gamma matrix *)
tns={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}[[1+leG]];(* choose tensor *)
rotl={{1,{}}};
psl={{5,{1,5}},{5,{6,10}},{2,{11,12}},{5,{12,16}}};(* {# for Bianchi, position of DF,DF,DD,DF} *)
sym={5,5};(* partition of (DF)2 *)
symrep={{1->1,2->2,3->3,4->4},{1->3,2->4,3->1,4->2}};
asym={};
pseoml={{5,{1,5}},{5,{6,10}},{5,{12,16}},{2,{17,18+leG}}};(* {# for eom, position of DF,DF,DF,γP2: 5 for DF, 2 for γP2, 3 for γDP2} *)

Adjacency matirx

(* generate adjacency matrix for p0th row *)
Genadjmat /: Genadjmat[matl0_,rank0_,p0_]:=
Module[{len,nb,t1,max,int,ptl,mat,matl},

len=Length[rank0];(* # of terms *)
nb=0;
Do[t1=matl0[[u]];
     max=Table[rank0[[j]]-Sum[t1[[k,j]],{k,1,p0-1}],{j,p0+1,len}];(* max value of j-th column *)
     int=rank0[[p0]]-Sum[t1[[p0,k]],{k,1,p0-1}];(* integer to be divided *)
     ptl=IntegerPartitions[int,{len-p0},Table[k,{k,0,Max[max]}]];(* partition of int into len-p0 *)
     ptl=Flatten[Table[Permutations[ptl[[v]]],{v,1,Length[ptl]}],1];(* duplicates are removed automatically *)
     Do[If[0≤Min[max-ptl[[v]]],
            Do[t1[[p0,p0+k]]=ptl[[v,k]];t1[[p0+k,p0]]=ptl[[v,k]],{k,1,len-p0}];
            mat[++nb]=t1;
     ];
   ,{v,1,Length[ptl]}];
,{u,1,Length[matl0]}];

matl=Table[mat[u],{u,1,nb}];
matl
];

matl=Genadjmat[{Table[0,{i,1,9},{j,1,9}]},{1,4,1,4,1,1,4,2,6},1]

(* normal ordered adjacency matrix by pm0 *)
Noradjmat /: Noradjmat[matl0_,pm0_]:=
Module[{len,matl,t1,t2,t3,chk},

len=Length[matl0[[1,1]]];(* # of terms *)
matl=matl0;
Do[t1=matl0[[u]];(* matrix which is ordered 1st *)
     Do[t2=matl0[[u]];(* reset t2 *)
          Do[t2[[pm0[[1,i]]]]=matl0[[u,pm0[[v,i]]]];,{i,1,Length[pm0[[1]]]}];(* exchange rows *)
          t3=t2;
          Do[t3[[All,pm0[[1,i]]]]=t2[[All,pm0[[v,i]]]];,{i,1,Length[pm0[[1]]]}];(* exchange columns *)
          chk=SelectFirst[Flatten[t1-t3],#=!=0&];(* check order *)
          If[chk<0,t1=t3;];
   ,{v,2,Length[pm0]}];
      matl[[u]]=t1;
,{u,1,Length[matl0]}];

Union[matl]
];

Noradjmat[matl,pm]

(* generate adjacency matrix *)
Adjmat /: Adjmat[rank0_,pm0_]:=
Module[{len,matl,nb,p0,t1,max,int,ptl,mv,mat},

len=Length[rank0];(* # of terms *)
matl={Table[0,{i,1,len},{j,1,len}]};(* initial adjacency matrix *)
Do[matl=Genadjmat[matl,rank0,u];,{u,1,len-4}];(* generate u-th row *)

nb=0;
p0=len-3;
Do[t1=matl[[u]];(* generate p0-th ~ len-th rows *)
     max=Table[rank0[[j]]-Sum[t1[[k,j]],{k,1,p0-1}],{j,p0+1,len}];
     int=rank0[[p0]]-Sum[t1[[p0,k]],{k,1,p0-1}];(* integer to be partitioned *)
     ptl=IntegerPartitions[int,{len-p0},Table[k,{k,0,int}]];
     ptl=Flatten[Table[Permutations[ptl[[v]]],{v,1,Length[ptl]}],1];(* partition list *)
     Do[If[0≤Min[max-ptl[[v]]],
            Do[t1[[p0,p0+k]]=ptl[[v,k]];t1[[p0+k,p0]]=ptl[[v,k]],{k,1,len-p0}];
            Do[mv[k]=rank0[[p0+k]]-Sum[t1[[p0+k,j]],{j,1,p0}];,{k,1,3}];
            t1[[p0+1,p0+2]]=(mv[1]+mv[2]-mv[3])/2;t1[[p0+2,p0+1]]=t1[[p0+1,p0+2]];
            t1[[p0+1,p0+3]]=(mv[1]-mv[2]+mv[3])/2;t1[[p0+3,p0+1]]=t1[[p0+1,p0+3]];
            t1[[p0+2,p0+3]]=(-mv[1]+mv[2]+mv[3])/2;t1[[p0+3,p0+2]]=t1[[p0+2,p0+3]];
            If[(0≤Min[t1[[p0+1,p0+2]],t1[[p0+1,p0+3]],t1[[p0+2,p0+3]]]),
                 mat[++nb]=t1(* take positive values *)
          ];
     ];
   ,{v,1,Length[ptl]}];
,{u,1,Length[matl]}];

matl=Table[mat[u],{u,1,nb}];
matl=Noradjmat[matl,pm0];
matl
];

adjl=Adjmat[rank,pm]

Print[Length[adjl]];adjl[[1]] // MatrixForm

Base list

(* transform adjacency matrix to tensor *)
MatInd /: MatInd[adjmat0_,fn0_,ind0_]:=
Module[{len,nb,ind,id},

len=Length[ind0];(* # of terms *)
nb=Table[0,{i,1,len}];(* number of indices inserted *)
ind=ind0;
id=1;
Do[Do[ind[[i,nb[[i]]+u]]=fn0[[id]];
         ind[[j,nb[[j]]+u]]=fn0[[id]];
         ++id;
   ,{u,1,adjmat0[[i,j]]}];
     nb[[i]]+=adjmat0[[i,j]];nb[[j]]+=adjmat0[[i,j]];
,{i,1,len-1},{j,i+1,len}];

ind=Flatten[ind];
ind
];

adjl[[1]] // MatrixForm

MatInd[adjl[[1]],fn,ind]

(* transform tensor to canonical form *)
CanForm /: CanForm[indf0_,fn0_,tns0_,rotl0_]:=
Module[{lfn,indpair,tc,hc,indeven,nb,cf},

lfn=Length[fn0];
indpair=Table[Position[indf0,fn0[[i]]],{i,1,Length[fn0]}];
indpair=Map[Flatten,indpair];

tc=1+Sum[rotl0[[u,1]]TensorReduce[TensorContract[tns0,indpair /. rotl0[[u,2]]]]/. {TensorContract[T1_,x1_]->H[x1]},{u,1,Length[rotl0]}];(* tc = 1+5H[x1]+7H[x2]+... *)

hc=0;
Do[cf=tc[[u]] /. {H[x_]->1};(* Coefficient of 5H[x1] : 5 *)
     indpair=Flatten[tc[[u]]/cf /. H->List,1];(* x1={{1,15},...} *)
     indpair=Flatten[Sort[Map[RotateLeft,indpair]]];(* exchange 2nd with 1st in x1 and sort *)
     indeven=Take[indpair,{2,2 lfn,2}];(* 2nd part of exchanged x1 *)
     nb=FromDigits[indeven-1,2lfn-1];(* 2lfn-2 = max# in (indeven-1) *)
     hc=hc+cf*H[nb];
,{u,2,Length[tc]}];

hc
];

tmp1=MatInd[adjl[[100]],fn,ind];Print[tmp1];
CanForm[tmp1,fn,tns,rotl]

(* indlist, canlist, hlist *)
Canlist /: Canlist[matl0_,tns0_,rotl0_]:=
Module[{rank,lfn,fn,ind,indl,canl,hlist,ps,lenh,hrdl,dig},

rank=Total[matl0[[1]]];(* rank *)
lfn=Total[rank]/2;(* length of fn *)
fn=Take[{f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22},lfn];(* fn *)
ind=TakeList[Table[0,{u,1,2lfn}],rank];(* index filled by 0 *)

(* generate indl and canl *)
indl=Table[MatInd[matl0[[u]],fn,ind],{u,1,Length[matl0]}];(* list of index *)
canl=Table[0,{u,1,Length[indl]}];
Do[canl[[u]]=CanForm[indl[[u]],fn,tns0,rotl0];,{u,1,Length[indl]}];(* list of canonical form *)
Do[If[canl[[u]]==0,indl[[u]]=Hold[Nothing];canl[[u]]=Hold[Nothing]];,{u,1,Length[matl0]}];
indl=ReleaseHold[indl];(* list of nonzero index *)
canl=ReleaseHold[canl];(* list of nonzero canonical form *)

hlist=Table[{},{u,1,Length[canl]}];
Do[ps=Position[canl[[u]],H[__]];
     hlist[[u]]=Extract[canl[[u]],ps];
,{u,1,Length[canl]}];
hlist=Union[Flatten[hlist]];
hlist=Table[hlist[[u,1]],{u,1,Length[hlist]}];(* list of #s in H[#] *)
lenh=Length[hlist];

Do[hrdl=Map[Mod[#,(10^i)]&,hlist];
     If[Length[Union[hrdl]]==lenh,hlist=hrdl;dig=i;Break[]];
,{i,IntegerLength[lenh]+1,Max[IntegerLength[hlist]]}];(* shorten #s in H[#] *)
hlist=Sort[hlist];(* sort *)
canl=canl /. {H[n_]->H[Mod[n,(10^dig)]]};(* shorten *)

{indl,canl,hlist}
];

{indl,canl,hlist}=Canlist[adjl,tns,rotl]

(* base list *)
Baselist /: Baselist[indl0_,canl0_,hlist0_,V0_]:=
Module[{mat2,indl,canl,hlist,dig,len,lenh,sum,eq,sol,psi,eq2,psh,eq3,chl,mat,hind,rep},

{indl,canl,hlist}={indl0,canl0,hlist0};
dig=Max[IntegerLength[hlist]];(* digit of hlist *)
len=Length[indl0];
lenh=Length[hlist0];

(* reduce indl and canl via Bianchi identity *)
sum=Sum[Expand[Ci[j]canl[[j]]],{j,1,len}];(* sum=Ci[1288]H[5479]+... *)
eq=Normal[CoefficientArrays[sum,Map[H,hlist]][[2]]];(* Coefficient of H[] *)
sol=Solve[eq==0];sol=sol[[1]];(* # solution = # independent index *)
psi=Table[{sol[[u,1,1]]},{u,1,Length[sol]}];(* positions of independent index *)
indl=Extract[indl,psi];(* base of index *)
canl=Extract[canl,psi];(* base of canonical form *)
len=Length[indl];(* redefine *)

(* eq is reevaluated since indl is reduced *)
sum=Sum[Expand[Ci[j]canl[[j]]],{j,1,len}];(* sum=-Ci[906]H[37604282243165]+... *)
eq=Normal[CoefficientArrays[sum,Map[H,hlist]][[2]]];(* Coefficient of H[] *)

(* reduce hlist *)
sum=Sum[Expand[Ch[i]eq[[i]]],{i,1,lenh}];(* sum=4Ch[69251]Ci[1]+... *)
eq2=Normal[CoefficientArrays[sum,Table[Ci[j],{j,1,len}]][[2]]];(* Coefficient of Ci[] *)
sol=Solve[eq2==0];sol=sol[[1]];(* # solution = # independent vector *)
psh=Table[{sol[[u,1,1]]},{u,1,Length[sol]}];(* position of independent vector *)
hlist=Extract[hlist,psh];
lenh=Length[hlist];

(* (Ci^T).I=(Ch^T).H *)
eq3=Extract[eq,psh];(* eqs for Ci[u] *)
chl=Table[Ch[u],{u,1,len}];
sol=Solve[eq3==chl,Table[Ci[u],{u,1,len}]];sol=sol[[1]];(* Ci=M.Ch *)
sol=Table[Expand[sol[[u,2]]],{u,1,len}];
mat=CoefficientArrays[sol,chl][[2]];(* M *)

hind=Table[V0[u],{u,1,len}].mat;(* (H^T)=(I^T).M *)
rep=Table[H[hlist[[u]]]->hind[[u]],{u,1,len}];(* {H[]->V0[]+...,...} *)

{indl,dig,rep}
];

basel=Baselist[indl,canl,hlist,V1]

(* classification of index via base list *)
ClassifyInd /: ClassifyInd[basel0_,indf0_,fn0_,tns0_,rotl0_]:=
Module[{indl,dig,rep,len,can},

{indl,dig,rep}=basel0;
len=Length[indl];(* lenh = Length[indl] *)
can=CanForm[indf0,fn0,tns0,rotl0] /.H[n_]->H[Mod[n,(10^dig)]] /.rep /. {H[_]->0};(* canonical form *)
can=Expand[can];

can
];

tmp1=basel[[1,1]];Print[tmp1];
ClassifyInd[basel,tmp1,fn,tns,rotl]

Print[Now];Do[tmp=ClassifyInd[basel,basel[[1,u]],fn,tns,rotl];Print[tmp];,{u,1,100}];Print[Now];

Bianchi identity

(* generate Bianchi identity for R,DF,DDF,DP2; psl0[[i,1]]={3,5,2,3} for R,DF,DD,DP2 *)
Bianchi /: Bianchi[basel0_,indf0_,rank0_,psl0_,fn0_,tns0_,rotl0_]:=
Module[{len,id,tm,ind,jg1,id2,ps,jg2,bidl,pm,bid,bid2,ind2,rep,sgn},

len=Length[psl0];
id=Table[{psl0[[i,1]],psl0[[i,2,2]]-psl0[[i,2,1]]+1},{i,1,len}];(* id={{3,4},{5,5},{2,2},{3,3}} for R,DF,DD,DP2 *)
tm=Table[Take[indf0,psl0[[i,2]]],{i,1,len}];(* R,DF,DD,DP2 of indf0 *)
ind=TakeList[indf0,rank0];(* index *)
jg1=Table[Apply[Or,Table[Length[Intersection[tm[[i]],ind[[u]]]]>=id[[i,1]],{u,1,Length[ind]}]],{i,1,len}];(* True if self Bianchi *)
id2=Table[If[jg1[[i]],id[[i,2]],0],{i,1,len}];(* extract id[[i,2]] for self Bianchi *)
ps=0;If[Max[id2]=!=0,ps=Position[id2,Max[id2]][[1,1]]];(* extract position of max id2, if max=!=0 *)
jg2=Table[Length[Union[tm[[i]]]]==id[[i,2]],{i,1,len}];(* True if fn's are different *)
Do[If[id[[i]]=={3,4},tm[[i]]=Take[tm[[i]],3]];,{i,1,len}];(* True for Riemann tensor *)

bidl={};
bid=ClassifyInd[basel0,indf0,fn0,tns0,rotl0];
If[ps=!=0,(* True if self Bianchi *)
     Which[id[[ps]]=={2,2},bid=bid-VBD[bid[[1]]],
                 id[[ps]]=={3,3},bid=bid-VBP[bid[[1]]];
   ];
     bidl=Join[bidl,{bid}];,
     Do[If[jg2[[i]],
          bid2=bid;
          Which[id[[i]]=={2,2},bid2=bid2-2VBD[bid[[1]]],
                      id[[i]]=={3,3},bid2=bid2-3VBP[bid[[1]]];
        ];
          pm=Table[RotateLeft[tm[[i]],u],{u,1,id[[i,1]]-1}];(* rotation for Bianchi id *)
          Do[rep=Table[psl0[[i,2,1]]+v-1->pm[[u,v]],{v,1,Length[pm[[u]]]}];
               sgn=If[OddQ[id[[i,1]]],1,(((-1))^u)];
               ind2=ReplacePart[indf0,rep];
               bid2=bid2+sgn ClassifyInd[basel0,ind2,fn0,tns0,rotl0];(* generate Bianchi identity *)
        ,{u,1,id[[i,1]]-1}];
          bidl=Join[bidl,{bid2}];
   ];
,{i,1,len}];
];

bidl=Union[bidl];
bidl
];

basel=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base6";

tmp1=basel[[1,40]];Print[tmp1];
Bianchi[basel,tmp1,rank,psl,fn,tns,rotl]

(* generate list of Bianchi identity for R,DF,DD,DP2 *)
Bianchilist /: Bianchilist[basel0_,rank0_,psl0_,fn0_,tns0_,rotl0_,V0_]:=
Module[{indfl,bidl,sum,eq,sol,ps},

indfl=basel0[[1]];
bidl=Table[Bianchi[basel0,indfl[[u]],rank0,psl0,fn0,tns0,rotl0],{u,1,Length[indfl]}];
bidl=Union[Flatten[bidl]];
sum=Sum[Expand[Cb[u]bidl[[u]]],{u,1,Length[bidl]}] /. {VBD[_]->0,VBP[_]->0};

eq=Normal[CoefficientArrays[sum,Table[V0[j],{j,1,Length[indfl]}]][[2]]];(* Coefficient of tns0 *)
sol=Solve[eq==0];sol=sol[[1]];
ps=Table[{sol[[u,1,1]]},{u,1,Length[sol]}];
bidl=Extract[bidl,ps];

bidl
];

bidl=Bianchilist[basel,rank,psl,fn,tns,rotl,V1]

bv1[6]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2bid6";(* leg=6, 9838 *)

Dimension dependent identity

generate multiplicity list to construct num list

(* generate multiplicity list *)
MpNumlistSet /: MpNumlistSet[lfn0_,rank0_]:=
Module[{ptl,mpl,id,nb,num,mpnuml},

ptl=2*IntegerPartitions[lfn0,{11},Table[i,{i,1,IntegerPart[Length[rank0]/2]}]];(* partition: {#0,…,#10}={4,2,2,2,2,2,2,2,2,2,2}, max#=[Length[rank0]/2] *)

mpl=Table[Normal[Counts[ptl[[u]]]],{u,1,Length[ptl]}];(* mpl={{4->1,2->10},...} *)
Do[id=0;
     Do[nb=mpl[[u,v,2]];(* number of integers with same multiplicity *)
          mpl[[u,v,2]]=mpl[[u,v,1]];
          mpl[[u,v,1]]=Table[w+id-1,{w,1,nb}];(* mplist={{{0}->4,{1,2,3,4,5,6,7,8,9,10}->2},...} *)
          id=id+nb;
  ,{v,1,Length[mpl[[u]]]}];
,{u,1,Length[mpl]}];

num=Table[{},{u,1,Length[rank0]}];
mpnuml=Table[{mpl[[u]],num},{u,1,Length[mpl]}];

mpnuml
];

mpnuml=MpNumlistSet[lfn,rank]

insert num into tensor without symmetry

(* Insert integers into p0th tensor of (DF^3)(Subscript[Dψ, 2])γ, mp0={{0}->4,{1,2,3,4,5,6,7,8,9,10}->2} *)
InsertNum /: InsertNum[rank0_,p0_,{mp0_,num0_}]:=
Module[{rank,dm,dml,mpnuml,gid,mp,num,ps,ct,nb1,mp1,nb2,mp2,res},

dm=Flatten[Table[u,{u,1,Length[mp0]},{v,1,Length[mp0[[u,1]]]}]];(* dammy=assign group id:{1,2,2,2,2,2,2,2,2,2,2} *)
dml=Union[Subsets[dm,{rank0[[p0]]}]];(* put dammy into p0th tensor *)
mpnuml={};(* multiplicity number list *)
Do[dm=dml[[u]];(* {1,2,2,2} for rank=4 *)
     gid=Union[dm];(* group id: {1,2} *)
     mp=Map[List,mp0];(* do this to replace mp[[ps]] below *)
     num=num0;
     Do[ps=gid[[v]];(* position in mp0 *)
          ct=Count[dm,ps];(* ct times taken from mp0[[ps,1]] *)
          {nb1,mp1}={Take[mp0[[ps,1]],ct],mp0[[ps,2]]-1};(* {nb1->mp1} *)
          {nb2,mp2}={Complement[mp0[[ps,1]],nb1],mp0[[ps,2]]};(* {nb2->mp2} *)
          Which[nb2=!={}&&mp1=!=0,mp[[ps]]={nb1->mp1,nb2->mp2};,
                      nb2=!={}&&mp1==0,mp[[ps]]={nb2->mp2};,
                      nb2=={}&&mp1=!=0,mp[[ps]]={nb1->mp1};,
                      nb2=={}&&mp1==0,mp[[ps]]={};
       ];
          num[[p0]]=Join[num[[p0]],nb1];
   ,{v,1,Length[gid]}];
     mp=Flatten[mp,1];
     mpnuml=Join[mpnuml,{{mp,num}}];
,{u,1,Length[dml]}];

mpnuml
];

InsertNum[rank,9,{{{0}->4,{1,2,3,4,5,6,7,8,9,10}->2},{{},{},{},{},{},{},{},{},{}}}]

(* insert num into tensors without symmetry. pl0={9,7,8} *)
MpNumlistWoS /: MpNumlistWoS[rank0_,mpnuml0_,pl0_]:=
Module[{len,mpnuml,ps},

len=Length[rank0];
mpnuml=mpnuml0;(* multiplicity number list *)
Do[mpnuml=Flatten[Table[InsertNum[rank0,pl0[[u]],mpnuml[[v]]],{v,1,Length[mpnuml]}],1],{u,1,Length[pl0]}];

ps={};
Do[Do[If[mpnuml[[u,1,v,2]]>len-Length[pl0],(* remove if mp is larger than #remaining *)
         ps=Join[ps,{{u}}];Break[]
  ];
,{v,1,Length[mpnuml[[u,1]]]}];
,{u,1,Length[mpnuml]}];
mpnuml=Delete[mpnuml,ps];

mpnuml
];

mpnuml0=MpNumlistSet[lfn,rank];
mpnuml1=MpNumlistWoS[rank,mpnuml0,{9,5,6,7,8}](* tensors without symmetry *)

Length[mpnuml1]

insert num into tensor with symmetry

(* insert num into tensor with symmetry. mp0={{0}->1,{2,3,4,5}->1,{6,7,8,9,10}->2}, pl0={1,2,3,4,5,6} *)
InsertNumSym /: InsertNumSym[rank0_,mp0_,pl0_,sym0_]:=
Module[{len,mpnuml,ps,numl,nbl,rep,dml},

len=Length[rank0];
mpnuml={{mp0,Table[{},{u,1,len}]}};(* multiplicity number list *)
Do[mpnuml=Flatten[Table[InsertNum[rank0,pl0[[u]],mpnuml[[v]]],{v,1,Length[mpnuml]}],1],{u,1,Length[pl0]}];
ps=Map[List,pl0];(* {{1},{2},{3},{4},{5},{6}} *)
numl=Table[Extract[mpnuml[[u,2]],ps],{u,1,Length[mpnuml]}];(* num list for symmetric tensor : {{{0},{2,3,4,5},{6},{6,7,8,9},{10},{7,8,9,10}},...} *)

nbl=Table[mp0[[u,1]],{u,1,Length[mp0]}];(* nblist: {{0},{2,3,4,5},{6,7,8,9,10}} *)
rep=Flatten[Table[Map[Rule[#,u]&,nbl[[u]]],{u,1,Length[nbl]}]];(* replace nblist: {0->1,2->2,3->2,4->2,5->2,6->3,7->3,8->3,9->3,10->3} *)
dml=Table[TakeList[Flatten[numl[[u]]],sym0],{u,1,Length[numl]}];(* (((DF))^3): {{{0,2,3,4,5},{6,6,7,8,9},{10,7,8,9,10}},...} *)
dml=dml /. rep;(* dammy list: {{{1,2,2,2,2},{3,3,3,3,3},{3,3,3,3,3}},...} *)

ps={};
Do[If[OrderedQ[dml[[u]]],ps=Join[ps,{{u}}];];,{u,1,Length[numl]}];(* check ordering of (((DF))^3) *)
mpnuml=Extract[mpnuml,ps];

mpnuml
];

mpl=Union[Table[mpnuml1[[u,1]],{u,1,Length[mpnuml1]}]]

tmp=mpl[[10]];Print[tmp];
InsertNumSym[rank,tmp,{1,2,3,4},sym]

(* Complete insertion of num into tensor with and without symmetry *)
MpNumlistSym /: MpNumlistSym[rank0_,mpnuml0_,pl0_,sym0_]:=(* mpnuml0=list of {mp,num} after insertion without symmetry *)
Module[{mpl,mpnuml1,plf,mpnuml,ps},

mpl=Union[Table[mpnuml0[[u,1]],{u,1,Length[mpnuml0]}]];(* union of mp list: {{{0}->1,{2,3,4,5}->1,{6,7,8,9,10}->2},...} *)
mpnuml1=Table[InsertNumSym[rank0,mpl[[u]],pl0,sym0],{u,1,Length[mpl]}];(* mpnuml1=list of {mp,num} after insertion with symmetry *)

plf=Complement[Table[u,{u,1,Length[rank0]}],pl0];(* position list filled with number *)
mpnuml=mpnuml0;
Do[ps=Position[mpl,mpnuml0[[u,1]]][[1,1]];
     mpnuml[[u]]=mpnuml1[[ps]];
     Do[mpnuml[[u,v,2]]=ReplacePart[mpnuml[[u,v,2]],Table[plf[[i]]->mpnuml0[[u,2,plf[[i]]]],{i,1,Length[plf]}]];
  ,{v,1,Length[mpnuml[[u]]]}];
     If[mpnuml[[u]]=={},mpnuml[[u]]=Hold[Nothing];];
,{u,1,Length[mpnuml0]}];

mpnuml=Flatten[ReleaseHold[mpnuml],1];
mpnuml
];

mpnuml2=MpNumlistSym[rank,mpnuml1,{1,2,3,4},sym]

numberl=Map[Flatten,Table[mpnuml2[[u,2]],{u,1,Length[mpnuml2]}],1]

Setup for inserting number into (DF)2DDFeP2

basel=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base6";
numberl=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2num6";

indl=Table[TakeList[basel[[1,u]],rank],{u,1,Length[basel[[1]]]}];(* divide index *)
numl=Table[TakeList[numberl[[u]],rank],{u,1,Length[numberl]}];(* divide number *)
Print[Length[indl]," " ,Length[numl]];
indpsl=Map[Makeindps[fn,#]&,indl];(* make list of indps *)
numpsl=Map[Makenumps,numl];(* make list of numps *)
Print[Length[indpsl]," " ,Length[numpsl]];
indidl=Table[FromDigits[Flatten[Sort[indpsl[[u]]]]],{u,1,Length[indpsl]}];
indkeyl=Association[Table[indidl[[u]]->u,{u,1,Length[indidl]}]];
dvl=Table[{},{u,1,4}];
Do[pm=Permutations[Take[{2,3,4,5,6,7,8,9,10},2u+1]];
       pm=Table[Join[{1},pm[[u]]],{u,1,Length[pm]}];
       dvl[[u]]=Table[Partition[pm[[u]],2],{u,1,Length[pm]}];
  dvl[[u]]=Union[Map[Sort,dvl[[u]],2]];(* dvl[[1]]={{{1,2},{3,4}},{{1,3},{2,4}},{{1,4},{2,3}}} *)
,{u,1,4}];

indpsl

Union[Map[OrderedQ,indpsl]]

indidl

indkeyl

dvl

Make indps and numps

(* make indps; positions of (f1,f2,...) in ind *)
Makeindps /: Makeindps[fn0_,ind0_]:=
Module[{indps,ps},

indps={};(* position of index: {{1,9},{2,4},{2,4},...,{8,9}} for {f1,f2,...,f12} *)
Do[ps=Position[ind0,fn0[[u]]];
     ps={ps[[1,1]],ps[[2,1]]};
     indps=Join[indps,{ps}];
,{u,1,Length[fn]}];

indps
];

tmp=indl[[1]];Print[tmp];
Makeindps[fn,tmp]

(* make numps: positions of (0,1,...,10) in num *)
Makenumps /: Makenumps[num0_]:=
Module[{numps,ps},

numps={};(* num and position: {{1,7,8,9},{8,9},...,{4,6}} for {0,1,...,10} *)
Do[ps=Position[num0,n-1];
     ps=Table[ps[[v,1]],{v,1,Length[ps]}];
     numps=Join[numps,{ps}];
,{n,1,11}];

numps
];

tmp=numl[[70]];Print[tmp];
Makenumps[tmp]

Insert number into index

(* transform tensor to canonical form *)
CfNum /: CfNum[fn0_,indl0_,indpsl0_,indkey0_,num0_,numps0_,symrep0_,asym0_,dvl0_]:=
Module[{num2,num468,len468,len,rep,pml,uni,numpsl,numidl,psl,ps,numodl,nb,nbl,repl,numl,sgn0,sgnl,cfl},

num2={};(* list of num2 *)
num468={};(* list of num4,6,8 *)
len468=0;
Do[len=Length[numps0[[u]]];
     If[len==2,
          num2=Join[num2,{numps0[[u]]}];,
          rep=Table[i->numps0[[u,i]],{i,1,len}];
          num468=Join[num468,{dvl0[[len/2-1]] /. rep}];
          ++len468;
   ];
,{u,1,Length[numps0]}];
num468=Join[{Join},num468,{1}];(* {Join,list1,list2,...,1} *)
num468=Flatten[Apply[Outer,num468],len468-1];
numpsl=Table[Join[num468[[u]],num2],{u,1,Length[num468]}];(* generate numps from numps0 *)

rep=symrep0;
pml=Table[Map[Sort,numps0 /. rep[[u]]],{u,1,Length[rep]}];(* permutations by symmetry *)
uni=Union[pml];(* union of permutations *)
Do[If[SubsetQ[uni,{pml[[u]]}],uni=Complement[uni,{pml[[u]]}];,rep[[u]]=Hold[Nothing]];,{u,1,Length[pml]}];(* remove multiplicity *)
rep=ReleaseHold[rep];
numpsl=Flatten[numpsl /. rep,1];(* generate numps from numps0 by symmetry *)
numpsl=Map[Sort,numpsl,{2}];(* sort at level 2 *)

numidl=Map[FromDigits,Map[Flatten,Map[Sort,numpsl]]];(* numps to numid *)
psl=Table[indkey0[numidl[[u]]],{u,1,Length[numidl]}];(* numid to position in indpsl *)
ps=Position[psl,Missing[_,_]];(* delete missing# *)
psl=Delete[psl,ps];
numpsl=Delete[numpsl,ps];

numodl=Map[Ordering,numpsl];(* ordering of numps compared with indps *)
nb=Flatten[Table[n-1,{n,1,11},{u,1,Length[numps0[[n]]]/2}]];(* {0,0,0,1,1,2,...,10} *)
nbl=Table[nb[[numodl[[u]]]],{u,1,Length[numpsl]}];(* sort by numod *)
repl=Table[Table[fn0[[i]]->nbl[[u,i]],{i,1,Length[fn0]}],{u,1,Length[numpsl]}];(* replace list of fn *)
numl=Table[indl0[[psl[[u]]]] /. repl[[u]],{u,1,Length[numpsl]}];(* replace *)
sgn0=Apply[Times,Map[Signature,num0]];
sgnl=sgn0*Table[(Apply[Times,Map[Signature,numl[[u]]]])*(Signature[Flatten[Extract[numl[[u]],asym0]]]),{u,1,Length[numpsl]}];(* signature. note Signature[{}]=1 *)
cfl=Table[sgnl[[u]]Apply[Times,Factorial[Counts[numpsl[[u]]]]],{u,1,Length[numpsl]}];(* n! *)

Sum[cfl[[u]]Ci[psl[[u]]],{u,1,Length[numpsl]}]
];

CfNum[fn,indl,indpsl,indkeyl,numl[[10]],numpsl[[10]],symrep,asym,dvl]

Solve coefficients

CfSolve /: CfSolve[eq0_,matcf0_,cfl0_]:=
Module[{leq,len,mateq,mat,mr,ar,nb,matcf},

leq=Length[eq0];
len=Length[matcf0];
mateq=CoefficientArrays[eq0,cfl0][[2]].matcf0;(* (coefficinet of cfl0).matcf0 *)
mateq=SparseArray[RowReduce[mateq]];
mr=MatrixRank[mateq];
If[mr==0,mateq=0;,mateq=Take[mateq,mr];];(* cut by rank of mateq *)

mat=SparseArray[{{i_,i_}->1},{len,len},0];
Do[ar=ArrayRules[-mateq[[u]]];
     nb=ar[[1,1,1]];
     ar=Delete[ar,1];
     mat[[nb]]=SparseArray[ar,{len}];
,{u,1,Length[mateq]}];

matcf=matcf0.mat;
matcf=SparseArray[matcf,{len,len}]
];

eq1={Zf[1]+2Zf[2]+Zf[4],Zf[1]+3Zf[4],Zf[3],3Zf[1]+6Zf[2]+3Zf[4],5Zf[2]+Zf[3]+Zf[5]};
eq2={Zf[6]+2Zf[2]+Zf[4],Zf[7]+3Zf[4],Zf[3],3Zf[1]+6Zf[2]+3Zf[4],5Zf[2]+Zf[3]+Zf[5]};Print[MatrixForm[CoefficientArrays[eq1][[2]]],MatrixForm[CoefficientArrays[eq2][[2]]]];

mat1=CfSolve[eq1,SparseArray[{i_,i_}->1,{7,7},0],Table[Zf[i],{i,1,7}]];
mat2=CfSolve[eq2,mat1,Table[Zf[i],{i,1,7}]];Print[MatrixForm[mat1],MatrixForm[mat2]];

Solve[eq1==0&&eq2==0]

{{Zf[2]->-((Zf[1])/3),Zf[3]->0,Zf[4]->-((Zf[1])/3),Zf[5]->((5 Zf[1])/3),Zf[6]->Zf[1],Zf[7]->Zf[1]}}

Gauss Eliminate

GaussEliminate /: GaussEliminate[mat0_]:=
Module[{len1,len2,ar,mat,mr,id,ps,ps1,ps2},

len1=Length[mat0];(* length of row *)
len2=Length[mat0[[1]]];(* length of column *)
ar=ArrayRules[mat0];
mat=SparseArray[ar,{len1,len2},0];
mr=MatrixRank[mat];

id=0;
Do[ps=Flatten[ArrayRules[mat[[All,j]]][[All,1]]];(* nonzero positions of jth column = {1,2,_} *)
     If[ps[[-1]]==_,ps=Delete[ps,-1]];(* delete _ *)
     If[MemberQ[ps,n_/;(n>id)],
          ps1=FirstCase[ps,n_/;(n>id)];(* 1st position in ps which is greater than id *)
          ps2=Complement[ps,{ps1}];(* ps\ps1 *)
          mat[[ps1]]=mat[[ps1]]/mat[[ps1,j]];(* normalize mat[[ps1]] *)
          Do[mat[[ps2[[i]]]]=mat[[ps2[[i]]]]-mat[[ps2[[i]],j]]*mat[[ps1]];(* eliminate mat[[ps2[[i]],j]] *)
               ar=ArrayRules[mat[[ps2[[i]]]]];
               mat[[ps2[[i]]]]=SparseArray[ar,{len2},0];
        ,{i,1,Length[ps2]}];
          id=id+1;
         {mat[[id]],mat[[ps1]]}={mat[[ps1]],mat[[id]]};(* exchange mat *)
          If[id==mr,Break[];];
   ];
,{j,1,len2}];

mat
];

GaussEliminate[{{0,0},{0,0}}] // MatrixForm

mat1=SparseArray[{{1,2}->1,{1,3}->0,{2,1}->1,{3,2}->0,{4,5}->3}];
mat2=GaussEliminate[mat1];Print[MatrixForm[mat1],MatrixForm[mat2]];

GaussEliminate2 /: GaussEliminate2[mat0_]:=
Module[{len1,len2,ar,mat,vec,mr,id,ps,ps1,ps2,matv},

len1=Length[mat0];(* length of row *)
len2=Length[mat0[[1]]];(* length of column *)
ar=ArrayRules[mat0];
mat=SparseArray[ar,{len1,len2},0];
vec=Table[Zf[i],{i,1,len1}];(* mat0.v+vec=0 *)
mr=MatrixRank[mat];

id=0;
Do[ps=Flatten[ArrayRules[mat[[All,j]]][[All,1]]];(* nonzero positions of jth column *)
     If[ps[[-1]]==_,ps=Delete[ps,-1]];(* delete _ *)
     If[MemberQ[ps,n_/;(n>id)],
          ps1=FirstCase[ps,n_/;(n>id)];(* 1st position in ps which is greater than id *)
          ps2=Complement[ps,{ps1}];(* ps\ps1 *)
          vec[[ps1]]=Expand[vec[[ps1]]/mat[[ps1,j]]];(* normalize vec[[ps1]] *)
          mat[[ps1]]=mat[[ps1]]/mat[[ps1,j]];(* normalize mat[[ps1]] *)
          Do[vec[[ps2[[i]]]]=Expand[vec[[ps2[[i]]]]-mat[[ps2[[i]],j]]*vec[[ps1]]];(* eliminate vec[[ps2[[i]],j]] *)
               mat[[ps2[[i]]]]=mat[[ps2[[i]]]]-mat[[ps2[[i]],j]]*mat[[ps1]];(* eliminate mat[[ps2[[i]],j]] *)
               ar=ArrayRules[mat[[ps2[[i]]]]];
               mat[[ps2[[i]]]]=SparseArray[ar,{len2},0];
        ,{i,1,Length[ps2]}];
          id=id+1;
         {vec[[id]],vec[[ps1]]}={vec[[ps1]],vec[[id]]};(* exchange vec *)
         {mat[[id]],mat[[ps1]]}={mat[[ps1]],mat[[id]]};(* exchange mat *)
          If[id==mr,Break[];];
   ];
,{j,1,len2}];

matv=CoefficientArrays[vec,Table[Zf[i],{i,1,len1}]][[2]];(* mat.v+matv.vec=0 -> mat=matv.mat0 *)

{mat,matv}
];

mat1=SparseArray[{{1,2}->1,{2,1}->1,{4,5}->3}];Print[MatrixForm[mat1]];

mat2=GaussEliminate2[mat1];Print[MatrixForm[mat2[[1]]],MatrixForm[mat2[[2]]],MatrixForm[mat2[[2]].mat1]];

EOM identity

EOMidentity /: EOMidentity[basel0_,rank0_,psl0_,var0_,V0_,VEOM0_]:=(*psl0={{P2},{DP2},{DFs}}*)
Module[{leg,indl,vareomlist},

leg=rank0[[-1]];
indl=Table[TakeList[basel0[[1,u]],rank0],{u,1,Length[basel0[[1]]]}];

Which[Length[psl0[[1]]]==1,(* check P2 and DFs *)
               vareomlist=P2toP2EOM[basel0,rank0,psl0[[3]],var0,V0,VEOM0];,
            Length[psl0[[2]]]==2,(* check DP2 and DFs *)
               vareomlist=DP2toDP2EOM[basel0,rank0,psl0[[3]],var0,V0,VEOM0];,
            Length[psl0[[3]]]>0,(* check DFs *)
               vareomlist=DFtoDFEOM[basel0,rank0,psl0[[3]],var0,V0,VEOM0];
];

vareomlist
];

XP2 to XP2 + EOM

P2toP2EOM /: P2toP2EOM[basel0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,indl,psDF,ind,sec,ps,sgn,gm0,cf,psi,gm,indf,fn,vcl,vclall,vcllist},

leg=rank0[[-1]];
indl=Table[TakeList[basel0[[1,u]],rank0],{u,1,Length[basel0[[1]]]}];
psDF=Flatten[psDF0];(* flattened position of DFs *)

Do[ind=indl[[u]];
     sec=Intersection[ind[[-2]],ind[[-1]]];
     Clear[indf,vcl,vclall];
     Switch[Length[sec],
        0,If[Apply[Or,Table[SubsetQ[ind[[psDF[[2i]]]],ind[[psDF[[2i-1]]]]],{i,1,Length[psDF0]}]],(* check EOM of DF *)
                vclall=V0[leg,u]+VEOM0[leg,u];,
                vclall=0(* identity *)
         ];,
        1,ps=Position[ind[[-1]],sec[[1]]][[1,1]];
              sgn=(((-1))^(leg-ps));
              gm0=ind[[-1]] /. {sec[[1]]->Nothing};
              Do[cf=-sgn(((-1))^(leg-i));
                   psi=ind[[-2]] /. {sec[[1]]->gm0[[i]]};
                   gm=Delete[gm0,i];
                   indf=Flatten[Join[Drop[ind,-2],{psi,gm}]];
                   fn=Union[indf];
                   vcl[i]=cf ClassifyInd[var0[leg-2],indf,fn,Tns[[1+leg-2]],rotl] /. {V0[x_]->V0[leg-2,x]};(* classify *)
            ,{i,1,leg-1}];
              vclall=V0[leg,u]+Sum[vcl[i],{i,1,leg-1}];,
        2,ps=Table[Position[ind[[-1]],ind[[-2,i]]][[1,1]],{i,1,2}];
              sgn=Signature[{ps[[1]],ps[[2]]}](((-1))^(ps[[1]]+ps[[2]]+1));
              gm0=ind[[-1]] /. {sec[[1]]->Nothing,sec[[2]]->Nothing};
              Do[cf=-2sgn(((-1))^(i+j));
                   psi={gm0[[i]],gm0[[j]]};
                   gm=Delete[gm0,{{i},{j}}];
                   indf=Flatten[Join[Drop[ind,-2],{psi,gm}]];
                   fn=Union[indf];
                   vcl[i,j]=cf ClassifyInd[var0[leg-4],indf,fn,Tns[[1+leg-4]],rotl] /. {V0[x_]->V0[leg-4,x]};(* classify *)
            ,{i,1,leg-3},{j,i+1,leg-2}];
              vclall=V0[leg,u]+Sum[vcl[i,j],{i,1,leg-3},{j,i+1,leg-2}];
  ];
   vcllist[u]=vclall;
,{u,1,Length[indl]}];

Table[vcllist[u],{u,1,Length[indl]}]
];

vv1[0]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base0";(* leg=0, 587 (587-495=92) *)
vv1[2]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base2";(* leg=2, 5427 (5427-4482=945) *)
vv1[4]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base4";(* leg=4, 11320 (11320-9491=1829) *)
vv1[6]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base6";(* leg=6, 11444 (11444-9838=1606) *)
vv1[8]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base8";(* leg=8, 7054 (7054-6290=764) *)
vv1[10]=<<"Mtheory/(DF)2DDFeP2data/(DF)2DDFeP2base10";(* leg=10, 2958 (2958-2729=229) *)
Print[Length[vv1[0][[1]]]," ",Length[vv1[2][[1]]]," ",Length[vv1[4][[1]]]," ",Length[vv1[6][[1]]]," ",Length[vv1[8][[1]]]," ",Length[vv1[10][[1]]]];

leG=10;
eomid=P2toP2EOM[vv1[leG],{1,4,1,4,1,1,4,2,leG},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

Check P2 ideitity

DecomposeP2EOM /: DecomposeP2EOM[indf0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,len,psDF,fn,vcl,ind,sec,dec,gm1,gm2,sgn,tmall,cf,indf,vclall},

leg=rank0[[-1]];
len=Total[rank0]-leg;
psDF=Flatten[psDF0];(* flattened position of DFs *)
fn=Union[indf0];
vcl[1]= ClassifyInd[var0[leg],indf0,fn,Tns[[leg+1]],rotl] /. {V0[n_]->V0[leg,n]};
ind=TakeList[indf0,rank0];(* {(DF)2DDFP2γ} *)
sec=Intersection[ind[[-2]],ind[[-1]]];
dec={};

If[Length[sec]==0,
     If[Apply[Or,Table[SubsetQ[ind[[psDF[[2i]]]],ind[[psDF[[2i-1]]]]],{i,1,Length[psDF0]}]],(* check EOM of DF *)
       dec={vcl[1]->-VEOM0[leg,vcl[1][[2]]]},
       dec={vcl[1]->vcl[1]};(* identity *)
 ];,
     gm1=ind[[-1]] /. Table[sec[[u]]->Nothing,{u,1,Length[sec]}];(* decompose γ into sgn{gm1,gm2} *)
     gm2=sec;
     sgn=Signature[ind[[-1]]]Signature[Join[gm1,gm2]];
     tmall=1-Expand[sgn F1[Flatten[Delete[ind,-1]]]GammaGamma[gm1,gm2]];
     Do[cf=tmall[[u]] /. {F1[_]->1,γ[_]->1};
         indf=Flatten[tmall[[u]]/cf /. {F1[n1_]γ[n2_]->{n1,n2},F1[n1_]->{n1}}];
         fn=Union[indf];
         leg=Length[indf]-len;
         vcl[u]=cf ClassifyInd[var0[leg],indf,fn,Tns[[leg+1]],rotl] /. {V0[n_]->V0[leg,n]};
    ,{u,2,Length[tmall]}];
    vclall=Sum[vcl[u],{u,1,Length[tmall]}];
    dec={vcl[1]->vclall};
];

dec
];

<<"Mtheory/GammaProduct.m";

DecomposeP2EOM[vv1[8][[1,341]],{1,4,1,4,1,1,4,2,8},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOM[vv1[8][[1,76]],{1,4,1,4,1,1,4,2,8},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOM[vv1[2][[1,88]],{1,4,1,4,1,1,4,2,2},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOM[vv1[2][[1,14]],{1,4,1,4,1,1,4,2,2},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOM2 /: DecomposeP2EOM2[indf0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{dec,leg,ind,sec,rank,tmall,cf,nb,indf,dec1},

dec=DecomposeP2EOM[indf0,rank0,psDF0,var0,V0,VEOM0];
leg=rank0[[-1]];
ind=TakeList[indf0,rank0];(* {(DF)2DDFP2γ} *)
sec=Intersection[ind[[-2]],ind[[-1]]];(* intersection of P2 and γ *)

If[Length[sec]>=2,
     Do[leg=leg-2;
       rank=rank0;
       rank[[-1]]=leg;
       tmall=1+dec[[1,2]] /. {V0[leg,n1_]->V0[n1],V0[_,_]->0,VEOM0[_,_]->0};
       Do[cf=tmall[[v]] /. {V0[_]->1};
            nb=tmall[[v]]/cf /. {V0[n1_]->n1};
            indf=var0[leg][[1,nb]];
            dec1=DecomposeP2EOM[indf,rank,psDF0,var0,V0,VEOM0];
            dec[[1,2]]=Expand[dec[[1,2]] /. dec1];
      ,{v,2,Length[tmall]}];
,{u,1,Length[sec]-1}];
];

dec
];

DecomposeP2EOM2[vv1[8][[1,76]],{1,4,1,4,1,1,4,2,8},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOM2[vv1[2][[1,3528]],{1,4,1,4,1,1,4,2,2},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOM2[vv1[2][[1,157]],{1,4,1,4,1,1,4,2,2},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

DecomposeP2EOMAll /: DecomposeP2EOMAll[basel0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,indfl,dec,vcllist},

leg=rank0[[-1]];
indfl=basel0[[1]];
Do[dec=DecomposeP2EOM2[indfl[[u]],rank0,psDF0,var0,V0,VEOM0];
     vcllist[u]=dec[[1,1]]-dec[[1,2]];
,{u,1,Length[indfl]}];

Table[vcllist[u],{u,1,Length[indfl]}]
];

leG=10;
eomidcheck=DecomposeP2EOMAll[vv1[leG],{1,4,1,4,1,1,4,2,leG},{{1,2},{3,4},{6,7}},vv1,V1,V1EOM]

eomidcheck==eomid

XDP2 to XDP2 + EOM

DP2toDP2EOM /: DP2toDP2EOM[basel0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,indl,psDF,ind,sec,id,ps,sgn,gm0,cf,psi,gm,dr,indf,fn,vcl,vclall,vcllist},

leg=rank0[[-1]];
indl=Table[TakeList[basel0[[1,u]],rank0],{u,1,Length[basel0[[1]]]}];
psDF=Flatten[psDF0];(* flattened position of DFs *)

Do[ind=indl[[u]];
     sec={Intersection[ind[[-2]],ind[[-1]]],Intersection[ind[[-3]],ind[[-1]]],Intersection[ind[[-3]],ind[[-2]]]};(* {ψγ,Dγ,Dψ} *)
     id=Map[Length,sec];
     Clear[indf,vcl,vclall]; 
     Switch[id,(* check EOM of Dψ *)
     {0,0,0},
            If[Apply[Or,Table[SubsetQ[ind[[psDF[[2i]]]],ind[[psDF[[2i-1]]]]],{i,1,Length[psDF0]}]],(* check EOM of DF *)
               vclall=V0[leg,u]+VEOM0[leg,u];,
               vclall=0(* identity *)
       ];,
     {1,0,0},
            ps=Position[ind[[-1]],sec[[1,1]]][[1,1]];
            sgn=(((-1))^(leg-ps));
            gm0=ind[[-1]] /. {sec[[1,1]]->Nothing};
            Do[cf=-sgn(((-1))^(leg-i));
                 psi=ind[[-2]] /. {sec[[1,1]]->gm0[[i]]};
                 gm=Delete[gm0,i];
                 indf=Flatten[Join[Drop[ind,-2],{psi,gm}]];
                 fn=Union[indf];
                 vcl[i]=cf ClassifyInd[var0[leg-2],indf,fn,Tns[[1+leg-2]],rotl] /. {V0[x_]->V0[leg-2,x]};(* classify *)
          ,{i,1,leg-1}];
            vclall=V0[leg,u]+Sum[vcl[i],{i,1,leg-1}]+VEOM0[leg,u];,
     {0,1,0},
            ps=Position[ind[[-1]],sec[[2,1]]][[1,1]];
            sgn=(((-1))^(leg-ps));
            gm0=ind[[-1]] /. {sec[[2,1]]->Nothing};
            Do[cf=-sgn(((-1))^(leg-i));
                 dr=ind[[-3]] /. {sec[[2,1]]->gm0[[i]]};
                 gm=Delete[gm0,i];
                 indf=Flatten[Join[Drop[ind,-3],{dr,ind[[-2]],gm}]];
                 fn=Union[indf];
                 vcl[i]=cf ClassifyInd[var0[leg-2],indf,fn,Tns[[1+leg-2]],rotl] /. {V0[x_]->V0[leg-2,x]};(* classify *)
          ,{i,1,leg-1}];
            vclall=V0[leg,u]+Sum[vcl[i],{i,1,leg-1}]+VEOM0[leg,u];,
     {0,0,1},
            vclall=V0[leg,u]+VEOM0[leg,u];(* no term for V2 *),
     {2,0,0},
            ps=Table[Position[ind[[-1]],ind[[-2,i]]][[1,1]],{i,1,2}];
            sgn=Signature[{ps[[1]],ps[[2]]}](((-1))^(ps[[1]]+ps[[2]]+1));
            gm0=ind[[-1]] /. {sec[[1,1]]->Nothing,sec[[1,2]]->Nothing};
            Do[cf=-2sgn(((-1))^(i+j));
                 psi={gm0[[i]],gm0[[j]]};
                 gm=Delete[gm0,{{i},{j}}];
                 indf=Flatten[Join[Drop[ind,-2],{psi,gm}]];
                 fn=Union[indf];
                 vcl[i,j]=cf ClassifyInd[var0[leg-4],indf,fn,Tns[[1+leg-4]],rotl] /. {V0[x_]->V0[leg-4,x]};(* classify *)
          ,{i,1,leg-3},{j,i+1,leg-2}];
            vclall=V0[leg,u]+Sum[vcl[i,j],{i,1,leg-3},{j,i+1,leg-2}]+VEOM0[leg,u];,
     {1,1,0},
            ps={Position[ind[[-1]],sec[[2,1]]][[1,1]],Position[ind[[-1]],sec[[1,1]]][[1,1]]};
            sgn=Signature[{ps[[1]],ps[[2]]}](((-1))^(ps[[1]]+ps[[2]]+1));
            gm0=ind[[-1]] /. {sec[[2,1]]->Nothing,sec[[1,1]]->Nothing};
            Do[cf=-sgn(((-1))^(i+j));
                 gm=Delete[gm0,{{i},{j}}];
                 indf[1]=Flatten[Join[Drop[ind,-3],{{gm0[[i]]},ind[[-2]] /. {sec[[1,1]]->gm0[[j]]},gm}]];
                 indf[2]=Flatten[Join[Drop[ind,-3],{{gm0[[j]]},ind[[-2]] /. {sec[[1,1]]->gm0[[i]]},gm}]];
                 fn=Union[indf[1]];
                 vcl[i,j]=cf ClassifyInd[var0[leg-4],indf[1],fn,Tns[[1+leg-4]],rotl]
                               -cf ClassifyInd[var0[leg-4],indf[2],fn,Tns[[1+leg-4]],rotl] /. {V0[x_]->V0[leg-4,x]};(* classify *)
          ,{i,1,leg-3},{j,i+1,leg-2}];
            vclall=V0[leg,u]+Sum[vcl[i,j],{i,1,leg-3},{j,i+1,leg-2}]+VEOM0[leg,u];,
     {1,0,1},
            vclall=V0[leg,u]+VEOM0[leg,u];,
     {2,1,0},
            ps={Position[ind[[-1]],ind[[-3,1]]][[1,1]],Position[ind[[-1]],ind[[-2,1]]][[1,1]],Position[ind[[-1]],ind[[-2,2]]][[1,1]]};
            sgn=Signature[{ps[[1]],ps[[2]],ps[[3]]}](((-1))^(leg+ps[[1]]+ps[[2]]+ps[[3]]+1));
            gm0=ind[[-1]] /. {sec[[2,1]]->Nothing,sec[[1,1]]->Nothing,sec[[1,2]]->Nothing};
            Do[cf=-2sgn(((-1))^(leg+i+j+k));
                 gm=Delete[gm0,{{i},{j},{k}}];
                 indf[1]=Flatten[Join[Drop[ind,-3],{{gm0[[i]]},{gm0[[j]],gm0[[k]]},gm}]];
                 indf[2]=Flatten[Join[Drop[ind,-3],{{gm0[[j]]},{gm0[[k]],gm0[[i]]},gm}]];
                 indf[3]=Flatten[Join[Drop[ind,-3],{{gm0[[k]]},{gm0[[i]],gm0[[j]]},gm}]];
                 fn=Union[indf[1]];
                 vcl[i,j,k]=cf ClassifyInd[var0[leg-6],indf[1],fn,Tns[[1+leg-6]],rotl]
                                    +cf ClassifyInd[var0[leg-6],indf[2],fn,Tns[[1+leg-6]],rotl]
                                    +cf ClassifyInd[var0[leg-6],indf[3],fn,Tns[[1+leg-6]],rotl] /. {V0[x_]->V0[leg-6,x]};(* classify *)
          ,{i,1,leg-5},{j,i+1,leg-4},{k,j+1,leg-3}];
            vclall=V0[leg,u]+Sum[vcl[i,j,k],{i,1,leg-5},{j,i+1,leg-4},{k,j+1,leg-3}]+VEOM0[leg,u];
  ];
    vcllist[u]=vclall;
,{u,1,Length[indl]}];

Table[vcllist[u],{u,1,Length[indl]}]
];

vv2[0]=<<"Mtheory/(DF)3eDP2data/(DF)3eDP2base0";(* leg=0, 194 (194-145=49) *)
vv2[2]=<<"Mtheory/(DF)3eDP2data/(DF)3eDP2base2";(* leg=2, 1809 (1809-1287=522) *)
vv2[4]=<<"Mtheory/(DF)3eDP2data/(DF)3eDP2base4";(* leg=4, 3774 (3774-2724=1050)*)
vv2[6]=<<"Mtheory/(DF)3eDP2data/(DF)3eDP2base6";(* leg=6, 3813 (3813-2909=904) *)
vv2[8]=<<"Mtheory/(DF)3eDP2data/(DF)3eDP2base8";(* leg=8, 2345 (2345-1887=458) *)
vv2[10]=<<"Mtheory/(DF)3eDP2data/(DF)3eDP2base10";(* leg=10, 981 (981-850=131) *)
Print[Length[vv2[0][[1]]]," ",Length[vv2[2][[1]]]," ",Length[vv2[4][[1]]]," ",Length[vv2[6][[1]]]," ",Length[vv2[8][[1]]]," ",Length[vv2[10][[1]]]];

leG=10;
eomid=DP2toDP2EOM[vv2[leG],{1,4,1,4,1,4,1,2,leG},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

Check DP2 ideitity

DecomposeDP2EOM /: DecomposeDP2EOM[indf0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,len,psDF,fn,vcl,ind,sec,id,dec,int,gm1,gm2,sgn,tmall,cf,indf,vclall},

leg=rank0[[-1]];
len=Total[rank0]-leg;
psDF=Flatten[psDF0];(* flattened position of DFs *)
fn=Union[indf0];
vcl[1]= ClassifyInd[var0[leg],indf0,fn,Tns[[leg+1]],rotl] /. {V0[n_]->V0[leg,n]};(* classify *)
ind=TakeList[indf0,rank0];(* {(DF)3DP2γ} *)
sec={Intersection[ind[[-2]],ind[[-1]]],Intersection[ind[[-3]],ind[[-1]]],Intersection[ind[[-3]],ind[[-2]]]};(* intersection *)
id=Map[Length,sec];
dec={};

Which[
   id=={0,0,0},
                If[Apply[Or,Table[SubsetQ[ind[[psDF[[2i]]]],ind[[psDF[[2i-1]]]]],{i,1,Length[psDF0]}]],(* check EOM of DF *)
                  dec={vcl[1]->0},
                  dec={vcl[1]->vcl[1]};(* identity *)
           ];,
   id=={0,0,1},
                dec={vcl[1]->0};,
   True,int=Intersection[Join[ind[[-3]],ind[[-2]]],ind[[-1]]];
                gm1=ind[[-1]] /. Table[int[[u]]->Nothing,{u,1,Length[int]}];(* decompose γ into sgn{gm1,gm2} *)
                gm2=int;
                sgn=Signature[ind[[-1]]]Signature[Join[gm1,gm2]];
                tmall=1-Expand[sgn F1[Flatten[Delete[ind,-1]]]GammaGamma[gm1,gm2]];
                Do[cf=tmall[[u]] /. {F1[_]->1,γ[_]->1};
                     indf=Flatten[tmall[[u]]/cf /. {F1[n1_]γ[n2_]->{n1,n2},F1[n1_]->{n1}}];
                     fn=Union[indf];
                     leg=Length[indf]-len;
                     vcl[u]=cf ClassifyInd[var0[leg],indf,fn,Tns[[leg+1]],rotl] /. {V0[n_]->V0[leg,n]};
              ,{u,2,Length[tmall]}];
               vclall=Sum[vcl[u],{u,1,Length[tmall]}];
               dec={vcl[1]->vclall};
];

dec
];

<<"Mtheory/GammaProduct.m";
Tns={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11};(* define (((DF))^3)(Subscript[Dψ, 2])γ *)
$Assumptions=Table[Tns[[2u+1]]∈Arrays[Table[11,{i,1,18+2u}],Reals,{Antisymmetric[{2,3,4,5}],Antisymmetric[{7,8,9,10}],Antisymmetric[{12,13,14,15}],Antisymmetric[{17,18}],Antisymmetric[Table[18+i,{i,1,2u}]],{Join[{6,7,8,9,10,1,2,3,4,5},Table[i,{i,11,18+2u}]],1},{Join[{1,2,3,4,5,11,12,13,14,15,6,7,8,9,10},Table[i,{i,16,18+2u}]],1}}],{u,0,5}];(* 2u = rank of gamma matrix *)
tns={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}[[1+leG]];(* choose tensor *)
rotl={{1,{}}};

DecomposeDP2EOM[vv2[0][[1,9]],{1,4,1,4,1,4,1,2,0},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM[vv2[2][[1,224]],{1,4,1,4,1,4,1,2,2},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM[vv2[8][[1,1]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM[vv2[8][[1,6]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM[vv2[8][[1,51]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM[vv2[8][[1,989]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM2 /: DecomposeDP2EOM2[indf0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{dec,leg,ind,sec,id,rank,tmall,cf,nb,indf,dec1},

dec=DecomposeDP2EOM[indf0,rank0,psDF0,var0,V0,VEOM0];
leg=rank0[[-1]];
ind=TakeList[indf0,rank0];(* {(DF)3DP2γ} *)
sec={Intersection[ind[[-2]],ind[[-1]]],Intersection[ind[[-3]],ind[[-1]]],Intersection[ind[[-3]],ind[[-2]]]};(* intersection *)
id=Map[Length,sec];

If[Total[id]>=2,
     Do[leg=leg-2;
        rank=rank0;
        rank[[-1]]=leg;
        tmall=1+dec[[1,2]] /. {V0[leg,n1_]->V0[n1],V0[_,_]->0};
        Do[cf=tmall[[v]] /. {V0[_]->1};
             nb=tmall[[v]]/cf /. {V0[n1_]->n1};
             indf=var0[leg][[1,nb]];
             dec1=DecomposeDP2EOM[indf,rank,psDF0,var0,V0,VEOM0];
             dec[[1,2]]=Expand[dec[[1,2]] /. dec1];
      ,{v,2,Length[tmall]}];
  ,{i,1,Total[id]-1}];
];

dec
];

DecomposeDP2EOM2[vv2[2][[1,224]],{1,4,1,4,1,4,1,2,2},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM2[vv2[8][[1,1]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM2[vv2[8][[1,6]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM2[vv2[8][[1,51]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOM2[vv2[8][[1,989]],{1,4,1,4,1,4,1,2,8},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

DecomposeDP2EOMAll /: DecomposeDP2EOMAll[basel0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,indfl,dec,nb,vcllist},

leg=rank0[[-1]];
indfl=basel0[[1]];
Do[dec=DecomposeDP2EOM2[indfl[[u]],rank0,psDF0,var0,V0,VEOM0];
     nb=dec[[1,1,2]];
     vcllist[u]=dec[[1,1]]-dec[[1,2]]+VEOM0[leg,nb];
     If[dec[[1,1]]==dec[[1,2]],vcllist[u]=0;];
,{u,1,Length[indfl]}];

Table[vcllist[u],{u,1,Length[indfl]}]
];

leG=10;
eomidcheck=DecomposeDP2EOMAll[vv2[leG],{1,4,1,4,1,4,1,2,leG},{{1,2},{3,4},{5,6}},vv2,V2,V2EOM]

eomid==eomidcheck

XDF to XDF + EOM

DFtoDFEOM /: DFtoDFEOM[basel0_,rank0_,psDF0_,var0_,V0_,VEOM0_]:=
Module[{leg,indl,psDF,ind,vcllist},

leg=rank0[[-1]];
indl=Table[TakeList[basel0[[1,u]],rank0],{u,1,Length[basel0[[1]]]}];
psDF=Flatten[psDF0];(* flattened position of DFs *)

Do[ind=indl[[u]];
     If[Apply[Or,Table[SubsetQ[ind[[psDF[[2i]]]],ind[[psDF[[2i-1]]]]],{i,1,Length[psDF0]}]],
          vcllist[u]=V0[leg,u]+VEOM0[leg,u];,
          vcllist[u]=0(* identity *)
   ];
,{u,1,Length[indl]}];

Table[vcllist[u],{u,1,Length[indl]}]
];