#OK to post homework.
#Jingyang Deng, Feb.8th, 2020, Optional Homework 1.

Help:=proc():print(`isintri(p,tri), isinter(line1,line2), isintertri(line,tri), isinhalfplane(p,line), partcolor(p,good), rotref(p,n), kncolor(p,n,good), colorize(n)`):end:

#To determine whether a point p is in a triangle area.
isintri:=proc(p,tri) local a,b,c,soln:
soln:=solve({a+b+c=1,a*tri[1][1]+b*tri[2][1]+c*tri[3][1]=p[1],a*tri[1][2]+b*tri[2][2]+c*tri[3][2]=p[2]},{a,b,c}):
soln:=evalf(soln):
assign(soln):
#print(a,b,c):
if a>=0 and b>=0 and c>=0 then
  RETURN(true):
else
  RETURN(false):
fi:
end:

#To determine whether two segments intersect with each other (excluding endpoints).
isinter:=proc(line1,line2) local t1,t2,l11,l12,l21,l22,soln:
t1:='t1':t2:='t2':
l11:=t1*line1[1][1]+(1-t1)*line1[2][1]:
l12:=t1*line1[1][2]+(1-t1)*line1[2][2]:
l21:=t2*line2[1][1]+(1-t2)*line2[2][1]:
l22:=t2*line2[1][2]+(1-t2)*line2[2][2]:
soln:=solve({l11=l21,l12=l22},{t1,t2}):
#print(soln):
soln:=evalf(soln):
#print(soln):
if soln=NULL then
  #print('NULL'):
  RETURN(false):
else
  assign(soln):
  if (not whattype(t1)=float) or (not whattype(t2)=float) then
    #print('inf'):
    RETURN(true):
  elif (t1>0 and t1<1 and t2>0 and t2<1) then
    #print('yes'):
    RETURN(true):
  else
    #print('no'):
    RETURN(false):
  fi:
fi:
end:

#To determine whether a line intersect (excluding there are only one common point) with a triangle area.
isintertri:=proc(line,tri):
if isinter(line,[tri[1],tri[2]])=true or isinter(line,[tri[2],tri[3]])=true or isinter(line,[tri[3],tri[1]])=true then
  #print(1):
  RETURN(true):
else
  #print(0):
  RETURN(false):
fi:
end:

#Given a point and a line. Determine whether the point is in the halfplane that contains (0,0)
isinhalfplane:=proc(p,line) local a,b,c,soln:
if line=[[1,0],[-1,0]] then
  RETURN(0):
fi:
a:='a':b:='b':c:='c':
soln:=solve({a*line[1][1]+b*line[1][2]+c=0,a*line[2][1]+b*line[2][2]+c=0},{a,b,c}):
#print(soln):
assign(soln):
if whattype(c)=symbol then
  if evalf(subs(c=1,c*(a*p[1]+b*p[2]+c)))>0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
elif whattype(b)=symbol then
  if evalf(subs(b=1,c*(a*p[1]+b*p[2]+c)))>0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
else
  if evalf(subs(a=1,c*(a*p[1]+b*p[2]+c)))>0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
fi:
end:

#Assign a value to each point whose argument is in [0,2pi/5] by 
#calculating the sum of some characteristic functions of halfplanes.
#good: a set of lines that intersect with triangle.
partcolor:=proc(p,good) local i,a:
#print(good):
a:=add(isinhalfplane(p,good[i]),i=1..nops(good)):
if a mod 2=0 then
  a:=-a:
fi:
a:
end:

#Rotation and reflection, which enable a point converted to 
#the point whose argument is in [0,2pi/5].
#This step ensures that the whole picture is symmetric.
rotref:=proc(p,n) local t,m,i,tri,tri1,a,b:
#print(p):
#print(p[1]):
a:=p[1]:
b:=p[2]:
tri:=[[0,0],[1,0],[cos(Pi/n)^2,cos(Pi/n)*sin(Pi/n)]]:
tri1:=[[0,0],[1,0],[cos(2*Pi/n),sin(2*Pi/n)]]:
#print(isintri(p,tri1)):
for i from 0 to (n-1) do
  t:=evalf(argument(a+b*I)):
  if t>=0 and t<=evalf(2*Pi/n) then
    break
  else
    m:=a:
    a:=cos(2*Pi/n)*a-sin(2*Pi/n)*b:
    b:=sin(2*Pi/n)*m+cos(2*Pi/n)*b:
  fi:
od:
if t>evalf(Pi/n) then
  m:=a:
  a:=cos(2*Pi/n-2*t)*a-sin(2*Pi/n-2*t)*b:
  b:=sin(2*Pi/n-2*t)*m+cos(2*Pi/n-2*t)*b:
fi:
[a,b]:
end:

#To assign a value to any point in the plane.
kncolor:=proc(p,n,good):
partcolor(rotref(p,n),good):
end:

#MAIN FUNCTION
#First, figure out the expressions of vertices, edges. 
#Define a specific triangle tri whose area is 1/2n of the polygon. 
#Assign values in this triangle. 
#Use rotation and reflection to define values in other regions.
#Make a contourplot.
###This algorithm is simple but expensive to run, so I just used listcontplot instead.
###However, a very long time is still needed.
###m denotes the precision of the outcome. Larger m means a more clear the colored graph is.
colorize:=proc(n,m) local i,j,v,edge,tri,tri1,good,kplot,lst,col:
v:=[seq([cos(2*Pi*i/n),sin(2*Pi*i/n)],i=0..(n-1))]:
#print(v):
edge:=[]:
for i from 1 to (n-1) do
  for j from (i+1) to n do
    edge:=[op(edge),[v[i],v[j]]]:
  od:
od:
tri:=[[0,0],[1,0],[cos(Pi/n)^2,cos(Pi/n)*sin(Pi/n)]]:
tri1:=[[0,0],[1,0],[cos(2*Pi/n),sin(2*Pi/n)]]:
#print(tri1):
good:=[]:
for i from 1 to nops(edge) do
  if isintertri(edge[i],tri)=true and (not edge=[[1,0],[-1,0]]) then
    good:=[op(good),edge[i]]:
  fi:
od:
#print(good):
lst:=[seq([seq([i/m,j/m],j=-m..m)],i=-m..m)]:
#print(lst):
col:=[]:
for i from 1 to (2*m+1) do
  for j from 1 to (2*m+1) do
    #print((i-1)*(2*m+1)+j,kncolor(lst[i][j],n,good)):
    col:=[op(col),kncolor(lst[i][j],n,good)]:
  od:
  #print(col[(i-1)*(2*m+1)+1..i*(2*m+1)]):
od:
col:=ArrayTools[Reshape](Array(col),[(2*m+1),(2*m+1)]):
plots[listcontplot](col,filledregions=true,coloring = ["Blue","Red"]):
end: