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 = \begin{cases}x_i + 1 (i = k),\\x_i (i \neq k)\end{cases}

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!}.