###############################################################################
###############################################################################
###  Maple procedures for the covariants appearing in: 
###  "Some covariants related to Steiner surfaces"
###  Franck Aries, Emmanuel Briand, Claude Bruchou 
###  01/31/2006 
### 
###  All of the programs are provided free for non-commercial, personal or non-profit use.
###  These programs are distributed in the hope that they will be useful, but WITHOUT ANY WARRANTY.
###############################################################################
###############################################################################

with(LinearAlgebra);

# In all what follows:
# x,lambda are lists of three unassigned names
# y,alpha are lists of four unassigned names
# f is a list of four ternary quadratic functions in the elements of x

## Auxiliary proocedures
# `covariant/aux/a`(f,x)
# `covariant/aux/a/matrix`(f,x)
# `covariant/aux/b`(f,x)
# `covariant/aux/b/matrix`(f,x)
# `covariant/aux/varXline/matrix`(lambda)
# `covariant/aux/ell/matrix`(f,x,y)
# `covariant/aux/matrixOfQuadraticForm`(F,x)
# `covariant/aux/HessianMatrix`(F,x,n)
# `covariant/aux/BorderedHessianMatrix`(F,x,vbord,hbord,n)
# `covariant/aux/HarmonicCovariantOfTwoConics`(C1,C2,x,lambda)

`covariant/Phi`:=table([
  1=`covariant/TangentPlane`,
  2=`covariant/PlaneSpanned`,
  3=`covariant/OtherLine`,
  4=`covariant/ImplicitEquation`,
  5=`covariant/AssociatedQuadric`,
  6=`covariant/PreimagePoint`,
  7=`covariant/CayleyCubic`,
  8=`covariant/TriplePoint`,
  9=`covariant/UnionOfTropes`,
 10=`covariant/Trihedron`,
 11=`covariant/ExceptionalTriangle`,
 12=`covariant/PlanePi`,
 13=`covariant/Conic`,
 14=`covariant/DualConic`,
 15=`covariant/Quadrilateral`
]):

## covariants
# `covariant/Discriminant`(f,x)
# `covariant/TangentPlane/matrix`(f,x,y)
# `covariant/TangentPlane`(f,x,y)
# `covariant/PlaneSpanned/matrix`(f,x,y,lambda)
# `covariant/PlaneSpanned`(f,x,y,lambda)
# `covariant/OtherLine/matrix`(f,x,lambda)
# `covariant/OtherLine`(f,x,lambda)
# `covariant/ImplicitEquation/matrix`(f,x,y)
# `covariant/ImplicitEquation`(f,x,y)
# `covariant/AssociatedQuadric`(f,x,y)
# `covariant/PreimagePoint/matrix`(f,x,y,lambda)
# `covariant/PreimagePoint`(f,x,y,lambda)
# `covariant/CayleyCubic/matrix`(f,x,alpha)
# `covariant/CayleyCubic`(f,x,alpha)
# `covariant/TriplePoint`(f,x,alpha)
# `covariant/UnionOfTropes`(f,x,y)
# `covariant/Trihedron`(f,x,y)
# `covariant/ExceptionalTriangle/matrix`(f,x)
# `covariant/ExceptionalTriangle`(f,x)
# `covariant/PlanePi`(f,x,y)
# `covariant/Conic`(f,x)
# `covariant/DualConic`(f,x,lambda)
# `covariant/Quadrilateral`(f,x)

###############################################################################
# Auxiliary Procedures
###############################################################################

`covariant/aux/a`:=proc(f,x)
  # The a-coeffs
  local a,i,j:
  a:=table([]):
  for i from 1 to 4 do
    for j from 1 to 3 do a[i,j]:=coeff(f[i],x[j],2) od:
  od:
  copy(a); 
end:

###

`covariant/aux/a/matrix`:=proc(f,x)
   # The matrix of the a-coeffs
   local a,i,j:
   a:=`covariant/aux/a`(f,x):
   Matrix(4,3,(i,j)->a[i,j]):
end:

###

`covariant/aux/b`:=proc(f,x)
  local b,L,i,j:
  # The b-coeffs
  b:=table([]):
  L:=[[2,3],[1,3],[1,2]]:
  for i from 1 to 4 do
    for j from 1 to 3 do 
      b[i,j]:=coeff(coeff(f[i],x[L[j][1]],1),x[L[j][2]],1)/2
    od:
  od:
  copy(b); 
end:

###

`covariant/aux/b/matrix`:=proc(f,x)
   # The matrix of the b-coeffs
   local b,i,j:
   b:=`covariant/aux/b`(f,x):
   Matrix(4,3,(i,j)->b[i,j]):
