Gustavo

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

The idea is to improve a recent result and prove that 4.5058(?)≤s(17) using these weights:

Trevor Green’s bound

The idea of the old proof (19+40*sqrt(2))/17≅4.4452…≤s(17) of Trevor Green is to pick 16 very interesting "unavoidable" points in a square of side 4.4452… and then he uses a lot of geometry to prove that any unit square must include at least one of them. So if we try to fit 17 unit squares there, at least two unit squares must share one of the 16 interesting points. The construction chooses 16 points out of a 4x6 grid.

I only found an image of the points in the old article, but I couldn't find the analytical definition. Looking at the formula for the side of the square, and using a rule, and some guessing, I think that the empty left/right margin is 0.5 and the empty top/bottom margin is sqrt(2)-1/2≅0.9142… With these choices, the diagonal segment in the original graphic has length 1, which is a very useful number to make triangles that have vertices that are unavoidable points. (I’d be glad to hear a confirmation.)

It uses a 6x4 grid with an empty margin of 0.9142… and 0.5000, and the total size of the grid is 2.6168… and 2.4452…

To compare the construction to the newer constructions, it’s better to symmetrize it. In this symmetrized version each unit square includes at least 4 points, but some points are thicker, and they count as double points (more details later).

Sam Burns’ bound

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 a finite number of directions and they use a program in Python to test “all” the possible “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 somewhat interesting points. (More details in the second post.) 

As I said, the weights are closer to the border than what I expected looking at the old examples, close to the second/penpenultimate row/column of the previous one. It also uses fewer weights, so I hope it’s easier to prove that it’s correct without a computer. I’d like to make a non symmetric version, that may be even better.

The program published by Sam Burns assumes that the empty margin is 0.5, so I had to modify it slightly to allow arbitrary borders with a variable M that is the double of the margin. The version with that modification, the new sizes and the new table of weight is at the bottom. Running that program and making the obvious changes to the explanation posted by Sam Burns proves(?) the new bound.

Conclusion and Future Work

• The distributions look quite discrete in the corner, but it has some strange bars near the center. It would be nice to increase the grid size and take a look. Also, the narrow empty margins appear to be useful.
• My search program in the second article is too slow (like 1 hour), so I avoided changing the size of the grid. It may be useful to explore other grid sizes in case there are some interesting coincidences. 
• Adding more digits takes only a few minutes, I didn't bother because it looks like refining the grid or using more directions for the rotations would make bigger changes.
• This result also automatically improves the lower bound of s(18), s(19) and s(20). But a more deep search for those values should provide even better bounds. I’ve seen too many cases where the total sum of the weight is 18. There is something interesting about 18.
• 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 version hopefully has like 1/8 of the weight and hopefully shows the almost equilateral triangles and is easier to understand without a computer. 

You may like to read the second post with details about how I got the new weights. 

Program to verify the bound

from __future__ import annotations

from bisect import bisect_left, bisect_right
from fractions import Fraction as F
import numpy as np

# Original version posted by Sam Burns 2026
# Modified by Gustavo Massaccesi 2026

# Proposed exact 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(45058, 10000)   # side of the square
M = F(15513, 10000)   # both empty borders
B = F(9973, 10000)
T = F(207107, 500000)
KMAX = 180
D = T / KMAX
WEIGHT_SCALE = 576    # min weight
NGRID = 29
LAST = NGRID - 1

# (i, j, w): every distinct D4 image of grid point (i,j) receives weight w/WEIGHT_SCALE.
CERT = [
    (0, 2, 165),
    (0, 11, 129),
    (1, 8, 36),
    (1, 10, 21),
    (1, 11, 15),
    (2, 2, 246),
    (2, 8, 129),
    (2, 9, 105),
    (2, 10, 36),
    (2, 11, 105),
    (5, 10, 36),
    (6, 10, 63),
    (6, 11, 12),
    (7, 10, 21),
    (8, 9, 33),
    (8, 11, 15),
    (9, 11, 75),
    (9, 14, 39),
    (10, 11, 25),
    (10, 12, 21),
    (10, 13, 24),
    (10, 14, 3),
    (11, 11, 16)
]


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 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]],
) -> int:
    """Return the exact minimum integer score 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.zeros((len(ue), len(ve)), dtype=np.int64)
    for u1, u2, v1, v2, w in rects:
        a, b = ui[u1], ui[u2]
        p, q = vi[v1], vi[v2]
        diff[a, p] += w
        diff[b, p] -= w
        diff[a, q] -= w
        diff[b, q] += w

    scores = diff.cumsum(axis=0).cumsum(axis=1)
    nu, nv = len(ue) - 1, len(ve) - 1

    best = 10**18
    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:
            row_min = int(scores[i, j0:j1 + 1].min())
            best = min(best, row_min)

    if best == 10**18:
        raise RuntimeError("center domain was not enumerated")
    return best


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 main() -> None:
    atoms = build_atoms()
    total = sum(w for _, _, w in atoms)

    print(f"atoms = {len(atoms)}")
    print(
        f"total_weight = {total}/{WEIGHT_SCALE}"
        f" = {total / WEIGHT_SCALE:.4f}"
    )
    #assert len(atoms) == 268
    #assert total == 169476
    assert total < 17 * WEIGHT_SCALE

    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 = 10**18
    argmin = -1
    for k, (c, s) in enumerate(net):
        m = verify_orientation(c, s, atoms)
        if m < global_min:
            global_min, argmin = m, k
        if k % 30 == 0 or k == KMAX:
            print(
                f"orientation {k:3d}/{KMAX}: "
                f"min={m}/{WEIGHT_SCALE}, "
                f"global={global_min}/{WEIGHT_SCALE}"
            )

    print(
        f"minimum_score = {global_min}/{WEIGHT_SCALE}"
        f" = {global_min / WEIGHT_SCALE:.4f} at k={argmin}"
    )
    assert global_min >= WEIGHT_SCALE

    print("CERTIFICATE CONDITIONS VERIFIED.")
    print(f"By the scaling argument: s(17) >= {L} = {L:.4f}.")

    
if __name__ == "__main__":
    main() 

Program to draw the images

#lang racket

(require racket/list)
(require metapict)


{define-syntax-rule (for/append clauses body ...)
  ; Todo: Add support for #:breack and #:final
  (append* (for/list clauses (begin body ...)))}

{define (mirror-x N atoms)
  (for/append ([a (in-list atoms)])
    (match a
      [(list x y w)
       (list (list x y w) (list (- N 1 x) y w))]))}

{define (mirror-y N atoms)
  (for/append ([a (in-list atoms)])
    (match a
      [(list x y w)
       (list (list x y w) (list x (- N 1 y) w))]))}

{define (mirror-d atoms)
  (for/append ([a (in-list atoms)])
    (match a
      [(list x y w)
       (list (list x y w) (list y x w))]))}

{define-values (L-Green atoms-Green)
  (let ()
    ; Todo: Confirm these are the correct lenghts.
    (define grid-nx 6)
    (define grid-ny 4)
    (define border-size-y (- (sqrt 2) 1/2))
    (define grid-size-y (/ (+ 12 (sqrt 8)) 17))
    (define border-size-x 1)
    (define L 
      (+ (* border-size-y 2) (* grid-size-y 3)))
    (define grid-size-x (/ (- L (* border-size-x 2)) 5))
    
    {define atoms/int '(#;()
                        (0 3 1) (1 3 1) (3 3 1) (5 3 1)
                        (0 2 1) (2 2 1) (4 2 1) (5 2 1)
                        (0 1 1) (1 1 1) (3 1 1) (5 1 1)
                        (0 0 1) (2 0 1) (4 0 1) (5 0 1))}
    (define min-weight 1)

    {define atoms (for/list ([a (in-list atoms/int)])
                    (match a
                      [(list x y w)
                       (list (+ border-size-x (* grid-size-x x))
                             (+ border-size-y (* grid-size-y y))
                             (/ w min-weight))]))}
    (values L atoms))}

{define-values (L-Green/S atoms-Green/S)
  (let ()
    ; Todo: Confirm these are the correct lenghts.
    (define grid-nx 6)
    (define grid-ny 4)
    (define border-size-y (- (sqrt 2) 1/2))
    (define grid-size-y (/ (+ 12 (sqrt 8)) 17))
    (define border-size-x 1)
    (define L 
      (+ (* border-size-y 2) (* grid-size-y (- grid-ny 1))))
    (define grid-size-x (/ (- L (* border-size-x 2)) (- grid-nx 1)))

    ; It's easier to calculate the overlaps by hand
    (define atoms/gen '(#;()
                        (0 1 2) (1 1 1) (2 1 1)
                        (0 0 2) (1 0 1) (2 0 1)))
    (define min-weight 4)

    (define atoms/int (remove-duplicates
                       (mirror-x grid-nx
                                 (mirror-y grid-ny
                                           atoms/gen))))
    (define atoms/one-dir (for/list ([a (in-list atoms/int)])
                            (match a
                              [(list x y w)
                               (list (+ border-size-x (* grid-size-x x))
                                     (+ border-size-y (* grid-size-y y))
                                     (/ w min-weight))])))
    (define atoms (mirror-d atoms/one-dir))
    (values L atoms))}

{define-values (L-Burns atoms-Burns)
  (let ()
    (define grid-n 29)
    (define L 44811/10000)
    (define M 1)
    (define grid-size (/ (- L M) grid-n))
    (define border-size (/ M 2))
    (define min-weight 10003)
    {define atoms/gen '(#;()
                        (1 11 107) (2 4 137) (2 9 214) (2 11 107) (2 12 137)
                        (3 4 3884) (3 7 214) (3 8 913) (3 9 214)
                        (3 10 214) (3 11 1234) (3 12 2189) (3 14 384)
                        (4 4 1961) (4 7 520) (4 8 214) (4 9 1413) (4 10 1234)
                        (4 11 1083) (4 13 137) (4 14 292)
                        (7 11 529) (7 12 33) (8 10 906) (8 11 384) (8 12 351)
                        (9 9 340) (9 10 180) (9 11 204) (9 12 549)
                        (10 12 879) (10 13 201) (10 14 378)
                        (11 11 396) (11 12 622) (11 13 204) (11 14 204))}

    (define atoms/int (remove-duplicates
                       (mirror-x grid-n
                                 (mirror-y grid-n
                                           (mirror-d
                                            atoms/gen)))))

    (define atoms (for/list ([a (in-list atoms/int)])
                    (match a
                      [(list x y w)
                       (list (+ border-size (* grid-size x))
                             (+ border-size (* grid-size y))
                             (/ w min-weight))])))
    (values L atoms))}

{define-values (L-Massaccesi atoms-Massaccesi)
  (let ()
    (define grid-n 29)
    (define L 45058/10000)
    (define M 15513/10000)
    (define grid-size (/ (- L M) grid-n))
    (define border-size (/ M 2))
    (define min-weight 576)
    {define atoms/gen '(#;()
                        (0 2 165) (0 11 129) (1 8 36) (1 10 21) (1 11 15)
                        (2 2 246) (2 8 129) (2 9 105) (2 10 36) (2 11 105)
                        (5 10 36) (6 10 63) (6 11 12) (7 10 21)
                        (8 9 33) (8 11 15) (9 11 75) (9 14 39)
                        (10 11 25) (10 12 21) (10 13 24) (10 14 3) (11 11 16))}

    (define atoms/int (remove-duplicates
                       (mirror-x grid-n
                                 (mirror-y grid-n
                                           (mirror-d
                                            atoms/gen)))))

    (define atoms (for/list ([a (in-list atoms/int)])
                    (match a
                      [(list x y w)
                       (list (+ border-size (* grid-size x))
                             (+ border-size (* grid-size y))
                             (/ w min-weight))])))
    (values L atoms))}

{define (draw-example L atoms #:gamma [gamma 2.0] #:scale [scale 0.07])
  [with-window (window -.1 (+ L .1) -.1 (+ L .1))
    (define big-fill-color "whitesmoke")
    (define big-border-color "black")
    (define dots-color (change-alpha "darkred" 0.75))
    (define big-square (curve (pt 0 0) -- (pt 0 L) -- (pt L L) -- (pt L 0) -- cycle))
    (draw (color big-fill-color (fill big-square))
          (penscale .1 (color big-border-color (draw big-square)))
          (draw* (for/list ([a (in-list atoms)])
                   (match a
                     [(list x y w)
                      (define s (* (sqrt (expt w (/ 1. gamma))) scale)) 
                      (penstyle 'transparent (color dots-color (filldraw (circle (pt x y) s))))]))))
    ]}

(scale 4 (draw-example L-Green atoms-Green #:gamma 2.0 #:scale .07))

(scale 4 (draw-example L-Green/S atoms-Green/S #:gamma 2.0 #:scale .07))

(scale 4 (draw-example L-Burns atoms-Burns #:gamma 2.0 #:scale .07))

(scale 4 (draw-example L-Massaccesi atoms-Massaccesi #:gamma 2.0 #:scale .07))

► Versión en Español 

The other day, I was reading yet another benchmark: “Lambda Land: Functional Languages Need Not Be Slow”. (Is it a benchmark? Nice article anyway.) It compares Racket, Python, Rust, Julia and JavaScript trying the Collatz conjecture up to 500,000. 

First, two disclaimers:

- There are lies, damned lies, statistics and benchmarks.

- The ultimate way of cheating in a benchmark is changing the compiler.

The idea is to try to improve the time of the Racket program following the spirit of the benchmark, and then extract some ideas to improve the compiler to automate the improvements when possible. The Racket compiler transforms the program to Chez Scheme code, so most of the optimization pass improvements will be actually in the Chez Scheme compiler and they will hopefully make programs in both languages faster. 

(I don’t know enough of the other languages to try something similar, so I won’t even try. I’d like to read any similar analysis.)

Each benchmark is different. In my opinion, the spirit of this benchmark is to write nice code and see what the compiler can do. I’m more used to benchmarks that try to squeeze every single millisecond and use curse primitives like #3%$unsafe-fl*+! so this one has a very different spirit, only allowing nice code.

With the changes, the run times are:

Program Version                                    Seconds

Original from Lambda Land       (Racket 9.0)          17.9

Faster now, but not so nice     (Racket 9.0)           5.3

Faster in the future and nice   (Racket 9.0)          17.8 

Faster in the future and nice   (Racket 9.3?)          4.9 (expected?)

(All times measured at the command line, outside DrRacket. By default DrRacket has “debugging” enabled and that adds a lot of additional time in exchange for better error reports.) 

As a baseline, this is the original code with two highlighted expressions that will be the slow ones:

(define (count-collatz n [cnt 1])

  (cond

    [(= n 1) cnt]

    [(even? n) (count-collatz (/ n 2) (+ 1 cnt))]

    [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))]))

Fast but not so nice code

My first step was to write not nice code that is fast. But not very ugly code, just a little ugly. The idea is to have a tradeoff between nice and fast code. I tried a lot of variants until I understood these two expressions were the most important to improve the speed.

The changes are

• Split the function in a part for fixnums and another for bignums.

• The division / is slow, so let’s replace it with unsafe-quotient. (green part)

• Also even? is slow, so let’s split the logic and use another unsafe (blue part)

(require racket/fixnum)

(require racket/unsafe/ops)

(define (count-collatz n [cnt 1])

  (if (fixnum? n)

      (cond

        [(= n 1) cnt]

        [(zero? (unsafe-fxremainder n 2)) (count-collatz (unsafe-fxquotient n 2) (+ 1 cnt))]

        [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))])

      (cond

        [(= n 1) cnt]

        [(even? n) (count-collatz (/ n 2) (+ 1 cnt))]

        [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))])))

Not so ugly in my opinion, and the run time goes from 17.9 seconds to 5.3 seconds on my computer. Note that the first branch for fixnums has a few tricks, but the second is just to the original code.

Division and quotient

We can measure the run time when we change the green expression:

Blue                              Green                    Seconds

(even? n)                         (/ n 2)                     17.9

(even? n)                         (quotient n 2)              17.8

(even? n)                         (fxquotient n 2)            17.3

(even? n)                         (unsafe-fxquotient n 2)      7.6

Let’s think of them like a chain of transformations. 

Changing / to quotient is super difficult for the compiler. The optimization pass doesn't know enough algebra to make this transformation. (I think someday this is possible in very specific cases like this one, that uses modulo 2 and a predicate, but just assume it’s impossible.) Some of the examples in the other languages use the integer division, so I think it’s “fair” to use quotient here.

After this first change the time improves a little and the difference with fxquotient is small too. But the last version with unsafe-fxquotient is much faster. The good news is that from the expression that uses quotient it’s possible for an improved version of the compiler to do the replacements during the compilation and get the faster code.

(I also tried arithmetic-shift. There are a few weird things to check there too.) 

The even? predicate

Now we can measure the run time when we change the blue expression:

Blue                              Green                    Seconds

(even? n)                         (unsafe-fxquotient n 2)      7.6

(zero? (remainder n 2))           (unsafe-fxquotient n 2)     20.8

(zero? (fxremainder n 2))         (unsafe-fxquotient n 2)     19.6

(zero? (unsafe-fxremainder n 2))  (unsafe-fxquotient n 2)      5.3

(zero? (bitwise-and n 1))         (unsafe-fxquotient n 2)      4.9

(zero? (fxand n 1))               (unsafe-fxquotient n 2)      4.9

(not (bitwise-bit-set? n 0))      (unsafe-fxquotient n 2)      5.4

Again, let’s think of them like a chain of transformations. 

Changing even? to (zero? (remainder _ 2)) is super bad, I’m not sure why. It’s worth checking in the future. Some of the other languages use similar expressions, so it’s a “fair” change but also it’s worrying that it’s so slow because someone may use it.

Anyway, after the initial change, it is possible for an improved version of the compiler to change the primitive during compilation to fxremainder and then to unsafe-fxremander and would make the code faster than the original version.

The next two versions with bitwise-and and fxand are slightly faster, but they are too different from the original code, so I don’t classify them as “fair”. In both cases the compiler changes them to use the equivalent of unsafe-fxand.

The last version with (not (bitwise-bit-set? _ 0)) is as fast as the version that uses (zero? (unsafe-fxremainder _ 2)), but it uses a predicate. This is more friendly for the optimization pass that understands better predicates than functions with binary results. This may be relevant in a distant future to allow the compiler to replace / with quotient.

Anyway, it’s better to make the optimization pass just transform even? to the unsafe version of the Chez Scheme primitive cs:fxeven? that is as fast as (zero? (fxand _ 1)).

A macro for the happy path

One of the nice features of Racket is that you can use macros to make weird transformations in the code. (If you don’t want your coworkers and future you to hate you, then use macros wisely, document them, use syntax-parse to get nice error messages and avoid using macros when there is another option.)

So we will define a macro happy-path

(define-syntax-rule (happy-path test 

                      body ...)

  (if test

     (begin body ...)

     (begin body ...)))

It’s a very silly macro. When test is true then it runs body ... and when test is false it runs body ... too! The interesting part is that the compiler can apply different optimizations to each branch. In particular, in our case test will be (fixnum? n) so an improved version of the compiler can make in the first branch all the changes we discussed before and so we get nice code that is fast. The final version is  

(define (count-collatz n [cnt 1])

  (happy-path (fixnum? n)

    (cond

      [(= n 1) cnt]

      [(even? n) (count-collatz (quotient n 2) (+ 1 cnt))]

      [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))])))

Conclusions

• If you know that all the arguments and the result are integers, use quotient instead of /. I think this is a good recommendation for all languages and compilers.

• Try happy-path to give an opportunity for the optimizer to make improvements in the most expected path, or both if you are very lucky. Probably (fixnum? n) and (flonum? n) are useful tests in general. It would be nice that the compiler can do that for you, but it’s difficult without making all the executable code twice as big and perhaps not getting any speed improvements.

• I expect these improvements to land on Racket 9.3 (mid 2026) so use your time machine and upgrade. I made a proof of concept branch [for Racket, for Chez Scheme], but the idea is to make more general versions of the improvements. (For example, if even? is magic, it’s nice that odd? is magic too, to keep the balance of the universe and avoid surprising users.) 

Code

Original from Lambda Lang Blog

#lang racket

(define (count-collatz n [cnt 1])

  (cond

    [(= n 1) cnt]

    [(even? n) (count-collatz (/ n 2) (+ 1 cnt))]

    [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))]))

(define (count-collatz-upto n)

  (for/fold

      ([max-seen 0])

      ([i (in-range 1 n)])

    (max max-seen (count-collatz i))))

(displayln (format "\nDone ~a" (count-collatz-upto 5000000)))

Faster now, but not so nice (Racket 9.0)

#lang racket

(require racket/fixnum)

(require racket/unsafe/ops)

(define (count-collatz n [cnt 1])

  (if (fixnum? n)

      (cond

        [(= n 1) cnt]

        [(zero? (unsafe-fxremainder n 2)) (count-collatz (unsafe-fxquotient n 2) (+ 1 cnt))]

        [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))])

      (cond

        [(= n 1) cnt]

        [(even? n) (count-collatz (/ n 2) (+ 1 cnt))]

        [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))])))

Faster in the future and nice (Racket 9.3?)

#lang racket

(define-syntax-rule (happy-path test

                      body ...)

  (if test

     (begin body ...)

     (begin body ...)))

(define (count-collatz n [cnt 1])

  (happy-path (fixnum? n)

    (cond

      [(= n 1) cnt]

      [(even? n) (count-collatz (quotient n 2) (+ 1 cnt))]

      [else (count-collatz (+ (* 3 n) 1) (+ 1 cnt))])))

(define (count-collatz-upto n)

  (for/fold

      ([max-seen 0])

      ([i (in-range 1 n)])

    (max max-seen (count-collatz i))))

(displayln (format "\nDone ~a" (count-collatz-upto 5000000)))