#
#
# Procedures for the search of triangles with many common transversals
#
#################################################################################
#
#  Data representation
#
#   Lines here are always given by two points. 
#   (Thus easily transformed into parametrizations).
#   As much as it is possible, everything is parametrized.
################################################################################
#
#  Pluecker(Line)            Returns the pluecker vector of a given line
#  Meet(Line1,Line2)         Evaluates the equation for the lines to meet
#  ParameterizeLine(Line,t)  Parameterizes the line with t=0 the first point
#                             and t=1 the second.
#  MakeRegulus(Line1, Line2, Line3, t)   Gives the linked parametrizations of 
#                            Lines 1 and 2 for the regulus defined by the 3 lines
#                            This parametrizes the hyperboloid, returning
#                            `false' if there is some problem with the
#                            parametrization 
#  Wiggle(Triangle)          This wiggles the triangle
#  numberOfTransversals(Triangle,Num)
#                            This coomputes the numnber of transversals.  It tries
#                            to handle all errors; when one is encountered, it
#                            Wiggles(Triangle), increments Num and calls itself.
#                            The third call results in an automatic end of the
#                            procedure. 
################################################################################
Pluecker := proc(Line)
#
#  Computes the Pluecker vector of a line
#  Returns [p01, p02, p03, p12, p13, p23]
#
return([Line[2][1]-Line[1][1], Line[2][2]-Line[1][2], Line[2][3]-Line[1][3] 
       ,Line[1][1]*Line[2][2]-Line[1][2]*Line[2][1]
       ,Line[1][1]*Line[2][3]-Line[1][3]*Line[2][1]
       ,Line[1][2]*Line[2][3]-Line[1][3]*Line[2][2]]):
end proc:
################################################################################
Meet := proc(Line1,Line2)
#
#  Evaluates the equation for Line1 to meet Line2 in Pluecker coordinates
#
 local L1, L2:
 L1:=Pluecker(Line1):
 L2:=Pluecker(Line2):
 return(L1[1]*L2[6]-L1[2]*L2[5]+L1[3]*L2[4]
       +L1[4]*L2[3]-L1[5]*L2[2]+L1[6]*L2[1]):
end proc:
################################################################################
ParameterizeLine := proc(Line,t)
#
#  Given a line by two points [ [x1,y1,z1], [x2,y2,z2] ]
#  Parameterizes the line with t=0 the first point and t=1 the second
#  Returns [x(t),y(t),z(t)]
#
 return([Line[1][1]+t*(Line[2][1]-Line[1][1])
        ,Line[1][2]+t*(Line[2][2]-Line[1][2])
        ,Line[1][3]+t*(Line[2][3]-Line[1][3])]):
end proc:
################################################################################
MakeRegulus := proc(Line1,Line2, Line3,t)
#
#  This gives parametrized points on Line1 and Line2 that are collinear with
#  Line 3.  This defines the regulus of lines meeting all three.  It also 
#  serves to parametrize the hyperboloid, via the command
#    ParameterizeLine(MakeRegulus(Line1, Line2, Line3,t), s):
#  Returns the boolean variable `false' if for some reason it cannot
#  make the Regulus.
#
 local Equation:
 Equation:=Meet([ParameterizeLine(Line1,t),ParameterizeLine(Line2,s)],Line3);
 if indets(Equation) = {s,t} then 
   return([ParameterizeLine(Line1,t),
    subs(s=solve(Equation=0,s),ParameterizeLine(Line2,s))]):
  else
   return(evalb(0=1)):
 end if:
end proc:
################################################################################
Wiggle:=proc(Old)
#
#  This wiggles a given configuration, multiplying all 
# coordinates by 1000 and adding a random element of {-20..20}
#
local New, i, j, k, die:
global _seed:
 _seed:=1:
 New:=Old:
 die := rand(-20..20):
 for i from 1 to 4 do 
  for j from 1 to 3 do
   for k from 1 to 3 do 
    New[i][j][k]:=1000*Old[i][j][k]+die():
   od:
  od:
 od:
 return(New):
end proc:
##############################################################################
numberOfTransversals := proc(T,Num)
#
#  Procedure to compute the number of transversals of a given
#  triangle T.  If it detects a four-tuple of lines whose common
#  transversals are either incorrect in number or dimension, then
#  it wiggles the triangle and calls itself to recompute the number.
#  This is troubling, but I'll put in a trip switch: it can only 
#  do this thrice.
#
#   Num keeps track of how often it has been called
#
#################################################################
 local file, NumberOfTransversals, Line, L, Regulus, 
       Equation, Points, j, P, M, N, SOL, L1, L2, EE:
 ################################################################
 NumberOfTransversals:=0:

 if Num>1 then 
  file:=fopen("CheckMe", APPEND):
  fprintf(file, "Data:=%a:\n\n",T):
  fclose(file):
 if Num >4 then quit: end if:
 end if:

 ################################################################
