"""Exact verification of the secp256k1 domain parameters and algebraic structure, cross-checked with PARI/GP.

Writes results/baselines/secp256k1_facts.json. Every fact is computed with exact integer arithmetic twice:
once in pure Python (src/wp34) and once in PARI/GP (independent implementation), and the two must agree.
"""
import argparse, json, subprocess, 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


def gp(script: str) -> list[str]:
    out = subprocess.run(["gp", "-q", "-f"], input="default(parisize,128000000);\n" + script + "\n",
                         capture_output=True, text=True, timeout=1200)
    if out.returncode != 0:
        raise RuntimeError(out.stderr)
    return out.stdout.split("\n")


def main():
    ap = argparse.ArgumentParser(); ap.add_argument("--out", default=str(ROOT / "results/baselines/secp256k1_facts.json"))
    args = ap.parse_args()
    t0 = time.time()
    facts = k1.verify_parameters()
    p, n, beta, lam = k1.p, k1.n, k1.beta, k1.lam
    twist = facts["twist_order"]
    script = f"""
p={p}; n={n}; E=ellinit([0,7],p); G=[{k1.Gx},{k1.Gy}];
print("p_is_prime=",isprime(p)); print("n_is_prime=",isprime(n)); print("ellcard=",ellcard(E));
print("G_on_curve=",ellisoncurve(E,G)); print("nG_is_infinity=",ellmul(E,G,n)==[0]);
print("trace_t=",ellap(E)); print("j_invariant=",lift(E.j));
print("beta_cubed_is_1_mod_p=",Mod({beta},p)^3==1); print("lambda_cubed_is_1_mod_n=",Mod({lam},n)^3==1);
print("phi_G_equals_lambda_G=",ellmul(E,G,{lam})==[Mod({beta},p)*G[1], G[2]]);
tw={twist}; f=factor(tw,10^7); print("twist_small_primes=",f[,1]~); print("twist_small_exponents=",f[,2]~);
c=tw; for(i=1,#f[,1], if(f[i,1]<10^7, c=c/f[i,1]^f[i,2])); print("twist_cofactor_is_prime=",isprime(c)); print("twist_cofactor_bits=",#binary(c));
print("embedding_degree=",znorder(Mod(p,n)));
print("n_minus_1_small_primes=",factor(n-1,10^7)[,1]~); print("p_minus_1_small_primes=",factor(p-1,10^7)[,1]~);
t=ellap(E); d=t^2-4*p; print("frobenius_discriminant=",d); print("disc_is_minus3_times_square=",issquare(-d/3));
z=sqrtn(Mod(1,p),3,&zz); print("cube_roots_unity_mod_p=",[lift(zz),lift(zz^2)]);
y=sqrtn(Mod(1,n),3,&yy); print("cube_roots_unity_mod_n=",[lift(yy),lift(yy^2)]);
print("kronecker_7_p=",kronecker(7,p)); print("num_roots_x3_plus_7=",#polrootsmod(x^3+7, p));
"""
    g = {}
    for line in gp(script):
        if "=" in line:
            key, val = line.split("=", 1)
            g[key.strip()] = val.strip()
    def I(k): return int(g[k])
    def B(k): return g[k] == "1"
    gp_facts = {
        "p_is_prime": B("p_is_prime"), "n_is_prime": B("n_is_prime"), "ellcard": I("ellcard"), "G_on_curve": B("G_on_curve"),
        "nG_is_infinity": B("nG_is_infinity"), "trace_t": I("trace_t"), "j_invariant": I("j_invariant"),
        "beta_cubed_is_1_mod_p": B("beta_cubed_is_1_mod_p"), "lambda_cubed_is_1_mod_n": B("lambda_cubed_is_1_mod_n"),
        "phi_G_equals_lambda_G": B("phi_G_equals_lambda_G"),
        "twist_small_primes": g["twist_small_primes"], "twist_small_exponents": g["twist_small_exponents"],
        "twist_cofactor_is_prime": B("twist_cofactor_is_prime"), "twist_cofactor_bits": I("twist_cofactor_bits"),
        "embedding_degree": I("embedding_degree"),
        "n_minus_1_small_primes": g["n_minus_1_small_primes"], "p_minus_1_small_primes": g["p_minus_1_small_primes"],
        "frobenius_discriminant": I("frobenius_discriminant"), "disc_is_minus3_times_square": B("disc_is_minus3_times_square"),
        "cube_roots_unity_mod_p": [int(v) for v in g["cube_roots_unity_mod_p"].strip("[]").split(",")],
        "cube_roots_unity_mod_n": [int(v) for v in g["cube_roots_unity_mod_n"].strip("[]").split(",")],
        "kronecker_7_p": I("kronecker_7_p"), "num_roots_x3_plus_7": I("num_roots_x3_plus_7"),
    }
    agree = {
        "p_prime": facts["p_is_prime"] == gp_facts["p_is_prime"] == True,
        "n_prime": facts["n_is_prime"] == gp_facts["n_is_prime"] == True,
        "order_equals_n": gp_facts["ellcard"] == n,
        "trace": facts["trace_t"] == gp_facts["trace_t"],
        "j": facts["j_invariant"] == gp_facts["j_invariant"] == 0,
        "phi_lambda": facts["phi_G_equals_lambda_G"] == gp_facts["phi_G_equals_lambda_G"] == True,
        "twist_cofactor_prime": facts["twist_cofactor_is_prime"] == gp_facts["twist_cofactor_is_prime"],
        "twist_cofactor_bits": facts["twist_cofactor_bits"] == gp_facts["twist_cofactor_bits"],
        "embedding_degree_gt_100": facts["embedding_degree_le_100"] is None and gp_facts["embedding_degree"] > 100,
        "beta_in_gp_roots": beta in gp_facts["cube_roots_unity_mod_p"],
        "lambda_in_gp_roots": lam in gp_facts["cube_roots_unity_mod_n"],
        "x0_points": facts["x0_points_exist"] == (gp_facts["kronecker_7_p"] == 1),
        "cm_discriminant": facts["disc_is_minus3_times_square"] == gp_facts["disc_is_minus3_times_square"] == True,
        "frobenius_discriminant": facts["frobenius_discriminant"] == gp_facts["frobenius_discriminant"],
    }
    # The GLV basis can be re-derived from (n, λ) by the extended Euclidean algorithm (GLV 2001, §4).
    r0, r1 = n, lam
    s0, s1 = 1, 0
    t0_, t1 = 0, 1
    seq = []
    while r1 != 0:
        q = r0 // r1
        r0, r1 = r1, r0 - q * r1
        s0, s1 = s1, s0 - q * s1
        t0_, t1 = t1, t0_ - q * t1
        seq.append((r0, t0_))
    import math
    sq = math.isqrt(n)
    idx = next(i for i, (r, t) in enumerate(seq) if r < sq)
    # candidate vectors (r_i, -t_i) and (r_{i+1}, -t_{i+1}) or (r_{i-1}, -t_{i-1})
    cands = [(seq[idx - 1][0], -seq[idx - 1][1]), (seq[idx][0], -seq[idx][1]), (seq[idx + 1][0], -seq[idx + 1][1])]
    pub = {(k1.a1, k1.b1), (k1.a2, k1.b2)}
    derived = set()
    for a, b in cands:
        for sgn in (1, -1):
            if (sgn * a, sgn * b) in pub:
                derived.add((sgn * a, sgn * b))
    facts["glv_basis_rederived_from_euclid"] = len(derived) == 2
    facts["glv_euclid_candidates"] = [(str(a), str(b)) for a, b in cands]
    measurements = {"python": {k: (str(v) if isinstance(v, int) and abs(v) > 2**63 else v) for k, v in facts.items()},
                    "pari_gp": {k: (str(v) if isinstance(v, int) and abs(v) > 2**63 else v) for k, v in gp_facts.items()},
                    "agreement": agree, "all_agree": all(agree.values())}
    record("secp256k1_parameter_verification", {"gp_version": subprocess.run(["gp", "--version-short"], capture_output=True, text=True).stdout.strip() or "2.19"},
           measurements, out=args.out, started=t0)
    print(json.dumps(agree, indent=1)); print("ALL AGREE:", all(agree.values()))
    print("twist cofactor bits", facts["twist_cofactor_bits"], "prime", facts["twist_cofactor_is_prime"], "| small factors", facts["twist_small_factors"])
    print("embedding degree (gp):", gp_facts["embedding_degree"] if gp_facts["embedding_degree"] < 10**6 else "huge")
    print("trace t =", facts["trace_t"], "conductor f =", facts["conductor_f"])
    print("GLV basis re-derived from Euclid:", facts["glv_basis_rederived_from_euclid"])
    print("x=0 points exist:", facts["x0_points_exist"], "| #roots of x^3+7:", gp_facts["num_roots_x3_plus_7"])
    print("n-1 small factors", facts["n_minus_1_small_factors"], "| p-1 small factors", facts["p_minus_1_small_factors"])


if __name__ == "__main__":
    main()