# How to solve a system of linear equations with binary variables mod 2 with a constraint that is not mod 2?

I am trying to solve a system of $$H \leq 4^n$$ polynomial equations of degree $$H$$, where the variables $$(x_1, x_2, \ldots, x_{2n}) \in \mathbb{Z}_2^{2n}$$ are binary. The problem was that these equations had a mixture of regular sums and sums mod 2. I have used linearisation to convert my non linear problem of $$2n$$ variables to a linear problem of up to $$2^n -1$$ variables:

As $$x_i$$ is binary, $$x_i^k = x_i$$ and I have let $$u_{ij} = x_i x_j$$ subject to $$-x_i + u_{ij} \leq 0$$, $$-x_j + u_{ij} \leq 0$$ and $$x_i + x_j + u_{ij} \leq 1$$. For products of 3 or more variables, I have analogously let $$u_{ijk \ldots} = x_i x_j x_k \ldots$$ subject to similar constraints. This brings me back to my original problems was that there were mixtures of sums mod 2 and sums that were not mod 2. This linearisation scheme allowed me to form:

$$A\underline{u} ~~~\text{mod 2}= \underline{0}$$

where $$\underline{u} = (u_1, u_2, \ldots, u_{11}, \ldots, u_{ijk \ldots}, \ldots)^T$$ where $$u_{ijk \ldots} = x_i x_j x_k \ldots$$, $$A$$ is an $$m \times (2^n -1)$$ matrix (known), and $$\underline{0}$$ is an $$m$$-dimensional vector of $$0$$s. All the linear equations used to form the matrix equation must be done mod 2. However the constraints from linearisation are not mod 2, and are of the form:

$$B\underline{u} \leq \underline{c}$$

where $$B$$ is a $$(2^n-1) \times (2^n-1)$$ matrix and $$\underline{c}$$ is a $$2^n -1$$ dimensional vector. My question is really how would you go about solving this? I have heard of the Euclid-Wallis algorithm from searching and then I could use those solutions to check if they satisfy the constraints. Does anyone know of any solvers that may be able to solve this but are there any known solvers (preferably in Python) that would be able to do this?