1"""Index-calculus applicability probe on prime-field curves (Axis A, hypothesis on Semaev decomposition).23Question: with a factor base F = {P ∈ E(F_p) : x(P) < B}, what fraction of random points R decomposes as R = P1 + P24(P_i ∈ F), and what is the implied cost of a relation versus the generic √n cost?5Theory (Semaev 2004; Gaudry 2009; Diem 2011): over a prime field there is no subfield structure; a random point is a6sum of m factor-base points with probability ≈ |F|^m / (m!·n), so with |F| = n^{1/m}·c the relation-collection cost7is ≈ m!·n^{1−1/m}·(cost of decomposition), never below n^{1/2} for m = 2 and worse for m ≥ 3 once the decomposition8cost (solving S_{m+1} = 0) is included. We measure the m = 2 decomposition rate empirically on 16–28-bit curves.9"""10import argparse, math, random, sys, time, pathlib, statistics as st11ROOT = pathlib.Path(__file__).resolve().parents[2]12sys.path.insert(0, str(ROOT / "src"))13from wp34.ec import random_prime_order_curve, sqrt_mod14from wp34.provenance import record151617def run(bits, frac_exp, trials, rng):18 E = random_prime_order_curve(bits, rng)19 n, p = E.n, E.p20 B = int(p ** frac_exp) # factor base: x < B, |F| ≈ B (two points per x on average… half the x are on-curve ⇒ ≈ B)21 # enumerate the factor base (points with x < B); store by x for O(1) membership22 F = {}23 for x in range(B):24 rhs = (x * x * x + E.a * x + E.b) % p25 y = sqrt_mod(rhs, p)26 if y is not None:27 F[x] = y28 size = sum(1 if y == 0 else 2 for y in F.values())29 hits = 030 dec_cost = 031 for _ in range(trials):32 R = E.random_point(rng)33 # R = P1 + P2 ⇔ R − P1 ∈ F for some P1 ∈ F: try all P1 (|F| subtractions) — the naive decomposition34 found = False35 for x, y in F.items():36 for s in (y, p - y) if y else (0,):37 dec_cost += 138 S = E.sub(R, (x, s))39 if S is None or S[0] < B:40 found = True41 break42 if found:43 break44 hits += found45 rate = hits / trials46 pred = 1 - math.exp(-size * size / (2 * n)) # m = 2 prediction: P[R ∈ F + F] ≈ 1 − exp(−|F|²/(2n))47 return {"bits": bits, "p": str(p), "n": str(n), "B": B, "factor_base_size": size, "trials": trials, "hits": hits,48 "decomposition_rate": rate, "predicted_rate": pred,49 "relations_needed": size, "expected_cost_group_ops": (size / rate) * size if rate else None,50 "sqrt_n": math.sqrt(n), "cost_over_sqrt_n": ((size / rate) * size / math.sqrt(n)) if rate else None,51 # even with a FREE decomposition oracle: |F| relations at 1/rate trials each, then sparse linear algebra ~ |F|²52 "oracle_relation_cost": (size / rate) if rate else None,53 "oracle_relation_cost_over_sqrt_n": (size / rate / math.sqrt(n)) if rate else None,54 "linear_algebra_cost_over_sqrt_n": size * size / math.sqrt(n)}555657def main():58 ap = argparse.ArgumentParser(); ap.add_argument("--seed", type=int, default=34)59 ap.add_argument("--out", default=str(ROOT / "results/synthetic_ecdlp/index_calculus_probe.json"))60 a = ap.parse_args(); rng = random.Random(a.seed); t0 = time.time()61 rows = []62 for bits, trials in [(16, 400), (20, 300), (24, 200), (28, 120)]:63 for fe in (0.5, 0.6):64 if bits * fe > 15.5: # keep |F| ≤ ~50k so the naive decomposition stays affordable65 continue66 r = run(bits, fe, trials, rng); r["frac_exp"] = fe; rows.append(r)67 print(f"{bits}b |F|=p^{fe}: size {r['factor_base_size']:6d} rate {r['decomposition_rate']:.4f} (pred {r['predicted_rate']:.4f}) naive cost/√n ≈ {round(r['cost_over_sqrt_n'],1)} oracle-relations/√n ≈ {round(r['oracle_relation_cost_over_sqrt_n'],2)} linalg/√n ≈ {round(r['linear_algebra_cost_over_sqrt_n'],1)}")68 record("index_calculus_decomposition_probe", {"m": 2, "sizes": [16, 20, 24, 28]}, {"rows": rows}, seed=a.seed, out=a.out, started=t0)697071if __name__ == "__main__":72 main()73