######################################################################
##RichardSerge.txt: Save this file as  RichardSerge.txt              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read RichardSerge.txt<Enter>                       #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: May, 2014
 
print(`Created:  May 2014`):
print(` This is RichardSerge.txt `):
print(`It is one of the packages that accompany the article `):
print(`  Enumerative Geometrical Genealogy (Or: The Sex Life of Points and Lines)  `):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`--------------------------------------------------------------------------`):
 print(`For a list of the SUPPORTINGprocedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------------------`):

print(`--------------------------------------------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------------------`):



with(combinat):
with(linalg):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` ApplyPT, ApplyPT1, Du11, Du1, FprT, IsPE1, Le, Milim1, Pt, Rep1, Rep2 `):

else
ezra(args):
fi:

end:

ezra:=proc(): 

if nargs=0 then
print(`Contains procedures:`):
print(` FindSTW1, FindSTW2, FindSTW3, GCP, GIP, IsC, IsI,  IsPE, MatI, MatI1, Nisim1, Nisim2, RCP, RIP, STwords, Tk, Tw  `):

elif nargs=1 and  args[1]=ApplyPT1 then
print(`ApplyPT1(A,Pt): applies a projective transformation, given by a 3 by 3 matrix A`):
print(`to the point Pt, try:`):
print(`ApplyPT1([[1,2,3],[2,3,4],[3,4,5]],[3,7]);`):


elif nargs=1 and  args[1]=ApplyPT then
print(`ApplyPT(A,L): applies a projective transformation, given by a 3 by 3 matrix A`):
print(`to the list of points L, try:`):
print(`ApplyPT([[1,2,3],[2,3,4],[3,4,5]],[[3,7],[6,3]]);`):

elif nargs=1 and  args[1]=Du1 then
print(`Du1(P): the dual of the polygon P, given in terms of lines`):
print(`Try: `):
print(` Du1([x+y+1,2*x+4*y+1]); `):

elif nargs=1 and  args[1]=Du11 then
print(`Du11(L): the dual of the line L`):
print(`Try: `):
print(` Du11(x+y+1); `):

elif nargs=1 and  args[1]=FindSTW1 then
print(`FindSTW1(n,k,R): Finds all Schwartz-Tabachnikov words of the First kind (preserving being inscribed) for an  n-gon`):
print(`by repeating a word of even length k up to R times. Try:`):
print(`FindSTW1(6,2,6);`):

elif nargs=1 and  args[1]=FindSTW2 then
print(`FindSTW2(n,k,R): Finds all Schwartz-Tabachnikov words of the Second kind (leading to projective equivalence`):
print(`by repeating a word of even length k up to R times. Try:`):
print(`FindSTW2(6,2,6);`):

elif nargs=1 and  args[1]=FindSTW3 then
print(`FindSTW3(n,R): Finds all Schwartz-Tabachnikov words of the First Kind(preserving being inscribed)`):
print(`by repeating a word of the form [1,r] with r between 2 and n/2. Try:`):
print(`FindSTW3(6,4);`):


elif nargs=1 and  args[1]=FprT then
print(`FprT(L1,L2): are the list of points L1, L2, projectively equivalent?, finds a matrix`):
print(`A such that L1 goes to L2. Try:`):
print(`FprT([[1,1],[2,2],[3,3],[4,4],[5,5],[6,6]],[[1,1],[2,2],[3,3],[4,4],[5,5],[6,6]]):`):

elif nargs=1 and  args[1]=IsC then
print(` IsC(L): inputs a list of lines and decides whether they are circumscribed `):
print(`Try:`):
print(`IsC(Du1(RIP(7,100)));`):

elif nargs=1 and  args[1]=IsI then
print(` IsI(L): inputs a list of points and decides whether they are inscribed `):
print(`Try:`):
print(`IsI(RIP(7,100));`):


elif nargs=1 and  args[1]=IsPE then
print(`IsPE(P,L): is some circular shift of  list of lines L projectively equivalent to P?`):
print(`if it is, it returns the shift, otherwise FAIL. Try:`):
print(`P:=RIP(6,40); IsPE(P,Tk(P,2)); `):


elif nargs=1 and  args[1]=IsPE1 then
print(`IsPE1(P,L): is the list of lines L projectively equivalent to P`):

elif nargs=1 and  args[1]=GCP then
print(`GCP(t,k): a generic CIRCUMSCRIBED k-gon parametrized by t[1],t[2],.., t[k] inscribed in the parabola y=x^2`):
print(`GCP(t,8);`):

elif nargs=1 and  args[1]=GIP then
print(`GIP(t,k): a generic INSCRIBED k-gon parametrized by t[1],t[2],.., t[k] inscribed in the parabola y=x^2`):
print(`GIP(t,8);`):

elif nargs=1 and  args[1]=Le then
print(`Le(P1,P2), the line joining P1 and P2`):

elif nargs=1 and  args[1]=MatI1 then
print(`MatI1(L): inputs a hexagon (given as a list of 6 points) and outputs the matrix that would decide whether they are inscribed in a conic`):
print(`Try:`):
print(`MatI1(GIP(t,6));`):

elif nargs=1 and  args[1]=MatI then
print(`MatI(L): inputs a polygon (given as a list of points, at least 6) and outputs the matrices that would decide whether they are inscribed in a conic`):
print(`Try:`):
print(`MatI(GIP(t,7));`):


elif nargs=1 and  args[1]=Milim1 then
print(`Milim1(n,r): all the words without double-letters  of length n in {1,...,r}.`):
print(`Try:`):
print(`Milim1(6,4);`):

elif nargs=1 and  args[1]=Nisim1 then
print(`Nisim1(N,k,R): all the Schwartz-Tabachnikov miracles of the first kind for`):
print(`inscribed n-gons from n=6 to n=N, by repeating a word of length k up to R times. Try:`):
print(`Nisim1(12,2,6);`):

elif nargs=1 and  args[1]=Nisim2 then
print(`Nisim2(N,k,R): all the Schwartz-Tabachnikov miracles of the second kind for`):
print(`inscribed n-gons from n=6 to n=N, by repeating a word of length k up to R times. Try:`):
print(`Nisim2(12,2,6);`):

elif nargs=1 and  args[1]=Pt then
print(`Pt(Le1,Le2), the point of intersection of Le1 and Le2`):

elif nargs=1 and  args[1]=Rep1 then
print(`Rep1(n,w,R): inputs a pos. integers n and R, and a word in {1,2,..., n/2}, of even length `):
print(` finds whether `):
print(` repeating it r times <=R, has property that Tw preserves the property of being `):
print(` an inscibed polygon. If successful, it returns (the smallest) r, otherwise it returns FAIL `):
print(` Try: `):
print(`Rep1(6,[2,1],5);`):


elif nargs=1 and  args[1]=Rep2 then
print(`Rep2(n,w,R): inputs a pos. integers n and R, and a word in {1,2,..., n/2}, of even length `):
print(` finds whether `):
print(` repeating it r times <=R, and dropping the last letter`):
print(` has property that Tw maps an inscribed n-gon to something projectively equivalent `):
print(` If successful, it returns (the smallest) r, otherwise it returns FAIL `):
print(` Try: `):
print(`Rep2(6,[2,1],5);`):

elif nargs=1 and  args[1]=RCP then
print(`RCP(k,K): a random CIRCUMSCRIBED k-gon of the parabola y=x^2 with parameters drawn randomly from 1 to K`):
print(`RCP(8,100);`):

elif nargs=1 and  args[1]=RIP then
print(`RIP(k,K): a random k-gon  inscribed in the parabola y=x^2, using values from 1 to K`):
print(`Try:`):
print(`RIP(8,100);`):


elif nargs=1 and  args[1]=STwords then
print(`STwords(n,r,s): all the good Schwartz-Tabachnikov words of even length r with largest letter s that`):
print(`map an inscribed n-gon to another one. `):
print(` Try: `):
print(` STwords(6,2,2); `):

elif nargs=1 and  args[1]=Tk then
print(`Tk(L,k): the map T_k(P), Try:`):
print(`Tk([[0,0],[1,0],[2,1],[4,2],[3,1],[4,0]],2);`):

elif nargs=1 and  args[1]=Tw then
print(`Applies the word w to the polygon L, try:`):
print(`P:=RIP(8,20): Tw(P,[2,1,2,1,2]);`):


else
 print(`There is no ezra for`,  args ):
fi:

end:





#The point of intersection of lines Le1 and Le2
Pt:=proc(Le1,Le2) local q:q:=solve(
{numer(normal(Le1)),numer(normal(Le2))},{x,y}):
[normal(simplify(subs(q,x))),normal(simplify(subs(q,y)))]:end:


#Def(Area of triangle ABC)
AREA:=proc(A,B,C):normal(expand((B[1]*C[2]-B[2]*C[1]-A[1]*C[2]+A[2]*C[1]
-B[1]*A[2]+B[2]*A[1])/2)):end:

#The eq. of the line joining A and B
Le:=proc(A,B) AREA(A,B,[x,y]):end:

#T1k(L,k): the map T_k(P), Try:
#T1k([[0,0],[1,0],[2,1],[4,2],[3,1],[4,0]]);
T1k:=proc(L,k) local i,M,n:

n:=nops(L):
M:=[]:

for i from 1 to n-k do
 
 M:=[op(M),Le(L[i],L[i+k])]:
od:

for i from n-k+1 to  n do
 M:=[op(M),Le(L[i],L[i+k-n])]:
od:
M:
end:

#T2k(L,k): the map T_k(P), Try:
#T2k([[0,0],[1,0],[2,1],[4,2],[3,1],[4,0]]);
T2k:=proc(L,k) local i,M,n:

n:=nops(L):
M:=[]:

for i from 1 to n-k do
 M:=[op(M),Pt(L[i],L[i+k])]:
od:

for i from n-k+1 to  n do
 M:=[op(M),Pt(L[i],L[i+k-n])]:
od:
M:
end:


#Tk(L,k): the map T_k(P), Try:
#Tk([[0,0],[1,0],[2,1],[4,2],[3,1],[4,0]]);
Tk:=proc(L,k) local i,M,n,kama:

if nops(convert(L,set))<>nops(L) then
 RETURN(FAIL):
fi:

if type(L[1],list) then
 T1k(L,k):
else
 T2k(L,k):
fi:

end:



#IsI1(L): inputs 6 points and decides whether they are inscribed
IsI1:=proc(L) local i:

if nops(L)<>6 then
 RETURN(FAIL):
fi:

evalb(normal(det([seq([L[i][1]^2,L[i][1]*L[i][2],L[i][2]^2,L[i][1],L[i][2],1],i=1..6)]))=0):

end:

#IsI(L): inputs a list of points and decides whether they are inscribed
IsI:=proc(L) local i:

if not( type(L,list) and {seq(type(L[i],list),i=1..nops(L))}={true} and 
 {seq(nops(L[i]),i=1..nops(L))}={2}) then
print(`Bad input`):
  RETURN(FAIL):
fi:

if nops(L)<6 then
 true:
else
 evalb({seq(IsI1([op(1..5,L), L[i] ]),i=6..nops(L))}={true}):
fi:



end:





#GIP(t,k): a generic inscribed k-gon parametrized by t[1],t[2],.., t[k] inscribed in the parabola y=x^2
GIP:=proc(t,k) local i:
[seq( [t[i],t[i]^2] ,i=1..k)]:
end:

#GCP(t,k): a generic circumscibed k-gon parametrized by t[1],t[2],.., t[k] inscribed in the parabola y=x^2
GCP:=proc(t,k) local i:
[seq(2*t[i]*x-y-t[i]^2 ,i=1..k)]:
end:

#RIP(k,K): a random inscribed k-gon  drawn with coeff. between 1 and K
#RIP(6,10);
RIP:=proc(k,K) local i,ra,gu,a:
ra:=rand(1..K):
gu:=[]:
for i from 1 to k do
a:=ra()+I*ra():
gu:=[op(gu),[a,a^2]]:
od:
gu:

end:

#RCP(k,K): a random circumscribed k-gon  drawn with coeff. between 1 and K
#RCP(6,10);
RCP:=proc(k,K) local i,ra,gu,a:
ra:=rand(1..K):
gu:=[]:
for i from 1 to k do
a:=ra()+I*ra():
gu:=[op(gu), evalc(2*a*x-y-a^2)]:
od:
gu:

end:


#Applies the word w to the polygon L, try:
#P:=RIP(8,100): Tw(P,[2,1,2,1,2]);
Tw:=proc(P,w) local P1,i:
P1:=P:

for i from 1 to nops(w) do
P1:=Tk(P1,w[i]):

if P1=FAIL then
 RETURN(FAIL):
fi:

od:

P1:

end:


#Milim1(n,r): all the words without double-letters  of length n in {1,...,r}.
#Try:
#Milim1(6,4);
Milim1:=proc(n,r) local gu,i,mu,j,w:
if n=0 then
RETURN({[]}):
fi:

if n=1 then
 RETURN({seq([i],i=1..r)}):
fi:

gu:={}:
mu:=Milim1(n-1,r):

for w in mu do
 for j from 1 to r do
  
if w[n-1]<>j then
  gu:=gu union {[op(w),j]}:
fi:
 od:
od:

gu:
end:

#STwords(n,r,s): all the good Schwartz-Tabachnikov words of even length r with largest letter s that
#map an inscribed n-gon to another one.
#Try:
#STwords(6,2,2);
STwords:=proc(n,r,s) local mu,gu,w,P:

if r mod 2<>0 then
 print(`r must be even`):
 RETURN(FAIL):
fi:



mu:=Milim1(r,s):
gu:={}:

for w in mu do

 P:=RIP(n,100):

 if IsI(Tw(P,w)) then

  P:=RIP(n,1000):

    if IsI(Tw(P,w)) then
     gu:=gu union {w}:
   fi:
 fi:

od:

gu:

end:

#Du11(L): the dual of the line L
#Try:
#Du11(x+y+1):
Du11:=proc(L) local c:


c:=subs({x=0,y=0},L):

if c=0 then 
 RETURN(FAIL):
fi:

[coeff(L,x,1)/c,coeff(L,y,1)/c]:

end:


#Du1(P): the dual of the polygon P, given in terms of lines
Du1:=proc(P) local i,gu,mu:
gu:=[]:

if (coeff(P[1],x,1)=0 or coeff(P[1],y,1)=0) then
 RETURN(FAIL):
fi:

for i from 1 to nops(P) do

mu:=Du11(P[i]):

if mu=FAIL then
  RETURN(FAIL):
else
gu:=[op(gu),mu]:
fi:

od:

gu:

end:


#IsC(P): inputs a list of lines and decides whether they are circumscribed
IsC:=proc(P) local P1:
P1:=Du1(P):
if P1=FAIL then
 RETURN(FAIL):
else

IsI(P1):

fi:

end:


#ApplyPT1(A,P1): applies a projective transformation, given by a 3 by 3 matrix A
#to the point Pt, try:
#ApplyPT1([[1,2,3],[2,3,4],[3,4,5]],[3,7]);
ApplyPT1:=proc(A,P1) local gu,i:
gu:=normal([seq(A[i][1]*P1[1]+A[i][2]*P1[2]+A[i][3],i=1..3)]):

if gu[3]=0 then
  RETURN(FAIL):
else
normal([gu[1]/gu[3],gu[2]/gu[3]]):
fi:

end:


#ApplyPT(A,L): applies a projective transformation, given by a 3 by 3 matrix A
#to the list of points L, try:
#ApplyPT1([[1,2,3],[2,3,4],[3,4,5]],[3,7]);
ApplyPT:=proc(A,L) local i;
[seq(ApplyPT1(A,L[i]),i=1..nops(L))];
end:




#FprT(L1,L2): are the list of points L1, L2, projectively equivalent?, finds a matrix
#A such that L1 goes to L2. Try:
#FprT([[1,1],[2,2],[3,3],[4,4],[5,5],[6,6]],[[1,1],[2,2],[3,3],[4,4],[5,5],[6,6]]):
FprT:=proc(L1,L2) local A,eq,var,L1a,a,i,j,var1,t:
if nops(L1)<>nops(L2) then
 RETURN(FAIL):
fi:

A:=[seq([seq(a[i,j],j=1..3)],i=1..3)]:
var:={seq(seq(a[i,j],j=1..3),i=1..3)}:
L1a:=ApplyPT(A,L1):

eq:={seq(numer(L1a[i][1]-L2[i][1]),i=1..nops(L2)),seq(numer(L1a[i][2]-L2[i][2]),i=1..nops(L2))}:

var:=solve(eq,var):
A:= subs(var,A):

if A=[[0$3]$3] then
 RETURN(FAIL):
fi:

var1:={}:


for t in var do
if op(1,t)=op(2,t) then
 var1:=var1 union {op(1,t)}:
fi:
od:


if nops(var1)<>1 then
RETURN(FAIL):
else
A:=subs({seq(t=1,t in var1)},A):
RETURN(A):
fi:


end:

#IsPE1(P,L): is the list of lines L projectively equivalent to P
IsPE1:=proc(P,L) local P1:
P1:=Du1(L):


if FprT(P,P1)=FAIL then
false:
else
 true:
fi:

end:

IsPE:=proc(P,L) local i,L1:

for i from 1 to nops(P) do
 L1:=[op(i..nops(L),L),op(1..i-1,L)]:
 if IsPE1(P,L1) then
  RETURN(i):
 fi:
od:

FAIL:

end:



#Rep1(n,w,R): inputs a pos. integers n and R, and a word in {1,2,..., n/2}, of even length 
#finds whether
#repeating it r times <=R has property that Tw preserves the property of being
#an inscibed polygon. If successful, it returns (the smallest) r, otherwise it returns FAIL
#Try:
#Rep1(6,[2,1],5);
Rep1:=proc(n,w,R) local r,P,Q,w1:
if nops(w) mod 2=1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

P:=RIP(n,200):
w1:=w:
Q:=Tw(P,w):

if IsI(Q) then
 RETURN(w1):
fi:

for r from 2 to R do
 w1:=[op(w1),op(w)]:
 Q:=Tw(Q,w):

  if Q=FAIL then
     RETURN(FAIL):
 fi:
 if IsI(Q) then
   
 Q:=Tw(RIP(n,100),w1):
   if IsI(Q) then
    RETURN(w1):
   fi:

 fi:
od:

FAIL:

end:



#Rep2(n,w,R): inputs a pos. integers n and R, and a word in {1,2,..., n/2}, of even length 
#finds whether
#repeating it r times <=R, and chopping the last letter
# has property that Tw preserves the property of being prog. equiv.
#an inscibed polygon. If successful, it returns (the smallest) r, otherwise it returns FAIL
#Try:
#Rep2(6,[2,1],5);
Rep2:=proc(n,w,R) local r,P,Q,w1:
if nops(w) mod 2=1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

P:=RIP(n,200):
w1:=w:
Q:=Tw(P,[op(1..nops(w)-1,w)]):

if IsPE(P,Q)<>FAIL then
 RETURN([op(1..nops(w)-1,w)]):
fi:

