(*Version 2023/03/23 by Kuo-Chang Chen,Ke-Ming Chang at National Tsing Hua University,Taiwan*)

(* In Albouy-Kaloshin's paper (Ann. Math. 2012) about central configurations (c.c.), the concept of singular sequences (s.s.) was introduced to investigate finiteness of c.c. The following is a program which consists of 3 algorithms. The first one automates the identification of diagrams for s.s., the second one find possible orders or primary variables, and the third one generate mass relations. The program involves exact computations, no floating point involved, so outputs are precise and rigorous. *)

(* Algorithm I: Find zw-diagrams. 
The following function 'zwmatrices' is used to finding the diagrams of singular sequences of normorlized central configurations with n bodies. 
input: the unmber of bodies;  output: diagrams *)

zwmatrices[nn_]:= (

(*BEGIN SELECTION OF $z$-MATRICES*)
ClearAll[fmx,fmx0,smx,fmxcnn,fmxdnn,smxcnn,FMXCNN,FMXDNN,FMX,FMXR,pm];
n=nn;
initialtime=AbsoluteTime[];  (* Current time. Not necessary. *)
allind=Table[k1,{k1,1,n}];
pm=Permutations@Table[k1,{k1,n}];(* Permutations of {1,2,...,n}. *)
Fstyle={Bold,RGBColor[0,0.0,0.7],13};
ConnectedMatrixQ[x2_]:=Count[Flatten[IdentityMatrix[n]+Sum[MatrixPower[x2,k2],{k2,1,n-1}]],0]==0;

zpart[{zmx4_,wmx4_}]:=zmx4;
wpart[{zmx4_,wmx4_}]:=wmx4;
trianglepart[x2_]:=Table[If[i2≤j2,x2[[i2,j2]],0],{i2,1,n},{j2,1,n}];

mxsort[x2_?ArrayQ]:=x2[[##]] & @@ (Ordering[x2~Total~{#}] &/@ {1,1}); 
(* Sort n*n symmetric matrices according to column sums first so that (many) matrices conjugate by permutation matrices are mepped to the same matrix. *)
mxreorder[x2_?ArrayQ]:=x2[[##]] & @@ (Ordering[Diagonal[x2]] &/@ {1,1}); 
(* reorder according to circled or not*)

A1=Tuples[{0,1},n(n-1)/2];a=Table[k1(k1-1)/2,{k1,1,n+1}];
T1={0}~Join~Table[k1,{k1,2,n}];
V1=Array[If[#1>#2,0,1]&,{n+1,n}];
S1=Table[If[i1>=j1,1,0],{i1,n},{j1,n}];
S1=Position[Flatten@S1,1];
S1=Table[{S1[[k1,1]]-k1+1},{k1,Length@S1}];
od={};
Do[od=AppendTo[od,n-k1];Do[od=AppendTo[od,od[[-1]]+n-j1],{j1,2,n-k1}],{k1,1,n-1}];
od=n(n-1)/2-od+1;
TransToMx[x2_]:=(A2=x2[[od]];
A2=Insert[A2,0,S1];
B2=Table[A2[[n*j2+k2]],{j2,0,n-1},{k2,1,n}];
B2=B2+Transpose[B2];
Return[B2]
);
fmx0=Table[TransToMx[A1[[k1]]],{k1,1,Length[A1]}];
fmx=Table[Table[DiagonalMatrix[V1[[n-j1+1]]]+fmx0[[k1]],{k1,Length@fmx0}],{j1,T1}];
(*Build all n*n symmetric matrices up to permutations with specified order on diagonal, and classfy them according to their trace *)
Str1=Style[" building all n*n ordered symmetric matrices",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k]]],{k,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

Do[C1=fmx[[j1]];
p1=Table[C1[[k1]][[Ordering@Total@C1[[k1]],Ordering@Total@C1[[k1]]]],{k1,1,Length@C1}];
fmx[[j1]]=DeleteDuplicates[p1];,{j1,1,n}];
Do[C1=fmx[[j1]];
p1=Table[C1[[k1]][[Ordering@Diagonal@C1[[k1]],Ordering@Diagonal@C1[[k1]]]],{k1,1,Length@C1}];
fmx[[j1]]=DeleteDuplicates[p1];,{j1,1,n}];
Str1=Style[" the first sorting",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k]]],{k,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Rule of Column Sums *)
RuleCSQ[x2_]:=Count[Total@x2,1]==0 &&Count[Total@x2,0]<n; (* A "Column-sum rule": No column sum is 1, column sums can't be all zero. *)
Do[fmx[[k1]]=Select[fmx[[k1]],RuleCSQ],{k1,1,n}];
Str1=Style[" First Rule of Column Sums",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* First Rule of Trace-2 Matrices - part 1. *)
RuleTr2MxPart1Q[x2_]:=x2[[n,n-1]]==1;
fmx[[2]]=Select[fmx[[2]],RuleTr2MxPart1Q];
Str1=Style[" First Rule of Trace-2 Matrices - part 1",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Second Rule of Triangles *)
prtriangle3[x3_]:=Flatten[Table[ x3[[i3,k3]]+x3[[j3,k3]]+x3[[k3,k3]],{k3,n-K1+1,n},{j3,1,n-K1},{i3,j3+1,n-K1}],2];
RuleTriangle2Q[x2_]:=Count[prtriangle3@x2,3]==0;(*Anoter "Rule of 3*3 Principal Minors".*)
Do[K1=k1;fmx[[k1]]=Select[fmx[[k1]],RuleTriangle2Q],{k1,2,n-2}];
Str1=Style[" Second Rule of Triangles",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* First Rule of Quadrilaterals. *)
foursum[x3_]:=Flatten[Table[x3[[i3,j3]]+x3[[i3,k3]]+x3[[i3,l3]]+x3[[j3,k3]]+x3[[j3,l3]]+x3[[k3,l3]] ,{l3,4,n},{k3,3,l3-1},{j3,2,k3-1},{i3,1,j3-1}]];
RuleQuad1Q[x2_]:= Count[foursum[x2],5]==0; (*rule 2f implies that the number of z-strokes between any given four vertices can't be 5*)
Do[fmx[[k1]]=Select[fmx[[k1]],RuleQuad1Q],{k1,1,n}];
Str1=Style[" First Rule of Quadrilaterals",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* First Rule of Triangles. *)
threelinesSP[x3_]:=Flatten[Table[(x3[[i3,j3]]+x3[[i3,k3]])*(x3[[j3,i3]]+1+x3[[j3,k3]])*(x3[[k3,i3]]+x3[[k3,j3]]+1),{k3,n-K1+2,n},{j3,n-K1+1,k3-1},{i3,1,n-K1}]];
RuleTriangle1Q[x2_]:= Count[threelinesSP[x2],8]==0; 
Do[K1=k1;fmx[[k1]]=Select[fmx[[k1]],RuleTriangle1Q],{k1,2,n-1}]; 
Str1=Style[" First Rule of Triangles",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* First Rule of Trace-2 Matrices - part 2. *)
twolinesSP[x3_]:=Flatten[Table[(x3[[i3,n]]+x3[[i3,n-1]]+x3[[n-1,n]]),{i3,1,n-2}]];
RuleTr2MxPart2Q[x2_]:= ( Count[twolinesSP[x2],2]==0); 
fmx[[2]]=Select[fmx[[2]],RuleTr2MxPart2Q];
Str1=Style[" First Rule of Trace-2 Matrices - part 2",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Rule of Trace-3 Matrices, used only for n≥6.*)
If[n>5,(
RuleTr3MxQ[x2_]:= (MXZ2=x2;ii2=n-3;While [ii2>0,If[(MXZ2[[ii2,n-2]]+MXZ2[[ii2,n-1]]+MXZ2[[ii2,n]])*Max[MXZ2[[ii2,n]]+MXZ2[[n-2,n]]+MXZ2[[n-1,n]],MXZ2[[ii2,n-1]]+MXZ2[[n-2,n-1]]+MXZ2[[n-1,n]],MXZ2[[ii2,n-2]]+MXZ2[[n-2,n-1]]+MXZ2[[n-2,n]]]==3,Return[False]];ii2--];Return[True]);
fmx[[3]]=Select[fmx[[3]],RuleTr3MxQ] ;(*this rule delete diagram with exactly 3 z-circled vertices. *)
Str1=Style[" Rule of Trace-3 Matrices",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];)];

(* Rule of Trace-0 Principal Minors. *)
Rule0DiagSubmxQ[x2_]:=(S2=Table[j2,{j2,1,n-K1}];If[K1==1,S2=S2~Join~{n}];P2=Subsets[S2];P2=Select[P2,Length@#>2&];kk2=1;While[kk2≤Length@P2,s2a=P2[[kk2]];s2b=Complement[allind,s2a];If[Count[Flatten@(x2[[s2a,s2b]]),1]==1,Return[False]];kk2++];Return[True]);
Do[K1=k1;fmx[[k1]]=Select[fmx[[k1]],Rule0DiagSubmxQ],{k1,1,n-1}];
Str1=Style[" Rule of Trace-0 Principal Minors",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Delete duplicates *)
Do[
dupind={};
Do[MXA=fmx[[i1,k1]];
TMXA=Table[MXA[[pm[[j1]],pm[[j1]]]],{j1,1,Length[pm]}];
If[Length@Intersection[TMXA,fmx[[i1]][[1;;k1-1]]]>0,dupind=Append[dupind,{k1}]],{k1,2,Length[fmx[[i1]]]}];
fmx[[i1]]=Delete[fmx[[i1]],dupind],{i1,1,n}]; 
Str1=Style[" deleting duplicates",Fstyle];
Str2=Style[ToString@Table[Length[fmx[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible ordered z-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(*END SELECTION OF $z$-MATRICES*)


(*BEGIN SELECTION OF $zw$-MATRICES*)

(*Produce all possible $w$-matrices,classify them according to traces.*)
smx=Table[{},{k1,1,n}];
Do[Do[smx[[k1]]=Union[smx[[k1]],Table [  fmx[[k1]][[i1]][[pm[[j1]],pm[[j1]]]],{j1,1,Length[pm]}]],{i1,1,Length[fmx[[k1]]]}],{k1,1,n}];
smxcnn=smx[[1]];
Do[smx[[n-k1]]=smx[[n-k1]]~Join~smx[[n-k1+1]],{k1,1,n-1}];
Str1=Style[" deleting duplicates",Fstyle];
Lsmx=Table[Length[smx[[k1]]]-Length[smx[[k1+1]]],{k1,1,n-1}];
Lsmx=Append[Lsmx,Length@smx[[n]]];
Str1=Style[" producing w-matrices",Fstyle];
Str2=Style[ToString@Lsmx,Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible w-matrices in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Classify $z$-matrices according to connectedness *)
fmxcnn=Table[{},{k1,1,n}];fmxdnn=fmxcnn;
Do[fmxcnn[[k1]]=Select[fmx[[k1]],ConnectedMatrixQ],{k1,1,n}];(*Split fmx into different categories.Begin with those having connected graphs.*)
Do[fmxdnn[[k1]]=Complement[fmx[[k1]],fmxcnn[[k1]]],{k1,1,n}];
Str1=Style[" classfying z-matrices according to connectedness",Fstyle];
Str2=Style[ToString@Table[Length[fmxcnn[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString@Table[Length[fmxdnn[[k1]]],{k1,1,n}],Fstyle];
Str4=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible connected ordered z-matrices in each class: ", Str2, ",
number of possible disconnected ordered z-matrices in each class: ", Str3, ",
cumulative time: "Str4
];Print[];

(*Produce all possible $zw$-matrices whose $z$-trace is less than or equal to $w$-trace.*)
FMXCNN={};
Do[A1=fmxcnn[[1]][[i1]];
TA1=Table[{A1,smxcnn[[j1]]},{j1,Length[smxcnn]}];
FMXCNN=FMXCNN~Join~TA1,{i1,Length[fmxcnn[[1]]]}];

FMXDNN=Table[{},{k1,1,n}];
Do[
Do[
A1=fmxdnn[[k1]][[i1]];
TA1=Table[{A1,smx[[k1]][[j1]]},{j1,Length[smx[[k1]]]}];
FMXDNN[[k1]]=FMXDNN[[k1]]~Join~TA1,{i1,Length[fmxdnn[[k1]]]}]
,{k1,1,n}];
Str1=Style[" Matching z- and w-matrices",Fstyle];
Str2=Style[ToString@Length[FMXCNN],Fstyle];
Str3=Style[ToString@Table[Length[FMXDNN[[k1]]],{k1,1,n}],Fstyle];
Str4=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices with connected z-part: ", Str2, ",
number of possible zw-matrices with disconnected z-part in each class: ", Str3, ",
cumulative time: "Str4
];Print[];

zwmxtime=AbsoluteTime[]; 

(* Rule of Circling - part 1. *)
K1=1;
RuleCirclingWQ[x2_]:=Total[Flatten@wpart[x2][[1;;n-K1,n-K1+1;;n]]];
Do[K1=k1;FMXDNN[[k1]]=Select[FMXDNN[[k1]],RuleCirclingWQ@#==0&],{k1,2,n-1}]; 
Str1=Style[" Rule of Circling - part 1",Fstyle];
Str2=Style[ToString@Table[Length[FMXDNN[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices with disconnected z-part in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Rule of Circling - part 2. *)
RuleCirclingZQ[x2_]:=(dW2=Diagonal@wpart[x2];pos0=Flatten@Position[dW2,0];pos1=Complement[allind,pos0];s2=Total[Flatten@zpart[x2][[pos0,pos1]]];Return[s2]);
Do[FMXDNN[[k1]]=Select[FMXDNN[[k1]],RuleCirclingZQ@#==0&],{k1,1,n-1}]; 
Str1=Style[" Rule of Circling - part 2",Fstyle];
Str2=Style[ToString@Table[Length[FMXDNN[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices with disconnected z-part in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Second Rule of Trace-2 Matrices. *)
Rule2ndTr2MxQ[x2_]:=(MXZ2=zpart[x2];If[Total[MXZ2[[n]]]==2&& Tr[wpart[x2]]==n-3,Return[False]]; If[n==5, MXW2=wpart[x2];If[Tr[MXW2]==2,MXW2=mxreorder[MXW2];If[MXW2[[5,4]]==1&&Total[MXW2[[5]]]==2&& Tr[MXZ2]==2,Return[False]]]];   Return[True]);
FMXDNN[[2]]=Select[FMXDNN[[2]],Rule2ndTr2MxQ];(*This rule only deletes ZW-matrices with exactly two z-circles*)
Str1=Style[" Second Rule of Trace-2 Matrices",Fstyle];
Str2=Style[ToString@Table[Length[FMXDNN[[k1]]],{k1,1,n}],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices with disconnected z-part in each class: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Combine all classes *)
FMX=FMXCNN;Do[FMX=FMX~Join~FMXDNN[[k1]],{k1,1,n}];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Third Rule of Triangles *)
RuleTriangle3Q[x2_]:=(X2=zpart[x2]+wpart[x2];kk2=n;While[kk2≥3,
jj2=2;While[jj2<kk2,
ii2=1;While[ii2<jj2,discri2=X2[[ii2,jj2]]+X2[[jj2,kk2]]+X2[[kk2,ii2]];If[discri2==4||discri2==5,Return[False]];ii2++];jj2++];kk2--];Return[True]);
FMX=Select[FMX,RuleTriangle3Q];
Str1=Style[" Third Rule of Triangles",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Rule of Two Column Sums *)
Rule2CSQ[x2_]:=(MXZ2=zpart[x2];MXZ2=MXZ2-DiagonalMatrix[Diagonal[MXZ2]];MXW2=wpart[x2];MXW2=MXW2-DiagonalMatrix[Diagonal[MXW2]];X2=MXZ2+MXW2;pos2=Position[X2,2];If[Length[pos2]==2,Return[False]];kk2=1;While[kk2≤Length[pos2],If[(Total[MXZ2[[pos2[[kk2,1]]  ]]]+Total[MXZ2[[pos2[[kk2,2]]  ]]]-2)*(Total[MXW2[[pos2[[kk2,1]]  ]]]+Total[MXW2[[pos2[[kk2,2]]  ]]]-2)==0,Return[False]];kk2++];Return[True]);
FMX=Select[FMX,Rule2CSQ];
Str1=Style[" Rule of Two Column Sums",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Fourth Rule of Triangles *)
RuleTriangle4Q[x2_]:=(MXZ2=zpart[x2];MXW2=wpart[x2];X2=MXZ2+MXW2;kk2=n;While[kk2≥3,
jj2=2;While[jj2<kk2,
ii2=1;While[ii2<jj2,s2=X2[[ii2,jj2]]*X2[[jj2,kk2]]*X2[[kk2,ii2]]*Max[(MXZ2[[ii2,jj2]]+MXZ2[[kk2,jj2]]+MXZ2[[ii2,kk2]]),(MXW2[[ii2,jj2]]+MXW2[[kk2,jj2]]+MXW2[[ii2,kk2]])];If[  s2!=0&&s2!=3&&  s2!=24 ,Return[False]];ii2++];jj2++];kk2--];Return[True]);
FMX=Select[FMX,RuleTriangle4Q];
Str1=Style[" Fourth Rule of Triangles",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Rule of Connected Components. *)
componentsZ[x3_]:= (MXZ3=zpart[x3];ZC3=MatrixPower[MXZ3,n-1];zc3=Flatten[Tuples[{{0}},n-1],1];Pz3=Table[k3,{k3,1,n}];kk3=1;While[Length[Pz3]>0, zc3[[kk3]]=Union[{Pz3[[1]]},Select[Pz3,ZC3[[Pz3[[1]],#]]>0&]];Pz3=Complement[Pz3,zc3[[kk3]]] ;kk3++];zc3=Select[zc3,Length[#]>1&];Return[zc3]);

componentsW[x3_]:= (MXW3=wpart[x3];WC3=MatrixPower[MXW3,n-1];wc3=Flatten[Tuples[{{0}},n-1],1];Pw3=Table[k3,{k3,1,n}];kk3=1;While[Length[Pw3]>0, wc3[[kk3]]=Union[{Pw3[[1]]},Select[Pw3,WC3[[Pw3[[1]],#]]>0&]];Pw3=Complement[Pw3,wc3[[kk3]]] ;kk3++];wc3=Select[wc3,Length[#]>1&];Return[wc3]);

RuleConnCompZQ[x2_]:=(MXZ2=zpart[x2];MXW2=wpart[x2];X2=MXZ2+MXW2;X2=X2-DiagonalMatrix[Diagonal[X2]];If[Tr[MXZ2]<3,Return[True]];zc2=Select[componentsZ[x2],Length[#]≥3&];kk2=1;While[kk2≤Length[zc2],ll2=Length[zc2[[kk2]]];subX2=X2[[zc2[[kk2]],zc2[[kk2]]]];subZ2=MXZ2[[zc2[[kk2]],zc2[[kk2]]]];If[Tr[subZ2]==1,Return[False]];If[Tr[subZ2]==ll2&&Count[Flatten[subX2],2]==ll2(ll2-1),Return[False]];kk2++];Return[True]);

RuleConnCompWQ[x2_]:=(MXW2=wpart[x2];MXZ2=zpart[x2];X2=MXZ2+MXW2;X2=X2-DiagonalMatrix[Diagonal[X2]];If[Tr[MXW2]<3,Return[True]];wc2=Select[componentsW[x2],Length[#]≥3&];kk2=1;While[kk2≤Length[wc2],ll2=Length[wc2[[kk2]]];subX2=X2[[wc2[[kk2]],wc2[[kk2]]]];subW2=MXW2[[wc2[[kk2]],wc2[[kk2]]]];If[Tr[subW2]==1,Return[False]];If[Tr[subW2]==ll2&&Count[Flatten[subX2],2]==ll2(ll2-1),Return[False]];kk2++];Return[True]);

FMX=Select[FMX,RuleConnCompZQ];
FMX=Select[FMX,RuleConnCompWQ];
Str1=Style[" Rule of Connected Components",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Second Rule of Quadrilaterals *)
odminortest[x3_]:=Flatten@Table[{If[l3==j3||l3==k3||(zpart[x3][[i3,j3]]+wpart[x3][[i3,j3]])(zpart[x3][[i3,k3]]+wpart[x3][[i3,k3]])(zpart[x3][[l3,j3]]+wpart[x3][[l3,j3]])(zpart[x3][[l3,k3]]+wpart[x3][[l3,k3]])==0,{},{zpart[x3][[i3,j3]]+zpart[x3][[l3,k3]],zpart[x3][[i3,k3]]+zpart[x3][[l3,j3]],wpart[x3][[i3,j3]]+wpart[x3][[l3,k3]],wpart[x3][[i3,k3]]+wpart[x3][[l3,j3]]}]},{k3,3,n},{i3,1,k3-1},{l3,i3+1,n},{j3,i3+1,k3-1}];
RuleQuad2Q[x2_]:=Count[odminortest[x2],1]==0;
FMX=Select[FMX,RuleQuad2Q];
Str1=Style[" Second Rule of Quadrilaterals",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

(* Rule of Pentagons *)
RulePentagonQ[x2_]:=(MXZ2=zpart[x2];MXW2=wpart[x2];X2=MXZ2+2MXW2;ii2=1;While[ii2<n-3,jj2=ii2+1;
While[jj2≤n,kk2=ii2+1;
While[kk2≤n,ll2=ii2+1;
While[ll2≤n,mm2=jj2+1;
While[mm2≤n,
If[(jj2-kk2)(kk2-ll2)(ll2-mm2)(jj2-ll2)(kk2-mm2)(jj2-mm2)≠0,side2={X2[[ii2,jj2]],X2[[jj2,kk2]],X2[[kk2,ll2]],X2[[ll2,mm2]],X2[[mm2,ii2]]};v2={Count[side2,0],Count[side2,1],Count[side2,2],Count[side2,3]};If[v2[[1]]==0&&(v2[[2]]==1||v2[[3]]==1||v2=={0,4,0,1}||v2=={0,0,4,1}),Return[False]]];mm2++;];ll2++;];kk2++;];jj2++;];ii2++;];Return[True]);
FMX=Select[FMX,RulePentagonQ];
Str1=Style[" Rule of Pentagons",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

mxruletime=AbsoluteTime[]; 

(* Delete duplicates *)
duplicateindices={};
Do[A1=FMX[[k1]];
TA1=Flatten[Table[{{A1[[1]][[pm[[j1]],pm[[j1]]]],A1[[2]][[pm[[j1]],pm[[j1]]]]},{A1[[2]][[pm[[j1]],pm[[j1]]]],A1[[1]][[pm[[j1]],pm[[j1]]]]}},{j1,1,Length[pm]}],1];
If[Length@Intersection[TA1,FMX[[1;;k1-1]]]>0,duplicateindices=Append[duplicateindices,{k1}]],{k1,2,Length@FMX}];
FMX=Delete[FMX,duplicateindices];

Str1=Style[" deleting duplicate",Fstyle];
Str2=Style[ToString@Length[FMX],Fstyle];
Str3=Style[ToString[AbsoluteTime[]-initialtime],Fstyle];
Print["After" ,Str1,",
number of possible zw-matrices: ", Str2, ",
cumulative time: "Str3
];Print[];

zwmxtime2=AbsoluteTime[]; 

(* END SELECTION OF zw-MATRICES *)


(* BEGIN DRAWING DIAGRAMS *)

FMXR=Table[{trianglepart[FMX[[k1,1]]],trianglepart[FMX[[k1,2]]]},{k1,Length[FMX]}]; (* Remove the lower left triangular parts for square matrices on the left and on the right *)

rededges[k2_]:=Rule@@@Position[zpart@FMXR[[k2]],1]; (* Collect red edges; i.e. z-edges. *)
bluedges[k2_]:=Rule@@@Position[wpart@FMXR[[k2]],1];(* Collect blue edges; i.e. w-edges. *)
alledges[k2_]:=Flatten[Append[rededges[k2],bluedges[k2]]]; (* Collect all edges. *)
Replace[Table[GraphPlot[DeleteDuplicates[Table[1->j1,{j1,2,n}]~Join~alledges[k1]],VertexLabeling->False,(* Delete duplicate edges, then determine color of edges. *)EdgeRenderingFunction->(Which[MemberQ[rededges[k1],First[#2]->Last[#2]]==True&&MemberQ[bluedges[k1],First[#2]->Last[#2]]==False,{RGBColor[1,0.3,0.3],Thickness[0.025],Line[#1]},MemberQ[rededges[k1],First[#2]->Last[#2]]==False&&MemberQ[bluedges[k1],First[#2]->Last[#2]]==True,{RGBColor[0.2,0.7,1],Thickness[0.025],Line[#1]},MemberQ[rededges[k1],First[#2]->Last[#2]]==True&&MemberQ[bluedges[k1],First[#2]->Last[#2]]==True,{Thickness[0.025],RGBColor[0.5,0.3,0.9],Line[#1]},MemberQ[rededges[k1],First[#2]->Last[#2]]==False&&MemberQ[bluedges[k1],First[#2]->Last[#2]]==False,{Thickness[0.025]}]&),(* When the edge comes from the left side, not from the right side, then it is red; if from the right side, not from the left side, then it is blue; if from both sides, then it is purple (zw-edge). *)VertexCoordinateRules->Flatten[Table[{k1->{Cos[2(k1-1)π/n],Sin[2(k1-1)π/n]}},{k1,n}],1],(* Place vertices at n-th roots of unity, in the counterclockwise order. *)VertexRenderingFunction->({#2/.Table[Which[FMXR[[k1,1]][[i1]][[i1]]==1&&FMXR[[k1,2]][[i1]][[i1]]==0,i1->RGBColor[1,0.3,0.3],FMXR[[k1,1]][[i1]][[i1]]==0&&FMXR[[k1,2]][[i1]][[i1]]==1,i1->RGBColor[0.2,0.7,1],FMXR[[k1,1]][[i1]][[i1]]==1&&FMXR[[k1,2]][[i1]][[i1]]==1,i1->RGBColor[0.5,0.3,0.9],FMXR[[k1,1]][[i1]][[i1]]==0&&FMXR[[k1,2]][[i1]][[i1]]==0,i1->White],{i1,n}],EdgeForm[Black],Disk[#1,.24],Black,FontSize->28,FontWeight->Bold,Text[#2,#1],Opacity[0.5],RGBColor[0,0.5,0],FontSize->24,Text[k1,{1.1,-1}]}&),SelfLoopStyle->None,ImageSize->{240,240}],{k1,1,Length[FMXR]}],List->Print,1,Heads->True];(* Color vertices: z-circles red, w-circles blue, zw-circles purple *)

graphingtime=AbsoluteTime[]; 

(* END DRAWING DIAGRAMS *)

Print["Number of possibilities: ",Length[FMXR],". Total computation time in seconds: ",AbsoluteTime[]-initialtime,". Time used for Z-matrix sorting, finding possible companion W-matrices, match Z and W matrices to ZW-matrices, applying matrix rules, duplicate ZW-matrix eliminations, and graphing:",{zwmxtime-initialtime,mxruletime-zwmxtime,zwmxtime2-mxruletime,graphingtime-zwmxtime2}]  (*Show the number of possibilities,and total computaion time in seconds.*)
)

zwmatrices[6]

(* Algorithm II: Find possible orders of primary variables *)

(* Algorithm II must be compiled after Algorithm I. One may use command FindOrder[n] or FindOrder[n,t] to generate the optimal $zw$-order matrices of Type~1, and the collections C_2 of Type-2 order matrices,the collection C_3 of Type-3 order matrices, for diagram $n$. A default time constraint is set to prevent from enumerating all Type~2 or Type~3 $zw$-order matrices
when there are too many possibilities. The optional argument $t$ is to change this time constraint. The optimal $zw$-order matrices of Type~1 is denoted {zo1,wo1}, and the collection C_2 is denoted zwotype2 in the program. These will be used in applications of Algorithm III. *)

(* Generate $r$-order matrices from $zw$-order matrices *)
Dis={{    {}  ,       {}             ,        {}            ,         {}          ,{}   },
        {   {}   ,       {}             ,        {}            ,{2,3,4}  ,{4}  },
        {   {}  ,        {}            ,        {}          ,{2,3,4}  ,{4}  },
        {   {}  ,{2,3,4}    ,{2,3,4}   ,{2,3,4}  ,{4}  },
        {  {} ,      {4}           ,        {4}         ,        {4}         ,{5}  }};
DisEdge= 
        {{    {}  ,       {}            ,       {}          ,         {}        ,{3}   },
        {   {}   ,       {}           ,        {}          ,        {2}      ,{}  },
        {   {}  ,        {}          ,        {1}        ,         {}        ,{}  },
        {   {}  ,      {2}         ,         {}          ,         {}        ,{}  },
        {  {3} ,      {}           ,        {}          ,          {}        ,{}  }};
Dis2[zorderset6_,worderset6_]:=DeleteDuplicates[Flatten[Table[Dis[[zorderset6[[i6]],worderset6[[j6]]]],{i6,1,Length[zorderset6]},{j6,1,Length[worderset6]}]]];
Dis2Edge[zorderset6_,worderset6_]:=DeleteDuplicates[Flatten[Table[DisEdge[[zorderset6[[i6]],worderset6[[j6]]]],{i6,1,Length[zorderset6]},{j6,1,Length[worderset6]}]]];
Dis3[{zomatrix5_,womatrix5_}]:=(dx5=Table[{},{i5,1,n},{j5,1,n}];Do[Do[If[XX1[[i5,j5]]==0,dx5[[i5,j5]]=Dis2[zomatrix5[[i5,j5]],womatrix5[[i5,j5]]],dx5[[i5,j5]]=Dis2Edge[zomatrix5[[i5,j5]],womatrix5[[i5,j5]]]];dx5[[j5,i5]]=dx5[[i5,j5]],{i5,j5+1,n}],{j5,1,n-1}];Return[dx5]);

(* PRINCIPLES FOR UPDATING ORDERS *)

(* Proposition 4.5:triangle inequalities *) 
TriIneqQ[order5a_,order5b_,order5c_]:=(If[Count[{order5a,order5b,order5c},Max[order5a,order5b,order5c]]==1,Return[False],Return[True]]);
TriIneqQ2[order4_,orderset4A_,orderset4B_]:=(ii4=1;While[ii4≤Length[orderset4A],jj4=1;While[jj4≤Length[orderset4B],If[TriIneqQ[order4,orderset4A[[ii4]],orderset4B[[jj4]]]==True,Return[True]];jj4++];ii4++];Return[False]);
TriIneq[{zomatrix3_,womatrix3_}]:=(Clear[zo3,wo3];zo3=zomatrix3;wo3=womatrix3;
Do[Do[
zo3[[i3,i3]]=Select[zo3[[i3,i3]],TriIneqQ2[#,zo3[[j3,j3]],zo3[[i3,j3]]]&];
zo3[[i3,j3]]=Select[zo3[[i3,j3]],TriIneqQ2[#,zo3[[j3,j3]],zo3[[i3,i3]]]&];
zo3[[j3,j3]]=Select[zo3[[j3,j3]],TriIneqQ2[#,zo3[[i3,i3]],zo3[[i3,j3]]]&];
zo3[[j3,i3]]=zo3[[i3,j3]];
wo3[[i3,i3]]=Select[wo3[[i3,i3]],TriIneqQ2[#,wo3[[j3,j3]],wo3[[i3,j3]]]&];
wo3[[i3,j3]]=Select[wo3[[i3,j3]],TriIneqQ2[#,wo3[[j3,j3]],wo3[[i3,i3]]]&];
wo3[[j3,j3]]=Select[wo3[[j3,j3]],TriIneqQ2[#,wo3[[i3,i3]],wo3[[i3,j3]]]&];
wo3[[j3,i3]]=wo3[[i3,j3]];,{i3,j3+1,n}],{j3,1,n-1}];
Do[Do[Do[
zo3[[i3,j3]]=Select[zo3[[i3,j3]],TriIneqQ2[#,zo3[[j3,k3]],zo3[[i3,k3]]]&];
zo3[[j3,k3]]=Select[zo3[[j3,k3]],TriIneqQ2[#,zo3[[i3,j3]],zo3[[i3,k3]]]&];
zo3[[i3,k3]]=Select[zo3[[i3,k3]],TriIneqQ2[#,zo3[[i3,j3]],zo3[[j3,k3]]]&];
zo3[[j3,i3]]=zo3[[i3,j3]];
zo3[[k3,i3]]=zo3[[i3,k3]];
zo3[[k3,j3]]=zo3[[j3,k3]];
wo3[[i3,j3]]=Select[wo3[[i3,j3]],TriIneqQ2[#,wo3[[j3,k3]],wo3[[i3,k3]]]&];
wo3[[j3,k3]]=Select[wo3[[j3,k3]],TriIneqQ2[#,wo3[[i3,j3]],wo3[[i3,k3]]]&];
wo3[[i3,k3]]=Select[wo3[[i3,k3]],TriIneqQ2[#,wo3[[i3,j3]],wo3[[j3,k3]]]&];
wo3[[j3,i3]]=wo3[[i3,j3]];
wo3[[k3,i3]]=wo3[[i3,k3]];
wo3[[k3,j3]]=wo3[[j3,k3]];,{i3,j3+1,n}],{j3,k3+1,n-1}],{k3,1,n-2}];
Return[{zo3,wo3}];)

(* Proposition 4.6 *)
Rule2e[{zomatrix3_,womatrix3_}]:=(Clear[zo3,wo3];zo3=zomatrix3;wo3=womatrix3;Do[If[i3≠j3&&i3≠k3&&j3≠k3&&(XX1[[i3,j3]]*XX1[[j3,k3]]*XX1[[k3,i3]]≠0),zo3[[i3,j3]]=Intersection[zo3[[i3,j3]],zo3[[j3,k3]],zo3[[k3,i3]]];zo3[[j3,k3]]=zo3[[i3,j3]];zo3[[k3,i3]]=zo3[[i3,j3]];wo3[[i3,j3]]=Intersection[wo3[[i3,j3]],wo3[[j3,k3]],wo3[[k3,i3]]];wo3[[j3,k3]]=wo3[[i3,j3]];wo3[[k3,i3]]=wo3[[i3,j3]]];,{i3,1,n},{j3,1,n},{k3,1,n}];Return[{zo3,wo3}])

(* Proposition 4.7(a) *)
Estimate[{zomatrix3_,womatrix3_}]:=(Clear[zo3,wo3];zo3=zomatrix3;wo3=womatrix3;Do[If[i3≠j3,If[XX1[[i3,j3]]>0,zo3[[i3,j3]]=Intersection[zo3[[i3,j3]],Table[6-wo3[[i3,j3]][[k3]],{k3,1,Length[wo3[[i3,j3]]]}]];
wo3[[i3,j3]]=Intersection[wo3[[i3,j3]],Table[6-zo3[[i3,j3]][[k3]],{k3,1,Length[zo3[[i3,j3]]]}]]
];If[XX1[[i3,j3]]==0,zo3[[i3,j3]]=Select[zo3[[i3,j3]],#≥6-2*Floor[Max[wo3[[i3,j3]]]/2]&];wo3[[i3,j3]]=Select[wo3[[i3,j3]],#≥6-2*Floor[Max[zo3[[i3,j3]]]/2]& ]]],{i3,1,n},{j3,1,n}];Return[{zo3,wo3}]);

(* Proposition 4.7(b) *)
CenterOfMass[{zomatrix3_,womatrix3_}]:=(Clear[zo3,wo3];zo3=zomatrix3;wo3=womatrix3;MZ3=Max@Flatten[Table[Table[zo3[[i3,j3]],{i3,j3+1,n}],{j3,1,n-1}]];MW3=Max@Flatten[Table[Table[wo3[[i3,j3]],{i3,j3+1,n}],{j3,1,n-1}]];Do[zo3[[i3,i3]]=Select[zo3[[i3,i3]],#≤MZ3&];wo3[[i3,i3]]=Select[wo3[[i3,i3]],#≤MW3&],{i3,1,n}];Return[{zo3,wo3}];)

(* Determine connected components of $z$-matrix,$w$-matrix *)
componentsZ[x3_]:= (MXZ3=zpart[x3];ZC3=MatrixPower[MXZ3,n-1];zc3=Flatten[Tuples[{{0}},n-1],1];Pz3=allind;kk3=1;While[Length[Pz3]>0, zc3[[kk3]]=Union[{Pz3[[1]]},Select[Pz3,ZC3[[Pz3[[1]],#]]>0&]];Pz3=Complement[Pz3,zc3[[kk3]]] ;kk3++];zc3=Select[zc3,Length[#]>1&];Return[zc3]);

componentsW[x3_]:= (MXW3=wpart[x3];WC3=MatrixPower[MXW3,n-1];wc3=Flatten[Tuples[{{0}},n-1],1];Pw3=allind;kk3=1;While[Length[Pw3]>0, wc3[[kk3]]=Union[{Pw3[[1]]},Select[Pw3,WC3[[Pw3[[1]],#]]>0&]];Pw3=Complement[Pw3,wc3[[kk3]]] ;kk3++];wc3=Select[wc3,Length[#]>1&];Return[wc3]);

(* Proposition 4.9 *)
Prop1[{zomatrix3_,womatrix3_}]:=(zo3=zomatrix3;wo3=womatrix3;Do[If[i3≠j3,kk3=1;Z3A={};W3A={};
While[kk3≤n,
If[kk3≠i3&&kk3≠j3,Z3A=Union[Z3A,zo3[[i3,kk3]]];W3A=Union[W3A,wo3[[i3,kk3]]]];
kk3++];
Z3B=Union[Z3A,zo3[[i3,j3]]];W3B=Union[W3A,wo3[[i3,j3]]];
If[(Z3B=={1}||Z3B=={3}||Z3B=={5})&&Max[wo3[[i3,j3]]]<Min[W3A],wo3[[i3,i3]]=Select[wo3[[i3,i3]],#≤Max@wo3[[i3,j3]]& ];wo3[[j3,j3]]=Select[wo3[[j3,j3]],#≤Max@wo3[[i3,j3]]& ];wo3[[i3,j3]]=Select[wo3[[i3,j3]],#≥Min@(wo3[[i3,i3]]~Join~wo3[[j3,j3]])& ];];If[(W3B=={1}||W3B=={3}||W3B=={5})&&Max[zo3[[i3,j3]]]<Min[Z3A],zo3[[i3,i3]]=Select[zo3[[i3,i3]],#≤Max@zo3[[i3,j3]]& ];zo3[[j3,j3]]=Select[zo3[[j3,j3]],#≤Max@zo3[[i3,j3]]& ];zo3[[i3,j3]]=Select[zo3[[i3,j3]],#≥Min@(zo3[[i3,i3]]~Join~zo3[[j3,j3]])& ];]],{i3,1,n},{j3,1,n}];Return[{zo3,wo3}])

(* Proposition 4.10 *)
Prop2[{zomatrix3_,womatrix3_}]:=(zo3=zomatrix3;wo3=womatrix3;Do[dis3={};jj3=1;While[jj3≤n,If[jj3≠i3,If[XX1[[i3,jj3]]==0,dis3=Union[dis3,Dis2[zo3[[i3,jj3]],wo3[[i3,jj3]]]],dis3=Union[dis3,Dis2Edge[zo3[[i3,jj3]],wo3[[i3,jj3]]]]]];jj3++];If[Complement[dis3,{4,5}]=={},ZO3={};WO3={};jj3=1;While[jj3≤n,If[jj3≠i3,ZO3=ZO3~Join~zo3[[jj3,jj3]];WO3=WO3~Join~wo3[[jj3,jj3]]];jj3++];If[Max[ZO3]==2||Max[ZO3]==4,zo3[[i3,i3]]=Select[zo3[[i3,i3]],#<=Max[ZO3]&],zo3[[i3,i3]]=Select[zo3[[i3,i3]],#<Max[ZO3]&]];If[Max[WO3]==2||Max[WO3]==4,wo3[[i3,i3]]=Select[wo3[[i3,i3]],#<=Max[WO3]& ],wo3[[i3,i3]]=Select[wo3[[i3,i3]],#<Max[WO3]& ]]],{i3,1,n}];Return[{zo3,wo3}]);

(* Compute the number of possible levels of orders (1) for all off-diagonal entries;(2) for all entries *)
NumberofType2[{zomatrix3_,womatrix3_}]:=(kk3=1;Do[Do[kk3=kk3*Length[zomatrix3[[i3,j3]]]*Length[womatrix3[[i3,j3]]],{i3,j3+1,n}],{j3,1,n-1}];Return[kk3]);
NumberofType3[{zomatrix3_,womatrix3_}]:=(kk3=1;Do[Do[kk3=kk3*Length[zomatrix3[[i3,j3]]]*Length[womatrix3[[i3,j3]]],{i3,j3,n}],{j3,1,n}];Return[kk3]);

(* Repeat above processes until possible levels becomes invariant *)
UpdateOrder[{zomatrix2_,womatrix2_}]:=(zo2=zomatrix2;wo2=womatrix2;l2=NumberofType3[{zo2,wo2}];L2=l2+1;While[l2≠L2,
l2=NumberofType3[{zo2,wo2}];
{zo2,wo2}=TriIneq[{zo2,wo2}];
{zo2,wo2}=Rule2e[{zo2,wo2}];
{zo2,wo2}=Estimate[{zo2,wo2}];
{zo2,wo2}=CenterOfMass[{zo2,wo2}];
{zo2,wo2}=Prop1[{zo2,wo2}];
{zo2,wo2}=Prop2[{zo2,wo2}];
L2=NumberofType3[{zo2,wo2}];
];
Return[{zo2,wo2}]);

UpdateOrderPart[{zomatrix2_,womatrix2_}]:=(zo2=zomatrix2;wo2=womatrix2;l2=NumberofType3[{zo2,wo2}];L2=l2+1;While[l2≠L2,
l2=NumberofType3[{zo2,wo2}];
{zo2,wo2}=TriIneq[{zo2,wo2}];
L2=NumberofType3[{zo2,wo2}];
];
Return[{zo2,wo2}]);

(* PosOffDiag is to find a position of an off-diagonal entry that contains more than one element; PosDiag is similar, but for diagonal entries *)
PosOffDiag[x2_]:=(ii2=1;While[ii2≤n-1,jj2=ii2+1;While[jj2≤n,L2=Length[x2[[ii2,jj2]]];If[L2>1,Return[{ii2,jj2,L2}]];jj2++];ii2++];Return[{1,1,0}])
PosDiag[x2_]:=(ii2=1;While[ii2≤n,L2=Length[x2[[ii2,ii2]]];If[L2>1,Return[{ii2,L2}]];ii2++];Return[{1,0}])


(* CRITERIA FOR FEASIBLE ORDER MATRICES *)

(* Proposition 4.11 *)
Rule2hQ[{zomatrix2_,womatrix2_}]:=(zo2=zomatrix2;wo2=womatrix2;dx2=Flatten[Dis3[{zo2,wo2}],1] ;If[Count[dx2,{2}]==2&&Count[dx2,{1}]==0,Return[False],Return[True]])

(* Compare orders,find candidates for max order terms among Z_{ij}'s, W_{ij}'s *)

(* Compare z_{ij}'s or w_{ij}'s *)
CompareOrderzw[ordermatrix5_,{ii5_,jj5_},{kk5_,ll5_}]:=(o5=ordermatrix5;

If[Max@o5[[ii5,jj5]]<Min@o5[[kk5,ll5]],Return[{-1}]];
If[Min@o5[[ii5,jj5]]>Max@o5[[kk5,ll5]],Return[{1}]];
If[Length[o5[[ii5,jj5]]]==1&&o5[[ii5,jj5]]==o5[[kk5,ll5]]&&Complement[o5[[ii5,jj5]],{1,3,5}]==0,Return[{0}]];

If[ii5≠jj5&&kk5≠ll5,

If[ii5==kk5&&Min@o5[[ii5,jj5]]≥Max@o5[[jj5,ll5]]&&Min@o5[[kk5,ll5]]≥Max@o5[[jj5,ll5]],If[Min@o5[[ii5,jj5]]>Max@o5[[jj5,ll5]]||Min@o5[[kk5,ll5]]>Max@o5[[jj5,ll5]]||MemberQ[{1,3,5},Max@o5[[jj5,ll5]]]==True,Return[{0}]]];
If[jj5==ll5&&Min@o5[[ii5,jj5]]≥Max@o5[[ii5,kk5]]&&Min@o5[[kk5,ll5]]≥Max@o5[[ii5,kk5]],If[Min@o5[[ii5,jj5]]>Max@o5[[ii5,kk5]]||Min@o5[[kk5,ll5]]>Max@o5[[ii5,kk5]]||MemberQ[{1,3,5},Max@o5[[ii5,kk5]]]==True,Return[{0}]]];
If[ii5==ll5&&Min@o5[[ii5,jj5]]≥Max@o5[[jj5,kk5]]&&Min@o5[[kk5,ll5]]≥Max@o5[[jj5,kk5]],If[Min@o5[[ii5,jj5]]>Max@o5[[jj5,kk5]]||Min@o5[[kk5,ll5]]>Max@o5[[jj5,kk5]]||MemberQ[{1,3,5},Max@o5[[jj5,kk5]]]==True,Return[{0}]]];
If[kk5==jj5&&Min@o5[[ii5,jj5]]≥Max@o5[[ll5,ii5]]&&Min@o5[[kk5,ll5]]≥Max@o5[[ll5,ii5]],If[Min@o5[[ii5,jj5]]>Max@o5[[ll5,ii5]]||Min@o5[[kk5,ll5]]>Max@o5[[ll5,ii5]]||MemberQ[{1,3,5},Max@o5[[ll5,ii5]]]==True,Return[{0}]]]];

If[ii5==jj5&&kk5≠ll5,
If[ii5==kk5&&Min@o5[[ii5,jj5]]≥Max@o5[[ll5,ll5]]&&Min@o5[[kk5,ll5]]≥Max@o5[[ll5,ll5]],If[Min@o5[[ii5,jj5]]>Max@o5[[ll5,ll5]]||Min@o5[[kk5,ll5]]>Max@o5[[ll5,ll5]]||MemberQ[{1,3,5},Max@o5[[ll5,ll5]]]==True,Return[{0}]]];
If[ii5==ll5&&Min@o5[[ii5,jj5]]≥Max@o5[[kk5,kk5]]&&Min@o5[[kk5,ll5]]≥Max@o5[[kk5,kk5]],If[Min@o5[[ii5,jj5]]>Max@o5[[kk5,kk5]]||Min@o5[[kk5,ll5]]>Max@o5[[kk5,kk5]]||MemberQ[{1,3,5},Max@o5[[kk5,kk5]]]==True,Return[{0}]]]];
Return[{1,0,-1}];
)

(* Compare r_{ij}'s or compare r_{ij} with 1 *)
CompareOrderr[rordermatrix5_,{ii5_,jj5_},{kk5_,ll5_}]:=(
ro5=rordermatrix5;

If[ii5≠jj5&&kk5≠ll5,
If[Max@ro5[[ii5,jj5]]<Min@ro5[[kk5,ll5]],Return[{-1}]];
If[Min@ro5[[ii5,jj5]]>Max@ro5[[kk5,ll5]],Return[{1}]];
If[Length[ro5[[ii5,jj5]]]==1&&ro5[[ii5,jj5]]==ro5[[kk5,ll5]]&&Complement[ro5[[ii5,jj5]],{1,3,5}]==0,Return[{0}]]];

If[ii5==jj5&&kk5≠ll5,
If[Min@ro5[[kk5,ll5]]≥4,Return[{-1}]];
If[Max@ro5[[kk5,ll5]]≤2,Return[{1}]];
If[ro5[[kk5,ll5]]=={3},Return[{0}] ]];
Return[{1,0,-1}];
)

(* Compare Z_{ij}'s according to Z_{ij}=z_{ij}^{-1/2} w_{ij}^{-3/2}=z_{ij}r_{ij}^{-3}=w_{ij}^{-1}r_{ij}^{-1}=W_{ij}^{1/3}w_{ij}^{-4/3} *)
CompareOrderZ1[{zomatrix4_,womatrix4_},{ii4_,jj4_},{kk4_,ll4_}]:=(zo4=zomatrix4;wo4=womatrix4;dis=Dis3[{zo4,wo4}];

If[ZX1[[ii4,jj4]]==1&&ZX1[[kk4,ll4]]==0,Return[{1}]];
If[ZX1[[ii4,jj4]]==1&&ZX1[[kk4,ll4]]==1,Return[{0}]];
If[ZX1[[ii4,jj4]]==0&&ZX1[[kk4,ll4]]==1,Return[{-1}]];

Comparez=CompareOrderzw[zo4,{ii4,jj4},{kk4,ll4}];
Comparew=CompareOrderzw[wo4,{ii4,jj4},{kk4,ll4}];
Comparer=CompareOrderr[dis,{ii4,jj4},{kk4,ll4}];

If[Comparez=={1}&&Comparew=={1},Return[{-1}]];
If[Comparez=={1}&&Comparew=={0},Return[{-1}]];
If[Comparez=={0}&&Comparew=={1},Return[{-1}]];
If[Comparez=={-1}&&Comparer=={1},Return[{-1}]];
If[Comparez=={0}&&Comparer=={1},Return[{-1}]];
If[Comparez=={-1}&&Comparer=={0},Return[{-1}]];
If[Comparew=={1}&&WX1[[kk4,ll4]]==1,Return[{-1}]];
If[Comparez=={0}&&Comparew=={0},Return[{0}]];
If[Comparez=={0}&&Comparer=={0},Return[{0}]];
If[Comparew=={0}&&WX1[[ii4,jj4]]==1&&WX1[[kk4,ll4]]==1,Return[{0}]];
If[Comparez=={-1}&&Comparew=={-1},Return[{1}]];
If[Comparez=={-1}&&Comparew=={0},Return[{1}]];
If[Comparez=={0}&&Comparew=={-1},Return[{1}]];
If[Comparez=={1}&&Comparer=={-1},Return[{1}]];
If[Comparez=={0}&&Comparer=={-1},Return[{1}]];
If[Comparez=={1}&&Comparer=={0},Return[{1}]];
If[Comparew=={-1}&&WX1[[ii4,jj4]]==1,Return[{1}]];
Return[{1,0,-1}];
)

(* Compare W_{ij}'s according to W_{ij}=w_{ij}^{-1/2} z_{ij}^{-3/2}=w_{ij}r_{ij}^{-3}=z_{ij}^{-1}r_{ij}^{-1}=Z_{ij}^{1/3}z_{ij}^{-4/3} *)
CompareOrderW1[{zomatrix4_,womatrix4_},{ii4_,jj4_},{kk4_,ll4_}]:=(zo4=zomatrix4;wo4=womatrix4;dis=Dis3[{zo4,wo4}];

If[WX1[[ii4,jj4]]==1&&WX1[[kk4,ll4]]==0,Return[{1}]];
If[WX1[[ii4,jj4]]==1&&WX1[[kk4,ll4]]==1,Return[{0}]];
If[WX1[[ii4,jj4]]==0&&WX1[[kk4,ll4]]==1,Return[{-1}]];

Comparez=CompareOrderzw[zo4,{ii4,jj4},{kk4,ll4}];
Comparew=CompareOrderzw[wo4,{ii4,jj4},{kk4,ll4}];
Comparer=CompareOrderr[dis,{ii4,jj4},{kk4,ll4}];

If[Comparez=={1}&&Comparew=={1},Return[{-1}]];
If[Comparez=={1}&&Comparew=={0},Return[{-1}]];
If[Comparez=={0}&&Comparew=={1},Return[{-1}]];
If[Comparew=={-1}&&Comparer=={1},Return[{-1}]];
If[Comparew=={0}&&Comparer=={1},Return[{-1}]];
If[Comparew=={-1}&&Comparer=={0},Return[{-1}]];
If[Comparez=={1}&&ZX1[[kk4,ll4]]==1,Return[{-1}]];
If[Comparez=={0}&&Comparew=={0},Return[{0}]];
If[Comparew=={0}&&Comparer=={0},Return[{0}]];
If[Comparez=={0}&&ZX1[[ii4,jj4]]==1&&ZX1[[kk4,ll4]]==1,Return[{0}]];
If[Comparez=={-1}&&Comparew=={-1},Return[{1}]];
If[Comparez=={-1}&&Comparew=={0},Return[{1}]];
If[Comparez=={0}&&Comparew=={-1},Return[{1}]];
If[Comparew=={1}&&Comparer=={-1},Return[{1}]];
If[Comparew=={0}&&Comparer=={-1},Return[{1}]];
If[Comparew=={1}&&Comparer=={0},Return[{1}]];
If[Comparez=={-1}&&ZX1[[ii4,jj4]]==1,Return[{1}]];
Return[{1,0,-1}];
)

(* Compare Z_{ij} and z_k *)
CompareOrderZ2[{zomatrix4_,womatrix4_},{ii4_,jj4_},{kk4_,ll4_}]:=(zo4=zomatrix4;wo4=womatrix4;dis=Dis3[{zo4,wo4}];

If[ZX1[[ii4,jj4]]==1&&ZX1[[kk4,ll4]]==0,Return[{1}]];
If[ZX1[[ii4,jj4]]==1&&ZX1[[kk4,ll4]]==1,Return[{0}]];
If[ZX1[[ii4,jj4]]==0&&ZX1[[kk4,ll4]]==1,Return[{-1}]];

Comparez=CompareOrderzw[zo4,{ii4,jj4},{kk4,ll4}];
Comparew=CompareOrderzw[wo4,{ii4,jj4},{kk4,ll4}];
Comparer=CompareOrderr[dis,{ii4,jj4},{kk4,ll4}];

If[Comparez=={-1}&&Comparer=={1},Return[{-1}]];
If[Comparez=={0}&&Comparer=={1},Return[{-1}]];
If[Comparez=={-1}&&Comparer=={0},Return[{-1}]];
If[Comparez=={0}&&Comparer=={0},Return[{0}]];
If[Comparez=={1}&&Comparer=={-1},Return[{1}]];
If[Comparez=={0}&&Comparer=={-1},Return[{1}]];
If[Comparez=={1}&&Comparer=={0},Return[{1}]];
Return[{-1,0,1}];
)

(* Compare W_{ij} and w_k *)
CompareOrderW2[{zomatrix4_,womatrix4_},{ii4_,jj4_},{kk4_,ll4_}]:=(zo4=zomatrix4;wo4=womatrix4;dis=Dis3[{zo4,wo4}];

If[WX1[[ii4,jj4]]==1&&WX1[[kk4,ll4]]==0,Return[{1}]];
If[WX1[[ii4,jj4]]==1&&WX1[[kk4,ll4]]==1,Return[{0}]];
If[WX1[[ii4,jj4]]==0&&WX1[[kk4,ll4]]==1,Return[{-1}]];

Comparez=CompareOrderzw[zo4,{ii4,jj4},{kk4,ll4}];
Comparew=CompareOrderzw[wo4,{ii4,jj4},{kk4,ll4}];
Comparer=CompareOrderr[dis,{ii4,jj4},{kk4,ll4}];

If[Comparew=={-1}&&Comparer=={1},Return[{-1}]];
If[Comparew=={0}&&Comparer=={1},Return[{-1}]];
If[Comparew=={-1}&&Comparer=={0},Return[{-1}]];
If[Comparew=={0}&&Comparer=={0},Return[{0}]];
If[Comparew=={1}&&Comparer=={-1},Return[{1}]];
If[Comparew=={0}&&Comparer=={-1},Return[{1}]];
If[Comparew=={1}&&Comparer=={0},Return[{1}]];
Return[{-1,0,1}];
)

(* Determine possible maximal terms in z_i-equation *)
MaxOrderZ[{zomatrix3_,womatrix3_},kk3_]:=(zo3=zomatrix3;wo3=womatrix3;max3=allind;s3=Complement[allind,{kk3}];ii3=1;While[ii3≤n-1,jj3=ii3+1;While[jj3≤n-1,
If[CompareOrderZ1[{zo3,wo3},{kk3,s3[[ii3]]},{kk3,s3[[jj3]]}]=={1},max3=Complement[max3,{s3[[jj3]]}];];
If[CompareOrderZ1[{zo3,wo3},{kk3,s3[[ii3]]},{kk3,s3[[jj3]]}]=={-1},max3=Complement[max3,{s3[[ii3]]}];];jj3++
];
If[CompareOrderZ2[{zo3,wo3},{kk3,kk3},{kk3,s3[[ii3]]}]=={1},max3=Complement[max3,{s3[[ii3]]}]];
If[CompareOrderZ2[{zo3,wo3},{kk3,kk3},{kk3,s3[[ii3]]}]=={-1},max3=Complement[max3,{kk3}];];
ii3++];Return[max3])

(* Determine possible maximal terms in w_i-equation *)
MaxOrderW[{zomatrix3_,womatrix3_},kk3_]:=(zo3=zomatrix3;wo3=womatrix3;max3=allind;s3=Complement[allind,{kk3}];ii3=1;While[ii3≤n-1,jj3=ii3+1;While[jj3≤n-1,
If[CompareOrderW1[{zo3,wo3},{kk3,s3[[ii3]]},{kk3,s3[[jj3]]}]=={1},max3=Complement[max3,{s3[[jj3]]}];];
If[CompareOrderW1[{zo3,wo3},{kk3,s3[[ii3]]},{kk3,s3[[jj3]]}]=={-1},max3=Complement[max3,{s3[[ii3]]}];];jj3++
];
If[CompareOrderW2[{zo3,wo3},{kk3,kk3},{kk3,s3[[ii3]]}]=={1},max3=Complement[max3,{s3[[ii3]]}]];
If[CompareOrderW2[{zo3,wo3},{kk3,kk3},{kk3,s3[[ii3]]}]=={-1},max3=Complement[max3,{kk3}];];
ii3++];Return[max3])

(* Proposition 4.12 *)
Rule1cQ[{zomatrix2_,womatrix2_}]:=(
ii2=1;While[ii2≤Length[zc1],zc2=Intersection[zc1[[ii2]],Flatten@Position[Diagonal[ZX1],1]];If[Length[zc2]>0&&Count[Flatten@zomatrix2[[zc2,zc2]],5]==Length[zc2],Return[False]];ii2++];ii2=1;While[ii2≤Length[wc1],wc2=Intersection[wc1[[ii2]],Flatten@Position[Diagonal[WX1],1]];If[Length[wc2]>0&&Count[Flatten@womatrix2[[wc2,wc2]],5]==Length[wc2],Return[False]];ii2++];Return[True])

(* Proposition 4.13 *)
OrderQ1[{zomatrix2_,womatrix2_}]:=(kk2=1;While[kk2≤n,
If[Length[MaxOrderZ[{zomatrix2,womatrix2},kk2]]==1||Length[MaxOrderW[{zomatrix2,womatrix2},kk2]]==1,Return[False]];kk2++;
];Return[True])

(* Proposition 4.14 *)
OrderQ2[{zomatrix2_,womatrix2_}]:=(dis2=Dis3[{zomatrix2,womatrix2}];ii2=1;While[ii2≤n-1,jj2=ii2+1;While[jj2≤n,
If[MaxOrderZ[{zomatrix2,womatrix2},ii2]=={ii2,jj2}&&MaxOrderZ[{zomatrix2,womatrix2},jj2]=={ii2,jj2}&&Intersection[dis2[[ii2,jj2]],{3}]=={},Return[False]];If[MaxOrderW[{zomatrix2,womatrix2},ii2]=={ii2,jj2}&&MaxOrderW[{zomatrix2,womatrix2},jj2]=={ii2,jj2}&&Intersection[dis2[[ii2,jj2]],{3}]=={},Return[False]];
jj2++;
];ii2++];Return[True])

(* Proposition 4.15 *)
OrderQ3[{zomatrix2_,womatrix2_}]:=(
maxZ=Table[MaxOrderZ[{zomatrix2,womatrix2},i2],{i2,n}];
maxW=Table[MaxOrderW[{zomatrix2,womatrix2},i2],{i2,n}];
ii2=1;
While[ii2≤n-2,jj2=ii2+1;
While[jj2≤n-1,kk2=jj2+1;
While[kk2≤n,

If[ZX1[[ii2,jj2]]ZX1[[jj2,kk2]]≠0 && ZX1[[kk2,ii2]]==0,
If[maxW[[jj2]]=={ii2,kk2},Return[False]];
If[Complement[maxW[[ii2]]⋃maxW[[jj2]]⋃maxW[[kk2]],{ii2,jj2,kk2}]=={},Return[False]];
];
If[WX1[[ii2,jj2]]WX1[[jj2,kk2]]≠0 && WX1[[kk2,ii2]]==0,
If[maxZ[[jj2]]=={ii2,kk2},Return[False]];
If[Complement[maxZ[[ii2]]⋃maxZ[[jj2]]⋃maxZ[[kk2]],{ii2,jj2,kk2}]=={},Return[False]];
];
If[ZX1[[jj2,kk2]]ZX1[[kk2,ii2]]≠0 && ZX1[[ii2,jj2]]==0,
If[maxW[[kk2]]=={ii2,jj2},Return[False]];
If[Complement[maxW[[ii2]]⋃maxW[[jj2]]⋃maxW[[kk2]],{ii2,jj2,kk2}]=={},Return[False]];
];
If[WX1[[jj2,kk2]]WX1[[kk2,ii2]]≠0 && WX1[[ii2,jj2]]==0,
If[maxZ[[kk2]]=={ii2,jj2},Return[False]];
If[Complement[maxZ[[ii2]]⋃maxZ[[jj2]]⋃maxZ[[kk2]],{ii2,jj2,kk2}]=={},Return[False]];

];
If[ZX1[[kk2,ii2]]ZX1[[ii2,jj2]]≠0 && ZX1[[jj2,kk2]]==0,
If[maxW[[ii2]]=={jj2,kk2},Return[False]];
If[Complement[maxW[[ii2]]⋃maxW[[jj2]]⋃maxW[[kk2]],{ii2,jj2,kk2}]=={},Return[False]];
];
If[WX1[[kk2,ii2]]WX1[[ii2,jj2]]≠0 && WX1[[jj2,kk2]]==0,
If[maxZ[[ii2]]=={jj2,kk2},Return[False]];
If[Complement[maxZ[[ii2]]⋃maxZ[[jj2]]⋃maxZ[[kk2]],{ii2,jj2,kk2}]=={},Return[False]];

];
kk2++];
jj2++];
ii2++];Return[True])

outputform1[ordermatrix2_]:=(OMA=ordermatrix2;OMB=Table[If[Length@OMA[[i2,j2]]>1,{ToString[Min@OMA[[i2,j2]]]<>"~"<>ToString[Max@OMA[[i2,j2]]]},OMA[[i2,j2]]],{i2,1,n},{j2,1,n}];Return[OMB])

(* --------------------- MAIN PROGRAM FOR FINDING ORDERS ----------------------------- *)
FindOrder[x_]:=FindOrder[x,60];
FindOrder[x_,time1_]:=(
t1=AbsoluteTime[];
N1=x;X1=FMX[[x]];Clear[zwofamily,zwotype2,zwotype3]; constrainttime=time1;TF1=0;Print[Style["finding zw-order matrices covering ",Bold,16],Style["FMX[["<>ToString[N1]<>"]]",Bold,16,RGBColor[0.8,0.4,0]],Style["...",Bold,16]];
Print[];

zo1=Table[{1,2,3,4,5},{i1,1,n},{j1,1,n}];
wo1=Table[{1,2,3,4,5},{i1,1,n},{j1,1,n}];
ZX1=zpart[X1];WX1=wpart[X1];XX1=ZX1+2*WX1;
zc1=componentsZ[X1];wc1=componentsW[X1];pmzw={};Do[If[ZX1[[pm[[i1]],pm[[i1]]]]==ZX1&&WX1[[pm[[i1]],pm[[i1]]]]==WX1,pmzw=pmzw~Join~{pm[[i1]]}],{i1,1,Length[pm]}];

(* Proposition 4.4(a) *)
Do[If[ZX1[[i1,i1]]==1,zo1[[i1,i1]]={5},zo1[[i1,i1]]={0,1,2,3,4}];If[WX1[[i1,i1]]==1,wo1[[i1,i1]]={5},wo1[[i1,i1]]={0,1,2,3,4}],{i1,1,n}];

(* Proposition 4.8 *)
Do[L1=Length[zc1[[i1]]];Pos1=Position[XX1[[zc1[[i1]],zc1[[i1]]]],1];Pos1=Complement[Pos1,Table[{j1,j1},{j1,1,n}]];If[Length[Pos1]==2&&Total[Diagonal[ZX1[[zc1[[i1]],zc1[[i1]]]]]]>0,Do[zo1[[zc1[[i1]][[Pos1[[j1]][[1]]]],zc1[[i1]][[Pos1[[j1]][[2]]]]]]={5};wo1[[zc1[[i1]][[Pos1[[j1]][[1]]]],zc1[[i1]][[Pos1[[j1]][[2]]]]]]={1},{j1,1,2}]],{i1,1,Length[zc1]}];
Do[L1=Length[wc1[[i1]]];Pos1=Position[XX1[[wc1[[i1]],wc1[[i1]]]],2];Pos1=Complement[Pos1,Table[{j1,j1},{j1,1,n}]];If[Length[Pos1]==2&&Total[Diagonal[WX1[[wc1[[i1]],wc1[[i1]]]]]]>0,Do[wo1[[wc1[[i1]][[Pos1[[j1]][[1]]]],wc1[[i1]][[Pos1[[j1]][[2]]]]]]={5};zo1[[wc1[[i1]][[Pos1[[j1]][[1]]]],wc1[[i1]][[Pos1[[j1]][[2]]]]]]={1},{j1,1,2}]],{i1,1,Length[wc1]}];

(* Proposition 4.4(b) *)
Do[If[i1≠j1,If[XX1[[i1,j1]]==1,zo1[[i1,j1]]=Intersection[zo1[[i1,j1]],{4,5}];wo1[[i1,j1]]=Intersection[wo1[[i1,j1]],{1,2}]]],{i1,1,n},{j1,1,n}];
Do[If[i1≠j1,If[XX1[[i1,j1]]==2,wo1[[i1,j1]]=Intersection[wo1[[i1,j1]],{4,5}];zo1[[i1,j1]]=Intersection[zo1[[i1,j1]],{1,2}]]],{i1,1,n},{j1,1,n}];
Do[If[i1≠j1,If[XX1[[i1,j1]]==3,wo1[[i1,j1]]=Intersection[wo1[[i1,j1]],{3}];zo1[[i1,j1]]=Intersection[zo1[[i1,j1]],{3}]]],{i1,1,n},{j1,1,n}];

(* Determined the order of z_{ij} and w_{ij} if there is a maximal z-or w-stroke *)
Do[If[i1≠j1,If[XX1[[i1,j1]]==1&&(Total[ZX1[[i1]]]+Total[ZX1[[j1]]]-ZX1[[i1,i1]]-ZX1[[j1,j1]])==2,zo1[[i1,j1]]=Intersection[zo1[[i1,j1]],{5}];wo1[[i1,j1]]=Intersection[wo1[[i1,j1]],{1}]]],{i1,1,n},{j1,1,n}];
Do[If[i1≠j1,If[XX1[[i1,j1]]==2&&(Total[WX1[[i1]]]+Total[WX1[[j1]]]-WX1[[i1,i1]]-WX1[[j1,j1]])==2,zo1[[i1,j1]]=Intersection[zo1[[i1,j1]],{1}];wo1[[i1,j1]]=Intersection[wo1[[i1,j1]],{5}]]],{i1,1,n},{j1,1,n}];

(* Proposition 4.4(b) *)
Do[If[i1≠j1,If[ZX1[[i1,j1]]==0,wo1[[i1,j1]]=Intersection[wo1[[i1,j1]],{2,3,4,5}]]],{i1,1,n},{j1,1,n}];
Do[If[i1≠j1,If[WX1[[i1,j1]]==0,zo1[[i1,j1]]=Intersection[zo1[[i1,j1]],{2,3,4,5}]]],{i1,1,n},{j1,1,n}];

(* Refine $zw$-order matrices by principles for updating orders *)
{zo1,wo1}=UpdateOrder[{zo1,wo1}];

Print[Style["result of setp 1:",Bold,15]];
Print["cumulative time: ",AbsoluteTime[]-t1 ];Print[];
Print["{z-order, w-order, r-order} matrices :"];
Print[{Style[MatrixForm[outputform1@zo1],Blue],Style[MatrixForm[outputform1@wo1],Blue],Style[MatrixForm[outputform1@Dis3[{zo1,wo1}]],RGBColor[0,0.5,0]]}];
Str="There are "<>ToString[NumberofType2[{zo1,wo1}]]<>" possible zw-order matrices of Type 2.";
Print[];Print[Str];
Str="There are "<>ToString[NumberofType3[{zo1,wo1}]]<>" possible zw-order matrices of Type 3.";
Print[Str];Print[];Print[];

(* ------------------------ Step 2:Type 2 order matrices --------------------------- *)
zwotype2={};zwotype3={};
TimeConstrained[zwofamily={{zo1,wo1}};
While[Length[zwofamily]>0,
L1=Length[zwofamily];Lz1=0;Lw1=0;Del={};
Do[If[NumberofType2[zwofamily[[i1]]]>1,
Del=Del~Join~{{i1}};
{Pzi,Pzj,Lz1}=PosOffDiag[zwofamily[[i1]][[1]]];
If[Lz1==0,
{Pwi,Pwj,Lw1}=PosOffDiag[zwofamily[[i1]][[2]]];
Do[
zwofamily=zwofamily~Join~{zwofamily[[i1]]};
zwofamily[[-1]][[2]][[Pwi,Pwj]]={zwofamily[[i1]][[2]][[Pwi,Pwj]][[j1]]};
zwofamily[[-1]][[2]][[Pwj,Pwi]]={zwofamily[[i1]][[2]][[Pwi,Pwj]][[j1]]};
zwofamily[[-1]]=UpdateOrder[zwofamily[[-1]]];
,{j1,1,Lw1}]
,Do[
zwofamily=zwofamily~Join~{zwofamily[[i1]]};
zwofamily[[-1]][[1]][[Pzi,Pzj]]={zwofamily[[i1]][[1]][[Pzi,Pzj]][[j1]]};
zwofamily[[-1]][[1]][[Pzj,Pzi]]={zwofamily[[i1]][[1]][[Pzi,Pzj]][[j1]]};
zwofamily[[-1]]=UpdateOrder[zwofamily[[-1]]];
,{j1,1,Lz1}]
]
],{i1,1,L1}];
zwofamily=Delete[zwofamily,Del];
zwotype2=zwotype2~Join~Select[zwofamily,NumberofType2[#]==1&];
zwofamily=Select[zwofamily,NumberofType2[#]>1& ];
];
(* End of While loop *)

zwotype2=Select[zwotype2,Rule2hQ];
zwotype2=Select[zwotype2,OrderQ1];
zwotype2=Select[zwotype2,OrderQ2];
zwotype2=Select[zwotype2,OrderQ3];
zwotype2=Select[zwotype2,Rule1cQ];

del={};ii1=2;While[ii1≤Length[zwotype2],jj1=1;tf1=0;While[jj1<ii1&&tf1==0,kk1=1;While[kk1≤Length[pmzw],If[zwotype2[[ii1]][[1]]==zwotype2[[jj1]][[1]][[pmzw[[kk1]],pmzw[[kk1]]]]&&zwotype2[[ii1]][[2]]==zwotype2[[jj1]][[2]][[pmzw[[kk1]],pmzw[[kk1]]]],del=del~Join~{{ii1}};tf1=1];kk1++];jj1++];ii1++];zwotype2=Delete[zwotype2,del];

Print[];
Print[Style["result of setp 2:",Bold,15]];
Print["cumulative time: ",AbsoluteTime[]-t1 ];Print[];
Print["There are ",Style[Length[zwotype2],Bold,RGBColor[0,0,0.7],14]," zw-order matrices of Type 2."];
Print["{z-order, w-order, r-order} matrices :"];Print[Table[{Style[MatrixForm[outputform1@zwotype2[[j1,1]]],Blue],Style[MatrixForm[outputform1@zwotype2[[j1,2]]],Blue],Style[MatrixForm[outputform1@Dis3[zwotype2[[j1]]]],RGBColor[0,0.5,0]]},{j1,1,Length[zwotype2]}]];
Str="There are "<>ToString[Sum[NumberofType3[zwotype2[[i1]]],{i1,1,Length[zwotype2]}]]<>" possible zw-order matrices of Type 3.";
Print[];Print[Str];Print[];Print[];

(* ----------------------- Step 3:Type 3 order matrices --------------------------- *)

zwofamily=zwotype2;
While[Length[zwofamily]>0,
L1=Length[zwofamily];Lz1=0;Lw1=0;Del={};
Do[
If[NumberofType3[zwofamily[[i1]]]>1,
Del=Del~Join~{{i1}};
{Pzd,Lz1}=PosDiag[zwofamily[[i1]][[1]]];
If[Lz1==0,
{Pwd,Lw1}=PosDiag[zwofamily[[i1]][[2]]];
Do[
zwofamily=zwofamily~Join~{zwofamily[[i1]]};
zwofamily[[-1]][[2]][[Pwd,Pwd]]={zwofamily[[i1]][[2]][[Pwd,Pwd]][[j1]]};
zwofamily[[-1]]=UpdateOrderPart[zwofamily[[-1]]];
,{j1,1,Lw1}],
Do[zwofamily=zwofamily~Join~{zwofamily[[i1]]};
zwofamily[[-1]][[1]][[Pzd,Pzd]]={zwofamily[[i1]][[1]][[Pzd,Pzd]][[j1]]};
zwofamily[[-1]]=UpdateOrderPart[zwofamily[[-1]]];
,{j1,1,Lz1}]
]],{i1,1,L1}];
zwofamily=Delete[zwofamily,Del];
zwotype3=zwotype3~Join~Select[zwofamily,NumberofType3[#]==1&];
zwofamily=Select[zwofamily,NumberofType3[#]>1&];];

del={};ii1=2;While[ii1≤Length[zwotype3],jj1=1;tf1=0;While[jj1<ii1&&tf1==0,kk1=1;While[kk1≤Length[pmzw],If[zwotype3[[ii1]][[1]]==zwotype3[[jj1]][[1]][[pmzw[[kk1]],pmzw[[kk1]]]]&&zwotype3[[ii1]][[2]]==zwotype3[[jj1]][[2]][[pmzw[[kk1]],pmzw[[kk1]]]],del=del~Join~{{ii1}};tf1=1];kk1++];jj1++];ii1++];zwotype3=Delete[zwotype3,del];

zwotype3=Select[zwotype3,OrderQ1];
zwotype3=Select[zwotype3,OrderQ2];
zwotype3=Select[zwotype3,Rule1cQ];

l1=Length[zwotype3];
Print[Style["result of setp 3:",Bold,15]];
Print["cumulative time: ",AbsoluteTime[]-t1 ];Print[];
Print["There are ",Style[l1,Bold,RGBColor[0,0,0.7],14]," zw-order matrices of Type 3."];

(* -------------------------- Output results of FindOrder ----------------------------- *)

If[l1≤10,Print["{z-order, w-order, r-order} matrices :"];Print[Table[{Style[MatrixForm[outputform1@zwotype3[[j1,1]]],Blue],Style[MatrixForm[outputform1@zwotype3[[j1,2]]],Blue],Style[MatrixForm[outputform1@Dis3[zwotype3[[j1]]]],RGBColor[0,0.5,0]]},{j1,1,Length[zwotype3]}]]];

If[Length[zwotype3]==0,Print[];Str="matrix FMX[["<>ToString[N1]<>"]] is impossible"; Print[Str]];TF1=1;Print[];Print[];
, constrainttime];

If[TF1==0,Str="stop running FMX[["<>ToString[N1]<>"]]";Print[Str];Print[];Print[];])

(* ----------------------------------- End of FindOrder ---------------------------------- *)

FindOrder[21]

(* Algorithm III: Find Mass Relations. Define functions after implementing "FindOrder[k]", then input "FindMassRel1[...]", "FindMassRel2[...]", or "FindMassRel3[...]".  Details about the inputs of these functions are given below.  *)

(*Set variables*)
Z=Table[ToExpression["Z"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
W=Table[ToExpression["W"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
u=Table[ToExpression["u"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
v=Table[ToExpression["v"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
zz=Table[ToExpression["z"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
ww=Table[ToExpression["w"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
rr=Table[ToExpression["r"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
delta=Table[ToExpression["δ"<>ToString[j1]<>ToString[k1]],{j1,1,n},{k1,1,n}];
z=Table[ToExpression["z"<>ToString[j1]],{j1,1,n}];
w=Table[ToExpression["w"<>ToString[j1]],{j1,1,n}];
zeta=Table[ToExpression["ζ"<>ToString[j1]],{j1,1,n}];
omega=Table[ToExpression["ω"<>ToString[j1]],{j1,1,n}];
M=Table[ToExpression["m"<>ToString[j1]],{j1,1,n}];
mu=Table[ToExpression["μ"<>ToString[j1]],{j1,1,n}];
Vardis=Flatten[Join[zz,ww,delta]];
Var=Flatten[Join[zz,ww,delta,z,w,M]];
zpart[{zmx4_,wmx4_}]:=zmx4;
wpart[{zmx4_,wmx4_}]:=wmx4;

(* Determine whether the possible max. order terms in a (Subscript[z, i])-equation or a (Subscript[w, i])-equation are indeed max. order *)
CertainMaxOrderZQ[{zomatrix3_,womatrix3_},k3_]:=(zo3=zomatrix3;wo3=womatrix3;
max3=MaxOrderZ[{zo3,wo3},k3];
If[Length[max3]==2,Return[True]];
If[Length[max3]>2,
Max3={};
Do[
Do[If[max3[[i3]]==k3&&max3[[i3]]≠max3[[j3]],
If[CompareOrderZ2[{zo3,wo3},{k3,max3[[i3]]},{k3,max3[[j3]]}]=={0},Max3=Max3~Join~{{max3[[i3]],max3[[j3]]}}];
];
If[max3[[j3]]==k3&&max3[[i3]]≠max3[[j3]],
If[CompareOrderZ2[{zo3,wo3},{k3,max3[[j3]]},{k3,max3[[i3]]}]=={0},Max3=Max3~Join~{{max3[[i3]],max3[[j3]]}}];
];
If[max3[[i3]]!=k3&&max3[[j3]]!=k3&&max3[[i3]]≠max3[[j3]],
If[CompareOrderZ1[{zo3,wo3},{k3,max3[[j3]]},{k3,max3[[i3]]}]=={0},Max3=Max3~Join~{{max3[[i3]],max3[[j3]]}}];
];
,{j3,i3+1,Length[max3]}];
,{i3,1,Length[max3]-1}];

If[Length[Max3]>1,
Do[Do[Do[If[Intersection[Max3[[i3]],Max3[[j3]]]≠{},Max3[[i3]]=Max3[[i3]]⋃Max3[[j3]];
Max3[[j3]]={};];
,{j3,i3+1,Length[Max3]}]
,{i3,1,Length[Max3]}]
,{ii3,1,Length[Max3]}];];
Max3=Complement[Max3,{{}}];

If[Length[Max3]==1&&Flatten[Max3]==max3,Return[True],Return[False]]

;]
)

CertainMaxOrderWQ[{zomatrix3_,womatrix3_},k3_]:=(zo3=zomatrix3;wo3=womatrix3;
max3=MaxOrderW[{zo3,wo3},k3];
If[Length[max3]==2,Return[True]];
If[Length[max3]>2,
Max3={};
Do[
Do[If[max3[[i3]]==k3&&max3[[i3]]≠max3[[j3]],
If[CompareOrderW2[{zo3,wo3},{k3,max3[[i3]]},{k3,max3[[j3]]}]=={0},Max3=Max3~Join~{{max3[[i3]],max3[[j3]]}}];
];
If[max3[[j3]]==k3&&max3[[i3]]≠max3[[j3]],
If[CompareOrderW2[{zo3,wo3},{k3,max3[[j3]]},{k3,max3[[i3]]}]=={0},Max3=Max3~Join~{{max3[[i3]],max3[[j3]]}}];
];
If[max3[[i3]]!=k3&&max3[[j3]]!=k3&&max3[[i3]]≠max3[[j3]],
If[CompareOrderW1[{zo3,wo3},{k3,max3[[j3]]},{k3,max3[[i3]]}]=={0},Max3=Max3~Join~{{max3[[i3]],max3[[j3]]}}];
];
,{j3,i3+1,Length[max3]}];
,{i3,1,Length[max3]-1}];

If[Length[Max3]>1,
Do[Do[Do[If[Intersection[Max3[[i3]],Max3[[j3]]]≠{},Max3[[i3]]=Max3[[i3]]⋃Max3[[j3]];
Max3[[j3]]={};];
,{j3,i3+1,Length[Max3]}]
,{i3,1,Length[Max3]}]
,{ii3,1,Length[Max3]}];];
Max3=Complement[Max3,{{}}];

If[Length[Max3]==1&&Flatten[Max3]==max3,Return[True],Return[False]]

;]
)

(* Determine the (Subscript[z, i])-equations corresponding to zw-order matrices for each i *)
EqnMZO[{zomatrix2_,womatrix2_},vertexSet2_]:=(
zo2=zomatrix2;wo2=womatrix2;eqnset2={};
Do[If[CertainMaxOrderZQ[{zo2,wo2},i2]==True,eqn2=0;max2=MaxOrderZ[{zo2,wo2},i2];Do[If[MemberQ[max2,j2]==True,eqn2=eqn2+M[[j2]]*Z[[i2,j2]]],{j2,1,n}];
eqnset2=eqnset2⋃{eqn2};]
,{i2,vertexSet2}];Return[eqnset2];
)

(* Determine the (Subscript[w, i])-equations corresponding to zw-order matrices for each i *)
EqnMWO[{zomatrix2_,womatrix2_},vertexSet2_]:=(
zo2=zomatrix2;wo2=womatrix2;eqnset2={};
Do[If[CertainMaxOrderWQ[{zo2,wo2},i2]==True,eqn2=0;max2=MaxOrderW[{zo2,wo2},i2];Do[If[MemberQ[max2,j2]==True,eqn2=eqn2+M[[j2]]*W[[i2,j2]]],{j2,1,n}];
eqnset2=eqnset2⋃{eqn2};]
,{i2,vertexSet2}];Return[eqnset2];
)

(* Determine the cluster equations according zw-order matrices *)
ClusterZ[{zomatrix2_,womatrix2_}]:=(VV2={};Do[Do[If[Min[zomatrix2[[i2,i2]]~Join~zomatrix2[[j2,j2]]]>Max[zomatrix2[[i2,j2]]],VV2=VV2~Join~{z[[i2]]-z[[j2]]}],{i2,j2+1,n}],{j2,1,n-1}];Return[VV2];)
ClusterW[{zomatrix2_,womatrix2_}]:=(VV2={};Do[Do[If[Min[womatrix2[[i2,i2]]~Join~womatrix2[[j2,j2]]]>Max[womatrix2[[i2,j2]]],VV2=VV2~Join~{w[[i2]]-w[[j2]]}],{i2,j2+1,n}],{j2,1,n-1}];Return[VV2];)

(* Define functions for searching mass relations *)

(* Replace (Subscript[m, i])(Subscript[Z, ii])'s , (Subscript[m, i])(Subscript[W, ii])'s by (Subscript[z, i]), (Subscript[w, i]) resp.  *)
Relation1[x2_]:=(y2=x2;Do[Do[y2[[j2]]=y2[[j2]]/.{M[[i2]]*Z[[i2,i2]]->z[[i2]],M[[i2]]*W[[i2,i2]]->w[[i2]]},{i2,1,n}],{j2,1,Length[y2]}];Return[y2])

(* For each pair of (i,j) with i>j, replace (Subscript[Z, ij]) , (Subscript[W, ij]) by -(Subscript[Z, ji]) , -(Subscript[W, ji]) resp.  *)
Relation2[x2_]:=(y2=x2;Do[Do[Do[y2[[k2]]=y2[[k2]]/.{Z[[i2,j2]]->-Z[[j2,i2]],W[[i2,j2]]->-W[[j2,i2]]},{i2,j2+1,n}],{j2,1,n-1}],{k2,1,Length[y2]}];Return[y2])

(* For every i,j, replace (Subscript[Z, ij]) , (Subscript[W, ij]) by ((Subscript[z, ij])^(-1/2))((Subscript[w, ij])^(-3/2)) , ((Subscript[w, ij])^(-1/2))((Subscript[z, ij])^(-3/2)) resp.  *)
Relation3ver1[x2_]:=(y2=x2;
Do[Do[Do[y2[[k2]]=y2[[k2]]/.{Z[[i2,j2]]->zz[[i2,j2]]^(-1/2)*ww[[i2,j2]]^(-3/2),W[[i2,j2]]->zz[[i2,j2]]^(-1/2)*ww[[i2,j2]]^(-3/2)},{j2,i2+1,n}],{i2,1,n-1}],{k2,1,Length[y2]}];Return[y2])

(* For every i,j, replace (Subscript[Z, ij]) , (Subscript[W, ij]) by (Subscript[u, ij])((Subscript[v, ij])^3) , (Subscript[v, ij])((Subscript[u, ij])^3) resp.  *)
Relation3ver2[x2_]:=(y2=x2;
Do[Do[Do[y2[[k2]]=y2[[k2]]/.{Z[[i2,j2]]->u[[i2,j2]]^(1)*v[[i2,j2]]^(3),W[[i2,j2]]->v[[i2,j2]]^(1)*u[[i2,j2]]^(3)},{j2,i2+1,n}],{i2,1,n-1}],{k2,1,Length[y2]}];Return[y2])

(* Change variable according to the setting (Subscript[z, i])=((Subscript[ζ, i])^2) , (Subscript[w, i])=((Subscript[ω, i])^2)  *)
ChangeVar1[x2_]:=(y2=x2;Do[Do[y2[[j2]]=y2[[j2]]/.{z[[i2]]->zeta[[i2]]^(2),w[[i2]]->omega[[i2]]^(2)},{i2,1,n}],{j2,1,Length[y2]}];Return[y2])
ChangeVar1Inv[x2_]:=(y2=x2;Do[Do[y2[[j2]]=y2[[j2]]/.{zeta[[i2]]->z[[i2]]^(1/2),omega[[i2]]->w[[i2]]^(1/2)},{i2,1,n}],{j2,1,Length[y2]}];Return[y2])

(* Change variable according to the setting (Subscript[m, i])=((Subscript[μ, i])^4)  *)
ChangeVar2[x2_]:=(y2=x2;Do[Do[y2[[j2]]=y2[[j2]]/.{M[[i2]]->mu[[i2]]^4},{i2,1,n}],{j2,1,Length[y2]}];Return[y2])
ChangeVar2Inv[x2_]:=(y2=x2;Do[Do[y2[[j2]]=y2[[j2]]/.{mu[[i2]]->M[[i2]]^(1/4)},{i2,1,n}],{j2,1,Length[y2]}];Return[y2])

(* Change variable according to the setting (Subscript[u, ij])=((Subscript[z, ij])^(-1/2)) ,  (Subscript[v, ij])=((Subscript[w, ij])^(-1/2)) *)
ChangeVar3[x2_]:=(y2=x2;Do[Do[Do[y2[[k2]]=Simplify@Expand[y2[[k2]]/.{zz[[i2,j2]]->u[[i2,j2]]^(-2),ww[[i2,j2]]->v[[i2,j2]]^(-2)}],{j2,i2+1,n}],{i2,1,n-1}],{k2,1,Length[y2]}];Return[y2])
ChangeVar3Inv[x2_]:=(y2=x2;Do[Do[Do[y2[[k2]]=Simplify@Expand[y2[[k2]]/.{u[[i2,j2]]->zz[[i2,j2]]^(-1/2),v[[i2,j2]]->ww[[i2,j2]]^(-1/2)}],{j2,i2+1,n}],{i2,1,n-1}],{k2,1,Length[y2]}];Return[y2])

(* If (Subscript[z, i])≻(Subscript[z, ij])  and  (Subscript[z, j])≻(Subscript[z, ij]), then (Subscript[z, i])~(Subscript[z, j]) , therefore we can replace (Subscript[z, j]) by (Subscript[z, i]) , similar for (Subscript[w, i])'s *)
ClusterReplace[{zomatrix2_,womatrix2_},x2_]:=(y2=x2;Do[Do[Do[If[Min[zomatrix2[[i2,i2]]~Join~zomatrix2[[j2,j2]]]>Max[zomatrix2[[i2,j2]]],y2[[k2]]=Simplify[y2[[k2]]/.{z[[i2]]->z[[j2]]}]];
If[Min[womatrix2[[i2,i2]]~Join~womatrix2[[j2,j2]]]>Max[womatrix2[[i2,j2]]],y2[[k2]]=y2[[k2]]/.{w[[i2]]->w[[j2]]}],{i2,j2+1,n}],{j2,1,n-1}],{k2,1,Length[y2]}];Return[y2])

 (* If a polynomial p has a factor ((Subscript[var, i])^k) for some variable (Subscript[var, i]) , then we change p to p/((Subscript[var, i])^k) *)
CancelFactor[x4_]:=(y4=x4;
Do[varset4=Variables[y4[[k4]]];Do[Ex4=Exponent[y4[[k4]],varset4[[j4]]^(-1)];If[Ex4<0,y4[[k4]]=Simplify[Expand[y4[[k4]]*varset4[[j4]]^(Ex4)]]
],{j4,1,Length[varset4]}],{k4,1,Length[y4]}];Return[y4])

(* If a rational function f=q/p and p has a factor ((Subscript[var, i])^k) for some variable (Subscript[var, i]) , then we change g to g*((Subscript[var, i])^k) *)
CancelDenominator[x3_]:=(y3=x3;
Do[varset3=Variables[y3[[k3]]];Do[Ex3=Exponent[y3[[k3]],varset3[[j3]]^(-1)];If[Ex3>0,y3[[k3]]=Simplify[Expand[y3[[k3]]*varset3[[j3]]^(Ex3)]]
],{j3,1,Length[varset3]}],{k3,1,Length[y3]}];Return[y3])

(* Determine whether the polynomial gives a mass relation. If it does, then we collect this mass relation. *)
MassRelationQ[x2_]:=(If[Count[Table[Complement[Variables[x2[[k2]]],M~Join~mu],{k2,1,Length[x2]}],{}]>0,Return[True],Return[False]])
TakeMassRela[x2_]:=(y2=x2;y2=Select[y2,Complement[Variables[#],M~Join~mu]=={}&];Return[y2])

(* Determine whether the pair of zw-order matrices gives any wedge equations. If it does, then we collect these wedge equations *)
FindWedgeV[{I2_,J2_}]:=(SS2=Select[Table[k,{k,1,n}],Max[Diagonal[I2][[#]],Diagonal[J2][[#]]]==5&];S2=Select[Subsets[SS2],Length[#]>0&&Length[#]<n&];jj2=1;While[jj2≤Length[S2],V2=S2[[jj2]];V2C=Complement[SS2,V2];U2=Complement[allind,V2~Join~V2C];If[Complement[Flatten[I2[[V2,V2]]],{5}]=={}&&Max[Flatten[I2[[V2C,V2C]]]]<5&&Complement[Flatten[I2[[V2,V2C]]],{5}]=={}&&Max[Flatten[J2[[V2,V2]]]]<5&&Complement[Flatten[J2[[V2C,V2C]]],{5}]=={}&&Complement[Flatten[J2[[V2,V2C]]],{5}]=={}&&Length@U2≤1&&Max@(Flatten@(I2[[U2,U2]])~Join~Flatten@(J2[[U2,U2]]))<=1,Return[{V2,V2C}]];jj2++];Return[{{},{}}])

(* Determine how many terms invoving given variables are there in the polynomial *)
NumberofInv[S4_,var4_]:=(ii4=0;Do[If[Intersection[Variables[S4[[j4]]],{var4}]≠{},ii4++],{j4,1,Length[S4]}];Return[ii4])

(* For a collection {(Subscript[f, i])} of polynomials, if there is some (Subscript[f, i])=g*var, where var is a variable and g is a polynomial independent of var, then we make substitution f/g for var to each (Subscript[f, j]) in the collection *)
simpleSubstitution[S3_,varset3a_]:=(s3=S3;varset3b=Select[varset3a,NumberofInv[s3,#]>1& ];ii3=1;tf2=0;While[ii3≤Length[s3]&&tf2==0,If[Length[ExpandAll@(s3[[ii3]])]==2,varset3c=Intersection[Variables[s3[[ii3]]],varset3b];jj3=1;While[jj3≤Length[varset3c]&&tf2==0,If[Exponent[s3[[ii3]],varset3c[[jj3]]]==1,substi3=Solve[s3[[ii3]]==0,varset3c[[jj3]]][[1]];Do[s3[[k3]]=Simplify@Expand[s3[[k3]]/.substi3],{k3,1,Length[s3]}];tf2=1];jj3++]];ii3++];Return[s3];)

(* Delete duplicates from a collection of polynomials of some variables *)
DelDupPoly[polyset3a_,varset3_]:=(polyset3b=polyset3a;
del={};
Do[If[Intersection[Variables[polyset3b[[i3]]],varset3]≠{},var3=Intersection[Variables[polyset3b[[i3]]],varset3][[1]];exp3=Exponent[polyset3b[[i3]],var3];
Do[If[Simplify[polyset3b[[i3]]*Coefficient[polyset3b[[j3]],var3^exp3]-polyset3b[[j3]]*Coefficient[polyset3b[[i3]],var3^exp3]]==0,del=del~Join~{{j3}}],{j3,i3+1,Length[polyset3b]}]],{i3,1,Length[polyset3b]-1}];
polyset3b=Delete[polyset3b,del];Return[polyset3b];
)

(* For a collection {(Subscript[f, i])} of polynomails and a set V of variables, if (Subscript[f, i])=(Subscript[f, j])*g and i<j, where g is independent of variables in V, then we delete (Subscript[f, j]) from the collection *)
DelDupPoly2[polyset3a_,varset3_]:=(polyset3b=polyset3a;
del={};
Do[
Do[If[PolynomialQ[Simplify@(Expand@(polyset3b[[i3]]/polyset3b[[j3]])),varset3]==True&&PolynomialQ[Simplify@(Expand@(polyset3b[[j3]]/polyset3b[[i3]])),varset3]==True,del=del~Join~{{j3}}],{j3,i3+1,Length[polyset3b]}],{i3,1,Length[polyset3b]-1}];
polyset3b=Delete[polyset3b,del];Return[polyset3b];
)

(* If g is a multiple of f, then we delete g *)
DelDupPoly3[polyset3a_,varset3_]:=(polyset3b=polyset3a;
del={};
Do[
Do[tf3=0;If[PolynomialQ[Simplify@(Expand@(polyset3b[[i3]]/polyset3b[[j3]])),varset3]==True,
del=del~Join~{{i3}};tf3=1];
If[tf3==0&&PolynomialQ[Simplify@(Expand@(polyset3b[[j3]]/polyset3b[[i3]])),varset3]==True,del=del~Join~{{j3}}],{j3,i3+1,Length[polyset3b]}],{i3,1,Length[polyset3b]-1}];
polyset3b=Delete[polyset3b,del];Return[polyset3b];
)

(* Apply above “simplesubstitution”, change each member to a polynomial, and delete constant polynomials and duplicates.*)
SimpleSubstitution[S2_,varset2_]:=(s2=S2;tf2=1;While[tf2==1,s2=simpleSubstitution[s2,varset2];
s2=CancelDenominator[s2];
s2=CancelFactor[s2];
s2=Complement[s2,{1,0,-1}];
s2=DelDupPoly[s2,varset2];
];Return[s2])

(* If f=((Subscript[g, 1])^(Subscript[k, 1]))((Subscript[g, 2])^(Subscript[k, 2]))···((Subscript[g, n])^(Subscript[k, n])) , then replace f by (Subscript[g, 1])(Subscript[g, 2])···(Subscript[g, n])*)
MinimizeExpofFac[S4_]:=(s4=S4;L4=Length[s4];If[L4>0,Do[
list4=FactorList[s4[[k4]]];If[Length[list4]>1, s4[[k4]]=1;
Do[poly4=list4[[j4,1]];
If[Length@poly4>1,
s4[[k4]]=ExpandAll@(s4[[k4]]*poly4)],{j4,2,Length[list4]}]]
,{k4,1,L4}]];Return[s4];)

(* If f=((Subscript[g, 1])^(Subscript[k, 1]))((Subscript[g, 2])^(Subscript[k, 2]))···((Subscript[g, n])^(Subscript[k, n]))((Subscript[h, 1])^(Subscript[l, 1]))((Subscript[h, 2])^(Subscript[l, 2]))···((Subscript[h, m])^(Subscript[l, m])) and all terms in (Subscript[h, i])'s are of the same sign , then replace f by (Subscript[g, 1])(Subscript[g, 2])···(Subscript[g, n])*)
MinimizeExpofFacP[S3_]:=(s3=S3;L3=Length[s3];If[L3>0,Do[
list3=FactorList[s3[[k3]]];If[Length[list3]>1, s3[[k3]]=1;Do[
If[Length@StringPosition[ToString@(list3[[j3,1]]),"-"]>0&&Length@StringPosition[ToString@(-list3[[j3,1]]),"-"]>0,s3[[k3]]=s3[[k3]]*list3[[j3,1]]],
{j3,2,Length[list3]}]]
,{k3,1,L3}]];Return[s3];)

(* Elimination process: apply some algorithm on a collection of polynomial to elimi a given variable.*)
Elimi[S2_,var2_]:=(s2=Select[S2,Intersection[Variables[#],{var2}]≠{}& ];
s2b=Complement[S2,s2];
exp2c=2*Max[Table[Exponent[s2[[k2]],var2],{k2,1,Length[s2]}]];kk2=0;possfac2=1;
While[(Length@s2>1)&&(kk2<exp2c+2),
s2=SortBy[s2,Length[#]& ];
s2=SortBy[s2,Exponent[#,var2]& ];
exp2a=Exponent[s2[[1]],var2];
CoeList2=CoefficientList[s2[[1]],var2];
TF2=1;
CoeFactor2=Replace[CoeList2,List->PolynomialGCD,1,Heads->True];
Do[If[Exponent[CoeFactor2,var2b]>0,TF2=0];,{var2b,VarSet1b}];
If[TF2==1,
s2[[1]]=ExpandAll@Simplify[s2[[1]]/CoeFactor2];
possfac2=possfac2*CoeFactor2;
If[Total@Exponent[CoeFactor2,VarSet1a]>0,
Print[Style[" new common factor of coefficients : ",RGBColor[0,0.5,0]],Style[CoeFactor2,RGBColor[0,0.5,0]],Style["     ( ",RGBColor[0,0.5,0]],Style[Length@CoeFactor2,RGBColor[0,0.5,0]],Style[" terms )",RGBColor[0,0.5,0]]];
Print[];]];
Do[
exp2b=Exponent[s2[[j2]],var2];
If[exp2b>=exp2a,
poly2=Expand[s2[[1]]*var2^(exp2b-exp2a)];
Coe2a=Coefficient[poly2,var2^exp2b];
Coe2b=Coefficient[s2[[j2]],var2^exp2b];
s2[[j2]]=ExpandAll@(s2[[j2]]*Coe2a-poly2*Coe2b);
],{j2,2,Length[s2]}
];
s2=CancelFactor[s2];
s2=MinimizeExpofFac[s2];
s2b=s2b~Join~Select[s2,Intersection[Variables[#],{var2}]=={}&];
s2=Select[s2,Intersection[Variables[#],{var2}]!={}&];
kk2++;];
s2=SortBy[s2,Exponent[#,var2]& ];
Return[{s2b,possfac2}];
)

(* for each v in the varset,delete the polynomial (if it exists) in S5 that is the unique polynomial involving variable v.*)
DeleteVar[S2_,varset2_]:=(s2=S2;If[Length[s2]>1,
Do[If[Length[Select[s2,(Variables[#]⋂{varset2[[k2]]})≠{}& ]]==1,del={};
Do[If[(Variables[s2[[j2]]]⋂{varset2[[k2]]})≠{},
del=del~Join~{{j2}}],{j2,1,Length[s2]}];s2=Delete[s2,del]],{k2,1,Length[varset2]}]];Return[s2];
)


(* First algorithm for finding mass relations *)
FindMassRel1[]:=(
t1=AbsoluteTime[];
Str="finding the mass relation of FMX[["<>ToString[N1]<>"]]";Print[Str];
Print[];
TF1=0;
EqnZ=(ZX1*Z).M;
EqnZ=Relation1@Complement[EqnZ,{0}];
EqnW=(WX1*W).M;
EqnW=Relation1@Complement[EqnW,{0}];

Polyset1=Union[EqnW,EqnZ];

Print[Style["Collected equations: ",Bold,15,Blue]];
If[Length@EqnZ+Length@EqnW>0,Print[Style["Main equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",EqnZ[[k1]]],{k1,1,Length@EqnZ}];
Do[Print[0," = ",EqnW[[k1]]],{k1,1,Length@EqnW}];
];
If[Length@ClusterZ[{zo1,wo1}]+Length@ClusterW[{zo1,wo1}]>0,Print[Style["Cluster equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",ClusterZ[{zo1,wo1}][[k1]]],{k1,1,Length@ClusterZ[{zo1,wo1}]}];
Do[Print[0," = ",ClusterW[{zo1,wo1}][[k1]]],{k1,1,Length@ClusterW[{zo1,wo1}]}];
];Print[];

Polyset1=Relation1[Polyset1];
Polyset1=Relation2[Polyset1];
Polyset1=ClusterReplace[{zo1,wo1},Polyset1];
Polyset1=Complement[Polyset1,{1,0,-1}];
Polyset1b=GroebnerBasis[Polyset1,M~Join~mu];Polyset1=Polyset1b;
Polyset1=CancelFactor[Polyset1];
Polyset1=Complement[Polyset1,{1,0,-1}];

If[MassRelationQ[Polyset1]==True,TF1=1;Polyset1=TakeMassRela[Polyset1];Print["Finished, we have mass relation ",Style[Polyset1,Blue]," for this diagram."]];

If[TF1==0,Print["We don't find an relation between only masses. We have the relations:"];Print[Polyset1]];

Print[];Print["time spent (seconds): ",AbsoluteTime[]-t1];Print[];Print[];Print[];
);


(* Second algorithm for finding mass relations *)
FindMassRel2[]:=(TF1=0;
t1=AbsoluteTime[];
Str="finding the mass relation of FMX[["<>ToString[N1]<>"]]...";Print[Str];
Print[];
If[Length[zwotype3]>0,TF1=1;jj1=1;
While[jj1≤Length[zwotype3]&&TF1==1,
{S1a,S1b}=FindWedgeV[zwotype3[[jj1]]];If[S1a=={},TF1=0;Print["some zw-order matrices don't satisfy the hypothesis that the system contains a wedge equation."];];jj1++];

If[TF1==1,Print["all zw-order matrices satisfy the hypothesis that the system contains a wedge equation."];Print[];massrellcm=1;PossFaclcm=1;jj1=1;
While[jj1≤Length[zwotype3],
{S1a,S1b}=FindWedgeV[zwotype3[[jj1]]];

WedgeV1a=Flatten@Table[M[[j1]]M[[i1]](zeta[[i1]]omega[[j1]])^(-1),{i1,S1a},{j1,S1b}];sign1=Tuples[{1,-1},{Length[S1a],Length[S1b]}];
WedgeV1b=Table[Dot[WedgeV1a,Flatten[sign1[[j1]],1]],{j1,1,Length[sign1]}];
VarSet1a={};VarSet1a=Variables[WedgeV1a];
WedgeV1b=WedgeV1b[[1;;2^(Length[S1b]*Length[S1a]-1)]];
Do[WedgeV1b[[j1]]=ExpandAll@WedgeV1b[[j1]],{j1,1,Length@WedgeV1b}];
index1=1;

Circledzc=Table[Intersection[zc1[[i1]],Flatten[Position[Diagonal[ZX1],1]]],{i1,1,Length[zc1]}];
Circledwc=Table[Intersection[wc1[[i1]],Flatten[Position[Diagonal[WX1],1]]],{i1,1,Length[wc1]}];
MassCenRela=Flatten[Table[ Sum[M[[ Circledzc[[i1,j1]]      ]]*z[[ Circledzc[[i1,j1]]      ]],  {j1,1,Length[Circledzc[[i1]]]}    ] ,{i1,1,Length[Circledzc]}]]~Join~Flatten[Table[ Sum[M[[Circledwc[[i1,j1]]      ]]*w[[ Circledwc[[i1,j1]]      ]],  {j1,1,Length[Circledwc[[i1]]]}    ] ,{i1,1,Length[Circledwc]}]];

Str="Finding mass relation with zwotype3[["<>ToString[jj1]<>"]] ";Print[Str];
Print[];
Print["The order matrices are "];Print[{MatrixForm[zwotype3[[jj1,1]]],MatrixForm[zwotype3[[jj1,2]]]}];
Print[];
Print[Style["Collected equations: ",Bold,15,Blue]];
If[Length@MassCenRela>0,Print[Style["Conditions of MassCenter: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",MassCenRela[[k1]]],{k1,1,Length@MassCenRela}];
];
If[Length@ClusterZ[{zo1,wo1}]+Length@ClusterW[{zo1,wo1}]>0,Print[Style["Cluster equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",ClusterZ[{zo1,wo1}][[k1]]],{k1,1,Length@ClusterZ[{zo1,wo1}]}];
Do[Print[0," = ",ClusterW[{zo1,wo1}][[k1]]],{k1,1,Length@ClusterW[{zo1,wo1}]}];
];
Print[Style["Wedge equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Replace[Prepend[Insert[WedgeV1a," ± ",Table[{k1},{k1,2,Length@WedgeV1a}]],"one of     0 = "],List->Print,1,Heads->True];
Print[];

While[index1≤Length[WedgeV1b],

Print[];
Polyset1=MassCenRela;
Polyset1=Complement[Polyset1,{1,0,-1}];
Polyset1=Polyset1~Join~{WedgeV1b[[index1]]};
Str="running wedge equation "<>ToString[index1]<>"...";Print[Str];
Print[];
Print["polynomials: "];
Print[Polyset1];Print[];
PossFac=1;
Polyset1=Relation1[Polyset1];
Polyset1=Relation2[Polyset1];
Polyset1=ChangeVar1[Polyset1];
Polyset1=ChangeVar2[Polyset1];
Polyset1=CancelDenominator[Polyset1];
Polyset1=ClusterReplace[zwotype3[[jj1]],Polyset1];

Polyset1=Table[Simplify[Polyset1[[j1]]],{j1,1,Length[Polyset1]}];
Polyset1=Complement[Polyset1,{1,0,-1}];

VarSet1a={};VarSet1a=Variables[Polyset1];
VarSet1b=Complement[VarSet1a,Union[M,mu]];
VarSet1b=SortBy[VarSet1b,Max[Table[Exponent[Polyset1[[k1]],#],{k1,1,Length[Polyset1]}]]&];
Polyset1=MinimizeExpofFac[Polyset1];

While[Length[VarSet1b]>0&&Length[Polyset1]>1&&MassRelationQ[Polyset1]==False,

VA1=Elimi[Polyset1,VarSet1b[[1]]];
Polyset1=VA1[[1]];
PossFac=PossFac*VA1[[2]];
Polyset1=Complement[Polyset1,{0}];
Polyset1=CancelFactor[Polyset1];
If[Length[Polyset1]>1,
Polyset1=DeleteVar[Polyset1,VarSet1b];
Polyset1=DelDupPoly2[Polyset1,VarSet1b];
];
VarSet1a={};Do[VarSet1a=VarSet1a⋃Variables[Polyset1[[k1]]],{k1,1,Length[Polyset1]}];
VarSet1b=Complement[VarSet1a,Union[M,mu]];
VarSet1b=SortBy[VarSet1b,Max[Table[Exponent[Polyset1[[k1]],#],{k1,1,Length[Polyset1]}]]&];
];

Polyset1=CancelFactor[Polyset1];
Polyset1=Complement[Polyset1,{1,0,-1}];

If[MassRelationQ[Polyset1]==True,
Polyset1=Simplify[TakeMassRela[Polyset1]];
Polyset1=MinimizeExpofFacP[Polyset1];
Polyset1=ChangeVar2Inv[Polyset1];
PossFac=ChangeVar2Inv[{PossFac}][[1]];
PossFac=MinimizeExpofFacP[{PossFac}][[1]];

Str="for zwotype3[["<>ToString[jj1]<>"]] and WedgeV[["<>ToString[index1]<>"]], we have mass relations";Print[Str];Print[Style[Polyset1,Blue]];Do[massrellcm=PolynomialLCM[massrellcm,poly1];,{poly1,Polyset1}];
PossFaclcm=PolynomialLCM[PossFaclcm,PossFac];
Print[" or  ",Style[PossFac,RGBColor[0,0.5,0]]," (the part obtained by factors of coefficients )"];Print[];Print[];];

If[MassRelationQ[Polyset1]==False,TF1=0;Str="We don't obtain any mass relation for zwotype3[["<>ToString[jj1]<>"]] and WedgeV[["<>ToString[index1]<>"]], we have relations";Print[Str];Print[Polyset1];];
Print[];Print[];
index1++;];jj1++

];];
If[TF1==1,Print["Finished, for this diagram , we have mass relations ",Style[massrellcm,Blue], " or ",Style[PossFaclcm,RGBColor[0,0.5,0]]," (the part obtained by factors of coefficients )" ],Print["Finished, we don't obtain any mass relations for this diagram."]];

];
Print[];Print["time spent (seconds): ",AbsoluteTime[]-t1];Print[];Print[];Print[];
)


(* Third algorithm for finding mass relations *)
FindMassRel3[]:=FindMassRel3[{zo1,wo1}];
FindMassRel3[{zomatrix1_,womatrix1_}]:=FindMassRel3[{zomatrix1,womatrix1},allind,allind];
FindMassRel3[{zomatrix1_,womatrix1_},vertexSet1a_,vertexSet1b_]:=FindMassRel3[{zomatrix1,womatrix1},vertexSet1a,vertexSet1b,False,{}];
FindMassRel3[{zomatrix1_,womatrix1_},vertexSet1a_,vertexSet1b_,RandomTF_,AddEqnSet_]:=(
t1=AbsoluteTime[];

ODind1=Flatten[Table[Table[{j1,i1},{i1,j1+1,n}],{j1,1,n-1}],1];
EdgePos1=Select[ODind1,XX1[[#[[1]],#[[2]]]]>0&];
TriangleSet={};
Do[Do[Do[If[Complement[{ {i1,j1},{j1,k1},{i1,k1}},EdgePos1]=={},
AppendTo[TriangleSet,{i1,j1,k1}];
];
,{i1,1,j1-1}];
,{j1,2,k1-1}];
,{k1,3,n}];

(*  (Subscript[z, ij])+(Subscript[z, jk])=(Subscript[z, ik]) ,  (Subscript[w, ij])+(Subscript[w, jk])=(Subscript[w, ik]) *)
Idenzz=Table[zz[[s1[[1]],s1[[2]]]]+zz[[s1[[2]],s1[[3]]]]-zz[[s1[[1]],s1[[3]]]],{s1,TriangleSet}];
Idenww=Table[ww[[s1[[1]],s1[[2]]]]+ww[[s1[[2]],s1[[3]]]]-ww[[s1[[1]],s1[[3]]]],{s1,TriangleSet}];
(*  Substitute (Subscript[u, ij])'s , (Subscript[v, ij])'s for ((Subscript[z, ij])^(-1/2)) , ((Subscript[w, ij])^(-1/2)) *)
Idenu=Table[u[[s1[[1]],s1[[2]]]]^(-2)+u[[s1[[2]],s1[[3]]]]^(-2)-u[[s1[[1]],s1[[3]]]]^(-2),{s1,TriangleSet}];
Idenv=Table[v[[s1[[1]],s1[[2]]]]^(-2)+v[[s1[[2]],s1[[3]]]]^(-2)-v[[s1[[1]],s1[[3]]]]^(-2),{s1,TriangleSet}];

Str="finding the mass relation of FMX[["<>ToString[N1]<>"]] with zw-order matrix ";Print[Str];
Print[{MatrixForm[outputform1@zomatrix1],MatrixForm[outputform1@womatrix1]}];
Print[];
Print[];
TF1=0;
PossFac=1;
EqnZ=Relation1@EqnMZO[{zomatrix1,womatrix1},vertexSet1a];
EqnW=Relation1@EqnMWO[{zomatrix1,womatrix1},vertexSet1b];
Polyset1=EqnZ⋃EqnW⋃AddEqnSet;
Polyset1b={};

Print[Style["Collected equations: ",Bold,15,Blue]];
If[Length@EqnZ+Length@EqnW>0,Print[Style["Main equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",EqnZ[[k1]]],{k1,1,Length@EqnZ}];
Do[Print[0," = ",EqnW[[k1]]],{k1,1,Length@EqnW}];
];
If[Length@ClusterZ[{zomatrix1,womatrix1}]+Length@ClusterW[{zomatrix1,womatrix1}]>0,Print[Style["Cluster equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",ClusterZ[{zomatrix1,womatrix1}][[k1]]],{k1,1,Length@ClusterZ[{zomatrix1,womatrix1}]}];
Do[Print[0," = ",ClusterW[{zomatrix1,womatrix1}][[k1]]],{k1,1,Length@ClusterW[{zomatrix1,womatrix1}]}];
];
VarSet1c=Variables[Relation3ver1[Relation2@Polyset1]];
Polyset1c={};
Do[If[Complement[Variables[Idenzz[[k1]]],VarSet1c]=={},Polyset1c=Polyset1c~Join~{Idenzz[[k1]]}],{k1,1,Length[Idenu]}];
Do[If[Complement[Variables[Idenww[[k1]]],VarSet1c]=={},Polyset1c=Polyset1c~Join~{Idenww[[k1]]}],{k1,1,Length[Idenv]}];
If[Length@Polyset1c>0,
Print[Style["Relations between separations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",Polyset1c[[k1]]],{k1,1,Length@Polyset1c}];
];
If[Length@AddEqnSet>0,
Print[Style["additional equations: ",Bold,13,RGBColor[0.8,0.4,0]]];
Do[Print[0," = ",AddEqnSet[[k1]]],{k1,1,Length@AddEqnSet}];
];Print[];

Polyset1=Relation1[Polyset1];
Polyset1=Relation2[Polyset1];
Polyset1=Relation3ver2[Polyset1];
VarSet1a=Variables[Polyset1];
VarSet1b=Complement[VarSet1a,Union[M,mu]];
Do[If[Complement[Variables[Idenu[[k1]]],VarSet1a]=={},Polyset1=Polyset1~Join~{Idenu[[k1]]}],{k1,1,Length[Idenu]}];
Do[If[Complement[Variables[Idenv[[k1]]],VarSet1a]=={},Polyset1=Polyset1~Join~{Idenv[[k1]]}],{k1,1,Length[Idenv]}];
Polyset1=ClusterReplace[{zomatrix1,womatrix1},Polyset1];
Polyset1=CancelDenominator[Polyset1];
Polyset1=CancelFactor[Polyset1];
Polyset1=Complement[Polyset1,{1,0,-1}];
Polyset1=ChangeVar3[Polyset1];
Polyset1=CancelDenominator[Polyset1];
Polyset1=CancelFactor[Polyset1];
Polyset1=Complement[Polyset1,{1,0,-1}];
VarSet1a={};Do[VarSet1a=VarSet1a⋃Variables[Polyset1[[k1]]],{k1,1,Length[Polyset1]}];
VarSet1b=Complement[VarSet1a,Union[M,mu]];
VarSet1b=SortBy[VarSet1b,Max[Table[Exponent[Polyset1[[k1]],#],{k1,1,Length[Polyset1]}]]&];
str="Change of variables ";
Print[str];
Print[Polyset1];Print[ToString[Length[Polyset1]]<>" polynomials, "<>ToString[Length[VarSet1b]]<>" variables "<>ToString[VarSet1b]];Print[];

str="Simple Substitutions ";
Print[str];
Polyset1=SimpleSubstitution[Polyset1,VarSet1b];
Polyset1=CancelFactor[Polyset1];
Polyset1=Complement[Polyset1,{1,0,-1}];
Polyset1=DelDupPoly2[Polyset1,VarSet1b];
Polyset1b=Polyset1b⋃TakeMassRela[Polyset1];
Polyset1=Complement[Polyset1,Polyset1b];
VarSet1a={};Do[VarSet1a=VarSet1a⋃Variables[Polyset1[[k1]]],{k1,1,Length[Polyset1]}];
VarSet1b=Complement[VarSet1a,Union[M,mu]];
VarSet1b=SortBy[VarSet1b,Max[Table[Exponent[Polyset1[[k1]],#],{k1,1,Length[Polyset1]}]]&];
Print[Polyset1b~Join~Polyset1];Print[ToString[Length[Polyset1b~Join~Polyset1]]<>" polynomials, "<>ToString[Length[VarSet1b]]<>" variables "<>ToString[VarSet1b]];Print[];
Polyset1=MinimizeExpofFac[Polyset1];

While[Length[VarSet1b]>0&&Length[Polyset1]>1&&MassRelationQ[Polyset1]==False,
If[RandomTF==True,Var1=RandomChoice[VarSet1b],Var1=VarSet1b[[1]]];
str="Eliminating variable "<>ToString[Var1]<>"...";
Print[str];Print[];
VA1=Elimi[Polyset1,Var1];
Polyset1=VA1[[1]];
PossFac=PossFac*VA1[[2]];
Polyset1=Complement[Polyset1,{0}];
Polyset1=CancelFactor[Polyset1];
Polyset1b=Polyset1b⋃TakeMassRela[Polyset1];
Polyset1=Complement[Polyset1,Polyset1b];

If[Length[Polyset1]>1,
Polyset1=DeleteVar[Polyset1,VarSet1b];
Polyset1=DelDupPoly2[Polyset1,VarSet1b];
];
VarSet1a={};Do[VarSet1a=VarSet1a⋃Variables[Polyset1[[k1]]],{k1,1,Length[Polyset1]}];
VarSet1b=Complement[VarSet1a,Union[M,mu]];
VarSet1b=SortBy[VarSet1b,Max[Table[Exponent[Polyset1[[k1]],#],{k1,1,Length[Polyset1]}]]&];
str="after eliminating "<>ToString[Var1];
Print[str];
NV1=Table[Length[f1],{f1,Polyset1b}];
NV1=NV1~Join~Table[Length[f1],{f1,Polyset1}];
If[Max[NV1]≤100,Print[Polyset1b~Join~Polyset1]];
Print[ToString[Length[Polyset1b~Join~Polyset1]]<>" polynomials, "<>ToString[Length[VarSet1b]]<>" variables "<>ToString[VarSet1b]];
Print["numbers of terms : ",NV1];Print[];Print[];
];

If[Polyset1b≠{},
Polyset1b=DelDupPoly3[Polyset1b,VarSet1b];
Polyset1b=Complement[Polyset1b,{1,-1}];
Polyset1b=Simplify[Polyset1b];
Polyset1b=MinimizeExpofFacP[Polyset1b];
Str="we have mass relations";Print[Str];Print[Style[Polyset1b,Blue]];
If[PossFac=!=1,Print[" or  ",Style[PossFac,RGBColor[0,0.5,0]]," (the part obtained by factors of coefficients )"]];
Print[];Print[];

Polyset1b=ChangeVar2[Polyset1b];
Polyset1b=MinimizeExpofFacP[Polyset1b];
Polyset1b=ChangeVar2Inv[Polyset1b];
PossFac=ChangeVar2[{PossFac}][[1]];
PossFac=MinimizeExpofFacP[{PossFac}][[1]];
PossFac=ChangeVar2Inv[{PossFac}][[1]];
Print["another expression:"];
Print[Style[Polyset1b,Blue]];
If[PossFac=!=1,Print[" or  ",Style[PossFac,RGBColor[0,0.5,0]]," (the part obtained by factors of coefficients )"]];
Print[];,

If[Length[Polyset1]>1,
FacFamilyV1=Table[{},{k1,1,Length[Polyset1]}];
Print[FacFamilyV1];Print[];
Do[FacFamilyV1[[k1]]=Table[FactorList[Polyset1[[k1]]][[j1,1]],{j1,1,Length@FactorList[Polyset1[[k1]]]}],{k1,1,Length@FacFamilyV1}];
Print[FacFamilyV1];Print[];
Do[Polyset1[[k1]]=1;
Do[Do[TF1b=1;If[PolynomialQ[Simplify@(Expand[Polyset1[[j1]]/FacFamilyV1[[k1]][[i1]]]),Variables[Polyset1[[j1]]]]==True,TF1b=TF1b,TF1b=TF1b*0];
,{j1,1,Length[Polyset1]}];If[TF1b==1,Polyset1[[k1]]=Polyset1[[k1]]*FacFamilyV1[[k1,i1]]];
,{i1,1,Length[FacFamilyV1[[k1]]]}];
,{k1,1,Length[Polyset1]}];

Polyset1=CancelDenominator[Polyset1];
Polyset1=CancelFactor[Polyset1];
Polyset1=DeleteVar[Polyset1,VarSet1b];
Polyset1=DelDupPoly2[Polyset1,VarSet1b];
];

If[Length[Polyset1]==1,Polyset1b=FactorList[Polyset1[[1]]];Polyset1[[1]]=1;
Do[Polyset1[[1]]=Polyset1[[1]]*Polyset1b[[i1]][[1]],{i1,1,Length[Polyset1b]}]];

Print[Polyset1];Print[ToString[Length[Polyset1]]<>" polynomials"];Print[];];

Print[];Print["time spent (seconds): ",AbsoluteTime[]-t1];Print[];Print[];Print[];
)

FindOrder[107];FindMassRel1[]

FindOrder[77];Do[FindMassRel3[zwotype2[[k]],{1,2},{1,2,3,4,5}],{k,1,Length@zwotype2}]

FindOrder[7,2];FindMassRel3[{zo1,wo1},{5,6},{5,6}]

FindOrder[31,2];FindMassRel3[{zo1,wo1},{4,5},{4,5}]

FindOrder[20];Do[FindMassRel3[zwotype2[[k]],{1,2,3,4,5,6},{1,2,3,4,5,6}],{k,1,Length@zwotype2}]