#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)`):
print(`dist(p,q), isinline(p,line), intrs(line1,line2), isre(p,q), findint(n), countv(n), countedge(n), countregion(n)`):end:

###PART1 colorize Kn
#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:

###PART 2 find V, E and R
#calculate the distance between p and q
dist:=proc(p,q):
evalf(sqrt((p[1]-q[1])^2+(p[2]-q[2])^2)):
end:

#determine whether a point is in a line.
isinline:=proc(p,line) local s,area,h,d1,d2,d3:
d1:=dist(p,line[1]):
d2:=dist(p,line[2]):
d3:=dist(line[1],line[2]):
s:=(d1+d2+d3)/2:
area:=sqrt(s*(s-d1)*(s-d2)*(s-d3)):
h:=2*area/d3:
if h<10^(-5) then
  RETURN(true):
else
  RETURN(false):
fi:
end:

#calculate the intersect point of two lines.
intrs:=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({}):
else
  assign(soln):
  if (not whattype(t1)=float) or (not whattype(t2)=float) then
    #print('inf'):
    RETURN({}):
  elif (t1>0 and t1<1 and t2>0 and t2<1) then
    #print('yes'):
    RETURN({[t1*line1[1][1]+(1-t1)*line1[2][1],t1*line1[1][2]+(1-t1)*line1[2][2]]}):
  else
    #print('no'):
    RETURN({}):
  fi:
fi:
end:

#determine whether two points are the same.
isre:=proc(p,q):
if evalf(sqrt((p[1]-q[1])^2+(p[2]-q[2])^2))<10^(-5) then
  RETURN(true):
else
  RETURN(false):
fi:
end:

#find all the interior vertices and edges.
findint:=proc(n) local v,edge,i,j,k,intv,p,f:
Digits:=30:
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:
intv:={}:
for i from 1 to (nops(edge)-1) do
  for j from (i+1) to nops(edge) do
    p:=evalf(intrs(edge[i],edge[j])):
    #print(p):
    if evalb(op(p)=NULL)=true then
      :
    else
      if intv={} then
        intv:=intv union p:
      else
        f:=0:
        for k from 1 to nops(intv) do
          #print(p,intv[k]):
          #print(p,intv[k]):
          if isre(op(p),intv[k])=true then
            f:=1:
            #print('same'):
            break
          fi:
        od:
        #print(f):
        if f=0 then
          intv:=intv union p:
          #print(intv):
        fi:
      fi:
    fi:
    #print(intv):
  od:
od:
[intv,edge]:
end:

#count # of vertices
countv:=proc(n) local a:
a:=findint(n):
nops(a[1])+n:
end:

#count # of edges
countedge:=proc(n) local a,edge,intv,i,j,co:
a:=findint(n):
intv:=a[1]:
edge:=a[2]:
co:=nops(edge):
for i from 1 to nops(edge) do
  for j from 1 to nops(intv) do
    if isinline(intv[j],edge[i])=true then
      co:=co+1:
    fi:
  od:
od:
co:
end:

#calculate # of regions
countregion:=proc(n):
1+countedge(n)-countv(n):
end: