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

# This program works for two different sequences:
Problem = "A322345" # number of interior points 
#Problem = "A298562" # "polygonal" Helly numbers

# "volume" always stands for half the area ("normalized area")

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

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

"""
Adaptation of a dynamic programming program for A070911

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

"""

lattice_width_factor  = 2/3 + 0.00001
# +0.00001 is because of the strict inequality (except for triangles.)

def lower_volume_bound(w):
    "lower bound on vol (=2*area) for polygon of lattice width w "
    LB = ceil( w**2 * lattice_width_factor )
    return LB

def min_width(vol):
    "upper bound on the lattice width for polygons of volume 'vol' "
    for w in count(1):
        if lower_volume_bound(w)>vol:
            return w-1
    

"""

Vol = 2*area
k = number of vertices
E = number of non-vertex points ("edge points") on the boundary
I = number of interior points

Pick's formula: Vol = 2I + k + E - 2 

Both problems can be investigated under the same framework by
and appropriate definition of a quantity that we call the "measure" M.
In both problems, M+2-k is specified, and we want to maximize k. 

A322345
=======
Maximal number of vertices of a convex lattice polygon containing n lattice
points in its interior. 

I=n is fixed, k should be maximized.

2n = 2I = Vol - E + 2 - k is fixed

Define "measure" M as Vol - E
(Warning: For n=0, this is unbounded.)

A298562:
=======
a(n) = g(Z^2,n) is the maximum integer k > 0 such that there exists
a lattice polygon with k vertices containing n+k lattice points.
(or in other words, n lattice points in the interior on on edges.)

I+E+k=n+k, n is given, k should be maximized.

2n = 2I+2E = Vol + E + 2 - k is fixed.

Define "measure" M as Vol + E

=======
In both problems: M = 2n+k-2

"""


A070911_value = {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,
    # ... more values are known
))}

# We can compute a rigorous upper bound "max_k" on max_result (=k)
# as follows. (for both problems)

# The smallest-area k-gon contains no points on the boundary except vertices.
# The smallest volumes V_k of k-gons are known, see OEIS A063984.
#
# If the smallest lattice k-gon has vol V_k = k+2I-2, this means that
# any k-gon of volume V_k or larger contains at least (V_k-k+2)/2 interior points.
# (and at least so many interior points plus points on edges).
#
#    A063984(k) = A070911(k)/2 - k/2 + 1
#               = min { n : A322345(n)>=k }
#
# In other words:
#
#    if n < A063984(k) = A070911(k)/2 - k/2 + 1
#         then A322345(n) < k and A298562(n) < k
    
for k in count(3):
    if A070911_value[k]-k+2 > 2*N_target:####!!! 2*
        # If there are ("not enough data from A070911"),
        # this will lead to an IndexError exception.
        max_k = k-1
        break

max_k += 5
# error margin for sloppy thinking, to be on the safe side.
maxVOL = A070911_value[max_k]
maxHEIGHT = min_width(maxVOL)

    
maxM = 2*N_target + max_k - 2 # M = 2n + k - 2
explanation = f"The program assumed that 2n+k is at most {maxM+2}."

if Problem == "A322345":
    Problemtext = "n interior points"
    def measure(vol,E):
        return vol-E
    # https://oeis.org/A322345 (n=0..30)
    known = [4,6,6,6,8,7,8,9,8,8,10,9,9,10,10,10,10,11,10,12,
             12,12,11,11,12,12,12,13,12,12,13]
    
    # We expect something like the inverse of A070911_value (plus a safety margin)    
    # estimated upper bound max_k on max_result (=k):
    rough_estimated_result = ceil((N_target*(1/3)*4)) + 5
    # not used, replaced by the rigorous derivation of max_k
        
elif Problem == "A298562":
    Problemtext = "Helly numbers, n points in interier or on edges"
    def measure(vol,E):
        return vol+E
    known = [4,6,6,6,8,7,8,9,8,8,10,9,9,10,10,10,10,11,11,12,
             12,12,11,11,12,12,12,13,12,12,13]

    
print(f"estimated={max_k} maxM={maxM} maxVOL={maxVOL} maxHEIGHT={maxHEIGHT}")
    
def yardstick(M,k): # the parameter "n"
    # M=measure
    # valid for both problems
    assert (M-k) % 2 == 0
    return (M+2-k)//2



