#OK to post homework.
#Jingyang Deng, Feb.9th, 2020, Optional Homework 2.

###PART1 colorize m*n
#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:

#determine whether a point is in a line.
isinline1:=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:

#To determine whether a line intersect (excluding there are only one common point) with a triangle area.
isintertri:=proc(line,tri,m,n):
if (line[1][2]=0 and line[2][2]=0) or (line[1][1]=m/2 and line[2][1]=m/2) or (isinline1([m/5,n/5],line)=true and line[1][1]=0) then
  RETURN(false):
fi:
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:

#To determine whether a line intersect (excluding there are only one common point) with a rectangle area.
isinterrect:=proc(line,rect,m,n):
if (line[1][1]=0 and line[2][1]=0) or (line[1][2]=0 and line[2][2]=0) or (line[1][1]=m/2 and line[2][1]=m/2) or (line[1][2]=n/2 and line[2][2]=n/2) then
  RETURN(false):
fi:
if (isinline1([m/5,n/5],line)=true and line[1][1]=0) then
  RETURN(true):
fi:
if isinter(line,[rect[1],rect[2]])=true or isinter(line,[rect[2],rect[3]])=true or isinter(line,[rect[3],rect[4]])=true or isinter(line,[rect[4],rect[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 (-1/4,1/4)
isinhalfplane:=proc(p,line) local a,b,c,soln:
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,(-a/4+b/4+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,(-a/4+b/4+c)*(a*p[1]+b*p[2]+c)))>0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
else
  if evalf(subs(a=1,(-a/4+b/4+c)*(a*p[1]+b*p[2]+c)))>0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
fi:
end:

#Assign a value to each point which is in tri/rect 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
rotref:=proc(p,m,n) local t,a,b,mid:
a:=evalf(p[1]-m/2):
b:=evalf(p[2]-n/2):
#print(a,b):
if a>0 then
  if b>0 then
    a:=-a:b:=-b:
  else
    a:=-a:
  fi:
else
  if b>0 then
    b:=-b:
  fi:
fi:
#print(a,b):
if m=n then
  t:=evalf(argument(a+b*I)):
  if t<evalf(-3*Pi/4) then
    mid:=a:
    a:=cos(-2*t-3*Pi/2)*a-sin(-2*t-3*Pi/2)*b:
    b:=sin(-2*t-3*Pi/2)*mid+cos(-2*t-3*Pi/2)*b:
  fi: 
fi:
a:=a+m/2:
b:=b+n/2:
evalf([a,b]):
end:

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

#MAIN FUNCTION
#prec: precision. 
#prec>=100 is recommended.
colorize:=proc(m,n,prec) local i,j,v,edge,tri,rect,good,kplot,lst,col,a,b,mid:
a:=min(m,n):
b:=max(m,n):
v:=[]:
for i from 0 to a do
  v:=[op(v),seq([i,j],j=0..b)]:
od:
edge:=[]:
for i from 1 to (nops(v)-1) do
  for j from (i+1) to nops(v) do
    edge:=[op(edge),[v[i],v[j]]]:
  od:
od:
good:=[]:
if a=b then
  tri:=[[0,0],[a/2,0],[a/2,b/2]]:
  for i from 1 to nops(edge) do
    if isintertri(edge[i],tri,a,b)=true then
      good:=[op(good),edge[i]]:
    fi:
  od:
else
  rect:=[[0,0],[a/2,0],[a/2,b/2],[0,b/2]]:
  for i from 1 to nops(edge) do
    if isinterrect(edge[i],rect,a,b)=true then
      good:=[op(good),edge[i]]:
    fi:
  od:
fi:
#print(good):
lst:=[seq([seq([i/prec,j/prec],j=0..b*prec)],i=0..a*prec)]:
col:=[]:
for i from 1 to (a*prec+1) do
  for j from 1 to (b*prec+1) do
    col:=[op(col),mncolor(lst[i][j],a,b,good)]:
  od:
  #print(col[((i-1)*(a*prec+1)+1)..((i-1)*(a*prec+1)+(b*prec+1))]):
od:
col:=ArrayTools[Reshape](Array(col),[(b*prec+1),(a*prec+1)]):
plots[listcontplot](col,filledregions=true,coloring = ["White","DarkViolet"]):
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,eps) 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^(-eps) then
  RETURN(true):
else
  RETURN(false):
fi:
end:

#calculate the intersect point of two lines.
intrs:=proc(line1,line2,m,n) 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({}):
  else
    if t1*line1[1][1]+(1-t1)*line1[2][1]>0 and t1*line1[1][1]+(1-t1)*line1[2][1]<m and t1*line1[1][2]+(1-t1)*line1[2][2]>0 and t1*line1[1][2]+(1-t1)*line1[2][2]<n then
      RETURN({[t1*line1[1][1]+(1-t1)*line1[2][1],t1*line1[1][2]+(1-t1)*line1[2][2]]}):
    else
      RETURN({}):
    fi:
  fi:
fi:
end:

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

#determine whether two infinite lines are the same.
islinere:=proc(line1,line2) local vec1,vec2,n1,n2,a,soln:
a:='a':
vec1:=evalf(line1[1]-line1[2]):
vec2:=evalf(line2[1]-line2[2]):
if not vec1[2]=0 then
  n1:=evalf(vec1[1]/vec1[2]):
else
  n1:=a:
fi:
if not vec2[2]=0 then
  n2:=evalf(vec2[1]/vec2[2]):
else
  n2:=a:
fi:
#print(n1,n2):
#print(evalb(n1=n2)):
if evalb(n1=n2)=false then
  RETURN(false):
elif evalb(n1=n2)=true then
  if evalb(n1=a)=true then
    #print(1):
    if not line1[1][2]=line2[1][2] then
      RETURN(false):
    else
      RETURN(true):
    fi:
  else
    soln:=solve({line1[1][1]+t*vec1[1]=line2[1][1],line1[1][2]+t*vec1[2]=line2[1][2]},{t}):
    #print(line1,line2):
    #print(soln):
    if soln=NULL then
      RETURN(false):
    else
      #print(2):
      RETURN(true):
    fi:
  fi:
else
  RETURN(`??`):
fi:
end:

#find all the interior vertices and edges.
#eps: epsilon
#digit: Digits
findint:=proc(m,n,eps,digit) local v,edge,i,j,k,intv,p,f,edgere:
Digits:=digit:
v:=[]:
for i from 0 to m do
  v:=[op(v),seq([i,j],j=0..n)]:
od:
#print(v):
edgere:=[]:
for i from 1 to (nops(v)-1) do
  for j from (i+1) to nops(v) do
    edgere:=[op(edgere),[v[i],v[j]]]:
  od:
od:
#print(edgere):
edge:=[edgere[1]]:
#print(edge[1]):
for i from 2 to (nops(edgere)) do
  for j from 1 to nops(edge) do
    #print(islinere(edgere[i],edge[j])):
    f:=1:
    if islinere(edgere[i],edge[j])=true then
      f:=0:
      break
    fi:
  od:
  if f=1 then
    edge:=[op(edge),edgere[i]]:
  fi:
  #print(edge):
od:
#print(edge):
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],m,n)):
    #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],eps)=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:
evalf([intv,edge]):
end:

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

#count # of edges
countedge:=proc(m,n,eps,digit) local a,edge,intv,i,j,co:
a:=findint(m,n,eps,digit):
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],eps)=true then
      co:=co+1:
    fi:
  od:
od:
co+2*m+2*n-4:
end:

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

###There are some mistakes for specific m and n when counting the number of vertices, edges and regions.
###I believe that my idea is right, but actually I do not know why the outcome sometimes is different from the correct answer...
###Perhaps because the numerical computation is not accurate enough