Skip to content

Instantly share code, notes, and snippets.

@kusano
Created August 12, 2022 06:56
Show Gist options
  • Save kusano/fa8f3b8d71b3a146fa5bc940e0e944b1 to your computer and use it in GitHub Desktop.
Save kusano/fa8f3b8d71b3a146fa5bc940e0e944b1 to your computer and use it in GitHub Desktop.
from fractions import Fraction
import math
import random
# 内積
def prod(u, v):
return sum(x*y for x, y in zip(u, v))
# Gram-Schmidtの直交化
# GSOベクトルのノルムの2乗 [|b0|^2, |b1|^2, ...] とGSO係数 mu を返す。
def GSO(M):
n = len(M)
mu = [[0]*n for _ in range(n)]
B = [m[:] for m in M]
for i in range(n):
mu[i] = [0]*i
for j in range(i):
mu[i][j] = Fraction(prod(M[i], B[j]), prod(B[j], B[j]))
for k in range(n):
B[i][k] -= mu[i][j]*B[j][k]
return [prod(b, b) for b in B], mu
# LLL基底簡約
def LLL(M, delta=Fraction(3, 4)):
n = len(M)
M = [m[:] for m in M]
B, mu = GSO(M)
k = 1
while k<n:
#print(k, [int(b).bit_length() for b in B])
# サイズ簡約
for i in range(k)[::-1]:
if abs(mu[k][i])>=Fraction(1, 2):
q = round(mu[k][i])
for l in range(n):
M[k][l] -= q*M[i][l]
# mu の更新
for l in range(i):
mu[k][l] -= q*mu[i][l]
mu[k][i] -= q
# Lovasz条件
if B[k]>=(delta-mu[k][k-1]**2)*B[k-1]:
k += 1
else:
# M[k] と M[k-1] の交換
M[k], M[k-1] = M[k-1], M[k]
# B の更新
b_old = B[k-1]
B[k-1] = B[k] + mu[k][k-1]**2*b_old
B[k] *= b_old/B[k-1]
# mu の更新
mu_old = mu[k][k-1]
mu[k][k-1] *= b_old/B[k-1]
for i in range(k-1):
t = mu[k][i]
mu[k][i] = mu[k-1][i]
mu[k-1][i] = t
for i in range(k+1, n):
t = mu[i][k]
mu[i][k] = mu[i][k-1] - mu_old*t
mu[i][k-1] = t + mu[k][k-1]*mu[i][k]
k = max(1, k-1)
return M
# p(x) = Σ[0<=i<=k]A[i]*x^i = 0 mod N について、
# |x0|<(1/2)*N^(1/k)-e を満たす解 x0 を返す。
def coppersmith(A, N, e=Fraction(1, 8)):
# パラメタの決定
k = len(A)-1
h = math.ceil(max(7/k, (k+e*k-1)/(e*k**2)))
X = int(N**(1/k-e)/2)
print(f"{k=}")
print(f"{h=}")
print(f"{X=}")
# 基底行列の構築
M = [[0]*(2*h*k-k) for _ in range(2*h*k-k)]
# 左上
for i in range(h*k):
M[i][i] = X**(h*k-i-1)
# 右上
A2 = [1] # (p(x))^j の係数
for j in range(h-1):
T = [0]*(len(A)+len(A2)-1)
for i0 in range(len(A2)):
for i1 in range(len(A)):
T[i0+i1] += A2[i0]*A[i1]
A2 = T
for l in range(k):
for i in range(len(A2)):
M[i+l][h*k+j*k+l] = A2[i]
# 右下
for j in range(h-1):
for l in range(k):
M[h*k+j*k+l][h*k+j*k+l] = N**(j+1)
# 右上を零行列、右下を単位行列にする。
M = M[h*k:]+M[:h*k]
for i in range(2*h*k-k):
for j in range(h*k, 2*h*k-k):
if j!=i:
t = M[i][j]
for l in range(2*h*k-k):
M[i][l] -= t*M[j][l]
# 左上の取り出し
M = [m[:h*k] for m in M[:h*k]]
# LLL基底簡約
B = LLL(M)
n = h*k
# B[:-1]*F = 0 となる F を求める。
B = B[:-1]
for i in range(n-1):
for j in range(n):
B[i][j] = Fraction(B[i][j]//X**(n-j-1))
for i in range(n-1):
for j in range(i, n-1):
if M[j][i]!=0:
B[i], B[j] = B[j], B[i]
break
assert B[i][i]!=0
t = B[i][i]
for j in range(n):
B[i][j] /= t
for j in range(n-1):
if j!=i:
t = B[j][i]
for l in range(n):
B[j][l] -= t*B[i][l]
g = 1
for i in range(n-1):
d = B[i][n-1].denominator
g = g*d//math.gcd(g, d)
F = [int(-B[i][n-1]*g) for i in range(n-1)]+[g]
print(f"{F=}")
# ΣF[i]*x^i = 0 を解く。
p = lambda x: F[0]+sum(F[i]*x**i for i in range(1, n))
while True:
x1 = random.randint(-X, X+1)
x2 = random.randint(-X, X+1)
if p(x1)*p(x2)<0:
while abs(x1-x2)>1:
m = (x1+x2)//2
if p(m)*p(x1)>=0:
x1 = m
else:
x2 = m
if p(x1)==0:
# 元の方程式の解となっていることも確認する。
if (A[0]+sum(A[i]*x1**i for i in range(1, k+1)))%N==0:
x0 = x1
break
return x0
N = 135896252889455179306044561856073657176589282065060485982049884024159935440451648268809083278093186378968141115107175023955352285741931510695513569631865587980780116850562946289647504541616229544460216616423916073902516352634221984329472372337896406341980543263900139265196004232936560031956301070039887638201
e = 3
c = 104064417885709974759351196219773048348376455313832829209947346446874683031677099263660660465743321265675101447205286367005536581370429707934588639596542365006675867115121072960551019978628425500909321182750343178517807287529652070367005206132362187414189170498876416717871986142239844250496022317643124516521
m0 = """
+>-<+>-<+>-<+>-<+>-<+>-<+>-<+>-<+>-<+>-<+
| FLAG{""" + "\0"*25 + """} |
+>-<+>-<+>-<+>-<+>-<+>-<+>-<+>-<+>-<+>-<+
"""
m0 = int.from_bytes(m0.encode(), "big")
# A[i]: (m0 + x * 256^49)^e - c の x^i の係数
A = [
m0**3 - c,
3 * m0**2 * 256**49,
3 * m0 * (256**49)**2,
(256**49)**3,
]
# A[3]を1にする。
for i in range(4):
A[i] = A[i] * pow((256**49)**3, -1, N) % N
flag = coppersmith(A, N)
print(flag.to_bytes(25, "big").decode())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment