Book Ramsey Numbers: Verification Notebook

This notebook is a self-contained mathematical and computational companion to the Book Ramsey numbers section in the paper.

The task templates, including the prompts given to the agents, are available for the first Station and the second Station.

The outline is as follows:

1. Problem Definitions
2. S1. A conference-graph family.
3. S2. A doubled Legendre family.
4. S3. A Yamada–Pott family.
5. Additional findings
6. Appendix: coverage for $1\le n\le200$

The required packages are:

NumPy

1. Problem Definitions

Given graphs $G_1,G_2$ , the Ramsey number $R(G_1,G_2)$ is the smallest $n$ such that every red-blue edge coloring of $K_n$ forces either a red copy of $G_1$ or a blue copy of $G_2$ . Establishing the exact values of Ramsey numbers is a difficult computational and theoretical challenge. The most famous Ramsey numbers are those where $G_1$ and $G_2$ are complete graphs, but many other choices have been studied extensively. The book graph $B_k$ consists of $k$ triangles that share a common edge. An open problem is whether

R(B_{n-1},B_n)=4n-1

holds for every positive integer $n$ . Rousseau and Sheehan (1978) established the upper bound, proving $R(B_{n-1},B_n)\leq4n-1$ for all $n$ . It therefore remains to prove the matching lower bound. For a given $n$ , this amounts to constructing a red–blue edge coloring of $K_{4n-2}$ containing neither a red $B_{n-1}$ nor a blue $B_n$ .

Let $G$ be the red graph on $N=4n-2$ vertices. Its edge-positive Seidel matrix is $S$ , where $S_{ij}=1$ on red edges, $S_{ij}=-1$ on blue edges, and $S_{ii}=0$ .

Lemma 1.1 (Seidel certificate). If $S\mathbf 1=-\mathbf 1$ , then $G$ is $(2n-2)$ -regular. For a red edge $ij$ having $a$ common red neighbors and a blue pair $ij$ having $b$ common blue neighbors,

(S^2)_{ij}=4(a-(n-2))\quad\text{and}\quad (S^2)_{ij}=4(b-(n-1)),

respectively. Hence nonpositive off-diagonal entries of $S^2$ certify that $G$ contains no red $B_{n-1}$ and its complement contains no blue $B_n$ .

Proof. The row-sum equation gives $2\deg_G(i)-(N-1)=-1$ , hence $\deg_G(i)=2n-2$ . For a red pair, partition the other $N-2$ vertices according to their two incident colors. There are $a$ of type red–red, $2n-3-a$ of each mixed type, and $a+2$ of type blue–blue. Adding $+1$ for equal signs and $-1$ for unequal signs gives the first formula. For a blue pair, the red–red and blue–blue classes have the same size $b$ , and each mixed class has size $2n-2-b$ , giving the second formula. $\square$

These shared routines are used by all verification cells below. The next cell authenticates the certificate bundle and defines exact adjacency, codegree, and Seidel checks.

Show code

Code cell 3 · In [1]

from __future__ import annotations

import hashlib
import itertools
from dataclasses import dataclass
from math import isqrt
from pathlib import Path
from typing import Iterable

import numpy as np

DATA_FILE = Path("book_ramsey_certificates.npz")
assert DATA_FILE.is_file(), f"Missing public certificate bundle: {DATA_FILE}"
CERTIFICATES = np.load(DATA_FILE, allow_pickle=False)

def matrix_to_masks(matrix: np.ndarray) -> tuple[int, ...]:
    matrix = np.asarray(matrix, dtype=np.uint8)
    assert matrix.ndim == 2 and matrix.shape[0] == matrix.shape[1]
    assert np.array_equal(matrix, matrix.T)
    assert not np.any(np.diag(matrix))
    assert set(np.unique(matrix)).issubset({0, 1})
    return tuple(
        sum(int(matrix[i, j]) << j for j in range(len(matrix)))
        for i in range(len(matrix))
    )

def masks_to_matrix_public(masks: tuple[int, ...]) -> np.ndarray:
    return np.array(
        [[(masks[i] >> j) & 1 for j in range(len(masks))] for i in range(len(masks))],
        dtype=np.uint8,
    )

def matrix_adjacency_string(matrix: np.ndarray) -> str:
    return "".join(
        "1" if matrix[i, j] else "0"
        for j in range(1, len(matrix)) for i in range(j)
    )

def verify_matrix(matrix: np.ndarray, n: int) -> dict[str, object]:
    masks = matrix_to_masks(matrix)
    N = len(masks)
    assert N == 4 * n - 2
    universe = (1 << N) - 1
    blue = tuple(universe ^ (1 << i) ^ masks[i] for i in range(N))
    degrees = tuple(row.bit_count() for row in masks)
    red_max = blue_max = -1
    red_bad = blue_bad = 0
    square_values = set()
    for j in range(1, N):
        for i in range(j):
            red_common = (masks[i] & masks[j]).bit_count()
            blue_common = (blue[i] & blue[j]).bit_count()
            mixed = N - 2 - red_common - blue_common
            square_values.add(red_common + blue_common - mixed)
            if (masks[j] >> i) & 1:
                red_max = max(red_max, red_common)
                red_bad += int(red_common > n - 2)
            else:
                blue_max = max(blue_max, blue_common)
                blue_bad += int(blue_common > n - 1)
    return {
        "n": n,
        "vertices": N,
        "degree_values": sorted(set(degrees)),
        "seidel_row_sums": sorted({2 * d - (N - 1) for d in degrees}),
        "seidel_square_off_diagonal": sorted(square_values),
        "red_max": red_max,
        "red_cap": n - 2,
        "blue_max": blue_max,
        "blue_cap": n - 1,
        "bad_pairs": red_bad + blue_bad,
        "valid": red_bad == blue_bad == 0,
        "sha256": hashlib.sha256(
            matrix_adjacency_string(matrix).encode("ascii")
        ).hexdigest(),
    }


def parse_adjacency_string(text: str, expected_vertices: int) -> tuple[int, ...]:
    """Decode an upper-triangular adjacency string into bit masks."""
    bits = "".join(str(text).split())
    assert set(bits) <= {"0", "1"}
    assert len(bits) == expected_vertices * (expected_vertices - 1) // 2
    masks = [0] * expected_vertices
    index = 0
    for j in range(1, expected_vertices):
        for i in range(j):
            if bits[index] == "1":
                masks[i] |= 1 << j
                masks[j] |= 1 << i
            index += 1
    return tuple(masks)

print("Loaded public arrays:", sorted(CERTIFICATES.files))

Saved output 1

Loaded public arrays: ['conference_n46', 'conference_n6', 'conference_n66', 'conference_source_q45', 'conference_source_q65', 'doubled_n6', 'doubled_n70', 'finite_n67', 'yamada_pott_n11', 'yamada_pott_n127', 'yamada_pott_n86']

2. S1. A conference-graph family.

A conference graph of order $q$ is a strongly regular graph with parameters

\left(q,\frac{q-1}{2},\frac{q-5}{4},\frac{q-1}{4}\right).

Let $A$ be its adjacency matrix and put $C=2A-J+I$ . The defining intersection numbers give

C\mathbf1=0,\qquad C^2=qI-J. \tag{1}

Theorem 2.1 (conference lift). Every conference graph of order $q$ produces a graph on $4q+2$ vertices with no red $B_q$ and no blue $B_{q+1}$ . Consequently

R(B_q,B_{q+1})=4q+3.

Proof. Form four $q$ -vertex chambers and the bulk matrix

H=\begin{pmatrix} C&I-C&I+C&-I-C\\ I-C&-C&-I-C&-I-C\\ I+C&-I-C&-C&-I+C\\ -I-C&-I-C&-I+C&C \end{pmatrix}.

With $a=(-\mathbf1,-\mathbf1,\mathbf1,\mathbf1)^\top$ and $b=(-\mathbf1,\mathbf1,-\mathbf1,\mathbf1)^\top$ , attach two endpoints by

S=\begin{pmatrix}H&a&b\\a^\top&0&-1\\b^\top&-1&0\end{pmatrix}.

Every off-diagonal entry of $S$ is $\pm1$ . Equation (1) gives $S\mathbf1=-\mathbf1$ . More explicitly, if

D=(4q+3)I-2J+2C,\quad E=-2(I+J+C),\quad F=2(I-J+C),

direct block multiplication yields

S^2=\begin{pmatrix} D&0&0&E&0&0\\ 0&D&F&0&-4\mathbf1&0\\ 0&F&D&0&0&-4\mathbf1\\ E&0&0&D&0&0\\ 0&-4\mathbf1^\top&0&0&4q+1&0\\ 0&0&-4\mathbf1^\top&0&0&4q+1 \end{pmatrix}. \tag{2}

