N_target = 44
EVEN_N = True # look only for even n-gons
#EVEN_N = False # looking at all n

# A different target value of N can be given on the command line.
# "even" or "all" can also be given on the command line.

"""
Lemma 1: Every edge is a primitive vector.

Lemma 2: If n is even, then P can be assumed to be centrally symmetric.

The correctness of the program (knowing when one can stop) depends
crucially on two lemmas that relate area to lattice width. The first lemma is about centrally symmetric polygons (arising when n is even), and
the second lemma about arbitrary polygons.

Definition:
The *lattice width* of a lattice polygon P *in a given lattice direction*
(a,b) in Z² is one less than the number of lattice lines orthogonal to (a,b)
intersecting P.
The *lattice width* of P is the smallest number lattice width over
all lattice directions.

Lemma 3: A centrally symmetric convex polygon P of lattice width w has area
at least w²/2.

Proof:
[ A weaker bound of w²/4 is given in the following paper
(p.175, before remark 2):
Imre Bárány and Norihide Tokushige, The minimum area of convex lattice n-gons.
Combinatorica, 24 (No. 2, 2004), 171-185.
https://www.renyi.hu/~barany/cikkek/94.pdf
Our proof uses the same setup but refines the argument in the end.
The improvement to w²/2 speeds up the computations by a factor 5-6. ]


1. We assume that the lattice width direction is vertical.
Thus P touches two horizontal lines H+: y=b and H-: y=-b, with b=w/2.
Let (a,0) be the intersection of P with the positive x-axis.

2. If a>=b, then the area is at least 2b²=w²/2, and we are done.
Thus, let us assume a<b.

3. By horizontal shearing, we assume that the tangent at (a,0)
(or some supporting line) is not vertical and has slope s at least +1.
(Since horizontal lattice-preserving shearing changes the "inverse
slope" 1/s in increments of 1, the inverse slope can always be brought into
the interval (0,1].)

4. The lattice width in direction (1,0) is b'>=b. Thus
P touches two vertical lines V+: x=b' and V-: x=-b'

5. The lattice width in direction (1,-1) is b">=b. Thus
P touches two diagonal lines D+: y-x=b" and D-: y-x=-b"
(Draw a picture!)

6. If (a,0) had a tangent of slope 1, P could not touch D-, since a<b".

Thus the slope s is strictly larger than 1, and P must touch V+
in the upper half-plane, and P must touch D- in the lower half-plane.

7. If the intersection of V+ and D- is below the x-axis,
we swap the roles of V+ and D- by a horizontal shearing transformation
combined with a vertical reflection: Each point (x,y) is mapped to (x-y,-y).
This turns V+ into a line of slope 1 and makes D- vertical.
Point on the x-axis are fixed.

Thus we can assume that the intersection of V+ and D- is above or
on the x-axis, at (b',c) with c>=0.

8. Pick a point where P touches each of the six lines, in a symmetric way.
Denote them by B,BR,R (Bottom, Bottom-Right, and Right) and the reflected
points by T,TL,L. (Draw a picture!)
Replace P by the convex hull of these points (making the area smaller).

9. Holding L,B,R, and T fixed, we can move the points BR and TL on their
respective lines D- and D+, keeping them symmetric.
The area depends linearly on the movement.

If the are would decrease when BR is moved towards R and towards V+,
we can stop this movement at the point (b",0) where D- intersects
the x-axis. The convex hull of (b",0),(-b",0),T,B has area 2b"b>w²/2.

10. Thus we are left with the case that we can decrease the area by
moving BR towards H-, reaching the point B'=(d,-b) with d := b"-1>=0.
Then we consider the convex hull of B', R, and their two
symmetric points. R has coordinates R=(b',e) with e>=c>=0.
The determinant of (d,-b) and (b',e) is bb'+de>=bb'>=b², and we are done.

QED.


Lemma 4: A (not necessarily symmetric) convex polygon P of lattice width w
has area at least w²/3.
The inequality is strict except when P is a triangle.

"""
# Since we use "volume", which is TWICE the area, the following constants
# are twice as big as the constants in the lemmas.

