Prime Number Theorem: Verification Notebook

This notebook is a self-contained mathematical and computational companion to the Prime number theorem 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. S1. A score of $0.980681$ valid for every $x$ .
3. S2. Why a direct Möbius cutoff fails.
4. Additional findings

The required packages are:

NumPy
GCC with OpenMP support for the complete-period scan

1. Problem Definitions

The prime number theorem describes the asymptotic density of the primes. This notebook studies a finite Chebyshev certificate that turns a globally valid floor-sum inequality into an explicit lower bound for that density. The original benchmark considers both lower and upper certificates; this notebook concerns its lower-bound side. If $\pi(x)$ counts the primes at most $x$ , the theorem states that

\lim_{x\to\infty}\frac{\pi(x)}{x/\log x}=1.

Let $f$ be a finitely supported real function on the positive integers such that

\sum_k \frac{f(k)}{k}=0.

Define

F_f(x)=\sum_k f(k)\left\lfloor\frac{x}{k}\right\rfloor, \qquad A(f)=-\sum_k\frac{f(k)\log k}{k}.

Because every jump of $F_f$ occurs at an integer, checking $F_f(x)\leq1$ at every positive integer is equivalent to checking it for every real $x\geq1$ .

The next cell imports the libraries used throughout the notebook.

Show code

Code cell 3 · In [1]

from pathlib import Path
from decimal import Decimal, ROUND_CEILING, ROUND_FLOOR, localcontext
from fractions import Fraction
import ctypes
import hashlib
import itertools
import math
import os
import subprocess
import tempfile
import time

import numpy as np

2. S1. A score of $0.980681$ valid for every $x$ .

2.1 Prime-counting consequence

The next lemma is the classical finite Chebyshev argument; see Diamond (1982).

Lemma 2.1 (finite Chebyshev weight). If $F_f(x)\le 1$ for every real $x\ge 1$ , then

\liminf_{x\to\infty}\frac{\pi(x)}{x/\log x}\ge A(f).

Proof. Write $\Lambda$ for the von Mangoldt function and $\psi(x)=\sum_{m\le x}\Lambda(m)$ . The identity

\log(N!)=\sum_{m\le N}\Lambda(m)\left\lfloor\frac{N}{m}\right\rfloor

gives, after interchanging two finite sums,

\sum_k f(k)\log\!\left(\left\lfloor\frac{x}{k}\right\rfloor!\right) =\sum_{m\le x}\Lambda(m)F_f(x/m) \le \psi(x).

Stirling’s formula and the balance condition show that the left side equals $xA(f)+O_f(\log x)$ : the terms proportional to $x\log x$ and to $x$ cancel. Hence

\liminf_{x\to\infty}\frac{\psi(x)}x\ge A(f).

Prime powers contribute only $O(\sqrt{x}\log x)$ to $\psi(x)-\vartheta(x)$ , while $\vartheta(x)\le\pi(x)\log x$ . It follows that

\pi(x)\ge\frac{\psi(x)-O(\sqrt{x}\log x)}{\log x},

which proves the lemma. $\square$

2.2 A 255-term all- $x$ certificate

Theorem 2.2 (finite all- $x$ construction). There exists a rational function $f$ , supported on $255$ positive integers, for which

\sum_k\frac{f(k)}k=0, \qquad F_f(x)\le1\quad(x\ge1),

and

0.9806811331448391184649165900079202241771553360836383429397186 <A(f)< 0.9806811331448391184649165900079202241771553360836383429397187.

Verification. The bundled data encode one such function. They store the nonanchor coefficients as integers $b_k$ on the grid $Q=10^{12}$ :

f(k)=\frac{b_k}{Q}\quad(k>1), \qquad f(1)=-\sum_{k>1}\frac{b_k}{Qk}.

The denominator $1$ is called the anchor because its value is determined by the balance equation rather than independently chosen.

The next cell authenticates the bundled coefficients, reconstructs the anchor coefficient, and checks exact reciprocal balance.

Show code

Code cell 5 · In [2]

DATA_PATH = Path("prime_number_theorem_construction.npz")
EXPECTED_SHA256 = "391b92f93ce0e1a3f65bfbee09009f9f043fac97f908dafad25da1cf1c4aeb0b"
actual_sha256 = hashlib.sha256(DATA_PATH.read_bytes()).hexdigest()
assert actual_sha256 == EXPECTED_SHA256

with np.load(DATA_PATH, allow_pickle=False) as data:
    keys = data["nonanchor_keys"].astype(np.int64)
    b_values = data["grid_numerators"].astype(np.int64)
    Q = int(data["grid"])
    PERIOD = int(data["period"])

assert len(keys) == len(b_values) == 255
assert len(set(map(int, keys))) == 255
assert np.all(keys > 1)
assert all(PERIOD % int(k) == 0 for k in keys)

anchor_numerator = -sum(
    int(b) * (PERIOD // int(k)) for k, b in zip(keys, b_values)
)
COMMON_DENOMINATOR = Q * PERIOD
common_numerators = {1: anchor_numerator}
common_numerators.update(
    {int(k): int(b) * PERIOD for k, b in zip(keys, b_values)}
)

balance_numerator = anchor_numerator + sum(
    common_numerators[int(k)] // int(k) for k in keys
)
assert balance_numerator == 0
assert all(-10 * Q <= int(b) <= 10 * Q for b in b_values)

print("data SHA-256:", actual_sha256)
support_size = int(np.count_nonzero(b_values)) + int(anchor_numerator != 0)
assert support_size == 255
print("listed denominators, including the anchor:", len(keys) + 1)
print("nonzero support size:", support_size)
print("common period:", f"{PERIOD:,}")
print("exact reciprocal-balance numerator:", balance_numerator)
print("coefficient grid: 1 /", Q)
print("\nFirst twelve nonanchor coefficients b_k/Q:")
for k, b in list(zip(keys, b_values))[:12]:
    print(f"  k={int(k):>4}: {int(b)}/{Q}")

Saved output 1

data SHA-256: 391b92f93ce0e1a3f65bfbee09009f9f043fac97f908dafad25da1cf1c4aeb0b
listed denominators, including the anchor: 256
nonzero support size: 255
common period: 5,354,228,880
exact reciprocal-balance numerator: 0
coefficient grid: 1 / 1000000000000

First twelve nonanchor coefficients b_k/Q:
  k=   2: -999999645495/1000000000000
  k=   3: -999999645495/1000000000000
  k=   5: -999999645495/1000000000000
  k=   6: 999999645495/1000000000000
  k=   7: -999999645495/1000000000000
  k=  10: 999999645495/1000000000000
  k=  11: -999999645495/1000000000000
  k=  13: -999999645495/1000000000000
  k=  14: 999999645495/1000000000000
  k=  15: 999999645495/1000000000000
  k=  17: -999999645495/1000000000000
  k=  19: -999999645495/1000000000000

Complete-period verification

Every listed denominator divides

L=5{,}354{,}228{,}880.

The balance equation therefore gives

F_f(x+L)-F_f(x) =\sum_k f(k)\frac{L}{k} =L\sum_k\frac{f(k)}k=0.

Thus $F_f$ is periodic, and it suffices to check the integer points $1,\ldots,L$ . The increment at $n$ is

F_f(n)-F_f(n-1)=\sum_{k\mid n}f(k).

Since every $k$ divides $L$ , the increment depends only on $\gcd(n,L)$ . The factorization

L=2^4 3^2\cdot5\cdot7\cdot11\cdot13\cdot17\cdot19\cdot23

has only

\tau(L)=5\cdot3\cdot2^7=1920

divisors. The next cell constructs the $1{,}920$ divisor-state increments used by the complete-period scan.

Show code

Code cell 7 · In [3]

def factor_integer(value):
    remaining = int(value)
    factors = []
    p = 2
    while p * p <= remaining:
        exponent = 0
        while remaining % p == 0:
            remaining //= p
            exponent += 1
        if exponent:
            factors.append((p, exponent))
        p = 3 if p == 2 else p + 2
    if remaining > 1:
        factors.append((remaining, 1))
    return factors


def state_description(factors):
    radices = [exponent + 1 for _, exponent in factors]
    strides = []
    total = 1
    for radix in radices:
        strides.append(total)
        total *= radix
    divisors = np.ones(total, dtype=np.int64)
    for state in range(total):
        rest = state
        value = 1
        for (prime, _), radix in zip(factors, radices):
            value *= prime ** (rest % radix)
            rest //= radix
        divisors[state] = value
    return strides, divisors


def gcd_state_codes(start, stop, factors, strides):
    codes = np.zeros(stop - start + 1, dtype=np.uint16)
    for (prime, exponent), stride in zip(factors, strides):
        power = 1
        for _ in range(exponent):
            power *= prime
            first = ((start + power - 1) // power) * power
            if first <= stop:
                codes[first - start :: power] += np.uint16(stride)
    return codes


factors = factor_integer(PERIOD)
strides, state_divisors = state_description(factors)
assert factors == [(2, 4), (3, 2), (5, 1), (7, 1), (11, 1),
                   (13, 1), (17, 1), (19, 1), (23, 1)]
assert len(state_divisors) == 1920

items = tuple(common_numerators.items())
increment_table = [
    sum(numerator for key, numerator in items if int(divisor) % key == 0)
    for divisor in state_divisors
]

LIMB_BASE = 10**18
table_high = np.asarray([value // LIMB_BASE for value in increment_table], dtype=np.int64)
table_low = np.asarray([value % LIMB_BASE for value in increment_table], dtype=np.uint64)
wheel_codes = np.ascontiguousarray(
    gcd_state_codes(1, 720_720, factors[:6], strides[:6]), dtype=np.uint16
)
assert len(wheel_codes) == 720_720
print("factorization:", factors)
print("gcd states:", len(state_divisors))
print("wheel states prepared:", len(wheel_codes))

Saved output 1

factorization: [(2, 4), (3, 2), (5, 1), (7, 1), (11, 1), (13, 1), (17, 1), (19, 1), (23, 1)]
gcd states: 1920
wheel states prepared: 720720

The next cell compiles the exact scanner, traverses the complete period, and returns the maximum, minimum, and endpoint floor-sum numerators.

Show code

Code cell 9 · In [4]

SCANNER_SOURCE = r'''
#include <stdint.h>
#include <stdlib.h>
#ifdef _OPENMP
#include <omp.h>
#endif

typedef __int128 i128;
#define WHEEL 720720U
#define STATES 1920U
#define LIMB_BASE ((i128)1000000000000000000LL)

typedef struct {
    i128 total, maximum, minimum;
    uint64_t argmax, argmin;
} block_stats;

static inline i128 join_limbs(int64_t high, uint64_t low) {
    return ((i128)high) * LIMB_BASE + (i128)low;
}

static inline void split_limbs(i128 value, int64_t *high, uint64_t *low) {
    i128 q = value / LIMB_BASE;
    i128 r = value % LIMB_BASE;
    if (r < 0) { q -= 1; r += LIMB_BASE; }
    *high = (int64_t)q;
    *low = (uint64_t)r;
}

int exact_prefix_scan(
    const int64_t *table_high,
    const uint64_t *table_low,
    const uint16_t *wheel_codes,
    uint64_t period,
    uint64_t block_size,
    int threads,
    int64_t *out_high,
    uint64_t *out_low,
    uint64_t *out_n
) {
    if (block_size < 1) return 1;
    uint64_t blocks = (period + block_size - 1) / block_size;
    block_stats *stats = (block_stats *)calloc((size_t)blocks, sizeof(block_stats));
    if (!stats) return 2;

#ifdef _OPENMP
    if (threads > 0) omp_set_num_threads(threads);
#endif
#pragma omp parallel for schedule(static)
    for (uint64_t block = 0; block < blocks; ++block) {
        uint64_t start = block * block_size + 1;
        uint64_t stop = start + block_size - 1;
        if (stop > period) stop = period;
        uint32_t wi = (uint32_t)((start - 1) % WHEEL);
        uint32_t c17 = (uint32_t)((start - 1) % 17U);
        uint32_t c19 = (uint32_t)((start - 1) % 19U);
        uint32_t c23 = (uint32_t)((start - 1) % 23U);
        i128 running = 0, best = 0, worst = 0;
        uint64_t best_n = 0, worst_n = 0;
        int initialized = 0;

        for (uint64_t n = start; n <= stop; ++n) {
            uint16_t state = wheel_codes[wi];
            if (++c17 == 17U) { state += 240U; c17 = 0; }
            if (++c19 == 19U) { state += 480U; c19 = 0; }
            if (++c23 == 23U) { state += 960U; c23 = 0; }
            if (++wi == WHEEL) wi = 0;
            running += join_limbs(table_high[state], table_low[state]);
            if (!initialized || running > best) { best = running; best_n = n; }
            if (!initialized || running < worst) { worst = running; worst_n = n; }
            initialized = 1;
        }
        stats[block].total = running;
        stats[block].maximum = best;
        stats[block].minimum = worst;
        stats[block].argmax = best_n;
        stats[block].argmin = worst_n;
    }

    i128 offset = 0, global_max = 0, global_min = 0;
    uint64_t global_argmax = 0, global_argmin = 0;
    int initialized = 0;
    for (uint64_t block = 0; block < blocks; ++block) {
        i128 candidate_max = offset + stats[block].maximum;
        i128 candidate_min = offset + stats[block].minimum;
        if (!initialized || candidate_max > global_max) {
            global_max = candidate_max; global_argmax = stats[block].argmax;
        }
        if (!initialized || candidate_min < global_min) {
            global_min = candidate_min; global_argmin = stats[block].argmin;
        }
        initialized = 1;
        offset += stats[block].total;
    }
    split_limbs(global_max, out_high, out_low);
    split_limbs(global_min, out_high + 1, out_low + 1);
    split_limbs(offset, out_high + 2, out_low + 2);
    out_n[0] = global_argmax;
    out_n[1] = global_argmin;
    free(stats);
    return 0;
}
'''


def join_limbs(high, low):
    return int(high) * LIMB_BASE + int(low)


def exhaustive_exact_scan(block_size=5_000_000, threads=None):
    thread_count = int(threads or min(48, os.cpu_count() or 1))
    with tempfile.TemporaryDirectory(prefix="pnt_exact_scan_") as temp_dir:
        source_path = Path(temp_dir) / "exact_scan.c"
        library_path = Path(temp_dir) / "exact_scan.so"
        source_path.write_text(SCANNER_SOURCE)
        subprocess.run(
            ["gcc", "-O3", "-march=native", "-fopenmp", "-shared", "-fPIC",
             str(source_path), "-o", str(library_path)],
            check=True,
        )
        library = ctypes.CDLL(str(library_path))
        function = library.exact_prefix_scan
        function.argtypes = [
            ctypes.POINTER(ctypes.c_int64),
            ctypes.POINTER(ctypes.c_uint64),
            ctypes.POINTER(ctypes.c_uint16),
            ctypes.c_uint64,
            ctypes.c_uint64,
            ctypes.c_int,
            ctypes.POINTER(ctypes.c_int64),
            ctypes.POINTER(ctypes.c_uint64),
            ctypes.POINTER(ctypes.c_uint64),
        ]
        function.restype = ctypes.c_int
        out_high = np.zeros(3, dtype=np.int64)
        out_low = np.zeros(3, dtype=np.uint64)
        out_n = np.zeros(2, dtype=np.uint64)
        started = time.perf_counter()
        status = function(
            table_high.ctypes.data_as(ctypes.POINTER(ctypes.c_int64)),
            table_low.ctypes.data_as(ctypes.POINTER(ctypes.c_uint64)),
            wheel_codes.ctypes.data_as(ctypes.POINTER(ctypes.c_uint16)),
            PERIOD,
            block_size,
            thread_count,
            out_high.ctypes.data_as(ctypes.POINTER(ctypes.c_int64)),
            out_low.ctypes.data_as(ctypes.POINTER(ctypes.c_uint64)),
            out_n.ctypes.data_as(ctypes.POINTER(ctypes.c_uint64)),
        )
        elapsed = time.perf_counter() - started
    if status != 0:
        raise RuntimeError(f"exact scanner returned status {status}")
    return {
        "maximum": join_limbs(out_high[0], out_low[0]),
        "minimum": join_limbs(out_high[1], out_low[1]),
        "final": join_limbs(out_high[2], out_low[2]),
        "argmax": int(out_n[0]),
        "argmin": int(out_n[1]),
        "threads": thread_count,
        "elapsed": elapsed,
    }


scan = exhaustive_exact_scan()
EXPECTED_MAXIMUM = 5_354_226_994_850_107_576_341
EXPECTED_ARGMAX = 210_039_343
EXPECTED_MINIMUM = -475_454_681_777_932_352_262_768
EXPECTED_ARGMIN = 5_144_189_536

assert scan["maximum"] == EXPECTED_MAXIMUM
assert scan["argmax"] == EXPECTED_ARGMAX
assert scan["minimum"] == EXPECTED_MINIMUM
assert scan["argmin"] == EXPECTED_ARGMIN
assert scan["final"] == 0
assert scan["maximum"] < COMMON_DENOMINATOR

direct_at_argmax = sum(
    numerator * (scan["argmax"] // key)
    for key, numerator in common_numerators.items()
)
direct_at_endpoint = sum(
    numerator * (PERIOD // key)
    for key, numerator in common_numerators.items()
)
assert direct_at_argmax == scan["maximum"]
assert direct_at_endpoint == 0

print("complete-period rows checked:", f"{PERIOD:,}")
print("worker threads:", scan["threads"])
print("scan time (seconds):", f"{scan['elapsed']:.3f}")
print("exact maximum numerator:", scan["maximum"])
print("exact maximum location:", scan["argmax"])
print("roof denominator:", COMMON_DENOMINATOR)
print("exact roof margin:", COMMON_DENOMINATOR - scan["maximum"])
print("period endpoint numerator:", scan["final"])

Saved output 1

complete-period rows checked: 5,354,228,880
worker threads: 48
scan time (seconds): 0.468
exact maximum numerator: 5354226994850107576341
exact maximum location: 210039343
roof denominator: 5354228880000000000000
exact roof margin: 1885149892423659
period endpoint numerator: 0

Logarithmic objective

The next cell uses directed decimal rounding to enclose $A(f)$ in a rigorous interval.

Show code

Code cell 11 · In [5]

def objective_interval(precision=80):
    lower = Decimal(0)
    upper = Decimal(0)
    widest_log_interval = Decimal(0)
    for key, b in zip(map(int, keys), map(int, b_values)):
        if b == 0:
            continue
        with localcontext() as context:
            context.prec = precision
            rounded_log = Decimal(key).ln(context=context)
            log_lower = rounded_log.next_minus(context=context)
            log_upper = rounded_log.next_plus(context=context)
        widest_log_interval = max(widest_log_interval, log_upper - log_lower)

        # The coefficient of log(k) is -b/(Q*k).
        # Choose endpoints according to the coefficient sign so the sums
        # remain rigorous lower and upper bounds.
        lower_log = log_lower if b < 0 else log_upper
        upper_log = log_upper if b < 0 else log_lower
        denominator = Decimal(Q * key)
        with localcontext() as context:
            context.prec = precision
            context.rounding = ROUND_FLOOR
            lower += Decimal(-b) * lower_log / denominator
        with localcontext() as context:
            context.prec = precision
            context.rounding = ROUND_CEILING
            upper += Decimal(-b) * upper_log / denominator
    return lower, upper, widest_log_interval


A_lower, A_upper, widest_log_interval = objective_interval()
TARGET = Decimal("0.980681")
EXPECTED_LOWER_PREFIX = Decimal(
    "0.9806811331448391184649165900079202241771553360836383429397186"
)
EXPECTED_UPPER_PREFIX = Decimal(
    "0.9806811331448391184649165900079202241771553360836383429397187"
)
assert A_lower > EXPECTED_LOWER_PREFIX
assert A_upper < EXPECTED_UPPER_PREFIX
assert A_lower > TARGET
assert A_lower <= A_upper

print("A lower endpoint:", A_lower)
print("A upper endpoint:", A_upper)
print("interval width:", A_upper - A_lower)
print("widest individual log enclosure:", widest_log_interval)
print("certified A(f) > 0.980681:", A_lower > TARGET)

Saved output 1

A lower endpoint: 0.98068113314483911846491659000792022417715533608363834293971860151222815392876850
A upper endpoint: 0.98068113314483911846491659000792022417715533608363834293971860151222815392878068
interval width: 1.218E-77
widest individual log enclosure: 2E-78
certified A(f) > 0.980681: True

Proof of Theorem 2.2. Balance makes $F_f$ periodic, the complete-period scan proves $F_f(x)\leq1$ , and the interval calculation proves the displayed score. $\square$

Station Reference. The construction and proof method come from Archive Room paper #32 and Research Center Evals #553, #578, #581, and #598. Eval #581 supplied the complete-period scan, and Eval #598 supplied the directed-rounding objective check.

Related Work. Diamond (1982) surveys the broader Chebyshev method and records Sylvester’s classical lower bound of approximately $0.956$ . Diamond and Erdős (1980) showed that sufficiently long finite Möbius data yield Chebyshev-type estimates arbitrarily close to $1$ . Georgiev et al. (2025) introduced this finite benchmark and reported $0.938$ . The explicit 255-term bounded rational all- $x$ certificate appears novel within the benchmark’s direct-roof class.

3. S2. Why a direct Möbius cutoff fails.

Let $\mu$ be the Möbius function, let $M(D)=\sum_{d\le D}\mu(d)$ , and put

g_D(n)=\sum_{2\le d\le D}\mu(d) \left(\left\{\frac nd\right\}-\frac12\right), \qquad L_D=\operatorname{lcm}(1,\ldots,D).

The balanced direct cutoff uses $f_D(d)=\mu(d)$ for $2\le d\le D$ and chooses $f_D(1)$ from $\sum f_D(d)/d=0$ .

Theorem 3.1 (Möbius-cutoff obstruction). If

R_D=\max_n F_{f_D}(n),

then

R_D=\Omega\!\left(\frac{D}{\log^2D}\right).

Moreover the unscaled objective $A(f_D)$ tends to $1$ . Consequently, dividing $f_D$ by the factor needed to impose the roof $F\le1$ gives score

O\!\left(\frac{\log^2D}{D}\right),

which tends to zero.

Proof. Define the incomplete Möbius divisor sum

E_D(n)=\sum_{\substack{d\mid n\\d\le D}}\mu(d).

For every $d\le D$ , we have $d\mid L_D$ and

\left\{\frac nd\right\}+\left\{\frac{L_D-n}{d}\right\} =1-\mathbf 1_{d\mid n}.

After subtracting the two half-unit centerings and summing, this gives the exact reflection identity

g_D(n)+g_D(L_D-n)=1-E_D(n).

If $O(D)=\max_m(-g_D(m))$ , at least one member of every reflected pair is no greater than their average, so

O(D)\ge\frac{\max_nE_D(n)-1}{2}. \tag{1}

Now let $n_*$ be the product of all primes in

D^{1/3}<p\le D^{1/2},

and let $r$ be the number of these primes. A divisor of $n_*$ that is at most $D$ uses exactly zero, one, or two band primes: every two-prime product is at most $D$ , while every three-prime product exceeds $D$ . Hence

E_D(n_*)=1-r+\binom r2=\frac{(r-1)(r-2)}2. \tag{2}

The standard estimates $\pi(x)=\Theta(x/\log x)$ give $r=\Theta(\sqrt D/\log D)$ . Equations (1)–(2) therefore imply

O(D)=\Omega(D/\log^2D). \tag{3}

The balanced cutoff path is

F_{f_D}(n)=-\sum_{d=2}^D\mu(d)\left\{\frac nd\right\} =-g_D(n)-\frac{M(D)-1}{2}. \tag{4}

A classical zero-free-region estimate for the zeta function gives

M(D)=o(D/\log^2D).

Combining this with (3)–(4) proves the asserted lower bound for $R_D$ .

Finally,

A(f_D)=-\sum_{d\le D}\frac{\mu(d)\log d}{d}\longrightarrow1.

Indeed, for $s>1$ ,

-\sum_{d\ge1}\frac{\mu(d)\log d}{d^s} =\left(\frac1{\zeta(s)}\right)',

and the right side tends to $1$ as $s\downarrow1$ because $1/\zeta(s)\sim s-1$ . The same zero-free-region estimate just used supplies convergence at $s=1$ . Thus $A(f_D)=1+o(1)$ , whereas the roof rescaling is at least a constant multiple of $D/\log^2D$ . This proves the score bound. $\square$

The next cell checks the reflection identity and prime-band witness for small cutoffs using rational arithmetic.

Show code

Code cell 14 · In [6]

def mobius_sieve(limit):
    mu = [1] * (limit + 1)
    is_prime = [True] * (limit + 1)
    primes = []
    mu[0] = 0
    for n in range(2, limit + 1):
        if is_prime[n]:
            primes.append(n)
            mu[n] = -1
        for p in primes:
            if n * p > limit:
                break
            is_prime[n * p] = False
            if n % p == 0:
                mu[n * p] = 0
                break
            mu[n * p] = -mu[n]
    return mu, primes


def lcm_upto(D):
    value = 1
    for n in range(1, D + 1):
        value = math.lcm(value, n)
    return value


def centered_g(D, n, mu):
    return sum(
        Fraction(mu[d]) * (Fraction(n % d, d) - Fraction(1, 2))
        for d in range(2, D + 1)
    )


mu_small, _ = mobius_sieve(12)
reflection_checks = 0
for D in (6, 8, 10, 12):
    L = lcm_upto(D)
    for n in range(L):
        E = sum(mu_small[d] for d in range(1, D + 1) if n % d == 0)
        assert centered_g(D, n, mu_small) + centered_g(D, L - n, mu_small) == 1 - E
        reflection_checks += 1

mu_band, primes = mobius_sieve(1_000_000)
band_rows = []
for D in (1_000, 10_000, 100_000, 1_000_000):
    band = [p for p in primes if p**3 > D and p**2 <= D]
    assert all(p * q <= D for i, p in enumerate(band) for q in band[i + 1:])
    assert all(
        band[i] * band[j] * band[k] > D
        for i in range(len(band))
        for j in range(i + 1, len(band))
        for k in range(j + 1, len(band))
    )
    r = len(band)
    E_closed = (r - 1) * (r - 2) // 2
    E_subset = 1 - r + r * (r - 1) // 2
    assert E_closed == E_subset
    band_rows.append((D, r, E_closed, Fraction(E_closed - 1, 2)))

print("exact reflection identities checked:", f"{reflection_checks:,}")
print("\nPrime-band finite witnesses:")
print("       D   band primes       E_D(n*)      reflection lower bound")
for D, r, E, lower_bound in band_rows:
    print(f"{D:>8,} {r:>13} {E:>15,} {str(lower_bound):>27}")

Saved output 1

exact reflection identities checked: 31,140

Prime-band finite witnesses:
       D   band primes       E_D(n*)      reflection lower bound
   1,000             7              15                           7
  10,000            17             120                       119/2
 100,000            51           1,225                         612
1,000,000           143          10,011                        5005

Station Reference. The reflection identity, prime-band witness, and cutoff obstruction come from Archive Room paper #18, building on paper #13. Supporting records are Research Center Evals #247, #265, #270, #272, #281, and #282 and Question Room discussion #4.

Related Work. Letendre (2020) studies the same incomplete Möbius convolution, and the prime-band extremal mechanism is classical. Esser, Tao, Totaro, and Wang (2022) study one-sided sawtooth inequalities in a different coefficient class. The reflection bridge from incomplete divisor sums to the direct-cutoff obstruction appears novel.

4. Additional findings

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

4.1 Prime-power-only constructions stop at $\log2$

We prove that every balanced certificate supported only on $1$ and prime powers satisfies the sharp bound

A(f)\le\log2,

with equality for $f(1)=1$ and $f(2)=-2$ . The Chinese remainder theorem separates the contributions from distinct primes, while an exact digit-moment argument gives the optimal contribution of each prime-power tower. Thus adding depth within isolated prime towers cannot approach the unit score; the $0.980681$ construction succeeds because its composite denominators couple different primes.

Call a support prime-power-only if every denominator other than $1$ is $p^j$ for some prime $p$ .

Theorem 4.1 (sharp prime-power ceiling). Every balanced prime-power-only construction satisfying $F_f(n)\le1$ for all $n$ obeys

A(f)\le\log2.

Equality is attained by $f(1)=1$ , $f(2)=-2$ .

Proof. For each prime $p$ , let $J_p$ be the largest exponent used and write

B_p(r)=-\sum_{j=1}^{J_p}f(p^j)\left\{\frac{r}{p^j}\right\}.

Balance gives $F_f(n)=\sum_pB_p(n)$ . The prime-power moduli $p^{J_p}$ are pairwise coprime, so the Chinese remainder theorem lets their residues vary independently. Therefore

\max_nF_f(n)=\sum_p c_p, \qquad c_p:=\max_rB_p(r)\ge0. \tag{5}

We first prove the sharp single-prime inequality

A_p:=-\sum_{j=1}^{J_p}\frac{f(p^j)j\log p}{p^j} \le c_p\frac{\log p}{p-1}. \tag{6}

Put $a_j=-f(p^j)$ and write a residue $r\pmod {p^J}$ in base $p$ as $r=\sum_{i=0}^{J-1}d_ip^i$ . Let $\lambda$ be $\log p/(p-1)$ times the following probability measure on digit strings:

$d_0=p-1$ with probability $1$ ;
independently for $i\ge1$ , $d_i=1$ with probability $(p-1)/p^i$ and $d_i=0$ otherwise.

Then $\lambda$ is nonnegative, has total mass $\log p/(p-1)$ , and

\int d_i\,d\lambda=\frac{\log p}{p^i}.

Since

\left\{\frac r{p^j}\right\}=\sum_{i<j}d_ip^{i-j},

we obtain the exact moment identity

\int\left\{\frac r{p^j}\right\}d\lambda =\sum_{i<j}p^{i-j}\frac{\log p}{p^i} =\frac{j\log p}{p^j}.

Consequently,

A_p=\int B_p(r)\,d\lambda(r) \le c_p\frac{\log p}{p-1},

which proves (6). It is sharp: using only $p$ with $-f(p)=c_pp/(p-1)$ makes $\max B_p=c_p$ and attains equality.

Finally, the roof and (5) give $\sum_pc_p\le1$ . The function $\log x/(x-1)$ decreases for $x>1$ , so its largest value at a prime is $\log2$ . Summing (6) yields

A(f)\le\sum_pc_p\frac{\log p}{p-1}\le\log2.

For $f(1)=1$ , $f(2)=-2$ , the path equals $1$ at odd integers and $0$ at even integers, and $A(f)=\log2$ . $\square$

The result shows why the successful $0.980681$ construction needs denominators that couple several primes: adding depth inside isolated prime towers cannot pass even the elementary binary value.

The next cell checks the base- $p$ digit moments used in the single-prime dual argument.

Show code

Code cell 18 · In [7]

def digit_dual_moments(p, J):
    probabilities = []
    for digits in itertools.product(range(p), repeat=J):
        if digits[0] != p - 1:
            continue
        probability = Fraction(1)
        for i in range(1, J):
            success = Fraction(p - 1, p**i)
            if digits[i] == 1:
                probability *= success
            elif digits[i] == 0:
                probability *= 1 - success
            else:
                probability = Fraction(0)
                break
        if probability:
            r = sum(digit * p**i for i, digit in enumerate(digits))
            probabilities.append((r, probability))
    assert sum(probability for _, probability in probabilities) == 1
    moments_without_log = []
    for j in range(1, J + 1):
        moment = sum(
            probability * Fraction(r % (p**j), p**j)
            for r, probability in probabilities
        ) / (p - 1)
        moments_without_log.append(moment)
    return moments_without_log


dual_cases = []
for p in (2, 3, 5, 7):
    for J in range(1, 6):
        moments = digit_dual_moments(p, J)
        expected = [Fraction(j, p**j) for j in range(1, J + 1)]
        assert moments == expected
        dual_cases.append((p, J))

print("exact digit-dual moment systems verified:", len(dual_cases))
print("cases: p in {2,3,5,7}, depth J in {1,2,3,4,5}")

Saved output 1

exact digit-dual moment systems verified: 20
cases: p in {2,3,5,7}, depth J in {1,2,3,4,5}

Station Reference. The prime-power decomposition, single-prime frontier, and $\log2$ ceiling come from Archive Room paper #31 and Research Center Evals #565 and #589.

Related Work. The two-term construction $f(1)=1$ , $f(2)=-2$ , attaining $A(f)=\log 2$ , is a classical example recorded by Diamond (1982). Esser, Tao, Totaro, and Wang (2022) solve an adjacent one-sided sawtooth optimization problem with a different coefficient and support class. The sharp upper bound for every prime-power-only construction in the signed Chebyshev certificate class appears novel.

1. Problem Definitions

2. S1. A score of 0.9806810.9806810.980681 valid for every xxx.

2.1 Prime-counting consequence

2.2 A 255-term all-xxx certificate

Complete-period verification

Logarithmic objective

3. S2. Why a direct Möbius cutoff fails.

4. Additional findings

4.1 Prime-power-only constructions stop at log⁡2\log2log2

2. S1. A score of $0.980681$ valid for every $x$ .

2.2 A 255-term all- $x$ certificate

4.1 Prime-power-only constructions stop at $\log2$