end:

###

`covariant/aux/varXline/matrix`:=proc(lambda)
   # The matrix of the x_i*<lambda,x>'s
   Matrix(3,6,[
    [lambda[1], 0       , 0        , 0        , lambda[3], lambda[2]],
    [0        ,lambda[2], 0        , lambda[3], 0        , lambda[1]],
    [0        , 0       , lambda[3], lambda[2], lambda[1], 0        ]
   ]):
end:

###

`covariant/aux/ell/matrix`:=proc(f,x,y)
   local T,C,i,j:
   # the matrix fo the ell's
   T:=`covariant/TangentPlane`(f,x,y):
   C:=[diff(T,x[1],x[1]),diff(T,x[2],x[2]),diff(T,x[3],x[3]),
       diff(T,x[2],x[3]),diff(T,x[1],x[3]),diff(T,x[1],x[2])]:
   # each entry of L corresponds to a column in the final matrix
   Matrix(3,6,(i,j)->1/6*diff(C[j],x[i])):
end:

###
 
`covariant/aux/matrixOfQuadraticForm`:=proc(F,x,n)
  local M,i,j:
  M:=Matrix(n,n,shape=symmetric):
  for i from 1 to n do
      M[i,i]:=coeff(F,x[i],2):
      for j from 1 to i-1 do
          M[i,j]:=coeff(coeff(F,x[j],1),x[i],1)/2
      od: 
  od:
  M:
end:

###

`covariant/aux/HessianMatrix`:=proc(F,x,n)
   # beware the factor 1/2.
   local i,j:
   Matrix(n,n,shape=symmetric,(i,j)->diff(F,x[i],x[j])/2):
end;

###

`covariant/aux/BorderedHessianMatrix`:=proc(F,x,vbord,hbord,n)
    local M1,M2,M3,M4,i,j:
    M1:=`covariant/aux/HessianMatrix`(F,x,n):
    M2:=Matrix(n,1,(i,j)->vbord[i]):
    M3:=Matrix(1,n,(i,j)->hbord[j]):
    M4:=Matrix(1,1,0):
    Matrix([[M1,M2],[M3,M4]]):
end:

###

`covariant/aux/HarmonicCovariantOfTwoConics`:=proc(C1,C2,x,lambda)
   # the covariant of the the conics C1(x) and C2(x) 
   # R(C1,C2,x,lambda)=0 (identically wrt x)
   # is the condition for the line <x,lambda> to be harmonically intercepted 
   # by the the two conics  
   # cf (Salmon, conics parag. 377)
  local u,v,H,h,R;
  H:=`covariant/aux/BorderedHessianMatrix`(u*C1+v*C2,x,lambda,lambda,3):
  h:=Determinant(H):
  R:=coeff(coeff(h,u,1),v,1):
  collect(R,lambda,distributed):
end:
   
###############################################################################
# The discriminant (Delta).
###############################################################################
 
`covariant/Discriminant`:=proc(f,x)
  local Q,P,R,y,alpha,i:
  Q:=`covariant/AssociatedQuadric`(f,x,y):
  P:=`covariant/TriplePoint`(f,x,alpha):
  R:=subs(seq(y[i]=coeff(P,alpha[i],1),i=1..4),Q):
end:

###############################################################################
# The Tangent Plane (Phi_1).
###############################################################################

`covariant/TangentPlane/matrix`:=proc(f,x,y)
   local M1,M2,i,j:
   # matrix in formula (12)
   M1:=Matrix(4,3,(i,j)->diff(f[i],x[j])):
   M2:=Matrix(4,1,(i,j)->y[i]):
   Matrix([[M1,M2]]):
end;

###

`covariant/TangentPlane`:=proc(f,x,y)
   local R:
   R:=1/8*Determinant(`covariant/TangentPlane/matrix`(f,x,y));
   collect(R,y,distributed):
end; # 

###############################################################################
# The plane spanned by the image of a line (Phi_2)
###############################################################################

`covariant/PlaneSpanned/matrix`:=proc(f,x,y,lambda) 
  # the matrix in formula (14) 
  local Ma,Mb,Mup,Mlambda,My,M0,i,j:

  Ma:=`covariant/aux/a/matrix`(f,x):
  Mb:=`covariant/aux/b/matrix`(f,x):

  Mup:=Matrix([[Ma,2*Mb]]):
  Mlambda:=`covariant/aux/varXline/matrix`(lambda):
  My:=Matrix(4,1,(i,j)->y[i]):
  M0:=Matrix(3,1,0):

  Matrix([[Mup,My],[Mlambda,M0]]);
