# file name: lines.sage
# compute the number of lines on a cubic surface
# usage: sage lines.sage 2 3 5 7

import sys

K = QQ
R.<a,b,c,d> = PolynomialRing(K, order = "lex")
S.<x,y,z,w> = PolynomialRing(K)


# moduli parameters:
p1, p2, p3, p4 = map( lambda x: int(x), sys.argv[1:])

# define a cubic form F:
F = x^3 + y^3 + z^3 + w^3 + p1*y*z*w + p2*x*z*w + p3*x*y*w + p4*x*y*z

J = ideal(F.jacob())
pd = J.dimension() - 1
print 
print "F =", F
print "Projective dimension of singular locus =", pd

# Set up a parametrized line L(t) depending on four 
# indeterminates.  The indeterminates a, b, c, d
# are parameters for the Grassmanian of lines in P^3, G(2,4):
L = lambda t: (1,t,a+b*t,c+d*t)

# compose F with L(t):
FL = lambda t: F(L(t))

# The ideal J defines the set of points X in an 
# affine cell of G(2,4) corresponding to lines
# which lie on the surface F = 0.  This is because
# a line which intersects a cubic surface in four
# points lies on the cubic surface.
J = ideal( FL(0), FL(1), FL(2), FL(3) )

# Compute a Groebner basis of J
B = J.groebner_basis()

# Check to see on what variables the elements
# of the Groebner basis depend:
print "Variables on which g[i] depends:"
for i in range(0,len(B)):
  print i, B[i].variables()

# Sure enough, the polyonmial g = B[0] is univariate.
# It cuts out the projection of the set X
# onto the d-axis.  Its degree is the 
# degree of the projection, which is generically
# the degree of X.

g = B[0]; dg = g.degree(); nf = len(factor(g))
print "degree(g) =", g.degree()
print "number of factors of g =", len(factor(g))
if nf == 1 and dg == 27:
	print "cubic surface has 27 lines"

print