Skip to content

Instantly share code, notes, and snippets.

@n-WN

n-WN/test.py

Created March 14, 2026 08:24
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save n-WN/86d7ce992676f5869ce265b01bc451be to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save n-WN/86d7ce992676f5869ce265b01bc451be to your computer and use it in GitHub Desktop.
Download ZIP
AliCTF 2026 Final - RSA factorization via genus-2 Jacobian GCD leak
Raw
test.py
#!/usr/bin/env python3
import argparse
import time
from math import gcd, isqrt
"""Attack the special RSA where q = ((p - 1)^2 + 4x^2) / 8 and p = x^2 + 2y^2.
The key observation behind the attack is that the explicit divisor
D = (X^2 + 3X + 2, w(X + 1)), w^2 = -30
has order dividing 2q modulo the prime p generated by test.py, so [2n]D is
the identity modulo p and a gcd-leak appears during Cantor arithmetic over Z/nZ.
"""
KNOWN_P = int(
"46402297179264908470708471623273486751352579716103704797281944801460358524678"
"42279018921962833172524978958523306002263641947191072311981468164777376340163"
"09248109423490734066891943610586766143148559064501553780842031247144580440009"
"38194474983786306069139448770032043723951601909376213368723549683354239163507"
)
KNOWN_Q = int(
"26914664793910201015089600129623785390294811149557238224416199450647622361540"
"61458290133058210358061211802203736734517014228513818848859863320201028662594"
"08636909405263960880006579031003383880970259284258277447905002584560808016151"
"95599654382682078087103707560460633655990368809688841580499721704151332137274"
"51938932866495598613452778552658798768008138970942831931260605440957375646498"
"93934369707231143720201745545288105250504997847404937401018339096895509727111"
"86962060225378214303250799105314406275478888857499845714130299510525486044637"
"9794254941560890382535289658268169964153224826572350999211110389811740960617"
)
KNOWN_R = int(
"24388034157995205768053575087784566752894503365408073207022062557399294688638"
"49871187552042732895106626518198027167886195098589690365603086574529209472200"
"96949588170313491329321346278296852048084345607744818834743398321194013425847"
"36079181920031140709129205873116058110573979367198811758488627834078439211237"
"80070866680332283061603290014343206557668754257839454399733908519690656633727"
"38942851736654217461175512542145828887781883388886484361840892076528161284480"
"83334698430779030641077014553115044418997193362038721689196980116230439648408"
"256003405800149003946832004060056815330449886786417953912090190244491671682413"
)
KNOWN_N = int(
"30458271324341532949475846017727469451149884244455456848418344623613451533435"
"79307435382746063628125152906089520106401768106957001062097608543296273342600"
"73494232932631354400410633128749140076171158983190907973593575053497854417048"
"45070155389191120194965141592685269642219235223032729325019714071736443084339"
"80477948134183047844592974766663701461571247677079016119150849583378377004941"
"01778736872707088950191066536845315590172909690731773340519248676713609734650"
"33186389749487953940534414613347699119330329207595957401343988936792722602584"
"65824747012697595376626703975988882515595764957311907919027365749197002820887"
"81656589748596542720171740310760955369300856300367678388501209155989865070239"
"46443201683467311037291475027254181164572349694445942115684287264474876679889"
"31325516339110486866743445852066095669998920438968216485232424791922223861252"
"34242699387282215677275222361951085429715508244131809462847100374697050628398"
"79709527194763646717676930379663192565578694000082841476832027375526271828242"
"91133445659115219846651737385850400768144769901213446461201922926263027143565"
"52210558900333779120514161371252950169057659365777899851550703510494190464197"
"66493585403363271907314371735096265003807865198468046761716803922538443995949"
"23125523614125873231662726372159350593472047446355770877880057916630068295202"
"17920732447287933284467288077074320117165209921696169342156202641657329841597"
"67065660140195565607999242430922731963562167060407187425962779655725641005967"
"8199505555708003269487619386599472185153902440223444128863861928158832935247"
)
class FactorFound(RuntimeError):
def __init__(self, factor):
super().__init__(f"non-invertible denominator leaked factor {factor}")
self.factor = factor
class QuadRing:
__slots__ = ("n", "d", "zero", "one", "w")
def __init__(self, n, d=30):
self.n = n
self.d = d
self.zero = (0, 0)
self.one = (1 % n, 0)
self.w = (0, 1 % n)
def coerce(self, value):
if isinstance(value, tuple):
return (value[0] % self.n, value[1] % self.n)
return (value % self.n, 0)
def is_zero(self, x):
return x[0] % self.n == 0 and x[1] % self.n == 0
def is_one(self, x):
return x[0] % self.n == 1 and x[1] % self.n == 0
def add(self, x, y):
return ((x[0] + y[0]) % self.n, (x[1] + y[1]) % self.n)
def sub(self, x, y):
return ((x[0] - y[0]) % self.n, (x[1] - y[1]) % self.n)
def neg(self, x):
return ((-x[0]) % self.n, (-x[1]) % self.n)
def mul(self, x, y):
a0, a1 = x
b0, b1 = y
return (
(a0 * b0 - self.d * a1 * b1) % self.n,
(a0 * b1 + a1 * b0) % self.n,
)
def inv(self, x):
a0, a1 = x
norm = (a0 * a0 + self.d * a1 * a1) % self.n
g = gcd(norm, self.n)
if g != 1:
if g == self.n:
raise FactorFound(g)
raise FactorFound(g)
inv_norm = pow(norm, -1, self.n)
return ((a0 * inv_norm) % self.n, (-a1 * inv_norm) % self.n)
def div(self, x, y):
return self.mul(x, self.inv(y))
def square(self, x):
return self.mul(x, x)
def format(self, x):
a0, a1 = x
if a1 == 0:
return str(a0)
if a0 == 0:
return f"{a1}*w"
return f"{a0} + {a1}*w"
class Poly:
__slots__ = ("ctx", "coeffs")
def __init__(self, ctx, coeffs):
self.ctx = ctx
self.coeffs = [ctx.coerce(c) for c in coeffs]
self._trim()
def _trim(self):
while len(self.coeffs) > 1 and self.ctx.is_zero(self.coeffs[-1]):
self.coeffs.pop()
@classmethod
def zero(cls, ctx):
return cls(ctx, [ctx.zero])
@classmethod
def one(cls, ctx):
return cls(ctx, [ctx.one])
def copy(self):
return Poly(self.ctx, list(self.coeffs))
def degree(self):
if self.is_zero():
return -1
return len(self.coeffs) - 1
def lc(self):
return self.coeffs[-1]
def is_zero(self):
return len(self.coeffs) == 1 and self.ctx.is_zero(self.coeffs[0])
def is_one(self):
return len(self.coeffs) == 1 and self.ctx.is_one(self.coeffs[0])
def __eq__(self, other):
return self.coeffs == other.coeffs
def __add__(self, other):
size = max(len(self.coeffs), len(other.coeffs))
out = []
for i in range(size):
a = self.coeffs[i] if i < len(self.coeffs) else self.ctx.zero
b = other.coeffs[i] if i < len(other.coeffs) else self.ctx.zero
out.append(self.ctx.add(a, b))
return Poly(self.ctx, out)
def __sub__(self, other):
size = max(len(self.coeffs), len(other.coeffs))
out = []
for i in range(size):
a = self.coeffs[i] if i < len(self.coeffs) else self.ctx.zero
b = other.coeffs[i] if i < len(other.coeffs) else self.ctx.zero
out.append(self.ctx.sub(a, b))
return Poly(self.ctx, out)
def __neg__(self):
return Poly(self.ctx, [self.ctx.neg(c) for c in self.coeffs])
def scale(self, coeff):
coeff = self.ctx.coerce(coeff)
return Poly(self.ctx, [self.ctx.mul(c, coeff) for c in self.coeffs])
def shift(self, places):
if self.is_zero():
return Poly.zero(self.ctx)
return Poly(self.ctx, [self.ctx.zero] * places + self.coeffs)
def __mul__(self, other):
out = [self.ctx.zero] * (len(self.coeffs) + len(other.coeffs) - 1)
for i, a in enumerate(self.coeffs):
if self.ctx.is_zero(a):
continue
for j, b in enumerate(other.coeffs):
if self.ctx.is_zero(b):
continue
out[i + j] = self.ctx.add(out[i + j], self.ctx.mul(a, b))
return Poly(self.ctx, out)
def monic(self):
if self.is_zero():
return self.copy()
return self.scale(self.ctx.inv(self.lc()))
def divmod(self, other):
if other.is_zero():
raise ZeroDivisionError("polynomial division by zero")
if self.degree() < other.degree():
return Poly.zero(self.ctx), self.copy()
rem = self.copy()
q_coeffs = [self.ctx.zero] * (self.degree() - other.degree() + 1)
inv_lc = self.ctx.inv(other.lc())
while not rem.is_zero() and rem.degree() >= other.degree():
coeff = self.ctx.mul(rem.lc(), inv_lc)
shift = rem.degree() - other.degree()
q_coeffs[shift] = self.ctx.add(q_coeffs[shift], coeff)
rem = rem - other.scale(coeff).shift(shift)
return Poly(self.ctx, q_coeffs), rem
def __floordiv__(self, other):
q, r = self.divmod(other)
if not r.is_zero():
raise ArithmeticError("polynomial division was not exact")
return q
def __mod__(self, other):
_, r = self.divmod(other)
return r
def xgcd(self, other):
r0, r1 = self.copy(), other.copy()
s0, s1 = Poly.one(self.ctx), Poly.zero(self.ctx)
t0, t1 = Poly.zero(self.ctx), Poly.one(self.ctx)
while not r1.is_zero():
q, r2 = r0.divmod(r1)
r0, r1 = r1, r2
s0, s1 = s1, s0 - q * s1
t0, t1 = t1, t0 - q * t1
if r0.is_zero():
return r0, s0, t0
inv_lc = self.ctx.inv(r0.lc())
return r0.scale(inv_lc), s0.scale(inv_lc), t0.scale(inv_lc)
def __repr__(self):
terms = []
for i, coeff in enumerate(self.coeffs):
if self.ctx.is_zero(coeff):
continue
term = self.ctx.format(coeff)
if i == 0:
terms.append(term)
elif i == 1:
terms.append(f"({term})*x")
else:
terms.append(f"({term})*x^{i}")
return " + ".join(terms) if terms else "0"
class Divisor:
__slots__ = ("ctx", "f", "u", "v")
def __init__(self, ctx, f, u, v):
self.ctx = ctx
self.f = f
self.u = u
self.v = v
@classmethod
def zero(cls, ctx, f):
return cls(ctx, f, Poly.one(ctx), Poly.zero(ctx))
def is_zero(self):
return self.u.is_one() and self.v.is_zero()
def __neg__(self):
if self.is_zero():
return self
return Divisor(self.ctx, self.f, self.u, -self.v)
def __eq__(self, other):
return self.u == other.u and self.v == other.v
def _compose(self, other):
a1, b1 = self.u, self.v
a2, b2 = other.u, other.v
if a1 == a2 and b1 == b2:
d, _, h3 = a1.xgcd(b1.scale(2))
a = (a1 // d)
a = a * a
b = (b1 + h3 * ((self.f - (b1 * b1)) // d)) % a
else:
d0, _, h2 = a1.xgcd(a2)
if d0.is_one():
a = a1 * a2
b = (b2 + h2 * a2 * (b1 - b2)) % a
else:
d, l, h3 = d0.xgcd(b1 + b2)
a = (a1 * a2) // (d * d)
b = (
b2
+ l * h2 * (b1 - b2) * (a2 // d)
+ h3 * ((self.f - (b2 * b2)) // d)
) % a
return a.monic(), b
def _reduce(self, a, b):
a2 = ((self.f - (b * b)) // a).monic()
b2 = (-b) % a2
if a2.degree() > 2:
return self._reduce(a2, b2)
return a2, b2
def add(self, other):
if self.is_zero():
return other
if other.is_zero():
return self
a, b = self._compose(other)
if a.degree() > 2:
a, b = self._reduce(a, b)
return Divisor(self.ctx, self.f, a, b)
def scalar_mul(self, k, progress=False):
out = Divisor.zero(self.ctx, self.f)
addend = self
total_bits = k.bit_length()
for bit_index in range(total_bits):
if progress and bit_index and bit_index % 256 == 0:
print(f"[+] processed {bit_index}/{total_bits} bits")
if (k >> bit_index) & 1:
out = out.add(addend)
if bit_index + 1 != total_bits:
addend = addend.add(addend)
return out
def __repr__(self):
return f"({self.u}, {self.v})"
def is_probable_prime(n):
if n < 2:
return False
small_primes = (
2,
3,
5,
7,
11,
13,
17,
19,
23,
29,
31,
37,
41,
43,
47,
53,
59,
61,
67,
71,
)
for p in small_primes:
if n == p:
return True
if n % p == 0:
return False
d = n - 1
s = 0
while d % 2 == 0:
s += 1
d //= 2
for a in small_primes:
if a >= n:
continue
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = pow(x, 2, n)
if x == n - 1:
break
else:
return False
return True
def tonelli_shanks(a, p):
a %= p
if a == 0:
return 0
if p == 2:
return a
if pow(a, (p - 1) // 2, p) != 1:
raise ValueError("not a square")
if p % 4 == 3:
return pow(a, (p + 1) // 4, p)
q = p - 1
s = 0
while q % 2 == 0:
s += 1
q //= 2
z = 2
while pow(z, (p - 1) // 2, p) != p - 1:
z += 1
m = s
c = pow(z, q, p)
t = pow(a, q, p)
r = pow(a, (q + 1) // 2, p)
while t != 1:
i = 1
t2i = pow(t, 2, p)
while t2i != 1:
t2i = pow(t2i, 2, p)
i += 1
if i == m:
raise ValueError("tonelli failure")
b = pow(c, 1 << (m - i - 1), p)
m = i
c = pow(b, 2, p)
t = (t * c) % p
r = (r * b) % p
return r
def cornacchia_x2_plus_2y2(p):
z = tonelli_shanks(-2, p)
for root in (z, p - z):
a, b = p, root
while b * b > p:
a, b = b, a % b
x = abs(b)
t = p - x * x
if t >= 0 and t % 2 == 0:
y2 = t // 2
y = isqrt(y2)
if y * y == y2:
return x, y
raise ValueError("no x^2 + 2y^2 representation found")
def integer_cuberoot(n):
if n < 0:
raise ValueError("cube root only supports nonnegative integers")
if n < 8:
return 0 if n == 0 else 1
x = 1 << ((n.bit_length() + 2) // 3)
while True:
y = (2 * x + n // (x * x)) // 3
if y >= x:
while (x + 1) ** 3 <= n:
x += 1
while x ** 3 > n:
x -= 1
return x
x = y
def build_curve_data(n):
ctx = QuadRing(n)
f = Poly(ctx, [0, -1, 0, 0, 0, 1])
u = Poly(ctx, [2, 3, 1])
v = Poly(ctx, [ctx.w, ctx.w])
return ctx, f, Divisor(ctx, f, u, v)
def try_complete_from_p(n, candidate):
if candidate <= 1 or n % candidate != 0:
return None
if candidate % 8 != 3 or not is_probable_prime(candidate):
return None
try:
x, y = cornacchia_x2_plus_2y2(candidate)
except ValueError:
return None
q = (((candidate - 1) * (candidate - 1)) + 4 * x * x) // 8
if q <= 1 or n % (candidate * q) != 0:
return None
r = n // (candidate * q)
if not is_probable_prime(q) or not is_probable_prime(r):
return None
return candidate, q, r, x, y
def try_complete_from_q(n, candidate_q):
if candidate_q <= 1 or n % candidate_q != 0 or not is_probable_prime(candidate_q):
return None
root = isqrt(8 * candidate_q)
for delta in range(-4, 5):
candidate_p = root + delta
result = try_complete_from_p(n, candidate_p)
if result is not None and result[1] == candidate_q:
return result
return None
def try_complete_from_pq_product(n, candidate_pq):
if candidate_pq <= 1 or n % candidate_pq != 0:
return None
root = integer_cuberoot(8 * candidate_pq)
for delta in range(-8, 9):
candidate_p = root + delta
result = try_complete_from_p(n, candidate_p)
if result is not None and result[0] * result[1] == candidate_pq:
return result
return None
def resolve_factor_leak(n, leaked_factor):
candidates = [leaked_factor, n // leaked_factor]
for candidate in candidates:
result = try_complete_from_p(n, candidate)
if result is not None:
return result
for candidate in candidates:
result = try_complete_from_q(n, candidate)
if result is not None:
return result
for candidate in candidates:
result = try_complete_from_pq_product(n, candidate)
if result is not None:
return result
return None
def factor_special_rsa(n, progress=False):
_, _, base = build_curve_data(n)
scalar = 2 * n
try:
base.scalar_mul(scalar, progress=progress)
except FactorFound as leak:
g = leak.factor
if progress:
print(f"[+] leaked gcd = {g}")
result = resolve_factor_leak(n, g)
if result is not None:
return result
raise RuntimeError("got a nontrivial gcd, but could not turn it into p")
raise RuntimeError("scalar multiplication finished without leaking a factor")
def verify_known_instance():
x, y = cornacchia_x2_plus_2y2(KNOWN_P)
q = (((KNOWN_P - 1) * (KNOWN_P - 1)) + 4 * x * x) // 8
print("p*q*r == n ?", KNOWN_P * KNOWN_Q * KNOWN_R == KNOWN_N)
print("p = x^2 + 2y^2 ?", KNOWN_P == x * x + 2 * y * y)
print("recovered q == known q ?", q == KNOWN_Q)
print("n bits =", KNOWN_N.bit_length())
print("p bits =", KNOWN_P.bit_length())
print("q bits =", KNOWN_Q.bit_length())
print("r bits =", KNOWN_R.bit_length())
def run_small_demo():
p, q, r = 11, 17, 23
n = p * q * r
print("[*] demo modulus =", n)
got_p, got_q, got_r, x, y = factor_special_rsa(n, progress=False)
print("[+] recovered p =", got_p)
print("[+] recovered q =", got_q)
print("[+] recovered r =", got_r)
print("[+] x, y =", x, y)
def parse_args():
parser = argparse.ArgumentParser(description="Attack the special RSA in test.py")
parser.add_argument("--n", type=lambda s: int(s, 0), default=KNOWN_N)
parser.add_argument("--check", action="store_true", help="verify the known instance")
parser.add_argument("--demo", action="store_true", help="run a tiny toy example")
parser.add_argument("--progress", action="store_true", help="print scalar-mul progress")
return parser.parse_args()
def main():
args = parse_args()
if args.check:
verify_known_instance()
return
if args.demo:
run_small_demo()
return
start = time.time()
p, q, r, x, y = factor_special_rsa(args.n, progress=args.progress)
elapsed = time.time() - start
print("p =", p)
print("q =", q)
print("r =", r)
print("x =", x)
print("y =", y)
print("elapsed =", round(elapsed, 3), "seconds")
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment