sage: #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-20 by the authors.  We attempt to follow the notation and description of the paper.  The notebook was written using Sage 3.2.3.
sage: #Define the Minkowski form, and a normalized version which outputs floats.
sage: def form(x,y): return -x[0]*y[0] + x[1]*y[1] + x[2]*y[2] + x[3]*y[3] + x[4]*y[4]
sage: def nform(x,y): 
...       xn = x.n()
...       yn = y.n()
...       return (form(xn,yn)/sqrt(form(xn,xn)*form(yn,yn)))
...
sage: def perp_q(x,y): return (form(x,y)==0)
sage: #Define the golden ratio and its reciprocal in the number field K.
sage: K.<alpha>=NumberField((x^2 - 1)^2 - 5,'alpha',embedding=1.79)
sage: tau = alpha^2 /2; taui = 1/2 * alpha^2 - 1
sage: #Define the 96 space-like normal vectors n_i corresponding to the walls of the 120-cell not in the subset C.  See the beginning of Section 3.
sage: A=[x.matrix() for x in AlternatingGroup(4)]
sage: pm_lst=[vector(K,[x,y,z,0]) for x in [-1,1] for y in [-1,1] for z in [-1,1]]
sage: n_short=[vector(K,[tau, 1, taui, 0]).pairwise_product(x) * y for y in A for x in pm_lst]
sage: n=[vector(K,[alpha] + x.list()) for x in n_short]
sage: len(n)
96
sage: #Each vector of n should have norm 3 - sqrt(5).
sage: set(form(x,x) for x in n)
set([-alpha^2 + 4])
sage: bool(-alpha^2 + 4 == 3 - sqrt(5))
True
sage: #Each vector of n should be Minkowski-orthogonal to 9 others.  See Prop.3.1.
sage: set(len([x for x in n if (form(x,y)==0)]) for y in n)
set([9])
sage: #The vectors of n should be points out of the 120-cell, meaning on distinct walls nform should never be positive.  If not orthogonal, walls are disjoint.  The means nform should not take values between -1 and 0.  See equation (1) of the paper.  Show the following should output floats very close to 0.0.  Warning: This command could be slow on your computer.
sage: n_pairs = [(n[j],n[k]) for j in range(len(n)-1) for k in range(j+1,len(n))]
sage: time set(nform(x,y) for (x,y) in n_pairs if (nform(x,y) > -1.0))
set([0.000000000000000, -1.45330080678069e-16, 2.90660161356137e-16, 2.90660161356137e-16])
Time: CPU 37.61 s, Wall: 37.95 s
sage: #Suppose we have an infinitesimal deformation v of n.  In the paper v is also denoted by n with a dot.
sage: R = PolynomialRing(K,480,"v")
sage: var('v'); v=R.gens()
sage: v_part = [v[k:k+5] for k in range(0,len(n)*5,5)]
sage: perp_pairs = [[j,k] for j in range(len(n)-1) for k in range(j+1,len(n))\
...       if perp_q(n[j],n[k])]
...
sage: len(perp_pairs)
432
sage: #The 432 linear polynomials p are satisfied by the infinitesimal deformation v.
sage: p=[form(v_part[j],n[k]) + form(n[j],v_part[k]) for [j,k] in perp_pairs]
sage: #Finally define the matrix of algebraic numbers Aalg by peeling off the correct coefficient of each p.  Unlike in Mathematica, here we compute the matrix rank of Aalg directly in the number field K.  The answer should be 374.
sage: M=MatrixSpace(K,432,480)
sage: Aalg = M([poly.monomial_coefficient(x) for poly in p for x in v])
sage: time Aalg.rank()
374
Time: CPU 68.80 s, Wall: 69.18 s