## solution

If you can compute $s_1$ from $s_1^e, {(as_1+b)}^e, {(a^2s_1+ab+b)}^e \cdots$, you can recover $s_2^e, {(as_2+b)}^e, {(a^2s_2+ab+b)}^e \cdots$. Now you can compute $s_2$, you also can recover $s_3^e, {(as_3+b)}^e, {(a^2s_3+ab+b)}^e \cdots$. By repeating this, you can compute all of the paramators(and can encrypt any desired plaintext).

So, let's focus on compute $s_1$ from $s_1^e, {(as_1+b)}^e, {(a^2s_1+ab+b)}^e \cdots$. This can be computed by using polynomial gcd, as compute $\rm{gcd}(s_1^e-\rm{debugout}_1, {(as_1+b)}^e-\rm{debugout}_2)$. The fastest way to do this (as far as I know) is using Half-GCD algorithm. But considering that this challenge has timelimit of 150 seconds and you should do this 30 times in a row, doing Half-GCD of polynomials that have degrees of 65535 looks quite infeasible.

But you can speed up this process. To understand, we need to talk about the inside of the half-gcd algorithm. This part is translation of [NyaanNyaan's article](https://nyaannyaan.github.io/library/fps/polynomial-gcd.hpp.html).

> There is a polynomial $a(x),b(x)$. For any polynomial $h$ such that $h\mid a,h\mid b$ holds $g\mid h$, let $\gcd(a,b)$ be the GCD of the polynomial.
> The $\gcd(a,b)$ can be obtained by Euclidian algorithm in the same way as for ordinary gcd, but since the order is reduced by at least one per division, the computational complexity is about $\mathrm{O}(N^2 \log N)$ when calculated naively. Therefore, an algorithm called the Half-GCD method is used to speed up the calculation.
>
> A single operation in the usual Euclidean reciprocal division scheme, expressed as a matrix, is shown below
> $$\left(
\begin{array}{cc}
0 & 1 \\
1 & -\mathrm{quo}(a,b)
\end{array}
\right)
\left(
\begin{array}{c}
a \\
b
\end{array}
\right)=
\left(
\begin{array}{c}
b \\
a - b\cdot \mathrm{quo}(a,b)
\end{array}
\right)$$
> Similarly, Euclidian algorithm is expressed as follows.
> $$
M_i M_{i-1}\cdots M_2M_1\left(
\begin{array}{c}
a \\
b
\end{array}
\right)=
\left(
\begin{array}{c}
\gcd(a,b) \\
0
\end{array}
\right)$$
> Therefore, if the matrices $M_1,M_2,M_3,\ldots,M_i$ can be obtained, then $\gcd(a,b)$ also can be obtained.
> The Half-GCD algorithm is an algorithm that takes advantage of the following facts to find approximately the front half of a matrix column $M$ with less computation.
> - Let $a = a_0 + a_1x^k, b=b_0+b_1x^k$. ($k \lt \deg(a) \le \deg(b))$ Here, the matrices $L_1,L_2,\dots$ and $M_1,M_2\ldots$ obtained when applying Euclidian algorithm to $(a_1,b_1)$ are partly identical(upto $\lfloor\frac{\deg(b)-k}{2}\rfloor$), starting from the beginning.

The last fact can be useful for our situation, since $s_1^e-\rm{debugout}_1$ and ${(as_1+b)}^e-\rm{debugout}_2$ always has the same coefficients upto $x^1$. This means we can pre-compute half of matricies.

By using this matricies, we can reduce the degree of polynomials to the half, which is 32768. Dispite that we reduced the degree, Applying Half-GCD algorithm to polynomial that has 32768 degree has also takes too much time.

But, what if we take in to the account the other polynomials? Long story short, if there is a $n$ polynomials that can be notated as $p_i = p_{i,0} + p_{i,1}x^k$, there is a matrix that reduce the sequence of $p_i$ into the degree of $\lfloor\frac{\max(\deg(p))-k}{n}\rfloor$ no matter what the value of $p_{i,0}$ is.

Below is the sagemath code that compute the matrix. If you are familiar with the implementation of Half-GCD, you will notice that the broad structure is similar.

G = Zmod(n)
F.<x> = PolynomialRing(G)

def reduce(polys: List[Polynomial], mat: Matrix):
  n = len(polys)

  start_deg = polys[0].degree()
  for i in range(n):
    if polys[i].degree() < start_deg:
      polys[i] += polys[0]
  assert all([poly.degree() == start_deg for poly in polys])

  def reduce_poly(i, j):
    q = polys[i] // polys[j]
    assert q.degree() <= 1
    mat[i] -= q * mat[j]
    polys[i] -= q * polys[j]

  for divisor_ind in range(n-1):
    for dividend_ind in range(n):
      if divisor_ind == dividend_ind: continue
      reduce_poly(dividend_ind, divisor_ind)

  reduce_poly(-2, -1)
  mat[-2] += mat[-1]
  polys[-2] += polys[-1]

  end_deg = polys[0].degree()
  assert start_deg - end_deg == n - 1
  assert all([poly.degree() == end_deg for poly in polys])

def bulk_gcd(polys: List):
  n = len(polys)
  
  if polys[0].degree() <= n + 1:
    return matrix.identity(F, n)

  m  = polys[0].degree() // 2

  M_1 = bulk_gcd([poly // x^m for poly in polys])
  V = M_1 * matrix(polys).T
  new_polys = list(V.T[0])

  if new_polys[0].degree() <= n + 1:
    return M_1
  
  reduce(new_polys, M_1)

  if new_polys[0].degree() <= n + 1:
    return M_1

  reduced_deg = polys[0].degree() - new_polys[0].degree()
  gained_pollution = ceil(reduced_deg / (n - 1))

  M_2 = bulk_gcd([poly // x^gained_pollution for poly in new_polys])
  return M_2 * M_1

Using this `bulk_gcd` function, we can compute the matrix for the $s_1^e, {(as_1+b)}^e, {(a^2s_1+ab+b)}^e \cdots$ as follows. `GATHER_COUNT` is the number of the polynomials. Here, we will assume that `GATHER_COUNT` is 20. In this case, this computation takes 30 minutes in my environment(ThinkPad X1 Carbon Gen9, Intel Core i7-1165G7).

# GATHER_COUNT = 10 # precalc:  8 mins
GATHER_COUNT = 20   # precalc: 30 mins
# GATHER_COUNT = 40 # precalc: 98 mins

prev_s = x
fs_prime = [prev_s^e]
while len(fs_prime) < GATHER_COUNT:
  prev_s = a * prev_s + b
  fs_prime.append(prev_s^e)

M = bulk_gcd(fs_prime)

Now, we can reduce the degree of the $s_1^e, {(as_1+b)}^e, {(a^2s_1+ab+b)}^e \cdots$ into 3277. GCD of these polynomials can be easily computed by Half-GCD algorithm. Moreover, we only need two of these polynomials to apply Half-GCD algorithms. Below is the code to compute the s using this method.

# rs: debut outputs
fs = [f - r for f, r in zip(fs_prime, rs)]
V = matrix(fs).T
V2 = M[:2] * V

u, v = V2[0, 0], V2[1, 0]
res = pgcd(u, v).monic()
s = -res.constant_coefficient()

But, now we have the problem. Computing `M[:2] * V` is too slow. This can be resolved using the fact that u and v has small degrees compare to fs. Since the larger element of the fs is not affecting to the essential part of u and v, we can trim `fs`.

fs = [f - 1337 for f in fs_prime]
M2 = M[:2]
V2 = M2 * matrix(fs).T
required_degree = max(V2[0,0].degree(), V2[1,0].degree())+1

fs_reduced = [f[:required_degree+1] for f in fs_prime]


fs = [f - r for f, r in zip(fs_reduced, rs)]
V = matrix(fs).T
V2 = M2 * V
u, v = V2[0, 0][:required_degree+1], V2[1, 0][:required_degree+1]
res = pgcd(u, v).monic()
s = -res.constant_coefficient()

Now, we can officialy compute the gcd of $s_1^e-\rm{debugout}_1, {(as_1+b)}^e-\rm{debugout}_2, {(a^2s_1+ab+b)}^e-\rm{debugout}_3 \cdots$.

The time constraints are a little tighter, as we definitely did not want to be solved with a simple Half-GCD. If there were teams that were able to get to the solution but could not solve it, I sincerely apologise.

Below is the full solution script. `pgcd.sage` is taken from [here](https://scrapbox.io/crypto-writeup-public/half-GCD).

import pickle
from typing import List
from ptrlib import remote
from tqdm import tqdm
from sage.all import *

load("pgcd.sage")

LOCAL = False

def get_conn():
    return remote("nc crypto.2023.zer0pts.com 7002")

conn = get_conn()
a = int(conn.recvline().strip().split()[2])
b = int(conn.recvline().strip().split()[2])
e = int(conn.recvline().strip().split()[2])
n = int(conn.recvline().strip().split()[2])
conn.close()

print(f'[+] {a=}')
print(f'[+] {b=}')
print(f'[+] {e=}')
print(f'[+] {n=}')

G = Zmod(n)
F.<x> = PolynomialRing(G)

def reduce(polys: List[Polynomial], mat: Matrix, debug=False):
  n = len(polys)

  start_deg = polys[0].degree()
  for i in range(n):
    if polys[i].degree() < start_deg:
      polys[i] += polys[0]
  assert all([poly.degree() == start_deg for poly in polys])

  def reduce_poly(i, j):
    q = polys[i] // polys[j]
    assert q.degree() <= 1
    mat[i] -= q * mat[j]
    polys[i] -= q * polys[j]

  for divisor_ind in range(n-1):
    PROGRESS.update()
    for dividend_ind in range(n):
      if divisor_ind == dividend_ind: continue
      reduce_poly(dividend_ind, divisor_ind)

  reduce_poly(-2, -1)
  mat[-2] += mat[-1]
  polys[-2] += polys[-1]

  end_deg = polys[0].degree()
  if debug: print(f'{end_deg=} {[poly.degree() for poly in polys]=}')
  assert start_deg - end_deg == n - 1
  assert all([poly.degree() == end_deg for poly in polys])

  return mat

def bulk_gcd(polys: List):
  n = len(polys)
  
  if polys[0].degree() <= n + 1:
    return matrix.identity(F, n)

  m  = polys[0].degree() // 2

  M_1 = bulk_gcd([poly // x^m for poly in polys])
  V = M_1 * matrix(polys).T
  new_polys = list(V.T[0])

  if new_polys[0].degree() <= n + 1:
    return M_1
  
  reduce(new_polys, M_1)

  if new_polys[0].degree() <= n + 1:
    return M_1

  reduced_deg = polys[0].degree() - new_polys[0].degree()
  gained_pollution = ceil(reduced_deg / (n - 1))

  M_2 = bulk_gcd([poly // x^gained_pollution for poly in new_polys])
  return M_2 * M_1

from time import time

ROUND = 30

# GATHER_COUNT = 10 # precalc:  8 mins, gcd: 7.3 secs/round (recovering: 0.6 secs, half-gcd: 6.8 secs)
GATHER_COUNT = 20   # precalc: 30 mins, gcd: 3.5 secs/round (recovering: 0.6 secs, half-gcd: 2.8 secs)
# GATHER_COUNT = 40 # precalc: 98 mins, gcd: 2.1 secs/round (recovering: 0.6 secs, half-gcd: 1.4 secs)

prev_s = x
fs_prime = [prev_s^e]
while len(fs_prime) < GATHER_COUNT:
  prev_s = a * prev_s + b
  fs_prime.append(prev_s^e)

PICKLE_FILE = f"data_{e}_{GATHER_COUNT}.pickle"

if os.path.exists(PICKLE_FILE):
  data = pickle.load(open(PICKLE_FILE, "rb"))
  print("[+] data loaded from cache")
  assert data["a"] == a
  assert data["b"] == b
  assert data["e"] == e
  assert data["n"] == n
else:
  print("[+] calculating matrix...")
  start = time()
  PROGRESS = tqdm(total=int(e * (GATHER_COUNT - 1) // GATHER_COUNT), smoothing=0)
  M = bulk_gcd(fs_prime)
  PROGRESS.close()
  end = time()
  print(f"[+] calculation takes: {end - start:.3f}")
  data = {}
  data["a"] = a
  data["b"] = b
  data["e"] = e
  data["n"] = n
  data["M"] = [[list(poly) for poly in row] for row in M]
  pickle.dump(data, open(PICKLE_FILE, "wb"))
  print(f"[+] saved to {PICKLE_FILE}")

M = matrix(F, [[F(poly) for poly in row] for row in data["M"]])

print("[+] checking required degrees")
fs = [f - 1337 for f in fs_prime]
M2 = M[:2]
V2 = M2 * matrix(fs).T
required_degree = max(V2[0,0].degree(), V2[1,0].degree())+1

print(f'[+] {required_degree=}')

fs_reduced = [f[:required_degree+1] for f in fs_prime]

recovering_ts = []
half_gcd_ts = []

def solve_round(rs):
  fs = [f - r for f, r in zip(fs_reduced, rs)]
  start = time()
  V = matrix(fs).T
  V2 = M2 * V
  end = time()
  print(f' | recovering takes: {end - start:.3f} secs'); recovering_ts.append(end - start)
  u, v = V2[0, 0][:required_degree+1], V2[1, 0][:required_degree+1]
  print(f' | {u.degree()=} {v.degree()=}')
  start = time()
  res = pgcd(u, v).monic()
  end = time()
  print(f' \\ half-gcd takes: {end - start:.3f} secs'); half_gcd_ts.append(end - start)
  assert res.degree() == 1
  s = -res.constant_coefficient()
  assert pow(s, e, n) == rs[0]
  return s

conn = get_conn()
start = time()
def elapsed(): return round(time() - start, 3)

for i in range(GATHER_COUNT):
  conn.sendline(b'00')
encrypted_rs_list = []
for i in range(GATHER_COUNT):
  conn.recvuntil(b'> ')
  leaks = [0]
  for j in range(ROUND):
    leaks.append(int(conn.recvlineafter(": ")))
  encrypted_rs = [a ^^ b for a, b in zip(leaks[:-1], leaks[1:])]
  assert(len(encrypted_rs) == ROUND)
  encrypted_rs_list.append(encrypted_rs)


result_key = 0
prev_keys = [0] * GATHER_COUNT

for i in range(ROUND):
  print(f'[+] {elapsed()=} round {i+1} / {ROUND}...')
  rs = [encrypted_rs[i] ^^ prev_keys[j] for j, encrypted_rs in enumerate(encrypted_rs_list)]
  s = int(solve_round(rs))
  prev_keys = []
  while len(prev_keys) < GATHER_COUNT:
    assert pow(s, e, n) == rs[len(prev_keys)]
    prev_keys.append(s)
    s = (a * s + b) % n
  result_key ^^= pow(s, e, n)
  result_key ^^= s

payload = int(result_key ^^ int.from_bytes(b"Give me the flag!", "big")).to_bytes(128, "big")

print(f'[+] {elapsed()=} send payload: {payload.hex()=}')
conn.sendlineafter("> ", payload.hex())

print(f'[+] {elapsed()=} all done! {elapsed() / ROUND:.3f} secs / round (recovering: {sum(recovering_ts) / ROUND:.2f} secs, half-gcd: {sum(half_gcd_ts) / ROUND:.2f} secs)')
print(conn.recvlineafter("The flag is: ", timeout=5).decode())

conn.close()

<!--
$$\left(
\begin{array}{cc}
1 & 0 & \cdots & -\rm{quo}(p_0, p_i) & \cdots & 0 & \cdots & 0 \\
0 & 1 & \cdots & -\rm{quo}(p_1, p_i) & \cdots & 0 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 1 & \cdots & 0 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & -\rm{quo}(p_j, p_i) & \cdots & 1 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & -\rm{quo}(p_n, p_i) & \cdots & 0 & \cdots & 1
\end{array}
\right)
\left(
\begin{array}{c}
p_1 \\
p_2 \\
\vdots \\
p_i \\
\vdots \\
p_j \\
\vdots \\
p_n
\end{array}
\right)=
\left(
\begin{array}{c}
p_1 - p_i \cdot \rm{quo}(p_1, p_i) \\
p_2 - p_i \cdot \rm{quo}(p_2, p_i) \\
\vdots \\
p_i \\
\vdots \\
p_j - p_i \cdot \rm{quo}(p_j, p_i) \\
\vdots \\
p_n - p_i \cdot \rm{quo}(p_n, p_i)
\end{array}
\right)$$

o: known
x: unknown

X: oooooxx <- deg(x)==6, lower two coefficients polynomial X has 

A: aooooox
B: booooox
C: cooooox
D: dooooox

↓

A: aeoooox
B:  foooox (-=b/a*A)
C:  goooox (-=c/a*A)
D:  hoooox (-=d/a*A)

↓

A:   oooxx (-=a/f*B*x-e/f*B)
B:  ooooox
C:   oooox (-=g/f*B)
D:   ooooo (-=h/f*B)

↓

A:    ooxx
B:    ooxx
C:   oooox
D:    ooox

↓

A:    ooxx
B:    ooxx
C:    ooxx
D:    ooox

If we implement the

The time constraints are a little tighter, as we definitely did not want to be solved with a simple half-gcd. If there were teams that were able to get to the solution but could not solve it, we sincerely apologise.
-->