This is the second post to improve a recent result and prove that 4.5058(?)≤s(17), but let’s first define s(17). Quoting the old article about the topic
The high level and historic details are in the first post. It also has a few nice drawings of the weights that appear in the previous and new proof. Let’s go now to the low level details.
Sam Burns' bound
(In case you didn’t read the first post, let’s repeat this part that is relevant.)
The idea to prove 4.4811(?)≤s(17) posted by Sam Burns using ChatGPT picks 268 somewhat interesting points in a square of side 4.4811 The points have different weights, and the total weight is only 16.9476. After some reductions, it’s only necessary to test 181 directions using “almost-unit” (actually .9973) squares and verify that the sum of weight inside each one of them is at least 1 (actually 1.0003). So if we try to fit 17 unit squares there, at least two unit squares must share at least one of the 268 interesting points.
To test this, they use a program in Python. This method has false negatives. If it verifies a solution then it’s surely correct, but if the program fails there is a tiny chance that it’s a mistake. This is fine to ensure the weight proves a lower bound.
It uses a 29x29 grid with an empty margin of 0.5000 (hardcoded) and the total size of the grid is 3.4811.
Changes for the new bound
My idea was to try different combinations of the margin and internal grid size, so the first trivial step was to add a new parameter M for the double of the length of the margin, so M=1.0000 in the example of Sam Burns. .
As I said, it’s not clear how the weights were selected in the example of Sam Burns. So for each fixed size, I decided to use linear programming to find them using linprog in scipy.
(Like 30 years ago someone told me that linear programming can solve more problems than expected. It took me like 20 years to understand that advice.)
Each possible weight distribution can be made symmetric, this reduces a lot of the search space, so we only keep “1/8" of the grid (0≤x≤y≤14) as in Sam Burns´ certificate, but now each one of them get a variable for the linear programming method. The variables are copied to the other “8” parts of the grid using the symmetry, and we want to minimize the sum of these variables counted with multiplicity.
If we get less than 17, we get a distribution of weight for the new lower bound of s(17). In general, a value that is less than n is useful for other s(n) that may be useful in the future.
We only need to discover the set of linear inequalities. The sum in each square must be at least 1, so if for example, if v0, v1, and v2 are in a possible almost unit square and v0, v3, and 2 copies of v7 are inside another, and v1 and v2 are in a third one then we have a system like
v0 + v1 + v2 ≥ 1
v0 + v3 + 2 v7 ≥ 1
v1 + v2 ≥ 1
In the post by Sam Burns they have explicit values of the weight. In each direction, after rotating the grid if for every point with rotated coordinates (u, v) they place four charges in (u±.9973/2, v±.9973/2) with alternating signs. Then they calculate the cumulative sums in the horizontal and vertical directions to get the sum of the weights that are inside a square of side .9973. As a final part, they select the smaller inside the rotated big square.
Now we don’t have the exact value of each weight, so in the first step we add a tuple with the index of the variable and the sign. Then we accumulate them in a dict and then we accumulate the results in another dict. As a final step, we make a set with all those dicts that are inside the big square. For the inclusion in the set in the last step the “dicts” must be immutable. I should have used a frozendict, but they will be available next month, so for now I’m using a tuple after sorting the keys. In conclusion, the result of each direction is a set, that removes the duplication of the tuples/frozendicts automatically.
I collect all the 181 sets in an even bigger set, that removes even more duplication, but transforming it into an array to use linprog uses too much memory. So I remove the “obvious” inclusions, that is when a tuple/frozendict can be obtained from another removing exactly one count of one variable. So in the example we can drop the first inequality and get
v0 + v3 + 2 v7 ≥ 1
v1 + v2 ≥ 1
This does not remove all inclusions, but it’s a very common pattern in this problem, it’s easy, fast enough and enough to solve the memory problem.
(Actually, linprog requires
- v0 - v3 - 2 v7 ≥ -1
- v1 - v2 ≥ -1
so expect to see a few extra minus here and there in the conversion.)
As the final result, the program shows a rounded version of the weight (using ceil). I think round is fine too but rounding up ensures the inequalities are not broken by rounding. Just in case it shows the new total after ceiling the weights.
Sorry that my Python is not so good but the code to calculate the weight and get the certificate is at the bottom. Beware I’m using L=4.5000 as the side of the square and 0.775 as the margin so M=1.5500 instead of the optimal values. (For some reason, this saves memory. probably not so long integers inside the fractions.)
It takes like an hour to run, so go and grab a coffee. There are probably many optimizations to improve that. IWIMM.
Searching
How I got the final side and margin is more complicated…
- I was feeling lucky, so I fixed L=4.5000 that is bigger that the previous bound 4.4811
- I wrapped the program in a big for, and tried with M from 0.0000 to 1.9500 (inclusive) using 0.0500 steps.
- Only M=1.5500 got a result that was less than 17!
- I looked at the weight, they have too many ??333 and ??5 so I multiplied them by 2 and 3 until I got almost integer numbers, and I rounded them to drop noise.
- So I have to multiply by 576, and those are the numbers reported in the previous post.
- The original program that uses explicit weight with numbers is much faster, so I used it to add some digits. The solution is good if the weights don't cross an imaginary almost-unit square boundary, so a small increase should not break it sometimes.
- I’m too lazy for binary search, so I used decimal search to make it more fun. It’s possible to run the 200+ case to add a digit in a few minutes, but I actually reduce the search using that the results can be compared if both the length of the margin and the length of the grid are bigger or smaller.
- After repeating it a few times, I got L=4.45058 and M=1.5513 that are the values reported in the previous post that includes the program to verify them.
Conclusion and Future Work
- No idea if my code is optimal, probably not. In particular, using frozendict may be nice.
- For an exploratory phase I guess the fractions can be replaced by floating point numbers but the program needs a lot of modifications.
- Changing the size of the grid is trivial, and I expect some interesting better cases there. My version takes like one hour to run so I’d not try it but I expect it to improve the precision and show if the optimal weights are punctual or a continuous distribution.
- Changing the number of the 181 directions is possible, but it needs some adjustment of the constants here and there. It’s left as an exercise for the reader. Anyway, this is necessary to get more digits and reduce the number of false negatives
- I’d like to find the non-symmetrical version. I have some ideas to try, so check again in a few days. A non-symmetrical hopefully has like 1/8 of the weight and hopefully shows the almost equilateral triangles and is easier to understand without a computer.
Program to find the weights
from __future__ import annotations
from bisect import bisect_left, bisect_right
from fractions import Fraction as F
from math import ceil
import numpy as np
from scipy.sparse import coo_array
from scipy.optimize import linprog
# Original version posted by Sam Burns 2026
# Modified by Gustavo Massaccesi 2026
# Try to find a lower-bound certificate for packing 17 unit squares in a square.
#
# All geometric quantities and predicates are rational. NumPy is used only for
# integer range-addition and cumulative sums; no floating-point geometry is used.
L = F(45000, 10000) # side of the square
M = F(15500, 10000) # both empty borders
B = F(9973, 10000)
T = F(207107, 500000)
KMAX = 180
D = T / KMAX
NGRID = 29
LAST = NGRID - 1
# (i, j, w): every distinct D4 image of grid point (i,j). Use # of variable instead of weight
def build_cert_vars() -> list[tuple[F, F, int]]:
by_index: dict[tuple[int, int], int] = {}
v = 0
for i in range((NGRID + 1)//2):
for j in range(i, (NGRID + 1)//2):
by_index[(i, j)] = v
v += 1
return [(i, j, w) for (i, j), w in sorted(by_index.items())]
v_CERT = build_cert_vars()
print(v_CERT)
def orbit(i: int, j: int) -> set[tuple[int, int]]:
n = LAST
return {
(i, j), (n - i, j), (i, n - j), (n - i, n - j),
(j, i), (n - j, i), (j, n - i), (n - j, n - i),
}
def build_atoms() -> list[tuple[F, F, int]]:
step = (L - M) / LAST
coord = [M / 2 + step * i for i in range(NGRID)]
by_index: dict[tuple[int, int], int] = {}
for i, j, w in v_CERT:
for ij in orbit(i, j):
if ij in by_index:
raise ValueError(f"duplicate orbit assignment at {ij}")
by_index[ij] = w
return [
(coord[i], coord[j], w)
for (i, j), w in sorted(by_index.items())
]
# Clip a convex rational polygon against U >= bound or U <= bound.
def clip_u(
poly: list[tuple[F, F]],
bound: F,
keep_ge: bool,
) -> list[tuple[F, F]]:
if not poly:
return []
out: list[tuple[F, F]] = []
def inside(p: tuple[F, F]) -> bool:
return p[0] >= bound if keep_ge else p[0] <= bound
prev = poly[-1]
prev_in = inside(prev)
for cur in poly:
cur_in = inside(cur)
if cur_in != prev_in:
u1, v1 = prev
u2, v2 = cur
if u2 == u1:
v = v1
else:
lam = (bound - u1) / (u2 - u1)
v = v1 + lam * (v2 - v1)
out.append((bound, v))
if cur_in:
out.append(cur)
prev, prev_in = cur, cur_in
return out
def center_domain(c: F, s: F) -> list[tuple[F, F]]:
# A B-square at orientation (c,s) lies in [0,L]^2 exactly when its
# center lies in [h,L-h]^2, with h=B(c+s)/2.
# Transform that square to the B-square's (U,V) frame.
h = B * (c + s) / 2
lo, hi = h, L - h
corners_xy = [(lo, lo), (hi, lo), (hi, hi), (lo, hi)]
return [(c * x + s * y, -s * x + c * y) for x, y in corners_xy]
def verify_orientation(
c: F,
s: F,
atoms: list[tuple[F, F, int]],
) -> set:
"""Return the all the combination of variables for one rational orientation."""
half = B / 2
dom = center_domain(c, s)
u_dom_min = min(u for u, _ in dom)
u_dom_max = max(u for u, _ in dom)
v_dom_min = min(v for _, v in dom)
v_dom_max = max(v for _, v in dom)
rects: list[tuple[F, F, F, F, int]] = []
u_events = {u_dom_min, u_dom_max}
v_events = {v_dom_min, v_dom_max}
# In center coordinates, atom membership is an axis-aligned rectangle.
for x, y, w in atoms:
pu = c * x + s * y
pv = -s * x + c * y
u1, u2 = pu - half, pu + half
v1, v2 = pv - half, pv + half
rects.append((u1, u2, v1, v2, w))
u_events.add(u1)
u_events.add(u2)
v_events.add(v1)
v_events.add(v2)
ue = sorted(u_events)
ve = sorted(v_events)
ui = {x: i for i, x in enumerate(ue)}
vi = {x: i for i, x in enumerate(ve)}
# Exact integer 2D difference array. Scores are constant in every open
# event cell. NumPy performs only integer arithmetic here.
diff = np.empty((len(ue), len(ve)), dtype=object)
for u1, u2, v1, v2, w in rects:
a, b = ui[u1], ui[u2]
p, q = vi[v1], vi[v2]
assert(diff[a, p] is None)
assert(diff[b, q] is None)
diff[a, p] = (w,1)
diff[b, p] = (w, -1)
diff[a, q] = (w, -1)
diff[b, q] = (w, 1)
#scores = diff.cumsum(axis=0).cumsum(axis=1)
# We almost can use a set here, but we need a dict for variables near the center
diffcum0 = np.empty((len(ue), len(ve)), dtype=object)
for v in range(len(ve)):
s = {}
for u in range(len(ue)):
if diff[u,v] is not None:
s = s.copy()
w, sg = diff[u,v]
if w in s:
t = s[w] + sg
if t == 0:
s.pop(w)
else:
s[w] = t
else:
s[w] = sg
diffcum0[u,v] = s
scores = np.empty((len(ue), len(ve)), dtype=object)
for u in range(len(ue)):
s = {}
# This should have been an empty frozendict
s_t = ()
for v in range(len(ve)):
if diffcum0[u,v] is not None and diffcum0[u,v]:
s = s.copy()
for w, sg in diffcum0[u,v].items():
if w in s:
t = s[w] + sg
if t == 0:
s.pop(w)
else:
s[w] = t
else:
s[w] = sg
# This should have been a frozendict
s_t = tuple(sorted(s.items()))
scores[u,v] = s_t
nu, nv = len(ue) - 1, len(ve) - 1
all = set()
for i in range(nu):
u0, u1 = ue[i], ue[i + 1]
if u1 <= u_dom_min or u0 >= u_dom_max:
continue
slab = clip_u(dom, u0, True)
slab = clip_u(slab, u1, False)
if not slab:
continue
vlo = min(v for _, v in slab)
vhi = max(v for _, v in slab)
if vhi <= vlo:
continue
# This may examine a superset of feasible event cells, which is
# conservative for a lower-bound verification.
j0 = max(0, bisect_right(ve, vlo) - 1)
j1 = min(nv - 1, bisect_left(ve, vhi) - 1)
if j0 <= j1:
for j in range(j0, j1):
all.add(scores[i, j])
return all
def angle_net() -> list[tuple[F, F]]:
out: list[tuple[F, F]] = []
for k in range(KMAX + 1):
t = T * k / KMAX
den = 1 + t * t
c = (1 - t * t) / den
s = 2 * t / den
assert c * c + s * s == 1
out.append((c, s))
# The final adjacent pair brackets pi/4.
assert out[-2][1] < out[-2][0]
assert out[-1][1] >= out[-1][0]
# If psi_k=2 arctan(t_k), half an adjacent angular gap is
# arctan(t_{k+1})-arctan(t_k), whose tangent is
# D/(1+t_k*t_{k+1}) <= D. Therefore every angle in [0,pi/4]
# is within an error epsilon < D of a net direction.
for k in range(KMAX):
t0 = T * k / KMAX
t1 = T * (k + 1) / KMAX
tan_half_gap = (t1 - t0) / (1 + t0 * t1)
assert tan_half_gap <= D
return out
def my_ceil(x):
return ceil(10**10 * x) / 10**10
def main() -> None:
print(f"L = {L} = {L:.4f}")
print(f"M = {M} = {M:.4f}")
atoms = build_atoms()
print(f"atoms = {len(atoms)}")
net = angle_net()
# For an orientation error epsilon <= D,
# cos(epsilon)+sin(epsilon) <= 1+epsilon <= 1+D.
contain = B * (1 + D)
print(f"angle_net_size = {len(net)}")
print(f"b*(1+d) = {contain} = {float(contain):.12f} < 1")
assert contain < 1
global_min = set()
for k, (c, s) in enumerate(net):
m = verify_orientation(c, s, atoms)
global_min.update(m)
if k % 1 == 0 or k == KMAX:
print(
f"orientation {k:3d}/{KMAX}: "
f"size={len(m)}, "
f"global={len(global_min)}"
)
print(f"size = {len(global_min)}")
global_min_red = set()
for t in global_min:
for w, r in t:
s = dict(t)
if r>1:
s[w]=r-1
else:
s.pop(w)
s_t = tuple(sorted(s.items()))
if s_t in global_min:
break
else:
global_min_red.add(t)
print(f"size (reduced) = {len(global_min_red)}")
c = np.zeros(len(v_CERT))
for x,y,w in atoms:
c[w] += 1
b_ub = -np.ones(len(global_min_red))
A_rows = []
A_cols = []
A_data = []
for k, t in enumerate(global_min_red):
for w, r in t:
A_rows.append(k)
A_cols.append(w)
A_data.append(-r)
A_ub = coo_array((A_data, (A_rows, A_cols)), shape=(len(global_min_red), len(v_CERT)))
result = linprog(c=c, A_ub=A_ub, b_ub=b_ub)
print(f"L = {L} = {L:.4f}")
print(f"M = {M} = {M:.4f}")
print(f"success = {result.success}")
if result.success:
print(f"min sum = {result.fun}")
print("values = ", result.x)
new_CERT = [(x, y, float(result.x[w])) for x,y,w in v_CERT if result.x[w] != 0.0]
#print(f"min sum = {c @ result.x}")
#print(f"new CERT = {new_CERT}")
new_r_CERT = [(x, y, float(my_ceil(w))) for x,y,w in new_CERT]
print(f"min sum (rounded) = {c @ list(map(my_ceil, result.x))}")
print(f"new CERT (rounded) = {new_r_CERT}")
if __name__ == "__main__":
main()