end:

###

`covariant/PlaneSpanned`:=proc(f,x,y,lambda)
   local R:
   R:=1/2*Determinant(`covariant/PlaneSpanned/matrix`(f,x,y,lambda));
   collect(R,y,distributed):
end;

###############################################################################
# The other line whose image spans the same plane (Phi_3)
###############################################################################

`covariant/OtherLine/matrix`:=proc(f,x,lambda)
  # the matrix in Formula  (16)
  local Ma,Mb,Mup,Mlambda,M0,Mx,i,j:

  Ma:=`covariant/aux/a/matrix`(f,x):
  Mb:=`covariant/aux/b/matrix`(f,x):

  Mup:=Matrix([[Ma,2*Mb]]):
  Mlambda:=`covariant/aux/varXline/matrix`(lambda):
  M0:=Matrix(4,1,0):
  Mx:=Matrix(3,1,(i,j)->x[i]):

  Matrix([[Mup,M0],[Mlambda,Mx]]);
end:

###

`covariant/OtherLine`:=proc(f,x,lambda)
   local R:
   R:=Determinant(`covariant/OtherLine/matrix`(f,x,lambda));
   collect(R,x,distributed):
end:

###############################################################################
# The implicit equation (Phi_4)
###############################################################################

`covariant/ImplicitEquation/matrix`:=proc(f,x,y)
  # The matrix in Formula (18)
  local Ma,Mb,Mup,Mell,My,M0,i,j:

  Ma:=`covariant/aux/a/matrix`(f,x):
  Mb:=`covariant/aux/b/matrix`(f,x):

  Mup:=Matrix([[Ma,Mb]]):

  Mell:=`covariant/aux/ell/matrix`(f,x,y):
  My:=Matrix(4,1,(i,j)->y[i]):
  M0:=Matrix(3,1,0):

  Matrix([[Mup,My],[Mell,M0]]);
end:

###

`covariant/ImplicitEquation`:=proc(f,x,y)
   local R:
   R:=6^3*Determinant(`covariant/ImplicitEquation/matrix`(f,x,y)):
   collect(R,y,distributed):
end:

###############################################################################
# The Associated Quadric (Phi_5)
###############################################################################


`covariant/AssociatedQuadric`:=proc(f,x,y)
  local T,P,S,lambda,i,j,k,R; 
  T:=`covariant/TangentPlane`(f,x,y); 
  P:=`covariant/PlaneSpanned`(f,x,y,lambda); 
  R:=0;
  for i from 1 to 3 do
   for j from 1 to 3 do
     for k from 1 to 3 do
       R:=R+diff(T,x[i],x[j],x[k])*diff(P,lambda[i],lambda[j],lambda[k]);
     od;
   od;
  od;
  collect(R/6,lambda,distributed);
end;

###############################################################################
# Preimage of a point (Phi_6)
###############################################################################

`covariant/PreimagePoint/matrix`:=proc(f,x,y,lambda)
  # The matrix in Formula (20)
  local Ma,Mb,Mup,Mell,Mlambda,My,M0,i,j:

  Ma:=`covariant/aux/a/matrix`(f,x):
  Mb:=`covariant/aux/b/matrix`(f,x):

  Mup:=Matrix([[Ma,Mb]]):

  Mell:=`covariant/aux/ell/matrix`(f,x,y):
  M0:=Matrix(4,1,0):
  Mlambda:=Matrix(3,1,(i,j)->lambda[i]):

  Matrix([[Mup,M0],[Mell,Mlambda]]);
end:

###

`covariant/PreimagePoint`:=proc(f,x,y,lambda)
   local R:
   R:=6^3*Determinant(`covariant/PreimagePoint/matrix`(f,x,y,lambda)):
   collect(R,lambda,distributed):
end:

###############################################################################
# The Cayley Cubic Surface (Phi_7)
###############################################################################

`covariant/CayleyCubic/matrix`:=proc(f,x,alpha)
  # The matrix in Formula (21)
  local F,i:
  F:=expand(add(alpha[i]*f[i],i=1..4)):
  `covariant/aux/matrixOfQuadraticForm`(F,x,3):
end:

###

`covariant/CayleyCubic`:=proc(f,x,alpha)
   local R:
   R:=Determinant(`covariant/CayleyCubic/matrix`(f,x,alpha)):
   collect(R,alpha,distributed)
end:

###############################################################################
# The Triple Point (Phi_8)
###############################################################################

