I have a set of $N$ equations and $N$ unknown variables $x_{i}$ like $0 = x_{0}x_{1} - 1 - x_{0}x_{2}$. That's an example. Realistically the equations will be much larger. Each $x_{i} \in \{0, 1\}$. Let $X$ be the set of all $x_{i}$.
The $i^{th}$ equation takes the form $C_{i} = a_{i,0} + a_{i,1} + \ldots + a_{i,m}$. The number of terms $m$ is not necessarily the same for each equation. For each equation, $C_{i} \in \{0, 1\}$ and each term $a_{i,j}$ is either $0$, $-1$, $1$, or a product of some subset of terms from $X$. This product may also be negative, like $-x_{0}x_{2}$ in the example at the beginning.
I'm looking for some ideas to solve this system of equations (which will have a solution by the way). Both theoretical and practical ideas are welcome. I've tried using sympy and symengine to express the equations, but both are extremely slow to even simplify/expand these equations, let alone solve them. The size of the problem I'm dealing with is $N$ ~ 200, $m$ ~ 1,000,000. Probably I will end up coding something in C or C++, but the idea of writing a non-linear equation solver for a problem of this complexity is daunting. SAT solvers like ortools cannot really handle a problem of this size, so I believe I need a critical insight to simplify the problem.
Also when I say "solve" the problem, I'm really looking for any solution, not all solutions.
You can linearize any products of binary variables. See https://or.stackexchange.com/questions/37/how-to-linearize-the-product-of-two-binary-variables for the details.