-- This is Macaulay2 code from the paper "A small frame and certificate of its injectivity" (arXiv:1502.04656) 

needsPackage "realroots2"
-- If this doesn't work, some useful commands are 
-- installPackage "realroots2"  
-- and 
-- loadPackage "realroots2"


-- First we define the ring of 16 variables we'll be working in
S = QQ[X_(1,1),X_(1,2), X_(1, 3), X_(1, 4), X_(2,2), X_(2, 3), X_(2, 4), X_(3,3),  X_(4,4), Y_(1,2), Y_(1, 3), Y_(1, 4), Y_(2, 3), Y_(2, 4),X_(3,4), Y_(3,4)];

-- In order to write our Hermitian matrix, we temporarily work in the field extension QQ(i)
kk = frac(QQ[I]/ideal(I^2+1))
R = kk[X_(1,1),X_(1,2), X_(1, 3), X_(1, 4), X_(2,2), X_(2, 3), X_(2, 4), X_(3,3),  X_(4,4), Y_(1,2), Y_(1, 3), Y_(1, 4), Y_(2, 3), Y_(2, 4),X_(3,4), Y_(3,4)];
con = map(R,R,{I=>-I})

-- Our 4x4 Hermitian matrix 
Q = matrix(R,
    {{X_(1,1), X_(1, 2) + I*Y_(1, 2), X_(1, 3) + I*Y_(1, 3), X_(1, 4) + I*Y_(1, 4)}, 
    {X_(1, 2) - I*Y_(1, 2), X_(2,2), X_(2, 3) + I*Y_(2, 3), X_(2, 4) + I*Y_(2, 4)}, 
    {X_(1, 3) - I*Y_(1, 3), X_(2, 3) - I*Y_(2, 3), X_(3,3), X_(3, 4) + I*Y_(3, 4)}, 
    {X_(1, 4) - I*Y_(1, 4), X_(2, 4) - I*Y_(2, 4), X_(3, 4) - I*Y_(3, 4),X_(4,4)}})

-- We take the ideal of 3x3 minors and lift it to the original ring S (with rational coefficients)
rk2R = minors(3,Q);
rk2 = ideal(sub(mingens rk2R, S));
rk2R==sub(rk2,R)

-- Now we define our frame
V= {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, 
     {1, 9*I, -5 - 7*I, -6 - 7*I}, 
     {1, 1 - I, -5 - 2*I, -1 - 8*I}, 
     {1, -2 + 4*I, -4 - 2*I, 3 + 8*I}, 
     {1, -3 + I, 1 - 8*I, 7 - 6*I}, 
     {1, 3 - 3*I, -8 + 7*I, -6 - 2*I}, 
     {1, -3 + 5*I, 5 + 6*I, 2*I}, 
     {1, -3 + 8*I, 5 - 5*I, -6 - 4*I}};


-- and turn the frame into linear equations on the matrix Q
LinR = V/(v->(matrix{{v_(0),con(v_(1)),con(v_(2)),con(v_(3))}}*Q*transpose(matrix{v}))_(0,0));
use S
L = ideal(sub(matrix{LinR},S));

-- Joining the rank-2 condition with our linear space gives the solution set of interest
pts = rk2 + L;

-- We check that the solution set has dimension 1 and degree 20. 
dim pts
degree pts 

-- We then eliminate all but two of the variables to obtain the polynomial f in X_(3,4), Y_(3,4) vanishing on our original solutions
use S
projPts = eliminate({X_(1,1), X_(1,2), X_(1, 3), X_(1, 4), X_(2,2), X_(2, 3), X_(2, 4), X_(3,3),  X_(4,4), Y_(1,2), Y_(1, 3), Y_(1, 4), Y_(2, 3), Y_(2,4)}, pts);
f = (mingens projPts)_(0,0)

-- Then we make the substitution x=X_(3,4) and Y_(3,4)=1 to obtain a univariate polynomial F in the variable x
T = QQ[x]
mm = map(T,S,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,x,1})
F = mm(f)

-- Now we use Sturm sequences (implemented in the realroots package) to check that F has no real roots
numRealSturm(F)


-- Our last computation checks that the ideal pts has no non-zero solutions with Y_(3,4)=0. 
--We do this by adjoining Y_(3,4) and X-1 to our equations for each variable X=X_(j,k), Y_(j,k_ and computing a Groebner basis.  
--For each variable, we find a Groebner basis 1, certifying that there are no such solutions.  
xyvar = flatten entries vars S;
xyvar/(X->gens(gb(pts + ideal(Y_(3,4), X-1))))