Define Function M: $M(x_0,\ldots,x_{n-1}) = \sum\limits_{i=0}^{n-1} \sum\limits_{j=i}^{n-1} q_{i, j}x_i x_j + \sum\limits_{i=0}^{n-1} u_ix_i + c $. ($q$ is for quad, $u$ is for uni. Notice that $q_{i,j}=q_{j,i}$)

The code gives $M(input)$ and  $M(input + flag)$. Let’s think about $F(x, y)=M(x + y) - M(x)$.

$$F(x, y)=M(x + y) - M(x) = \sum\limits_{i=0}^{n-1} \sum\limits_{j=i}^{n-1} q_{i, j} (x_i y_j + y_i x_j + y_i y_j) + \sum\limits_{i=0}^{n-1} u_i y_i $$

For some integer $k$ that $0 \leq k \leq n-1$, define $x’$:

$$ x’_i = x_i + 1 (\text{when}\ i = k), x_i (\text{when}\ i \neq k) $$

Then,

$$F(x’, y)-F(x, y)=\sum\limits_{i=0}^{k} q_{i, k} y_i + \sum\limits_{j=k}^{n-1} q_{k, j} y_j = q_{k,k} y_k+\sum\limits_{i=0}^{n-1}q_{i, k}y_i (\because q_{i,j}=q_{j,i})$$

Note that $q_{k,k}$ appears twice.

Therefore, we can get $n$ linear equations of $flag$ by getting a value of $F(x’, flag)-F(x, flag)$ for each $k$, and it’s clear that they are solvable.

So, do the row reduction and get the flag. The code is here:

from pwn import *

r = remote('mq.eatpwnnosleep.com', 12345)

mq_str = r.recvuntil('\n')
mq_var_str = mq_str.replace(' ', '').split('+')

quad = [ [ 0 for _ in range(32) ] for _ in range(32) ]
uni = [ 0 for _ in range(32) ]
c = 0
p = 131

# Parse MQpoly
for var in mq_var_str:
    tmp = var.split('x')
    if len(tmp) == 1:
        c = int(tmp[0])
    elif len(tmp) == 2:
        uni[ int(tmp[1])-1 ] = int(tmp[0])
    else:
        x, y, z = int(tmp[0]), int(tmp[1]), int(tmp[2])
        quad[y-1][z-1] = x
        quad[z-1][y-1] = x

# Get values
def get_value(s):
    r.send(s)
    tmp = r.recv()
    print(tmp[:4])
    return (int(tmp[2:4], 16) - int(tmp[:2], 16)) % p

base_input = [ 'a' for _ in range(32) ]
base_value = get_value( ''.join(base_input) )
value = [ 0 for _ in range(32) ]

for i in range(32):
    base_input[i] = 'b'
    value[i] = get_value( ''.join(base_input) )
    value[i] = (value[i] - base_value) % p
    base_input[i] = 'a'

print("value:", value)

# Solve (Gauss Elimination)
def xgcd(b, a):
    x0, x1, y0, y1 = 1, 0, 0, 1
    while a != 0:
        q, b, a = b // a, a, b % a
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return  b, x0, y0

def mulinv(b, n):
    g, x, _ = xgcd(b, n)
    if g == 1:
        return x % n

def swap(i, j):
    for k in range(32):
        quad[i][k], quad[j][k] = quad[j][k], quad[i][k]

def mult(i, m):
    for k in range(32):
        quad[i][k] *= m
        quad[i][k] %= p

    value[i] *= m
    value[i] %= p

def reduc(i, j, m):
    for k in range(32):
        quad[j][k] -= quad[i][k] * m
        quad[j][k] %= p

    value[j] -= value[i] * m
    value[j] %= p

for i in range(32):
    quad[i][i] *= 2 # q_{k,k} added twice in the eq
    quad[i][i] %= p

for i in range(32):
    for j in range(i, 32):
        if quad[j][i] != 0:
            break
    if i != j:
        swap(i, j)

    pivot_inv = mulinv(quad[i][i], p)
    mult(i, pivot_inv)

    for j in range(i+1, 32):
        reduc(i, j, quad[j][i] % p)

for i in range(31, -1, -1):
    for j in range(i-1, -1, -1):
        reduc(i, j, quad[j][i] % p)

print(quad)
print(value)

flag = ''.join([ chr(v) for v in value ])
print(flag)

The flag is SCTF{Try_MQ_be_happy!}.