#####################################################################################################
# This file contains a function designed to run under Maple(TM) that finds the lines in the complex 
# 3-space covered by a given parametrization of a rational surface. To be used, one must run before:
#           > read "/path_to_file/Lines_in_surface_public.mpl":
#
# Let X = [x_1(t,s), x_2(t,s), x_3(t,s)] be a list of rational functions with rational coefficients. 
# Let w be a variable. Then
#           > LinesSurfaceParam(X,s,t,w); 
# returns a list [r_1,. . ., r_n] such that any entry has 
# the form [aw+b,cw+d,ew+f], where a, b, c, d, e and f are algebraic numbers. Any line in the complex
# 3-space having infinite points in the image of the complex plane by X, has a parametrization among
# the r_i.
#
# Authors: Juan G. Alcazar and Jorge Caravantes
#
# Reference: Juan G. Alcazar and Jorge Caravantes. "On the computation of the straight lines 
# contained in a rational surface". arXiv:1603.03959 [math.AG]
# 
# License: Creative commons
#####################################################################################################


LinesSurfaceParam:=proc(X,s,t,w)
# X is the parametrization (three rational functions in the variables s, t).
# w is the variable we are using for the lines
local X_s,X_t,X_ss,X_st,X_tt, #partial derivatives of the parametrization
   E, F, G, # Coefficients of the 1st fundamental form
   E_s, E_t, F_s, F_t, G_s, G_t, # partial derivatives of E, F, G
   Gamma_111, Gamma_112, Gamma_121, # Christoffel symbols
   Gamma_122, Gamma_221, Gamma_222,  # moreChristoffel symbols
   e_, f_, g_,  # Coefficients of the 2nd fundamental form
   n_, c, ct,   # normal vector to the surface, gcd of all the coordinates of n_, relevant part of its content.
   M_, delta2M,  #Condition to be asymptotic, its content wrt 'w'
   N, delta2N,   #Condition to be geodesic combined with M_w, its content wrt 'w'
   chi, delta1, delta2, delta3, delta4, eta,   # Those of the paper
   cont, # content of a polynomial
   gamma, cc,    # dw/dt, gcd of all coordinates of X_s and X_t
   den, den0,     # denominator of the parametrization, specialized denominator
   rectas_candidatas,    #Output
   lista,  #  [chi, delta1, delta2, delta3]
   lista2, factores, #two more lists of polynomials with factors of the elements of lista
   l, l0, l0_s,   #a member of lista, specialization of the member, partial derivative
   l_t0,           # specialization of partial derivative
   a0, a,     # 1st coordinates od pts, 1st coordinate of a pt, all in parameter plane
   p0, p, v0,  #  2nd coordinates of pts, 2nd coord of a pt, tangent vectors, all in parameter plane
   P, V,    # Image point, tangent vector to the surface at P
   pol,     # condition to find a line
   i, j, rpt;	# Aux variables
    # We begin by converting the data to what we need:
    # Coefficients of the first fundamental form and their derivatives:
den := lcm(denom(X[1]),lcm(denom(X[2]),denom(X[3]))):
X_s := simplify(diff(X,s)): 
X_t := simplify(diff(X,t)):
X_ss := simplify(diff(X_s,s)):
X_st := simplify(diff(X_s,t)):
X_tt := simplify(diff(X_t,t)):
E := simplify(X_t[1]^2+X_t[2]^2+X_t[3]^2):
E_s := simplify(diff(E,s)):
E_t := simplify(diff(E,t)):
F := simplify(X_t[1]*X_s[1]+X_t[2]*X_s[2]+X_t[3]*X_s[3]):
F_s := simplify(diff(F,s)):
F_t := simplify(diff(F,t)):
G := simplify(X_s[1]^2+X_s[2]^2+X_s[3]^2):
G_s := simplify(diff(G,s)):
G_t := simplify(diff(G,t)):
    # Christoffel's symbols:
Gamma_111 := simplify( (G*E_t-2*F*F_t+F*E_s)/(2*(E*G-F^2)) ): 
Gamma_112 := simplify( (2*E*F_t-E*E_s-F*E_t)/(2*(E*G-F^2)) ):
Gamma_121 := simplify( (G*E_s-F*G_t)/(2*(E*G-F^2)) ):
Gamma_122 := simplify( (E*G_t-F*E_s)/(2*(E*G-F^2)) ):
Gamma_221 := simplify( (2*G*F_s-G*G_t-F*G_s)/(2*(E*G-F^2)) ):
Gamma_222 := simplify( (E*G_s-2*F*F_s+F*G_t)/(2*(E*G-F^2)) ):
    # We now compute the normal vector, but without normaiz
n_ := simplify(LinearAlgebra[CrossProduct](X_t,X_s)):
    #We now remove common factors of all the coordinates of n_
c := gcd(numer(n_[1]),gcd(numer(n_[2]),numer(n_[3]))):
n_ := [seq(simplify(n_[i]/c),i=1..3)]:
    # We now compute the coefficiets of the second fundamental form, up to such normalization
e_ := simplify(X_tt[1]*n_[1]+X_tt[2]*n_[2]+X_tt[3]*n_[3]):
f_ := simplify(X_st[1]*n_[1]+X_st[2]*n_[2]+X_st[3]*n_[3]):
g_ := simplify(X_ss[1]*n_[1]+X_ss[2]*n_[2]+X_ss[3]*n_[3]):
# Now let us try vertical lines in the parameter space
eta := simplify(gcd(numer(g_),numer(Gamma_221))):
if eta = 0 then error "The surface is ruled": fi:
cont := content(eta,s):
cont := simplify(cont/gcd(cont,diff(cont,t))):
ct := gcd(cont,c):#maybe this is not necessary
cc := gcd(gcd(numer(X_s[1]),numer(X_s[2])),numer(X_s[3])):
cc := gcd(gcd(gcd(numer(X_t[1]),numer(X_t[2])),numer(X_t[3])),cc):
cont := simplify(numer(simplify(cont/(den*cc)))):
a0 := [solve(cont,t)]:
rectas_candidatas := []:
for a in a0 do
  p := 0:
  if simplify(subs(t=a,X_s))<>[0,0,0] then
    while subs(t=a,s=p,den)=0 or subs(t=a,s=p,[X_s[1],X_s[2],X_s[3]])=[0,0,0] do
      p:=p+1: 
	  od:
    P := subs(t=a,s=p,X):
    V := subs(t=a,s=p,[X_s[1],X_s[2],X_s[3]]):
    rectas_candidatas := [op(rectas_candidatas), [P[1]+V[1]*w,P[2]+V[2]*w,P[3]+V[3]*w]]:
	fi:
  od:  
# Now singular lines given by ct
# First we remove lines that contract to a point:
ct := numer(simplify(ct/gcd(numer(X_s[1]),gcd(numer(X_s[2]),numer(X_s[3]))))):
a0 := solve(ct,t):
for a in a0 do
  p := 0:
  while subs(t=a,s=p,X_s)=[0,0,0] do p := p+1: od:
  P := subs(t=a,s=p,X):
  V := subs(t=a,s=p,[X_s[1],X_s[2],X_s[3]]):
  pol := simplify(LinearAlgebra[CrossProduct](simplify(subs(t=a,X)-P),V)):
  if pol =[0,0,0] then
    rectas_candidatas := [op(rectas_candidatas), [P[1]+V[1]*w,P[2]+V[2]*w,P[3]+V[3]*w]]:
	fi:
  od:
# We go now for other lines:
    #First equation:
M_ := numer(simplify(e_+2*f_*w+g_*w^2)):
delta2M := content(M_,w):#gcd(coeff(M_,w,2),gcd(coeff(M_,w,1),coeff(M_,w,0))):
M_ := simplify( M_/delta2M ):
    # Second equation:
gamma := simplify(Gamma_221*w^3+(2*Gamma_121-Gamma_222)*w^2+(Gamma_111-2*Gamma_122)*w-Gamma_112):
N := numer( simplify(diff(M_,w)*gamma + diff(M_,s)*w + diff(M_,t))):
delta2N := content(N,w):#gcd(gcd(gcd(coeff(N,w,3),coeff(N,w,2)),coeff(N,w,1)),coeff(N,w,0)): 
N := simplify( N/delta2N ):
chi := simplify(resultant(M_,N,w)):
# WARNING: Maybe, as it is now, the algorithm won't detect ruled surfaces
# when the lines are vertical in the parameter plane.
if chi = 0 then
  error "The surface is ruled":
  fi:
    # removing denominators:
delta1 := gcd(coeff(diff(M,w),w,1),coeff(diff(M,w),w,0)):
delta2 := delta2M*delta2N:
delta3 := numer(e_*g_-f_^2):
delta4 := numer(n_[1]^2+n_[2]^2+n_[3]^2):
lista:=[chi,delta1,delta2,delta3,c]:
# Now we eliminate redundant info:
# This means we get pairwise coprime square free polynomials 
for i from 1 to nops(lista) do 
  if lista[i]<>0 then 
    l := lista[i]:
    l := simplify(l/gcd(l,diff(l,s))):
    lista[i] := numer (simplify( l/den )):
  else
    lista[i] :=1:  
    fi:
  od:
