"""Stage 2–3 automation: exhaustive search for exact integer / modular identities among the secp256k1 constants.

The generator enumerates candidate relations of bounded complexity and the checker evaluates each one with exact
arithmetic. Every surviving identity is then classified against the literature (CM theory of j = 0 curves,
GLV lattice structure). This is a deliberately simple "discovery engine": its purpose is to measure how many of the
known structural identities an automated search rediscovers, and whether it finds any that are NOT explained by theory.

Search spaces
  Q1  binary quadratic forms: u² + c·u·v + d·v² = target  for (u, v) among constant pairs, |c|,|d| ≤ 6, target ∈ {p, n, 4p, 4n, …}
  Q2  modular polynomial relations: x² + c·x + d ≡ 0 (mod m) for x ∈ {β, λ, …}, m ∈ {p, n}, |c|,|d| ≤ 6
  Q3  linear relations among the GLV basis entries: a2 = e1·a1 + e2·b1, b2 = e3·a1 + e4·b1, |e_i| ≤ 3
  Q4  "λ ≡ polynomial(β)" style cross-modulus relations (expected to fail: different moduli) — negative control
  Q5  small-coefficient relations between Gx, Gy, β, p (expected to fail — the generator is not special) — negative control
Writes results/automated_discovery/identity_search.json.
"""
import itertools, json, math, sys, time, pathlib
ROOT = pathlib.Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT / "src"))
from wp34 import secp256k1 as k1
from wp34.provenance import record