for r from 2 to R do
 w1:=[op(w1),op(w)]:
 Q:=Tw(P,[op(1..nops(w1)-1,w1)]):

  if Q=FAIL then
     RETURN(FAIL):
  fi:
 if  IsPE(P,Q)<>FAIL then
 RETURN([op(1..nops(w1)-1,w1)]):

 fi:
od:

FAIL:

end:






#MatI1(L): inputs 6 points and outputs the matrix that would decide whether they are inscribed in a conic
MatI1:=proc(L) local i:

if nops(L)<>6 then
 RETURN(FAIL):
fi:

[seq([L[i][1]^2,L[i][1]*L[i][2],L[i][2]^2,L[i][1],L[i][2],1],i=1..6)]:

end:

#MatI(L): inputs a list of points and outputs the matrices that would decide whether they are inscribed in a conic
MatI:=proc(L) local i:

if not( type(L,list) and {seq(type(L[i],list),i=1..nops(L))}={true} and 
 {seq(nops(L[i]),i=1..nops(L))}={2}) then
print(`Bad input`):
  RETURN(FAIL):
fi:


[seq( MatI1([op(1..5,L), L[i] ]),i=6..nops(L))]:




end:




#FindSTW1(n,k,R): finds all Schwartz-Tabachnikov words of the First kind (preserving being inscribed) for an  n-gon
#by repeating a word of even length k up to R times. Try:
#FindSTW1(6,2,6);
FindSTW1:=proc(n,k,R) local gu,mu,w,ha:

if not k mod 2=0 then
 print(`second argument must be an even pos. integer`):
 RETURN(FAIL):
fi:

mu:=Milim1(k,trunc((n-1)/2)):

gu:={}:

for w in mu do
 ha:=Rep1(n,w,R):

  if ha<>FAIL then
    gu:=gu union {ha}:
 fi:
od:

gu:

end:

#FindSTW2(n,k,R): finds all Schwartz-Tabachnikov words of the Second kind (leading to projective equivalence)
#by repeating a word of even length k up to R times. Try:
#FindSTW2(6,2,6);
FindSTW2:=proc(n,k,R) local gu,mu,w,ha:

if not k mod 2=0 then
 print(`second argument must be an even pos. integer`):
 RETURN(FAIL):
fi:

mu:=Milim1(k,trunc((n-1)/2)):

gu:={}:

for w in mu do
 ha:=Rep2(n,w,R):

  if ha<>FAIL then
    gu:=gu union {ha}:
 fi:
od:

gu:

end:


#Nisim1(N,k,R): all the Schwartz-Tabachnikov miracles of the first kind for
#inscribed n-gons from n=6 to n=N, by repeating a word of length k up to R times. Try:
#Nisim1(12,2,6);
Nisim1:=proc(N,k,R) local gu,n:
gu:=[]:

for n from 6 to N do

 gu:=[op(gu), [n,FindSTW1(n,k,R)]]:

od:

gu:

end:


#Nisim2(N,k,R): all the Schwartz-Tabachnikov miracles of the second kind for
#inscribed n-gons from n=6 to n=N, by repeating a word of length k up to R times. Try:
#Nisim2(12,2,6);
Nisim2:=proc(N,k,R) local gu,n:
gu:=[]:

for n from 6 to N do

 gu:=[op(gu), [n,FindSTW2(n,k,R)]]:

od:

gu:

end:




#FindSTW3(n,R): finds all Schwartz-Tabachnikov words of the third kind (preserving being inscribed) for an  n-gon
#by repeating [r,1] for r from 2 to trunc((n-1)/2)  up to R times. Try:
#FindSTW3(6,6);
FindSTW3:=proc(n,R) local gu,mu,w,ha,r:


mu:={ seq( [1,r], r=2..trunc((n-1)/2) )}:

gu:={}:

for w in mu do
 ha:=Rep1(n,w,R):

  if ha<>FAIL then
    gu:=gu union {ha}:
 fi:
od:

gu:

end: