(* This code verifies that Rudin's triangulation is 2-regular as stated
 in the article "A result on flip-graph connectivity" by L. Pournin. 
The computation done with mathematica is formal, showing that this
  result actually holds for the given set of height functions w1, w2 *)

(* Clears all variables *)
Clear["Global`*"];

(* Here 'epsilon' denotes the length of the shortest side of a triangle
  whose longest side is of length 1 and two of whose angles are 1° and 
  61°, as defined in the original article by Rudin (Rudin, Bulletin of
the AMS 64, 1958) *)
epsilon=Sin[Pi/180]/Sin[61*Pi/180];

(* Definition of the coordinates of the 14 vertices with same notations
  as in the original article. Vertices Xi are defined explicitely, and 
  all the other vertices are defined as functions of the Xi *)
X1={0,0,0};
X2={1,0,0};
X3={1/2,Sqrt[3]/2,0};
X4={1/2,1/(2*Sqrt[3]),Sqrt[2/3]};
U1=(X1+X3)/2;
U2=(X2+X4)/2;
Y1=(Cos[Pi/180]*epsilon*Sqrt[3]/2)*(X3-X1)+(Sin[Pi/180]*epsilon)*(X4-U1)+X1;
Y2=(Cos[Pi/180]*epsilon*Sqrt[3]/2)*(X4-X2)+(Sin[Pi/180]*epsilon)*(X1-U2)+X2;
Y3=(Cos[Pi/180]*epsilon*Sqrt[3]/2)*(X1-X3)+(Sin[Pi/180]*epsilon)*(X2-U1)+X3;
Y4=(Cos[Pi/180]*epsilon*Sqrt[3]/2)*(X2-X4)+(Sin[Pi/180]*epsilon)*(X3-U2)+X4;
Z1=(2*Sin[Pi/180]*Cos[59*Pi/180]/Sqrt[3])*(X2-X1)+(4*Sin[Pi/180]*Sin[59*Pi/180]/3)*(X4-(X1+X2)/2)+X1;
Z2=(2*Sin[Pi/180]*Cos[59*Pi/180]/Sqrt[3])*(X3-X2)+(4*Sin[Pi/180]*Sin[59*Pi/180]/3)*(X1-(X2+X3)/2)+X2;
Z3=(2*Sin[Pi/180]*Cos[59*Pi/180]/Sqrt[3])*(X4-X3)+(4*Sin[Pi/180]*Sin[59*Pi/180]/3)*(X2-(X3+X4)/2)+X3;
Z4=(2*Sin[Pi/180]*Cos[59*Pi/180]/Sqrt[3])*(X1-X4)+(4*Sin[Pi/180]*Sin[59*Pi/180]/3)*(X3-(X4+X1)/2)+X4;

(* The rows of the matrix 'vertices' correspond to the vertices of the
  triangulation with the following convention:
 Row 1: X1;
   Row 2: X2;
   Row 3: X3;
   Row 4: X4;
   Row 5: U1;
   Row 6: U2;
   Row 7: Y1;
   Row 8: Y2;
   Row 9: Y3;
   Row 10: Y4;
   Row 11: Z1;
   Row 12: Z2;
   Row 13: Z3;
   Row 14: Z4; *)

Vertices=FullSimplify[{X1,X2,X3,X4,U1,U2,Y1,Y2,Y3,Y4,Z1,Z2,Z3,Z4}];

(* With the notations of the original paper, the categories (1) to (11) 
  are now defined. Category (i) is represented by a matrix 'Categoryi'
  whose rows correspond to each tetrahedron and whose columns
  correspond to the vertices of those tetrahedron. Each of those
  vertices is represented as the index of its row (between 1 and 14)
  in the matrix 'Vertices' *)
Category1=(({{1, 11, 2, 7}, {2, 12, 3, 8}, {3, 13, 4, 9}, {4, 14, 1, 10}}));
Category2=(({{1, 12, 2, 7}, {2, 13, 3, 8}, {3, 14, 4, 9}, {4, 11, 1, 10}}));
Category3=(({{11, 12, 2, 7}, {12, 13, 3, 8}, {13, 14, 4, 9}, {14, 11, 1, 10}}));
Category4=(({{11, 12, 2, 8}, {12, 13, 3, 9}, {13, 14, 4, 10}, {14, 11, 1, 7}}));
Category5=(({{11, 12, 13, 14}}));
Category6=(({{11, 12, 7, 13}, {12, 13, 8, 14}, {13, 14, 9, 11}, {14, 11, 10, 12}}));
Category7=(({{1, 12, 7, 13}, {2, 13, 8, 14}, {3, 14, 9, 11}, {4, 11, 10, 12}}));
Category8=(({{1, 12, 9, 13}, {2, 13, 10, 14}, {3, 14, 7, 11}, {4, 11, 8, 12}}));
Category9=(({{1, 5, 9, 13}, {2, 6, 10, 14}, {3, 5, 7, 11}, {4, 6, 8, 12}}));
Category10=(({{1, 5, 7, 13}, {2, 6, 8, 14}, {3, 5, 9, 11}, {4, 6, 10, 12}}));
Category11=(({{11, 5, 7, 13}, {12, 6, 8, 14}, {13, 5, 9, 11}, {14, 6, 10, 12}}));

(* The matrix 'T1' contains the subset of the tetrahedra of Rudin's
  triangulation which will be shown to belong to the subdivision
  realized by height function w1. In this matrix, each row corresponds
  to a tetrahedron and each column to a vertex of those tetrahedra,
   with the same convention as above: each vertex is represented as the
index of its row (between 1 and 14) in the matrix 'Vertices' *)
T1={Category1[[2]],Category1[[4]],Category2[[2]],Category2[[4]],Category3[[2]],Category3[[4]],Category4[[2]],Category4[[4]],Category11[[2]],Category11[[4]],Category6[[1]],Category6[[3]],Category7[[1]],Category7[[3]],Category8[[1]],Category8[[3]],Category9[[1]],Category9[[3]],Category10[[1]],Category10[[3]],Category5[[1]]};

(* The matrix 'T2' contains the subset of the tetrahedra of Rudin's
  triangulation which will be shown to belong to the subdivision
  realized by height function w2. In this matrix, each row corresponds
  to a tetrahedron and each column to a vertex of those tetrahedra,
   with the same convention as above: each vertex is represented as the
index of its row (between 1 and 14) in the matrix 'Vertices' *)
T2={Category1[[1]],Category1[[3]],Category2[[1]],Category2[[3]],Category3[[1]],Category3[[3]],Category4[[1]],Category4[[3]],Category11[[1]],Category11[[3]],Category6[[2]],Category6[[4]],Category7[[2]],Category7[[4]],Category8[[2]],Category8[[4]],Category9[[2]],Category9[[4]],Category10[[2]],Category10[[4]],Category5[[1]]};

(* Matrix A1 of the inequation system A1.w1>0 is now built. Each 
 contraint of this system corresponds to a pair constituted of a
 tetrahedron (in matrix 'T1') and a vertex of the triangulation
 (in matrix 'Vertices'). Each constraint checks that the vertex,
   lifted by height function w1, is above the affine hull of the
  tetrahedron, lifted by the same height function in dimension 4.
 The number of constraints therefore the number of tetrahedra in
 T1 times the number of vertices of the triangulation that do not
 belong to the tetrahedron being checked (i.e. 10 vertices). As
 T1 contains 21 tetrahedra The system A1.w1>0 has 210 rows. Vector
 w1 has 14 rows (one per vertex). Hence A1 has 14 columns. The
 matrix 'A1' is first initialized with zeros and then filled with
 its actual coefficients. Those coefficients are subdeterminant of
 a determinant described by Takeushi (LCNS 2098, 330-338, 2001).
In the paper, one also finds a comprehensive description of the
 system A.w>0. *)
A1=ConstantArray[0,{210,14}];

(* k is the index on the row of matrix A that is being filled.
It is initialized as 1. *)  
k=1;

(* This is a first loop on the tetrahedra *)
For[i=1,i<22,i++, (* The tetrahedron is taken from the matrix T1. Here, a, b, c, and d
  are is vertices with the above convention: each vertex is
represented as the index of its row (between 1 and 14) in the
matrix 'Vertices'*)
a=T1[[i]][[1]];
b=T1[[i]][[2]];
c=T1[[i]][[3]];
d=T1[[i]][[4]];
(* A first subdeterminant is calculated. It is called 'Dx' and
  only depends on the tetrahedron. *)
Clear[Dx];Dx=Det[{Vertices[[a]]-Vertices[[c]],Vertices[[b]]-Vertices[[c]],Vertices[[d]]-Vertices[[c]]}];

(* This is a second, nested loop on the vertices *)
For[x=1,x<15,x++,
(* Only builds a conatraint for the vertices that are 
  not vertices of the tetrahedron *)
If[x≠a&&x≠b&&x≠c&&x≠d,
(* Prints which is the contraint that is being built *)
Print["Constraint ",k,"; tetrahedron:",i,"; vertex:",x];
(* A second subdeterminant is calculated. It is called 'Da' and
  depends on both the tetrahedron and the vertex. *)
Clear[Da];
Da=Det[{Vertices[[b]]-Vertices[[d]],Vertices[[c]]-Vertices[[d]],Vertices[[x]]-Vertices[[d]]}];
(* A third subdeterminant is calculated. It is called 'Db' and
  depends on both the tetrahedron and the vertex. *)
Clear[Db];
Db=Det[{Vertices[[a]]-Vertices[[d]],Vertices[[c]]-Vertices[[d]],Vertices[[x]]-Vertices[[d]]}];
(* A fourth subdeterminant is calculated. It is called 'Dc' and
  depends on both the tetrahedron and the vertex. *)
Clear[Dc];
Dc=Det[{Vertices[[a]]-Vertices[[d]],Vertices[[b]]-Vertices[[d]],Vertices[[x]]-Vertices[[d]]}];
(* A last subdeterminant is calculated. It is called 'Dd' and
  depends on both the tetrahedron and the vertex. *)
Clear[Dd];
Dd=Det[{Vertices[[a]]-Vertices[[c]],Vertices[[b]]-Vertices[[c]],Vertices[[x]]-Vertices[[c]]}];
(* Writes the non-zero coefficients of matrix A in row k. Those
 coeeficients depend on the sign of Dx. *)
If[Dx<=0,
A1=ReplacePart[A1,{k,a}->-Da];
A1=ReplacePart[A1,{k,b}->Db];
A1=ReplacePart[A1,{k,c}->-Dc];
A1=ReplacePart[A1,{k,d}->Dd];
A1=ReplacePart[A1,{k,x}->-Dx];,
A1=ReplacePart[A1,{k,a}->Da];
A1=ReplacePart[A1,{k,b}->-Db];
A1=ReplacePart[A1,{k,c}->Dc];
A1=ReplacePart[A1,{k,d}->-Dd];
A1=ReplacePart[A1,{k,x}->Dx];];
k++;
]
]
]

(* Matrix A2 of the inequation system A2.w2>0 is now built. Each 
 contraint of this system corresponds to a pair constituted of a
 tetrahedron (in matrix 'T2') and a vertex of the triangulation
 (in matrix 'Vertices'). Each constraint checks that the vertex,
   lifted by height function w2, is above the affine hull of the
  tetrahedron, lifted by the same height function in dimension 4.
 The number of constraints therefore the number of tetrahedra in
 T2 times the number of vertices of the triangulation that do not
 belong to the tetrahedron being checked (i.e. 10 vertices). As
 T2 contains 21 tetrahedra The system A2.w2>0 has 210 rows. Vector
 w2 has 14 rows (one per vertex). Hence A2 has 14 columns. The
 matrix 'A2' is first initialized with zeros and then filled with
 its actual coefficients. Those coefficients are subdeterminant of
 a determinant described by Takeushi (LCNS 2098, 330-338, 2001).
In the paper, one also finds a comprehensive description of the
 system A.w>0. *)
A2=ConstantArray[0,{210,14}];

(* k is the index on the row of matrix A that is being filled.
It is initialized as 1. *)  
k=1;

(* This is a first loop on the tetrahedra *)
For[i=1,i<22,i++, (* The tetrahedron is taken from the matrix T1. Here, a, b, c, and d
  are is vertices with the above convention: each vertex is
represented as the index of its row (between 1 and 14) in the
matrix 'Vertices'*)
a=T2[[i]][[1]];
b=T2[[i]][[2]];
c=T2[[i]][[3]];
d=T2[[i]][[4]];
(* A first subdeterminant is calculated. It is called 'Dx' and
  only depends on the tetrahedron. *)
Clear[Dx];Dx=Det[{Vertices[[a]]-Vertices[[c]],Vertices[[b]]-Vertices[[c]],Vertices[[d]]-Vertices[[c]]}];

(* This is a second, nested loop on the vertices *)
For[x=1,x<15,x++,
(* Only builds a conatraint for the vertices that are 
  not vertices of the tetrahedron *)
If[x≠a&&x≠b&&x≠c&&x≠d,
(* Prints which is the contraint that is being built *)
Print["Constraint ",k,"; tetrahedron:",i,"; vertex:",x];
(* A second subdeterminant is calculated. It is called 'Da' and
  depends on both the tetrahedron and the vertex. *)
Clear[Da];
Da=Det[{Vertices[[b]]-Vertices[[d]],Vertices[[c]]-Vertices[[d]],Vertices[[x]]-Vertices[[d]]}];
(* A third subdeterminant is calculated. It is called 'Db' and
  depends on both the tetrahedron and the vertex. *)
Clear[Db];
Db=Det[{Vertices[[a]]-Vertices[[d]],Vertices[[c]]-Vertices[[d]],Vertices[[x]]-Vertices[[d]]}];
(* A fourth subdeterminant is calculated. It is called 'Dc' and
  depends on both the tetrahedron and the vertex. *)
Clear[Dc];
Dc=Det[{Vertices[[a]]-Vertices[[d]],Vertices[[b]]-Vertices[[d]],Vertices[[x]]-Vertices[[d]]}];
(* A last subdeterminant is calculated. It is called 'Dd' and
  depends on both the tetrahedron and the vertex. *)
Clear[Dd];
Dd=Det[{Vertices[[a]]-Vertices[[c]],Vertices[[b]]-Vertices[[c]],Vertices[[x]]-Vertices[[c]]}];
(* Writes the non-zero coefficients of matrix A in row k. Those
 coeeficients depend on the sign of Dx. *)
If[Dx<=0,
A2=ReplacePart[A2,{k,a}->-Da];
A2=ReplacePart[A2,{k,b}->Db];
A2=ReplacePart[A2,{k,c}->-Dc];
A2=ReplacePart[A2,{k,d}->Dd];
A2=ReplacePart[A2,{k,x}->-Dx];,
A2=ReplacePart[A2,{k,a}->Da];
A2=ReplacePart[A2,{k,b}->-Db];
A2=ReplacePart[A2,{k,c}->Dc];
A2=ReplacePart[A2,{k,d}->-Dd];
A2=ReplacePart[A2,{k,x}->Dx];];
k++;
]
]
]

(* Initializes height functions w1 and w2. *)
w1={12,10000,12,10000,11,10,10,10,10,10,0,0,0,0};
w2={10000,12,10000,12,10,11,10,10,10,10,0,0,0,0};
(* Formally calculates the products A1.w1 and A2.w2 and stores them
  as vectors 'b1' and 'b2' respectively *)
b1=A1.w1;
b2=A2.w2;
(* the integer n1 will be the number of positive elements in vector b1
  and the integer n2 will be the number of positive elements in
  vector b2.  *)
n1=0;
n2=0;
(* counts the number of positive (formally) elements in vector b *)
For[i=1,i<211,i++,If[b1[[i]]>0,n1++];If[b2[[i]]>0,n2++];]
(* prints the result *)
Print["Number of positive element in vector b1: ",n1," out of 210"];
Print["Number of positive element in vector b2: ",n2," out of 210"];
If[n1==210&&n2==210,Print["Rudin's triangulation in 2-regular!"]]