#!/usr/bin/env python3
"""Verify the paper's finite arithmetic witnesses using the standard library.

All logical checks use integers or fractions.Fraction. Floating-point numbers
are used only to format the final human-readable interval-bound summary.

The pure-product lower bounds at 7, 9 and 10 factors are theorems of Maltby,
recorded explicitly by Cipu (2004, p. 39). They are literature inputs. This
program checks the explicit products attaining those bounds and every rational
inequality used to extend their separation from ideal PTE witnesses to a real
parameter interval. It does not claim to reproduce Maltby's exhaustive search.

Run from any directory:
    python /path/to/supplement/verify_pte.py
"""

from __future__ import annotations

import argparse
from fractions import Fraction
import json
from pathlib import Path


def require(condition: bool, message: str) -> None:
    if not condition:
        raise AssertionError(message)


def multiply(a: list[int], b: list[int]) -> list[int]:
    out = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            out[i + j] += x * y
    return out


def pure_product(lengths: list[int]) -> list[int]:
    out = [1]
    for length in lengths:
        require(isinstance(length, int) and length > 0,
                "Every product length must be a positive integer")
        factor = [1] + [0] * (length - 1) + [-1]
        out = multiply(out, factor)
    return out


def divide_one_minus_z(q: list[int]) -> list[int]:
    require(len(q) >= 2, "A nonzero constant is not divisible by 1-z")
    acc = 0
    partial_sums = []
    for coefficient in q:
        acc += coefficient
        partial_sums.append(acc)
    require(partial_sums[-1] == 0, "Nonzero remainder on division by 1-z")
    quotient = partial_sums[:-1]
    require(multiply(quotient, [1, -1]) == q, "Reconstruction failed")
    return quotient


def incidence(record: dict) -> list[int]:
    a, b = record["A"], record["B"]
    require(all(isinstance(x, int) and x >= 0 for x in a + b),
            f"{record['id']}: exponents must be nonnegative integers")
    require(len(a) == len(set(a)) and len(b) == len(set(b)),
            f"{record['id']}: a repeated element was found")
    require(not set(a).intersection(b), f"{record['id']}: supports intersect")
    q = [0] * (max(a + b) + 1)
    for x in a:
        q[x] = 1
    for x in b:
        q[x] = -1
    return q


def verify_pair(record: dict) -> dict:
    label, m = record["id"], record["m"]
    a, b = record["A"], record["B"]
    require(len(a) == len(b) == m, f"{label}: wrong cardinality")
    for r in range(m):
        require(sum(x ** r for x in a) == sum(x ** r for x in b),
                f"{label}: moment {r} does not cancel")
    require(sum(x ** m for x in a) != sum(x ** m for x in b),
            f"{label}: unexpected cancellation of moment {m}")
    q = incidence(record)
    require(q[-1] == (-1) ** m, f"{label}: wrong leading sign")
    norm = sum(x * x for x in q)
    require(norm == 2 * m, f"{label}: squared coefficient norm is not 2m")
    p = q
    for _ in range(m):
        p = divide_one_minus_z(p)
    require(p[-1] == 1, f"{label}: quotient is not monic")
    require(sum(p) != 0, f"{label}: zero order exceeds m")
    require(len(p) - 1 == record["expected_degree"], f"{label}: wrong degree")
    if "expected_coefficients" in record:
        require(p == record["expected_coefficients"], f"{label}: coefficient mismatch")
    if record.get("require_nonpalindromic"):
        require(p != p[::-1], f"{label}: expected a nonpalindromic quotient")
    if "pure_product_lengths" in record:
        require(pure_product(record["pure_product_lengths"]) == q,
                f"{label}: pure-product expansion mismatch")
    return {"id": label, "m": m, "degree": len(p) - 1,
            "norm_squared": norm, "palindromic": p == p[::-1]}


def fourier_majorant(delta: Fraction, max_index: int) -> list[Fraction]:
    """rho[n] = c_n(-delta)/c_0(-delta), n >= 1, and rho[0] = 1.

    For every |beta| <= delta < 1/2, |c_n(beta)|/c_0(beta) <= rho[n].
    """
    require(0 < delta < Fraction(1, 2), "delta must lie in (0,1/2)")
    rho = [Fraction(1)]
    if max_index == 0:
        return rho
    rho.append(delta / (1 - delta))
    for n in range(1, max_index):
        rho.append(rho[-1] * (n + delta) / (n + 1 - delta))
    require(all(rho[n + 1] < rho[n] for n in range(1, max_index)),
            "Fourier majorant must decrease")
    return rho


