#!/usr/bin/python3

"""The algorithms from Appendix A-C of the manuscript
"The Largest Contained Quadrilateral and
the Smallest Enclosing Parallelogram of a Convex Polygon"
by Günter Rote, arXiv:1905.11203v2 [cs.CG], June 2019

The programs refer to a global polygon P for their input.
They only return the area of the smallest enclosing parallelogram
or largest contained quadrilateral, without any identifying information
about these objects.

If started as a main program, the program reads polygons in counterclockwise
order from the file "max-triangle/test-polygons", which is supposed to
contain one polygon per line, in the format

       x1 y1 x2 y2 ... xn yn

Rational coordinates are allowed in the format "p/q"; otherwise they
must be integers.
The programs use exact arithmetic.
Such example input polygons can be found with the source
files of the arXiv preprint of Kallus [Kal] and at
https://github.com/ykallus/max-triangle/releases/tag/v1.0

[Kal] Yoav Kallus. A linear-time algorithm for the maximum-area
      inscribed triangle in a convex polygon. Preprint, June
      2017. arXiv:1706.03049.
"""

import fractions

class Point2D:
    "also used as Vector2D"
    def __init__(self,x,y):
        self.x = x
        self.y = y
    def __sub__(self,other):
        return Point2D(self.x-other.x, self.y-other.y)
    def __str__(self): # for printing
        return "({},{})".format(self.x,self.y)

class Polygon:
    def __init__(self,P):
        self.P = P
        self.n = len(P)
    def __len__(self):
        return self.n
    def __getitem__(self,i):
        "indices are taken modulo n"
        return self.P[i % self.n]

def read_polygon(line):
    p2 = line.split()
    P = []

    if "/" in line:
        type = fractions.Fraction
    else:
        type = int
    for i in range(0,len(p2),2):
        x,y = p2[i:i+2]
        P.append(Point2D(type(x),type(y)))
    return Polygon(P)

def area(a,b,c):
    # a,b,c counterclockwise
    return det(P[b]-P[a],P[c]-P[a])/2
def det(u,v):
    return u.x*v.y-u.y*v.x

###########################################
def both(): # Appendix A
  n = len(P)
  a0 = a = 1
  c = 2
  while area(a,a+1,c+1) > area(a,a+1,c) :
    c0 = c = c + 1 #find the point pc with supporting line parallel to (a,a+1)
  next_AC = "A" #The corner A slides on the edge (a,a+1).
  uAC = P[c] - P[a+1] #the direction u where A hits the next vertex)
  b = a
  while area(c,a,b+1) > area(c,a,b):
    b = b + 1 #find the point pb with supporting line parallel to (a,c).
  d = c
  while area(a,c,d+1) > area(a,c,d):
    d = d + 1 #find the other point pd with supporting line parallel to (a,c).
  if area(b,b+1,d+1) <= area(b,b+1,d):
    next_BD = "B" #The parallelogram side b hits an edge of P before d does.
    uBD = P[b+1] - P[b]
  else:
    next_BD = "D" #The parallelogram side d hits an edge of P before b does.
    uBD = P[d] - P[d+1]
  maxarea = 0 #the area of the largest contained quadrilateral
  minarea = None #the area of the smallest enclosing parallelogram
  while True:
#     print (a,b%n,c%n,d%n,a0,c0,n,next_AC,next_BD,minarea)
      if det(uBD , uAC ) >= 0:
          #The parallelogram side b or d touches an edge of P.
          if next_AC == "A":
              v = P[a+1]-P[a]
          else:
              v = P[c+1]-P[c]
          newarea = abs(det(v,P[c]-P[a])*det(uBD,P[d]-P[b])/det(v,uBD))
#         print (newarea, a,b%n,c%n,d%n,n,next_AC,next_BD)
          if minarea is None or newarea<minarea:
              minarea = newarea
          if next_BD == "B": b = b + 1
          else:              d = d + 1
          if area(b,b+1,d+1) <= area(b,b+1,d):
              next_BD = "B" #The parallelogram side b hits an edge of P before d does.
              uBD = P[b+1] - P[b]
          else:
              next_BD = "D" #The parallelogram side d hits an edge of P before b does.
              uBD = P[d] - P[d+1]
              
      else: #The sliding corner A or C reaches a vertex of P.
          if next_AC == "A": a = a + 1
          else:              c = c + 1
          maxarea = max(maxarea, det(P[c]-P[a],P[d]-P[b])/2)
          if area(a,a+1,c+1) <= area(a,a+1,c): #Which of A and C slides on an edge of P ?)
              next_AC = "A" #The corner A slides on the edge (a,a+1).
              uAC = P[c] - P[a+1] #the direction u where A hits the next vertex
          else:
              next_AC = "C" #The corner C slides on the edge (c,c+1).
              uAC = P[c+1] - P[a] #the direction u where C hits the next vertex
      if (a%n, c%n) == (c0%n, a0%n): return (minarea,maxarea)
  
###########################################
def contained4gon(): # Appendix B
  n = len(P)
  _,_,a0 = min((p.y,p.x,i) for i,p in enumerate(P.P))
  a = b = a0
  _,_,c0 = max((p.y,p.x,i) for i,p in enumerate(P.P))
  c = d = c0
  maxarea = 0
  while True:
#      print (a,b,c,d,len(P))
      while area(c,a,b+1) > area(c,a,b):
          b = b + 1 #find the point b with supporting line parallel to (a,c)
      while area(a,c,d+1) > area(a,c,d):
          d = d + 1 #find the other point d with supporting line parallel to (a,c)
      maxarea = max(maxarea, det(P[c]-P[a],P[d]-P[b])/2)
      if area(a,a+1,c+1) <= area(a,a+1,c): #advance (a, c) to the next antipodal pair
          a = a + 1
      else:
          c = c + 1
      if (a%n, c%n) == (c0 , a0): return maxarea
  
###########################################
def enclosing_parallelogram(): # Appendix C  
  c = 2; d = 2; a = 3; n = len(P)
  while area(c,c+1,a+1) > area(c,c+1,a) : #initialization
      a = a + 1 #find the opposite point a with supporting line parallel to (c,c+1)
  minarea = None
  for b in range(1,n+1):
      while area(b,b+1,d+1) > area(b,b+1,d):
          d = d + 1 #update the point d opposite to (b,b+1)
      while area(b,b+1,a) > area(b,b+1,c+1):
          c = c + 1 #search for antipodal pair AC parallel to u, with C on an edge)
          while (area(c,c+1,a+1) > area(c,c+1,a) or
                (area(c,c+1,a+1)== area(c,c+1,a) and area(b,b+1,a+1)>=area(b,b+1,c))):
              a = a + 1 #update the point a opposite to (c,c+1)
      if area(b,b+1,a)>=area(b,b+1,c):
          v = P[c+1]-P[c]
          u = P[b+1]-P[b]
          newarea = abs(det(v,P[c]-P[a])*det(u,P[d]-P[b])/det(v,u))
          if minarea is None or newarea<minarea:
              minarea=newarea
  return minarea

if __name__ == "__main__":
  for N0,line in enumerate( open("max-triangle/test-polygons")):
    P = read_polygon(line)
    print("#"+str(N0+1)+': n='+str(len(P))
          +', P='+",".join(str(x) for x in P.P))
    ma2 = contained4gon()
    mi2 = enclosing_parallelogram()
    mi,ma = both()
    print (mi, mi2, ma, ma2)
    if ma!=ma2 or mi!=mi2:
        print ("error")
        exit(1)