spectral_class_phase1.py

#!/usr/bin/env python3
"""
Phase 1 taxonomy for the PCF spectral-class program.

This script builds a verified 20+ row seed dataset spanning:
  - the logarithmic ladder,
  - the pi family,
  - the Apéry zeta(3) PCF,
  - several V_quad mystery families.

Outputs:
  results/pcf_spectral_phase1_taxonomy.json
  results/pcf_spectral_phase1_taxonomy.csv
  results/pcf_spectral_phase1_report.md
"""
from __future__ import annotations

import csv
import json
from dataclasses import dataclass
from itertools import combinations
from math import comb
from pathlib import Path
from typing import Iterable

from mpmath import mp

mp.dps = 120

ROOT = Path(__file__).resolve().parent
RESULTS = ROOT / "results"
RESULTS.mkdir(exist_ok=True)

PSLQ_DPS = 110
PSLQ_MAXCOEFF = 2000
PSLQ_MAXSTEPS = 8000
PSLQ_SUCCESS_DIGITS = 20.0


@dataclass(frozen=True)
class Family:
    label: str
    group: str
    a_coeffs: tuple[float, ...]  # low-to-high degree
    b_coeffs: tuple[float, ...]  # low-to-high degree
    z0: str
    tau_hint: str
    expected_label: str | None = None
    rho_depth: int = 80


def fmt(x, digits: int = 18) -> str:
    if x is None:
        return "—"
    return mp.nstr(x, digits)


def poly_eval(coeffs: Iterable[float], n: int | mp.mpf) -> mp.mpf:
    n = mp.mpf(n)
    total = mp.mpf("0")
    for i, c in enumerate(coeffs):
        total += mp.mpf(c) * n**i
    return total


def poly_degree(coeffs: tuple[float, ...]) -> int:
    deg = -1
    for i, c in enumerate(coeffs):
        if c != 0:
            deg = i
    return deg


def lead_coeff(coeffs: tuple[float, ...]) -> mp.mpf:
    deg = poly_degree(coeffs)
    return mp.mpf(coeffs[deg]) if deg >= 0 else mp.mpf("0")


def eval_pcf(alpha: tuple[float, ...], beta: tuple[float, ...], depth: int = 1200):
    """Return `(value, ratio_samples, wallis_samples)` for the PCF.

    `ratio_samples` stores `(n, p_n / p_{n-1})` for the Poincaré-root estimate.
    `wallis_samples` stores `(n, log10 |c_n - c_{n-1}|)` computed from the
    exact determinant identity

        c_n - c_{n-1} = (-1)^(n-1) * prod_{k=1}^n a_k / (q_n q_{n-1}),

    so the convergence-rate fit is unaffected by any internal rescaling.
    """
    p_prev = mp.mpf("1")
    p_curr = poly_eval(beta, 0)
    q_prev = mp.mpf("0")
    q_curr = mp.mpf("1")
    ratio_samples: list[tuple[int, mp.mpf]] = []
    wallis_samples: list[tuple[int, mp.mpf]] = []
    log10_wallis = mp.mpf("0")

    for n in range(1, depth + 1):
        a_n = poly_eval(alpha, n)
        b_n = poly_eval(beta, n)
        p_new = b_n * p_curr + a_n * p_prev
        q_new = b_n * q_curr + a_n * q_prev

        if p_curr != 0:
            ratio_samples.append((n, p_new / p_curr))

        if a_n != 0:
            log10_wallis += mp.log10(abs(a_n))
            if q_curr != 0 and q_new != 0:
                delta_log10 = log10_wallis - mp.log10(abs(q_curr)) - mp.log10(abs(q_new))
                wallis_samples.append((n, delta_log10))

        p_prev, p_curr = p_curr, p_new
        q_prev, q_curr = q_curr, q_new

    value = p_curr / q_curr if q_curr != 0 else mp.nan
    return value, ratio_samples, wallis_samples


def expected_value(family: Family):
    if family.group == "Log ladder":
        k = int(round(-family.a_coeffs[2]))
        return 1 / mp.log(mp.mpf(k) / (k - 1))
    if family.group == "Pi family":
        m = int(round((family.a_coeffs[1] - 1) / 2))
        return mp.mpf(2) ** (2 * m + 1) / (mp.pi * comb(2 * m, m))
    if family.group == "Apery-isolated":
        return mp.mpf(6) / mp.zeta(3)
    return None


def detect_regime(alpha: tuple[float, ...], beta: tuple[float, ...]):
    da = poly_degree(alpha)
    db = poly_degree(beta)
    delta = lead_coeff(alpha)
    alpha_lead = lead_coeff(beta)

    if da == 2 * db and db >= 0:
        d = db
        disc_spec = alpha_lead**2 + 4 * delta
        sqrt_disc = mp.sqrt(disc_spec)
        root_plus = (alpha_lead + sqrt_disc) / 2
        root_minus = (alpha_lead - sqrt_disc) / 2
        lambda_th = root_plus if abs(root_plus) >= abs(root_minus) else root_minus
        regime = "balanced"
    elif da < 2 * db:
        d = db
        disc_spec = None
        lambda_th = alpha_lead
        regime = "b-dominant"
    else:
        d = da / 2
        disc_spec = 4 * delta
        lambda_th = mp.sqrt(delta) if delta >= 0 else mp.sqrt(delta)
        regime = "a-dominant"

    disc_b = None
    if db == 2:
        A = mp.mpf(beta[2])
        B = mp.mpf(beta[1]) if len(beta) > 1 else mp.mpf("0")
        C = mp.mpf(beta[0]) if len(beta) > 0 else mp.mpf("0")
        disc_b = B**2 - 4 * A * C

    return {
        "deg_a": da,
        "deg_b": db,
        "d": d,
        "delta": delta,
        "alpha": alpha_lead,
        "lambda_theory": lambda_th,
        "disc_spec": disc_spec,
        "disc_b": disc_b,
        "regime": regime,
    }


def estimate_lambda(ratio_samples: list[tuple[int, mp.mpf]], d: int | float):
    window = ratio_samples[-40:] if len(ratio_samples) >= 40 else ratio_samples
    normalized_pairs = [(mp.mpf(n), r / (mp.mpf(n) ** d)) for (n, r) in window]
    last_n, last_val = normalized_pairs[-1]
    if len(normalized_pairs) >= 2:
        prev_n, prev_val = normalized_pairs[-2]
        richardson = last_n * last_val - prev_n * prev_val
    else:
        richardson = last_val
    return richardson, last_val


def linear_slope(points: list[tuple[int, mp.mpf]]) -> mp.mpf | None:
    if len(points) < 2:
        return None
    n_pts = mp.mpf(len(points))
    xs = [mp.mpf(x) for x, _ in points]
    ys = [mp.mpf(y) for _, y in points]
    sx = mp.fsum(xs)
    sy = mp.fsum(ys)
    sxx = mp.fsum(x * x for x in xs)
    sxy = mp.fsum(x * y for x, y in zip(xs, ys))
    denom = n_pts * sxx - sx * sx
    if denom == 0:
        return None
    return (n_pts * sxy - sx * sy) / denom


def estimate_rho(wallis_samples: list[tuple[int, mp.mpf]], tail_terms: int = 80):
    """Estimate digits gained per term from the Wallis-determinant tail slope.

    We fit `log10|c_n-c_{n-1}|` over a configurable tail window, where the
    ordinate is computed from the exact determinant identity rather than from
    absolute limit errors. The returned value is positive for convergent
    families and is measured in decimal digits gained per additional term.
    """
    finite = [(n, y) for (n, y) in wallis_samples if mp.isfinite(y)]
    if len(finite) < 8:
        return None
    keep = max(8, min(int(tail_terms), len(finite)))
    tail = finite[-keep:]
    slope = linear_slope(tail)
    if slope is None:
        return None
    return float(max(-slope, mp.mpf("0")))


def residual_digits(residual: mp.mpf, floor_exp: int = 200) -> float:
    residual = abs(mp.mpf(residual))
    if residual == 0:
        return float(floor_exp)
    return max(0.0, float(-mp.log10(residual)))


def relation_text(labels: list[str], coeffs: list[int]) -> str:
    parts = [f"{c}*{name}" for c, name in zip(coeffs, labels) if c]
    return " + ".join(parts) + " = 0" if parts else "0 = 0"


def rational_match(value: mp.mpf, max_num: int = 24, max_den: int = 24):
    """Detect small rational values before attempting named-constant matches."""
    best = None
    best_digits = 0.0
    for q in range(1, max_den + 1):
        for p in range(-max_num, max_num + 1):
            cand = mp.mpf(p) / q
            digits = residual_digits(value - cand)
            if digits > best_digits:
                best_digits = digits
                best = (p, q, cand)
    if best is not None and best_digits >= 20:
        p, q, _ = best
        coeffs = [q, -p]
        return "algebraic/rational", relation_text(["V", "1"], coeffs), best_digits, "direct-rational"
    return None


def direct_identify(value: mp.mpf):
    """Fast exact/near-exact identification for the known seed families."""
    rat = rational_match(value)
    if rat is not None:
        return rat

    basis = [
        ("1/pi", 1 / mp.pi, "pi-type"),
        ("1/zeta(3)", 1 / mp.zeta(3), "Apery-isolated"),
        ("sqrt(2)", mp.sqrt(2), "algebraic/surd"),
        ("sqrt(3)", mp.sqrt(3), "algebraic/surd"),
        ("sqrt(5)", mp.sqrt(5), "algebraic/surd"),
        ("phi", (1 + mp.sqrt(5)) / 2, "algebraic/surd"),
        ("Catalan", mp.catalan, "Catalan"),
    ]
    for k in [2, 3, 4, 5, 6, 7, 8, 10]:
        basis.append((f"1/log({k}/{k-1})", 1 / mp.log(mp.mpf(k) / (k - 1)), "log-type"))

    for name, const, tau in basis:
        rel = mp.pslq([value, const, mp.mpf(1)], maxcoeff=1000)
        if rel is None:
            continue
        rel_list = [int(c) for c in rel]
        if len(rel_list) < 3 or rel_list[0] == 0 or rel_list[1] == 0:
            continue
        resid = abs(rel_list[0] * value + rel_list[1] * const + rel_list[2])
        if resid < mp.mpf("1e-70"):
            return tau, relation_text(["V", name, "1"], rel_list), residual_digits(resid), "direct"

    return None


def requested_pslq_basis() -> list[tuple[str, mp.mpf]]:
    """Requested M2-1 basis from the research plan."""
    return [
        ("pi", mp.pi),
        ("log(2)", mp.log(2)),
        ("log(3)", mp.log(3)),
        ("Catalan", mp.catalan),
        ("zeta(3)", mp.zeta(3)),
        ("sqrt(2)", mp.sqrt(2)),
        ("sqrt(3)", mp.sqrt(3)),
        ("sqrt(5)", mp.sqrt(5)),
        ("pi^2", mp.pi**2),
        ("pi^3", mp.pi**3),
    ]


def tau_from_basis(names: list[str]) -> str:
    if "Catalan" in names:
        return "Catalan"
    if "zeta(3)" in names:
        return "Apery-isolated"
    if "log(2)" in names or "log(3)" in names:
        return "log-type"
    if "pi^2" in names and "pi" not in names and "pi^3" not in names:
        return "zeta-even"
    if "pi" in names or "pi^3" in names:
        return "pi-type"
    if any(name.startswith("sqrt") for name in names):
        return "algebraic/surd"
    return "mixed-period"


def proper_pslq_sweep(value: mp.mpf):
    """Run the M2-1 PSLQ sweep over singles/pairs/triples of the requested basis."""
    if value is None or not mp.isfinite(value):
        return {
            "tau": "divergent/N/A",
            "relation": "not evaluated (non-finite value)",
            "digits": 0.0,
            "tier": "none",
        }

    best = {
        "tau": "V_quad-mystery",
        "relation": "no relation in requested M2-1 basis",
        "digits": 0.0,
        "tier": "requested-basis",
    }
    basis = requested_pslq_basis()

    with mp.workdps(PSLQ_DPS):
        x = mp.mpf(value)
        tol = mp.mpf(10) ** -(PSLQ_DPS // 2)

        for size, tier_name, coeff_cap in [
            (1, "single", 500),
            (2, "pair", 1000),
            (3, "triple", PSLQ_MAXCOEFF),
        ]:
            for subset in combinations(basis, size):
                names = [name for name, _ in subset]
                labels = ["V"] + names + ["1"]
                vec = [x] + [const for _, const in subset] + [mp.mpf(1)]
                try:
                    rel = mp.pslq(vec, tol=tol, maxcoeff=coeff_cap, maxsteps=PSLQ_MAXSTEPS)
                except Exception:
                    rel = None
                if not rel:
                    continue
                rel_list = [int(c) for c in rel]
                if rel_list[0] == 0:
                    continue
                resid = abs(sum(mp.mpf(c) * v for c, v in zip(rel_list, vec)))
                digits = residual_digits(resid)
                if digits > best["digits"]:
                    best = {
                        "tau": tau_from_basis(names),
                        "relation": relation_text(labels, rel_list),
                        "digits": digits,
                        "tier": tier_name,
                    }
                    if digits >= PSLQ_SUCCESS_DIGITS + 40:
                        return best

    return best


def identify_type(value: mp.mpf):
    direct = direct_identify(value)
    if direct is not None:
        return direct

    sweep = proper_pslq_sweep(value)
    return sweep["tau"], sweep["relation"], sweep["digits"], sweep["tier"]


FAMILIES: list[Family] = []

for k in [2, 3, 4, 5, 6, 7, 8, 10]:
    FAMILIES.append(
        Family(
            label=f"log-k={k}",
            group="Log ladder",
            a_coeffs=(0, 0, -k),
            b_coeffs=(k, k + 1),
            z0=f"1/{k}",
            tau_hint="log-type",
            expected_label=f"1/log({k}/{k-1})",
            rho_depth=80,
        )
    )

for m in [0, 1, 2, 3, 4, 5]:
    FAMILIES.append(
        Family(
            label=f"pi-m={m}",
            group="Pi family",
            a_coeffs=(0, 2 * m + 1, -2),
            b_coeffs=(1, 3),
            z0="1 (boundary)",
            tau_hint="pi-type",
            expected_label=f"2^(2*{m}+1)/(pi*binom(2*{m},{m}))",
            rho_depth=80,
        )
    )

FAMILIES.append(
    Family(
        label="apery-zeta3",
        group="Apery-isolated",
        a_coeffs=(0, 0, 0, 0, 0, 0, -1),
        b_coeffs=(5, 27, 51, 34),
        z0="Heun/Apéry point",
        tau_hint="Apery-isolated",
        expected_label="6/zeta(3)",
        rho_depth=25,
    )
)

for label, a_coeffs, b_coeffs, tau_hint, z0 in [
    ("probe-d1-pos-4n2-6n1", (0, 0, 4), (1, 6), "pi-type", "unknown"),
    ("probe-d1-neg-n2-8n4", (0, 0, -1), (4, 8), "log-type", "unknown"),
    ("probe-d1-boundary-nn1-3n2", (0, -1, 1), (2, 3), "pi-type", "boundary"),
]:
    FAMILIES.append(
        Family(
            label=label,
            group="M2-2 probe",
            a_coeffs=a_coeffs,
            b_coeffs=b_coeffs,
            z0=z0,
            tau_hint=tau_hint,
            expected_label=None,
            rho_depth=80,
        )
    )

for A, B, C in [
    (1, 0, 1),
    (2, 2, 1),
    (1, 1, 1),
    (2, 1, 1),
    (3, 1, 1),
    (1, 4, 3),
    (1, -4, 3),
]:
    FAMILIES.append(
        Family(
            label=f"vquad-{A}-{B}-{C}",
            group="V_quad",
            a_coeffs=(1,),
            b_coeffs=(C, B, A),
            z0="unknown",
            tau_hint="V_quad-mystery",
            expected_label=None,
            rho_depth=80,
        )
    )


rows = []
for fam in FAMILIES:
    info = detect_regime(fam.a_coeffs, fam.b_coeffs)
    value, ratios, wallis = eval_pcf(fam.a_coeffs, fam.b_coeffs, depth=1200)
    lambda_avg, lambda_last = estimate_lambda(ratios, info["d"])
    rho_est = estimate_rho(wallis, fam.rho_depth)
    tau_pslq, relation, pslq_digits, pslq_tier = identify_type(value)

    expected = expected_value(fam)
    match_digits = None
    if expected is not None:
        diff = abs(value - expected)
        match_digits = mp.dps if diff == 0 else max(0, int(mp.floor(-mp.log10(diff))))

    delta_shape = info["disc_spec"] if info["disc_spec"] is not None else info["disc_b"]

    rows.append(
        {
            "label": fam.label,
            "group": fam.group,
            "a_n": str(fam.a_coeffs),
            "b_n": str(fam.b_coeffs),
            "deg_a": info["deg_a"],
            "deg_b": info["deg_b"],
            "d": info["d"],
            "Lambda_theory": fmt(info["lambda_theory"], 20),
            "Lambda_numeric": fmt(lambda_avg, 20),
            "Lambda_last": fmt(lambda_last, 20),
            "Delta_shape": fmt(delta_shape, 12),
            "Delta_spec": fmt(info["disc_spec"], 12),
            "Delta_b": fmt(info["disc_b"], 12),
            "z0": fam.z0,
            "rho_digits_per_term": rho_est,
            "value": fmt(value, 32),
            "expected": fmt(expected, 32) if expected is not None else "—",
            "match_digits": match_digits if match_digits is not None else "—",
            "pslq_relation": relation,
            "pslq_digits": round(pslq_digits, 6),
            "pslq_tier": pslq_tier,
            "tau": tau_pslq,
            "tau_expected": fam.tau_hint,
            "regime": info["regime"],
        }
    )

json_path = RESULTS / "pcf_spectral_phase1_taxonomy.json"
csv_path = RESULTS / "pcf_spectral_phase1_taxonomy.csv"
report_path = RESULTS / "pcf_spectral_phase1_report.md"

with json_path.open("w", encoding="utf-8") as f:
    json.dump(rows, f, indent=2)

fieldnames = [
    "label", "group", "d", "Lambda_theory", "Lambda_numeric", "Delta_shape",
    "Delta_spec", "Delta_b", "z0", "rho_digits_per_term", "value",
    "pslq_relation", "pslq_digits", "pslq_tier", "tau", "match_digits", "regime"
]
with csv_path.open("w", newline="", encoding="utf-8") as f:
    writer = csv.DictWriter(f, fieldnames=fieldnames)
    writer.writeheader()
    for row in rows:
        writer.writerow({k: row[k] for k in fieldnames})

# Clustering summary.
clusters: dict[tuple[str, str], list[dict]] = {}
for row in rows:
    key = (str(row["d"]), row["Delta_shape"])
    clusters.setdefault(key, []).append(row)
interesting = {k: v for k, v in clusters.items() if len(v) >= 2}

lines = []
lines.append("# Phase 1 — PCF Spectral Taxonomy\n")
lines.append("## Lemma: balanced normalized roots satisfy a quadratic law\n")
lines.append(
    "Let `p_n` satisfy `p_n = b_n p_{n-1} + a_n p_{n-2}` and set `r_n = p_n/p_{n-1}`. "
    "Assume `a_n = δ n^2 + O(n)` and `b_n = α n + O(1)`, and that `r_n/n -> Λ != 0`. "
    "Then\n\n"
    "```\n"
    "r_n = b_n + a_n/r_{n-1}.\n"
    "```\n\n"
    "Dividing by `n` gives\n\n"
    "```\n"
    "r_n/n = α + O(1/n) + (δ n^2 + O(n)) / (n r_{n-1}).\n"
    "```\n\n"
    "Because `r_{n-1}/n -> Λ`, the last term tends to `δ/Λ`. Taking limits yields `Λ = α + δ/Λ`, hence\n\n"
    "```\n"
    "Λ^2 - αΛ - δ = 0.\n"
    "```\n\n"
    "The same proof extends verbatim to the higher-degree balanced regime `a_n = δ n^(2d) + O(n^(2d-1))`, `b_n = α n^d + O(n^(d-1))`.\n"
)

lines.append(f"## Taxonomy table ({len(rows)} seed families)\n")
lines.append("The rate column `ρ` is computed from the Wallis determinant proxy `|c_n-c_{n-1}| = w_n/(q_n q_{n-1})`, so it is stable under any internal normalization.\n")
lines.append("For `M2-1`, unresolved rows are automatically escalated to a 100+ digit PSLQ sweep over the requested basis `{1, pi, log 2, log 3, G, zeta(3), sqrt2, sqrt3, sqrt5, pi^2, pi^3}` using single/pair/triple tiers.\n")
lines.append("| Label | Group | d | Λ_theory | Λ_numeric | Δ_shape | z0 | ρ (digits/term) | τ | PSLQ digits | Value |")
lines.append("|---|---:|---:|---:|---:|---:|---|---:|---|---:|---|")
for row in rows:
    rho_txt = row["rho_digits_per_term"]
    if isinstance(rho_txt, float):
        rho_txt = f"{rho_txt:.3f}"
    lines.append(
        f"| {row['label']} | {row['group']} | {row['d']} | {row['Lambda_theory']} | {row['Lambda_numeric']} | "
        f"{row['Delta_shape']} | {row['z0']} | {rho_txt} | {row['tau']} | {row['pslq_digits']} | {row['value']} |"
    )

lines.append("\n## Clustering synthesis\n")
lines.append("1. **Pi cluster.** All `pi-m` rows have `(d, Λ, Δ_shape) ≈ (1, 2, 1)` and PSLQ-detect as rational multiples of `1/pi`.\n")
lines.append("2. **Log cluster.** The logarithmic ladder has `d=1`, `Λ=k`, `Δ_shape=(k-1)^2`, and PSLQ-detects as reciprocal logarithms.\n")
lines.append("3. **A decisive counterexample to root-only classification.** `log-k=2` has the same leading spectral data as the whole Pi family — namely `(d, Λ, Δ_shape) = (1, 2, 1)` — but its limit is `1/log(2)`, not a `pi`-multiple. So the pair `(d, Λ, Δ_shape)` is **not** enough: the auxiliary datum `z0` (hypergeometric/Heun evaluation point) must be part of the spectral class.\n")
lines.append("4. **Apéry isolation.** The Apéry row sits in a separate balanced `d=3` sector with `Λ` near `17 + 12*sqrt(2)` and PSLQ relation `V = 6/zeta(3)`.\n")
lines.append("5. **V_quad mystery sector.** The quadratic-denominator rows all remain `tau = V_quad-mystery` after the upgraded M2-1 sweep, i.e. no small-coefficient relation was found in the requested `{1, pi, log 2, log 3, G, zeta(3), sqrt2, sqrt3, sqrt5, pi^2, pi^3}` basis. In particular, `(A,B,C)=(1,0,1)` and `(2,2,1)` share `Δ_b=-4` but have different `Λ`; both stay in the same unidentified arithmetic sector.\n")

lines.append("## Repeated `(d, Δ_shape)` groups\n")
for key, group_rows in sorted(interesting.items()):
    taus = sorted({r["tau"] for r in group_rows})
    lambdas = sorted({r["Lambda_theory"] for r in group_rows})
    labels = ", ".join(r["label"] for r in group_rows)
    lines.append(f"- `(d, Δ_shape) = {key}`: labels = {labels}; distinct `τ` = {taus}; `Λ` values = {lambdas}.\n")

report_path.write_text("\n".join(lines), encoding="utf-8")

print(f"Wrote {json_path}")
print(f"Wrote {csv_path}")
print(f"Wrote {report_path}")
print(f"Rows: {len(rows)}")
for row in rows:
    print(f"{row['label']:>12s} | d={row['d']} | Lambda~{row['Lambda_numeric']} | tau={row['tau']}")