1"""Exact verification of the secp256k1 domain parameters and algebraic structure, cross-checked with PARI/GP.23Writes results/baselines/secp256k1_facts.json. Every fact is computed with exact integer arithmetic twice:4once in pure Python (src/wp34) and once in PARI/GP (independent implementation), and the two must agree.5"""6import argparse, json, subprocess, sys, time, pathlib7ROOT = pathlib.Path(__file__).resolve().parents[2]8sys.path.insert(0, str(ROOT / "src"))9from wp34 import secp256k1 as k110from wp34.provenance import record111213def gp(script: str) -> list[str]:14 out = subprocess.run(["gp", "-q", "-f"], input="default(parisize,128000000);\n" + script + "\n",15 capture_output=True, text=True, timeout=1200)16 if out.returncode != 0:17 raise RuntimeError(out.stderr)18 return out.stdout.split("\n")192021def main():22 ap = argparse.ArgumentParser(); ap.add_argument("--out", default=str(ROOT / "results/baselines/secp256k1_facts.json"))23 args = ap.parse_args()24 t0 = time.time()25 facts = k1.verify_parameters()26 p, n, beta, lam = k1.p, k1.n, k1.beta, k1.lam27 twist = facts["twist_order"]28 script = f"""29p={p}; n={n}; E=ellinit([0,7],p); G=[{k1.Gx},{k1.Gy}];30print("p_is_prime=",isprime(p)); print("n_is_prime=",isprime(n)); print("ellcard=",ellcard(E));31print("G_on_curve=",ellisoncurve(E,G)); print("nG_is_infinity=",ellmul(E,G,n)==[0]);32print("trace_t=",ellap(E)); print("j_invariant=",lift(E.j));33print("beta_cubed_is_1_mod_p=",Mod({beta},p)^3==1); print("lambda_cubed_is_1_mod_n=",Mod({lam},n)^3==1);34print("phi_G_equals_lambda_G=",ellmul(E,G,{lam})==[Mod({beta},p)*G[1], G[2]]);35tw={twist}; f=factor(tw,10^7); print("twist_small_primes=",f[,1]~); print("twist_small_exponents=",f[,2]~);36c=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));37print("embedding_degree=",znorder(Mod(p,n)));38print("n_minus_1_small_primes=",factor(n-1,10^7)[,1]~); print("p_minus_1_small_primes=",factor(p-1,10^7)[,1]~);39t=ellap(E); d=t^2-4*p; print("frobenius_discriminant=",d); print("disc_is_minus3_times_square=",issquare(-d/3));40z=sqrtn(Mod(1,p),3,&zz); print("cube_roots_unity_mod_p=",[lift(zz),lift(zz^2)]);41y=sqrtn(Mod(1,n),3,&yy); print("cube_roots_unity_mod_n=",[lift(yy),lift(yy^2)]);42print("kronecker_7_p=",kronecker(7,p)); print("num_roots_x3_plus_7=",#polrootsmod(x^3+7, p));43"""44 g = {}45 for line in gp(script):46 if "=" in line:47 key, val = line.split("=", 1)48 g[key.strip()] = val.strip()49 def I(k): return int(g[k])50 def B(k): return g[k] == "1"51 gp_facts = {52 "p_is_prime": B("p_is_prime"), "n_is_prime": B("n_is_prime"), "ellcard": I("ellcard"), "G_on_curve": B("G_on_curve"),53 "nG_is_infinity": B("nG_is_infinity"), "trace_t": I("trace_t"), "j_invariant": I("j_invariant"),54 "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"),55 "phi_G_equals_lambda_G": B("phi_G_equals_lambda_G"),56 "twist_small_primes": g["twist_small_primes"], "twist_small_exponents": g["twist_small_exponents"],57 "twist_cofactor_is_prime": B("twist_cofactor_is_prime"), "twist_cofactor_bits": I("twist_cofactor_bits"),58 "embedding_degree": I("embedding_degree"),59 "n_minus_1_small_primes": g["n_minus_1_small_primes"], "p_minus_1_small_primes": g["p_minus_1_small_primes"],60 "frobenius_discriminant": I("frobenius_discriminant"), "disc_is_minus3_times_square": B("disc_is_minus3_times_square"),61 "cube_roots_unity_mod_p": [int(v) for v in g["cube_roots_unity_mod_p"].strip("[]").split(",")],62 "cube_roots_unity_mod_n": [int(v) for v in g["cube_roots_unity_mod_n"].strip("[]").split(",")],63 "kronecker_7_p": I("kronecker_7_p"), "num_roots_x3_plus_7": I("num_roots_x3_plus_7"),64 }65 agree = {66 "p_prime": facts["p_is_prime"] == gp_facts["p_is_prime"] == True,67 "n_prime": facts["n_is_prime"] == gp_facts["n_is_prime"] == True,68 "order_equals_n": gp_facts["ellcard"] == n,69 "trace": facts["trace_t"] == gp_facts["trace_t"],70 "j": facts["j_invariant"] == gp_facts["j_invariant"] == 0,71 "phi_lambda": facts["phi_G_equals_lambda_G"] == gp_facts["phi_G_equals_lambda_G"] == True,72 "twist_cofactor_prime": facts["twist_cofactor_is_prime"] == gp_facts["twist_cofactor_is_prime"],73 "twist_cofactor_bits": facts["twist_cofactor_bits"] == gp_facts["twist_cofactor_bits"],74 "embedding_degree_gt_100": facts["embedding_degree_le_100"] is None and gp_facts["embedding_degree"] > 100,75 "beta_in_gp_roots": beta in gp_facts["cube_roots_unity_mod_p"],76 "lambda_in_gp_roots": lam in gp_facts["cube_roots_unity_mod_n"],77 "x0_points": facts["x0_points_exist"] == (gp_facts["kronecker_7_p"] == 1),78 "cm_discriminant": facts["disc_is_minus3_times_square"] == gp_facts["disc_is_minus3_times_square"] == True,79 "frobenius_discriminant": facts["frobenius_discriminant"] == gp_facts["frobenius_discriminant"],80 }81 # The GLV basis can be re-derived from (n, λ) by the extended Euclidean algorithm (GLV 2001, §4).82 r0, r1 = n, lam83 s0, s1 = 1, 084 t0_, t1 = 0, 185 seq = []86 while r1 != 0:87 q = r0 // r188 r0, r1 = r1, r0 - q * r189 s0, s1 = s1, s0 - q * s190 t0_, t1 = t1, t0_ - q * t191 seq.append((r0, t0_))92 import math93 sq = math.isqrt(n)94 idx = next(i for i, (r, t) in enumerate(seq) if r < sq)95 # candidate vectors (r_i, -t_i) and (r_{i+1}, -t_{i+1}) or (r_{i-1}, -t_{i-1})96 cands = [(seq[idx - 1][0], -seq[idx - 1][1]), (seq[idx][0], -seq[idx][1]), (seq[idx + 1][0], -seq[idx + 1][1])]97 pub = {(k1.a1, k1.b1), (k1.a2, k1.b2)}98 derived = set()99 for a, b in cands:100 for sgn in (1, -1):101 if (sgn * a, sgn * b) in pub:102 derived.add((sgn * a, sgn * b))103 facts["glv_basis_rederived_from_euclid"] = len(derived) == 2104 facts["glv_euclid_candidates"] = [(str(a), str(b)) for a, b in cands]105 measurements = {"python": {k: (str(v) if isinstance(v, int) and abs(v) > 2**63 else v) for k, v in facts.items()},106 "pari_gp": {k: (str(v) if isinstance(v, int) and abs(v) > 2**63 else v) for k, v in gp_facts.items()},107 "agreement": agree, "all_agree": all(agree.values())}108 record("secp256k1_parameter_verification", {"gp_version": subprocess.run(["gp", "--version-short"], capture_output=True, text=True).stdout.strip() or "2.19"},109 measurements, out=args.out, started=t0)110 print(json.dumps(agree, indent=1)); print("ALL AGREE:", all(agree.values()))111 print("twist cofactor bits", facts["twist_cofactor_bits"], "prime", facts["twist_cofactor_is_prime"], "| small factors", facts["twist_small_factors"])112 print("embedding degree (gp):", gp_facts["embedding_degree"] if gp_facts["embedding_degree"] < 10**6 else "huge")113 print("trace t =", facts["trace_t"], "conductor f =", facts["conductor_f"])114 print("GLV basis re-derived from Euclid:", facts["glv_basis_rederived_from_euclid"])115 print("x=0 points exist:", facts["x0_points_exist"], "| #roots of x^3+7:", gp_facts["num_roots_x3_plus_7"])116 print("n-1 small factors", facts["n_minus_1_small_factors"], "| p-1 small factors", facts["p_minus_1_small_factors"])117118119if __name__ == "__main__":120 main()121