p, n, beta, lam = k1.p, k1.n, k1.beta, k1.lam
t = p + 1 - n
f = math.isqrt((4 * p - t * t) // 3)
consts = {"p": p, "n": n, "t": t, "f": f, "a1": k1.a1, "b1": k1.b1, "a2": k1.a2, "b2": k1.b2,
          "beta": beta, "lambda": lam, "Gx": k1.Gx, "Gy": k1.Gy}
KNOWN = {
    "Q1:t^2+0*t*f+3*f^2=4p": "CM theory: 4p = t² + 3f² for j = 0 (Frobenius π = (t + f√−3)/2 ∈ Z[ω])",
    "Q1:a1^2-1*a1*b1+1*b1^2=n": "n = N(a1 + b1·ω): the GLV basis vector has Eisenstein norm n (GLV 2001 / Z[ω] structure)",
    "Q1:a2^2-1*a2*b2+1*b2^2=n": "same for the second basis vector",
    "Q3:a2=1*a1-1*b1": "second GLV vector = ω·(first) up to sign: (a1,b1)·ω rotates the lattice (Eisenstein units)",
    "Q3:b2=1*a1+0*b1": "same",
    "Q2:beta^2+1*beta+1=0 mod p": "β is a primitive cube root of unity (β² + β + 1 = 0)",
    "Q2:lambda^2+1*lambda+1=0 mod n": "λ is a primitive cube root of unity in Z/nZ",
    "Q1:t^2+0*t*b1+3*b1^2=4p": "b1 = −f: the GLV vector's second coordinate is the CM conductor (4p = t² + 3f²)",
}
# Closed form found by the search and then explained (CM theory): the GLV lattice {(k1,k2): k1 + k2·λ ≡ 0 (mod n)} is the
# principal ideal (π̄ − 1)·Z[ω] of Z[ω] where π̄ = (t − f√−3)/2 is the Frobenius conjugate mapping ω ↦ λ, hence
#   a1 = (t − f)/2 − 1,  b1 = −f,  a2 = (t + f)/2 − 1 = a1 − b1,  b2 = a1,   and  a1² − a1·b1 + b1² = N(π̄ − 1) = n.
CLOSED_FORM = {"a1": (t - f) // 2 - 1, "b1": -f, "a2": (t + f) // 2 - 1, "b2": (t - f) // 2 - 1}
EISENSTEIN_NORM_FAMILY = "quadratic form u² ± uv + v² = n among {a1, b1, a2, b2, f}: all are Eisenstein norms of unit multiples / conjugates of π̄ − 1 (CM theory, GLV lattice = ideal (π̄ − 1))"


def search():
    found, tested = [], 0
    names = list(consts)
    targets = {"p": p, "n": n, "4p": 4 * p, "4n": 4 * n, "2p": 2 * p, "3p": 3 * p, "3n": 3 * n, "p-1": p - 1, "n-1": n - 1}
    # Q1 binary quadratic forms (ordered pairs; forms u² + c uv + d v²)
    for u, v in itertools.permutations(names, 2):
        U, V = consts[u], consts[v]
        for c in range(-6, 7):
            for d in range(-6, 7):
                val = U * U + c * U * V + d * V * V
                tested += 1
                for tn, tv in targets.items():
                    if val == tv:
                        cs = f"{c:+d}" if c else ""
                        found.append({"space": "Q1", "id": f"Q1:{u}^2{c:+d}*{u}*{v}{d:+d}*{v}^2={tn}".replace("+", "+").replace("^2+0*", "^2+0*"),
                                      "relation": f"{u}^2 {c:+d}·{u}·{v} {d:+d}·{v}^2 = {tn}"})
    # Q2 modular quadratics
    for xname in ("beta", "lambda", "Gx", "Gy", "t", "f", "a1", "b1"):
        X = consts[xname]
        for mname, m in (("p", p), ("n", n)):
            for c in range(-6, 7):
                for d in range(-6, 7):
                    tested += 1
                    if (X * X + c * X + d) % m == 0:
                        found.append({"space": "Q2", "id": f"Q2:{xname}^2{c:+d}*{xname}{d:+d}=0 mod {mname}",
                                      "relation": f"{xname}^2 {c:+d}·{xname} {d:+d} ≡ 0 (mod {mname})"})
    # Q3 linear relations in the GLV basis
    for e1, e2 in itertools.product(range(-3, 4), repeat=2):
        tested += 2
        if k1.a2 == e1 * k1.a1 + e2 * k1.b1:
            found.append({"space": "Q3", "id": f"Q3:a2={e1:+d}*a1{e2:+d}*b1", "relation": f"a2 = {e1}·a1 {e2:+d}·b1"})
        if k1.b2 == e1 * k1.a1 + e2 * k1.b1:
            found.append({"space": "Q3", "id": f"Q3:b2={e1:+d}*a1{e2:+d}*b1", "relation": f"b2 = {e1}·a1 {e2:+d}·b1"})
    # Q4 negative control: λ ≡ c0 + c1 β + c2 β² (mod n) for small c_i and also mod p — must not hold (except trivialities)
    q4 = 0
    for c0, c1, c2 in itertools.product(range(-6, 7), repeat=3):
        tested += 2
        if (c0 + c1 * beta + c2 * beta * beta - lam) % n == 0 or (c0 + c1 * beta + c2 * beta * beta - lam) % p == 0:
            q4 += 1
            found.append({"space": "Q4", "id": f"Q4:lambda={c0}+{c1}beta+{c2}beta^2", "relation": "cross-modulus relation (unexpected)"})
    # Q5 negative control: Gy² ≡ Gx³ + 7 is the curve equation (expected, trivial) — exclude it; look for
    # small-coefficient linear relations c1·Gx + c2·Gy + c3·beta + c4 ≡ 0 (mod p), |c_i| ≤ 6, not all zero
    q5 = 0
    for c in itertools.product(range(-6, 7), repeat=4):
        if not any(c):
            continue
        tested += 1
        if (c[0] * k1.Gx + c[1] * k1.Gy + c[2] * beta + c[3]) % p == 0:
            q5 += 1
            found.append({"space": "Q5", "id": f"Q5:{c}", "relation": "small linear relation among Gx, Gy, beta (unexpected)"})
    return found, tested, q4, q5


def main():
    t0 = time.time()
    found, tested, q4, q5 = search()
    # normalise ids and classify
    out = []
    for h in found:
        hid = h["id"].replace("+0*", "+0*")
        known = None
        for k, v in KNOWN.items():
            # tolerant matching on structure
            kk = k.replace("+", "").replace("-", "").replace("*", "")
            hh = hid.replace("+", "").replace("-", "").replace("*", "")
            if kk == hh:
                known = v
        if known is None and h["space"] == "Q1" and h["relation"].endswith("= n") and all(v in h["relation"] for v in ["^2"]):
            known = EISENSTEIN_NORM_FAMILY
        h["classification"] = "known_reproduced" if known else ("unexpected_investigate" if h["space"] in ("Q4", "Q5") else "found_unclassified")
        h["explanation"] = known
        out.append(h)
    # The expected count of chance hits for Q5: 13^4 trials, each true with prob 1/p ≈ 2^-256 → 0.
    closed_ok = all(getattr(k1, k) == v for k, v in CLOSED_FORM.items())
    meas = {"candidates_tested": tested, "identities_found": len(out), "found": out,
            "closed_form_glv_basis": {k: str(v) for k, v in CLOSED_FORM.items()}, "closed_form_matches_libsecp256k1": closed_ok,
            "closed_form_statement": "a1 = (t−f)/2 − 1, b1 = −f, a2 = (t+f)/2 − 1, b2 = a1 with 4p = t² + 3f²; lattice = ideal (π̄−1)Z[ω]",
            "negative_controls": {"Q4_cross_modulus_hits": q4, "Q5_small_linear_hits": q5,
                                  "expected_chance_hits": tested / p},
            "unexplained": [h for h in out if h["classification"] != "known_reproduced"]}
    record("identity_search_secp256k1_constants", {"spaces": ["Q1", "Q2", "Q3", "Q4", "Q5"], "coeff_bound": 6}, meas,
           out=str(ROOT / "results/automated_discovery/identity_search.json"), started=t0)
    print(f"tested {tested} candidates, found {len(out)} identities; negative controls Q4={q4} Q5={q5}; closed-form GLV basis matches: {closed_ok}")
    for h in out:
        print(f"  [{h['classification']}] {h['relation']}  — {h['explanation'] or ''}")


if __name__ == "__main__":
    main()