def verify_exclusion(case: dict, pairs: dict[str, dict]) -> dict:
    m = case["m"]
    record = pairs[case["pair_id"]]
    require(record["m"] == m, "Exclusion witness has the wrong size")
    a, b = record["A"], record["B"]
    g = case["pure_norm_lower"]
    delta = Fraction(*case["delta"])
    eta = Fraction(*case["gap_lower"])
    require(g % 2 == 0 and g > 2 * m, "A positive even pure-product gap is required")
    cap = 2 ** (m - 1)
    rho = fourier_majorant(delta, max(cap, max(a + b) - min(a + b)))
    # H[p] is the sum of the p largest possible off-diagonal entries
    # in a row indexed by distinct integer exponents: each distance occurs
    # at most twice, so the maximizing distances are 1,1,2,2,... .
    h = [Fraction(0)]
    for p in range(1, cap + 1):
        h.append(h[-1] + rho[(p + 1) // 2])
    require(h[cap] < 1, "The uniform row bound must be below one")
    lower_candidates = [(2 * s * (1 - h[s]), s)
                        for s in range(g // 2, cap + 1)]
    lower, minimizing_s = min(lower_candidates)
    # For S > 2*cap the bound is S*(1-H[cap]), hence increases in S.
    cross = 2 * sum((rho[abs(x - y)] for x in a for y in b), Fraction(0))
    same = 2 * sum((rho[abs(v[i] - v[j])]
                    for v in (a, b) for i in range(len(v)) for j in range(i)),
                   Fraction(0))
    upper = 2 * m + max(cross, same)
    gap = lower - upper
    require(gap > eta, f"m={m}: the claimed uniform separation failed")
    lower_rational = Fraction(*case["lower_rational"])
    upper_rational = Fraction(*case["upper_rational"])
    require(lower > lower_rational and upper < upper_rational,
            f"m={m}: a displayed rational bound failed")
    require(lower_rational - upper_rational > eta,
            f"m={m}: the displayed bounds do not imply the claimed gap")
    require(all(bound - upper > eta for bound, _ in lower_candidates),
            f"m={m}: an intermediate even norm failed")
    return {"m": m, "delta": str(delta), "proved_gap": str(eta),
            "gap_decimal": float(gap), "lower_decimal": float(lower),
            "upper_decimal": float(upper), "minimizing_norm": 2 * minimizing_s,
            "rational_norm_checks": len(lower_candidates)}


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--pairs", type=Path,
                        default=Path(__file__).with_name("pte_pairs.json"))
    args = parser.parse_args()
    data = json.loads(args.pairs.read_text(encoding="utf-8"))
    require(data["schema_version"] == 1, "Unknown data schema")
    records = data["pairs"]
    ids = [record["id"] for record in records]
    require(len(ids) == len(set(ids)), "Duplicate record id")
    pairs = {record["id"]: record for record in records}
    for record in records:
        report = verify_pair(record)
        print(f"PTE {report['id']}: moments=OK, degree={report['degree']}, "
              f"norm_squared={report['norm_squared']}, monic=YES, "
              f"palindromic={'YES' if report['palindromic'] else 'NO'}")
    q12 = [-x for x in incidence(pairs["m12"])]
    p11 = q12
    for _ in range(11):
        p11 = divide_one_minus_z(p11)
    require(p11[-1] == 1 and sum(x * x for x in q12) == 24,
            "The size-twelve witness did not give the stated m=11 upper bound")
    print("m=11 upper-bound witness: monic quotient and norm_squared=24 verified")
    for record in data["pure_product_witnesses"]:
        require(len(record["lengths"]) == record["m"], "Wrong factor count")
        q = pure_product(record["lengths"])
        require(sum(abs(x) for x in q) == record["expected_l1"], "Pure-product l1 mismatch")
        require(sum(x * x for x in q) == record["expected_l2_squared"],
                "Pure-product squared l2 mismatch")
        print(f"Pure product m={record['m']}: l1={record['expected_l1']}, "
              f"norm_squared={record['expected_l2_squared']}, exact expansion=OK")
    print("Uniform exclusion bounds, using the cited pure-product lower bounds:")
    for case in data["product_exclusion_intervals"]:
        report = verify_exclusion(case, pairs)
        print(f"  m={report['m']}, delta={report['delta']}, "
              f"gap/c0>{report['proved_gap']}, "
              f"computed rational bound={report['gap_decimal']:.12f}, "
              f"norm cases={report['rational_norm_checks']}")
    print("All integer and rational checks passed.")


if __name__ == "__main__":
    main()