> restart: with(ListTools): with(LinearAlgebra): with(Groebner):

> Q:=Matrix([[x[1,1],x[1,2],x[1,3],x[1,4]],[x[1,2],x[2,2],x[2,3],x[2,4]],[x[1,3],x[2,3],x[3,3],x[3,4]],[x[1,4],x[2,4],x[3,4],x[4,4]]]):

> U:=Transpose(Matrix([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[1,0,-4,-3],[1,1,2,3]])):

> V:=Transpose(Matrix([[0,-5,2,-2],[5,0,-2,1],[-1,0,0,0],[-1,5,-3,0],[-17,4,-2,-3],[-5,4,2,-1]])):

> for t from 1 to 6 do

> NewU[1..4,t]:=U[1..4,t]/VectorNorm(U[1..4,t],2):

> NewV[1..4,t]:=V[1..4,t]/VectorNorm(V[1..4,t],2):

> end do:

> L:=[seq(simplify(Transpose(NewU[1..4,k]).Q.NewU[1..4,k]+Transpose(NewV[1..4,k]).Q.NewV[1..4,k]),k=1..6)]:

> S:=combinat[choose](4,3):

> mjk:=[seq(seq(Determinant(SubMatrix(Q,r,c)),r=S),c=S)]:

> Eq:=[op(L),op(mjk)]:

> Gb:=Basis(Eq,plex(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[3,4],x[4,4])):

> f0:=subs({x[3,4]=1,x[4,4]=y},Gb[1]):

> s := sturmseq(f0,y):

> NumberZero:=sturm(s, y, -infinity, +infinity):

> printf(" The number of nonzero real roots of f_0 is %d \n", NumberZero):

> for j from 1 to 4 do

> for k from j to 4 do

> NewEq:=[op(Eq),x[j,k]-1,x[3,4],x[4,4]]:

> NewGb:=Basis(NewEq,plex(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[3,4],x[4,4])):

> if NewGb[1]=1 then

> printf(" If x[%d,%d]=1 then the polynomial system has no roots \n", j,k):

> end if

> end do:

> end do:

The number of nonzero real roots of f_0 is 0

If x[1,1]=1 then the polynomial system has no roots

If x[1,2]=1 then the polynomial system has no roots

If x[1,3]=1 then the polynomial system has no roots

If x[1,4]=1 then the polynomial system has no roots

If x[2,2]=1 then the polynomial system has no roots

If x[2,3]=1 then the polynomial system has no roots

If x[2,4]=1 then the polynomial system has no roots

If x[3,3]=1 then the polynomial system has no roots

If x[3,4]=1 then the polynomial system has no roots

If x[4,4]=1 then the polynomial system has no roots