(******************* EVALUATE THE FOLLOWING CODE FIRST FOR BASIC USAGE **************)

(*Function to compute Bruhat Decomposition*)BruhatDecomposition[M_?MatrixQ]:=Module[{n,mat,U1,U2,P,T,availableRows,pivotRow,pivotVal,alpha,beta,pivots,j,i,c},n=Length[M];
mat=M;(*Working copy of the matrix*)(*Initialize matrices*)U1=IdentityMatrix[n];
U2=IdentityMatrix[n];
P=ConstantArray[0,{n,n}];
T=ConstantArray[0,{n,n}];
availableRows=Range[n];(*Track which rows are not yet used as pivots*)pivots=ConstantArray[0,n];(*Store pivot row index for each column*)(*Iterate through each column j from 1 to n*)Do[(*1. Find the PIVOT:Bottom-most non-zero entry in column j among available rows*)(*We use Simplify to handle symbolic expressions correctly*)pivotRow=SelectFirst[Reverse[availableRows],!TrueQ[Simplify[mat[[#,j]]==0]]&];
If[MissingQ[pivotRow],Print["Error: Matrix is singular or symbolic cancellation failed."];
Abort[]];
pivots[[j]]=pivotRow;
pivotVal=mat[[pivotRow,j]];
(*Mark this row as used*)availableRows=DeleteCases[availableRows,pivotRow];
(*2. ROW OPERATIONS:Eliminate entries ABOVE the pivot (i<pivotRow)*)(*Operation:Row_i=Row_i-alpha*Row_pivot*)(*Update U1:Col_pivot=Col_pivot+alpha*Col_i*)Do[If[i<pivotRow,alpha=Simplify[mat[[i,j]]/pivotVal];
If[alpha=!=0,mat[[i]]=Simplify[mat[[i]]-alpha*mat[[pivotRow]]];
U1[[All,pivotRow]]=Simplify[U1[[All,pivotRow]]+alpha*U1[[All,i]]];];],{i,1,n}];
(*3. COLUMN OPERATIONS:Eliminate entries to the RIGHT of the pivot (c>j)*)(*Operation:Col_c=Col_c-beta*Col_j*)(*Update U2:Row_j=Row_j+beta*Row_c*)Do[beta=Simplify[mat[[pivotRow,c]]/pivotVal];
If[beta=!=0,mat[[All,c]]=Simplify[mat[[All,c]]-beta*mat[[All,j]]];
U2[[j]]=Simplify[U2[[j]]+beta*U2[[c]]];],{c,j+1,n}];
(*4. Build P and T components*)P[[pivotRow,j]]=1;
T[[j,j]]=pivotVal;,{j,1,n}];
(*Return the list of matrices*){U1,P,T,U2}];

(*************************************** Path matrices ******************************)

(*Define the function P(z,k,n)*)P[z_,k_,n_]:=Module[{mat},(*Start with an n x n Identity Matrix*)mat=IdentityMatrix[n];
(*Modify the 4 specific entries*)mat[[k,k]]=0;
mat[[k,k+1]]=1;
mat[[k+1,k]]=-1;
mat[[k+1,k+1]]=z;
(*Return the resulting matrix*)mat];
Q[z_,k_,n_]:=Module[{mat},(*Start with an n x n Identity Matrix*)mat=IdentityMatrix[n];
(*Modify the 4 specific entries*)mat[[k,k]]=z;
mat[[k,k+1]]=-1;
mat[[k+1,k]]=1;
mat[[k+1,k+1]]=0;
(*Return the resulting matrix*)mat];
tildeW0[n_]:=SparseArray[{i_,j_}/;j==n-i+1:>(-1)^(i-1),{n,n}]//Normal
(*Example Input Variables*)
zVal=z;(*You can replace'z' with a number if desired*)kVal=2;(*Ensure k<n*)nVal=5;
(*Generate the matrix*)
resultMatrix=P[zVal,kVal,nVal];

tildeW0[n_]:=SparseArray[{i_,j_}/;j==n-i+1:>(-1)^(i-1),{n,n}]//Normal

leadingPrincipalMinors[m_?MatrixQ]:=Module[{n},(*Check if the matrix is square*)If[Dimensions[m][[1]]≠Dimensions[m][[2]],Return["Error: Matrix must be square"]];
n=Length[m];
(*Print["Leading Principal Minors:"];*)
(*Generate the determinant for each k x k submatrix*)Simplify[Table[Det[m[[1;;k,1;;k]]],{k,1,n}]]
]

(**********)
antiDiagonalPrefixProducts[m_?MatrixQ]:=Module[{n},n=Length[m];
(*Generate the list of anti-diagonal entries:(n,1),(n-1,2),...,(1,n)*)Table[Product[m[[n-j+1,j]],{j,1,k}],{k,1,n}]]
(*************)

Bruhats[M_]:=Module[{Msym2=M},
(*Print["--- Input Matrix M ---"];*)
Print["Morse Differential at this step:"];
Print[MatrixForm[Msym2//.Table[z[i]->Subscript[z,i],{i,0,20}]]];
{U1res,Pres,Tres,U2res}=BruhatDecomposition[Msym2];
Print["Its Bruhat decomposition components U1.Rook.U2 now follow."];
Print["U1 (First Upper Unipotent):"];
Print[MatrixForm[U1res//.Table[z[i]->Subscript[z,i],{i,0,20}]]];
(*Print["P (Permutation):"];
Print[MatrixForm[Pres//.Table[z[i]->Subscript[z,i],{i,0,20}]]];*)
Print["Rook matrix:"];
Print[MatrixForm[Expand[Pres.Tres]//.Table[z[i]->Subscript[z,i],{i,0,20}]]];Print["U2 (Second Upper Unipotent):"];
Print[MatrixForm[U2res//.Table[z[i]->Subscript[z,i],{i,0,20}]]];
antiDiagonalPrefixProducts[Expand[Pres.Tres]]//.Table[z[i]->Subscript[z,i],{i,0,20}]

(*Print["--- Verification ---"];*)
(*Print["M == U1 . P . T . U2 ?"];*)
(*recon=Simplify[U1res.Pres.Tres.U2res];*)
(*Print[recon==Msym2];*)
(*If[recon===Simplify[Msym2],Print["SUCCESS"]];*)]

BruhatAndClusterVariables[numList_List,n_]:=Module[{variables,Bmatrices,Pmatrices,Pproducts,Bproducts},

(*1. Create a list of indexed variables z_1,z_2,...*)
variables=Table[Subscript[z,k],{k,Length[numList]}];

(*2. Create the individual matrices P[n,i_k,z_k]*)
Pmatrices=MapThread[P[#1,#2,n]&,{variables,numList}];
Bmatrices=MapThread[Q[#1,#2,n]&,{variables,numList}];

(*3. Compute cumulative products in reverse order:P_k*...*P_1*)
(*FoldList applies the function to the previous result and the next element*)
Pproducts=FoldList[Dot[#2,#1]&,First[Pmatrices],Rest[Pmatrices]];
Bproducts=FoldList[Dot[#1,#2]&,First[Bmatrices],Rest[Bmatrices]];

(*4. Simplify and Print each product in MatrixForm*)
Do[Print["---------------- Data right after crossing number ",k,"----------------"];
(*Print[Simplify[Bruhats[Pproducts[[k]].tildeW0[n]]]];*)
Print["        ## Cerf theory data is displayed first, as follows ##"];
indk=n-numList[[k]];
(*Print[(-1)^(indk*(indk-1)/2+indk*(n-1))*Expand[Simplify[Bruhats[Pproducts[[k]].tildeW0[n]]][[n-numList[[k]]]]]];*)
BruhatProd=(-1)^(indk*(indk-1)/2+indk*(n-1))*Expand[Simplify[Bruhats[Pproducts[[k]].tildeW0[n]]][[n-numList[[k]]]]];
Print[Style["Product of corresponding Bruhat numbers for the Morse differential, extracted from rook matrix, with corrected sign (RHS of Comparison Lemma):",Blue]];
Print[Framed[BruhatProd,FrameStyle->Blue,Background->LightYellow]];
Print["         ## Cluster theory data for braid variety (equiv. augmentation variety now follows ##"];
Print["Braid matrix at this step:"];
Print[MatrixForm[Simplify[Bproducts[[k]]]]];
(*******)
Print[Style["Cluster variable associated to this crossing, extracted from braid matrix (LHS of Comparison Lemma):",Red]];
Print[Framed[Expand[Simplify[leadingPrincipalMinors[Bproducts[[k]]][[numList[[k]]]]]],FrameStyle->Red,Background->LightYellow]];
(*If[(-1)^(indk*(indk-1)/2+indk*(n-1))*Expand[Simplify[leadingPrincipalMinors[Bproducts[[k]]][[numList[[k]]]]]]==Expand[Simplify[Bruhats[Pproducts[[k]].tildeW0[n]]][[n-numList[[k]]]]],Print["SUCCESS"],Print["FAIL"]];*)
(*Print["test:"];*)
(*If[Expand[Simplify[leadingPrincipalMinors[Bproducts[[k]]][[numList[[k]]]]]]!=Expand[Simplify[Bruhats[Pproducts[[k]].tildeW0[n]]][[n-numList[[k]]]]],Print["FAIL"]]*)
Print["----------------------------------------------------------------"],{k,Length[Pproducts]}];
]

(******************* INSTRUCTIONS FOR BASIC USAGE **************)

(* Given an n-stranded positive braid word, expressed as an ordered list of crossings L={i1,i2,...,il}, with indices i1,i2,...,il between 1 and n-1, call function:
			
			BruhatAndClusterVariables[{i1,i2,...,il},n]

The output contains the information of the Morse differentials, rook matrices and associated Bruhat products, and then the braid matrices and associated cluster variables.*)

(******************* EXAMPLES OF USAGE **************)

(********* Example in 2-strands for braid word 1,1,1,1,1 ************)
BruhatAndClusterVariables[{1,1,1,1,1},2]

(********* Example in 3-strands for braid word 1,2,1,2,1,1,2,2,1 ************)
BruhatAndClusterVariables[{1,2,1,2,1,1,2,2,1},3]

(********* Example in 4-strands for braid word 1,2,3,1,3,2,1,3,2,2,1 ************)
BruhatAndClusterVariables[{1,2,3,1,3,2,1,3,2,2,1},4]