lattice_width_factor_even = 1
lattice_width_factor_odd  = 2/3 + 0.00001
# +0.00001 is because of the strict inequality (except for triangles, but
# those are found in the first iteration.)

def lower_bound(k,w):
    """lower bound on vol (=2*area) for k-gon of lattice width w """
    if k%2:
        lattice_width_factor = lattice_width_factor_odd
    else:
        lattice_width_factor = lattice_width_factor_even
    LB = ceil( w**2 * lattice_width_factor )
    if k%2 != LB%2:
        LB +=1 # vol must be odd for odd k and even for even k
    return LB

"""
Proof of Lemma 4:

The polygon P and the symmetrized version Q = 1/2 * (P-P)
(which is half the "difference body" P-P={ a-b | a,b in P }),
have the same width in every direction, and Q is centrally
symmetric. We can apply Lemma 3 to Q.
The area of Q is at most 3/2 times the area of P, and
equality holds only when P is a triangle.
See any of the following papers. QED

Hans Rademacher, Über den Vektorenbereich eines konvexen ebenen Bereiches.
Jahresbericht d. Deutschen Math.-Verein. 34 (1926): 64-78. http://eudml.org/doc/145705
Theodor Estermann, Zwei neue Beweise eines Satzes von Blaschke und Rademacher.
Jahresbericht d. Deutschen Math.-Verein. 36 (1927): 197-200. http://eudml.org/doc/145763
Theodor Estermann, Über den Vektorenbereich eines konvexen Körpers.
Math Z. 28 (1928), 471-475. https://doi.org/10.1007/BF01181177
C. A. Rogers and G. C. Shephard, The difference body of a convex body.
Arch. Math. 8 (1957), 220-233. doi:10.1007/BF01899997


I believe the true factor for Lemma 4 might be 3/8, and it is sharp for
the triangle (-w/2,-w/2),(w/2,0),(0,w/2), for the case of even w.
(By analogy with the ordinary width, not lattice width:
"It is well-known that of all convex sets of a given width, the
equilaterial triangle has the smallest area.")


In the proof, we would consider the 8-gon formed as the
intersection of a horizontal slab of width w, a vertical slab
of width w1>=1, and two slabs of directions 45° and -45° of
*vertical* widths w3>=w and w4>=w, respectively. (Their Euclidean
widths are >=w/sqrt(2).)

Some sides of the 8-gon may degenerate into
a point, in case several lines are concurrent.
The body must touch each of these edges (even if they are degenerated to a point).

Lemma 5.
A) For three successive vertices v,v',v" of P, the triangle vv'v" contains
no interior point.
B) If n is even, then there cannot even be points on the diagonal vv".

Proof:
A) Otherwise v' could be replaced by that point, reducing the area.

B) Bárány and Tokushige proved that the set E of edge vectors is
convex in the following sense: Every vector e in the convex hull of E
that is a primitive vector, must belong to E.
Let e=v->v' and e'=v'->v". If vv"=e+e' would have gcd d>2 the vector
(e+e')/d, which is in the convex hull of e,e',-e,-e', is missing
from E, a contradiction.

This Lemma is used to shortcut the propagation process.
Property B is really strong. It does not hold for the optimal 5-gons, 7-gons,
and 11-gons, but by preliminary experiments,
it seems that for larger k, it might hold also in the odd case.

Property B means the following. Consider a pair (q,f). Then the
outgoing edges q->p that need to be considered end on the next
lattice line parallel to f. The points q form an arithmetic progression
p0, p0+f, p0+2f, ....
When considering a range of diameter D (roughly the maximum height),
then these are about D/|f| successor points altogether (on all levels py).

Under Property A, we have to consider in addition all multiples of the vectors
(q-f)->(p0+t*f).
I.e.,
p = (q-f) + s*(p0+t*f-(q-f))
  = s*p0 - (s-1)*q + (s*t-s-1)*f
This would be D/|f| log(D*|f|) ~ D log D/|f| points.
Typically, the average length of the lists (the number of
vectors f for a point q) is very small (about 4 for height py=100),
and |f| is not always short.

(By contrast, currently, every point p gets O(D²) inputs, one from every other
point q (looking at the propagation from the incoming side).)

Lemma 6:
If n is even, then P has two edges parallel to the lattice width direction

Proof:
Suppose not.
Then there is a pair of extreme VERTICES where the strip of minimum lattice width
touches P. (as opposed to a pair of extreme EDGES.)
Let v be such a vertex, and let v- and v+ its neighbors.
As usual, we assume that the lattice width direction is horizontal,
as we assume that v=(0,0).
Let u- and u+ be the intersection of the edges v -> v- and v -> v+ with the line y=1.
The open interval (u-,u+) cannot contain an integer point, by Lemma 5B.
Hence, by horizontal shearing we can assume 0 <= u- < u+ <= 1.
P is contained in the wedge between v -> v- and v -> v+.
Hence the coordinates highest point v'=(x,y) (opposite to v) fulfill the inequality 0 <= x <= y.
If x<y, then the width in vertical direction is smaller than the horizontal width (which is y).
If x=0, then the width in direction (1,1) is smaller than y.
Both cases lead to a contradiction. QED

===================================================
Also in the algorithm (as in the proof of Lemma 3), we make the assumption
that the lattice width direction is vertical.
(The smallest number of parallel lattice lines intersecting P
are the horizontal lines.)

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

The approach is by dynamic programming. It is the
algorithm in the following paper:

Finding minimum area k-gons.
David Eppstein, Mark Overmars, Günter Rote und Gerhard Woeginger.
Discrete and Computational Geometry 7 (1992), 45-58.
https://doi.org/10.1007/BF02187823

see also
Counting convex polygons in planar point sets.
Joseph S. B. Mitchell, Günter Rote, Gopalakrishnan Sundaram,
and Gerhard Woeginger.
Information Processing Letters 56 (1995), 45-49
http://page.mi.fu-berlin.de/rote/Papers/abstract/Counting+k-subsets+and+convex+k-gons+in+the+plane.html

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

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.)

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

By horizontal lattice-preserving shearing transformations, it suffices to
store results only for the range 0 <= px < py
(The answer for any point p in the plane can then be figured out
by an appropriate shearing transformation.)

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

Main data structure:
====================
min_gon is a dictionary
min_gon[(px,py)] is a dictionary

min_gon[(px,py)][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.

[ If we want to enumerate ALL optimal solutions, we should allow
the "vol" values to be weakly increasing. We should then use the
weaker lower bound also for even k, in case we are interested in solutions
that are not centrally symmetric.
Another possibility is COUNTING the number of solutions. (For each height
separately. Every solution has up to 2k representations with at least
one flat edge, and infinitely many without flat edge.)
We should consider all solutions of height up until the value is confirmed. ]
[ We should have proceeded from left to right, swapping x and y.
Then the notion of slope could have been used more naturally. ]

"vol" is TWICE the area.

https://oeis.org/A070911
1,2,5,6,13,14,21,28,43,48,65,80,103,118, / 151,174,213,242,289,328,387,420,497
The odd values after the slash were unconfirmed according to OEIS, as of August 2023.

POTENTIAL IMPROVEMENTS:
In the n-gons computed, the left and right halves have always equal size +-1, for odd n?
So k needs to go only to n/2, in practice.
(This is definitely warranted for even-n-gons).
Can one show that the difference is at most 1?
Is the "best"-list for a given point a convex function of k? 
[ NO: This is very much not the case, not even in the middle of the range:
Example: (px,py) = (20,29), k in range(2, 31)
Best vols:[0,5,6,13,14,29,64,81,122,171,210,255,334,421,502,641,746,865,1072,1235,1438,1657,1924,2175,2440,2821,3214,3631,4096]
Differences are far from monotone:
5 1 7 1 15 35 17 41 49 39 45 79 87 81 139 105 119 207 163 203 219 267 251 265 381 393 417 465 ]

The program should have an option to evaluate only the even
ones; they are faster anyway because of the better constant factor in
the lower bound, AND for the above reason.

For computing best paths to the top vertex (without constraints on the
incoming slopes, when we take best1[0] in the final combination), one
could do a meet-in-the-middle divide-and-conquer approach.
Partition into k/2 + k/2, and try all meeting points.
-> reduce range of k by 1/2; factor-2 improvement

Another big opportunity for improvement would be to restrict
the range of outgoing points p that can be meaningfully reached from
a given q.

"""

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
    ))} # optional, just for guidance; not used by the algorithm
known_values = {}    

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

if len(sys.argv)>1:
    N_target = int(sys.argv[1])
for s in sys.argv[1:]:
    if s.lower() == "even":
        EVEN_N = True
    if s.lower() == "all":
        EVEN_N = False
    
    

if EVEN_N:
    b_file = open("b089187.txt","w")
    a_file = open("a089187.txt","w")
    output_so_far = 1
else:
    b_file = open("b070911.txt","w")
    a_file = open("a070911.txt","w")
    output_so_far = 2 # the values are not confirmed in order
    # because the even ones come out earlier.
    
    
min_gon = dict()
record_vol = [n**3 for n in range(N_target+1)] # loose upper bound as a start value
min_height = [0 for n in range(N_target+1)] # smallest height of a min-area polygon
max_height = [0 for n in range(N_target+1)]
# largest height of a min-area polygon with at least one horizontal edge,
# AS FOUND BY THE PROGRAM, not definite. The program does not look for the largest heights.

k_START = 2 if EVEN_N else 3 # the values below k_start are trivial

confirmed = [True]*k_START + [False]*(N_target+1-k_START) # checks which entries are confirmed.
how_achieved = [None]*(N_target+1)

# Auxiliary procedures to construct the solution, once the optimal
# area has been determined
def find_polygon(px,py,k,alpha0=0):
    # find the smallest k-gon ending in (px,py) by backtracing
    # 0<px<py. The startpoint (0,0) is not produced.
    alphaT = alpha0
    result = [(px+alphaT*py,py)]
    (fx,fy),vol = min_gon[px,py][k][0]
    for k in range (k-1,1,-1):
        qx,qy = px-fx,py-fy
        vol -= qx*py-qy*px
        alpha =  qx//qy
        px,py = qx-alpha*qy,qy
        alphaT += alpha # alpha is accumulated
        result.append((px+alphaT*py,py))
        for (fx,fy),vol2 in min_gon[px,py][k]:
            if vol2==vol: # this must be the right entry, recognition by vol is easiest
                break
        else:
            error("not found",k,px,py)
    return result

def print_solution(k,how_achieved,vol,for_b_file=False):
    # parameters k and vol are only for output, have no effect on the procedure.
    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 []
    polygon = [(0,0)]+p1+p2+last # counterclockwise
    if not for_b_file:
        print (f"--- a smallest {2*k if EVEN_N else k}-gon,",               
               f"area = {vol}:" if EVEN_N else f"area = {vol}/2:",
               polygon)
    else:
        result = vol # //2 if EVEN_N else vol
        print(k,result, file=b_file, flush=True)
        print(k,result, ",".join(f"({x},{y})" for x,y in polygon),
              file=a_file, flush=True) # compact format without spaces

## Start the computation

py_target_start = 1
py_target_end=20 # will be iteratively doubled when not sufficient
target_height = 1

while py_target_start<= target_height:
  collect=tuple([[[] for k in range(N_target+1)] for px in range(py)]
                for py in range(py_target_end+1)  )

  print(f"ONE ROUND STARTS: propagate from rows {1}-{py_target_end-1}",
        f"to rows {py_target_start}-{py_target_end}")
  for qy in range(py_target_end):
    # expand from row qy into rows in py_target_start <= py <= py_target_end
    py_lower = max(py_target_start,qy+1)  
    print(f"propagate from row {qy}",
        f"to rows {py_lower}-{py_target_end}", flush=True)

    for qx in range(qy):
        for k,sides in min_gon[qx,qy].items():
            if k>=N_target:
                continue
            for (fx,fy),vol in sides:
                # search first point q-d above qy-fy on the next lattice line
                # This could be done by the extended Euclidean algorithm.
                # Linear search is o.k. in this context.
                for dy in range(fy):
                    dx = fx*dy//fy + 1
                    # assert dx*fy-fx*dy>0
                    if dx*fy-fx*dy==1:
                        break
                #else: assert 0
                #print(f"{(dx,dy)=} {(fx,fy)=}")
                if EVEN_N:
                    px,py = qx+fx-dx,qy+fy-dy
                    while py<=py_target_end:
                        if py>=py_target_start:
                            Delta_x,Delta_y = px-qx,py-qy
                            # gcd(Delta_x,Delta_y)==1 by construction
                            alpha = px//py
                            collect[py][px-alpha*py][k].append(
                                (vol + qx*py-qy*px,
                                 Delta_x-alpha*Delta_y,Delta_y))
                        px,py = px+fx,py+fy
                else:
                    rx,ry = qx-fx,qy-fy
                    k0x,k0y = qx-dx,qy-dy # shoot from r in direction k0
                    while k0y<=py_target_end:
                        stepx,stepy = k0x-rx,k0y-ry
                        px,py = k0x,k0y
                        # the inner loop is necessary for odd n, as far as we know.
                        while py<=py_target_end:
                            if py>=py_lower:
                                Delta_x,Delta_y = px-qx,py-qy
                                if gcd(Delta_x,Delta_y)==1:
                                    alpha = px//py
                                    collect[py][px-alpha*py][k].append(
                                        (vol + qx*py-qy*px,
                                         Delta_x-alpha*Delta_y,Delta_y))
                            px,py = px+stepx,py+stepy
                        k0x,k0y = k0x+fx,k0y+fy
     
    py = qy+1
    # consolidate the lists in collect[py]
    if py<py_target_start:
        continue
    for px in range(py): # 2-gons are the base case.
        if gcd(px,py)==1:
            min_gon[px,py]={2:[((px,py),0)]}
        else:
            min_gon[px,py]={}
            
    work = num_work = 0 # statistics
    for px,l2 in enumerate(collect[py]):
      for k,triples in enumerate(l2):
          if not triples:
              continue
          work += len(triples)
          num_work += 1 # statistics
          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 optimal solutions:
                      # Then multiple vol-entries should be kept, but SORTED BY SLOPE.
                  prev_fx,prev_fy,prev_vol = fx,fy,vol
              
          min_gon[px,py][k+1]=result
          # print (f"{(px,py)=} {k=} {triples=}"); print(result)
          
    # put together solutions
    # update target_height
    # compute new record areas
    improved = [False]*(N_target+1)

    def check_record(k, vol, how, text=""):        
        if k<k_START or k>N_target:
            return
        if vol<record_vol[k]:
            #if text: print(text)
            if confirmed[k]:
                print(k,vol,how,text)
            assert not confirmed[k]

            record_vol[k] = vol
            how_achieved[k] = how
            min_height[k] = py
            improved[k] = True
        if vol==record_vol[k] and how[0] != "DIAGONAL":
            # DIAGONAL case not counted: Without a flat edge, height is unbounded.
            max_height[k] = py

    # Combine "right" k1-gons with "left" k2-gons
    for px in range(py):
        if EVEN_N:
            # even case is easy: insert a horizontal edge at the top: 
            ppx0 = (-px)%py  # unnecessarily complicated:
            # One could just turn the solution by 180 degrees.
            # It is possible that this produces an asymmetric solution.
            for k1,best1 in min_gon[px,py].items():
                _,vol1 = best1[0]
                check_record(k1, # a (2k)-gon
                             vol1+py, # actual area, not "volume"
                             ("2 FLAT",px,py,k1,k1))
        else:
          ppx0 = (-px)%py  # A) exact match along long edge (0,0)-(px,py)
          ppx1 = (1-px)%py # B) horizontal edge of length 1 at the top
          for k1,best1 in min_gon[px,py].items():
            _,vol1 = best1[0]
            for k2,best2 in min_gon[ppx0,py].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))
            for k2,best2 in min_gon[ppx1,py].items():
                _,vol2 = best2[0]
                check_record(k1+k2-1, vol1+vol2+py, ("1 FLAT",px,py,k1,k2))

    target_height = 0 # how high we still have to go
    for k in range(k_START,N_target+1):
        if not confirmed[k]:
            num_vertices = 2*k if EVEN_N else k
            vol_factor = 2 if EVEN_N else 1
            if lower_bound(num_vertices, py+1)>=record_vol[k]*vol_factor:
                # evaluated for increased py in the next iteration
                #print("LB",k,py+1,LB)
                confirmed[k] = True
                print_solution(k,how_achieved[k],record_vol[k])
            else:
                # find necessary target height to confirm the current record:
                t_k = py+2
                while lower_bound(num_vertices, t_k)<record_vol[k]*vol_factor:  t_k += 1
                target_height = max(target_height,t_k-1)
        if confirmed[k] and k==output_so_far+1:
            print_solution(k,how_achieved[k],record_vol[k], True) # print to a_file and b_file
            output_so_far += 1
                

    print("N=%d, height %d"%(N_target,py) +
        (".\nRESULT:" if target_height==0 else "->%d:"%target_height),
        ",".join(("*" if improved[k] else "")+
                 str(record_vol[k])+
                 ("" if confirmed[k] else "?")
            for k in range(k_START,N_target+1)),
        #
        # Meaning of the signs:
        # ?  still unconfirmed
        # *  has just been improved
        #
        "work=%d"%work, # work done in this step, total items in "collect"
        "" if num_work==0 else "avg=%.2f"%(work/num_work), # average list size in "collect"
        flush=True)
    if py%10==0:
      print("smallest heights:",*("%2d"%x for x in min_height[k_START:]))
      print("largest heights: ",*("%2d"%x for x in max_height[k_START:]))
    if py%10==0:
        # average list length rises very slowly:
        # 3.5 for height py=70, 3.9 for height py=100,
        # max length, on the other hand is 71 for py=100.
        # the long lists are mostly associated with steep vectors p,
        # and they contain many f-vectors of the form (-1,..) or (1,..)
        print("avg list length",
        sum(len(best1)
            for ppy in range (1,py+1)
            for ppx in range (ppy)
            for k1,best1 in min_gon[ppx,ppy].items()) /
        sum(1
            for ppy in range (1,py+1)
            for ppx in range (ppy)
            for k1,best1 in min_gon[ppx,ppy].items()))
        print("max list length",
        max(len(best1)
            for ppy in range (1,py+1)
            for ppx in range (ppy)
            for k1,best1 in min_gon[ppx,ppy].items()))

    if all(confirmed):
        break

  py_target_start = py_target_end+1
  py_target_end = py_target_end * 2 # make smaller steps if memory space is an issue.
  py_target_end = min(py_target_end,target_height)

a_file.close()
b_file.close()
                    
    
print("k:               ",*("%2d"%x for x in range(k_START,N_target+1)))
print("smallest heights:",*("%2d"%x for x in min_height[k_START:]))
print("largest heights: ",*("%2d"%x for x in max_height[k_START:]))
# "largest heights" are preliminary, since the algorithm stops
# too early to find the largest heights.
# "smallest heights" is definitely the smallest lattice width that
# a convex lattice n-gon of minimum area can have.

# There is quite a large variation of heights, in line with the
# result of Bárány and Tukushige that, in the limit,
# the shapes resemble ellipses that are more and more oblong.
# But already a single n-gon has necessarily a large variation of
# lattice widths, when considering all edge directions.
# This is because the lattice width in direction (a,b) with gcd(a,b)=1
# is the Euclidean width times the length of (a,b).
# Since the lengths of edge vectors vary a lot (a vector (1,0)
# may be adjacent to a vector (L,1) for large L, while the
# Euclidean width will hardly change when turning to the next edge,
# a large fluctuation is inherent.
# The limiting behavior predicted by Bárány and Tukushige
# may set in much later than how far the calculations go.