P63 optimal quadrature on the torus: n = 8, 13, 21 and 34 proved optimal

The values

E* is written as N/n⁴; in lowest terms n = 21 and 34 give 1457/64827 and 6345/668168. The exact integer scores of the grid minima are in §5.

0. Setting

A point set is P = {p₁, …, pₙ} ⊂ 𝕋², coordinates read modulo 1, coincidences allowed. B₂(t) = t² − t + 1/6 is the second Bernoulli polynomial. Let

k(t) = 1 + 6B₂(t) = 2 − 6t + 6t²  (0 ≤ t ≤ 1),  C(t₁, t₂) = k(t₁)·k(t₂),

with k extended with period 1. k(1 − t) = k(t), 1/2 ≤ k ≤ 2 on [0, 1] and k(0) = 2, so C(0) = 4. The squared error is

E(P) = (1/n²) ∑ᵢ ∑ⱼ C(pᵢ − pⱼ) − 1  (the n terms with i = j included)

The site takes S = 10⁹ and answer coordinates Xᵢ/S with integers 0 ≤ Xᵢ < S. For a grid difference t = |Xᵢ − Xⱼ| each factor is k(t/S) = (2S² − 6tS + 6t²)/S², so n²S⁴E is an integer: the score, smaller is better. The public display is √E rounded up to 12 decimals.

Fourier expansion. On [0, 1], B₂(t) = π⁻² ∑ cos(2πmt)/m² over the positive integers m, so

k(t) = ∑ r(m)·exp(2πimt)  (m ∈ ℤ),  r(0) = 1,  r(m) = 3/(π²m²)  (m ≠ 0),  Ĉ(h) = r(h₁)·r(h₂),

E(P) = ∑ Ĉ(h)·|(1/n) ∑ᵢ exp(2πi h·pᵢ)|²  (h ∈ ℤ² ∖ {0}).

So E ≥ 0, and E is the square of the worst-case error of the equal-weight rule (1/n)∑f(pᵢ) over the unit ball of the Hilbert space with reproducing kernel C; the public metric is its square root √E.

Why the 6 matters. In the family of kernels 1 + γB₂, γ is the weight of every nonzero frequency relative to the constant: r(m) = γ/(2π²m²). γ = 6 is the site’s fixed choice. It makes the nonzero weights sum to exactly 1 (k(0) = 1 + 6·(1/6) = 2), and it makes k a polynomial with integer coefficients, so the E of a grid answer is a rational number whose denominator divides n²S⁴ and the site can score it exactly as an integer. Changing γ shifts the balance between frequencies on the axes and off them; the optimal values change with it, and the optimal point sets may too. Every certificate below is for γ = 6 only and says nothing about other weights.

Symmetries. E is unchanged by translations, by relabelling, by t₁ ↦ −t₁, t₂ ↦ −t₂ and by swapping the coordinates. Translating by a grid point and reflecting X ↦ (S − X) mod S both keep an answer on the grid.

1. The continuous lower bound: an auxiliary function

For 0 ≤ p ≤ q take the symmetric cosine products

ψ(p, p)(t) = cos 2πpt₁ · cos 2πpt₂,  ψ(p, q)(t) = cos 2πpt₁ · cos 2πqt₂ + cos 2πqt₁ · cos 2πpt₂  (p < q),

and let b = ∑ a(p, q)·ψ(p, q) over 0 ≤ p ≤ q ≤ M, and g = C − b.

Lemma 1.1 (Fourier coefficients). b̂(0) = a(0, 0). For h ≠ 0: on the four axis frequencies (0, ±q), (±q, 0), b̂ = a(0, q)/2; on the four frequencies (±p, ±p), b̂ = a(p, p)/4; for 0 < p < q, on the eight frequencies (±p, ±q), (±q, ±p), b̂ = a(p, q)/4. Each frequency belongs to exactly one basis function, so b̂(h) ≥ 0 for every h ≠ 0 exactly when every non-constant coefficient a(p, q) is ≥ 0. The constant a(0, 0) has no sign condition.

Lemma 1.2. If every non-constant coefficient is nonnegative, then for any n points

∑ᵢ ∑ⱼ b(pᵢ − pⱼ) = ∑ b̂(h)·|∑ᵢ exp(2πi h·pᵢ)|² ≥ n²·a(0, 0)  (h ∈ ℤ²).

Lemma 1.3 (the bound). If moreover C ≥ b on all of 𝕋², then E(P) ≥ E♭ := a(0, 0) − b(0)/n + 4/n − 1. Precisely, there is the identity

n²·(E(P) − E♭) = ∑ g(pᵢ − pⱼ)  (i ≠ j)  + ∑ b̂(h)·|∑ᵢ exp(2πi h·pᵢ)|²  (h ≠ 0).

Proof: write out the n diagonal terms C(0) = 4; for i ≠ j write C = g + b; add and subtract the n diagonal values b(0) and use the expansion of Lemma 1.2. Both sums on the right are nonnegative. “i ≠ j” means distinct indices only; points may coincide, since C ≥ b holds at 0 too.

Lemma 1.4 (equality). Write ℓⱼ = (j/n, (kj mod n)/n). E(Lₙ) = E♭ exactly when g vanishes at every nonzero lattice difference ℓᵢ − ℓⱼ and ∑ⱼ exp(2πi h·ℓⱼ) = 0 for every h ≠ 0 with b̂(h) > 0. Since h·ℓⱼ ≡ j(h₁ + kh₂)/n, that sum is n when h₁ + kh₂ ≡ 0 (mod n) and 0 otherwise, so it suffices that the support of b avoid the dual lattice {h : h₁ + kh₂ ≡ 0 (mod n)}. As g ≥ 0 attains its minimum 0 at these points, ∇g = 0 there as well.

The certificates

Each coefficient is exactly A(ζ) + B(ζ)/π with ζ = exp(2πi/4n) and A, B rational polynomials reduced modulo the cyclotomic polynomial Φ₄ₙ (degree below φ(4n) = 24, 24, 64); for n = 8 every coefficient has the form a + (c + d√2)/π with rationals a, c, d. Every identity below is checked in ℚ(ζ)[q], with q a formal variable standing for 1/π; since π is transcendental they imply the real identities. Aₙ denotes the nonzero lattice differences folded into [0, 1/2]² (numerators over n below), the contact points.

For every n we checked exactly: (a) each A and B is invariant under ζ ↦ ζ⁻¹, so the coefficients are real; (b) at every point of Aₙ and at every unfolded nonzero lattice difference, g = 0 and ∂g/∂t₁ = ∂g/∂t₂ = 0 (using cos(2πj/n) = (ζ⁴ʲ + ζ⁻⁴ʲ)/2 and its companions, derivatives compared after dividing by 2π); (c) a(0, 0) − b(0)/n + 4/n − 1 = E*ₙ; (d) every non-constant coefficient is either exactly 0 in ℚ(ζ)[q] or proved strictly positive by interval evaluation of ∑ Aⱼ cos(2πj/4n) + π⁻¹ ∑ Bⱼ cos(2πj/4n); (e) no positive coefficient has a frequency on the dual lattice; (f) a direct rational sum gives E(Lₙ) = E*ₙ.

Every coefficient for n = 8

The n = 8 coefficients are short enough to list in full (q = 1/π). The other coefficients (0, 6), (1, 2), (1, 4), (1, 6), (2, 3), (2, 5), (3, 5), (3, 6), (4, 5) are exactly 0; no positive coefficient has a frequency on the dual lattice h₁ + 5h₂ ≡ 0 (mod 8).

Check: a(0, 0) − b(0)/8 + 1/2 − 1 = 0.903930 − 0.282348 − 0.5 = 0.121582… = 249/2048.

2. C ≥ b on the whole torus

Lemma 2.1 (reduction). C and b are invariant under t ↦ 1 − t in each coordinate, so it suffices to prove g ≥ 0 on Q = [0, 1/2]². On the closed square [0, 1]², C is the polynomial (2 − 6t₁ + 6t₁²)(2 − 6t₂ + 6t₂²), so g is smooth on Q and Taylor’s formula applies directly.

Lemma 2.2 (third derivatives). Expand b into cos·cos terms and write ω(p) = 2πp. Since k‴ = 0, k″ = 12 and |k′| ≤ 6, ∂³C/∂t₁³ = 0 and |∂³C/∂t₁²∂t₂| ≤ 72, which gives the global bounds

M₃₀ = ∑ |a|·ω(p)³,  M₂₁ = 72 + ∑ |a|·ω(p)²ω(q),  and M₁₂, M₀₃ by symmetry,  T₃ = M₃₀ + 3M₂₁ + 3M₁₂ + M₀₃.

Lemma 2.3 (contact squares). Let a ∈ Aₙ, r = 2⁻ᵏ with [a − r, a + r]² ⊂ (0, 1)², and μ ≥ 0. If

α = g₁₁(a) − (M₃₀ + M₂₁)r − μ > 0,  β = g₂₂(a) − (M₁₂ + M₀₃)r − μ > 0,  αβ > (|g₁₂(a)| + (M₂₁ + M₁₂)r)²,

then Hessian(g) ≥ μI on the whole square. With g(a) = 0 and ∇g(a) = 0, Taylor’s formula along segments inside the square gives g(t) ≥ (μ/2)|t − a|². The basic cover uses μ = 0, the quantitative cover (§3) μ = 2/1000.

Lemma 2.4 (leaf bound). For a square with centre c and half-side h, let λ = 0 if Hessian(g)(c) passes the 2×2 semidefiniteness test and otherwise λ = min(g₁₁(c) − |g₁₂(c)|, g₂₂(c) − |g₁₂(c)|, 0) (Gershgorin). Then on the whole square

g(t) ≥ g(c) − h·(|g₁(c)| + |g₂(c)|) + λh² − T₃h³/6.

Proof: second-order Taylor expansion at c with third-order remainder; the displacement d has |dᵢ| ≤ h and |d|² ≤ 2h², so for λ ≤ 0, ½dᵀH(c)d ≥ ½λ|d|² ≥ λh²; the remainder is at most T₃h³/6.

The cover. The submission supplies a dyadic quadtree partition of Q. We check that the leaves are distinct, not nested, and have total area equal to that of Q (so they partition it), then accept each leaf if it lies in a contact square, or if the bound of Lemma 2.4 is ≥ 0 (basic cover), or if it is ≥ an upper bound for dist(t, Aₙ)²/1000 on the leaf, namely min over a of [(|c₁ − a₁| + h)² + (|c₂ − a₂| + h)²]/1000 (quantitative cover). A leaf that fails is split into four by us, up to 8 further levels. All values, derivatives and third-derivative bounds are computed in fixed-point interval arithmetic at scale 2²⁰⁰ with outward rounding; π is bracketed by the alternating series of Machin’s formula, and sin, cos use Taylor series with alternating remainder on |θ| ≤ π/4.

Every leaf was accepted; none failed. Our leaf bound (Lemma 2.4) differs from the first-order bound the submission used, so for n = 21 and 34 a few leaves had to be split once more by us; the contact-square column counts the pieces that then lay in contact squares.

Proof of Theorem 1. Check (d) of §1 with Lemma 1.1 gives b̂(h) ≥ 0 for h ≠ 0; this section gives C ≥ b; by Lemma 1.3 and check (c), every n points have E ≥ E*ₙ. By check (f) (or by (b), (e) and Lemma 1.4), Lₙ attains it. ∎

The grid for n = 8. The coordinates of L₈ are multiples of 1/8 = 0.125, so L₈ is itself a nine-decimal answer, with score 64·10³⁶·249/2048 = 7781250000000000000000000000000000000; the site’s scoring formula gives the same integer. Every grid answer is also a continuous point set, so no score can be smaller: Theorem 2 holds for n = 8.

3. Localisation: a better grid answer must hug the lattice

From here on n ∈ {13, 21, 34}. The coordinates j/n of Lₙ are not nine-decimal, so the grid minimum is strictly above E*ₙ and needs a further argument.

Lemma 3.1 (quadratic growth). g(t) ≥ dist(t, Aₙ)²/1000 on Q. This is the quantitative cover of §2: in a contact square Lemma 2.3 with μ = 2/1000 gives g ≥ |t − a|²/1000 ≥ dist²/1000, and on every other leaf the lower bound for g is at least an upper bound for dist²/1000. Let Uₙ be all images of Aₙ under the coordinate reflections, modulo 1; we checked exactly that Uₙ = {(u/n, ±ku/n) : u = 1, …, n − 1}. If the folded image of t ∈ 𝕋² is within ε of Aₙ, the same signs put t within Euclidean distance ε (mod 1) of a point of Uₙ. So g(t) < ε²/1000 forces t within ε of Uₙ.

Lemma 3.2 (localisation). Let E₀ be the E of the current record answer (its score over n²S⁴) and Δ = E₀ − E*ₙ > 0. If a grid answer P has E(P) ≤ E₀, the identity of Lemma 1.3 (every term on the right nonnegative) gives g(pᵢ − pⱼ) ≤ n²(E(P) − E*ₙ) ≤ n²Δ for each i ≠ j. We checked exactly that 1000n²Δ < ε², so every difference lies within ε of Uₙ. Both coordinates of every point of Uₙ are nonzero multiples of 1/n, at distance at least 1/n > ε from 0, so no two points coincide.

Lemma 3.3 (residue rigidity). Assume gcd(k, n) = 1, 2k ≢ 0 (mod n) and 6ε < 1/n (checked for all three n). If P satisfies the conclusion of Lemma 3.2 then, after a translation by a grid point, a relabelling and possibly y ↦ −y, pⱼ = ℓⱼ + δⱼ with δ₀ = 0 and |δⱼ| < ε.

Proof: translate so that p₀ = 0. Each pᵢ (i ≠ 0) lies within ε of some (rᵢ/n, sᵢ/n) ∈ Uₙ with rᵢ ≢ 0 and sᵢ ≡ ±krᵢ. If rᵢ = rⱼ for i ≠ j, the x-component of pᵢ − pⱼ is within 2ε of 0 yet within ε of a nonzero u/n, giving 1/n < 3ε, a contradiction. So r takes each value 1, …, n − 1 once; relabel by r and put π(j) = sⱼ, π(0) = 0. For j ≠ l, pⱼ − pₗ lies within 2ε of ((j − l)/n, (π(j) − π(l))/n) and within ε of a point of Uₙ; these two points of (1/n)ℤ² are less than 3ε apart, while distinct points of (1/n)ℤ² (mod 1) are at least 1/n > 6ε > 3ε apart, so that point of Uₙ is ((j − l)/n, (π(j) − π(l))/n), i.e. π(j) − π(l) ≡ ±k(j − l) (mod n). Consecutive increments π(j + 1) − π(j) are k or −k; two adjacent increments of opposite sign would make a two-step increment 0, which must be ±2k ≢ 0. So all increments have the same sign and π(j) ≡ kj or −kj; y ↦ −y turns the second into the first. Hence pⱼ lies within ε of ℓⱼ. ∎

Lemma 3.4 (normalisation for n = 34). One may further require every δⱼ to have x-component ≤ 0. Take δᵣ with the largest x-component, translate by the grid point pᵣ and relabel j ↦ j − r. The lattice is cyclic, ℓⱼ₊ᵣ − ℓᵣ ≡ ℓⱼ, so the new perturbation is δ′ⱼ = δⱼ₊ᵣ − δᵣ, with x-component ≤ 0 and a priori |δ′ⱼ| < 2ε. Lemma 3.2 applied to the difference pⱼ₊ᵣ − pᵣ puts it within ε of a point of Uₙ; that point and ℓⱼ both lie in (1/n)ℤ² and are less than 3ε < 1/n apart, so it is ℓⱼ and |δ′ⱼ| < ε. This uses only a grid translation and a relabelling, so no answer is lost; it is used for n = 34 only, to make the enumeration finish. ∎

4. Reduction to a rational ellipsoid

The variables are δ = (δ₁, …, δₙ₋₁) ∈ ℝ²⁽ⁿ⁻¹⁾ with δ₀ = 0. For i ≠ j, ℓᵢ − ℓⱼ ≡ (x⁰, y⁰) with both coordinates in [1/n, 1 − 1/n]; the perturbation (u, v) = δᵢ − δⱼ has |u|, |v| < 2ε < 1/n, so x⁰ + u and y⁰ + v stay in (0, 1), where k is the polynomial 2 − 6t + 6t². Each ordered pair expands exactly as

k(x⁰+u)k(y⁰+v) = k(x⁰)k(y⁰) + [k′(x⁰)k(y⁰)u + k(x⁰)k′(y⁰)v] + [6k(y⁰)u² + k′(x⁰)k′(y⁰)uv + 6k(x⁰)v²] + 6[k′(x⁰)uv² + k′(y⁰)u²v] + 36u²v².

Summing and dividing by n² gives E(ℓ + δ) = E*ₙ + G·δ + Q(δ) + R(δ), where:

• the constant term is E(Lₙ) = E*ₙ;
• the linear term G vanishes: checked exactly, and it also follows from the lattice symmetry d ↦ −d with k′(1 − x) = −k′(x);
• the quadratic term is Q(δ) = δᵀMδ with an explicit rational matrix M of order 2(n − 1). An exact LDLᵀ factorisation of M − (λ/2)I has all pivots positive, so Q(δ) ≥ (λ/2)|δ|²;
• the remainder: with |k′| ≤ 6 and |u|, |v| ≤ 2ε, |6k′(x⁰)uv² + 6k′(y⁰)u²v| ≤ 36·2ε·(u² + v²) and 36u²v² ≤ 36·4ε²·min(u², v²) ≤ 72ε²(u² + v²), so each ordered pair contributes at most 72(ε + ε²)(u² + v²) in degrees three and four. Each unordered pair appears twice and ∑ᵢ<ⱼ |δᵢ − δⱼ|² = n∑|δᵢ|² − |∑δᵢ|² ≤ n|δ|², so

|R(δ)| ≤ (144/n)(ε + ε²)|δ|² ≤ η·Q(δ),  η = 288(ε + ε²)/(nλ) < 1.

Hence E(P) − E*ₙ ≥ (1 − η)Q(δ), and E(P) ≤ E₀ forces Q(δ) ≤ Δ/(1 − η). In integers: write Sℓ as an integer part base plus a fraction f ∈ [0, 1); the grid answer’s coordinates are base + z with z ∈ ℤ²⁽ⁿ⁻¹⁾ and Sδ = z − f (|Sδ| < 1500 while Sℓ is at least S/n from 0 and from S, so nothing wraps around). Therefore

(z − f)ᵀ M (z − f) ≤ T := S²Δ/(1 − η).

5. The complete integer enumeration

Factor M = LDLᵀ exactly (reconstructed entry by entry, all D positive). With w = z − f, (z − f)ᵀM(z − f) = ∑ₖ Dₖ(wₖ + ∑ L(t, k)w_t)² over t > k. Go depth-first from the last variable: at level k the fixed variables give a centre cₖ = fₖ − ∑ L(t, k)w_t, and the remaining budget ρ gives the integer interval |zₖ − cₖ| ≤ √(ρ/Dₖ). This lists every integer point of the ellipsoid. For n = 34 the x-coordinates also obey zₖ ≤ 0 (Lemma 3.4: Sδ ≤ 0 ⇔ z ≤ f ⇔ z ≤ 0, since 0 ≤ f < 1).

Our C++ implementation rounds D down and T up and replaces f and L by enclosing intervals, all at 40-bit fixed point, so it lists a superset of the exact set; arithmetic is signed 128-bit with a guard on every magnitude, and a violated guard aborts. Every listed vector is turned back into a nine-decimal answer and scored exactly with the site’s integer formula.

For n = 13 and 21 an exact rational enumeration in pure Python reproduces the same node counts and vectors. Every listed vector (13, 21 and 4 of them) scores exactly the record; we did not analyse whether they are symmetric images of one answer, and the proof does not need it.

Proof of Theorem 2 (n = 13, 21, 34). Suppose a grid answer P scores no more than the record, i.e. E(P) ≤ E₀. Lemmas 3.2–3.4 turn it, by grid translations, relabelling and reflection (all preserving the grid and the score), into Lₙ plus a perturbation whose integer offset z lies in the ellipsoid of §4, so it appears in the enumeration; every listed vector scores at least the record. So no grid answer scores lower, and the record answer itself is valid: the grid minimum is the record. It is strictly above ⌈n²S⁴E*⌉, yet both display the same. ∎

6. Conclusion and scope

All four rows are proved: the continuous minimum is E*ₙ, attained by Lₙ, and the grid minimum equals the current record (the two coincide for n = 8). These point sets are Fibonacci lattices (n = Fₘ with generator Fₘ₋₁ or its negative modulo n).

The proof does not address uniqueness or classification of the minimisers, and the submission itself explicitly makes no such claim; we do not assert that Lₙ is the only continuous minimiser, nor that the optimal grid answer is unique.

Our verification

The site replayed every certificate with two programs of its own: tools/p63-certificates.py (exact algebra, covers, localisation constants, Hessian and remainder, scoring) and tools/p63-ellipsoid-enum.cpp (integer enumeration). From the submission we read only the coefficients in certificate.json, the cover leaf lists and the answer points; its programs, their result claims and node counts were not used.

• algebra (all four n): checks (a)–(f) of §1. Real coefficients: 22, 40, 81, 219; contact and lattice-difference points checked: 9, 15, 25, 41; positive coefficients 12, 22, 45, 123 and exact zeros 9, 17, 35, 95; both E(Lₙ) and the bound identity give E*ₙ. Under a second per n.
• cover: every leaf in the table of §2, basic (all four n) and quantitative (n = 13, 21, 34), with no failures; 1–29 seconds each.
• grid-prepare: the record answer’s exact score, 1000n²Δ < ε², Uₙ = {(u, ±ku)}, gcd(k, n) = 1, 2k ≢ 0, 6ε < 1/n, zero linear term, the exact LDLᵀ of M − (λ/2)I, η < 1 and T, and writes the enumeration input.
• p63-ellipsoid-enum: the node and vector counts of §5 (about 22–34 seconds for n = 34); grid-finish scores every vector exactly. The 4 vectors for n = 34 are identical to the ones the submission lists.
• record: the submitted n = 8 answer (which attains E*₈ exactly) and the n = 13, 21, 34 grid answers (whose scores equal the current records) rescored with the site formula; scores and displays agree with the tables above.

python tools/p63-certificates.py --root ROOT algebra|record N
python tools/p63-certificates.py --root ROOT cover N --kind basic|quantitative
python tools/p63-certificates.py --root ROOT grid-prepare N --out enumN.txt
g++ -O2 -std=c++17 tools/p63-ellipsoid-enum.cpp -o p63-ellipsoid-enum
p63-ellipsoid-enum enumN.txt leavesN.txt
python tools/p63-certificates.py --root ROOT grid-finish N --leaves leavesN.txt
python tools/p63-certificates.py --root ROOT grid-python N

What is re-derived by hand rather than checked by a program: the reasoning of Lemmas 1.1–1.4 and 2.1–2.4, the reasoning of Lemmas 3.1–3.4 (the program checks only the numerical conditions they need), and the remainder estimate of §4 (the program checks η and T, plus a non-proof spot check at random points).

Differences from the submission: its reported enumeration node counts (n = 13: 9,581 and 7,782; n = 21: 247,286 and 168,411; n = 34: 117,256,345 and 117,256,191) differ from ours because counting conventions and rounding widths differ; the vector counts (13, 21, 4) agree. Its covers were built with a different leaf bound; we judged them with our own, which needed a few extra splits (§2). We found no mathematical error.

Credit

The adopted contribution earns zzzcy #308 one permanent +2 proof award. The grid records stay with their holders.

Open P63