# This is for OEIS A357888
# a(n) is the minimal squared length of the longest side of a strictly convex grid polygon of smallest area

# This is a modification of the program written for
# https://oeis.org/A070911

N_target = 60

# A different target value of N can be given on the command line.


"""

We restrict to edges of squared length at most L, for L=1,2,3,....

===================================================
We assume that the lowest point (or one of the two lowest points)
is the origin O=(0,0).

=========================================================

For each pair of points p and q and every k, we remember the smallest
counterclockwise k-gon (0,0), .... , q, p

Restriction:
* The successive points move strictly upwards.
  (Thus, p is the top point, and there is a "long edge" from O to p.)

* By symmetry, we assume that p is in the right half-plane or on the y-axis.

Polygons with one or two horizontal edges will be treated specially.

===========================

In the end we put together a counterclockwise k-gon ending in p
with a counterclockwise k'-gon starting in p or in p-(1,0) and
ending in O or in (-1,0).

Main data structure:
====================
min_gon is a 3-level dictionary.

min_gon[py][px][k] is a list of pairs (f,vol)=((fx,fy),vol), where
the incoming edge f=(fx,fy)=q-p is a primitive vector (pointing upwards, fy>0)
and vol = 2 * min-area of a k-gon ((0,0), .... ,q,p).
In the list, the vectors f are sorted clockwise by direction,
and the corresponding values "vol" are strictly increasing.
(If pairs don't fit this order, we can eliminate one of them, as DOMINATED.)

We fill this dictionary row by row, increasing py.

"vol" is TWICE the area.


"""

# for https://oeis.org/A070911
known_values = {i+3:v for i,v in enumerate((
    1,2,5,6,13,14,21,28,43,48,65,80,103,118,151,174,213,242,289,
    328,387,420,497,548,625,690,783,860,967,1046,1177,1264,1409,
    1498,1671,1780,1955,2078,2279,2444,2651,2824,3051,3240,3493,
    3676,3969,4176,4489,4714,5045,5302,5623,5906,6243,6556,6991,
    7224,7731,8040,8479,8878,9365, 9804,10313,10774,11323,11796,
    12359,12836,13459,13948,14615,15114,15811,16364,17073,17670,
    18431,19024,19833,20436,21321,21968,22901,23518,24527,25270,
    26283,27050,28105,28896,29993,30798,31959,32830,34019,34946,
    36145,37140 # up to n=102
    ))}

#a = ", ".join(str((known_values[n]-n+2)//2) for n in range(3,56))
#print(a);print(len(a));exit()
#b_file = open("b070911.txt","w")
#for k,vol in known_values.items(): print(k,vol,file=b_file)
#b_file = open("b089187.txt","w")
#for k,vol in known_values.items():
#    if k%2==0: print(k//2,vol//2,file=b_file)
#print(file=b_file); b_file.close(); exit()

from math import tan,pi,gcd,sqrt,ceil
import itertools
import sys

if len(sys.argv)>1:
    N_target = int(sys.argv[1])

min_gon = dict()

how_achieved = [None]*(N_target+1)
reached = [None]*(N_target+1)

def print_solution(k,how_achieved,vol):
    how,px,py,k1,k2 = how_achieved
    #print(f"{(k,how,px,py,k1,k2)=}")
    p1 = find_polygon(px,py,k1)
    p1.reverse()
    if how=="1 FLAT":
        top_point = px-1
    else:
        top_point = px
    ppx = (-top_point)%py
    alpha0 = (-top_point - ppx) // py # usually, alpha0==-1
    leftshift = 1 if how=="2 FLAT" else 0
    p2 = [(-x-leftshift,y) for x,y in find_polygon(ppx,py,k2,alpha0)]
    #print(f"-- {how=} {alpha0=}, {p2[0][0]=},{top_point=},{ppx=},{px=},{py=},{k=}")
    assert p2[0][0]==top_point-leftshift
    if how=="DIAGONAL":
        p2 = p2[1:] #eliminate the common point
    last = [(-1,0)] if how=="2 FLAT" else []
    print ("--- a smallest %d-gon,"%k, "area = %d/2:"%vol,
           [(0,0)]+p1+p2+last) # counterclockwise

## Start the computation

for L in itertools.count(1):
  if all(reached[4:]): break
  for py in itertools.count(1): # for py = 1,2,...
    # print(f"{L=} {py=}")
    # print(min_gon)
    # for k in min_gon.items(): print(k)
    min_gon[py] = dict()  
    # proceed row by row
    for px in itertools.count(0):
        if px**2+py**2>L:
            break
        min_gon[py][px]={2:[((px,py),0)]}

    found_something = False
    for k in range(2,N_target):
        collect={}
        # We will first *collect* the contributions to row py
        # for k+1 before sorting and processing them.         
        for qy in range(1,py):
            Delta_y = py-qy # >0
            if Delta_y**2 > L:
                continue
            for qx,di in min_gon[qy].items():
                #print(f"BB {qy=} {qx=} {di=}")
                if k not in di:
                    continue
                sides = di[k]
                assert(sides)
                px = 0
                row_py_finished = False
                for (fx,fy),vol in sides:
                    # sweep clockwise around q; vol is increasing.
                    # simultaneously sweep the point p on
                    # row py from left to right, clockwise around q.
                    while (# vector q->p is counterclockwise from (fx,fy):
                            # Delta_x/Delta_y < fx/fy 
                            (px-qx) * fy < fx * Delta_y ):
                        # store this:
                        Delta_x = px-qx
                        if Delta_x**2+Delta_y**2>L:
                            if px<qx:
                                px += 1
                                continue
                            else:
                                row_py_finished = True
                                break                                
                        if gcd(Delta_x,Delta_y)==1:
                            if px not in collect:
                                collect[px]=[]
                            collect[px].append(
                                (vol + qx*py-qy*px, Delta_x,Delta_y))
                        px += 1
                    if row_py_finished: break

        # print(f"AA {py=}",collect)
        found_something = found_something or collect
        for px,triples in collect.items():
            if not triples:
                continue

            # This is one more crucial piece of the program:
            # We should keep only those incoming edges that are not DOMINATED.
            # If (fx',fy') enters clockwise from (fx,fy) and has larger volume,
            # then it cannot be part of a minimum-area solution.
            # If it has equal volume, it could be useful, but in any situation it
            # could always be replaced by (fx,fy).
            #
            # We solve this but looking at the entries sorted by slope.
            #
            result = []
            prev_fx, prev_fy = -1,0 # ensure that test is true at the first iteration
            prev_vol = -1
            for vol,fx,fy in (sorted(triples)):
                # smallest area first <==> result should be ordered clockwise.
                if (# this is the first time OR (fx,fy) is clockwise from prev:
                    # fx/fy > prev_fx/prev_fy
                      fx*prev_fy > prev_fx*fy):
                    if vol!=prev_vol:                   
                        result.append(((fx,fy),vol))
                    else: # replace the last entry:
                        result[-1] = ((fx,fy),vol)
                        # adjust this treatment when generating ALL solutions:
                        # Then multiple vol-entries should be kept, but SORTED BY SLOPE.
                    prev_fx,prev_fy,prev_vol = fx,fy,vol

            if px not in min_gon[py]:
                min_gon[py][px]={}
            min_gon[py][px][k+1]=result
            if 0:
                print(f"p=({px},{py}) {k=}")
                for t in result: print("  ",t)
    if py>1 and not found_something: break        
    
    # compute new record areas
    improved = [False]*(N_target+1)

    def check_record(k, vol, how, text=""):        
        if k>N_target:
            return
        if vol>0 and vol==known_values[k]:
            if not reached[k]:
                reached[k] = L
            
                how_achieved[k] = how
#            improved[k] = True

    # Combine "right" k1-gons with "left" k2-gons
    for px,di in min_gon[py].items():
        for k1,best1 in di.items():
            _,vol1 = best1[0]
            # A) exact match along long edge (0,0)-(px,py) # turn 180°
            for k2,best2 in min_gon[py][px].items():
                _,vol2 = best2[0]
                check_record(k1+k2-2,
                             vol1+vol2, ("DIAGONAL",px,py,k1,k2))
                check_record(k1+k2,
                             vol1+vol2+2*py, ("2 FLAT",px,py,k1,k2))
            # B) horizontal edge of length 1 at the top
            #ppx1 = px-1
            if px>0:
              for k2,best2 in min_gon[py][px-1].items():
                _,vol2 = best2[0]
                check_record(k1+k2-1, vol1+vol2+py, ("1 FLAT",px,py,k1,k2))
                 #f"EDGE k={k1+k2-2} {k1=} {k2=} {vol1=} {vol2=} p=({px},{py})")

    if 0 and py==20:
        k1=5
        print(*( min_gon[px,py][k1][0][1] for px in range(py)))
        # this sequence is symmetric and varies between 9 and 42.
        exit()
        
    if 0 and py==30:
        li = sorted((k,best[0][1]) for k,best in ( min_gon[25,py].items()))
        for k,it in min_gon[25,py].items():  print(k,it)
        # most vectors are (necessarily) rather short in the y direction,
        # as soon as k becomes larger
        ks = set(k for k,v in li)
        vols = [v for k,v in li]
        kr= range(min(ks),max(ks)+1)
        assert set(kr)==ks
        print(k,kr,vols)
        print(*(v-vols[i] for i,v in enumerate(vols[1:])))
        exit()
        
    if not min_gon[py]: break  
  print("N=%d, L=%d, height %d"%(N_target,L,py) ,
          flush=True)
  print(reached)
  if all(reached[4:]): break

assert reached[3:29] == [  # according to OEIS
2, 1, 2, 2, 5, 2, 5, 5, 5, 5, 10, 5, 10, 5, 13, 10, 13, 10, 13, 13, 17, 13, 17, 13, 25, 17]
  
# http://antiton.de/PolygonalAreas/index.html?(0,0),(3,1),(5,2),(8,4),(9,5)