(*This notebook is a companion to the paper "Infinitesimal rigidity of a compact hyperbolic 4-orbifold with totally geodesic boundary", by Tarik Aougab and Peter Storm.  It was prepared 2009-02-13 by the authors.  We attempted to follow the notation and description of the paper.  The notebook was written for Wolfram Mathematica 7.0, and may not run on old versions.  No complaining about slow run times!  Go have a coffee while the computer's running.*)

(*Define the Minkowski form, and a normalized version.*)
form[x_, y_] := x.y - 2 x[[1]] * y[[1]]
nform[x_,y_]:= form[x,y] / Sqrt[form[x,x] form[y,y] ]
perpQ[x_,y_]:=(FullSimplify[form[x,y] ] == 0)

(*Define the golden ratio and its reciprocal.*)
τ = (1 + Sqrt[5])/2;
τi = -(1 - Sqrt[5])/2;

(*Define the 96 space-like normal vectors (Subscript[n, i]) corresponding to the walls of the 120-cell not in the subset C.  See the beginning of Section 3.*)
Needs["Combinatorica`"]
e=IdentityMatrix[4];
A=AlternatingGroup[4]/.{1->e[[1]],2->e[[2]],3->e[[3]],4->e[[4]]};
pmLst=Map[Append[#,0]&,Tuples[{1,-1},3]];
n=Map[
Prepend[#[[1]] . (#[[2]] * {τ,1,τi,0}), Sqrt[2τ] ] &,
Tuples[{A,  pmLst}] ];
Clear[A,pmLst];(*clean up the namespace*)

(*Do some tests.*)
(*There should be 96 vectors in n.*)
Length[n]
(*Each vector should have norm 3 - Sqrt[5].*)
Union[Map[FullSimplify[form[#,#] ]&,n] ]
(*Each vector should be orthogonal to 9 others.  See Prop.3.1.*)
Union[Map[Function[x,Length[Select[n,perpQ[#,x]& ] ] ],n] ]
(*Space-like vectors should be pointing out of the 120-cell, meaning on distinct walls nform should never be positive.  If not orthogonal, distinct walls are disjoint.  This means nform should never take values between -1 and 0.  See equation (1).  So the following should output {0.}.*)
Union[Select[Map[N[nform@@#]&,Subsets[n,{2}] ], #>-1.0&] ]

(*Suppose we have an infinitesimal deformation v of n.  In the paper v is also denoted by n with a dot.*)
(*V has another meaning to Mathematica.  So unprotect it.*)
Unprotect[V];
V=Array[v,480];
(*It will be convenient to partition V.*)
Vp=Partition[V,5];
(*resort to a bit of procedural programming here*)
perpPairs= Reap[Do[
If[perpQ[n[[i]],n[[j]] ],Sow[{i,j}] ], {i,Length[n]-1},{j,i+1,Length[n]}] ][[2,1]];
(*Should be 432 of them.*)
Length[perpPairs]==432
(*Define the 432 linear polynomials p.*)
p=Map[form[Vp[[#[[1]] ]],n[[#[[2]] ]] ] + form[n[[#[[1]] ]], Vp[[#[[2]] ]] ]&,perpPairs];

Aalg =Map[D[#[[1]],{#[[2]]}]&,Outer[List,p,V],{2}];
Length[Aalg]
Length[Transpose[Aalg]]
Aalg[[234,145]]
α=Sqrt[1 + Sqrt[5]];
(*Every entry of Aalg should be in Q(α).*)
algNumQ[x_]:=Module[{tmp},
tmp=ToNumberField[x,α];
Element[tmp,Rationals] || Head[tmp]===AlgebraicNumber]
And@@Map[algNumQ,Flatten[Aalg] ]

(*Define function replacing a number in Q(α) with a 4x4 rational matrix.  First make a slightly modified ToNumberField fct.*)
toNumberField[x_]:=Module[{tmp},
tmp=ToNumberField[x,α];
If[Head[tmp]===AlgebraicNumber,tmp[[2]],{tmp,0,0,0}] ] 
explodeNumber[x_]:=Transpose[Map[toNumberField,x*{1,α,α^2,α^3}] ] (*Check the paper's example and a few others.*)
Map[MatrixForm[explodeNumber[#]]&,{-1/2 * (1 + Sqrt[5]),1/2*(1-Sqrt[5]),5/4,0}]

(*Do the same for an entire matrix with entries in Q(α).  This fct is slight hack, but OK.*)
explodeRow[row_]:=Module[{tmp},
tmp = Map[explodeNumber, row];
{Join@@(#[[1]]&/@tmp),Join@@(#[[2]]&/@tmp),Join@@(#[[3]]&/@tmp),Join@@(#[[4]]&/@tmp)}]
explodeMatrix[mat_]:=Join@@Map[explodeRow,mat]

(*Computing AQ is very slow.  Start the computation and take a break.*)
AQ=explodeMatrix[Aalg];
Map[Length, {AQ,Transpose[AQ]}]
And@@Map[NumberQ,Flatten[AQ] ]

(*Uncomment the following slow command if you'd like to compute a graphical representation of AQ.*)
(*ArrayPlot[AQ]*)

(*Let's see how big the denominators if entries of AQ are.  The output should be simply {1,2}.*)
Union[Denominator/@Flatten[AQ] ]

(*Multiply each row of AQ by the LCM of the denominators of that row.  This doesn't alter the rank.*)
AZ=Map[# * LCM@@(Denominator/@#)&, AQ];
(*Verify that AZ is a matrix of integers.*)
And@@(IntegerQ/@Flatten[AZ])
(*Check that AZ has the correct dimensions.*)
Length/@{AZ,Transpose[AZ]}
AZ[[1323,768]]

(*Save AZ to a text file, in case you'd like to analyze it with other software.*)
Save["AZmatrix",AZ]

(*As in the final paragraph of the paper, compute the rank of the matrix AZ modulo the prime 3.  The answer should be 1496.*)MatrixRank[AZ,Modulus->3]

(*Finally, let's compute the rank over the integers.  On my laptop it took about 12 seconds, but this computation may take a long time.*)
MatrixRank[AZ]