# INPUT: A cayley octad in P^3(QQ): The first five points have to be in standard position. The two points o6, o7 can be set by the user. The eighth point is then computed automatically. # OUTPUT: A quartic, its 36 linear matrix representations, and its 63 Gram matrices kk=QQ #The first 5 points must be as follows for the formula for o8 to work correctly: o1=matrix([[1,0,0,0]]); o2=matrix([[0,1,0,0]]); o3=matrix([[0,0,1,0]]); o4=matrix([[0,0,0,1]]); o5=matrix([[1,1,1,1]]); # ACTUAL INPUT: First entries of o6 and o7 should be 1 o6=matrix([[1,3,5,-2]]); o7=matrix([[1,4,9,16]]); # These seven points determine a net of quadrics, all of whom vanish at 8 points. Here is a formula for the eigth point (found by Bernd) o8=matrix([[(o7[0,2]*o6[0,1]-o7[0,2]*o6[0,3]-o7[0,1]*o6[0,2]+o7[0,1]*o6[0,3]-o6[0,1]*o7[0,3]+o6[0,2]*o7[0,3])/(o7[0,2]*o7[0,1]*o6[0,1]*o6[0,3]-o7[0,2]*o7[0,1]*o6[0,2]*o6[0,3]+o7[0,2]*o6[0,1]*o6[0,2]*o7[0,3]-o7[0,2]*o6[0,1]*o6[0,3]*o7[0,3]-o7[0,1]*o6[0,1]*o6[0,2]*o7[0,3]+o7[0,1]*o6[0,2]*o6[0,3]*o7[0,3]),(o7[0,2]*o6[0,3]-o6[0,2]*o7[0,3]-o7[0,2]+o6[0,2]-o6[0,3]+o7[0,3])/(o7[0,2]*o6[0,2]*o6[0,3]-o7[0,2]*o6[0,2]*o7[0,3]+o7[0,2]*o6[0,3]*o7[0,3]-o6[0,2]*o6[0,3]*o7[0,3]-o7[0,2]*o6[0,3]+o6[0,2]*o7[0,3]),(o7[0,1]*o6[0,3]-o6[0,1]*o7[0,3]-o7[0,1]+o6[0,1]-o6[0,3]+o7[0,3])/(o7[0,1]*o6[0,1]*o6[0,3]-o7[0,1]*o6[0,1]*o7[0,3]+o7[0,1]*o6[0,3]*o7[0,3]-o6[0,1]*o6[0,3]*o7[0,3]-o7[0,1]*o6[0,3]+o6[0,1]*o7[0,3]),(o7[0,2]*o6[0,1]-o7[0,1]*o6[0,2]-o7[0,2]+o7[0,1]-o6[0,1]+o6[0,2])/(o7[0,2]*o7[0,1]*o6[0,1]-o7[0,2]*o7[0,1]*o6[0,2]+o7[0,2]*o6[0,1]*o6[0,2]-o7[0,1]*o6[0,1]*o6[0,2]-o7[0,2]*o6[0,1]+o7[0,1]*o6[0,2])]]) o=[o1,o2,o3,o4,o5,o6,o7,o8] R.=PolynomialRing(kk) # Now we find the net of quadrics vanishing at these 8 points. This can be simplified a bit since we have included the unit vectors in the octad (the diagonal entries of the matrices corresponding to these quadrics will all be zero). A1=transpose(o1)*o1; A2=transpose(o2)*o2; A3=transpose(o3)*o3; A4=transpose(o4)*o4; A5=transpose(o5)*o5; A6=transpose(o6)*o6; A7=transpose(o7)*o7; A8=transpose(o8)*o8; Q=matrix([[A1[0,0],2*A1[0,1],2*A1[0,2],2*A1[0,3],A1[1,1],2*A1[1,2],2*A1[1,3],A1[2,2],2*A1[2,3],A1[3,3]], [A2[0,0],2*A2[0,1],2*A2[0,2],2*A2[0,3],A2[1,1],2*A2[1,2],2*A2[1,3],A2[2,2],2*A2[2,3],A2[3,3]], [A3[0,0],2*A3[0,1],2*A3[0,2],2*A3[0,3],A3[1,1],2*A3[1,2],2*A3[1,3],A3[2,2],2*A3[2,3],A3[3,3]], [A4[0,0],2*A4[0,1],2*A4[0,2],2*A4[0,3],A4[1,1],2*A4[1,2],2*A4[1,3],A4[2,2],2*A4[2,3],A4[3,3]], [A5[0,0],2*A5[0,1],2*A5[0,2],2*A5[0,3],A5[1,1],2*A5[1,2],2*A5[1,3],A5[2,2],2*A5[2,3],A5[3,3]], [A6[0,0],2*A6[0,1],2*A6[0,2],2*A6[0,3],A6[1,1],2*A6[1,2],2*A6[1,3],A6[2,2],2*A6[2,3],A6[3,3]], [A7[0,0],2*A7[0,1],2*A7[0,2],2*A7[0,3],A7[1,1],2*A7[1,2],2*A7[1,3],A7[2,2],2*A7[2,3],A7[3,3]], [A8[0,0],2*A8[0,1],2*A8[0,2],2*A8[0,3],A8[1,1],2*A8[1,2],2*A8[1,3],A8[2,2],2*A8[2,3],A8[3,3]]]); # the column span of L is the net of quadrics and we choose a basis L1, L2, L3 (in 4x4 mtrix form) L=transpose(kernel(transpose(Q)).matrix()) L1=matrix([[L[0,0],L[1,0],L[2,0],L[3,0]],[L[1,0],L[4,0],L[5,0],L[6,0]],[L[2,0],L[5,0],L[7,0],L[8,0]],[L[3,0],L[6,0],L[8,0],L[9,0]]]) L2=matrix([[L[0,1],L[1,1],L[2,1],L[3,1]],[L[1,1],L[4,1],L[5,1],L[6,1]],[L[2,1],L[5,1],L[7,1],L[8,1]],[L[3,1],L[6,1],L[8,1],L[9,1]]]) L3=matrix([[L[0,2],L[1,2],L[2,2],L[3,2]],[L[1,2],L[4,2],L[5,2],L[6,2]],[L[2,2],L[5,2],L[7,2],L[8,2]],[L[3,2],L[6,2],L[8,2],L[9,2]]]) R.=PolynomialRing(kk) # N is the first linear matrix representation and det(N) will be the quartic we work with N=x*L1+y*L2+z*L3 f=det(N) f=f.base_extend(R) print("The quartic is") print(f) print("") Grad=ideal(f,diff(f,x),diff(f,y),diff(f,z)) if Grad.dimension()!=0: sys.exit("The quartic f is not smooth! Choose a more generic octad.") Z=block_matrix([o1,o2,o3,o4,o5,o6,o7,o8],nrows=8,subdivide=False) # From this we can make a 8x8 symmetric matrix, D, of the 28 bitangents (diagonal entries of D will be 0) D=Z*N*transpose(Z) # For each of the subsets {i,j,k,l} of {0,..,7} we can make a linear change of coordinates T so that o[i]=e1,o[j]=e2,o[k]=e3,o[l]=e4. Then by permuting the entries of T^tNT we get a inequivalent LDR of the quartic f. It turns out that a set {i,j,k,l} and its complement give equivalent LDRs, giving 35 new LDRs total. Thus is suffices to look at 4-subsets of {0..7} containing 0, given by 3-subsets of {1..7}. Subs=Subsets([1,2,3,4,5,6,7],3).list(); LMR=[[[],N]]; for S in Subs: T = block_matrix([o[0],o[S[0]],o[S[1]],o[S[2]]],nrows=4,subdivide=False) NN = T*N*transpose(T) M = matrix([[0,NN[2,3],NN[1,3],NN[1,2]],[NN[2,3],0,NN[0,3],NN[0,2]],[NN[1,3],NN[0,3],0,NN[0,1]],[NN[1,2],NN[0,2],NN[0,1],0]]) if det(M)!=det(T)^2*det(N): print "trouble at";print S LMR=LMR+[[S,M]] print("The 36 LMRs are stored in 'LMR'.") QuadraticMonomials=[x^2,x*y,x*z,y^2,y*z,z^2] QuarticMonomials=[x^4,x^3*y,x^3*z,x^2*y^2,x^2*y*z,x^2*z^2, x*y^3,x*y^2*z,x*y*z^2,x*z^3,y^4,y^3*z,y^2*z^2,y*z^3,z^4] Gram=[] # Produce a list L that contains all 63 Steiner complexes. This is a (very) nested list: Each entry consists of six pairs of bitangents, given as six pairs of pairs of indices, referring to the entries in the matrix D StComplexes=[] for j in [0..7]: for i in [0..j-1]: L0=[] for k in [0..7]: if k!=i and k!=j: L0=L0+[[[i,k],[j,k]]] StComplexes=StComplexes+[L0] for k in [3..7]: for j in [2..k-1]: for i in [1..j-1]: J=[1..7]; J.remove(i); J.remove(j); J.remove(k); StComplexes=StComplexes+[ [ [[0,i],[j,k]],[[0,j],[i,k]],[[0,k],[i,j]],[[J[0],J[1]],[J[2],J[3]]],[[J[0],J[2]],[J[1],J[3]]],[[J[0],J[3]],[J[1],J[2]]] ] ] # Now we will run through all the Steiner complexes and produce the corresponding quadratic determinantal representation. QuadraticMonomials=[x^2,x*y,x*z,y^2,y*z,z^2] Gram=[] for Duple in [0..62]: St=StComplexes[Duple] #We are going to work in the P^5 of conics: #First produce the coefficient vectors C0,...,C5 of the six reducible conics corresponding to the six pairs of bitangents in the Steiner complex StBT=[] for i in range(6): StBT=StBT+[map((D[St[i][0][0],St[i][0][1]]*D[St[i][1][0],St[i][1][1]]).monomial_coefficient,QuadraticMonomials)] C0=matrix(StBT[0]); C1=matrix(StBT[1]); C2=matrix(StBT[2]); C3=matrix(StBT[3]); C4=matrix(StBT[4]); C5=matrix(StBT[5]) #These sit in a big matrix StMatrix=block_matrix([C0,C1,C2,C3,C4,C5],nrows=6,subdivide=False) #C0,...,C5 span a 2-dimensional plane in P^5 Plane=block_matrix([C0,C1,C2],nrows=3,subdivide=False) # We dutifully check whether that is actually true if not rank(StMatrix)==3: sys.exit("ERROR: Six points that were expected to lie on a conic actually don't. We give up.") # We extend the basis of our plane by a basis of the orthogonal complement and compute the coordinates of the points C3,C4,C5 in this new basis; (C0,C1,C2 have of course become the first three unit vectors) Orth=kernel(transpose(Plane)).matrix() Coord=transpose(block_matrix([Plane,Orth],nrows=2,subdivide=False)) InverseCoord = Coord.inverse() B3=InverseCoord*transpose(C3) B4=InverseCoord*transpose(C4) B5=InverseCoord*transpose(C5) # We "know" that C0,...,C5 lie on a unique conic in our plane. We compute the matrix of that conic in the basis B0,B1,B2: FirstConic=matrix([[(B3*transpose(B3))[0,1],(B3*transpose(B3))[0,2],(B3*transpose(B3))[1,2]],[(B4*transpose(B4))[0,1], (B4*transpose(B4))[0,2],(B4*transpose(B4))[1,2]],[(B5*transpose(B5))[0,1],(B5*transpose(B5))[0,2],(B5*transpose(B5))[1,2]]]) SecondConic=kernel(transpose(FirstConic)).matrix() Conic=matrix([[0,SecondConic[0,0],SecondConic[0,1]],[SecondConic[0,0],0,SecondConic[0,2]],[SecondConic[0,1],SecondConic[0,2],0]]) # The dual will give us the Gram matrix we want: ConicDual=Conic.adjoint() # ConicDual is a non-degenerate quadratic form on our plane. Extending that form by 0 on the orthogonal complement gives the rank-3-quadratic form corresponding to the Gram matrix: ZeroMatrix=matrix([[0,0,0],[0,0,0],[0,0,0]]) BigConic=block_diagonal_matrix([ConicDual,ZeroMatrix],subdivide=False) # We just have to undo the coordinate change back to the usual monomial basis: GramSt=Coord*BigConic*transpose(Coord) mon=matrix(QuadraticMonomials) TestPol=(mon*GramSt*transpose(mon))[0,0] for m in QuarticMonomials: if f.monomial_coefficient(m)!=0: GramSt=(f.monomial_coefficient(m)/(TestPol.monomial_coefficient(m)))*GramSt TestPol=(f.monomial_coefficient(m)/(TestPol.monomial_coefficient(m)))*TestPol break # ...and it better be a multiple of our original quartic det(N): if TestPol!=f: print "WARNING: Something went wrong. The Gram matrix found for the Steiner complex "; print St; print "is incorrect!" # All went well. We produce a list whose entries are 63 pairs of a Steiner complex and the corresponding Gram matrix! Gram=Gram+[GramSt] print "The 63 Gram matrices are stored in 'Gram'."