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