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