`covariant/TriplePoint`:=proc(f,x,alpha)
local Q,C,y,R,i,j:
Q:=`covariant/AssociatedQuadric`(f,x,y):
C:=`covariant/CayleyCubic`(f,x,alpha):
R:=0:
  for i from 1 to 4 do
     for j from 1 to 4 do
       R:=R+diff(Q,y[i],y[j])*diff(C,alpha[i],alpha[j]);
     od;
  od;
collect(R,alpha,distributed):
end;

###############################################################################
# Union of the tropes (Phi_9) 
###############################################################################

`covariant/UnionOfTropes`:=proc(f,x,y)
  local E,Q,R:
  E:=`covariant/ImplicitEquation`(f,x,y);
  Q:=`covariant/AssociatedQuadric`(f,x,y);
  R:=E+Q^2:
  collect(R,y,distributed);
end:

###############################################################################
# The Trihedron of the singular lines (Phi_10)
###############################################################################

`covariant/Trihedron`:=proc(f,x)
   local y,E,P,R:
   E:=`covariant/ImplicitEquation`(f,x,y):
   P:=`covariant/TriplePoint`(f,x,alpha):
   R:=add(diff(E,y[i])*diff(P,alpha[i]),i=1..4):
   collect(R,x,distributed):
end:

###############################################################################
# The Exceptional Triangle (Phi_11)
###############################################################################

`covariant/ExceptionalTriangle/matrix`:=proc(f,x)
  local lambda,L:
  L:=`covariant/OtherLine`(f,x,lambda):
  `covariant/aux/matrixOfQuadraticForm`(L,lambda,3)
end:

###

`covariant/ExceptionalTriangle`:=proc(f,x)
  local M,R:
  M:=`covariant/ExceptionalTriangle/matrix`(f,x):
  R:=Determinant(M):
  collect(R,x,distributed)
end:

###############################################################################
# Polar plane Pi of the associated Quadric and the triple point (Phi_12) 
###############################################################################

`covariant/PlanePi`:=proc(f,x,y)
  local alpha,Q,P,R,i:
  Q:=`covariant/AssociatedQuadric`(f,x,y):
  P:=`covariant/TriplePoint`(f,x,alpha):
  R:=0;
  for i from 1 to 4 do 
    R:=R+coeff(P,alpha[i],1)*diff(Q,y[i]);
  od;
  collect(R,y,distributed);
end:

###############################################################################
# The conic, inverse image of Pi (Phi_13) 
###############################################################################

`covariant/Conic`:=proc(f,x)
   local y,Pi,i,R:
   Pi:=`covariant/PlanePi`(f,x,y):
   R:=subs(seq(y[i]=f[i],i=1..4),Pi):
   collect(R,x,distributed)
end:

###############################################################################
# The Dual Conic to the inverse image of Pi (Phi_14) 
###############################################################################

`covariant/DualConic`:=proc(f,x,lambda)
  local y,Q,H,R,i,j:
  Q:=`covariant/AssociatedQuadric`(f,x,y):
  H:=`covariant/aux/HarmonicCovariantOfTwoConics`:
  R:=0:
  for i from 1 to 4 do  
    for j from 1 to 4 do
       R:=R+diff(Q,y[i],y[j])*H(f[i],f[j],x,lambda):
    od;
  od;
  collect(R,lambda,distributed);
end;

###############################################################################
# The quadrilateral (Phi_15)
###############################################################################

`covariant/Quadrilateral`:=proc(f,x)
   local Q,R,i,y:
   Q:=`covariant/AssociatedQuadric`(f,x,y):
   R:=subs(seq(y[i]=f[i],i=1..4),Q):
   collect(R,x,distributed):
end:

###############################################################################

`covariant/ListPhi`:=proc(f,x,y,lambda,alpha)
# To get the covariants Phi_1...Phi_{15} as a list.
[
`covariant/TangentPlane`(f,x,y),
`covariant/PlaneSpanned`(f,x,y,lambda),
`covariant/OtherLine`(f,x,lambda),
`covariant/ImplicitEquation`(f,x,y),
`covariant/AssociatedQuadric`(f,x,y),
`covariant/PreimagePoint`(f,x,y,lambda),
`covariant/CayleyCubic`(f,x,alpha),
`covariant/TriplePoint`(f,x,alpha),
`covariant/UnionOfTropes`(f,x,y),
`covariant/Trihedron`(f,x),
`covariant/ExceptionalTriangle`(f,x),
`covariant/PlanePi`(f,x,y),
`covariant/Conic`(f,x),
`covariant/DualConic`(f,x,lambda),
`covariant/Quadrilateral`(f,x)
]
end:

###############################################################################
# END
###############################################################################