1"""Reference (pure Python, exact) generic discrete-logarithm algorithms on an elliptic-curve group.23These are the *validation* implementations: slow but transparent. The Rust engine in4experiments/synthetic_ecdlp/rho_engine is the benchmarking implementation and must agree with these.56All functions take a Curve E with prime order n (E.n), a generator G and a target Q = dG, and return (d, stats).7"""8from __future__ import annotations9import math, random10from .ec import Curve111213def bsgs(E: Curve, G, Q, n: int | None = None, m: int | None = None):14 """Baby-step giant-step. Memory m (default ⌈√n⌉) entries; time ≈ m + n/m group operations."""15 n = n or E.n16 m = m or math.isqrt(n) + 117 E.ops.reset()18 table = {}19 P = None20 for j in range(m): # baby steps jG, j = 0..m-121 key = P22 table.setdefault(key, j)23 P = E.add(P, G)24 mG_neg = E.neg(P) # −mG25 R = Q26 for i in range(0, n // m + 2): # giant steps Q − i·mG27 if R in table:28 d = (i * m + table[R]) % n29 return d, {"ops": E.ops.as_dict(), "table_entries": len(table), "giant_steps": i + 1, "memory_entries": m}30 R = E.add(R, mG_neg)31 raise ValueError("BSGS failed (Q not in <G>?)")323334class RAddingWalk:35 """Teske's r-adding walk f(P) = P + M_{h(P)}, M_j = u_j G + v_j Q, with scalar bookkeeping (a, b): P = aG + bQ."""3637 def __init__(self, E: Curve, G, Q, n: int, r: int, rng: random.Random):38 self.E, self.G, self.Q, self.n, self.r = E, G, Q, n, r39 self.u = [rng.randrange(1, n) for _ in range(r)]40 self.v = [rng.randrange(1, n) for _ in range(r)]41 self.M = [E.add(E.mul(u, G), E.mul(v, Q)) for u, v in zip(self.u, self.v)]4243 def partition(self, P) -> int:44 # Hash of the x-coordinate; a multiplicative hash avoids correlation with the low bits of x.45 return (P[0] * 0x9E3779B97F4A7C15 >> 7) % self.r4647 def step(self, P, a, b):48 j = self.partition(P)49 return self.E.add(P, self.M[j]), (a + self.u[j]) % self.n, (b + self.v[j]) % self.n505152def pollard_rho(E: Curve, G, Q, n: int | None = None, r: int = 20, seed: int = 0, max_iter: int | None = None):53 """Single-walk Pollard rho with Floyd cycle detection (reference). Returns (d, stats)."""54 n = n or E.n55 rng = random.Random(seed)56 E.ops.reset()57 walk = RAddingWalk(E, G, Q, n, r, rng)58 a0, b0 = rng.randrange(n), rng.randrange(n)59 X = E.add(E.mul(a0, G), E.mul(b0, Q))60 if X is None:61 return (-a0 * pow(b0, -1, n)) % n, {"ops": E.ops.as_dict(), "iterations": 0}62 E.ops.reset()63 tort = (X, a0, b0)64 hare = (X, a0, b0)65 it = 066 while True:67 tort = walk.step(*tort)68 hare = walk.step(*walk.step(*hare))69 it += 170 if max_iter and it > max_iter:71 raise RuntimeError("rho exceeded max_iter")72 if tort[0] is None or hare[0] is None:73 # landed on infinity: restart with new seed74 return pollard_rho(E, G, Q, n, r, seed + 1_000_003, max_iter)75 if tort[0] == hare[0]:76 a1, b1 = tort[1], tort[2]77 a2, b2 = hare[1], hare[2]78 if (b1 - b2) % n == 0:79 return pollard_rho(E, G, Q, n, r, seed + 1_000_003, max_iter) # useless collision80 d = (a2 - a1) * pow(b1 - b2, -1, n) % n81 assert E.mul(d, G) == Q82 return d, {"ops": E.ops.as_dict(), "iterations": it, "floyd_steps": 3 * it}838485def pollard_rho_dp(E: Curve, G, Q, n: int | None = None, r: int = 20, dp_bits: int = 4, walks: int = 4,86 seed: int = 0, neg_map: bool = False, max_steps: int | None = None):87 """Distinguished-point (van Oorschot–Wiener) rho with `walks` sequential pseudo-parallel walks.8889 If neg_map is True the walk acts on classes {P, −P} (canonical representative: smaller y) and the usual90 2-cycle look-ahead escape is used; fruitless cycles longer than 2 are broken by a bounded-length check.91 """92 n = n or E.n93 rng = random.Random(seed)94 E.ops.reset()95 walk = RAddingWalk(E, G, Q, n, r, rng)96 mask = (1 << dp_bits) - 197 table = {}98 stats = {"fruitless_2cycles": 0, "fruitless_long": 0, "dp_count": 0, "iterations": 0, "restarts": 0}99100 def canon(P, a, b):101 if P[1] > E.p - P[1]:102 return (P[0], E.p - P[1]), (-a) % n, (-b) % n103 return P, a, b104105 def new_walk():106 a, b = rng.randrange(n), rng.randrange(n)107 P = E.add(E.mul(a, G), E.mul(b, Q))108 while P is None:109 a, b = rng.randrange(n), rng.randrange(n)110 P = E.add(E.mul(a, G), E.mul(b, Q))111 return canon(P, a, b) if neg_map else (P, a, b)112113 state = [new_walk() for _ in range(walks)]114 steps_since_check = [0] * walks115 anchor = [None] * walks116 while True:117 for w in range(walks):118 P, a, b = state[w]119 if neg_map:120 j = walk.partition(P)121 P2 = E.add(P, walk.M[j])122 if P2 is None:123 state[w] = new_walk(); stats["restarts"] += 1; continue124 P2, a2, b2 = canon(P2, (a + walk.u[j]) % n, (b + walk.v[j]) % n)125 # 2-cycle look-ahead (Bos–Kleinjung–Lenstra 2010): if f(P2) == P, use the next partition126 j2 = walk.partition(P2)127 P3 = E.add(P2, walk.M[j2])128 if P3 is not None and canon(P3, 0, 0)[0] == P:129 stats["fruitless_2cycles"] += 1130 j = (j + 1) % r131 P2 = E.add(P, walk.M[j])132 if P2 is None:133 state[w] = new_walk(); stats["restarts"] += 1; continue134 P2, a2, b2 = canon(P2, (a + walk.u[j]) % n, (b + walk.v[j]) % n)135 P, a, b = P2, a2, b2136 # long fruitless cycles: compare with an anchor point every 4096 steps137 steps_since_check[w] += 1138 if steps_since_check[w] % 4096 == 0:139 anchor[w] = P140 elif anchor[w] is not None and P == anchor[w]:141 stats["fruitless_long"] += 1142 # escape by doubling143 P = E.double(P); a, b = 2 * a % n, 2 * b % n144 P, a, b = canon(P, a, b)145 anchor[w] = None146 else:147 P, a, b = walk.step(P, a, b)148 if P is None:149 state[w] = new_walk(); stats["restarts"] += 1; continue150 stats["iterations"] += 1151 if max_steps and stats["iterations"] > max_steps:152 raise RuntimeError("exceeded max_steps")153 if P[0] & mask == 0:154 stats["dp_count"] += 1155 if P in table:156 a1, b1 = table[P]157 if (b1 - b) % n != 0:158 d = (a - a1) * pow(b1 - b, -1, n) % n159 assert E.mul(d, G) == Q160 stats["ops"] = E.ops.as_dict()161 stats["table_entries"] = len(table)162 return d, stats163 # same walk merged with itself (useless); restart this walk164 state[w] = new_walk(); stats["restarts"] += 1; continue165 table[P] = (a, b)166 state[w] = new_walk() # standard vOW: after a DP, start a fresh walk167 continue168 state[w] = (P, a, b)169170171def kangaroo(E: Curve, G, Q, lo: int, hi: int, n: int | None = None, seed: int = 0):172 """Pollard's kangaroo (lambda) method for d ∈ [lo, hi]: one tame and one wild kangaroo, DP-free173 'van Oorschot–Wiener trap' variant (tame kangaroo lays a trap, wild follows). Expected ≈ 2√w operations."""174 n = n or E.n175 w = hi - lo + 1176 rng = random.Random(seed)177 E.ops.reset()178 k = max(1, int(math.log2(math.sqrt(w)))) # jumps 2^0 .. 2^k, mean ≈ √w/2 (Pollard 2000)179 jumps = [1 << i for i in range(k + 1)]180 J = [E.mul(j, G) for j in jumps]181 mean = sum(jumps) / len(jumps)182 N = int(2 * math.sqrt(w) / mean * 4) + 10 # tame path length (van Oorschot–Wiener: ~2√w total)183184 salt = rng.getrandbits(64) | 1185186 def h(P):187 if P is None: # the walk may pass through the point at infinity when the interval is the whole group188 return 0189 return ((P[0] ^ salt) * 0x9E3779B97F4A7C15 >> 11) % len(jumps)190191 # tame kangaroo starts at the upper end of the interval and lays a trap at every point it visits192 t = hi193 T = E.mul(t, G)194 traps = {}195 for _ in range(N):196 traps[T] = t197 j = h(T)198 T = E.add(T, J[j]); t += jumps[j]199 traps[T] = t200 tame_dist = t - hi201 # wild kangaroo starts at Q = dG; it is caught once it lands on a trap. It passes the whole tame path202 # after at most w + tame_dist distance, so beyond that the attempt has failed.203 Wp, wdist = Q, 0204 it = 0205 while wdist <= w + tame_dist:206 if Wp in traps:207 d = (traps[Wp] - wdist) % n208 if E.mul(d, G) == Q:209 return d, {"ops": E.ops.as_dict(), "iterations": N + it, "tame_steps": N, "wild_steps": it}210 j = h(Wp)211 Wp = E.add(Wp, J[j]); wdist += jumps[j]212 it += 1213 # failure: retry with a different seed/longer tame path214 return kangaroo(E, G, Q, lo, hi, n, seed + 1)215216217def pohlig_hellman(E: Curve, G, Q, N: int, seed: int = 0):218 """Pohlig–Hellman for a point G of composite order N (smooth N makes this cheap); uses BSGS in each prime-power subgroup."""219 from sympy import factorint220 from sympy.ntheory.modular import crt221 rems, mods = [], []222 for q, e in factorint(N).items():223 qe = q ** e224 Gq = E.mul(N // qe, G)225 Qq = E.mul(N // qe, Q)226 # solve in the subgroup of order q^e digit by digit227 x = 0228 Gq1 = E.mul(q ** (e - 1), Gq) # order q229 for k in range(e):230 Qk = E.mul(q ** (e - 1 - k), E.sub(Qq, E.mul(x, Gq)))231 dk, _ = bsgs(E, Gq1, Qk, n=q)232 x += dk * q ** k233 rems.append(x); mods.append(qe)234 d, _ = crt(mods, rems)235 return int(d) % N, {}236