https://sanya.sweetduet.info/ctfcoppersmith/

solve.py

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())