Motivation
Normally, given an equation \(f(x) = 0\), one can find all roots by WolframAlpha or some magical algorithms. However, the congruence version \(f(x) \equiv 0 \pmod n\) is a lot harder. In the context of CTF crypto challenges, a common trick is to somehow find another relation \(g(x) = 0\) and calculate the gcd of \(f\) and \(g\) (which has degree 1 most of the time). This can be generalized to multi-variable case: given \(n\) variables and \(n + 1\) congruence equations with the same modulus, it is possible to solve the equations in time constant to \(n\) (assuming each arithmetic operation takes \(O(1)\) time).
In the case of two variables, I found an idea to make things simpler. It might be a well-known trick, and in that case, this article serves as a tutorial :p
Naive method
Given \(f(x, y) = g(x, y) = h(x, y) = 0\), we can treat \(x\) as a constant and \(y\) as the variable. After calculating \(\gcd(f, g)\) and \(\gcd(g, h)\), we get two equations of \(x\), which reduces to the one-variable case.
However, the issue is that the coefficents grow exponentially in the process of gcd calculation. Specifically, consider two equations
\[A = a_n(x)y^n + a_{n-1}(x)y^{n-1} + \dots + a_1(x)y + a_0(x)\] \[B = b_n(x)y^n + b_{n-1}(x)y^{n-1} + \dots + b_1(x)y + b_0(x)\]
where all \(a_i\) and \(b_i\) has bounded degree \(d\). To cancel the \(y^n\) term, we calculate \(b_n(x)A - a_n(x)B\), and in the worst case, the coefficeints of the result have degree \(2d\). As there are \(n\) iterations, the final result has degree \(d2^n\), which is unacceptable.
Another approach
We still want to calculate the gcd, but we need something better than the Euclidean method's exponential-size result.
Given \(f(x, y) = g(x, y) = 0\), consider the polynomial ring divided by \(g(x, y)\). By division, we set \(y\) as the major variables, meaning that if the degree of \(g(x, y)\) with respect to \(y\) is \(d\), then every \(y^d\) (or greater) term must be substituted. The idea is that for \(0 \le i < d\), \[f(x, y)y^i = f^i_{d-1}(x)y^{d-1} + \dots + f^i_{1}(x)y + f^i_0(x)\] can be seen as a vector of polynomials \(v_i = (f^i_0, f^i_1, \dots, f^i_{d-1})\). These \(d\) vectors span a space, and if a vector \((t(x), 0, 0, \dots, 0)\) is in the space, then \(t(x)\) is the gcd we want. This motivates us to put these \(d\) vectors in a \(d \times d\) square matrix and perform Guassian elimination, but Guassian elimination has the same issue as Euclidean method – the coefficent growth. Instead, we calculate the determinant of the matrix, which is a multiple of \(t(x)\), by Barieiss algorithm to avoid exponential-size calculation.
The last detail is the calculation of \(f(x, y)y^i\). Notice that while \(f(x, y)\) has bounded degree \(d\), \(f(x, y)y\) can have higher degree because when the \(y^n\) term exists, we need to perform a division. The good news is that this process has linear degree growth, so \(f(x, y)y^i\) has degree at most \(id\), and the determinant will have degree \(O(d^2)\).
Remark
The above method utilizes a lot of polynomial arithmetic operation. I found sympy galois tools useful, but they probably didn't implement FFT because it is slower than what I thought. Sage might have relevant tools, but most of them require \(n\) to be a prime. Ironically, in RSA setting, if they found out that a number does not have a modular inverse modulo \(n\), that would be even better.
I also didn't think about the complexity of the above method but just make a rough estimation of \(O(d^4)\). FFT can improve the complexity, but the purpose of the article is just to share this idea. Also, it is still a problem to generalize this to \(n\) variables as now it only works with constant number of variables.