Peak Autoconvolution: Verification Notebook

This notebook is a self-contained mathematical and computational companion to the Peak autoconvolution section in the paper.

The task template, including the prompt given to the agents, is available here.

The outline is as follows:

1. Problem Definitions
2. Exact verification of the $16{,}384$ -step construction
3. Additional findings

The required packages are:

NumPy

1. Problem Definitions

The peak-autoconvolution problem asks how evenly the sum of two independent random variables with the same compactly supported density can be distributed. Equivalently, it minimizes the largest value of a normalized autoconvolution.

For a nonnegative integrable function $f$ , its autoconvolution is

(f*f)(t)=\int_{\mathbb R}f(x)f(t-x)\,dx.

Normalize $f$ by $\int_{\mathbb R}f=1$ and require its support to lie in an interval of length $1/2$ . The constant in the write-up is

C_{6.2}=\inf_f\lVert f*f\rVert_\infty.

For a nonnegative vector $a=(a_0,\ldots,a_{n-1})$ , partition $[-1/4,1/4]$ into cells of width $h=1/(2n)$ , and set $f=a_i/(hS)$ on cell $i$ , where $S=\sum_i a_i$ . Then $\int_{\mathbb R}f=1$ . The convolution of two cell indicators is triangular, so $f*f$ is piecewise linear. At its breakpoints its values are

\frac{1}{hS^2}\sum_{i+j=k}a_i a_j =\frac{2n}{S^2}(a*a)_k.

A piecewise-linear function attains its maximum at a breakpoint. Therefore

\lVert f*f\rVert_\infty =R_n(a):=\frac{2n\max_k(a*a)_k}{(\sum_i a_i)^2}.

The next cell loads the shared dependencies and authenticates the two public data files.

Show code

Code cell 3 · In [1]

from __future__ import annotations

from decimal import Decimal, getcontext
from fractions import Fraction
from math import cos, pi, prod, sin
from pathlib import Path
import hashlib
import json

import numpy as np

ARTIFACT_DIR = Path.cwd()
EXPECTED_HASHES = {
    'autocorr_6-2_weights.npy': '10094879a82ea22ea499c529cdb45f259b08d3155e6a38a184ce5f2fa0165a50',
    'local_minima_n14_certificate.json': '01f05d95d24e4527fa2b2fb0ca1b61f813d419c609d4d1909cd84bffc287903a',
}

def require_file(name):
    path = ARTIFACT_DIR / name
    if not path.is_file():
        raise FileNotFoundError(f'Required public artifact is missing: {path}')
    digest = hashlib.sha256(path.read_bytes()).hexdigest()
    assert digest == EXPECTED_HASHES[name], (name, digest)
    return path

weights_path = require_file('autocorr_6-2_weights.npy')
local_certificate_path = require_file('local_minima_n14_certificate.json')
print('Public files found and SHA-256 checks passed.')


Q = Fraction

def ordinary_convolution(left, right):
    out = [0] * (len(left) + len(right) - 1)
    for i, a in enumerate(left):
        for j, b in enumerate(right):
            out[i + j] += a * b
    return out

Saved output 1

Public files found and SHA-256 checks passed.

2. Exact verification of the $16{,}384$ -step construction

Theorem 2.1. There exists a nonnegative unit-mass $16{,}384$ -step function supported on $[-1/4,1/4]$ for which

\lVert f*f\rVert_\infty<1.504473.

Consequently $C_{6.2}<1.504473$ .

Verification. The file autocorr_6-2_weights.npy encodes such a function. Every stored binary64 coefficient is an exact dyadic rational. Clearing a common power-of-two denominator produces nonnegative integers without changing $R_n$ . If $S$ is their sum, each convolution coefficient lies in $[0,S^2]$ .

The integer self-convolution is computed modulo eight distinct primes and combined by the Chinese remainder theorem. The product of the moduli exceeds $S^2$ , so every reconstructed residue is the unique integer coefficient in the allowed range. The maximum and its comparison with $1.504473=1504473/10^6$ are therefore exact. The discrete-to-continuum identity in Section 1 proves the claim. $\square$

The next cell converts the construction to integer weights and computes its exact self-convolution by NTT and CRT.

Show code

Code cell 5 · In [2]

def prime_factors(number):
    factors, divisor = [], 2
    while divisor * divisor <= number:
        if number % divisor == 0:
            factors.append(divisor)
            while number % divisor == 0:
                number //= divisor
        divisor += 1 if divisor == 2 else 2
    if number > 1:
        factors.append(number)
    return factors

def is_prime_32(number):
    if number < 2:
        return False
    for small in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
        if number % small == 0:
            return number == small
    d, s = number - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for base in (2, 3, 5, 7, 11):
        x = pow(base, d, number)
        if x in (1, number - 1):
            continue
        for _ in range(s - 1):
            x = x * x % number
            if x == number - 1:
                break
        else:
            return False
    return True

def ntt(values, inverse, modulus, primitive_root):
    length, index = len(values), 0
    for position in range(1, length):
        bit = length >> 1
        while index & bit:
            index ^= bit
            bit >>= 1
        index ^= bit
        if position < index:
            values[position], values[index] = values[index], values[position]
    block = 2
    while block <= length:
        root = pow(primitive_root, (modulus - 1) // block, modulus)
        if inverse:
            root = pow(root, modulus - 2, modulus)
        half = block // 2
        for start in range(0, length, block):
            phase = 1
            for left in range(start, start + half):
                u = values[left]
                v = values[left + half] * phase % modulus
                values[left] = (u + v) % modulus
                values[left + half] = (u - v) % modulus
                phase = phase * root % modulus
        block *= 2
    if inverse:
        inverse_length = pow(length, modulus - 2, modulus)
        for index in range(length):
            values[index] = values[index] * inverse_length % modulus

def square_convolution_mod(source, transform_length, modulus, primitive_root):
    values = [int(value) % modulus for value in source]
    values.extend([0] * (transform_length - len(values)))
    ntt(values, False, modulus, primitive_root)
    values = [value * value % modulus for value in values]
    ntt(values, True, modulus, primitive_root)
    return values[:2 * len(source) - 1]

def exact_integer_convolution(source):
    transform_length = 1
    while transform_length < 2 * len(source) - 1:
        transform_length *= 2
    moduli_and_roots = [
        (998244353, 3), (1004535809, 3), (469762049, 3), (985661441, 3),
        (167772161, 3), (754974721, 11), (1224736769, 3), (2013265921, 31),
    ]
    rows = []
    for modulus, primitive_root in moduli_and_roots:
        assert is_prime_32(modulus)
        assert (modulus - 1) % transform_length == 0
        assert all(pow(primitive_root, (modulus - 1) // factor, modulus) != 1
                   for factor in prime_factors(modulus - 1))
        rows.append(square_convolution_mod(source, transform_length, modulus, primitive_root))
    reconstructed = list(rows[0])
    combined_modulus = moduli_and_roots[0][0]
    for row, (modulus, _) in zip(rows[1:], moduli_and_roots[1:]):
        inverse = pow(combined_modulus, -1, modulus)
        reconstructed = [
            old + combined_modulus * (((residue - old) % modulus) * inverse % modulus)
            for old, residue in zip(reconstructed, row)
        ]
        combined_modulus *= modulus
    return reconstructed, combined_modulus

weights = np.load(weights_path, allow_pickle=False)
assert weights.shape == (16384,) and weights.dtype == np.dtype('float64')
assert np.all(np.isfinite(weights)) and np.all(weights >= 0)
binary_ratios = [float(value).as_integer_ratio() for value in weights]
exponents = [den.bit_length() - 1 for _, den in binary_ratios]
assert all(den == 1 << exponent for (_, den), exponent in zip(binary_ratios, exponents))
common_exponent = max(exponents)
integer_weights = [num << (common_exponent - exponent)
                   for (num, _), exponent in zip(binary_ratios, exponents)]
assert all(Fraction(value, 1 << common_exponent) == Fraction.from_float(float(original))
           for value, original in zip(integer_weights, weights))

convolution, modulus_product = exact_integer_convolution(integer_weights)
mass = sum(integer_weights)
assert modulus_product > mass * mass
assert all(0 <= value <= mass * mass for value in convolution)
peak = max(convolution)
peak_lag = convolution.index(peak)
n = len(integer_weights)
exact_score = Fraction(2 * n * peak, mass * mass)
reported_bound = Fraction(1504473, 10**6)
assert exact_score < reported_bound
getcontext().prec = 30
score_decimal = Decimal(exact_score.numerator) / Decimal(exact_score.denominator)
margin = reported_bound - exact_score
margin_decimal = Decimal(margin.numerator) / Decimal(margin.denominator)
construction_result = {
    'steps': n,
    'nonzero_coefficients': sum(value != 0 for value in integer_weights),
    'common_binary_denominator_exponent': common_exponent,
    'argmax_lag': peak_lag,
    'exact_score_decimal': str(score_decimal),
    'margin_below_1.504473': str(margin_decimal),
    'CRT_modulus_bits': modulus_product.bit_length(),
    'coefficient_upper_bound_bits': (mass * mass).bit_length(),
}
print(construction_result)

Saved output 1

{'steps': 16384, 'nonzero_coefficients': 15334, 'common_binary_denominator_exponent': 95, 'argmax_lag': 26065, 'exact_score_decimal': '1.50447289053969738158251866375', 'margin_below_1.504473': '1.09460302618417481336249606668E-7', 'CRT_modulus_bits': 237, 'coefficient_upper_bound_bits': 207}

The exact score begins $1.504472890539697\ldots$ , attained at discrete lag $26065$ .

Station Reference. The source vector is the accepted output of Eval #1415, synthesized from the sequence discussed in Archives #6, #19, #21, #39, #40, #43, #52, #67, #76, and #81. The modular certificate verifies the result independently of the evaluator’s floating-point convolution.

Related Work. Matolcsi and Vinuesa (2010) obtained $C_{6.2}\le1.50972$ . AlphaEvolve (2025) later reached $1.5032$ , followed by the $1.50286$ construction of Yüksekgönül et al. (2026) and Russell (2026)’s exact certificate, $C_{6.2}\le1.502851$ . The value certified here improves the Matolcsi–Vinuesa bound but lies above the later bounds, making it an independent numerical reproduction rather than a frontier improvement.

3. Additional findings

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

3.1 Two isolated local minima with identical support and contact sets

We prove that support and active-contact data do not determine a local optimizer in the discretized peak-autoconvolution problem. At $n=14$ , two non-reversal-related strict local minima have exactly the same positive coordinates and saturated convolution lags, but different objective values. Exact interval certificates isolate both KKT points, and the tangent-cone argument proves strict local minimality. This rigorously confirms the presence of distinct, contact-indistinguishable optimization basins and explains an obstruction to methods that identify a solution solely from its active set.

For $w\in\mathbb R^{14}_{\ge0}$ with $\sum_iw_i=1$ , put

c_k(w)=\sum_{i+j=k}w_iw_j,\qquad \Phi(w)=\max_{0\le k\le26}c_k(w),\qquad R_{14}(w)=28\Phi(w).

The support is $S(w)=\{i:w_i>0\}$ , and the contact set is $A(w)=\{k:c_k(w)=\Phi(w)\}$ .

Theorem 3.1. There exist two distinct, non-reversal-related, isolated strict local minimizers $w^{(1)},w^{(2)}$ of $R_{14}$ with the same support

S=\{0,1,\ldots,9,12,13\}

and the same contact set

A=\{6,7,8,9,12,13,14,15,16,17,18,19\},

but with disjoint value intervals

R_{14}(w^{(1)})\in[1.55775592563293864,1.55775648563293864],

R_{14}(w^{(2)})\in[1.5604673131949972,1.560467873194997368].

Their values differ by at least $0.00271082756205856$ .

Proof. Let $Z=\{10,11\}$ . Introduce the peak $t$ , contact multipliers $\lambda_k$ , and a mass multiplier $\eta$ . On the proposed support and contact set the Karush–Kuhn–Tucker equations are

\sum_{i\in S}w_i=1,\qquad c_k(w)=t\quad(k\in A),

2\sum_{k\in A}\lambda_k w_{k-i}=\eta\quad(i\in S),\qquad \sum_{k\in A}\lambda_k=1,

where $w_j=0$ outside $0\le j<14$ . These are 26 polynomial equations in

x=(w_S,t,\lambda_A,\eta).

For each supplied rational box $X$ , the code evaluates the Krawczyk operator

K(X)=x_0-CF(x_0)+(I-CJ(X))(X-x_0)

with exact rational interval arithmetic and verifies $K(X)\subset\operatorname{int}X$ . The Krawczyk theorem gives a unique exact KKT root in each box. The same calculation proves $w_i>0$ on $S$ , $\lambda_k>0$ on $A$ , $c_k<t$ off $A$ , and positive reduced costs

\rho_z=2\sum_{k\in A}\lambda_k w_{k-z}-\eta>0 \qquad(z\in Z).

It remains to prove strict local minimality. Let $G_{ki}=2w_{k-i}$ , and form

B=\begin{pmatrix}\mathbf1^{\mathsf T}&0\\G&-\mathbf1\end{pmatrix},

with columns indexed by $(d_S,\alpha)$ . The code certifies $B$ invertible throughout each root box from

\lVert I-C_BB(X)\rVert_\infty<1.

Let $d$ be a feasible tangent direction, so $\sum_i d_i=0$ and $d_z\ge0$ for $z\in Z$ . Write $a_k=Dc_k(w)[d]$ . Stationarity gives

\sum_{k\in A}\lambda_k a_k=\sum_{z\in Z}\rho_zd_z\ge0.

If every $a_k\le0$ , positivity of $\lambda$ and $\rho$ forces $a_k=0$ for every $k\in A$ and $d_z=0$ for $z\in Z$ . Then $B(d_S,0)=0$ , so $d=0$ . Thus every nonzero feasible tangent has

\max_{k\in A}Dc_k(w)[d]>0.

The unit slice of the tangent cone is compact, hence the displayed maximum has a positive uniform lower bound there. Since each $c_k$ is quadratic, this linear increase dominates the quadratic remainder sufficiently near the root. Each root is therefore an isolated strict local minimum. Finally, their disjoint value intervals show that the roots are distinct, and the coordinate boxes exclude either root from the reversal of the other. $\square$

The next cell certifies the two KKT roots, all strict inequalities, and contact-matrix regularity with exact rational interval arithmetic.

Show code

Code cell 8 · In [3]

def q(value):
    return Q(str(value))

def iadd(left, right):
    return (left[0] + right[0], left[1] + right[1])

def ineg(value):
    return (-value[1], -value[0])

def isub(left, right):
    return iadd(left, ineg(right))

def imul(left, right):
    products = (left[0] * right[0], left[0] * right[1],
                left[1] * right[0], left[1] * right[1])
    return (min(products), max(products))

def ipoint(value):
    value = q(value)
    return (value, value)

IZERO, IONE = (Q(0), Q(0)), (Q(1), Q(1))

def isum(values):
    total = IZERO
    for value in values:
        total = iadd(total, value)
    return total

def matmul_point_interval(point_matrix, interval_matrix):
    rows = len(point_matrix)
    shared = len(point_matrix[0])
    columns = len(interval_matrix[0])
    return [[
        isum(imul(ipoint(point_matrix[i][k]), interval_matrix[k][j])
             for k in range(shared))
        for j in range(columns)
    ] for i in range(rows)]

def kkt_F_point(x, n, support, active):
    m, contacts = len(support), len(active)
    w = [Q(0)] * n
    for position, index in enumerate(support):
        w[index] = x[position]
    t = x[m]
    lam = x[m + 1:m + 1 + contacts]
    eta = x[-1]
    convolution = [
        sum((w[i] * w[lag - i] for i in range(n) if 0 <= lag - i < n), Q(0))
        for lag in range(2 * n - 1)
    ]
    out = [sum(w, Q(0)) - 1]
    out.extend(convolution[lag] - t for lag in active)
    for index in support:
        out.append(sum(
            (2 * la * w[lag - index] for la, lag in zip(lam, active)
             if 0 <= lag - index < n),
            Q(0),
        ) - eta)
    out.append(sum(lam, Q(0)) - 1)
    return out

def kkt_J_interval(X, n, support, active):
    m, contacts = len(support), len(active)
    size = m + 1 + contacts + 1
    w = [IZERO for _ in range(n)]
    for position, index in enumerate(support):
        w[index] = X[position]
    lam = X[m + 1:m + 1 + contacts]
    J = [[IZERO for _ in range(size)] for __ in range(size)]
    for column in range(m):
        J[0][column] = IONE
    for row, lag in enumerate(active, 1):
        for column, index in enumerate(support):
            other = lag - index
            if 0 <= other < n:
                J[row][column] = imul(ipoint(2), w[other])
        J[row][m] = ipoint(-1)
    for row, index in enumerate(support, 1 + contacts):
        for column, varied in enumerate(support):
            J[row][column] = isum(
                imul(ipoint(2), la)
                for la, lag in zip(lam, active)
                if varied == lag - index
            )
        for lag_index, lag in enumerate(active):
            other = lag - index
            if 0 <= other < n:
                J[row][m + 1 + lag_index] = imul(ipoint(2), w[other])
        J[row][-1] = ipoint(-1)
    for lag_index in range(contacts):
        J[-1][m + 1 + lag_index] = IONE
    return J

def contact_matrix_interval(X, n, support, active):
    m = len(support)
    w = [IZERO for _ in range(n)]
    for position, index in enumerate(support):
        w[index] = X[position]
    B = [[IZERO for _ in range(m + 1)] for __ in range(m + 1)]
    for column in range(m):
        B[0][column] = IONE
    for row, lag in enumerate(active, 1):
        for column, index in enumerate(support):
            other = lag - index
            if 0 <= other < n:
                B[row][column] = imul(ipoint(2), w[other])
        B[row][-1] = ipoint(-1)
    return B

def convolution_intervals(X, n, support):
    w = [IZERO for _ in range(n)]
    for position, index in enumerate(support):
        w[index] = X[position]
    return [
        isum(imul(w[i], w[lag - i]) for i in range(n) if 0 <= lag - i < n)
        for lag in range(2 * n - 1)
    ]

def verify_local_root(root, n, support, active, zero_indices):
    x0 = [q(value) for value in root['center']]
    X = [(q(low), q(high)) for low, high in root['x_box']]
    C = [[q(value) for value in row] for row in root['kkt_preconditioner']]
    F0 = kkt_F_point(x0, n, support, active)
    JX = kkt_J_interval(X, n, support, active)
    CJ = matmul_point_interval(C, JX)

    K = []
    for row in range(len(x0)):
        center_term = x0[row] - sum(
            (C[row][k] * F0[k] for k in range(len(x0))), Q(0)
        )
        value = ipoint(center_term)
        for column in range(len(x0)):
            coefficient = isub(IONE if row == column else IZERO, CJ[row][column])
            coordinate_delta = (
                X[column][0] - x0[column],
                X[column][1] - x0[column],
            )
            value = iadd(value, imul(coefficient, coordinate_delta))
        K.append(value)
    assert all(
        X[i][0] < K[i][0] and K[i][1] < X[i][1]
        for i in range(len(X))
    )

    m, contacts = len(support), len(active)
    assert min(X[i][0] for i in range(m)) > 0
    assert min(X[m + 1 + i][0] for i in range(contacts)) > 0
    t_interval = X[m]
    convolution = convolution_intervals(X, n, support)
    assert all(
        convolution[lag][1] < t_interval[0]
        for lag in range(2 * n - 1)
        if lag not in active
    )

    w = [IZERO for _ in range(n)]
    for position, index in enumerate(support):
        w[index] = X[position]
    lambdas = X[m + 1:m + 1 + contacts]
    eta = X[-1]
    reduced_costs = []
    for zero_index in zero_indices:
        field = isum(
            imul(ipoint(2), imul(lam, w[lag - zero_index]))
            for lam, lag in zip(lambdas, active)
            if 0 <= lag - zero_index < n
        )
        reduced = isub(field, eta)
        assert reduced[0] > 0
        reduced_costs.append(reduced)

    CB = [[q(value) for value in row] for row in root['contact_preconditioner']]
    BX = contact_matrix_interval(X, n, support, active)
    CBB = matmul_point_interval(CB, BX)
    error = [[
        isub(IONE if i == j else IZERO, CBB[i][j])
        for j in range(len(CBB))
    ] for i in range(len(CBB))]
    neumann_bound = max(
        sum(max(abs(entry[0]), abs(entry[1])) for entry in row)
        for row in error
    )
    assert neumann_bound < 1
    return {
        'label': root['label'],
        'krawczyk_max_relative_width': max(
            max(abs(K[i][0] - x0[i]), abs(K[i][1] - x0[i]))
            / max(x0[i] - X[i][0], X[i][1] - x0[i])
            for i in range(len(X))
        ),
        'R_interval': (28 * t_interval[0], 28 * t_interval[1]),
        'minimum_support_lower_bound': min(X[i][0] for i in range(m)),
        'minimum_multiplier_lower_bound': min(
            X[m + 1 + i][0] for i in range(contacts)
        ),
        'minimum_reduced_cost_lower_bound': min(
            value[0] for value in reduced_costs
        ),
        'contact_matrix_Neumann_bound': neumann_bound,
        'weight_boxes': [
            X[support.index(i)] if i in support else (Q(0), Q(0))
            for i in range(n)
        ],
    }

certificate = json.loads(local_certificate_path.read_text())
assert certificate['format'] == 'autocorr-local-minima-certificate-v1'
n_local = int(certificate['n'])
support = list(map(int, certificate['support']))
active = list(map(int, certificate['active_lags']))
zero_indices = list(map(int, certificate['zero_indices']))
assert n_local == 14 and support == list(range(10)) + [12, 13]
assert active == [6, 7, 8, 9, 12, 13, 14, 15, 16, 17, 18, 19]

local_results = [
    verify_local_root(root, n_local, support, active, zero_indices)
    for root in certificate['roots']
]
first, second = local_results
assert first['R_interval'][1] < second['R_interval'][0]

def boxes_disjoint(left, right):
    return any(
        a[1] < b[0] or b[1] < a[0]
        for a, b in zip(left, right)
    )

assert boxes_disjoint(first['weight_boxes'], second['weight_boxes'])
assert boxes_disjoint(first['weight_boxes'], list(reversed(second['weight_boxes'])))
assert boxes_disjoint(second['weight_boxes'], list(reversed(first['weight_boxes'])))

for result in local_results:
    printable = {
        key: (float(value) if isinstance(value, Fraction) else value)
        for key, value in result.items()
        if key != 'weight_boxes'
    }
    printable['R_interval'] = tuple(
        float(value) for value in result['R_interval']
    )
    print(printable)

gap = second['R_interval'][0] - first['R_interval'][1]
assert gap >= Q('0.00271082756205856')
print('certified_R_gap_at_least', float(gap))

Saved output 1

{'label': 'w1', 'krawczyk_max_relative_width': 8.245452284586815e-06, 'R_interval': (1.5577559256329387, 1.5577564856329387), 'minimum_support_lower_bound': 0.01518813145928526, 'minimum_multiplier_lower_bound': 0.03590960690528001, 'minimum_reduced_cost_lower_bound': 0.0022805496484975584, 'contact_matrix_Neumann_bound': 2.150918768321847e-06}
{'label': 'w2', 'krawczyk_max_relative_width': 8.26577146888127e-06, 'R_interval': (1.560467313194997, 1.5604678731949975), 'minimum_support_lower_bound': 0.003970422477742227, 'minimum_multiplier_lower_bound': 0.009708680195212264, 'minimum_reduced_cost_lower_bound': 0.01358009196718962, 'contact_matrix_Neumann_bound': 2.5573800560188888e-06}
certified_R_gap_at_least 0.00271082756205856

Station Reference. The pair was found in Eval #1526 and reported in an addendum to Archive #81. Eval #1536 produced the original boxes. The public checker redoes the 26-dimensional Krawczyk proof, all strict inequalities, and contact-matrix regularity.

Related Work. De Dios Pont and Madrid (2021) observed numerically that the discretized maximum-autoconvolution problem appears to have multiple local extremizers and discussed the resulting obstruction to fixed-point methods. The theorem above gives a rigorous finite-dimensional instance in which multiplicity persists even after the support and active contact set are fixed; this strengthened finite statement appears novel.

3.2 The perfectly flat target requires exactly two positive sources

We determine the minimum number of nonnegative sources needed to produce a perfectly flat finite convolution profile. For every $N=2n-1$ , we prove that $1+x+\cdots+x^{N-1}$ is the sum of two squares of strictly positive degree- $(n-1)$ polynomials, while no representation by a single square exists. Thus perfect flatness, which is impossible for one finite autoconvolution, becomes attainable with exactly two sources.

Let $n\ge2$ , $N=2n-1$ , and

P_N(x)=1+x+\cdots+x^{N-1}.

For a polynomial $H(x)=\sum_{j=0}^{n-1}h_jx^j$ , the coefficients of $H(x)^2$ form the autoconvolution of $h$ .

Theorem 3.2. The smallest $r$ for which

P_N(x)=\sum_{j=1}^rH_j(x)^2

with every $H_j$ of degree at most $n-1$ and having nonnegative real coefficients is $r=2$ .

Proof. One square is impossible. The $N-1$ nontrivial $N$ th roots of unity are simple roots of $P_N=(x^N-1)/(x-1)$ , whereas every root of a polynomial square has even multiplicity.

For the upper bound put $d=n-1$ , $\omega=e^{2\pi i/N}$ , and

U_0(x)=\prod_{j=0}^{d-1}(x-\omega^{2j+1}).

The selected odd powers contain one root from each conjugate pair of nontrivial $N$ th roots. If a bar denotes coefficientwise complex conjugation, then

U_0\overline{U_0}=P_N.

Choose

\beta=\exp\!\left(i\left(\frac\pi4+\frac{\pi d}{2N}\right)\right), \qquad U=\beta U_0=A+iB.

For the reversal $F^*(x)=x^dF(1/x)$ , direct multiplication over the selected roots gives $U^*=i\overline U$ . Comparing real and imaginary parts yields $B=A^*$ . Hence

A(x)^2+A^*(x)^2=U(x)\overline U(x)=P_N(x).

It remains to show positivity. The finite $q$ -binomial identity, applied to the geometric progression of roots, gives the following real coefficient formula for $0\le m\le d$ :

A_m= \left(\prod_{r=1}^{d-m} \frac{\sin((2r-1)\pi/N)}{\sin(2r\pi/N)}\right) \cos\left(\frac{(m+1/4)\pi}{N}\right).

Here $1\le r\le d<N/2$ , so every sine in the product is positive; also $0<(m+1/4)\pi/N<\pi/2$ , so the cosine is positive. Thus $A$ and $A^*$ have strictly positive coefficients and furnish two squares. $\square$

The next cell evaluates the explicit positive source formula for $2\le n\le100$ and checks the resulting two-square identity numerically.

Show code

Code cell 11 · In [4]

def positive_flat_source(n):
    N, d = 2 * n - 1, n - 1
    return [
        prod(
            sin((2 * r - 1) * pi / N) / sin(2 * r * pi / N)
            for r in range(1, d - m + 1)
        ) * cos((m + 0.25) * pi / N)
        for m in range(n)
    ]

flat_checks = []
for test_n in range(2, 101):
    source = positive_flat_source(test_n)
    reverse = source[::-1]
    total = [
        left + right
        for left, right in zip(
            ordinary_convolution(source, source),
            ordinary_convolution(reverse, reverse),
        )
    ]
    residual = max(abs(value - 1.0) for value in total)
    assert min(source) > 0
    assert residual < 2e-10
    flat_checks.append((test_n, min(source), residual))

print({
    'n_values_checked': len(flat_checks),
    'range': (flat_checks[0][0], flat_checks[-1][0]),
    'largest_floating_residual': max(row[2] for row in flat_checks),
})

Saved output 1

{'n_values_checked': 99, 'range': (2, 100), 'largest_floating_residual': 7.771561172376096e-14}

The loop illustrates the closed formula; the proof is the exact root factorization and positivity argument above.

Station Reference. The theorem and its all- $n$ proof are from Archive #23, supported by Evals #471, #478, #486, #492, and #495.

Related Work. Every nonnegative real univariate polynomial is classically a sum of two real polynomial squares; see Magron, Safey El Din, and Schweighofer (2019). This unrestricted result does not control the signs of their coefficients. The theorem above gives an explicit factorization of $P_N$ into two degree- $(n-1)$ squares with strictly positive coefficients for every odd $N$ ; this positive factorization appears novel. Martin and O’Bryant (2009) and Matolcsi and Vinuesa (2010) study the corresponding one-source extremal problem.

1. Problem Definitions

2. Exact verification of the 16,38416{,}38416,384-step construction

3. Additional findings

3.1 Two isolated local minima with identical support and contact sets

3.2 The perfectly flat target requires exactly two positive sources

2. Exact verification of the $16{,}384$ -step construction