# Generate ALL solutions

# More precisely, we classify two polygons as equivalent if there is
# a lattice preserving affine transformation that maps one
# polygon to the other, and we output one sample from each equivalence class.

"""
This is a modification of the program for A089187

In that program, because of the normalization that was used to limit the height (py), rotated and
reflected copies of the same solution may or may not be generated,
but it is guaranteed that at least one solution from every
equivalence class is produced.
Thus, in the current program, whenever an n-gon is found, we produce all 2n copies where
two consecutive edges (in either order) are mapped to (0,1),(0,0),(1,0),
and we select a representative copy (the lex-smallest one).

For n<=100, it turned out that the solution is unique most of the times, except
for n=7,11,15,20,22,45, when there are 2 different solutions. More than 2 solutions have never been found.

We also looked for counterexample to the property of Lemma 5B, which is not proved in the odd case.
The counterexamples that we found were the 5-gon, the 7-gon (both optimal solutions), and the 11-gon
(both optimal solutions).

We also computed the smallest and largest "height" when any of the polygon edges is mapped to
the segment (0,0),(1,0). (The lattice width in an edge direction).
Here are the results found:

k:                 3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
#different k-gons  1  1  1  1  2  1  1  1  2  1  1  1  2  1  1  1  1  2  1  2  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  2  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
smallest heights:  1  1  2  2  3  3  4  4  5  6  6  7  8  9 10 10 12 11 15 13 14 18 16 19 22 20 23 23 26 27 30 28 27 29 30 33 33 37 37 40 38 44 40 43 44 46 47 50 54 48 60 52 64 65 69 69 70 74 74 75 80 76 81 81 86 87 92 88 86 90 94 95 90 101 100 99 105 105 110 110 114 116 116 120 123 125 127 126 132 127 128 129 128 131 136 134 134 137
largest heights:   1  1  3  2  6  4  7  7 15 10 16 16 23 19 31 28 36 34 39 43 57 46 66 64 77 71 92 85 100 91 107 109 133 117 159 141 170 149 186 183 196 190 211 205 234 218 254 250 283 265 289 294 304 298 313 325 383 334 442 408 439 419 474 463 498 485 509 504 537 520 578 556 634 583 630 614 659 653 721 677 735 706 812 766 885 775 932 876 994 890 999 923 1017 951 1087 1016 1134 1058

One can observer quite a large variation of heights between the minimum and the maximum,
when considering the different edge direction.
Some of this variation is inherent,
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 necessarily vary (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 cannot be avoided.

The result of Bárány and Tukushige (see below) implies that, in the limit,
the shapes resemble ellipses that are more and more oblong.
The observed variation may be a sign of this effect,
or limiting behavior predicted by Bárány and Tukushige
may set in much later than how far the calculations go.
"""

N_target = 44

# A different target value of N can 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 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 seems that it might hold also
in the odd case, but I cannot prove it.

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
    ))} # for guidance;
#known_values = {}    

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

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

    
    
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,


num_kgons = [0 for n in range(N_target+1)]
Solutionlist = [set() for n in range(N_target+1)]
Exceptionlist = []

k_START = 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)


def check_Lemma5(p):
    for i,((x0,y0),(x1,y1),(x2,y2)) in enumerate(zip( p, p[1:]+p[:1], p[2:]+p[:2])):
        u1,v1 = x1-x0,y1-y0
        u2,v2 = x2-x0,y2-y0
        det = u1*v2-u2*v1
        #print(x0,y0,det,p)
        assert det>0
        if det>1:
            if tuple(p) not in Exceptionlist:
                Exceptionlist.append(tuple(p))
            print(f"Exception {len(p)}-gon {p} start at p[{i}]={(x0,y0)} det={det}")

def normalize(p):
    """The lex-smallest polygon in the first quadrant that
    has first edge (0,0) (1,0) and last vertex (0,1).

    update largest height as a side result.
    """
    k = len(p)
    max_height[k] = max(max_height[k],
                        max (y for polygon in images(p)
                             for x,y in polygon))
    return min(images(p))

def images(p):
    found = [False]*len(p)
    for i in range(len(p)):
        r = p[i:]+p[:i] # cyclic shift
        x0,y0 = r[0]
        s = list((x-x0, y-y0) for x,y in r)
        u1,v1 = s[1]  # 2nd point
        u2,v2 = s[-1] # last point
        det = u1*v2-u2*v1
#        print(p,i,s)
        assert det>0
        if det>1:
            continue # consider only triplets with det=1
        found[i] = True
        # map V1,V2 to (1,0),(0,1)
        #          ( v2  -u2 )
        # matrix = (         )
        #          ( -v1  u1 )
        yield tuple( (v2*x-u2*y, -v1*x+u1*y) for x,y in s)
        # map -V2,-V1 to (1,0),(0,1) # reversed order
        yield ((0,0),)+tuple((-v1*x+u1*y, v2*x-u2*y)
                     for x,y in reversed(s[1:]))
    if not all(found[i] or found[i+1] for i in range(-1,len(p)-1)):
        if len(p) not in (7,5):
            print("largest height not correct",p,found)
            error()

# Auxiliary procedures to construct the solution, once the optimal
# area has been determined
def find_all_polygons(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.
    (fx,fy),vol = min_gon[px,py][k][0]
    prev_fx, prev_fy = -1,0 # ensure that test is true at the first iteration
    return find_all_polygons_recursive(px,py,prev_fx,prev_fy,vol,k,alpha0,[])

def find_all_polygons_recursive(px,py,prev_fx,prev_fy,vol,k,alpha,result):
    result = result + [(px+alpha*py,py)] # don't append
    if k==2:
        #print("***",result)
        yield result # + [(0,0)]
        return
    
    #print(f" min_gon[{px},{py}][{k}] = {min_gon[px,py][k]} {vol=}"
    #                       f" {(prev_fx,prev_fy)=}")
    for (fx,fy),vol2 in min_gon[px,py][k]:
        if vol2>vol:
                break
        elif vol2==vol and (
                # (fx,fy) is clockwise from prev
                # fx/fy > prev_fx/prev_fy
        fx*prev_fy > prev_fx*fy):
            qx,qy = px-fx,py-fy
            new_vol = vol - (qx*py-qy*px)
            D_alpha =  qx//qy
            ffx = fx-D_alpha*fy
            qqx = qx-D_alpha*qy
            for r in (
           find_all_polygons_recursive(qqx,qy,ffx,fy,new_vol,k-1,alpha+D_alpha,result)):
                yield r

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_all_solutions(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)=}")
    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
    #print(f"STARTING *  {(px,py)=} {k1=}")
    for p1 in find_all_polygons(px,py,k1):
        #print(f"{(px,py)=} {k1=} {p1=}")
        p1 = list(reversed(p1))
        for p20 in find_all_polygons(ppx,py,k2,alpha0):
            p2 = [(-x-leftshift,y) for x,y in p20]
            #print(f"-- {how=} {alpha0=}, {p2[0][0]=},{top_point=},{ppx=},{px=},{py=},{k=}")
            print(p1,p2)
            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
            num_kgons[k] +=1
            Solutionlist[k].add(normalize(polygon))
            print (f"--- a {k}-gon,", f"area = {vol}/2:", polygon)
            check_Lemma5(polygon)
            
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 = list(reversed(p1))
    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 {k}-gon,", f"area = {vol}/2:", polygon)

## 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 True: # not EVEN_N:
                    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((vol,Fraction(fx,fy),fx,fy)
                                             for vol,fx,fy in triples)):
              # ordered clockwise.
              # For each slope, remember only the smallest volume,
              # and only if it is not smaller than all previous volumes.
              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):
                  result.append(((fx,fy),vol))
                  prev_fx,prev_fy = fx,fy
          for vol,fx,fy in []: # OLD VERSION (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

          # print(f"COLLECT {(py,py)=} {k=} {triples}")
          # print(result)
                  
          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

        if vol==known_values[k]:
            if (how[0] != "DIAGONAL" and (k not in (3,4) or py==1)):
                #print(f"NOW PRINTING {(k,vol,how)=}")
                print_all_solutions(k,how,vol)

    # Combine "right" k1-gons with "left" k2-gons
    for px in range(py):
          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 = k
            vol_factor = 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)

    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
  if py_target_end > 400:
      py_target_end += 25 # make smaller steps if memory space is an issue.
  else:
      py_target_end = min (35 + py_target_end, py_target_end * 2)
  # py_target_end = py_target_end * 2
  py_target_end = min(py_target_end,target_height)

Exceptionlist = []
a_file = open("a070911all.txt","w")    
for k,l in enumerate(Solutionlist):
    # print(f"{k=}")
    for polygon in l:
        print(k, ",".join(f"({x},{y})" for x,y in polygon))
        check_Lemma5(polygon)
        print(k,record_vol[k], ",".join(f"({x},{y})" for x,y in polygon),
              file=a_file) # compact format without spaces
a_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:]))

print("# k-gons found:  ",*("%2d"%x for x in num_kgons[k_START:]))
print("#different k-gons",*("%2d"%len(x) for x in Solutionlist[k_START:]))

print ("Exceptions to Property 5B")
for polygon in Exceptionlist.copy():
    check_Lemma5(polygon)