######################################################################
##GeneralizedMorley: Save this file as  GeometryMiracles             #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read GeneralizedMorley<Enter>                      #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: May-June, 2014
 
print(`Created:  May-June 2014`):
print(` This is GeneralizedMorley `):
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 MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------------------`):




ezra:=proc(): 

if nargs=0 then
print(`Contains procedures: DisGM, DisGMv , GMorleyDis `):

elif nargs=1 and args[1]=DisGM then
print(`DisGM(A12,A23,A31,r,d): tries to find all polynomials of degree <=d in the distance-squared of the r-sector generalized Morely triangle`):
print(`The generalized Morley (the original one is r=3) is obtained by drawing r-sectors at each of the three vertices`):
print(`and taking the points of intersection of these r-sectors that are closest to the edges of the original trianle`):
print(`and calling the lenght-sides squares, A12,A13,A23. To rediscover (and rigorously prove!) Morley's theorem `):
print(`Try: `):
print(`DisGM(A12,A23,A31,3,1); `):


elif nargs=1 and args[1]=DisGMv then
print(`DisGMv(A12,A23,A31,r,d): verbose version of DisGM(A12,A23,A31,r,d)`):
print(`Try: `):
print(`DisGMv(A12,A23,A31,3,1); `):

elif nargs=1 and args[1]=GMorleyDis then
print(`The length^2 of the sides of the general Morley triangle with r-sectors, in terms of m and n`):
print(`GMorleyDis(m,n,r)`):
print(`The length^2 of the sides of the general Morley triangle with r-sectors, in terms of m and n`):
print(`Try: `):
print(`GMorleyDis(m,n,3);`):
else
 print(`There is no ezra for`,  args ):
fi:

end:

#ez:=proc(): print(`Morley(), MorleyDis(m,n), GMorleyDis(m,n,r), GenPol(Resh,a,deg) `): 
#print(` DisGM1(A12,A23,A31,r,d), DisGM(A12,A23,A31,r,d) `):
#end:

 

#Def (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:tan(a_1+a_2+...) expressed in terms li[1]:=tan(a_1), li[2]:=tan(a_2) ..
TS:=proc(li) local i,t:if nargs=1 then RETURN(args[1]) else
t:=TS(seq(args[i],i=2..nargs)): RETURN((args[1]+t)/(1-t*args[1])):fi:end:

#Def (The square of the distance of points A and B)
DeSq:=proc(A,B):normal((A[1]-B[1])^2+(A[2]-B[2])^2): end:

#Def(Is the triangle ABC equilateral?)
ItIsEqui:=proc(A,B,C):evalb(normal
({DeSq(A,B)-DeSq(A,C),DeSq(B,C)-DeSq(C,A)})={0}):end:

Morley:=proc() local m,n,C,D,E,F:
C:=Pt(y-TS(m,m,m)*x,y+TS(n,n,n)*(x-1)):D:=Pt(y-m*x,y+n*x-n):
E:=Pt(y-TS(m,m)*x,y-C[2]-(x-C[1])*TS(m,m,-n,sqrt(3))):
F:=Pt(y+TS(n,n)*(x-1),y-C[2]+(x-C[1])*TS(n,n,-m,sqrt(3))):ItIsEqui(D,E,F):end:


#The length^2 of the sides of the Morley triangle in terms of m and n
MorleyDis:=proc(m,n) local C,D,E,F:
C:=Pt(y-TS(m,m,m)*x,y+TS(n,n,n)*(x-1)):D:=Pt(y-m*x,y+n*x-n):
E:=Pt(y-TS(m,m)*x,y-C[2]-(x-C[1])*TS(m,m,-n,sqrt(3))):
F:=Pt(y+TS(n,n)*(x-1),y-C[2]+(x-C[1])*TS(n,n,-m,sqrt(3))):
[DeSq(D,E),DeSq(E,F),DeSq(F,D)]:
end:


#The length^2 of the sides of the general Morley triangle with r-sectors, in terms of m and n
GMorleyDis:=proc(m,n,r) local C,D,E,F,ka:

if not type(r,integer) and r>1 then
 RETURN(FAIL):
fi:

if r=2 then
C:=Pt(y-TS(m$r)*x,y+TS(n$r)*(x-1)):D:=Pt(y-m*x,y+n*x-n):
E:=Pt(y-TS(m$(r-1))*x,y-C[2]+(x-C[1])/TS(m$(r-1),-n)):
F:=Pt(y+TS(n$(r-1))*(x-1),y-C[2]-(x-C[1])/TS(n$(r-1),-m)):
RETURN([DeSq(D,E),DeSq(E,F),DeSq(F,D)]):
else
ka:=tan(Pi/r):
C:=Pt(y-TS(m$r)*x,y+TS(n$r)*(x-1)):D:=Pt(y-m*x,y+n*x-n):
E:=Pt(y-TS(m$(r-1))*x,y-C[2]-(x-C[1])*TS(m$(r-1),-n,ka)):
F:=Pt(y+TS(n$(r-1))*(x-1),y-C[2]+(x-C[1])*TS(n$(r-1),-m,ka)):
[DeSq(D,E),DeSq(E,F),DeSq(F,D)]:
fi:

end:

 
 

#GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a,
#and a non-negative integer deg, outputs a generic polynomial
#of total degree deg,in the variables Resh, indexed by the letter a of
#followed by the set of coeffs. For example, try:
#GenPol([x,y],a,1);
GenPol:=proc(Resh,a,deg)
local var,POL1,i,x,Resh1,POL:

if not type(Resh,list) or nops(Resh)=0 then
 ERROR(`Bad Input`):
fi:

x:=Resh[1]:
if nops(Resh)=1 then
RETURN(convert([seq(a[i]*x^i,i=0..deg)],`+`),{seq(a[i],i=0..deg)}):
fi:

Resh1:=[op(2..nops(Resh),Resh)]:
var:={}:
POL:=0:

for i from 0 to deg do
POL1:=GenPol(Resh1,a[i],deg-i):
var:=var union POL1[2]:
POL:=expand(POL+POL1[1]*x^i):
od:
POL,var:
end:



#GenHomPol(Resh,a,deg): Given a list of variables Resh, a symbol a,
#and a non-negative integer deg, outputs a generic Homog. polynomial
#of total degree deg,in the variables Resh, indexed by the letter a of
#followed by the set of coeffs. For example, try:
#GenHomPol([x,y],a,1);
GenHomPol:=proc(Resh,a,deg)
local var,POL1,i,x,Resh1,POL:

if not type(Resh,list) or nops(Resh)=0 then
 ERROR(`Bad Input`):
fi:

x:=Resh[1]:
if nops(Resh)=1 then
RETURN(a[deg]*x^deg,{a[deg]}):
fi:

Resh1:=[op(2..nops(Resh),Resh)]:
var:={}:
POL:=0:

for i from 0 to deg do
POL1:=GenHomPol(Resh1,a[i],deg-i):
var:=var union POL1[2]:
POL:=expand(POL+POL1[1]*x^i):
od:
POL,var:
end:



#DisGM1(A12,A23,A31,r,d): tries to find a polynomial of degree d in the distance-squared of the r-sector generalized Morely triangle
#Try: DisGM1(A12,A23,A31,3,1); A12,A13,A23 are the distance-squared of that triangle
DisGM1:=proc(A12,A23,A31,r,d) local a,P,eq,var,lu,m,n,ra,m1,n1,eq1,var1,lu1,
tov,var11,tov1,fu,fu1,ku:

P:=GenPol([A12,A23,A31],a,d):
var:=P[2]:
P:=P[1]:


lu:=GMorleyDis(m,n,r):

if lu=[0,0,0] then
RETURN(lu):

fi:
eq:={}:

ra:=rand(100..200):

while nops(eq)<nops(var)+10 do

m1:=(ra()+1)/(ra()+3):
n1:=(ra()+1)/(ra()+4):
lu1:=subs({m=m1,n=n1},lu):
eq1:=expand(subs({A12=lu1[1],A23=lu1[2],A31=lu1[3]}, P)):
eq:=eq union {eq1}:
od:

var1:=solve(eq,var):

if var1=NULL or subs(var1,P)=0 then
 RETURN(FAIL):
fi:

tov:={}:

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

P:=subs(var1,P):

fu:={seq(coeff(P,tov1,1), tov1 in tov)}:
tov:={}:

for fu1 in fu do
 ku:=normal(subs({A12=lu[1],A23=lu[2],A31=lu[3]},fu1)):
 if ku<>0 then
  print(fu1, `did not work out`):
 else
  tov:=tov union {fu1}:
 fi:
od:

tov:   

end:



#DisGMhom1(A12,A23,A31,r,d): tries to find a Homg. polynomial of degree d in the distance-squared of the r-sector generalized Morely triangle
#Try: DisGMhom1(A12,A23,A31,3,1); A12,A13,A23 are the distance-squared of that triangle
DisGMhom1:=proc(A12,A23,A31,r,d) local a,P,eq,var,lu,m,n,ra,m1,n1,eq1,var1,lu1,tov,var11,tov1:

P:=GenHomPol([A12,A23,A31],a,d):
var:=P[2]:
P:=P[1]:


lu:=GMorleyDis(m,n,r):

if lu=[0,0,0] then
RETURN(lu):

fi:
eq:={}:

ra:=rand(100..1000):

while nops(eq)<nops(var)+10 do

m1:=ra()/ra():
n1:=ra()/ra():
lu1:=subs({m=m1,n=n1},lu):
eq1:=expand(subs({A12=lu1[1],A23=lu1[2],A31=lu1[3]}, P)):
eq:=eq union {eq1}:
od:


var1:=solve(eq,var):

if var1=NULL or subs(var1,P)=0 then
 RETURN(FAIL):
fi:

tov:={}:

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


P:=subs(var1,P):


{seq(coeff(P,tov1,1), tov1 in tov)}:

end:



#DisGM(A12,A23,A31,r,d): tries to find a polynomial of degree <=d in the distance-squared of the r-sector generalized Morely triangle
#Try: DisGM(A12,A23,A31,3,1); A12,A13,A23 are the distance-squared of that triangle
DisGM:=proc(A12,A23,A31,r,d) local d1,gu:

for d1 from 1 to d do

gu:=DisGM1(A12,A23,A31,r,d1):

if gu<>FAIL then
 RETURN(gu):
fi:

od:

FAIL:

end:



#DisGMhom(A12,A23,A31,r,d): tries to find a homog. polynomial of degree d in the distance-squared of the r-sector generalized Morely triangle
#Try: DisGMhom(A12,A23,A31,3,1); A12,A13,A23 are the distance-squared of that triangle
DisGMhom:=proc(A12,A23,A31,r,d) local d1,gu:

for d1 from 1 to d do

gu:=DisGMhom1(A12,A23,A31,r,d1):

if gu<>FAIL then
 RETURN(gu):
fi:

od:

FAIL:

end:






#DisGMv(A12,A23,A31,r,d): verbose version of DisGM(A12,A23,A31,r,d);
#Try: DisGMv(A12,A23,A31,3,1);
DisGMv:=proc(A12,A23,A31,r,d) local gu,t0,i:
t0:=time():
gu:=DisGM(A12,A23,A31,r,d):

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

if gu=[0,0,0] then 


print(`Theorem: Take any triangle, and  draw `, r ,  `-sectors at each vertex. `):
print(`The Generalized Morley triangle has as its vertices the points of intersection of`):
print(r, ` -sectors closest to the edges emanating from the two vertices of each edge `):
print(`These are concurrent.`):
else


print(`Theorem: Take any triangle, and  draw `, r ,  `-sectors at each vertex. `):
print(`The Generalized Morley triangle has as its vertices the points of intersection of`):
print(r, ` -sectors closest to the edges emanating from the two vertices of each edge `):
print(`Let `, A12, A23, A31 , `be the DISTANCES-SQUARED, of that new triangle`):
print(`divided by the square of the length  of the base( i.e. the side closest to the intersection`):
print(`of edges A12 and A31 of the generalized Morley  triangle. `):

if nops(gu)>1 then

print(`We have the following polynomial relations, of degree`, d, `relating them.`):

for i from 1 to nops(gu) do
lprint(numer(gu[i])=0):
od:

else

print(`We have the following polynomial relation, of degree`, d, `relating them.`):
lprint(numer(gu[1])=0):
fi:
fi:

print(`This ends this cute article that took`, time()-t0, `seconds to generate. `):
end: