1. constructs every \(P(n)\) through \(10^7+1\) by overwriting multiples of

increasing primes;

2. independently reconstructs every \(P(n)\) using a smallest-prime-factor sieve

and the recurrence

\(P(n)=\max(\operatorname{spf}(n),P(n/\operatorname{spf}(n)))\);

3. checks \(P(n)\) a third way by trial division for \(1\leq n\leq20000\);

4. checks hard-coded exact joint and marginal counts.

All checks passed. The final run used Python 3 standard library only, took 28.952

seconds, and peaked at roughly 130 MiB.

3.2 Exact counts

Each entry is the number of \(n\) with \(2\leq n\leq X\) satisfying the two

strict live-page inequalities.

| \(X\) | \((1/3,1/3)\) | \((1/3,1/2)\) | \((1/2,1/3)\) | \((1/2,1/2)\) | \((1/2,2/3)\) | \((2/3,1/2)\) | \((2/3,2/3)\) |

|---:|---:|---:|---:|---:|---:|---:|---:|

| \(10^3\) | 0 | 10 | 12 | 36 | 114 | 110 | 242 |

| \(10^4\) | 8 | 86 | 99 | 492 | 1,226 | 1,204 | 2,663 |

| \(10^5\) | 86 | 891 | 926 | 5,765 | 13,198 | 13,182 | 28,225 |

| \(10^6\) | 1,053 | 9,208 | 9,347 | 63,102 | 138,909 | 139,020 | 294,981 |

| \(10^7\) | 11,617 | 96,970 | 97,497 | 675,882 | 1,451,554 | 1,451,271 | 3,040,588 |

At \(X=10^7\), let \(J\) be the joint count, \(L\) and \(R\) the two finite

marginal counts on the common \(X-1\)-element sample, and define the

finite-marginal independence ratio \(J(X-1)/(LR)\). The exact marginal counts

for exponents \(1/3,1/2,2/3\) are respectively

on both sides for this endpoint.

| \((\alpha,\beta)\) | \(J/X\) | \(J(X-1)/(LR)\) | \(\rho(1/\alpha)\rho(1/\beta)\) |

|---:|---:|---:|---:|

| \((1/3,1/3)\) | 0.001161700000 | 0.608162689029 | 0.002362775412 |

| \((1/3,1/2)\) | 0.009697000000 | 0.815915566645 | 0.014915620996 |

| \((1/2,1/3)\) | 0.009749700000 | 0.820349798919 | 0.014915620996 |

| \((1/2,1/2)\) | 0.067588200000 | 0.914030618552 | 0.094158652798 |

| \((1/2,2/3)\) | 0.145155400000 | 0.958106973096 | 0.182434707832 |

| \((2/3,1/2)\) | 0.145127100000 | 0.957920177239 | 0.182434707832 |

| \((2/3,2/3)\) | 0.304058800000 | 0.979556118706 | 0.353471737677 |

These figures are exact counts followed by displayed decimal evaluations. They do

not establish monotonicity, convergence, or a counterexample. In particular,

convergence of even the one-variable smooth-number density is slow at these

sizes.

3.3 Computational ceiling

(c) The present method is \(O(X\log\log X)\) time and \(O(X)\) memory. A

literal Python extrapolation to \(10^9\) is about one core-hour and over 12 GiB;

a segmented compiled implementation would reduce memory and constants but still

only produce another finite prefix. No finite cutoff can remove the exceptional

scales in (2) or prove the uniform \(o(X)\) statement (11). The needed advance is

analytic, not a larger census.

4. Reproduction

Run:

python3 runs/erdos928_wave7t_reverify.py

Expected final lines include:

independent trial-division prefix 1..20000: OK
large-factor inclusion-exclusion/prime-quadruple prefix audit: OK
independent full smallest-factor reconstruction: OK
ALL CHECKS PASSED

The complete standalone source is in runs/erdos928_wave7t_reverify.py and is

reproduced below.

#!/usr/bin/env python3
"""
From-scratch exact finite verifier for Erdős problem 928.

For rational gamma=a/b the live problem's strict condition

    P(n) < n**gamma

is tested without floating point as

    P(n)**b < n**a.

The program builds the largest-prime-factor table by an Eratosthenes sieve,
checks the entire table against an independent smallest-prime-factor sieve
and recurrence, checks a prefix a third way by trial division, counts seven
(alpha,beta) pairs at the advertised cutoffs, and verifies the hard-coded
exact census.
Only the Python standard library is used.

Default cost (limit 10^7): about 130 MiB peak RAM and roughly half a minute on
the exe.dev VM used for the report.  A smaller exploratory run is available
with --limit, but the advertised-census assertions run only at covered
cutoffs.
"""

from __future__ import annotations

import argparse
from array import array
from dataclasses import dataclass
import gc
from math import log
from time import perf_counter


@dataclass(frozen=True, order=True)
class Exponent:
    numerator: int
    denominator: int

    def label(self) -> str:
        return f"{self.numerator}/{self.denominator}"


ONE_THIRD = Exponent(1, 3)
ONE_HALF = Exponent(1, 2)
TWO_THIRDS = Exponent(2, 3)
EXPONENTS = (ONE_THIRD, ONE_HALF, TWO_THIRDS)

PAIRS = (
    (ONE_THIRD, ONE_THIRD),
    (ONE_THIRD, ONE_HALF),
    (ONE_HALF, ONE_THIRD),
    (ONE_HALF, ONE_HALF),
    (ONE_HALF, TWO_THIRDS),
    (TWO_THIRDS, ONE_HALF),
    (TWO_THIRDS, TWO_THIRDS),
)

DEFAULT_CHECKPOINTS = (10**3, 10**4, 10**5, 10**6, 10**7)

# Filled from one run and then checked by every fresh run.  These are counts
# of n with 2 <= n <= X satisfying both strict inequalities.
EXPECTED_COUNTS: dict[int, tuple[int, ...]] = {
    # X: counts in the order PAIRS above
    10**3: (0, 10, 12, 36, 114, 110, 242),
    10**4: (8, 86, 99, 492, 1226, 1204, 2663),
    10**5: (86, 891, 926, 5765, 13198, 13182, 28225),
    10**6: (1053, 9208, 9347, 63102, 138909, 139020, 294981),
    10**7: (
        11617, 96970, 97497, 675882, 1451554, 1451271, 3040588
    ),
}

# At X=10^7, counts on the common sample 2 <= n <= X of the left
# condition at n and the right condition at n+1, in EXPONENTS order.
EXPECTED_FINAL_LEFT = (437056, 2719287, 5571397)
EXPECTED_FINAL_RIGHT = (437056, 2719287, 5571397)


def largest_prime_factors(limit: int) -> array:
    """Return P(n) for 0 <= n <= limit, with the harmless convention P(1)=1."""
    lpf = array("I", [0]) * (limit + 1)
    if limit >= 1:
        lpf[1] = 1
    for prime in range(2, limit + 1):
        if lpf[prime] == 0:
            # Primes are visited in increasing order, so the final divisor
            # assigned to each integer is exactly its largest prime divisor.
            for multiple in range(prime, limit + 1, prime):
                lpf[multiple] = prime
    return lpf


def largest_prime_factor_trial(n: int) -> int:
    """Independent trial-division implementation, used only on a prefix."""
    if n == 1:
        return 1
    remaining = n
    answer = 1
    divisor = 2
    while divisor * divisor <= remaining:
        while remaining % divisor == 0:
            answer = divisor
            remaining //= divisor
        divisor = 3 if divisor == 2 else divisor + 2
    if remaining > 1:
        answer = remaining
    return answer


def independent_full_lpf_check(reference: array) -> None:
    """
    Recompute every P(n) by a structurally different method.

    First find the smallest prime factor by marking only previously-unmarked
    composites starting at p^2.  Then use

        P(n) = max(spf(n), P(n / spf(n)))

    in increasing n order.  This does not reuse the prime list or the
    overwriting-by-increasing-primes logic of largest_prime_factors().
    """
    limit = len(reference) - 1
    spf = array("I", [0]) * (limit + 1)
    for prime in range(2, limit + 1):
        if spf[prime] == 0:
            spf[prime] = prime
            if prime * prime <= limit:
                for multiple in range(prime * prime, limit + 1, prime):
                    if spf[multiple] == 0:
                        spf[multiple] = prime

    alternative = array("I", [0]) * (limit + 1)
    if limit >= 1:
        alternative[1] = 1
    for n in range(2, limit + 1):
        smallest = spf[n]
        cofactor_lpf = alternative[n // smallest]
        alternative[n] = (
            smallest if smallest > cofactor_lpf else cofactor_lpf
        )

    for n, (first, second) in enumerate(zip(reference, alternative)):
        assert first == second, (n, first, second)


def smooth_exact(prime_factor: int, n: int, exponent: Exponent) -> bool:
    """Integer-only form of prime_factor < n**(numerator/denominator)."""
    return (
        prime_factor ** exponent.denominator
        < n ** exponent.numerator
    )


def build_flags(lpf: array, limit: int) -> dict[Exponent, bytearray]:
    flags = {exponent: bytearray(limit + 1) for exponent in {
        exponent for pair in PAIRS for exponent in pair
    }}
    for n in range(2, limit + 1):
        prime_factor = lpf[n]
        flags[ONE_THIRD][n] = prime_factor**3 < n
        flags[ONE_HALF][n] = prime_factor**2 < n
        flags[TWO_THIRDS][n] = prime_factor**3 < n**2
    return flags


def exact_census(
    flags: dict[Exponent, bytearray],
    checkpoints: tuple[int, ...],
) -> tuple[
    dict[int, tuple[int, ...]],
    dict[int, tuple[int, ...]],
    dict[int, tuple[int, ...]],
]:
    running_joint = [0] * len(PAIRS)
    running_left = [0] * len(EXPONENTS)
    running_right = [0] * len(EXPONENTS)
    joint_result: dict[int, tuple[int, ...]] = {}
    left_result: dict[int, tuple[int, ...]] = {}
    right_result: dict[int, tuple[int, ...]] = {}
    checkpoint_set = set(checkpoints)
    last = checkpoints[-1]
    third = flags[ONE_THIRD]
    half = flags[ONE_HALF]
    two_thirds = flags[TWO_THIRDS]
    for n in range(2, last + 1):
        left_third, left_half, left_two_thirds = (
            third[n], half[n], two_thirds[n]
        )
        right_third, right_half, right_two_thirds = (
            third[n + 1], half[n + 1], two_thirds[n + 1]
        )
        running_left[0] += left_third
        running_left[1] += left_half
        running_left[2] += left_two_thirds
        running_right[0] += right_third
        running_right[1] += right_half
        running_right[2] += right_two_thirds
        running_joint[0] += left_third & right_third
        running_joint[1] += left_third & right_half
        running_joint[2] += left_half & right_third
        running_joint[3] += left_half & right_half
        running_joint[4] += left_half & right_two_thirds
        running_joint[5] += left_two_thirds & right_half
        running_joint[6] += left_two_thirds & right_two_thirds
        if n in checkpoint_set:
            joint_result[n] = tuple(running_joint)
            left_result[n] = tuple(running_left)
            right_result[n] = tuple(running_right)
    return joint_result, left_result, right_result


def rho_for_used_exponents(exponent: Exponent) -> float:
    """
    Dickman rho at u=1/gamma for gamma in {1/2,2/3,1/3}.

    For 1 <= u <= 2, rho(u)=1-log(u).  At u=3, numerical integration
    of rho'(u)=-rho(u-1)/u gives the displayed value.  Composite Simpson
    with 2^18 panels is deterministic and far more accurate than the table
    precision used in the report.
    """
    u = exponent.denominator / exponent.numerator
    if u <= 2:
        return 1.0 - log(u)
    assert u == 3
    panels = 1 << 18
    step = 1.0 / panels

    def integrand(t: float) -> float:
        return (1.0 - log(t - 1.0)) / t

    total = integrand(2.0) + integrand(3.0)
    total += 4.0 * sum(integrand(2.0 + j * step)
                       for j in range(1, panels, 2))
    total += 2.0 * sum(integrand(2.0 + j * step)
                       for j in range(2, panels, 2))
    integral = total * step / 3.0
    return 1.0 - log(2.0) - integral


def verify_prefix(lpf: array, flags: dict[Exponent, bytearray], prefix: int) -> None:
    for n in range(1, prefix + 1):
        trial = largest_prime_factor_trial(n)
        assert lpf[n] == trial, (n, lpf[n], trial)
        if n >= 2:
            for exponent in (ONE_THIRD, ONE_HALF, TWO_THIRDS):
                assert bool(flags[exponent][n]) == smooth_exact(
                    trial, n, exponent
                )


def verify_large_factor_reduction_prefix(
    lpf: array,
    flags: dict[Exponent, bytearray],
    prefix: int,
) -> None:
    """Exhaustively audit the exact Section 2.2 reduction for gamma=2/3."""
    smooth = flags[TWO_THIRDS]
    live_joint = 0
    bad_left_count = 0
    bad_right_count = 0
    bad_joint = 0
    quadruples = 0
    for n in range(2, prefix + 1):
        bad_left = not smooth[n]
        bad_right = not smooth[n + 1]
        live_joint += smooth[n] & smooth[n + 1]
        bad_left_count += bad_left
        bad_right_count += bad_right
        bad_joint += bad_left & bad_right
        if bad_left and bad_right:
            p, q = lpf[n], lpf[n + 1]
            a, b = n // p, (n + 1) // q
            assert a * p == n and b * q == n + 1
            assert b * q - a * p == 1
            assert p**3 >= n**2 and q**3 >= (n + 1) ** 2
            assert a < p and b < q
            quadruples += 1

    sample_size = prefix - 1
    assert live_joint == (
        sample_size - bad_left_count - bad_right_count + bad_joint
    )
    assert quadruples == bad_joint


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--limit", type=int, default=10**7)
    parser.add_argument("--trial-prefix", type=int, default=20_000)
    parser.add_argument(
        "--skip-full-cross-check",
        action="store_true",
        help="skip the independent all-n smallest-factor reconstruction",
    )
    args = parser.parse_args()
    if args.limit < 2:
        raise SystemExit("--limit must be at least 2")

    checkpoints = tuple(x for x in DEFAULT_CHECKPOINTS if x <= args.limit)
    if not checkpoints or checkpoints[-1] != args.limit:
        checkpoints += (args.limit,)

    started = perf_counter()
    lpf = largest_prime_factors(args.limit + 1)
    sieve_done = perf_counter()
    if not args.skip_full_cross_check:
        independent_full_lpf_check(lpf)
        gc.collect()
    full_check_done = perf_counter()
    flags = build_flags(lpf, args.limit + 1)
    flags_done = perf_counter()
    verify_prefix(lpf, flags, min(args.trial_prefix, args.limit + 1))
    verify_large_factor_reduction_prefix(
        lpf, flags, min(args.trial_prefix, args.limit)
    )
    prefix_done = perf_counter()
    counts, left_counts, right_counts = exact_census(flags, checkpoints)
    census_done = perf_counter()

    for cutoff, expected in EXPECTED_COUNTS.items():
        if cutoff in counts:
            assert counts[cutoff] == expected, (
                cutoff, counts[cutoff], expected
            )
    if args.limit == 10**7 and EXPECTED_FINAL_LEFT:
        assert left_counts[args.limit] == EXPECTED_FINAL_LEFT
        assert right_counts[args.limit] == EXPECTED_FINAL_RIGHT

    print("pair order:", " ".join(
        f"({alpha.label()},{beta.label()})" for alpha, beta in PAIRS
    ))
    print("exact counts (2 <= n <= X):")
    for cutoff in checkpoints:
        row = counts[cutoff]
        print(f"X={cutoff}: " + " ".join(map(str, row)))

    print("densities count/X at final X:")
    final_counts = counts[args.limit]
    for (alpha, beta), count in zip(PAIRS, final_counts):
        rho_product = (
            rho_for_used_exponents(alpha) * rho_for_used_exponents(beta)
        )
        print(
            f"  ({alpha.label()},{beta.label()}): "
            f"{count / args.limit:.12f}  "
            f"rho-product={rho_product:.12f}  "
            f"difference={count / args.limit - rho_product:+.12f}"
        )

    print("finite-sample independence ratios at final X:")
    sample_size = args.limit - 1
    exponent_index = {exponent: i for i, exponent in enumerate(EXPONENTS)}
    for (alpha, beta), joint in zip(PAIRS, final_counts):
        left = left_counts[args.limit][exponent_index[alpha]]
        right = right_counts[args.limit][exponent_index[beta]]
        ratio = joint * sample_size / (left * right)
        residual = joint * sample_size - left * right
        print(
            f"  ({alpha.label()},{beta.label()}): "
            f"left={left} right={right} joint={joint} "
            f"ratio={ratio:.12f} residual={residual}"
        )

    print(
        f"independent trial-division prefix 1.."
        f"{min(args.trial_prefix, args.limit + 1)}: OK"
    )
    print(
        "large-factor inclusion-exclusion/prime-quadruple prefix audit: OK"
    )
    print(
        "independent full smallest-factor reconstruction: "
        f"{'SKIPPED' if args.skip_full_cross_check else 'OK'}"
    )
    print(
        "timings seconds: "
        f"sieve={sieve_done-started:.3f} "
        f"full-check={full_check_done-sieve_done:.3f} "
        f"flags={flags_done-full_check_done:.3f} "
        f"prefix={prefix_done-flags_done:.3f} "
        f"census={census_done-prefix_done:.3f} "
        f"total={census_done-started:.3f}"
    )
    print("ALL CHECKS PASSED")


if __name__ == "__main__":
    main()

PARTIAL: Tao--Teräväinen gives the predicted ordinary density outside a zero-log-density set of scales; the supplied independently checked census is exact through 10^7, and the remaining fixed-shift all-scales task is precisely the prime-quadruple discrepancy (11).