1"""Short Weierstrass elliptic curves over prime fields with exact integer arithmetic.23Points are affine tuples (x, y) or None for the point at infinity. Every group operation is counted in an4OpCounter so that experiments can report operation counts alongside wall-clock time.5"""6from __future__ import annotations7import random8from dataclasses import dataclass, field91011class OpCounter:12 __slots__ = ("add", "dbl", "inv", "mul")1314 def __init__(self):15 self.add = self.dbl = self.inv = self.mul = 01617 def reset(self):18 self.add = self.dbl = self.inv = self.mul = 01920 def group_ops(self) -> int:21 return self.add + self.dbl2223 def as_dict(self) -> dict:24 return {"add": self.add, "dbl": self.dbl, "inv": self.inv, "mul": self.mul}252627def legendre(a: int, p: int) -> int:28 return pow(a % p, (p - 1) // 2, p)293031def sqrt_mod(a: int, p: int):32 """Tonelli–Shanks square root modulo an odd prime p, or None if a is a non-residue."""33 a %= p34 if a == 0:35 return 036 if legendre(a, p) != 1:37 return None38 if p % 4 == 3:39 return pow(a, (p + 1) // 4, p)40 q, s = p - 1, 041 while q % 2 == 0:42 q //= 243 s += 144 z = 245 while legendre(z, p) != p - 1:46 z += 147 m, c, t, r = s, pow(z, q, p), pow(a, q, p), pow(a, (q + 1) // 2, p)48 while t != 1:49 i, t2 = 0, t50 while t2 != 1:51 t2 = t2 * t2 % p52 i += 153 b = pow(c, 1 << (m - i - 1), p)54 m, c, t, r = i, b * b % p, t * b * b % p, r * b % p55 return r565758@dataclass59class Curve:60 p: int61 a: int62 b: int63 n: int | None = None # order of the subgroup generated by G (prime in our instances)64 h: int = 1 # cofactor65 G: tuple | None = None66 name: str = ""67 ops: OpCounter = field(default_factory=OpCounter, repr=False, compare=False)6869 # ---- curve invariants -------------------------------------------------------------------70 def discriminant(self) -> int:71 return (-16 * (4 * pow(self.a, 3, self.p) + 27 * pow(self.b, 2, self.p))) % self.p7273 def j_invariant(self) -> int:74 a3 = 4 * pow(self.a, 3, self.p)75 den = (a3 + 27 * pow(self.b, 2, self.p)) % self.p76 return 1728 * a3 * pow(den, -1, self.p) % self.p # j = 1728 · 4a³ / (4a³ + 27b²)7778 def is_on_curve(self, P) -> bool:79 if P is None:80 return True81 x, y = P82 return (y * y - (x * x * x + self.a * x + self.b)) % self.p == 08384 # ---- group law ------------------------------------------------------------------------------85 def neg(self, P):86 return None if P is None else (P[0], (-P[1]) % self.p)8788 def add(self, P, Q):89 if P is None:90 return Q91 if Q is None:92 return P93 p = self.p94 x1, y1 = P95 x2, y2 = Q96 if x1 == x2:97 if (y1 + y2) % p == 0:98 return None99 return self.double(P)100 self.ops.add += 1101 self.ops.inv += 1102 lam = (y2 - y1) * pow(x2 - x1, -1, p) % p103 x3 = (lam * lam - x1 - x2) % p104 y3 = (lam * (x1 - x3) - y1) % p105 return (x3, y3)106107 def double(self, P):108 if P is None:109 return None110 p = self.p111 x1, y1 = P112 if y1 == 0:113 return None114 self.ops.dbl += 1115 self.ops.inv += 1116 lam = (3 * x1 * x1 + self.a) * pow(2 * y1, -1, p) % p117 x3 = (lam * lam - 2 * x1) % p118 y3 = (lam * (x1 - x3) - y1) % p119 return (x3, y3)120121 def sub(self, P, Q):122 return self.add(P, self.neg(Q))123124 def mul(self, k: int, P):125 """Left-to-right double-and-add. k may be any integer (negative allowed)."""126 if P is None or k == 0:127 return None128 if k < 0:129 k, P = -k, self.neg(P)130 R = None131 for bit in bin(k)[2:]:132 R = self.double(R)133 if bit == "1":134 R = self.add(R, P)135 return R136137 def mul_wnaf(self, k: int, P, w: int = 4):138 """Signed-window (wNAF) scalar multiplication; same result as mul, fewer additions."""139 if P is None or k == 0:140 return None141 if k < 0:142 k, P = -k, self.neg(P)143 # precompute odd multiples 1P, 3P, ..., (2^{w-1}-1)P144 twoP = self.double(P)145 table = [P]146 for _ in range((1 << (w - 2)) - 1):147 table.append(self.add(table[-1], twoP))148 naf = []149 while k > 0:150 if k & 1:151 d = k % (1 << w)152 if d >= (1 << (w - 1)):153 d -= 1 << w154 naf.append(d)155 k -= d156 else:157 naf.append(0)158 k >>= 1159 R = None160 for d in reversed(naf):161 R = self.double(R)162 if d > 0:163 R = self.add(R, table[(d - 1) // 2])164 elif d < 0:165 R = self.add(R, self.neg(table[(-d - 1) // 2]))166 return R167168 # ---- utilities ----------------------------------------------------------------------------169 def random_point(self, rng: random.Random):170 while True:171 x = rng.randrange(self.p)172 rhs = (x * x * x + self.a * x + self.b) % self.p173 y = sqrt_mod(rhs, self.p)174 if y is not None:175 return (x, y if rng.random() < 0.5 else (-y) % self.p)176177 def lift_x(self, x: int, even_y: bool = True):178 rhs = (x * x * x + self.a * x + self.b) % self.p179 y = sqrt_mod(rhs, self.p)180 if y is None:181 return None182 if (y % 2 == 0) != even_y:183 y = (-y) % self.p184 return (x, y)185186 def order_of_point(self, P, bound: int | None = None) -> int:187 """Naive order (only for tiny curves); uses n if known."""188 if self.n is not None and self.mul(self.n, P) is None:189 # smallest divisor of n·h that kills P190 N = self.n * self.h191 for d in sorted(_divisors(N)):192 if self.mul(d, P) is None:193 return d194 k, R = 1, P195 while R is not None:196 R = self.add(R, P)197 k += 1198 if bound and k > bound:199 raise ValueError("order bound exceeded")200 return k201202 def count_points_naive(self) -> int:203 """#E(F_p) by Legendre symbols — only for p < ~10^6."""204 total = 1 # infinity205 for x in range(self.p):206 rhs = (x * x * x + self.a * x + self.b) % self.p207 if rhs == 0:208 total += 1209 elif legendre(rhs, self.p) == 1:210 total += 2211 return total212213214def _divisors(N: int):215 from sympy import divisors216 return divisors(N)217218219# ---- synthetic curve generation via PARI/GP (SEA point counting) ---------------------------------------220def gp_ellcard_batch(curves: list[tuple[int, int, int]]) -> list[int]:221 """Orders #E(F_p) for a list of (p, a, b) using PARI/GP's ellcard (SEA)."""222 import subprocess223 lines = ["default(parisize, 64000000);"]224 for (p, a, b) in curves:225 lines.append(f"print(ellcard(ellinit([{a},{b}],{p})));")226 out = subprocess.run(["gp", "-q", "-f"], input="\n".join(lines) + "\n", capture_output=True, text=True,227 timeout=600)228 if out.returncode != 0:229 raise RuntimeError(out.stderr)230 return [int(s) for s in out.stdout.split()]231232233def random_prime_order_curve(bits: int, rng: random.Random, j0: bool = False, batch: int = 32) -> Curve:234 """Random curve y² = x³ + ax + b over a random `bits`-bit prime p with prime group order.235236 If j0 is True, a = 0 and p ≡ 1 (mod 3) so that the curve has the order-6 automorphism group of secp256k1.237 """238 from sympy import isprime, randprime239 while True:240 if j0:241 while True:242 p = randprime(1 << (bits - 1), 1 << bits)243 if p % 3 == 1:244 break245 else:246 p = randprime(1 << (bits - 1), 1 << bits)247 cands = []248 for _ in range(batch):249 a = 0 if j0 else rng.randrange(p)250 b = rng.randrange(1, p)251 if (4 * a ** 3 + 27 * b ** 2) % p == 0:252 continue253 cands.append((p, a, b))254 orders = gp_ellcard_batch(cands)255 for (p, a, b), N in zip(cands, orders):256 if isprime(N) and N != p:257 E = Curve(p, a, b, n=N, h=1, name=f"synthetic-{bits}b{'-j0' if j0 else ''}")258 E.G = E.random_point(rng)259 assert E.mul(N, E.G) is None260 return E261262263def cube_roots_of_unity(m: int, rng: random.Random | None = None) -> list[int]:264 """The two primitive cube roots of unity modulo a prime m ≡ 1 (mod 3)."""265 assert m % 3 == 1266 rng = rng or random.Random(0)267 while True:268 g = rng.randrange(2, m)269 beta = pow(g, (m - 1) // 3, m)270 if beta != 1:271 return sorted([beta, beta * beta % m])272273274def endomorphism_lambda(E: Curve, beta: int) -> int:275 """Return λ with (βx, y) = [λ]P on E (prime order n ≡ 1 mod 3, a = 0)."""276 assert E.a == 0 and E.n is not None and E.n % 3 == 1277 lams = cube_roots_of_unity(E.n)278 P = E.G279 phiP = (beta * P[0] % E.p, P[1])280 for lam in lams:281 if E.mul(lam, P) == phiP:282 return lam283 raise ValueError("no λ matches β")284