"""
===================================================
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 inspired by 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 possible measures
of the counterclockwise k-gons O,....,q,p, under the following constraint:
* The successive points move strictly upwards.
  (Thus, p is the top point, and there is a "long edge" from O to p.)

More precisely, for each p, for each value M of the measure, and for each k,
we store the most clockwise point q (as seen from p) for which
a k-gon (O,....,q,p) with that measure exists.

Reason: If we have two polygons with different points q and q', where
q is clockwise from q', we can always replace the second one by the first,

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

In evaluating the measure, we don't take into account the boundary points
on the "long" edge Op, because that might become an interior edge when
gluing another edge to it.

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 (-r,0) for some r>0 and
ending in p or in p-(r',0).

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

min_gon[(px,py)][k,M] is the most clockwise incoming edge f=(fx,fy)
(pointing upwards, fy>0), for which a k-gon with measure M exists.

We fill this dictionary row by row, increasing py.

"vol" is TWICE the area.



POTENTIAL IMPROVEMENTS:

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

"""

print(Problem, Problemtext+".", "target =", N_target)

min_gon = dict()
record = dict()

how_achieved = dict()

# Auxiliary procedures to construct the solution, once the optimal
# area has been determined
def find_polygon(px,py,k,M,alpha0=0):
    # find the smallest k-gon ending in (px,py) by backtracing
    # 0<px<py. The two points (0,0),(1,0) are not produced.
    alphaT = alpha0
    result = [(px+alphaT*py,py)]
    for k in range (k,2,-1):
        (fx,fy) = min_gon[px,py][k,M]
        qx,qy = px-fx,py-fy
        M -= measure(qx*py-qy*px, gcd(fx,fy)-1)
        alpha =  qx//qy
        px,py = qx-alpha*qy,qy
        alphaT += alpha # alpha is accumulated
        result.append((px+alphaT*py,py))
    return result

def print_solution(n,how_achieved,k):
    if how_achieved[0]=="LONG":
        error("not implemented")
    h_top,h_bottom,px,py,k1,k2,M1,M2 = how_achieved
    #print(f"{(k,how,px,py,k1,k2)=}")
    p1 = find_polygon(px,py,k1,M1)
    p1.reverse()
    top_point = px-(h_top-h_bottom)
    ppx = (-top_point)%py
    alpha0 = (-top_point - ppx) // py # usually, alpha0==-1
    leftshift = h_bottom
    p2 = [(-x-leftshift,y) for x,y in find_polygon(ppx,py,k2,M2,alpha0)]
    #print(f"-- {how=} {alpha0=}, {p2[0][0]=},{top_point=},{ppx=},{px=},{py=},{k=}")
    assert p2[0][0]==top_point-leftshift
    if h_top==0:
        p2 = p2[1:] #eliminate the common point
    last = [(-h_bottom,0)] if h_bottom>0 else []
    polygon = [(0,0)]+p1+p2+last # counterclockwise
    out = "["+",".join(f"({x},{y})" for x,y in polygon)+"]" # compact format without spaces
    print ("an optimal polygon for n=%d,"%n, "k = %d:"%k, out)
    print(n,k,out, file=a_file)
    print(n,k, file=b_file)


## Start the computation

def clockwise(f,g):
    "is vector f strictly clockwise from g?"
    fx,fy=f
    gx,gy=g
    # assume fy,gy>0
    # return fx/fy > gx/gy
    return fx*gy > gx*fy

for py in range(1,maxHEIGHT+1): #itertools.count(1): # for py = 1,2,...
    # proceed row by row
    for px in range(py): # 2-gons are the base case.
        result_list = min_gon[px,py]={(2,measure(0,gcd(px,py)-1)):(px,py)}

        work = num_work = 0 # statistics
        sum_alpha = num_alpha = 0
        for qy in range(1,py):
            # We look at each starting point (qx0,qy) with 0<=qx0<qy<py
            # Then we consider the series of affine images
            #    (qx0,qy), (qx0+qy,qy), (qx0+2*qy,qy), ..., (qx0+alpha*qy,qy), ...
            # until we are sure that there are no more contributions to row py
            # to be expected.
            #
            Delta_y = py-qy # >0
            for qx0 in range(qy):
                for alpha in itertools.count(0): # for alpha = 0,1,2,...
                    qx = qx0+alpha*qy # apply shearing by alpha
                    vol = qx*py-qy*px
                    if vol<=0:
                        continue
                    if vol>maxVOL:
                        break
                    Delta_x = px-qx
                    E = gcd(Delta_x,Delta_y)-1
                    for (k,M),(fx0,fy) in min_gon[qx0,qy].items():
                        fx = fx0+alpha*fy # apply the shearing by alpha
                        if k>=N_target-1:
                            continue
                        if not clockwise((fx,fy),(Delta_x,Delta_y)):
                            continue
                        M_new = M + measure(vol,E)
                        if M_new> maxM:
                            continue
                        if (k+1,M_new) not in result_list or (
                          clockwise((Delta_x,Delta_y), result_list[k+1,M_new])):
                            result_list[k+1,M_new] = (Delta_x,Delta_y)
            if 0:
                #print(f"p=({px},{py}) {k=}")
                for t in result: print("  ",t)

    # compute new record areas
    improved = [False]*(N_target+1)

    def check_record(k, M, how, text=""):
        if k==2:
            return
        n = yardstick(M,k)
        if n>N_target:
            return
        if n not in record or k>record[n]:
            record[n] = k
            how_achieved[n] = how
            #improved[n] = True

    # Combine "right" k1-gons with "left" k2-gons
    if py==1:
        continue # arbitrarily large no interior lattice points
    for px in range(py):
        #0) consider whole polygon with "long" edge Op.
        for k,M in min_gon[px,py].keys():
            M_new = M + measure(0, gcd(px,py)-1)
            check_record(k, M_new, ("LONG",px,py,k,M), text="")
        ppx0 = (-px)%py  # A) exact match along long edge (0,0)-(px,py)
        for k1,M1 in min_gon[px,py].keys():
            for k2,M2 in min_gon[ppx0,py].keys():
                M = M1+M2
                check_record(k1+k2-2,
                             M, (0,0,px,py,k1,k2,M1,M2))
                for hori in count(1):
                    MM = M + measure(2*hori*py,2*(hori-1))
                    # 2 flat edges of length hori
                    if MM>maxM: break
                    check_record(k1+k2,
                             MM, (hori,hori,px,py,k1,k2,M1,M2))
                 #f"DIAG k={k1+k2-3} {k1=} {k2=} {vol1=} {vol2=} p=({px},{py})")
            for delta_h in count(1):
                # B) horizontal edge of length at the top is
                # delta_h longer than at the bottom
                M1_plus = M1 + measure(delta_h*py, delta_h-1)
                if M1_plus > maxM: break
                ppx1 = (delta_h-px)%py
                for k2,M2 in min_gon[ppx1,py].keys():
                    M = M1_plus+M2
                    # B1) no edge at the bottom:
                    #print(Problem,f"{py=} {k1=} {M1=} {k2=} {M2=}  {M1_plus=} ")
                    check_record(k1+k2-1,
                             M, (delta_h,0,px,py,k1,k2,M1,M2))
                    for hori in count(1):
                        # B2) flat edges of length hori and hori+delta_h
                        MM = M + measure(2*hori*py,2*hori-1)
                           # hori at the top edge and hori-1 at the bottom
                        if MM>maxM: break
                        check_record(k1+k2,
                             MM, (delta_h+hori,hori,px,py,k1,k2,M1,M2))

    print(py,end=":")
    ks = set(record.keys())
    if max(record.values())<= max_k:
        ok = "OK"
    else:
        ok = "still bigger than expected"
    
    if ks == set(range(max(ks)+1)):
        records = list(record[k] for k in sorted(ks))
        print(",".join(str(r) for r in records),
              all((x==y) for x,y in zip(records,known)), ok, flush=True)
    else:
        print( list(sorted(record.items())))
                        
print("Sequence",Problem)

b_file = open(f"b{Problem[1:]}.txt","w")
a_file = open(f"a{Problem[1:]}.txt","w")
for n in sorted(ks):
    k=record[n]
    #print (f"{n=} {k=}")
    print_solution(n, how_achieved[n],k)
a_file.close()
b_file.close()

    
OK = max(record.values()) <= max_k
VERYGOOD = max(record.values()) <= max_k*2/3
print(f"The largest found value, {max(record.values())},",
      "is well below" if VERYGOOD else
      "did not exceed" if OK else "exceeded",
      f"the estimate {max_k}.")
if OK:
    print("The results should be safe.")
else:
    print("Run the program with a larger estimate.")
print(explanation)