Because $C$ is zero on the diagonal and $\pm1$ off it, (2) has diagonal $4q+1$ and every off-diagonal entry lies in $\{0,-4\}$ . Apply the Seidel certificate with $n=q+1$ , then combine the resulting lower bound with the universal upper bound. $\square$

Paley conference graphs supply every prime-power $q\equiv1\pmod4$ . Beyond the Paley case, Seberry and Whiteman constructed symmetric conference matrices of order $5\cdot9^{2t+1}+1$ for every $t\ge0$ using Mathon’s construction. The resulting conference graphs have order $q=5\cdot9^{2t+1}$ , so the theorem also proves the Book Ramsey conjecture whenever

n=5\cdot9^{2t+1}+1,\qquad t\ge0.

The first member has $q=45$ and $n=46$ . The public data file contains this source and the additional known conference graph of order $65$ , which supplies $n=66$ .

The next cell defines finite-field sources and the four-chamber conference lift.

Show code

Code cell 5 · In [2]

from __future__ import annotations

import hashlib
import itertools
from dataclasses import dataclass
from math import isqrt
from typing import Iterable

CHAMBERS = ('BB', 'BR', 'RB', 'RR')
CHAMBER_INDEX = {name: i for i, name in enumerate(CHAMBERS)}

CHAMBER_RULES = {
    'BB|BB': {'R': (1, 0), 'not_R': (-1, 1)},
    'BB|BR': {'diag': (1, 0), 'R': (-1, 0), 'not_R': (1, 0)},
    'BB|RB': {'diag': (1, 0), 'R': (1, 0), 'not_R': (-1, 0)},
    'BB|RR': {'diag': (-1, 1), 'R': (-1, 1), 'not_R': (1, 0)},
    'BR|BR': {'R': (-1, 0), 'not_R': (1, 1)},
    'BR|RB': {'diag': (-1, 0), 'R': (-1, 0), 'not_R': (1, 1)},
    'BR|RR': {'diag': (-1, 0), 'R': (-1, 0), 'not_R': (1, 0)},
    'RB|RB': {'R': (-1, 0), 'not_R': (1, 1)},
    'RB|RR': {'diag': (-1, 0), 'R': (1, 0), 'not_R': (-1, 0)},
    'RR|RR': {'R': (1, 0), 'not_R': (-1, 1)},
}

ENDPOINT_RULES = {
    'u': {'BB': (-1, 0), 'BR': (-1, 1), 'RB': (1, 0), 'RR': (1, 0)},
    'v': {'BB': (-1, 0), 'BR': (1, 0), 'RB': (-1, 1), 'RR': (1, 0)},
    'u|v': (-1, 0),
}

def is_prime(value: int) -> bool:
    if value < 2:
        return False
    if value % 2 == 0:
        return value == 2
    for divisor in range(3, isqrt(value) + 1, 2):
        if value % divisor == 0:
            return False
    return True

def prime_power_decomposition(value: int) -> tuple[int, int] | None:
    if value < 2:
        return None
    if is_prime(value):
        return value, 1
    for p in range(2, isqrt(value) + 1):
        if value % p != 0 or not is_prime(p):
            continue
        remaining = value
        exponent = 0
        while remaining % p == 0:
            remaining //= p
            exponent += 1
        if remaining == 1:
            return p, exponent
    return None

def _trim(poly: list[int]) -> list[int]:
    while len(poly) > 1 and poly[-1] == 0:
        poly.pop()
    return poly

def _poly_mod(poly: list[int], divisor: list[int], p: int) -> list[int]:
    out = [coeff % p for coeff in poly]
    divisor = _trim([coeff % p for coeff in divisor[:]])
    if divisor == [0]:
        raise ZeroDivisionError('zero polynomial divisor')
    while len(out) >= len(divisor) and out != [0]:
        coeff = out[-1] % p
        offset = len(out) - len(divisor)
        if coeff:
            for i, div_coeff in enumerate(divisor):
                out[offset + i] = (out[offset + i] - coeff * div_coeff) % p
        _trim(out)
    return out

def _is_irreducible_monic(poly: list[int], p: int) -> bool:
    degree = len(poly) - 1
    if degree <= 0 or poly[-1] % p != 1 or poly[0] % p == 0:
        return False
    for divisor_degree in range(1, degree // 2 + 1):
        for coeffs in itertools.product(range(p), repeat=divisor_degree):
            divisor = list(coeffs) + [1]
            if divisor[0] % p == 0:
                continue
            if _poly_mod(poly, divisor, p) == [0]:
                return False
    return True

def find_irreducible_polynomial(p: int, degree: int) -> tuple[int, ...]:
    if degree == 1:
        return (0, 1)
    for coeffs in itertools.product(range(p), repeat=degree):
        candidate = list(coeffs) + [1]
        if _is_irreducible_monic(candidate, p):
            return tuple(candidate)
    raise RuntimeError(f'No irreducible polynomial found over F_{p} of degree {degree}')

@dataclass(frozen=True)
class FiniteField:
    p: int
    degree: int
    modulus: tuple[int, ...]

    @classmethod
    def of_size(cls, q: int) -> FiniteField:
        decomposition = prime_power_decomposition(q)
        if decomposition is None:
            raise ValueError(f'{q} is not a prime power')
        p, degree = decomposition
        return cls(p, degree, find_irreducible_polynomial(p, degree))

    @property
    def size(self) -> int:
        return self.p ** self.degree

    def coeffs(self, value: int) -> tuple[int, ...]:
        out = []
        x = value
        for _ in range(self.degree):
            out.append(x % self.p)
            x //= self.p
        return tuple(out)

    def encode(self, coeffs: Iterable[int]) -> int:
        value = 0
        factor = 1
        for coeff in coeffs:
            value += (coeff % self.p) * factor
            factor *= self.p
        return value

    def sub(self, a: int, b: int) -> int:
        if self.degree == 1:
            return (a - b) % self.p
        a_coeffs = self.coeffs(a)
        b_coeffs = self.coeffs(b)
        return self.encode((a_coeffs[i] - b_coeffs[i]) % self.p for i in range(self.degree))

    def mul(self, a: int, b: int) -> int:
        if self.degree == 1:
            return (a * b) % self.p
        a_coeffs = self.coeffs(a)
        b_coeffs = self.coeffs(b)
        product = [0] * (2 * self.degree - 1)
        for i, ai in enumerate(a_coeffs):
            for j, bj in enumerate(b_coeffs):
                product[i + j] = (product[i + j] + ai * bj) % self.p
        for deg in range(len(product) - 1, self.degree - 1, -1):
            coeff = product[deg] % self.p
            if coeff:
                offset = deg - self.degree
                for r in range(self.degree):
                    product[offset + r] = (product[offset + r] - coeff * self.modulus[r]) % self.p
        return self.encode(product[: self.degree])

def paley_conference_relation(q: int) -> tuple[int, ...]:
    if q % 4 != 1:
        raise ValueError(f'q must be congruent to 1 mod 4, got {q}')
    field = FiniteField.of_size(q)
    squares = {field.mul(x, x) for x in range(1, q)}
    masks = [0] * q
    for j in range(1, q):
        for i in range(j):
            if field.sub(j, i) in squares:
                masks[i] |= 1 << j
                masks[j] |= 1 << i
    return tuple(masks)

KNOWN_NON_PRIME_POWER_CONFERENCE_MASKS = {}

KNOWN_NON_PRIME_POWER_CONFERENCE_SOURCES = {
    45: 'Mathon PC strongly regular graph with parameters (45,22,10,11)',
    65: 'Gritsenko strongly regular graph with parameters (65,32,15,16)',
}

def known_conference_relation(q: int) -> tuple[int, ...]:
    if q in KNOWN_NON_PRIME_POWER_CONFERENCE_MASKS:
        return KNOWN_NON_PRIME_POWER_CONFERENCE_MASKS[q]
    if q % 4 == 1 and prime_power_decomposition(q) is not None:
        return paley_conference_relation(q)
    raise ValueError(
        'automatic construction requires n - 1 to be either a prime power congruent to 1 mod 4 '
        'or one of the embedded non-prime-power conference cases q=45,65'
    )

def conference_relation_source(q: int) -> str:
    if q in KNOWN_NON_PRIME_POWER_CONFERENCE_SOURCES:
        return KNOWN_NON_PRIME_POWER_CONFERENCE_SOURCES[q]
    if q % 4 == 1 and prime_power_decomposition(q) is not None:
        return f'Paley conference graph over GF({q})'
    return 'not available in this notebook'

def conference_relation_stats(r_masks: tuple[int, ...]) -> dict[str, object]:
    q = len(r_masks)
    expected_degree = (q - 1) // 2
    expected_lambda = (q - 5) // 4
    expected_mu = (q - 1) // 4
    loop_count = 0
    symmetry_mismatches = 0
    degrees = []
    adjacent_common = []
    nonadjacent_common = []
    for i, mask_i in enumerate(r_masks):
        loop_count += (mask_i >> i) & 1
        degrees.append((mask_i & ~(1 << i)).bit_count())
        for j in range(i + 1, q):
            symmetry_mismatches += int(((r_masks[i] >> j) & 1) != ((r_masks[j] >> i) & 1))
    for j in range(1, q):
        for i in range(j):
            common = (r_masks[i] & r_masks[j]).bit_count()
            if (r_masks[i] >> j) & 1:
                adjacent_common.append(common)
            else:
                nonadjacent_common.append(common)
    ok = (
        q % 4 == 1
        and loop_count == 0
        and symmetry_mismatches == 0
        and set(degrees) == {expected_degree}
        and set(adjacent_common) == {expected_lambda}
        and set(nonadjacent_common) == {expected_mu}
    )
    return {
        'ok': ok,
        'q': q,
        'degree_values': sorted(set(degrees)),
        'adjacent_common_values': sorted(set(adjacent_common)),
        'nonadjacent_common_values': sorted(set(nonadjacent_common)),
        'expected': (expected_degree, expected_lambda, expected_mu),
        'loop_count': loop_count,
        'symmetry_mismatches': symmetry_mismatches,
    }

def _decode_vertex(index: int, q: int) -> tuple[str, str | None, int | None]:
    if index < 4 * q:
        return 'chamber', CHAMBERS[index // q], index % q
    return 'endpoint', 'u' if index == 4 * q else 'v', None

def sk_rule_for_pair(i: int, j: int, q: int, r_masks: tuple[int, ...]) -> tuple[int, int]:
    type_i, label_i, coord_i = _decode_vertex(i, q)
    type_j, label_j, coord_j = _decode_vertex(j, q)
    if type_i == 'endpoint' and type_j == 'endpoint':
        return ENDPOINT_RULES['u|v']
    if type_i == 'endpoint' or type_j == 'endpoint':
        endpoint = label_i if type_i == 'endpoint' else label_j
        chamber = label_j if type_i == 'endpoint' else label_i
        return ENDPOINT_RULES[endpoint][chamber]
    if CHAMBER_INDEX[label_i] > CHAMBER_INDEX[label_j]:
        label_i, label_j = label_j, label_i
        coord_i, coord_j = coord_j, coord_i
    key = f'{label_i}|{label_j}'
    if coord_i == coord_j:
        status = 'diag'
    else:
        status = 'R' if (r_masks[coord_i] >> coord_j) & 1 else 'not_R'
    return CHAMBER_RULES[key][status]

def adjacency_string_from_masks(masks: tuple[int, ...] | list[int]) -> str:
    return ''.join('1' if (int(masks[j]) >> i) & 1 else '0' for j in range(1, len(masks)) for i in range(j))

def build_sk_adjacency_from_relation(r_masks: tuple[int, ...], *, strict: bool = True) -> str:
    check = conference_relation_stats(r_masks)
    if strict and not check['ok']:
        raise ValueError(f'R does not have the required conference parameters: {check}')
    q = len(r_masks)
    vertex_count = 4 * q + 2
    masks = [0] * vertex_count
    for j in range(1, vertex_count):
        for i in range(j):
            sign, _slack = sk_rule_for_pair(i, j, q, r_masks)
            if sign == 1:
                masks[i] |= 1 << j
                masks[j] |= 1 << i
    return adjacency_string_from_masks(masks)

def construct_book_ramsey_witness(n: int) -> str:
    q = n - 1
    if q <= 1 or q % 4 != 1:
        raise ValueError('automatic construction requires n - 1 > 1 and congruent to 1 mod 4')
    return build_sk_adjacency_from_relation(known_conference_relation(q))

def accepted_n_values(max_n: int) -> list[int]:
    return [
        n for n in range(1, max_n + 1)
        if n - 1 > 1
        and (n - 1) % 4 == 1
        and (
            prime_power_decomposition(n - 1) is not None
            or (n - 1) in KNOWN_NON_PRIME_POWER_CONFERENCE_MASKS
        )
    ]

The next cell authenticates the two non-prime-power conference sources and verifies every conference-family value with $n\le200$ .

Show code

Code cell 7 · In [3]

# Replace the two embedded non-prime-power sources by the public data arrays,
# and verify their source parameters before they can be used.
for q in (45, 65):
    source_matrix = CERTIFICATES[f"conference_source_q{q}"]
    source_masks = matrix_to_masks(source_matrix)
    expected_source_hash = {
        45: "7eb2f051b14be393d037f1d6afc460f35eae1cab6b13ca113d71da754c497029",
        65: "a8d866eaf0fe0a79f2da94518cd40a6b0c717e17311fe4038e37f53e3d34c85a",
    }[q]
    assert hashlib.sha256(matrix_adjacency_string(source_matrix).encode()).hexdigest() == expected_source_hash
    stats = conference_relation_stats(source_masks)
    assert stats["ok"], stats
    KNOWN_NON_PRIME_POWER_CONFERENCE_MASKS[q] = source_masks

conference_values = accepted_n_values(200)
assert len(conference_values) == 30
conference_checks = {}
for n in conference_values:
    adjacency = construct_book_ramsey_witness(n)
    matrix = masks_to_matrix_public(parse_adjacency_string(adjacency, 4 * n - 2))
    result = verify_matrix(matrix, n)
    assert result["valid"]
    assert result["seidel_row_sums"] == [-1]
    assert set(result["seidel_square_off_diagonal"]).issubset({-4, 0})
    conference_checks[n] = result

assert conference_checks[46]["sha256"] == "383f12d0f85f5b689d1c7509828c8bf57dc67345e336f6a88fc9b0eb5867e904"
assert conference_checks[66]["sha256"] == "8bdce7706cb90625b64eb0941c94242092704b52d9dc0af440d0f6e300e84a71"
print(f"Conference family: verified all {len(conference_checks)} values n <= 200.")
print("Values:", conference_values)

Saved output 1

Conference family: verified all 30 values n <= 200.
Values: [6, 10, 14, 18, 26, 30, 38, 42, 46, 50, 54, 62, 66, 74, 82, 90, 98, 102, 110, 114, 122, 126, 138, 150, 158, 170, 174, 182, 194, 198]

Station Reference. Archive #41 found the Paley subfamily; Archive #89 gives the general conference lift; Archive #104 identified equivalent presentations. The principal checks are Evals #4148, #4156, #4161, and #4171.

Related Work. Conference graphs and Paley conference graphs are classical. Mathon (1978) and Seberry and Whiteman (1988) constructed the non-prime-power conference-matrix family used here, and Gritsenko (2021) supplied the order- $65$ graph. Concurrently, Dai and Lin (2026) used conference matrices in a different construction for diagonal book Ramsey numbers on $4n$ vertices. The four-chamber lift here treats the near-diagonal pair $(B_q,B_{q+1})$ on $4q+2$ vertices and appears to be new.

3. S2. A doubled Legendre family.

Let $Q>3$ be a prime power with $Q\equiv3\pmod8$ , let $K=(\mathbb F_Q^\times)^2$ , and put $q=|K|=(Q-1)/2$ . If $\chi$ is the quadratic character, define

x(t)=-\chi(1+t),\qquad y(1)=-1,\quad y(t)=-\chi(1-t)\ (t\ne1).

Lemma 3.1 (periodic Legendre source). The functions $x$ and $y$ have sums $1$ and $-1$ , $x(t^{-1})=x(t)$ , and, for $\delta\ne1$ in $K$ ,

\sum_{t\in K}\{x(t)x(\delta t)+y(t)y(\delta t)\}=-2. \tag{3}

Proof. Since every $t\in K$ is a square, $x(t^{-1})=-\chi((1+t)/t)=x(t)$ . The subgroup indicator gives

\sum_{t\in K}\chi(1+t) =\frac12\sum_{t\in\mathbb F_Q^\times}(1+\chi(t))\chi(1+t)=-1,

because the ordinary and character-weighted sums are both $-1$ . Hence $\sum_Kx=1$ . Similarly,

\sum_{t\in K}\chi(1-t) =\frac12\left(-1-\chi(-1)\right)=0;

changing the zero value at $t=1$ to $y(1)=-1$ gives $\sum_Ky=-1$ . For (3), the exceptional value in $y$ contributes two seam corrections, which cancel because $\chi(1-\delta)+\chi(1-\delta^{-1})=0$ . The remaining two polynomials are

P_\pm(t)=\delta t^2\pm(\delta+1)t+1.

Both have square leading coefficient and nonzero square discriminant $(\delta-1)^2$ , so $\sum_{t\in\mathbb F_Q^\times}\chi(P_\pm(t))=-2$ . Inserting $\mathbf1_K(t)=(1+\chi(t))/2$ , the two character-weighted sums cancel under $t\mapsto-t$ because $P_+(-t)=P_-(t)$ and $\chi(-1)=-1$ . The remaining half-sum is $-2$ , proving (3). $\square$

Enumerate $K$ cyclically, let $C_x,C_y$ be the developed matrices, and set

A=C_x-I,\quad B=C_y,\quad T=\begin{pmatrix}A&B\\B^\top&-A\end{pmatrix}.

The correlation identity says

C_xC_x^\top+C_yC_y^\top=(2q+2)I-2J.

Since $C_x$ is symmetric and the developed matrices commute, block multiplication has zero cross blocks and gives the square identity below. The two row sums give the stated actions on $j$ and $r$ . Thus, with $j=\mathbf1_{2q}$ , $r=(\mathbf1_q,-\mathbf1_q)^\top$ , and $R=\operatorname{diag}(I_q,-I_q)$ , (3) gives

Tj=-j,\quad Tr=r,\quad T^2=(2q-1)I-4\operatorname{diag}(W,W),\quad W=\frac{J-I+A}{2}. \tag{4}

Here $W$ is binary with zero diagonal; the cross blocks vanish because developed matrices over the abelian group $K$ commute.

Theorem 3.2 (two-endpoint lift). Let $Q>3$ be a prime power with $Q\equiv3\pmod8$ . There is a coloring on $2Q$ vertices with no red $B_{(Q-1)/2}$ and no blue $B_{(Q+1)/2}$ . Consequently,

R(B_{(Q-1)/2},B_{(Q+1)/2})=2Q+1.

Proof. Define

M=\begin{pmatrix}T&T+R\\T+R&-T\end{pmatrix},\quad u=\binom{-r}{-r},\quad v=\binom{j}{-j},\quad L=(u\ v),\quad F=\begin{pmatrix}0&-1\\-1&0\end{pmatrix},

and $S=\left(\begin{smallmatrix}M&L\\L^\top&F\end{smallmatrix}\right)$ . The diagonal correction $R$ makes every off-diagonal entry of $S$ signed, and (4) gives $S\mathbf1=-\mathbf1$ . The bulk-square calculation is

(M^2)_{\rm diag}=2T^2+TR+RT+I,\quad (M^2)_{\rm cross}=TR-RT,

TR+RT=2\operatorname{diag}(A,A),\qquad TR-RT=\begin{pmatrix}0&-2B\\2B^\top&0\end{pmatrix}.

After adding $LL^\top$ , each diagonal $q$ -block is $(4q-1)I-8W+2A+2J$ , and the cross contribution consists of $-2(B+J)$ and $2(B^\top-J)$ . All their off-diagonal entries are $0$ or $-4$ . The two bulk-to-endpoint columns of $S^2$ are $(-2(r+j),0)^\top$ and $(0,2(r-j))^\top$ , also $0/-4$ valued, while the endpoint block equals $(4q+1)I_2$ . The Seidel certificate with $n=q+1=(Q+1)/2$ completes the proof. $\square$

The next cell defines the periodic Legendre source and its two-endpoint lift.

Show code

Code cell 10 · In [4]

from __future__ import annotations

import hashlib
import itertools
from dataclasses import dataclass
from math import isqrt
from typing import Iterable

def is_prime(value: int) -> bool:
    if value < 2:
        return False
    if value % 2 == 0:
        return value == 2
    for divisor in range(3, isqrt(value) + 1, 2):
        if value % divisor == 0:
            return False
    return True

def prime_power_decomposition(value: int) -> tuple[int, int] | None:
    if value < 2:
        return None
    if is_prime(value):
        return value, 1
    for p in range(2, isqrt(value) + 1):
        if value % p != 0 or not is_prime(p):
            continue
        remaining = value
        exponent = 0
        while remaining % p == 0:
            remaining //= p
            exponent += 1
        if remaining == 1:
            return p, exponent
    return None

def _trim(poly: list[int]) -> list[int]:
    while len(poly) > 1 and poly[-1] == 0:
        poly.pop()
    return poly

def _poly_mod(poly: list[int], divisor: list[int], p: int) -> list[int]:
    output = [coefficient % p for coefficient in poly]
    divisor = _trim([coefficient % p for coefficient in divisor[:]])
    while len(output) >= len(divisor) and any(output):
        coefficient = output[-1]
        offset = len(output) - len(divisor)
        for i, divisor_coefficient in enumerate(divisor):
            output[offset + i] = (
                output[offset + i] - coefficient * divisor_coefficient
            ) % p
        _trim(output)
    return output

def _is_irreducible_monic(poly: list[int], p: int) -> bool:
    degree = len(poly) - 1
    if degree <= 0 or poly[-1] % p != 1 or poly[0] % p == 0:
        return False
    for divisor_degree in range(1, degree // 2 + 1):
        for coefficients in itertools.product(range(p), repeat=divisor_degree):
            if coefficients[0] == 0:
                continue
            divisor = list(coefficients) + [1]
            if _poly_mod(poly, divisor, p) == [0]:
                return False
    return True

def find_irreducible_polynomial(p: int, degree: int) -> tuple[int, ...]:
    if degree == 1:
        return (0, 1)
    for coefficients in itertools.product(range(p), repeat=degree):
        candidate = list(coefficients) + [1]
        if _is_irreducible_monic(candidate, p):
            return tuple(candidate)
    raise RuntimeError(f'No irreducible polynomial over F_{p} of degree {degree}.')

@dataclass(frozen=True)
class FiniteField:
    p: int
    degree: int
    modulus: tuple[int, ...]

    @classmethod
    def of_size(cls, size: int) -> 'FiniteField':
        decomposition = prime_power_decomposition(size)
        if decomposition is None:
            raise ValueError(f'{size} is not a prime power.')
        p, degree = decomposition
        return cls(p, degree, find_irreducible_polynomial(p, degree))

    @property
    def size(self) -> int:
        return self.p ** self.degree

    def coefficients(self, value: int) -> tuple[int, ...]:
        output = []
        for _ in range(self.degree):
            output.append(value % self.p)
            value //= self.p
        return tuple(output)

    def encode(self, coefficients: Iterable[int]) -> int:
        value = 0
        factor = 1
        for coefficient in coefficients:
            value += (coefficient % self.p) * factor
            factor *= self.p
        return value

    def add(self, a: int, b: int) -> int:
        return self.encode(
            (x + y) % self.p
            for x, y in zip(self.coefficients(a), self.coefficients(b))
        )

    def sub(self, a: int, b: int) -> int:
        return self.encode(
            (x - y) % self.p
            for x, y in zip(self.coefficients(a), self.coefficients(b))
        )

    def mul(self, a: int, b: int) -> int:
        if self.degree == 1:
            return (a * b) % self.p
        a_coefficients = self.coefficients(a)
        b_coefficients = self.coefficients(b)
        product = [0] * (2 * self.degree - 1)
        for i, ai in enumerate(a_coefficients):
            for j, bj in enumerate(b_coefficients):
                product[i + j] = (product[i + j] + ai * bj) % self.p
        for degree in range(len(product) - 1, self.degree - 1, -1):
            coefficient = product[degree]
            offset = degree - self.degree
            for i, modulus_coefficient in enumerate(self.modulus):
                product[offset + i] = (
                    product[offset + i] - coefficient * modulus_coefficient
                ) % self.p
        return self.encode(product[:self.degree])

    def pow(self, a: int, exponent: int) -> int:
        result = 1
        while exponent:
            if exponent & 1:
                result = self.mul(result, a)
            a = self.mul(a, a)
            exponent //= 2
        return result

def distinct_prime_divisors(value: int) -> tuple[int, ...]:
    divisors = []
    divisor = 2
    while divisor * divisor <= value:
        if value % divisor == 0:
            divisors.append(divisor)
            while value % divisor == 0:
                value //= divisor
        divisor += 1
    if value > 1:
        divisors.append(value)
    return tuple(divisors)

def primitive_element(field: FiniteField) -> int:
    order = field.size - 1
    divisors = distinct_prime_divisors(order)
    for candidate in range(2, field.size):
        if all(
            field.pow(candidate, order // divisor) != 1
            for divisor in divisors
        ):
            return candidate
    raise RuntimeError(f'No primitive element found in GF({field.size}).')

def doubled_skew_paley_source(Q: int) -> dict[str, object]:
    if Q <= 3 or Q % 8 != 3 or prime_power_decomposition(Q) is None:
        raise ValueError('Q must be a prime power greater than 3 and congruent to 3 mod 8.')

    field = FiniteField.of_size(Q)
    generator = primitive_element(field)
    square_generator = field.mul(generator, generator)
    q = (Q - 1) // 2

    subgroup = []
    value = 1
    for _ in range(q):
        subgroup.append(value)
        value = field.mul(value, square_generator)
    assert value == 1 and len(set(subgroup)) == q

    def chi(a: int) -> int:
        if a == 0:
            return 0
        return 1 if field.pow(a, q) == 1 else -1

    x = tuple(-chi(field.add(1, t)) for t in subgroup)
    y = tuple(-1 if t == 1 else -chi(field.sub(1, t)) for t in subgroup)
    return {
        'Q': Q,
        'q': q,
        'field': field,
        'primitive_element': generator,
        'square_subgroup_generator': square_generator,
        'x': x,
        'y': y,
    }

def periodic_legendre_audit(x: tuple[int, ...], y: tuple[int, ...]) -> dict[str, object]:
    q = len(x)
    correlations = [
        sum(
            x[t] * x[(t + shift) % q]
            + y[t] * y[(t + shift) % q]
            for t in range(q)
        )
        for shift in range(1, q)
    ]
    audit = {
        'q': q,
        'x0': x[0],
        'x_sum': sum(x),
        'y_sum': sum(y),
        'x_inverse_symmetric': all(x[d] == x[-d % q] for d in range(q)),
        'correlation_values': sorted(set(correlations)),
    }
    audit['valid'] = (
        audit['x0'] == 1
        and audit['x_sum'] == 1
        and audit['y_sum'] == -1
        and audit['x_inverse_symmetric']
        and audit['correlation_values'] == [-2]
    )
    return audit

def periodic_legendre_source_matrix(
    x: tuple[int, ...],
    y: tuple[int, ...],
) -> tuple[tuple[int, ...], ...]:
    q = len(x)
    A = tuple(
        tuple(x[(j - i) % q] - int(i == j) for j in range(q))
        for i in range(q)
    )
    B = tuple(
        tuple(y[(j - i) % q] for j in range(q))
        for i in range(q)
    )
    T = [[0] * (2 * q) for _ in range(2 * q)]
    for i in range(q):
        for j in range(q):
            T[i][j] = A[i][j]
            T[i][q + j] = B[i][j]
            T[q + i][j] = B[j][i]
            T[q + i][q + j] = -A[i][j]
    return tuple(tuple(row) for row in T)

def compile_two_endpoint_graph(
    T: tuple[tuple[int, ...], ...],
) -> tuple[int, ...]:
    q = len(T) // 2
    source_order = 2 * q
    vertex_count = 4 * q + 2
    grading = (1,) * q + (-1,) * q
    masks = [0] * vertex_count

    def sign(i: int, j: int) -> int:
        if j < 4 * q:
            super_i, source_i = divmod(i, source_order)
            super_j, source_j = divmod(j, source_order)
            if super_i == super_j:
                return T[source_i][source_j] if super_i == 0 else -T[source_i][source_j]
            correction = grading[source_i] if source_i == source_j else 0
            return T[source_i][source_j] + correction

        if i < 4 * q:
            super_i, source_i = divmod(i, source_order)
            if j == 4 * q:
                return -grading[source_i]
            return 1 if super_i == 0 else -1

        return -1

    for j in range(1, vertex_count):
        for i in range(j):
            if sign(i, j) == 1:
                masks[i] |= 1 << j
                masks[j] |= 1 << i
    return tuple(masks)

def adjacency_string_from_masks(masks: tuple[int, ...]) -> str:
    return ''.join(
        '1' if (masks[j] >> i) & 1 else '0'
        for j in range(1, len(masks))
        for i in range(j)
    )

def construct_doubled_skew_paley_witness(Q: int) -> str:
    source = doubled_skew_paley_source(Q)
    audit = periodic_legendre_audit(source['x'], source['y'])
    if not audit['valid']:
        raise AssertionError(f'Periodic Legendre audit failed: {audit}')
    T = periodic_legendre_source_matrix(source['x'], source['y'])
    return adjacency_string_from_masks(compile_two_endpoint_graph(T))

def family_n_values(
    minimum_n: int = 1,
    maximum_n_exclusive: int = 200,
) -> list[int]:
    values = []
    for Q in range(3, 2 * maximum_n_exclusive):
        if Q <= 3 or Q % 8 != 3 or prime_power_decomposition(Q) is None:
            continue
        n = (Q + 1) // 2
        if minimum_n <= n < maximum_n_exclusive:
            values.append(n)
    return values

KNOWN_NON_PRIME_POWER_CONFERENCE_ORDERS = {45, 65}
KNOWN_FINITE_OVERLAPS = {22, 34}

The next cell constructs and verifies every doubled-Legendre family member with $n\le200$ .

Show code

Code cell 12 · In [5]

doubled_values = family_n_values(1, 201)
assert len(doubled_values) == 21
doubled_checks = {}
for n in doubled_values:
    adjacency = construct_doubled_skew_paley_witness(2 * n - 1)
    matrix = masks_to_matrix_public(parse_adjacency_string(adjacency, 4 * n - 2))
    result = verify_matrix(matrix, n)
    assert result["valid"]
    assert result["seidel_row_sums"] == [-1]
    assert set(result["seidel_square_off_diagonal"]).issubset({-4, 0})
    doubled_checks[n] = result
print(f"Doubled Legendre family: verified all {len(doubled_checks)} values n <= 200.")
print("Values:", doubled_values)

Saved output 1

Doubled Legendre family: verified all 21 values n <= 200.
Values: [6, 10, 14, 22, 30, 34, 42, 54, 66, 70, 82, 90, 106, 114, 122, 126, 142, 154, 166, 174, 190]

Station Reference. Archive #50 supplies the periodic-Legendre source identity, and Archive #51 supplies the two-endpoint compiler and family. The principal executable records are Evals #1657, #1681, and #1702.

Related Work. Periodic Legendre pairs and their Hadamard-matrix applications are classical; see Fletcher, Gysin, and Seberry (2001). Turturean (2026) used a different lift from the same type of source, reaching $n=(Q+1)/4$ . The two-endpoint lift verified here reaches $n=(Q+1)/2$ and appears to be new. Its $n=70$ member was independently announced by Epoch AI (2026) later in July.

4. S3. A Yamada–Pott family.

Let $q\ge7$ be a prime power with $q\equiv3\pmod4$ , put $m=(q-1)/2$ , choose a square generator $g$ of $\mathbb F_q^\times$ , and let $\chi(0)=0$ be the quadratic character. On $\mathbb Z_m$ define

u_r=\chi(1+g^r),\quad h_r=\chi(1-g^r),\quad X=\{r:u_r=1\}.

Then $u_{-r}=u_r$ , $h_{-r}=-h_r$ , $\sum_ru_r=-1$ , and $|X|=(m-1)/2$ . For $a\ne0$ , put

H(a)=\sum_r h_rh_{r+a},\qquad K(a)=|X\cap(X+a)|.

Lemma 4.1 (Yamada–Pott correlation identity). For every $a\ne0$ ,

H(a)+4K(a)=m-4. \tag{5}

Proof. Since $u_r=2\mathbf1_X(r)-1$ , $\sum_ru_ru_{r+a}=2-m+4K(a)$ . With $t=g^a\ne1$ , the sum of this correlation and $H(a)$ is the sum over nonzero squares $x$ of

\chi((1+x)(1+tx))+\chi((1-x)(1-tx)).

Replace the square-subgroup sum by one half of the sum over $\mathbb F_q^\times$ weighted by $1+\chi(x)$ . Each unweighted quadratic has square leading coefficient, nonzero discriminant $(1-t)^2$ , and sum $-2$ over $\mathbb F_q^\times$ ; the two weighted sums cancel under $x\mapsto-x$ because $\chi(-1)=-1$ . The displayed sum is therefore $-2$ , which rearranges to (5). $\square$

Set $E=\{0\}\cup\{\pm g^x:x\in X\}$ and $D=\mathbb F_q\setminus E$ . Then $|E|=m$ , $|D|=m+1$ . The square and nonsquare parts of $E$ cannot intersect, so $|E\cap g^aE|=1+2K(a)$ . Inclusion–exclusion gives

I(a):=|D\cap g^aD|=q-2|E|+|E\cap g^aE|=2+2K(a).

Thus (5) is equivalent to

H(a)+2I(a)=m. \tag{6}

Work in the additive group $\mathbb Z_m\times\mathbb F_q$ . For a difference $(r,c)$ define

\begin{aligned} A&=\{(r,c):r\ne0,\ c=0\text{ or }h_r\chi(c)=1\},\\ B&=\{(r,c):r\ne0,\ c=0\text{ or }h_r\chi(c)=-1\},\\ C&=\{(r,c):cg^{-r}\in D\}. \end{aligned} \tag{7}

Let $\Gamma_q$ be the graph on $\{0,1\}\times(\mathbb Z_m\times\mathbb F_q)$ whose within-fibre difference sets are $A$ and $B$ in the two respective fibres and whose cross-fibre difference set is $C$ .

Theorem 4.2 (Yamada–Pott lift). With $n=(q^2-q+2)/4$ , the graph $\Gamma_q$ has $4n-2$ vertices and proves

R(B_{(q^2-q-2)/4},B_{(q^2-q+2)/4})=q^2-q+1.

Proof. Write $R_U(d)=|U\cap(U+d)|$ . For $r\ne0$ , the row indicators in (7) are

\mathbf1_A(r,c)=\frac{1+h_r\chi(c)+\mathbf1_{\{0\}}(c)}2,\qquad \mathbf1_B(r,c)=\frac{1-h_r\chi(c)+\mathbf1_{\{0\}}(c)}2.

Expanding their correlations and using $\sum_c\chi(c)\chi(c-b)=-1$ for $b\ne0$ gives, when $a,b\ne0$ ,

\begin{aligned} R_A(a,b)&=\frac{(m-2)(2m+3)-H(a)-2h_a\chi(b)}4,\\ R_B(a,b)&=\frac{(m-2)(2m+3)-H(a)+2h_a\chi(b)}4. \end{aligned} \tag{8}

To derive (9), write $x=cg^{-r}$ and $y=(c+b)g^{-(r+a)}$ . For $b\ne0$ , a pair $(x,y)\in D^2$ determines an admissible $r$ precisely when $g^ay-x$ lies in one prescribed nonzero quadratic-character coset. The involution $(x,y)\mapsto(-x,-y)$ exchanges the two cosets, while the excluded zero differences are exactly the $I(a)$ pairs with $x=g^ay$ . Therefore

R_C(a,b)=\frac{(m+1)^2-I(a)}2. \tag{9}

If $b=0$ , each of those $I(a)$ pairs has $m$ choices of $r$ ; if $a=0$ and $b\ne0$ , the same involution excludes the $|D|=m+1$ diagonal pairs. Together with the indicator expansion (writing the exceptional rows $r=0,-a$ separately), this gives

\begin{aligned} R_A(a,0)=R_B(a,0)&=\frac{m^2-4+mH(a)}2,&R_C(a,0)&=mI(a),\\ R_A(0,b)=R_B(0,b)&=\frac{m^2-1}{2},&R_C(0,b)&=\frac{m(m+1)}2. \end{aligned} \tag{10}

Equations (6), (8)–(10) show that a same-fibre red edge has at most $n-2$ common red neighbors, while a same-fibre red nonedge has at most $n-1$ .

For a cross-fibre difference $(a,b)$ , note that $\mathbf1_A(r,c)+\mathbf1_B(r,c) =\mathbf1_{\{r\ne0\}}(1+\mathbf1_{\{0\}}(c))$ . Every nonzero column of $C$ has $(m+1)/2$ entries and its zero column is empty. Convolving the displayed identity with $\mathbf1_C$ therefore gives the common-red count

m^2-1+\begin{cases}0,&b=0,\\ (m+1)/2-\mathbf1_C(a,b),&b\ne0. \end{cases}

Hence it is $n-2$ on a red cross edge and at most $n-1$ otherwise. Finally, $|A|=|B|=m^2-1$ and $|C|=m(m+1)$ , so the graph is regular of degree $mq-1=2n-2$ . For a red nonedge, regularity makes its red and blue common-neighbor counts equal. Thus every red and blue book bound holds. The universal upper bound gives equality. $\square$

The next cell defines the Yamada–Pott source, checks its correlation identity, and constructs the lifted graphs.

Show code

Code cell 15 · In [6]

from __future__ import annotations

import hashlib
from math import isqrt

def is_prime(value: int) -> bool:
    if value < 2:
        return False
    if value % 2 == 0:
        return value == 2
    return all(value % divisor for divisor in range(3, isqrt(value) + 1, 2))

def distinct_prime_divisors(value: int) -> tuple[int, ...]:
    divisors = []
    candidate = 2
    while candidate * candidate <= value:
        if value % candidate == 0:
            divisors.append(candidate)
            while value % candidate == 0:
                value //= candidate
        candidate += 1
    if value > 1:
        divisors.append(value)
    return tuple(divisors)

def yamada_pott_source(Q: int) -> dict[str, object]:
    """Construct the compact, search-free source over the field of order Q."""
    if Q < 7 or Q % 4 != 3 or prime_power_decomposition(Q) is None:
        raise ValueError('Q must be a prime power at least 7 and congruent to 3 mod 4.')

    field = FiniteField.of_size(Q)
    m = (Q - 1) // 2
    rho = primitive_element(field)
    g = field.mul(rho, rho)
    powers_list = [1]
    for _ in range(1, m):
        powers_list.append(field.mul(powers_list[-1], g))
    powers = tuple(powers_list)
    assert field.mul(powers[-1], g) == 1 and len(set(powers)) == m

    def chi(value: int) -> int:
        if value == 0:
            return 0
        return 1 if field.pow(value, m) == 1 else -1

    u = tuple(chi(field.add(1, value)) for value in powers)
    h = tuple(chi(field.sub(1, value)) for value in powers)
    X = {r for r, sign in enumerate(u) if sign == 1}
    E = {0}
    for r in X:
        E.update((powers[r], field.sub(0, powers[r])))
    D = set(range(Q)) - E
    return {
        'Q': Q,
        'field': field,
        'm': m,
        'carrier_size': m * Q,
        'n': (Q * Q - Q + 2) // 4,
        'rho': rho,
        'g': g,
        'powers': powers,
        'chi': chi,
        'u': u,
        'h': h,
        'X': X,
        'E': E,
        'D': D,
    }

def source_identity_audit(source: dict[str, object]) -> dict[str, object]:
    """Check every instance of the finite-field invariant (1)--(3)."""
    field, m, g = source['field'], source['m'], source['g']
    u, h, X, D = source['u'], source['h'], source['X'], source['D']
    rows = []
    for a in range(1, m):
        H = sum(h[r] * h[(r + a) % m] for r in range(m))
        shifted_X = {(x + a) % m for x in X}
        K = len(X & shifted_X)
        multiplier = field.pow(g, a)
        shifted_D = {field.mul(multiplier, d) for d in D}
        I = len(D & shifted_D)
        rows.append((a, H, K, I, H + 4 * K, H + 2 * I))
    return {
        'u_sum': sum(u),
        'X_size': len(X),
        'D_size': len(D),
        'u_symmetric': all(u[-r % m] == u[r] for r in range(m)),
        'h_skew': all(h[-r % m] == -h[r] for r in range(m)),
        'H_plus_4K_values': sorted({row[4] for row in rows}),
        'H_plus_2I_values': sorted({row[5] for row in rows}),
        'valid': (
            sum(u) == -1
            and len(X) == (m - 1) // 2
            and len(D) == m + 1
            and all(u[-r % m] == u[r] for r in range(m))
            and all(h[-r % m] == -h[r] for r in range(m))
            and all(row[4] == m - 4 and row[5] == m for row in rows)
        ),
    }

def affine_connection_sets(
    source: dict[str, object],
) -> tuple[set[tuple[int, int]], set[tuple[int, int]], set[tuple[int, int]]]:
    """Compile (4.3) into subsets of Z_m x the additive group of F_Q."""
    Q, field, m, g = source['Q'], source['field'], source['m'], source['g']
    h, D, chi = source['h'], source['D'], source['chi']
    inverse_powers = tuple(field.pow(g, (-r) % m) for r in range(m))
    A, B, C = set(), set(), set()
    for r in range(m):
        for c in range(Q):
            if r != 0:
                symbol = chi(c)
                if c == 0 or h[r] * symbol == 1:
                    A.add((r, c))
                if c == 0 or h[r] * symbol == -1:
                    B.add((r, c))
            if field.mul(c, inverse_powers[r]) in D:
                C.add((r, c))
    return A, B, C

def carrier_coordinates(source: dict[str, object]) -> tuple[tuple[int, int], ...]:
    """Order the carrier, preserving the published prime-field certificate order."""
    Q, field, m = source['Q'], source['field'], source['m']
    if field.degree == 1:
        carrier = tuple((x % m, x % Q) for x in range(m * Q))
    else:
        carrier = tuple((r, c) for r in range(m) for c in range(Q))
    assert len(carrier) == len(set(carrier)) == m * Q
    return carrier

def construct_graph_rows(Q: int) -> tuple[dict[str, object], tuple[int, ...]]:
    """Generate the full graph as one Python-integer bitset per vertex."""
    source = yamada_pott_source(Q)
    audit = source_identity_audit(source)
    if not audit['valid']:
        raise AssertionError(f'Finite-field source audit failed: {audit}')
    A, B, C = affine_connection_sets(source)
    field, m = source['field'], source['m']
    carrier = carrier_coordinates(source)
    carrier_index = {value: index for index, value in enumerate(carrier)}
    carrier_size = len(carrier)
    rows = [0] * (2 * carrier_size)

    def translate(value: tuple[int, int], difference: tuple[int, int]) -> tuple[int, int]:
        r, c = value
        dr, dc = difference
        return (r + dr) % m, field.add(c, dc)

    for x, value in enumerate(carrier):
        for difference in A:
            rows[x] |= 1 << carrier_index[translate(value, difference)]
        for difference in B:
            y = carrier_index[translate(value, difference)]
            rows[carrier_size + x] |= 1 << (carrier_size + y)
        for difference in C:
            y = carrier_index[translate(value, difference)]
            rows[x] |= 1 << (carrier_size + y)
            rows[carrier_size + y] |= 1 << x
    metadata = {
        'source': source,
        'source_audit': audit,
        'A_size': len(A),
        'B_size': len(B),
        'C_size': len(C),
    }
    return metadata, tuple(rows)

def family_parameters(maximum_n_exclusive: int) -> list[tuple[int, int]]:
    output = []
    Q = 7
    while (Q * Q - Q + 2) // 4 < maximum_n_exclusive:
        if Q % 4 == 3 and prime_power_decomposition(Q) is not None:
            output.append((Q, (Q * Q - Q + 2) // 4))
        Q += 1
    return output

The next cell constructs and verifies all five Yamada–Pott values with $n\le200$ .

Show code

Code cell 17 · In [7]

yamada_values = [n for _, n in family_parameters(201)]
assert yamada_values == [11, 28, 86, 127, 176]
yamada_checks = {}
for Q, n in family_parameters(201):
    _, rows = construct_graph_rows(Q)
    result = verify_matrix(masks_to_matrix_public(rows), n)
    assert result["valid"]
    assert result["seidel_row_sums"] == [-1]
    assert max(result["seidel_square_off_diagonal"]) <= 0
    yamada_checks[n] = result
print(f"Yamada--Pott family: verified all {len(yamada_checks)} values n <= 200.")
print("Values:", yamada_values)

Saved output 1

Yamada--Pott family: verified all 5 values n <= 200.
Values: [11, 28, 86, 127, 176]

Station Reference. Archive #50 and Evals #3206, #3208, and #3210 contain the periodic-correlation identity and the finite affine constructions for $n=11,28,86$ . An external expert recognized their shared Yamada–Pott structure and combined these Station results to establish the general theorem.

Related Work. Yamada–Pott designs are classical; the source identity is described in the design-theoretic framework of Arasu, Bulutoglu, and Hollon (2020). The lift from these designs to sharp book-Ramsey colorings appears to be new.

5. Additional findings

These are additional results that may be interesting but are not included in the spotlight.

5.1 An exact finite value at $n=67$

Here we report a finite construction that settles $n=67$ : a regular two-fibre coloring on $266$ vertices built from multiplier orbits over $\mathbb Z_{133}$ . Exact common-neighbor counts, together with the universal upper bound, prove

R(B_{66},B_{67})=267.

Theorem 5.1. $R(B_{66},B_{67})=267$ .

Proof. The public data file contains a binary adjacency matrix on $266$ vertices in which every red edge has at most $65$ common red neighbors and every blue edge has at most $66$ common blue neighbors. This gives $R(B_{66},B_{67})\ge267$ ; Rousseau and Sheehan’s upper bound proves equality. The SHA-256 digest fixes the labelled matrix. $\square$

The next cell verifies the stored $n=67$ certificate directly.

Show code

Code cell 20 · In [8]

n67 = verify_matrix(CERTIFICATES["finite_n67"], 67)
assert n67["sha256"] == "819dbb3122cbd14d543c2f60588e770d2e5ae00fcd0f98c1bb06df88b23a1cf1"
assert n67["red_max"] == 65 and n67["blue_max"] == 66 and n67["valid"]
n67

Saved output 1

{'n': 67,
 'vertices': 266,
 'degree_values': [132],
 'seidel_row_sums': [-1],
 'seidel_square_off_diagonal': [-8, -4, 0],
 'red_max': 65,
 'red_cap': 65,
 'blue_max': 66,
 'blue_cap': 66,
 'bad_pairs': 0,
 'valid': True,
 'sha256': '819dbb3122cbd14d543c2f60588e770d2e5ae00fcd0f98c1bb06df88b23a1cf1'}

Station Reference. The multiplier-orbit construction is Archive #28. Its search, construction, and checks are Evals #766, #768, #769, #772, #774, #776, #779, #780, #783, #784, and #786.

Related Work. Rousseau and Sheehan (1978) proved the matching universal upper bound. Relative to Wesley (2026), Lidický et al. (2025), Turturean (2026), and Epoch AI (2026), the verified $n=67$ construction lies outside the previously reported finite range and the three infinite families above, and this exact value appears to be new.

5.2 Reconstructing a symmetric set from its autocorrelation

Construction searches often produce a desired autocorrelation profile before establishing whether any binary set realizes it. For a centrally symmetric subset $S$ of an odd-order finite abelian group, we prove the exact parity rule

x\in S\quad\Longleftrightarrow\quad r(2x)\ \text{is odd},

where $r$ is the labelled periodic autocorrelation of $S$ . Within the centrally symmetric class, this rule recovers every membership bit and yields a complete quadratic-time realizability test. It therefore removes the reconstruction ambiguity within the symmetric source class used by the book-Ramsey searches.

Theorem 5.2 (parity decoder). Let $G$ be a finite abelian group of odd order, let $S=-S\subseteq G$ , and define $r(d)=|S\cap(S+d)|$ . Then

x\in S\quad\Longleftrightarrow\quad r(2x)\equiv1\pmod2.

Thus, among centrally symmetric subsets of $G$ , the labelled periodic autocorrelation determines $S$ uniquely. Given a proposed profile $r$ , reconstructing $S_r=\{x:r(2x)\text{ is odd}\}$ and comparing its autocorrelation with $r$ gives an exact test for realization by a centrally symmetric subset.

Proof. The elements counted by $r(2x)$ can be written as those $y$ for which $y\in S$ and $y+2x\in S$ . The involution $y\mapsto-y-2x$ preserves this set because $S=-S$ . All nonfixed elements occur in pairs. Since $G$ has odd order, the unique possible fixed point is $y=-x$ , and it occurs exactly when $x\in S$ . This proves the parity rule. Doubling is bijective on an odd-order group, so every membership bit is recovered. To test a proposed profile, form $S_r=\{x:r(2x)\text{ is odd}\}$ and compare its actual autocorrelation with $r$ . $\square$

Station Reference. Archive #7 gives the theorem and algorithm, checked in Evals #265, #280, and #284.

Related Work. Averkov (2009) proved a broader uniqueness theorem for centrally symmetric finite Euclidean sets from autocorrelation; Rosenblatt and Seymour (1982) give general homometry background. The theorem here is a periodic finite-group specialization with the explicit parity decoder used by the construction searches.

6. Appendix: coverage for $1\le n\le200$

The arithmetic conditions above give 30 conference values, 21 doubled-Legendre values, and five Yamada–Pott values for $1\le n\le200$ . Their union contains 43 values. Against the finite, Paley, and Legendre results publicly available when the Station discoveries were made in mid-July 2026, 28 of these values were open: 19 first supplied by the conference family, six more by the doubled family, and three more by the Yamada–Pott family. Epoch AI announced a concurrent finite construction at $n=70$ at the end of July. The count of 28 records the discovery-time impact of the Station families, not exclusivity relative to that later announcement.

The .npz file contains representative adjacency matrices, including the two non-prime-power conference sources and an independent finite witness at $n=67$ . Each array is binary, symmetric, zero-diagonal, and stored without Python objects.

The next cell fixes the expected digests, recomputes the family ledgers and discovery-time comparison, and verifies every portable array.

Show code

Code cell 25 · In [9]

EXPECTED_DIGESTS = {'conference_source_q45': '7eb2f051b14be393d037f1d6afc460f35eae1cab6b13ca113d71da754c497029', 'conference_n46': '383f12d0f85f5b689d1c7509828c8bf57dc67345e336f6a88fc9b0eb5867e904', 'conference_source_q65': 'a8d866eaf0fe0a79f2da94518cd40a6b0c717e17311fe4038e37f53e3d34c85a', 'conference_n66': '8bdce7706cb90625b64eb0941c94242092704b52d9dc0af440d0f6e300e84a71', 'conference_n6': '84ce382953bf24a148b0476e5b3179f373c7276b56768257c454df80f5e5dd90', 'doubled_n6': '29f0688397a1d4ac8f1453b7b172a6eaad529236293fe6e55ea97ce3852e98a3', 'doubled_n70': '91ca25491075a53ae4c5aec3fc27f9ecf70001dd6546ee70893cc35196e5a0ea', 'yamada_pott_n11': '22c2a5e7bb2ebee5b112a38f2b64f56a34a16eedca1542826a0985b855492506', 'yamada_pott_n86': '55a612cc15f50a96c28868ab112c88ed67a73af98610f81b2b39fd9f80e08c6b', 'yamada_pott_n127': 'c0ee4fd1b6fbcc95f89e44daf19a1866a3b8d336fd4e4734a80346466bda8f6b', 'finite_n67': '819dbb3122cbd14d543c2f60588e770d2e5ae00fcd0f98c1bb06df88b23a1cf1'}

def is_prime_public(value: int) -> bool:
    if value < 2:
        return False
    if value % 2 == 0:
        return value == 2
    return all(value % d for d in range(3, isqrt(value) + 1, 2))

def is_prime_power_public(value: int) -> bool:
    if value < 2:
        return False
    if is_prime_public(value):
        return True
    for p in range(2, isqrt(value) + 1):
        if not is_prime_public(p):
            continue
        remainder = value
        while remainder % p == 0:
            remainder //= p
        if remainder == 1:
            return True
    return False

finite = set(range(1, 22)) | {
    22, 23, 24, 26, 28, 29, 30, 32, 34, 36, 38, 39,
    40, 42, 43, 44, 46, 47, 48, 50, 52, 54, 56,
}
paley = {n for n in range(1, 201) if (2*n-1) % 4 == 1 and is_prime_power_public(2*n-1)}
legendre_prior = {n for n in range(1, 201) if (4*n-1) % 8 == 3 and is_prime_power_public(4*n-1)}
prior = finite | paley | legendre_prior
conference_set = set(conference_values)
doubled_set = set(doubled_values)
yamada_set = set(yamada_values)
station_union = conference_set | doubled_set | yamada_set
station_new = station_union - prior
assert len(station_union) == 43
assert len(station_new) == 28
assert len(conference_set - prior) == 19
assert len(doubled_set - (prior | conference_set)) == 6
assert len(yamada_set - (prior | conference_set | doubled_set)) == 3

reconstructed = {
    "conference_n6": masks_to_matrix_public(parse_adjacency_string(construct_book_ramsey_witness(6), 22)),
    "conference_n46": masks_to_matrix_public(parse_adjacency_string(construct_book_ramsey_witness(46), 182)),
    "conference_n66": masks_to_matrix_public(parse_adjacency_string(construct_book_ramsey_witness(66), 262)),
    "doubled_n6": masks_to_matrix_public(parse_adjacency_string(construct_doubled_skew_paley_witness(11), 22)),
    "doubled_n70": masks_to_matrix_public(parse_adjacency_string(construct_doubled_skew_paley_witness(139), 278)),
}
for Q in (7, 19, 23):
    meta, rows = construct_graph_rows(Q)
    reconstructed[f"yamada_pott_n{meta['source']['n']}"] = masks_to_matrix_public(rows)

data_checks = {}
for key in CERTIFICATES.files:
    matrix = CERTIFICATES[key]
    digest = hashlib.sha256(matrix_adjacency_string(matrix).encode("ascii")).hexdigest()
    assert digest == EXPECTED_DIGESTS[key]
    if key in reconstructed:
        assert np.array_equal(matrix, reconstructed[key])
    if key.startswith("conference_source_"):
        q = int(key.rsplit("q", 1)[1])
        assert conference_relation_stats(matrix_to_masks(matrix))["ok"]
    else:
        n = int(key.rsplit("n", 1)[1])
        result = verify_matrix(matrix, n)
        assert result["valid"]
        data_checks[key] = result

print("Union of the three new families below 200:", len(station_union))
print("Open at the mid-July Station discovery cutoff:", len(station_new), sorted(station_new))
print("All portable arrays passed hash and mathematical checks.")

Saved output 1

Union of the three new families below 200: 43
Open at the mid-July Station discovery cutoff: 28 [62, 66, 70, 74, 82, 86, 90, 98, 102, 106, 110, 114, 122, 126, 127, 138, 142, 150, 154, 158, 166, 170, 174, 176, 182, 190, 194, 198]
All portable arrays passed hash and mathematical checks.

6.1 Source ledger for every $n\le200$

Each entry lists every source family in this notebook that covers the parameter. The codes are: L, prior exact values for $n\le21$ ; H, prior finite heuristic constructions; P, Paley; Sz, the prior Legendre/Szekeres family; C, Station conference; D, Station doubled Legendre; YP, Station Yamada–Pott; and F, the standalone Station finite certificate at $n=67$ (Section 5.1). Station contributions are bold. Overlaps are retained. Open means only that none of these listed sources covers the parameter; it is not an exhaustive literature claim. The listed sources cover 132 of the 200 parameters. The three infinite Station families cover 43 values and newly resolve 28 relative to the stated prior baseline; the standalone value $n=67$ is not included in either headline count.

n	status / sources	n	status / sources	n	status / sources	n	status / sources
1	L, Sz	51	P	101	Open	151	Open
2	L	52	H	102	C	152	Open
3	L, P, Sz	53	Sz	103	Open	153	Open
4	L	54	H, C, D	104	Open	154	D
5	L, P, Sz	55	P	105	Sz	155	Sz
6	L, C, D	56	H	106	D	156	Open
7	L, P, Sz	57	P, Sz	107	Open	157	P
8	L	58	Open	108	Open	158	C
9	L, P	59	Open	109	Open	159	P
10	L, C, D	60	Open	110	C	160	Open
11	L, Sz, YP	61	P, Sz	111	Sz	161	Sz
12	L	62	C	112	Open	162	Open
13	L, P	63	P, Sz	113	Open	163	Open
14	L, C, D	64	Open	114	C, D	164	Open
15	L, P, Sz	65	Open	115	P	165	Sz
16	L	66	C, D	116	Open	166	D
17	L, Sz	67	F	117	P, Sz	167	Open
18	L, C	68	Open	118	Open	168	Open
19	L, P	69	P	119	Open	169	P
20	L	70	D	120	Open	170	C
21	L, P, Sz	71	Sz	121	P	171	Sz
22	H, D	72	Open	122	C, D	172	Open
23	H	73	Open	123	Sz	173	Sz
24	H	74	C	124	Open	174	C, D
25	P	75	P	125	Sz	175	P
26	H, C	76	Open	126	C, D	176	YP
27	P, Sz	77	Sz	127	YP	177	P
28	H, YP	78	Open	128	Open	178	Open
29	H	79	P	129	P	179	Open
30	H, C, D	80	Open	130	Open	180	Open
31	P	81	Open	131	Sz	181	P
32	H	82	C, D	132	Open	182	C
33	Sz	83	Sz	133	Open	183	Open
34	H, D	84	Open	134	Open	184	Open
35	Sz	85	P	135	P	185	Sz
36	H	86	YP	136	Open	186	Open
37	P	87	P, Sz	137	Sz	187	P
38	H, C	88	Open	138	C	188	Open
39	H	89	Open	139	P	189	Open
40	H	90	C, D	140	Open	190	D
41	P, Sz	91	P	141	P, Sz	191	Open
42	H, C, D	92	Open	142	D	192	Open
43	H	93	Open	143	Sz	193	Open
44	H	94	Open	144	Open	194	C
45	P, Sz	95	Sz	145	P	195	P
46	H, C	96	Open	146	Open	196	Open
47	H	97	P	147	P, Sz	197	Sz
48	H	98	C	148	Open	198	C
49	P	99	P	149	Open	199	P
50	H, C	100	Open	150	C	200	Open

1. Problem Definitions

2. S1. A conference-graph family.

3. S2. A doubled Legendre family.

4. S3. A Yamada–Pott family.

5. Additional findings

5.1 An exact finite value at n=67n=67n=67

5.2 Reconstructing a symmetric set from its autocorrelation

6. Appendix: coverage for 1≤n≤2001\le n\le2001≤n≤200

6.1 Source ledger for every n≤200n\le200n≤200

5.1 An exact finite value at $n=67$

6. Appendix: coverage for $1\le n\le200$

6.1 Source ledger for every $n\le200$