for i from 1 to nops(lista)-1 do  
  for j from i+1 to nops(lista) do
    lista[i]:=simplify(lista[i]/gcd(lista[j],lista[i])):
	od:
  od:
# It is preferable to treat the factors of the polynomials in lista, since they are more simple.
lista2 := []:
for l in lista do
  factores := evala(AFactors(l))[2]:
  lista2 := [op(lista2), seq(factores[i][1],i=1..nops(factores))]:
  od:
# To finish, we remove curves consisting of singular points
for i from 1 to nops(lista2) do
  lista2[i]:=numer(simplify(lista2[i]/cc)):
  od:
# Now we can compute lines
for l in lista2 do
# l is squarefree. Now we must find a vertical line without critical points of l
# also this line must be such that all pints in d have finite image via the parametrization
  if degree(l)<>0 then
    a:=0:
    l0 := simplify(subs(t=a, l)):
    den0 := simplify(subs(t=a,den)):
    l0_s := diff(l0,s):
	rpt:=false:
	if resultant(l0,l0_s,s)=0 or resultant(l0,den0,s)=0 then
	  rpt:=true:
# Now we get all the points of l in the critical line x=a
    elif degree(l0)=0 then
	  rpt:=true:
	else
	  p0 := [solve(l0,s)][1]: #we just need one solution
	  if subs(t=a,s=p0,X_s)=[0,0,0] or subs(t=a,s=p0,X_t)=[0,0,0] then
	    rpt:=true:
		fi:
      fi:
    while  rpt  do
# If the point is singular when restricting the parametrization of the curve, we must find a different one
      a := a+1:
      l0 := simplify(subs(t=a, l)):
      den0 := simplify(subs(t=a,den)):
      l0_s := diff(l0,s):
	  rpt:=false:
      if resultant(l0,l0_s,s)=0 or resultant(l0,den0,s)=0 then
  	    rpt:=true:
# Now we get all the points of l in the critical line x=a
      elif degree(l0)=0 then
	    rpt:=true:
	  else
	    p0 := [solve(l0,s)][1]: #we just need one solution
	    if subs(t=a,s=p0,X_s)=[0,0,0] or subs(t=a,s=p0,X_t)=[0,0,0] then
	      rpt:=true:
		  fi:
        fi:
      od:
# the tangent vectors have first coordinate 1 and second coordinate l_s(a,p0)
    l_t0 := subs(t=a,diff(l,t)):
    v0 := [-subs(s=p0,l0_s),subs(s=p0,l_t0)]:
# We have a candidate point plus vector (so line)
    V := subs(t=a,s=p0,[X_t[1]*v0[1]+X_s[1]*v0[2],X_t[2]*v0[1]+X_s[2]*v0[2],X_t[3]*v0[1]+X_s[3]*v0[2]]):
    P := subs(t=a,s=p0,X):
# Now, P+wV is a line we now test if it is contained in the surface 
# We suppose now that the surface is smooth:
    if V<>[0,0,0] then
      pol := LinearAlgebra[CrossProduct](V,X-P):
      pol := gcd(numer(pol[1]),gcd(numer(pol[2]),numer(pol[3]))):
      if rem(pol, l, s)=0 then
        rectas_candidatas := [op(rectas_candidatas), [P[1]+V[1]*w, P[2]+V[2]*w, P[3]+V[3]*w]]:
        fi:
    elif subs(t=a,s=p0,X_t)<>[0,0,0] then
      V := subs(t=a,s=p0,X_t):
      pol := LinearAlgebra[CrossProduct](V,X-P):
      pol := gcd(numer(pol[1]),gcd(numer(pol[2]),numer(pol[3]))):
      if rem(pol, l, s)=0 then
        rectas_candidatas := [op(rectas_candidatas), [P[1]+V[1]*w, P[2]+V[2]*w, P[3]+V[3]*w]]:
        fi:
    elif subs(t=a,s=p0,X_s)<>[0,0,0] then
      V := subs(t=a,s=p0,X_s):
      pol := simplify([V[2]*(X[3]-P[3])-V[3]*(X[2]-P[2]),V[3]*(X[1]-P[1])-V[1]*(X[3]-P[3]),V[1]*(X[2]-P[2])-V[2]*(X[1]-P[1])]):
      pol := gcd(numer(pol[1]),gcd(numer(pol[2]),numer(pol[3]))):
      if rem(pol, l, s)=0 then
        rectas_candidatas := [op(rectas_candidatas), [P[1]+V[1]*w, P[2]+V[2]*w, P[3]+V[3]*w]]:
        fi:
    else
      print("Too singular point:", [a,p0]):
      fi:
    fi:
  od:
return simplify(rectas_candidatas):
end proc: