1"""Analysis of the ECDLP studies: statistics with bootstrap CIs, empirical complexity fits, figures and LaTeX tables.23Inputs : results/synthetic_ecdlp/{scaling_generic,scaling_j0,sweep_r,parallel,tmto}_<host>.json4Outputs: results/synthetic_ecdlp/analysis.json, figures/ecdlp_*.pdf (+ copies in paper/figures), paper/tables/ecdlp_*.tex5"""6import glob, json, math, pathlib, re, sys7import numpy as np8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt11ROOT = pathlib.Path(__file__).resolve().parents[2]12R = ROOT / "results/synthetic_ecdlp"; F = ROOT / "figures"; PF = ROOT / "paper/figures"; T = ROOT / "paper/tables"13for d in (F, PF, T):14 d.mkdir(parents=True, exist_ok=True)15plt.rcParams.update({"font.family": "serif", "font.size": 9, "mathtext.fontset": "stix", "axes.grid": True, "grid.alpha": 0.3})16rng = np.random.default_rng(34)17ANALYSIS = {}18MACROS = {}192021def load(study):22 files = sorted(glob.glob(str(R / f"{study}_*.json")))23 files = [f for f in files if ".chunk" not in f and ".raw" not in f]24 if not files:25 return None, None26 d = json.loads(pathlib.Path(files[0]).read_text())27 runs = [r for r in d["measurements"]["runs"] if r]28 for r in runs:29 r["config"] = r["label"].split("/")[1]; r["inst"] = r["label"].split("/")[0]30 r["bits"] = int(re.match(r"[gj](\d+)b", r["inst"]).group(1))31 r["n"] = int(r["params"]["instance"]["n"]) if "instance" in r["params"] else int(r["params"]["n"])32 return d, runs333435def boot_mean(x, B=2000):36 x = np.asarray(x, float)37 if len(x) < 2:38 return float(x.mean()), float(x.mean()), float(x.mean())39 means = rng.choice(x, size=(B, len(x)), replace=True).mean(axis=1)40 return float(x.mean()), float(np.percentile(means, 2.5)), float(np.percentile(means, 97.5))414243def boot_ratio(a, b, B=2000):44 """Bootstrap CI for mean(a)/mean(b) with paired resampling when lengths match."""45 a, b = np.asarray(a, float), np.asarray(b, float)46 if len(a) == len(b):47 idx = rng.integers(0, len(a), size=(B, len(a)))48 r = a[idx].mean(axis=1) / b[idx].mean(axis=1)49 else:50 r = rng.choice(a, size=(B, len(a))).mean(axis=1) / rng.choice(b, size=(B, len(b))).mean(axis=1)51 return float(a.mean() / b.mean()), float(np.percentile(r, 2.5)), float(np.percentile(r, 97.5))525354def fit_power(ns, ys, with_log=False):55 """Least squares on log y = log a + α log n (+ β log log2 n). Returns (a, α, β, r²) and bootstrap CI for α."""56 ns, ys = np.asarray(ns, float), np.asarray(ys, float)57 X = np.column_stack([np.ones_like(ns), np.log(ns)] + ([np.log(np.log2(ns))] if with_log else []))58 y = np.log(ys)59 coef, *_ = np.linalg.lstsq(X, y, rcond=None)60 pred = X @ coef; r2 = 1 - ((y - pred) ** 2).sum() / ((y - y.mean()) ** 2).sum()61 alphas = []62 for _ in range(1000):63 idx = rng.integers(0, len(ns), len(ns))64 c, *_ = np.linalg.lstsq(X[idx], y[idx], rcond=None); alphas.append(c[1])65 return {"a": float(math.exp(coef[0])), "alpha": float(coef[1]), "beta": float(coef[2]) if with_log else 0.0, "r2": float(r2),66 "alpha_ci": [float(np.percentile(alphas, 2.5)), float(np.percentile(alphas, 97.5))]}676869def summarise(runs, key="group_ops"):70 out = {}71 for cfg in sorted({r["config"] for r in runs}):72 out[cfg] = {}73 for bits in sorted({r["bits"] for r in runs if r["config"] == cfg}):74 rs = [r for r in runs if r["config"] == cfg and r["bits"] == bits]75 n = float(np.mean([r["n"] for r in rs]))76 ops = [r[key] for r in rs]; mean, lo, hi = boot_mean(ops)77 its = [r.get("iterations", 0) for r in rs]78 out[cfg][bits] = {"n_mean": n, "count": len(rs), "ops_mean": mean, "ops_ci": [lo, hi], "ops_median": float(np.median(ops)),79 "ops_sd": float(np.std(ops)), "ops_over_sqrt_n": mean / math.sqrt(n), "iter_mean": float(np.mean(its)),80 "iter_over_sqrt_n": float(np.mean(its)) / math.sqrt(n), "wall_mean": float(np.mean([r["wall_s"] for r in rs])),81 "cpu_mean": float(np.mean([r["cpu_s"] for r in rs])), "rss_max_mb": max(r["peak_rss_bytes"] for r in rs) / 2**20,82 "table_entries_mean": float(np.mean([r.get("table_entries", 0) or 0 for r in rs])),83 "dp_count_mean": float(np.mean([r.get("dp_count", 0) or 0 for r in rs])),84 "fruitless_mean": float(np.mean([r.get("fruitless_cycles", 0) or 0 for r in rs])),85 "ok": all(r["ok"] for r in rs)}86 return out878889LABELS = {"bsgs": "BSGS", "rho_floyd": "rho, Floyd", "rho_brent": "rho, Brent", "rho_dp": "rho, DP ($r=20$)", "rho_dp_r128": "rho, DP ($r=128$)",90 "rho_dp_neg": "rho, DP + negation ($r=128$)", "kangaroo": "kangaroo (full interval)",91 "rho_dp_aut3": r"rho, DP + $\phi$-orbits", "rho_dp_neg_aut3": r"rho, DP + negation + $\phi$ ($r=128$)", "rho_dp_neg_aut3_r20": r"rho, DP + negation + $\phi$ ($r=20$)"}92THEORY = {"bsgs": lambda n: 2 * math.sqrt(n), "rho_floyd": lambda n: 3 * math.sqrt(math.pi * n / 2), "rho_brent": lambda n: 1.0 * math.sqrt(math.pi * n / 2) * 1.15,93 "rho_dp": lambda n: math.sqrt(math.pi * n / 2), "rho_dp_r128": lambda n: math.sqrt(math.pi * n / 2), "rho_dp_neg": lambda n: math.sqrt(math.pi * n / 4),94 "kangaroo": lambda n: 2 * math.sqrt(n), "rho_dp_aut3": lambda n: math.sqrt(math.pi * n / 6), "rho_dp_neg_aut3": lambda n: math.sqrt(math.pi * n / 12),95 "rho_dp_neg_aut3_r20": lambda n: math.sqrt(math.pi * n / 12)}969798def analyse_scaling():99 d, runs = load("scaling_generic")100 if not runs:101 return102 S = summarise(runs); ANALYSIS["scaling_generic"] = {"host": d["host"], "summary": S}103 # fits per config on group_ops104 fits = {}105 for cfg, per in S.items():106 ns = [v["n_mean"] for v in per.values()]; ys = [v["ops_mean"] for v in per.values()]107 if len(ns) >= 3:108 big = [(v["n_mean"], v["ops_mean"], v["iter_mean"]) for b, v in per.items() if b >= 32]109 fits[cfg] = {"power": fit_power(ns, ys), "power_log": fit_power(ns, ys, True),110 "power_ge32": fit_power([t[0] for t in big], [t[1] for t in big]) if len(big) >= 3 else None,111 "power_iter_ge32": fit_power([t[0] for t in big], [t[2] for t in big]) if len(big) >= 3 and all(t[2] > 0 for t in big) else None}112 ANALYSIS["scaling_generic"]["fits"] = fits113 # figure: runtime scaling (ops and wall) log-log with theory lines114 fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.3))115 for cfg, per in S.items():116 bits = sorted(per); ns = [per[b]["n_mean"] for b in bits]117 axes[0].errorbar(ns, [per[b]["ops_mean"] for b in bits], yerr=[[per[b]["ops_mean"] - per[b]["ops_ci"][0] for b in bits], [per[b]["ops_ci"][1] - per[b]["ops_mean"] for b in bits]], marker="o", ms=3, lw=1, capsize=2, label=LABELS.get(cfg, cfg))118 axes[1].plot(ns, [per[b]["wall_mean"] for b in bits], marker="o", ms=3, lw=1, label=LABELS.get(cfg, cfg))119 nn = np.logspace(math.log10(2**20), math.log10(2**56), 50)120 axes[0].plot(nn, np.sqrt(np.pi * nn / 2), "k--", lw=0.8, label=r"$\sqrt{\pi n/2}$")121 axes[0].plot(nn, 2 * np.sqrt(nn), "k:", lw=0.8, label=r"$2\sqrt{n}$")122 for ax in axes:123 ax.set_xscale("log", base=2); ax.set_yscale("log"); ax.set_xlabel("group order $n$")124 axes[0].set_ylabel("group operations (mean, 95% CI)"); axes[1].set_ylabel("wall-clock seconds (1 thread)")125 axes[0].legend(fontsize=6, ncol=1, loc="lower right"); fig.tight_layout()126 fig.savefig(F / "ecdlp_scaling.pdf"); fig.savefig(PF / "ecdlp_scaling.pdf"); plt.close(fig)127 # figure: normalised ops/√n with theory constants128 fig, ax = plt.subplots(figsize=(6.4, 3.2))129 for cfg, per in S.items():130 bits = sorted(per)131 ax.plot(bits, [per[b]["ops_over_sqrt_n"] for b in bits], marker="o", ms=3, lw=1, label=LABELS.get(cfg, cfg))132 for c, lab in [(2, r"$2$ (BSGS, kangaroo)"), (math.sqrt(math.pi / 2), r"$\sqrt{\pi/2}$ (rho)"), (math.sqrt(math.pi / 4), r"$\sqrt{\pi/4}$ (negation)")]:133 ax.axhline(c, color="gray", lw=0.6, ls="--"); ax.text(bits[-1] + 0.3, c, lab, fontsize=6, va="center")134 ax.set_xlabel("group size (bits)"); ax.set_ylabel(r"group operations / $\sqrt{n}$"); ax.set_yscale("log"); ax.legend(fontsize=6, ncol=2)135 fig.tight_layout(); fig.savefig(F / "ecdlp_constants.pdf"); fig.savefig(PF / "ecdlp_constants.pdf"); plt.close(fig)136 # table137 rows = [r"\begin{tabular}{@{}llrrrrrr@{}}", r"\toprule Algorithm & Bits & Inst. & Group ops (mean) & 95\% CI & ops$/\sqrt n$ & Theory & Wall (s) \\ \midrule"]138 for cfg, per in S.items():139 for i, (b, v) in enumerate(sorted(per.items())):140 if b not in (20, 28, 36, 44, 48, 52, 56):141 continue142 th = THEORY[cfg](v["n_mean"]) / math.sqrt(v["n_mean"])143 rows.append(f"{LABELS.get(cfg, cfg) if i == 0 else ''} & {b} & {v['count']} & {v['ops_mean']:,.0f} & [{v['ops_ci'][0]:,.0f}, {v['ops_ci'][1]:,.0f}] & {v['ops_over_sqrt_n']:.2f} & {th:.2f} & {v['wall_mean']:.3f} \\\\".replace(",", r"\,"))144 rows += [r"\bottomrule", r"\end{tabular}"]145 (T / "ecdlp_scaling.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")146 rows = [r"\begin{tabular}{@{}lrrrrrr@{}}", r"\toprule Algorithm & $\hat\alpha$ (all sizes, ops) & 95\% CI & $R^2$ & $\hat\alpha$ ($\ge 32$ bits, ops) & 95\% CI & $\hat\alpha$ ($\ge 32$ bits, iterations) \\ \midrule"]147 for cfg, ff in fits.items():148 p, pg, pi = ff["power"], ff["power_ge32"], ff["power_iter_ge32"]149 g = f"{pg['alpha']:.3f} & [{pg['alpha_ci'][0]:.3f}, {pg['alpha_ci'][1]:.3f}]" if pg else "--- & ---"150 it = f"{pi['alpha']:.3f} [{pi['alpha_ci'][0]:.3f}, {pi['alpha_ci'][1]:.3f}]" if pi else "---"151 rows.append(f"{LABELS.get(cfg, cfg)} & {p['alpha']:.3f} & [{p['alpha_ci'][0]:.3f}, {p['alpha_ci'][1]:.3f}] & {p['r2']:.4f} & {g} & {it} \\\\")152 rows += [r"\bottomrule", r"\end{tabular}"]153 (T / "ecdlp_fits.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")154 for cfg, ff in fits.items():155 k = cfg.replace('_', '').replace('128', 'big')156 MACROS[f"alpha{k}"] = f"{ff['power']['alpha']:.3f}"157 MACROS[f"alphaLo{k}"] = f"{ff['power']['alpha_ci'][0]:.3f}"; MACROS[f"alphaHi{k}"] = f"{ff['power']['alpha_ci'][1]:.3f}"158 if ff["power_ge32"]:159 MACROS[f"alphaG{k}"] = f"{ff['power_ge32']['alpha']:.3f}"; MACROS[f"alphaGLo{k}"] = f"{ff['power_ge32']['alpha_ci'][0]:.3f}"; MACROS[f"alphaGHi{k}"] = f"{ff['power_ge32']['alpha_ci'][1]:.3f}"160 alphas = [ff["power_ge32"]["alpha"] for ff in fits.values() if ff["power_ge32"]]161 MACROS["alphaGMin"] = f"{min(alphas):.3f}"; MACROS["alphaGMax"] = f"{max(alphas):.3f}"162 big = max(S["rho_dp"]); MACROS["scalingMaxBits"] = str(big); MACROS["scalingInstances"] = str(S["rho_dp"][big]["count"])163 MACROS["rhoDpConstBig"] = f"{S['rho_dp'][big]['ops_over_sqrt_n']:.2f}"; MACROS["rhoDpIterConstBig"] = f"{S['rho_dp'][big]['iter_over_sqrt_n']:.2f}"164 MACROS["rhoNegIterConstBig"] = f"{S['rho_dp_neg'][big]['iter_over_sqrt_n']:.2f}"165 MACROS["bsgsConstBig"] = f"{S['bsgs'][max(S['bsgs'])]['ops_over_sqrt_n']:.2f}"; MACROS["kangarooConstBig"] = f"{S['kangaroo'][max(S['kangaroo'])]['ops_over_sqrt_n']:.2f}"166 MACROS["floydConstBig"] = f"{S['rho_floyd'][max(S['rho_floyd'])]['ops_over_sqrt_n']:.2f}"; MACROS["brentConstBig"] = f"{S['rho_brent'][max(S['rho_brent'])]['ops_over_sqrt_n']:.2f}"167 MACROS["scalingAllOk"] = "all" if all(v["ok"] for per in S.values() for v in per.values()) else "NOT ALL"168 MACROS["scalingHost"] = d["host"]["cpu"]169 # speedup neg vs plain (paired by instance)170 sp = {}171 for b in sorted(S["rho_dp_neg"]):172 a = [r["iterations"] for r in runs if r["config"] == "rho_dp_r128" and r["bits"] == b]173 c = [r["iterations"] for r in runs if r["config"] == "rho_dp_neg" and r["bits"] == b]174 sp[b] = boot_ratio(a, c)175 ANALYSIS["scaling_generic"]["neg_speedup_iterations"] = sp176 MACROS["negSpeedupBig"] = f"{sp[big][0]:.2f}"; MACROS["negSpeedupBigLo"] = f"{sp[big][1]:.2f}"; MACROS["negSpeedupBigHi"] = f"{sp[big][2]:.2f}"177 PB = 44178 allsp = [r["iterations"] for r in runs if r["config"] == "rho_dp_r128" and r["bits"] >= PB]; allc = [r["iterations"] for r in runs if r["config"] == "rho_dp_neg" and r["bits"] >= PB]179 pooled = boot_ratio(np.array(allsp) / np.array([math.sqrt(r["n"]) for r in runs if r["config"] == "rho_dp_r128" and r["bits"] >= PB]),180 np.array(allc) / np.array([math.sqrt(r["n"]) for r in runs if r["config"] == "rho_dp_neg" and r["bits"] >= PB]))181 MACROS["pooledMinBits"] = str(PB); MACROS["pooledN"] = str(len(allsp))182 # plain rho constant pooled over ≥ PB bits, with CI, against √(π/2)183 c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] in ("rho_dp", "rho_dp_r128") and r["bits"] >= PB]184 m_, lo_, hi_ = boot_mean(c); MACROS["rhoConstPooled"] = f"{m_:.3f}"; MACROS["rhoConstPooledLo"] = f"{lo_:.3f}"; MACROS["rhoConstPooledHi"] = f"{hi_:.3f}"; MACROS["rhoConstPooledN"] = str(len(c))185 c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == "rho_dp_neg" and r["bits"] >= PB]186 m_, lo_, hi_ = boot_mean(c); MACROS["negConstPooled"] = f"{m_:.3f}"; MACROS["negConstPooledLo"] = f"{lo_:.3f}"; MACROS["negConstPooledHi"] = f"{hi_:.3f}"187 c = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == "bsgs" and r["bits"] >= 36]188 m_, lo_, hi_ = boot_mean(c); MACROS["bsgsConstPooled"] = f"{m_:.3f}"; MACROS["bsgsConstPooledLo"] = f"{lo_:.3f}"; MACROS["bsgsConstPooledHi"] = f"{hi_:.3f}"189 c = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == "kangaroo" and r["bits"] >= PB]190 m_, lo_, hi_ = boot_mean(c); MACROS["kangConstPooled"] = f"{m_:.3f}"; MACROS["kangConstPooledLo"] = f"{lo_:.3f}"; MACROS["kangConstPooledHi"] = f"{hi_:.3f}"191 ANALYSIS["scaling_generic"]["pooled_constants"] = {k: v for k, v in MACROS.items() if "Pooled" in k}192 MACROS["negSpeedupPooled"] = f"{pooled[0]:.2f}"; MACROS["negSpeedupPooledLo"] = f"{pooled[1]:.2f}"; MACROS["negSpeedupPooledHi"] = f"{pooled[2]:.2f}"193 ANALYSIS["scaling_generic"]["neg_speedup_pooled_ge36"] = pooled194195196def analyse_j0():197 d, runs = load("scaling_j0")198 if not runs:199 return200 S = summarise(runs); ANALYSIS["scaling_j0"] = {"host": d["host"], "summary": S}201 base = "rho_dp"202 ratios = {}203 for cfg in ("rho_dp_neg", "rho_dp_aut3", "rho_dp_neg_aut3", "rho_dp_neg_aut3_r20"):204 ratios[cfg] = {}205 for b in sorted(S[cfg]):206 a = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] == b]207 c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == cfg and r["bits"] == b]208 ratios[cfg][b] = boot_ratio(a, c)209 PB = 44210 a = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] >= PB]211 c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == cfg and r["bits"] >= PB]212 ratios[cfg]["pooled_ge36"] = boot_ratio(a, c) # key kept for compatibility; pooled over bits ≥ PB213 ratios[cfg]["pooled_min_bits"] = PB; ratios[cfg]["pooled_n"] = len(a)214 m_, lo_, hi_ = boot_mean(c); ratios[cfg]["const_pooled"] = [m_, lo_, hi_]215 ratios[cfg]["cpu_ns_per_iter_ge40"] = float(np.mean([1e9 * r["cpu_s"] / max(r["iterations"], 1) for r in runs if r["config"] == cfg and r["bits"] >= 40]))216 # wall-clock ratio (what the user actually gains)217 a = [r["wall_s"] for r in runs if r["config"] == base and r["bits"] >= 40]; c = [r["wall_s"] for r in runs if r["config"] == cfg and r["bits"] >= 40]218 ratios[cfg]["wall_pooled_ge40"] = boot_ratio(a, c)219 a = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] >= PB]; c = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == cfg and r["bits"] >= PB]220 ratios[cfg]["ops_pooled_ge36"] = boot_ratio(a, c)221 a = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] >= PB]222 m_, lo_, hi_ = boot_mean(a); ratios["_plain_const_pooled"] = [m_, lo_, hi_]223 ratios["_plain_cpu_ns_per_iter_ge40"] = float(np.mean([1e9 * r["cpu_s"] / max(r["iterations"], 1) for r in runs if r["config"] == base and r["bits"] >= 40]))224 ANALYSIS["scaling_j0"]["iteration_speedup_vs_plain"] = ratios225 fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.2))226 for cfg, lab, th in [("rho_dp_neg", "negation", math.sqrt(2)), ("rho_dp_aut3", r"$\phi$-orbits", math.sqrt(3)), ("rho_dp_neg_aut3", r"negation + $\phi$ ($r=128$)", math.sqrt(6)), ("rho_dp_neg_aut3_r20", r"negation + $\phi$ ($r=20$)", math.sqrt(6))]:227 bits = sorted(b for b in ratios[cfg] if isinstance(b, int))228 m = [ratios[cfg][b][0] for b in bits]; lo = [ratios[cfg][b][0] - ratios[cfg][b][1] for b in bits]; hi = [ratios[cfg][b][2] - ratios[cfg][b][0] for b in bits]229 l = axes[0].errorbar(bits, m, yerr=[lo, hi], marker="o", ms=3, lw=1, capsize=2, label=lab)230 axes[0].axhline(th, color=l[0].get_color(), lw=0.6, ls="--")231 axes[0].set_xlabel("group size (bits)"); axes[0].set_ylabel("iteration speed-up vs plain rho (paired, 95% CI)"); axes[0].legend(fontsize=6)232 axes[0].set_title(r"dashed: $\sqrt{2}$, $\sqrt{3}$, $\sqrt{6}$", fontsize=8)233 for cfg in ("rho_dp_neg", "rho_dp_aut3", "rho_dp_neg_aut3", "rho_dp_neg_aut3_r20"):234 bits = sorted(S[cfg]); axes[1].plot(bits, [S[cfg][b]["fruitless_mean"] / max(S[cfg][b]["iter_mean"], 1) * 1e3 for b in bits], marker="o", ms=3, lw=1, label=LABELS.get(cfg, cfg))235 axes[1].set_xlabel("group size (bits)"); axes[1].set_ylabel("fruitless cycles per 1000 iterations"); axes[1].legend(fontsize=6)236 fig.tight_layout(); fig.savefig(F / "ecdlp_endomorphism.pdf"); fig.savefig(PF / "ecdlp_endomorphism.pdf"); plt.close(fig)237 rows = [r"\begin{tabular}{@{}lrrrrr@{}}", r"\toprule Classes & Theory & Iterations$/\sqrt n$ (95\% CI) & Iteration speed-up (paired, 95\% CI) & ns per iteration & Wall-clock speed-up \\ \midrule",238 f"none (plain rho, $r=128$) & $\\sqrt{{\\pi/2}} = 1.253$ & {ratios['_plain_const_pooled'][0]:.3f} [{ratios['_plain_const_pooled'][1]:.3f}, {ratios['_plain_const_pooled'][2]:.3f}] & 1 & {ratios['_plain_cpu_ns_per_iter_ge40']:.0f} & 1 \\\\"]239 for cfg, lab, th in [("rho_dp_neg", r"$\{\pm P\}$ (negation)", r"$\sqrt2 = 1.414$"), ("rho_dp_aut3", r"$\{P,\phi P,\phi^2P\}$", r"$\sqrt3 = 1.732$"), ("rho_dp_neg_aut3", r"order-6 orbits, $r=128$", r"$\sqrt6 = 2.449$"), ("rho_dp_neg_aut3_r20", r"order-6 orbits, $r=20$", r"$\sqrt6 = 2.449$")]:240 p, o, w, cp = ratios[cfg]["pooled_ge36"], ratios[cfg]["ops_pooled_ge36"], ratios[cfg]["wall_pooled_ge40"], ratios[cfg]["const_pooled"]241 rows.append(f"{lab} & {th} & {cp[0]:.3f} [{cp[1]:.3f}, {cp[2]:.3f}] & {p[0]:.3f} [{p[1]:.3f}, {p[2]:.3f}] & {ratios[cfg]['cpu_ns_per_iter_ge40']:.0f} & {w[0]:.2f} [{w[1]:.2f}, {w[2]:.2f}] \\\\")242 rows += [r"\bottomrule", r"\end{tabular}"]243 (T / "ecdlp_j0.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")244 for cfg, key in [("rho_dp_neg", "Neg"), ("rho_dp_aut3", "Aut"), ("rho_dp_neg_aut3", "NegAut"), ("rho_dp_neg_aut3_r20", "NegAutRtwenty")]:245 p = ratios[cfg]["pooled_ge36"]; MACROS[f"jzero{key}Speedup"] = f"{p[0]:.2f}"; MACROS[f"jzero{key}SpeedupLo"] = f"{p[1]:.2f}"; MACROS[f"jzero{key}SpeedupHi"] = f"{p[2]:.2f}"246 w = ratios[cfg]["wall_pooled_ge40"]; MACROS[f"jzero{key}Wall"] = f"{w[0]:.2f}"; MACROS[f"jzero{key}WallLo"] = f"{w[1]:.2f}"; MACROS[f"jzero{key}WallHi"] = f"{w[2]:.2f}"247 MACROS[f"jzero{key}NsIter"] = f"{ratios[cfg]['cpu_ns_per_iter_ge40']:.0f}"; MACROS[f"jzero{key}Const"] = f"{ratios[cfg]['const_pooled'][0]:.3f}"248 MACROS["jzeroPlainNsIter"] = f"{ratios['_plain_cpu_ns_per_iter_ge40']:.0f}"; MACROS["jzeroPlainConst"] = f"{ratios['_plain_const_pooled'][0]:.3f}"249 MACROS["jzeroPooledN"] = str(ratios["rho_dp_neg"]["pooled_n"]); MACROS["jzeroPooledMinBits"] = str(ratios["rho_dp_neg"]["pooled_min_bits"])250 MACROS["jzeroMaxBits"] = str(max(S["rho_dp"])); MACROS["jzeroHost"] = d["host"]["cpu"]251 MACROS["jzeroAllOk"] = "all" if all(v["ok"] for per in S.values() for v in per.values()) else "NOT ALL"252 b = max(S["rho_dp_neg_aut3"]); MACROS["jzeroFruitlessPerK"] = f"{S['rho_dp_neg_aut3'][b]['fruitless_mean'] / max(S['rho_dp_neg_aut3'][b]['iter_mean'], 1) * 1e3:.2f}"253 MACROS["jzeroFruitlessPerKRtwenty"] = f"{S['rho_dp_neg_aut3_r20'][b]['fruitless_mean'] / max(S['rho_dp_neg_aut3_r20'][b]['iter_mean'], 1) * 1e3:.2f}"254255256def analyse_sweep():257 d, runs = load("sweep_r")258 if not runs:259 return260 res = {}261 for bits in sorted({r["bits"] for r in runs}):262 res[bits] = {}263 for cfg in sorted({r["config"] for r in runs}):264 rs = [r for r in runs if r["bits"] == bits and r["config"] == cfg]265 if not rs:266 continue267 x = [r["iterations"] / math.sqrt(r["n"]) for r in rs]; m, lo, hi = boot_mean(x)268 res[bits][cfg] = {"iter_over_sqrt_n": m, "ci": [lo, hi], "count": len(rs), "ops_over_sqrt_n": float(np.mean([r["group_ops"] / math.sqrt(r["n"]) for r in rs]))}269 ANALYSIS["sweep_r"] = {"host": d["host"], "results": res}270 fig, ax = plt.subplots(figsize=(6.4, 3.2))271 rs_ = [2, 3, 4, 8, 16, 20, 32, 64, 128]272 for bits, mk in zip(sorted(res), ("o", "s")):273 for mixed, ls in ((0, "-"), (2, "--"), (4, ":")):274 xs, ys, lo, hi = [], [], [], []275 for r in rs_:276 k = f"r{r}_m{mixed}"277 if k in res[bits]:278 xs.append(r); ys.append(res[bits][k]["iter_over_sqrt_n"]); lo.append(ys[-1] - res[bits][k]["ci"][0]); hi.append(res[bits][k]["ci"][1] - ys[-1])279 if xs:280 ax.errorbar(xs, ys, yerr=[lo, hi], marker=mk, ms=3, lw=1, ls=ls, capsize=2, label=f"{bits} bits, mixed={mixed or 'none'}")281 ax.axhline(math.sqrt(math.pi / 2), color="k", lw=0.8, ls="--"); ax.text(2.1, math.sqrt(math.pi / 2) * 1.02, r"random mapping $\sqrt{\pi/2}$", fontsize=7)282 ax.set_xscale("log", base=2); ax.set_xlabel("number of partitions $r$ of the adding walk"); ax.set_ylabel(r"iterations / $\sqrt{n}$ (mean, 95% CI)"); ax.legend(fontsize=6, ncol=2)283 fig.tight_layout(); fig.savefig(F / "ecdlp_sweep_r.pdf"); fig.savefig(PF / "ecdlp_sweep_r.pdf"); plt.close(fig)284 rows = [r"\begin{tabular}{@{}lrrrrrr@{}}", r"\toprule $r$ & \multicolumn{3}{c}{36 bits: iterations$/\sqrt n$} & \multicolumn{3}{c}{40 bits: iterations$/\sqrt n$} \\", r" & plain & mixed 2 & mixed 4 & plain & mixed 2 & mixed 4 \\ \midrule"]285 for r in rs_:286 cells = []287 for bits in (36, 40):288 for mixed in (0, 2, 4):289 k = f"r{r}_m{mixed}"; v = res.get(bits, {}).get(k)290 cells.append(f"{v['iter_over_sqrt_n']:.3f}" if v else "---")291 rows.append(f"{r} & " + " & ".join(cells) + r" \\")292 rows += [r"\midrule", r"random-mapping constant $\sqrt{\pi/2}$ & 1.253 & & & 1.253 & & \\", r"\bottomrule", r"\end{tabular}"]293 (T / "ecdlp_sweep_r.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")294 best = min(((bits, cfg, v["iter_over_sqrt_n"], v["ci"]) for bits in res for cfg, v in res[bits].items() if cfg != "r2_m0"), key=lambda t: t[2])295 MACROS["sweepBestCfg"] = best[1].replace("_", r"\_"); MACROS["sweepBestConst"] = f"{best[2]:.3f}"; MACROS["sweepBestLo"] = f"{best[3][0]:.3f}"; MACROS["sweepBestHi"] = f"{best[3][1]:.3f}"296 MACROS["sweepRtwo"] = f"{res[40]['r2_m0']['iter_over_sqrt_n']:.2f}" if "r2_m0" in res.get(40, {}) else "n/a"297 MACROS["sweepRthree"] = f"{res[40]['r3_m0']['iter_over_sqrt_n']:.2f}" if "r3_m0" in res.get(40, {}) else "n/a"298 MACROS["sweepRtwenty"] = f"{res[40]['r20_m0']['iter_over_sqrt_n']:.3f}"; MACROS["sweepRhundred"] = f"{res[40]['r128_m0']['iter_over_sqrt_n']:.3f}"299 MACROS["sweepConfigs"] = str(len({r["config"] for r in runs})); MACROS["sweepInstances"] = str(len({r["inst"] for r in runs}))300301302def analyse_parallel():303 d, runs = load("parallel")304 if not runs:305 return306 res = {}307 for t in sorted({r["threads"] for r in runs}):308 rs = [r for r in runs if r["threads"] == t]309 res[t] = {"wall_mean": float(np.mean([r["wall_s"] for r in rs])), "ops_mean": float(np.mean([r["group_ops"] for r in rs])),310 "ops_per_wall_s": float(np.mean([r["group_ops"] / r["wall_s"] for r in rs])), "count": len(rs), "cpu_mean": float(np.mean([r["cpu_s"] for r in rs]))}311 base = res[1]312 for t, v in res.items():313 v["speedup_throughput"] = v["ops_per_wall_s"] / base["ops_per_wall_s"]; v["efficiency"] = v["speedup_throughput"] / t314 v["ops_ratio_to_1thread"] = v["ops_mean"] / base["ops_mean"]315 ANALYSIS["parallel"] = {"host": d["host"], "results": res}316 fig, axes = plt.subplots(1, 2, figsize=(7.0, 3.0))317 ts = sorted(res); axes[0].plot(ts, [res[t]["speedup_throughput"] for t in ts], marker="o", ms=3, label="measured throughput speed-up"); axes[0].plot(ts, ts, "k--", lw=0.8, label="linear")318 axes[0].set_xlabel("threads"); axes[0].set_ylabel("throughput speed-up vs 1 thread"); axes[0].legend(fontsize=7)319 axes[1].plot(ts, [res[t]["ops_ratio_to_1thread"] for t in ts], marker="s", ms=3, color="C3"); axes[1].axhline(1, color="k", lw=0.8, ls="--")320 axes[1].set_xlabel("threads"); axes[1].set_ylabel("total group ops, relative to 1 thread")321 fig.tight_layout(); fig.savefig(F / "ecdlp_parallel.pdf"); fig.savefig(PF / "ecdlp_parallel.pdf"); plt.close(fig)322 rows = [r"\begin{tabular}{@{}rrrrrr@{}}", r"\toprule Threads & Wall (s) & Group ops ($\times 10^6$) & Mops/s & Speed-up & Efficiency \\ \midrule"]323 for t in ts:324 v = res[t]; rows.append(f"{t} & {v['wall_mean']:.2f} & {v['ops_mean']/1e6:.1f} & {v['ops_per_wall_s']/1e6:.0f} & {v['speedup_throughput']:.2f} & {v['efficiency']:.2f} \\\\")325 rows += [r"\bottomrule", r"\end{tabular}"]326 (T / "ecdlp_parallel.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")327 tmax = max(ts); MACROS["parMaxThreads"] = str(tmax); MACROS["parSpeedupMax"] = f"{res[tmax]['speedup_throughput']:.1f}"; MACROS["parEffMax"] = f"{res[tmax]['efficiency']:.2f}"328 t16 = 16 if 16 in res else tmax; MACROS["parSpeedupSixteen"] = f"{res[t16]['speedup_throughput']:.1f}"; MACROS["parEffSixteen"] = f"{res[t16]['efficiency']:.2f}"329 MACROS["parOpsRatioMax"] = f"{res[tmax]['ops_ratio_to_1thread']:.2f}"; MACROS["parHost"] = d["host"]["cpu"]; MACROS["parCores"] = str(d["host"]["cores"])330331332def analyse_tmto():333 d, runs = load("tmto")334 if not runs:335 return336 bs = {}337 for r in runs:338 if r["config"].startswith("bsgs_m"):339 m = int(r["config"][6:]); bs.setdefault(m, []).append(r)340 rh = {}341 for r in runs:342 if r["config"].startswith("rho_dp"):343 dpb = int(r["config"][6:]); rh.setdefault(dpb, []).append(r)344 B = {m: {"ops_mean": float(np.mean([r["group_ops"] for r in rs])), "table_entries": float(np.mean([r["table_entries"] for r in rs])), "table_mb": float(np.mean([r.get("table_bytes", 0) for r in rs])) / 2**20,345 "rss_mb": float(np.mean([r["peak_rss_bytes"] for r in rs])) / 2**20, "wall": float(np.mean([r["wall_s"] for r in rs])), "n": float(np.mean([r["n"] for r in rs])), "pred": float(np.mean([m + r["n"] / (2 * m) for r in rs]))} for m, rs in bs.items()}346 Rr = {dpb: {"ops_mean": float(np.mean([r["group_ops"] for r in rs])), "ops_ci": boot_mean([r["group_ops"] for r in rs])[1:], "dp_count": float(np.mean([r["dp_count"] for r in rs])), "dp_mb": float(np.mean([r.get("dp_bytes", 0) for r in rs])) / 2**20,347 "rss_mb": float(np.mean([r["peak_rss_bytes"] for r in rs])) / 2**20, "wall": float(np.mean([r["wall_s"] for r in rs])), "n": float(np.mean([r["n"] for r in rs])), "iter": float(np.mean([r["iterations"] for r in rs]))} for dpb, rs in rh.items()}348 ANALYSIS["tmto"] = {"host": d["host"], "bsgs": B, "rho": Rr}349 fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.2))350 ms = sorted(B); axes[0].plot([B[m]["table_entries"] for m in ms], [B[m]["ops_mean"] for m in ms], "o-", ms=3, label="BSGS measured (36 bits)")351 axes[0].plot([B[m]["table_entries"] for m in ms], [B[m]["pred"] for m in ms], "k--", lw=0.8, label=r"$m + n/(2m)$ (table keyed by $x$: $\pm j$ resolved)")352 axes[0].set_xscale("log", base=2); axes[0].set_yscale("log"); axes[0].set_xlabel("table entries (memory)"); axes[0].set_ylabel("group operations"); axes[0].legend(fontsize=7)353 ds = sorted(Rr); axes[1].errorbar([max(Rr[k]["dp_count"], 1) for k in ds], [Rr[k]["ops_mean"] for k in ds], yerr=[[Rr[k]["ops_mean"] - Rr[k]["ops_ci"][0] for k in ds], [Rr[k]["ops_ci"][1] - Rr[k]["ops_mean"] for k in ds]], fmt="s-", ms=3, capsize=2, label="rho, DP (40 bits)")354 for k in ds:355 axes[1].annotate(f"$2^{{-{k}}}$", (max(Rr[k]["dp_count"], 1), Rr[k]["ops_mean"]), xytext=(3, 3), textcoords="offset points", fontsize=6)356 axes[1].axhline(math.sqrt(math.pi * Rr[ds[0]]["n"] / 2), color="k", lw=0.8, ls="--", label=r"$\sqrt{\pi n/2}$")357 axes[1].set_xscale("log", base=2); axes[1].set_xlabel("stored distinguished points (memory)"); axes[1].set_ylabel("group operations"); axes[1].legend(fontsize=7)358 fig.tight_layout(); fig.savefig(F / "ecdlp_tmto.pdf"); fig.savefig(PF / "ecdlp_tmto.pdf"); plt.close(fig)359 rows = [r"\begin{tabular}{@{}lrrrrr@{}}", r"\toprule Algorithm / parameter & Memory (entries) & Memory (MB) & Group ops & Predicted & Wall (s) \\ \midrule"]360 for m in ms:361 v = B[m]; rows.append(f"BSGS, $m=2^{{{int(math.log2(m))}}}$ (36 bits) & {v['table_entries']:,.0f} & {v['table_mb']:.1f} & {v['ops_mean']:,.0f} & {v['pred']:,.0f} & {v['wall']:.3f} \\\\".replace(",", r"\,"))362 rows.append(r"\midrule")363 for k in ds:364 v = Rr[k]; rows.append(f"rho, DP density $2^{{-{k}}}$ (40 bits) & {v['dp_count']:,.0f} & {v['dp_mb']:.2f} & {v['ops_mean']:,.0f} & {math.sqrt(math.pi * v['n'] / 2):,.0f} & {v['wall']:.3f} \\\\".replace(",", r"\,"))365 rows += [r"\bottomrule", r"\end{tabular}"]366 (T / "ecdlp_tmto.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")367 MACROS["tmtoBsgsMinM"] = f"2^{{{int(math.log2(min(ms)))}}}"; MACROS["tmtoBsgsMaxM"] = f"2^{{{int(math.log2(max(ms)))}}}"368 dev = max(abs(B[m]["ops_mean"] / B[m]["pred"] - 1) for m in ms); MACROS["tmtoBsgsMaxDev"] = f"{100*dev:.1f}"369 MACROS["tmtoRhoDpMin"] = str(min(ds)); MACROS["tmtoRhoDpMax"] = str(max(ds))370 MACROS["tmtoRhoOpsRange"] = f"{min(Rr[k]['ops_mean'] for k in ds)/1e6:.2f}--{max(Rr[k]['ops_mean'] for k in ds)/1e6:.2f}"371372373def main():374 analyse_scaling(); analyse_j0(); analyse_sweep(); analyse_parallel(); analyse_tmto()375 (R / "analysis.json").write_text(json.dumps(ANALYSIS, indent=1, default=str))376 (T / "numbers_ecdlp.tex").write_text("% generated by experiments/synthetic_ecdlp/analyze_ecdlp.py\n" + "\n".join(f"\\newcommand{{\\{k}}}{{{v}}}" for k, v in MACROS.items()) + "\n")377 print("studies analysed:", list(ANALYSIS), "| macros:", len(MACROS))378379380if __name__ == "__main__":381 main()382