# CIGISICGICGICG - SECCON CTF 2023

## problem.py

import os
from functools import reduce
from secrets import randbelow


flag = os.getenvb(b"FLAG", b"FAKEFLAG{THIS_IS_FAKE}")
p1 = 21267647932558653966460912964485513283
a1 = 6701852062049119913950006634400761786
b1 = 19775891958934432784881327048059215186
p2 = 21267647932558653966460912964485513289
a2 = 10720524649888207044145162345477779939
b2 = 19322437691046737175347391539401674191
p3 = 21267647932558653966460912964485513327
a3 = 8837701396379888544152794707609074012
b3 = 10502852884703606118029748810384117800


def prod(x: list[int]) -> int:
    return reduce(lambda a, b: a * b, x, 1)


def xor(x: bytes, y: bytes) -> bytes:
    return bytes([xi ^ yi for xi, yi in zip(x, y)])


class ICG:
    def __init__(self, p: int, a: int, b: int) -> None:
        self.p = p
        self.a = a
        self.b = b
        self.x = randbelow(self.p)

    def _next(self) -> int:
        if self.x == 0:
            self.x = self.b
            return self.x
        else:
            self.x = (self.a * pow(self.x, -1, self.p) + self.b) % self.p
            return self.x


class CIG:
    L = 256

    def __init__(self, icgs: list[ICG]) -> None:
        self.icgs = icgs
        self.T = prod([icg.p for icg in self.icgs])
        self.Ts = [self.T // icg.p for icg in self.icgs]

    def _next(self) -> int:
        ret = 0
        for icg, t in zip(self.icgs, self.Ts):
            ret += icg._next() * t
            ret %= self.T
        return ret % 2**self.L

    def randbytes(self, n: int) -> bytes:
        ret = b""
        block_size = self.L // 8
        while n > 0:
            ret += self._next().to_bytes(block_size, "big")[: min(n, block_size)]
            n -= block_size
        return ret


if __name__ == "__main__":
    random = CIG([ICG(p1, a1, b1), ICG(p2, a2, b2), ICG(p3, a3, b3)])
    enc_flag = xor(flag, random.randbytes(len(flag)))
    leaked = random.randbytes(300)
    print(f"{enc_flag = }")
    print(f"{leaked = }")

## output.txt

enc_flag = b'\xd5 \xc3b\xa3\xa1\xd6\xe5Sv\xe7%n\xd6\xd6UcQNYU\x1arR\xdes\xb4\x12\xc9\xed\x1a\xc6^=\xe1\xe3p@\xe65\x19\x18S\x80\xa4TE\x7f\x92\x07&"\xdf\xc9\xe1\xbd@QL\xcf\x90\x98\xd9C$\xcb\xb4U'
leaked = b"^\xed>\x03\xad\x8c\x1d\xa1\xe29\x83\x92\xbdm\xefL\xe5\xe5\xab\xc9\xffZ\xbd8\x95\x97\xa3i/k\xb1\x8dSD\x1e\x92;\x87\xa7\x16\xdc\x98\x15\x1ba\xc3fQ\xa9\t\xe8ak.0\xe3\x93\xba\x82\xb2%\xc2\x88]u\xeb\xfctKw\xe5\xcc\xd2\xce\xa7\x8c\xd6T\xe3\xfa$\xec\xca\xcc\x1a\x08\xbd3\xdd!D\xc8\xa7}\xeb\xd2=\xfb\x96\xeek\xdef>\xedm\t\x12\xe6\xeeO\xc5\xbe\xcev\x9aB\x90\x84\x981j'\xb18\xbb\x08\x93\xbd\xf9\xb1>/\x81\x83]\x93C\x84D\x9b5\xd0l\xcfQa\xe3\x1ev !\xd6W\xbc\x9b\xccV\xd65\x84\t\xd2\xdd\xde\xffs\xcc\x80\x16\x9cg\xcf\xa4&l\x8f\x82J\x16\xc7qNN\x90\x89\xef\xa6\xb8\x8c\xcb\xf8q\x0f.)\xa7 \x8b\x14\x83\xca-\x7fvP\x1a\x08\xb6^\x18\xd5\x9b\x01\xfa[\xdf3J\xc0\x85\x02\xe3\x16\\\x93\x17B\xd6\x8e2\xabia\xf1hT+][\x19c<\x06\xea%m\xc0\x01\xc6'\x95t\xf3\xf4\xd7\xe1f\xcd\x8f\xb0\xa3\\\xcfv\xa8\xfb\xb6\x03\xc4R\xe0\x10\xbb\xcb>\x0e\x94H8\xbe\x0c\xf6\x9c\xbf\xa1^\x178\t1\xda\xd4\xc3cm\x84}\x9d\x84"

## CIG to singular ICG

We can turn the CIG into a single combined ICG by multiplying the output by the CRT inverse components, since the CIG construction is very similar to CRT, example here: https://hackmd.io/@Utaha1228/By1GclUAh.

T = p1 * p2 * p3
M = crt([pow(p2 * p3, -1, p1), pow(p1 * p3, -1, p2), pow(p1 * p2, -1, p3)], [p1, p2, p3])
A = crt([a1, a2, a3], [p1, p2, p3])
B = crt([b1, b2, b3], [p1, p2, p3])

Using the above values we can define a new ICG $u_i = M*\text{CIG}_i$:
$$u_i = \frac{A}{u_{i - 1}} + B \bmod T\\
M*\text{CIG}_i = \frac{A}{M*\text{CIG}_{i - 1}} + B \bmod T$$

## Multivariate coppermsith given truncated output
According to the paper: https://www.iacr.org/archive/pkc2012/72930609/72930609.pdf, it should be possible, given enough bits of $\text{CIG}_i$, to solve for the original value of $\text{CIG}_i$ using multivariate Coppersmith.

We're given a leaked output which corresponds to the lower 256 bits of $\text{CIG}_i$, so we can define:

$$\text{CIG}_i = r_i + 2^{256}*s_i$$

The $r_i$ is known, and the $s_i$ is unknown. All $p_i < 2^{125}$, and $T < 2^{373}$, so $|s_i| < 2^{117}$.

## Quadratic LCG to ICG

Paper claims that the explicit shift polynomials used for the Quadratic LCG can be changed to the ICG, but stay the same, they're just rearranged.

Let $v_0$ and $v_1$ be two consectuive outputs from the QLCG, and we know $k$ LSB of the output, so let $w_0$ and $w_1$ be the $k$ LSBs of $v_0$ and $v_1$.

$$v_0 = 2^kx_0 + w_0\\
v_1 = 2^kx_1 + w_1$$
---

## Neopro cooking

from sage.all import *
from problem import *

exec(open('output.txt').read()) # read enc_flag and leaked
leaked = leaked[:64] # truncate for awesomeness

def flatter(M):
    # compile https://github.com/keeganryan/flatter and put it in $PATH
    # if flatter is not available we can just return M.LLL() instead
    from subprocess import check_output
    from re import findall
    z = "[[" + "]\n[".join(" ".join(map(str, row)) for row in M) + "]]"
    ret = check_output(["flatter"], input=z.encode())
    return matrix(M.nrows(), M.ncols(), map(int, findall(b"-?\\d+", ret)))

# solve xy+Ax+By+C=0 (mod N), with |x|,|y|<Δ and using m shifts
# m can be determined automatically if not provided
# follows this paper: https://www.iacr.org/archive/pkc2012/72930609/72930609.pdf
# only returns x because I am lazy
def solve_x(A, B, C, N, Δ, m=None):
    F,(x,y) = ZZ['x,y'].objgens()
    f = x*y + A*x + B*y + C
    
    if m is None:
        t = log(Δ,N)
        m = ceil((1-6*t)/(6*t-2))
        print(f'Automatically picked {m=}')
    
    def make_poly(i, j):
        f_pow = min(i, j)
        x_pow = max(0, i - j)
        y_pow = max(0, j - i)
        N_pow = m - f_pow
        return (f**f_pow * x**x_pow * y**y_pow * N**N_pow)(x = Δ*x, y = Δ*y)

    mat,c = Sequence([make_poly(i,j) for i in range(m+1) for j in range(m+1)]).coefficient_matrix(sparse=False)
    lll = flatter(mat) / diagonal_matrix(vector(c(x=Δ,y=Δ))) * c
    print('reduced', flush=True)
    
    hs = [F(p) for p in lll]
    for j in range(m+1):
        for i in range(j):
            for root,_ in hs[i].resultant(hs[j],y).univariate_polynomial().roots():
                return root
                
# everything below this line is specific to this challenge

N = ZZ(p1 * p2 * p3)
R = crt([a1*(N/p1)**2, a2*(N/p2)**2, a3*(N/p3)**2], [p1, p2, p3])
S = crt([b1*N/p1, b2*N/p2, b3*N/p3], [p1, p2, p3])
# CIG maps x |-> R/x + S

# our equation is:
# (x*2^256 + x_lo) * (y*2^256 + y_lo - S) = R (mod N)
x_lo, y_lo = [int.from_bytes(z,'big') for z in [leaked[:32], leaked[32:]]]
x,y = Zmod(N)['x,y'].gens()
f = (x*2**256+x_lo)*(y*2**256+y_lo-S)-R
f /= 2**512

x = solve_x(*f.change_ring(ZZ).coefficients()[1:], N, 2**116)
print(f'Winner! {x=}')

W0 = mod(x*2**256 + x_lo, N)

W1 = R/(W0-S)
W2 = R/(W1-S)
W3 = R/(W2-S)
xorpad = b''.join((int(w) % 2**256).to_bytes(32,'big') for w in [W3,W2,W1])
print(xor(xorpad, enc_flag))

---

# Bivariate coppersmith of a specific form

We are interested in finding the roots of $f(x,y) = xy + Ax + By + C \pmod{N}$, where $|x|,|y| < \Delta$.

The naive lattice corresponding to this polynomial is

$$
\begin{pmatrix}
\Delta^2 & A\Delta & B\Delta & C \\
0 & N\Delta & 0 & 0\\
0 & 0 & N\Delta & 0\\
0 & 0 & 0 & N
\end{pmatrix},
$$

which has determinant $N^3 \Delta^4$, and thus a shortest vector of norm roughly $N^{3/4}\Delta$. For this to give a solution in $\mathbb{Z}$, we need this norm to be less than $N$, resulting in $\Delta < N^{1/4}$.

Now the aim is to add Coppersmith-style shifted polynomials to the lattice so that we can go up to $\Delta = N^{116/372} \approx N^{0.312}$. The paper claims that we can asympotically achieve $N^{1/3}$, and the purpose of this treatise is to work out how this works.

Let's first attempt to multiply by $x$, i.e. add the polynomial $x^2y + Ax^2 + Bxy + Cx$ to the matrix.

Let's do this ordering: $x^2y, x^2, xy, x, y, 1$.

$$
\begin{pmatrix}
\Delta^3 & A\Delta^2 & B\Delta^2 & C\Delta & 0 & 0 \\
0 & N\Delta^2 & 0 & 0 & 0 & 0 \\
0 & 0 & \Delta^2 & A\Delta & B\Delta & C\\
0 & 0 & 0 & N\Delta & 0 & 0 \\
0 & 0 & 0 & 0 & N\Delta & 0 \\
0 & 0 & 0 & 0 & 0 & N
\end{pmatrix},
$$

which has determinant $N^4\Delta^9$, and thus a shortest vector of norm roughly $N^{2/3}\Delta^{3/2}$. This gives us $\Delta < N^{2/9}$ which is worse than the naive version...

Alright, let's go all-in and do everything up to $x^2y^2$.

We get $f^2 = x^2y^2 + 2Ax^2y + 2Bxy^2 + A^2x^2 + 2(AB+C)xy + B^2y^2 + 2ACx + 2BCy + C^2$.

$$
\begin{pmatrix}
\Delta^4 & 2A\Delta^3 & 2B\Delta^3 & A^2\Delta^2 & 2(AB+C)\Delta^2 & B^2\Delta^2 & 2AC\Delta & 2BC\Delta & C^2\\
0 & \Delta^3 & 0 & A\Delta^2 & B\Delta^2 & 0 & C & 0 & 0\\
0 & 0 & \Delta^3 & 0 & A\Delta^2 & B\Delta^2 & 0 & C & 0\\
0 & 0 & 0 & N\Delta^2 &0&0&0&0&0\\
0 & 0 & 0 & 0 & \Delta^2 & 0 & A\Delta & B\Delta & C\\
0&0&0&0&0&N\Delta^2&0&0&0\\
0&0&0&0&0&0&N\Delta&0&0\\
0&0&0&0&0&0&0&N\Delta&0\\
0&0&0&0&0&0&0&0&N\\
\end{pmatrix},
$$

which has determinant $N^5\Delta^{18}$ for $\Delta < N^{2/9}$ which is the same as before.

But... if we add in $xyf$ in there, we can remove a factor of $N$, which makes it $N^4\Delta^{18}<N^9 \implies N^{5/18}$. (This corresponds to the $m=2$ case.)

Yeah, this seems to work pretty well. If we push this up to $m=7$ as per the paper, then we get a 64x64 matrix which will give us $\Delta < N^{5/16}$.