#
#  Checks to make sure that the lines in different triangles do not meet
#
 for L1 in combinat[choose](T[1],2) do 
  EE :=  [combinat[choose](T[2],2)[], 
            combinat[choose](T[3],2)[],
            combinat[choose](T[4],2)[]]:
  for L2 in EE do
   if Meet(L1,L2)=0 then 
    return(numberOfTransversals(Wiggle(T),Num+1)):
   end if:
  end do:
 end do:
 for L1 in combinat[choose](T[2],2) do 
  EE :=  [combinat[choose](T[3],2)[],
           combinat[choose](T[4],2)[]]:
  for L2 in EE do
   if Meet(L1,L2)=0 then 
    return(numberOfTransversals(Wiggle(T),Num+1)):
   end if:
  end do:
 end do:
 for L1 in combinat[choose](T[3],2) do 
  for L2 in combinat[choose](T[4],2) do
   if Meet(L1,L2)=0 then 
    return(numberOfTransversals(Wiggle(T),Num+1)):
   end if:
  end do:
 end do:

 ################################################################
 NumberOfTransversals:=0:
 ################################################################
 for Line[1] in combinat[choose](T[1],2) do 
  for  Line[2] in combinat[choose](T[2],2) do 
   L[2]:=ParameterizeLine(Line[2],t):
   for  Line[3] in combinat[choose](T[3],2) do 
    L[3]:=ParameterizeLine(Line[3],t):
    for  Line[4] in combinat[choose](T[4],2) do 
     L[4]:=ParameterizeLine(Line[4],t):

     Regulus:=MakeRegulus(Line[1], Line[2], Line[3], t):
     if type(Regulus,boolean) then 
     #
     #  This is the error condition for the parametrization of the Regulus
     #
     return(numberOfTransversals(Wiggle(T),Num+1)):
     #
     end if:
     Equation:=numer(Meet(Line[4],Regulus)):
     #
     if indets(Equation)<>{t} then 
     #
     return(numberOfTransversals(Wiggle(T),Num+1)):
     #
     end if:
     Points:=[solve(Equation=0,t)]:
     #
     if nops(Points)<>2 or indets(Points)<>{} then 
     #
     return(numberOfTransversals(Wiggle(T),Num+1)):
     #
     end if:
     #
     if not type(Points[1],nonreal) then 
      for j from 1 to 2 do 
       P:=Points[j]:
       if 0<=evalf(P) and evalf(P)<=1 then 
        try
          M[j]:=subs(t=P,Regulus)
         catch:
          return(numberOfTransversals(Wiggle(T),Num+1)):
        end try: 
#       
# By the parametrization of the regulus, the second coordinate is on Line[2]
#       
        try
          P:=evalf(solve(L[2][1]=M[j][2][1],t))
         catch:
          return(numberOfTransversals(Wiggle(T),Num+1)):
        end try: 

        if indets(P)<>{} then 
          return(numberOfTransversals(Wiggle(T),Num+1)):
        end if:

        if 0<=P and P<=1 then
#        
# The lines 3 and 4 are slightly different
#       

         N[j]:=ParameterizeLine(M[j],s):

         try
           SOL:=solve({L[3][1]=N[j][1],L[3][2]=N[j][2],L[3][3]=N[j][3]},{s,t})
          catch:
           return(numberOfTransversals(Wiggle(T),Num+1)):
         end try: 

         try
           P:=evalf(subs(op(SOL),t))
          catch:
           return(numberOfTransversals(Wiggle(T),Num+1)):
         end try: 
 
         if indets(P)<>{} then 
           return(numberOfTransversals(Wiggle(T),Num+1)):
         end if:

         if 0<=P and P<=1 then

          try
            SOL:=solve({L[4][1]=N[j][1],L[4][2]=N[j][2],L[4][3]=N[j][3]},{s,t})
           catch:
            return(numberOfTransversals(Wiggle(T),Num+1)):
          end try: 

          try
            P:=evalf(subs(op(SOL),t))
           catch:
            return(numberOfTransversals(Wiggle(T),Num+1)):
          end try: 
   
          if indets(P)<>{} then 
            return(numberOfTransversals(Wiggle(T),Num+1)):
          end if:

          if 0<=P and P<=1 then
           NumberOfTransversals:=NumberOfTransversals+1:     
          end if:
         end if:
        end if:
       end if:
      end do:
     end if:
##################################################################
    end do:
   end do:
  end do:
 end do:
 return(NumberOfTransversals):
end proc: