Limitations of the Macaulay matrix approach for using the HHL algorithm to solve multivariate polynomial systems

our improved algorithm using a reduction employing the Valiant-Vazirani aﬃne hashing method, and also extend the result to polynomial systems over F q improving on subsequent work by Chen, Gao and Yuan [CGY18]. We also suggest a new approach for extracting the solution of the Boolean polynomial system via a generalization of the quantum coupon collector problem [ABC + 20].


Introduction
Solving systems of multivariate polynomial equations is a fundamental problem that is NPcomplete even when the polynomials are restricted over F 2 . The problem can be reduced to solving an exponential number of linear equations via the so-called Macaulay matrix, which holds coefficients of linear equations that come from the input polynomials, and multiples of them (multiplying each polynomial by each monomial up to a certain degree). Each monomial is represented by a new variable, recasting the polynomial equations and their multiples as linear equations. The usual classical approach to solve a polynomial system is based on computing the Gröbner basis of the corresponding polynomial ideal by triangularizing the Macaulay matrix. There is a vast literature on characterizing and improving the complexity of solving various types of polynomial systems using the Macaulay matrix [AFI + 04, CG17, CKPS04, DS13, Die04, Per16, WW15].
In quantum computing, the HHL [HHL09] Quantum Linear System (QLS) algorithm outputs a quantum state |x⟩ such that µA |x⟩ = |b⟩ for an exponentially large matrix A with certain properties, and a quantum state |b⟩, in time O κ 2 s 2 , 1 where κ is the condition number of A, µ is a normalization factor, and s is the sparsity of the matrix A, while state-of-the-art QLS algorithms [Amb12, CKS17, CGJ19, GSLW19, SSO19, LT20] have complexity O (κs). Although, the QLS algorithm is BQPcomplete [HHL09], meaning that it captures all essential features of quantum computing, a natural "killer-application" is still to be discoveredshowing the difficulty of finding a practically interesting problem instance that satisfies all stringent conditions. For example, to efficiently solve the classical equation Ax = b using the original HHL algorithm, where implicit access is given to an exponentially large matrix A and b, the following must be satisfied: the state |b⟩ can be efficiently prepared, the sought data can be efficiently extracted from the output state |x⟩, and the matrix A should be sparse and wellconditioned [Aar15].
Chen and Gao [CG21] made an interesting connection between the exponential size Macaulay matrix and the HHL algorithm. While they use Gröbner bases in their proof of correctness, they do not explicitly compute the Gröbner basis and instead use the HHL algorithm to solve the exponentially large system of linear equations resulting from the Macaulay matrix. They show that the access requirements that usually cause so much trouble, can all be resolved for this application, namely: they can efficiently compute the entries of an appropriate sparse matrix A, prepare |b⟩, and extract the answer from |x⟩. However, a major question was left open: what is the condition number of the matrices, driving the running time? Intuitively, for worst case instances of polynomial systems, the condition number of the resulting matrix should be large because the approach can solve NP-complete problems. This being said, the analysis of the condition number was left open, both in general, and for special cases such as breaking cryptosystems which have distributions over specific problem instances that might be easier than the worst case. Therefore, the algorithm of Chen and Gao [CG21] together with the follow-up work of Chen, Gao, and Yuan [CGY18] presented a potential quantum threat on multivariate cryptosystems. However, there was no consensus on the strength of this potential quantum attack, as its cryptanalysis was wide open.
In this paper we prove an exponential lower bound on the condition number κ of the matrix A related to the Boolean polynomial system, which shows that the quantum algorithm in [CG21] takes exponential time in the worst case. We also give a Grover-based brute-force search algorithm that outperforms their quantum algorithm for solving Boolean polynomial systems when there is a unique solution or all solutions have the same Hamming-weight. Specifically, in the unique solution case we give a simple proof that the condition number κ is Ω((3n) h/2 ), where h is the Hamming weight of the solution to the original n-variable Boolean polynomial system. Meanwhile, a simple Grover-based brute-force search algorithm over the possible assignments to the variables takes time O In fact, we give "robust" lower bounds on the condition number by also considering "truncated" QLS algorithms [HHL09, GSLW19]. Namely, if the singular-values of A are only inverted on a well-conditioned subspace and the overlap of the solution x with such a subspace is large enough, then a "truncated" QLS algorithm can provide a sufficiently accurate solutionx. In order to give a bound on the performance of such "truncated" versions of the QLS algorithm, we define the concept of the truncated QLS condition number κ b (A) := ∥A∥ A + b / ∥b∥, 2 which is also a lower bound on κ = ∥A∥ A + . All of our lower bounds also apply to the truncated QLS condition number, ruling out further improvement by truncated QLS algorithms. These results provide strong evidence that the quantum algorithm of [CG21] (at least in its original form) does not present a fatal cryptanalytic threat, and give generic tools for analyzing the strength of individual cryptosystems against this type of quantum attack.
Finally, we refine Chen and Gao's algorithm to the point that our lower bound does not always rule out the possibility of a superpolynomial quantum speedup even for unique solu-tions. In particular, the lower bound changes from (3n) h/2 to 2 h/2 on our refined algorithm, so for h = Θ(log n) the lower bound is only a polynomial, while the brute-force algorithm takes quasipolynomial time. Thus, it is conceivable that the condition number is upper bounded by poly(n) for some set of interesting input equations, potentially yielding a superpolynomial speedup. We leave it open to find a problem instance whose associated Macaulay matrix has a small enough condition number so that the running time of our refined quantum algorithm gives a speedup over the best classical or Grover-based algorithm. Such an example could result in a new type of quantum speedup and one that uses the HHL algorithm in a novel way.
The core ingredient of our refined algorithm is to show that the Macaulay matrix can be simplified to what we call the Boolean Macaulay matrix over C by exploiting that the input consists of quadratic polynomials over C, but restricted to 0/1 solutions. The Boolean Macaulay matrix is a submatrix of the original Macaulay matrix that can be obtained via Gaussian elimination over C. This construction of the Boolean Macaulay matrix over C is different from the Boolean Macaulay matrix over F 2 as defined in [BFSS13] since they are over different fields. This matrix preserves the solution set while its size is much smaller compared to the original Macaulay matrix -ultimately leading to a smaller lower bound Ω(2 h/2 ) on the condition number.
The correctness of our refined algorithm can be shown via the equivalence between the Boolean Macaulay linear system and the Macaulay linear system, where the correctness of the latter has been proven in [CG21]. For completeness, we also provide a simple self-contained proof of correctness for our improved algorithm in Appendix A, which is more elementary than that of the original algorithm proposed by Chen and Gao [CG21], since our proof does not require Gröbner bases. Instead, our proof combines a special case of the reduction in [CGY18] with the Valiant-Vazirani affine hashing method, reducing any Boolean polynomial system with more than one solution to one that has a unique solution. For a Boolean Macaulay linear system that has a unique solution, we also provide an alternative approach for extracting the Boolean solution of the corresponding Boolean polyno-mial system from the output quantum state, i.e., the normalized monomial solution vector of the Boolean Macaulay linear system. Specifically, we reformulate this problem as a generalization of the quantum coupon collector problem, and prove that O(log n) iterations suffice for extracting a solution, whereas Chen and Gao's algorithm uses O(n) iterations. On the other hand, the affine hashing reduction introduces O (n) extra rounds, so the total number of iterations in our algorithm can be bounded as O (n log n).
2 Quantum algorithms for solving polynomial systems Chen and Gao [CG21] proposed using the HHL algorithm to solve a Boolean polynomial system via solving an exponentially large linear system of equations. Now we discuss the two main parameters that appear in the complexity of the QLS algorithm, and their relevance in our case.  only bound the complexity of standard QLSP solvers, but also lower-bounds the complexity of fine-tuned truncated variants.
s : In order to efficiently solve the Quantum Linear System Problem (QLSP), we need a succinct representation of the vector b and the matrix A. In our case both b and A are s-sparse for some polynomially large s, i.e., they have at most s nonzero entries in every column and row. In order to utilize their sparsity, we also need to be able to efficiently compute the locations and the values of the nonzero entries; this is easy to do in our case, since the vector b is sparse and the (Boolean) Macaulay matrix A has a quasi-Toeplitz structure.
More precisely, it suffices to have efficient (quantum) circuits computing the locations and values of the nonzero elements of A, allowing us to perform the transformations where i labels the row indices of the matrixĀ standing for A and A T , k labels the nonzero entries ofĀ (which is assumed to be s-sparse), c(i, k) represents the column index of the k-th nonzero entry of the matrixĀ in row i, and |0⟩ in (2) represents a (large enough) ancillary system in which the matrix elementĀ ij can be stored (as a bit-string with sufficient precision so that errors can be neglected Thereby, the QLSP problem can also be efficiently solved for non-sparse b or A, whenever we can build efficient quantum circuits C b , C A . For example, if we have both b and A stored in a quantum random access memory (qRAM) using an appropriate data structure, then we can build the efficient quantum circuits C b , C A as described in [KP17, CGJ19, GSLW19].
For completeness, let us mention two additional types of quantum algorithms that have been proposed for solving polynomial systems:  [BFSS13], that applies Grover search on both the exhaustive search and the consistency check subroutines. Under certain assumptions, the running time of this quantum algorithm is slightly better than O 2 n/2 , which is the running time of trivial Grover search algorithm. Similarly, Bernstein and Yang [BY18] gave a quantum algorithm (GroverXL) for random polynomial systems over a finite field While QAOA is only a heuristic algorithm and Grover search only represents a quadratic speedup, the QLS algorithm -being BQPcomplete -captures the power of quantum computation and promises rigorous superpolynomial speedups, giving strong motivation for the study of the Macaulay matrix approach. 3 Reducing polynomial system solving over a finite field F q to polynomial system solving over C First, we present the F q = F 2 special case of Chen, Gao and Yuan's approach [CGY18]. In this special case there is a bijection between the 4 The generalization of the Boolean Macaulay matrix method in [BFSS13] is equivalent to the reduced XL approach that first appeared in [CKPS04] and then defined in [Die04]. solution sets of the corresponding polynomial systems over F 2 and C.
Problem 3.1. Solve a system of n-variate quadratic polynomials with Boolean variables over F 2 .
Problem 3.2. Solve a system of n-variate quadratic polynomials over C, together with the F 2 field equations that force variables to be Boolean.
Let #F denote the maximum number of nonzero terms in any polynomial in F. Also let #f be a shorthand for #{f }.

Lemma 3.3.
There is a polynomial-time reduction from Problem 3.1 on n variables and a set of m equations F to Problem 3.2 on n + m · ⌊log 2 #F⌋ variables and n + m · (⌊log 2 #F⌋ + 1) equations.
Proof. Solving Problem 3.1 is equivalent to solving the following polynomial system in variables x 1 , . . . , x n , z 1 , . . . , z m over C: ∀i ∈ [m] : ∀j ∈ [n] : The field equations (5) force each x j to be 0 or 1, and therefore each f i evaluates to an even or odd integer. Also, it is easy to see that each integer z i in equation (4) is in [0, #f i ] for i = 1, . . . m. and can be treated as a polynomial. Equations (3)-(4) force f i to evaluate to an even integer, so it is 0 mod 2.
Chen, Gao and Yuan [CGY18] then represent each z 1 , . . . , z m by the bits in its binary expansion. For each variable z i ∈ [0, #f i ] a polynomial is introduced with Boolean variables y ib to represent its value in binary, i.e., z i = ⌊log 2 #f i ⌋ b=1 2 b y ib , and y 2 ib − y ib = 0 for b = 1, . . . , ⌊log 2 #f i ⌋ .
Substituting the polynomials and Boolean constraints corresponding to each z i into the polynomial equations (3)-(5), we get a following polynomial system F over C: ∀j ∈ [n] : It is easy to see that there is a bijection between the set of solutions of F ⊆ F 2 [x 1 , x 2 . . . , x n ] and the set of solutions of (6)-(8) over C. On one hand, given a solution (s 1 , . . . , s n ) of F ⊆ F 2 [x 1 , x 2 . . . , x n ], evaluating f i (s 1 , . . . , s n ) over C gives an even number z i , and its binary expansion gives the values of the y ib variables. On the other hand, let (s j ), (t ib ) be a solution to (6)-(8). For each j, s j ∈ {0, 1} by (8), and for each i, f i (s 1 , . . . , s n ) = ⌊log 2 #F ⌋ b=1 2 b y ib by (6). Due to (7) this implies that f i (s 1 , . . . , s n ) ≡ 0 mod 2, so (s j ) is a solution of F.
Given a set of polynomials F ⊆ F q [X], Chen, Gao and Yuan [CGY18] propose an analogous approach for reducing the polynomial system F to a polynomial system in the form of Problem 3.2. However, in general even if F has a unique solution, the constructed polynomial systems F C may have multiple solutions. Now we present two additional reduction steps that make the polynomial systems in Problem 3.2 easier to handle: Red1: In order to ensure that the polynomial system F ⊆ C[X] in Problem 3.2 has no more than one solution, we employ the Valiant-Vazirani affine hashing method [VV86]. Suppose that the polynomial system F has S ∈ [2 n ] different solutions. The main idea of the affine hashing method is the follow- , then they isolate a unique solution with probability at least 1 8 . Even if we don't know the number of solutions a priori, we can loop over all possible values of ⌊log 2 (S)⌋ ∈ {0, 1, . . . , n}; making O (ln(1/ε)) trials for all possible choices of ⌊log 2 (S)⌋ gives at least one system F R = 0 with a unique solution with probability at least 1−ε. This amounts to O (n log(1/ε)) different polynomial systems to check. Remember that F R ⊆ F 2 [X] whereas F ⊆ C[X]. By Lemma 3.3, we can reduce the new polynomials F R to polynomials F RC ⊆ C[X, Y ], where Y is the set of new variables introduced during the reduction. Finally, if F R isolated a unique solution of F, then the polynomial system F RC ∪ F has a unique solution. Thus, without loss of generality, we can always assume that Problem 3.2 has a unique solution. Red2: . . , f ′ m have no constant terms. In case no polynomial in F has a constant term, the all-zero vector is a trivial solution. Otherwise, let c i denote the constant term of f i , and let us assume without loss of generality that c 1 ̸ = 0. Then we can simply set f ′ 1 := −f 1 /c 1 , and . . , m}. The above two reductions increase the parameters considered in this paper only moderately. Indeed, Red1 introduces at most O (n log(n)) new equations and variables, while Red2 only affects the number of nonzero terms in the polynomial system. Moreover, Red2 increases #F and the total number of nonzero terms #f i by at most a factor of 2 (for the latter we shall choose f 1 to be the polynomial with a nonzero constant term that also has the fewest nonzero coefficients).

Macaulay linear systems and their tQLScn
In this section we define the Macaulay linear system of a set of polynomials F ⊆ C[x 1 , . . . , x n ] and show that when F has a unique solution, the condition number of the matrix is Ω( n h ), where h is the Hamming weight of the solution. We show that the lower bound also holds when using max degree instead of total degree in the definition, as was done in [CG21], showing that their proposed quantum algorithm for solving polynomial equations by using the QLS algorithm to solve a Macaulay linear system in general takes time Ω((3n) h/2 ). We also show that if there are t different solutions, but they have the same Hamming weight h, then the above lower bound can reduce by at most a factor of √ t. Finally, we give a formula that can be used for giving a lower bound on the condition number for any number of solutions and present computational evidence that this analytical lower bound is exponentially large in terms of the smallest Hamming weight among the solutions.

Macaulay linear systems
There is a well-known approach for solving polynomial systems by linearizing them with the help of introducing new latent auxiliary variables. The advantage of this approach is that the problem becomes linear, but the downside is that the new problem is exponentially large. The matrix of the resulting linear system is called the Macaulay matrix.
is the matrix where each row is labeled by a pair of polynomials (m, f ) and contains the corresponding coefficient vector of the polynomialmf . The rows range over all f ∈ F and monomialsm such that mf has degree at most d. The columns are labeled by the set of monomials in x 1 , . . . , x n of degree at most d and are ordered with respect to a specified monomial ordering. The element in the row corresponding to (m, f ) and the column corresponding to the monomialm ′ is the coefficient ofm ′ in the polynomialmf .
In the above definition one can interpret the degree as either the total degree or the max degree (the maximum degree of any variable) of multivariate polynomials resulting in different notions of the Macaulay matrix. When it is necessary we will always clarify which definition is being used. For example, [CG21] uses the max degree, so all references to that paper refer to the max degree version of this definition.
In the classical setting, the goal is to compute the Gröbner basis from the Macaulay matrix, where the Gröbner basis is a set of polynomials G = {g 1 , g 2 , . . . , g r } such that for the leading term of any polynomial f in the ideal I = (F), there exists a polynomial g k ∈ G such that LM (g k )|LM (f ) 5 . Note that the size of the Macaulay matrix depends on the selected degree. 5 Here LM (f ) is the leading term of the polynomial f For a set of m quadratic polynomials with n variables, the degree is approximately lower bounded by n √ m [CKPS04]. When m = αn, the degree is upper bounded by c α n for some constant c α [BFSS13]. In the quantum setting, the goal is to compute the monomials up to a certain degree. In this paper, we only provide the upper bound of the degree that applies to any quadratic polynomials. It might be interesting to check the degree of some special polynomial systems.
Row operations on the matrixM correspond to polynomial addition, subtraction, and scalar multiplication in the polynomial ideal ⟨F⟩ [Bat13, B + 18], and these operations preserve the common roots of the system. Classically, Gaussian elimination can be performed on this matrix, and the entries can be read out from the row-reduced matrix [Die04, CKPS04]. However, in the quantum case, we cannot directly do Gaussian elimination on this matrix and look at the row-reduced matrix. Instead, Chen and Gao showed that the QLS algorithm can be used for sampling from nonzero solutions of the following related linear system. In other words, the Macaulay matrix is the augmented matrix from the Macaulay linear system M⃗ y = ⃗ b. Due to reduction Red2 in the previous section, it can be assumed that exactly one of the input polynomials has a nonzero constant term. So we may assume without loss of generality that ⃗ b = [1 0] T , where the vector ⃗ b corresponds to the column vector indexed by the degree 0 monomial 1.
Chen and Gao's algorithm applies the QLS algorithm to output the quantum state |ŷ⟩ that can be measured in order to sample from monomials with nonzero value in a valid assignment corresponding to a solution. We will lower bound the condition number of M, which will in turn lower bound the running time of the proposed algorithm.
Chen and Gao proposed to set the max degree to 3n in the Macaulay linear system, and they showed with this choice if a set of polynomials F has a unique solution, then the linear system also has a unique solution [CG21, Lemma 4.1]. The output state |ŷ⟩ of the QLS algorithm then corresponds to this unique solution of the linear system, and the solution of F can be efficiently obtained from measuring the state |ŷ⟩.
Classically, the solving degree of F is used, which is at most n + 2 as shown by Caminata and Gorla [CG17,Theorem 3.26]. This allows computing the Gröbner basis of the polynomial ideal ⟨F⟩ via Gaussian elimination of the linear system, and the solution of F can be obtained from the Gröbner basis. However, for some polynomial systems, the affine subspace of all solutions of the linear system has no well-understood structure even though F has a unique solution. In this case, the QLS algorithm outputs a state |ŷ⟩ that corresponds to the smallest ℓ 2 -norm solution of the linear system. In general we don't know how to extract the solution of F from such states |ŷ⟩.
When F has more than one solution and the max degree is set to be 3n, the dimension of the affine subspace of all solutions of the linear system equals the number of solutions of F minus 1. For each solution a ∈ {0, 1} n of F, there is a corresponding solutionŷ a of the linear system. Those solutionsŷ a are linearly independent and any solution of the linear system is an affine combination of the solutionsŷ a [CG21, Lemma 3.18]. Again, the QLS algorithm outputs a state |ŷ⟩ corresponding to the smallest ℓ 2 norm vectorŷ in this affine subspace.
Chen and Gao [CG21] showed that M is an O (m · #F)-sparse, row computable matrix and ⃗ b can be efficiently prepared as a quantum state. In particular, assuming |F| = O (poly(n)), they show that the QLS algorithm can be run in time O (poly(n)κ(M)). There is strong complexity theoretic evidence that in general running the QLS algorithm requires time Ω(κ(M)) [HHL09], so a lower bound on the condition number also lower bounds the running time.

Lower bound on the truncated QLS condition number κ ⃗ b (M)
In this section we give a lower bound on the tQLScn κ ⃗ b (M). Since κ ⃗ b (M) ≤ κ(M), this also implies a lower bound on the time complexity of Chen and Gao's [CG21] algorithm.
In order to prove a lower bound on κ ⃗ b (M), it suffices to lower bound the length of the solution vector ⃗ y = M + ⃗ b, since Here, the first equality is the definition of κ ⃗ b , and the inequality follows from ⃗ b = 1, and because ∥M∥ ≥ 1, as M has at least one matrix element which has absolute value at least 1. 6 In order to understand the length ∥⃗ y∥ of the solution vector ⃗ y = M + ⃗ b, let us first study the monomial solution vector ⃗ y (a) corresponding to a binary solution a. For a degree-d Macaulay linear system, a monomial exponent 0 ̸ = e ∈ N n is a "valid" coordinate of ⃗ y (a) if e ∈ {0, 1, . . . , d} n (and e i ≤ d for total degree), moreover the solution vector satisfies ⃗ y If the Hamming weight of a is h and M is constructed with max degree, then the number of such non-zero coordinates (mono- When M is constructed with total degree, the number of such non-zero coordinates (monomi- Suppose a 1 , a 2 , · · · , a t ∈ {0, 1} n are the t solu- The t solutions a 1 , a 2 , · · · , a t of F must be nonzero because the first equation has constant term b 1 = 1. Let ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t 7 be the 0/1 solution vectors of the linear system M⃗ y = ⃗ b under the assignments a 1 , a 2 , · · · , a t respectively. When the max degree of the linear system is set to be 3n, the affine subspace of all solutions of the linear system is spanned by the monomial solution vectors ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t [CG21, Theorem 3.21 and Lemma 6 For the Macaulay matrix construction, let f be x 2 −x, the row (1, f ) will have a matrix element of magnitude 1. 7 Note that the monomial solution vector corresponding to a binary solution ai is ⃗ y a i . Here and in the following discussion of multiple solutions case, we write ⃗ y a i as ⃗ yi for simplicity. 4.1], but this property might also hold for lower degrees. From now on we assume that degree d is such that the linear system has this property, then ⃗ y = M + ⃗ b has the minimum ℓ 2 norm in the affine subspace spanned by ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t .
If all the t solutions a 1 , a 2 , · · · , a t ∈ {0, 1} n have the same Hamming-weight, we can lower bound the length ∥⃗ y∥ of the solution vector ⃗ y = M + ⃗ b by the following lemma. Lemma 4.3. Suppose ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t are vectors with 0 and 1 entries such that their Hamming weights are equal. Then every vector in their (complex) affine hull A has length (ℓ 2 norm) at Proof. Every entry of the vectors ⃗ y i is either 0 or 1, therefore ⟨⃗ y i , 1 1 1⟩ = ∥⃗ y i ∥ 1 (we denote by 1 1 1 the all-1 vector). Since the vectors have the same Hamming weight we also have ∥⃗ . Let ⃗ y be any vector in A, then ⟨⃗ y, 1 1 1⟩ = ⟨⃗ y 1 , 1 1 1⟩, and in particular ∥⃗ y∥ 1 ≥ ⟨⃗ y, If the minimum ℓ 2 -norm solution ⃗ y =M + ⃗ b happens to be a convex combination ⃗ y = t i=1 w i ⃗ y i of the (possibly differing Hamming weight) solution vectors ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t , then we can similarly lower bound the length ∥⃗ y∥.
Lemma 4.4. Suppose ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t are 0/1 vectors and ⃗ y 1 has the minimum Hamming weight. Then every vector in their convex hull A has The second inequality is true because ⟨⃗ y i , ⃗ y i ⟩ ≥ ⟨⃗ y 1 , ⃗ y 1 ⟩ for any i ∈ [t] as ⃗ y i are a 0/1 vectors and ⃗ y 1 has the minimum Hamming weight. The third inequality follows from Cauchy-Schwarz.
Theorem 4.5. Suppose a 1 , a 2 , · · · , a t ∈ {0, 1} n are the t solutions of F = F 1 ∪ F 2 , where F 1 = {f 1 , . . . , f m } and F 2 = {x 2 1 −x 1 , . . . , x 2 n −x n }, and let h be the minimum Hamming weight of the t solutions a 1 , a 2 , · · · , a t . Let d be the selected degree on constructing the Macaulay linear system M⃗ y = ⃗ b and let ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t be the corresponding solution vectors of the Macaulay linear system M⃗ y = ⃗ b under the assignments a 1 , a 2 , · · · , a t respectively.
If all the t solutions a 1 , a 2 , · · · , a t have the same Hamming weight h or the minimum ℓ 2norm solution vector ⃗ y = M + ⃗ b is in the convex hull of ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t , then the tQLScn of M of F in the Macaulay linear system In particular in the setup in [CG21], using max Now we give a lower bound in terms of the smallest Hamming weight in the binary solution set. For this we use the following purely geometrical lemma.
Note that the problem of finding the length of the shortest vector in an affine subspace can also be reformulated as the following SDP: Tr(Gρ) subject to: Tr(1 1 1·1 1 1 T ρ) = 1, where without loss of generality we can assume that the above optimizer ρ has rank 1.
By the weak duality of SDPs, and utilizing the dual of the above problem, we get: subject to: G − γ1 1 1 · 1 1 1 T ⪰ 0.
In fact the following proof of Lemma 4.6 shows that the above inequality is tight, i.e., strong duality always holds for this SDP.
Proof. The shortest vector s in the affine subspace is orthogonal to all vectors of the form v i − v j , therefore we must have that ⟨s, v i ⟩ is constant for all i ∈ [k], i.e., V † s ∝ 1 1 1. Since s is in the column space of V , we can write it in the form If s ̸ = 0 then this implies that 1 1 1 is in the column space of G, consequently 0 ̸ = 1 1 1, G + 1 1 1 and thereby s = V · w for w = G + 1 1 1/ 1 1 1, G + 1 1 1 . 8 From this we can conclude γ * = ∥s∥ 2 = ⟨w, Gw⟩ = 1/ 1 1 1, G + 1 1 1 . Conversely, if 1 1 1 is in the column space of G then s = V · w for w = G + 1 1 1/ 1 1 1, G + 1 1 1 is a nonzero vector, which is the shortest vector in A due to the fact that it is orthogonal to all vectors of the form v i − v j . We can conclude that γ * ̸ = 0 iff 1 1 1 is in the column space of G.
Suppose that the Boolean solution set of the polynomial system is S = {a 1 , a 2 , · · · , a t }. Let A S be the affine subspace corresponding to the solution set S, spanned by the monomial solution vectors y 1 , y 2 , · · · , y t of the linear system M⃗ y = ⃗ b corresponding to the Boolean solutions a 1 , a 2 , · · · , a t respectively. We wish to lower bound the length of the shortest vector in A S ; for this it suffices to find the length of the shortest vector in an enlarged affine subspace A S ′ ⊇ A S .
Let S ′ be the symmetrized solution set of S by applying all possible permutations of the variables of a i 's and taking their union. Let A S ′ be the affine subspace corresponding to the solution set S ′ and let v be the shortest vector in A S ′ . For each Hamming weight h that appears in S ′ , there is a symmetrized monomial solution vector v h . This v h equals the average over all monomial solution vectors that are associated to Boolean solutions of Hamming weight h.
Next we will argue that the minimum ℓ 2 -norm vector v ∈ A S ′ is an affine combination of the symmetrized monomial solution vectors v h . If we apply an induced permutation on the coordinates of v according to a permutation of the Boolean variables, then the ℓ 2 -norm of the resulting vector u ∈ A S ′ is equal to the ℓ 2 -norm of v. Because v has the minimum ℓ 2 -norm in A S ′ we have u+v 2 2 ≥ ∥v∥, and due to ∥u∥ = ∥v∥ by the triangle inequality u+v 2 2 ≤ ∥v∥ (equality holds if and only if u = v), the resulting vector u is equal to v. Therefore, the shortest vector v is invariant under all the possible induced permutations, therefore we can conclude that v is an affine combination of the symmetrized monomial solution vectors v h . Now, we can lower bound M + F b by finding the lowest ℓ 2 norm of a vector in the affine subspace spanned by the symmetrized vectors v h corresponding to the Hamming weights h that appear in S. This can be achieved by considering the Gram matrix as explained in Lemma 4.6.
In order to compute this Gram matrix, we need to understand the symmetrized vectors v h . For this, let us introduce the following orthonormal vector system (b s ), corresponding to the set of monomials m d s that contain exactly s vari- and the Gram matrix G of the symmetrized vectors has matrix elements Together with Lemma 4.6 this enables us to give a lower bound on the smallest ℓ 2 -norm solution in terms of the minimal Hamming weight appearing in the solution set as follows.

Theorem 4.7. Suppose that F is a Boolean polynomial system with n Boolean variables where each solution have Hamming weight at least h, and d ≥ n.
Recall that tQLScn, defined as κ ⃗ b (M) = ∥M∥ ∥M +⃗ b∥ ∥ ⃗ b∥ , is also lower bounded by the smallest ℓ 2 -norm solution by equation (9). Then the degree-d Macaulay linear system's tQLScn is lower bounded by where G ∈ R n×n is the Gram matrix whose (i, j) matrix element is c d s = d s for max degree, while c d s = d s for total degree, and finally G (h) is the bottom-right (n − h + 1) by (n − h + 1) minor of G. 9 The expression in (15) is difficult to bound analytically, but it appears to be exponentially large in terms of h for large enough d. In particular, we could verify 10 that for max degree d = 3n Equation (15) is lower bounded by h h /2 for every h ∈ [n] up to n = 300.

Comparison to brute-force search
In case there is a unique Boolean solution we showed that the lower bound of the running time of the quantum algorithm using the HHL algorithm is exponential in the Hamming weight of the unique Boolean solution, and we provided strong evidence that this is also true when there are multiple solutions.
It is useful to compare the HHL-based approach to classical brute-force search and also to using Grover's algorithm. In case we know that the unique solution has Hamming weight h, we can simply classically search through all the n h different Hamming-weight-h assignments of the original polynomial system. We can also use Grover search to find such an assignment with O When d ≥ n, due to the triangular shape of the nonzero coefficients of the vectors in (14) it is easy to see that G has full rank, i.e., it is positive definite. It follows, that all principal submatrix of G are also positive definite, i.e., have full rank, therefore G (h) is invertible.  evaluations of the polynomials. By comparing this to the lower bounds of Theorem 4.5 one can see that in case d+h ≥ n Grover's algorithm performs at least as good as the HHL based algorithm (up to some potential lower order correction 4 √ h), and the algorithm of Chen and Gao [CG21] where the max degree d = 3n is definitely outperformed by Grover search.
In case there are multiple solutions, but all their Hamming weights are the same, Theorem 4.5 ensures that we do not get a bigger reduction in the condition number than the analogous speedup we can already achieve by plain Grover search. So the above Grover-based algorithm still performs just as competitively.
In the general case of having multiple solutions with different Hamming weights, the situation is harder to analyze, but we could still obtain an exponential lower bound on tQLScn in terms of the smallest Hamming weight solution up to n = 300, providing strong evidence for Chen and Gao's algorithm [CG21] having a best case complexity that is exponentially large in terms of the minimal Hamming weight of a solution, making it unlikely that their algorithm would give a substantial improvement over brute-force Grover search.

The Boolean Macaulay linear system and its tQLScn
In this section we give an equivalent but more efficient way to represent the Macaulay matrix using the fact that we are only searching for 0/1 solutions in C. This results in a smaller lower bound on the tQLScn of size Ω(2 h/2 ). While the quantum algorithm's running time is still exponentially large for larger Hamming weight so-

The Boolean Macaulay matrix over C
In this section we again include the field polynomials F 2 = {x 2 1 − x 1 , . . . , x 2 n − x n } for the field F 2 together with the input polynomials F 1 . Solving the system F = F 1 ∪ F 2 forces the roots to be effectively Boolean even though the underlying field is C. This allows all monomials in an equation to be replaced with equivalent multilinear versions and a reduced Macaulay matrix will be defined that has a more compact form. This was done in [BFSS13] for finite fields, where the extra equations prevented solutions from being in field extensions. We derive the analogous matrix when the solutions are forced to be Boolean but the arithmetic is over C. Additionally, in our case (similarly to Chen and Gao's original construction [CG21]) the structure of the Boolean solutions makes it possible to extract the Boolean solutions from measuring the quantum state corresponding the the solution vectors (over C).
Let ψ : R → R map a monomial to its mul- . Lemma 5.3 will show that having max degree higher than 1 becomes redundant, so the following definition only has rows up to max degree 1, and in this section we will set the total degree d = n, the number of variables. For notation, let m andm ′ denote monomials, and let m, m ′ , and m ′′ denote multilinear monomials (i.e., monomials with max degree at most 1). Note that compared to the Macaulay matrix, in addition to forcing answers to be Boolean, the Boolean Macaulay matrix is reduced in a certain way, by eliminating polynomials with max degree at least 2. Next we will show that the Boolean Macaulay matrix can be obtained as a submatrix of the Macaulay matrixM of max degree d corresponding to the set of polynomials F = F 1 ∪ F 2 after Gaussian reduction on the rows.
First we consider the special case when F 1 = ∅ and perform the row reduction on the Macaulay matrix of F 2 to show that the field equations F 2 take a special form.
Lemma 5.2. Let the Macaulay matrixM 2 of max degree d of F 2 have its columns ordered such that they are partitioned into two parts the following way: the labels on the right side are multilinear monomials (including the degree 0 monomial 1) and ordered in ascending order with respect to the integer represented by the exponent vector of the multilinear monomial, and let the left side columns be labeled by nonmultilinear monomials and ordered under any monomial order.
Then using row operations, where I 2 is the identity matrix of dimension (d + 1) n − 2 n with rows and columns labeled by nonmultilinear monomials, and rows with zeros are removed.
Proof. The rows ofM 2 are indexed by pairs of polynomials (m, x 2 j − x j ), wherem is a monomial and maxdegm(x 2 j − x j ) ≤ d. The approach is to first change each row, which starts with coefficients for a polynomial Π i . At this point the left side of the matrix has at most one 1 in each row. The second step is to zero out the bottom rows.
For the first step, work in descending total degree of the polynomials, starting at degree nd. Let the current row have the coefficients of i is the highest degree term have not changed yet, and therefore, one of the rows has the coefficients of i , which has decreased the total degree of the second term by one while keeping the set of variables the same. This is repeated until the row has the coefficients At the end of the first step, each row in the left side (i.e., columns indexed by nonmultilinear monomials) has at most one 1. This is in fact a constructive argument showing that Consider any two rows indexed bym( . Ifmx 2 i =m ′ x 2 j , they have the same set of variables, and therefore ψ(mx i ) = ψ(m ′ x j ), so the rows are equal and one can be eliminated (zeroed out). Keep doing this until for every leading nonmultilinear monomial there is only one row where the corresponding coefficient is nonzero.
Because for every column in the left part there is a unique nonzero row with the corresponding leading monomial, the matrix can be written (up to permutation of the rows) as Proof. By Lemma 5.2, using row operations on the F 2 submatrix ofM we get a matrix L 1 R 1 I 2 B 2 , where zero rows from F 2 are removed.
Row operations utilizing I 2 can then be used to zero out the top left, resulting in From the polynomial perspective, this maps all the nonmultilinear polynomials to their corresponding multilinear polynomials under ψ, i.e., for each monomial n i=1 x a i i , the map encodes the coefficient vector of ψ(( n i=1 x a i i )f j ), 1 ≤ j ≤ m into the Macaulay matrix as a row vector.
Recall that the Macaulay matrix has rows labeled by pairs; observe that rows (m, f i ) and (m ′ , f i ) will be equal when ψ(m) = ψ(m ′ ). In particular for any nonmultiliear monomialm the rows (m, f i ) and (ψ(m), f i ) will be equal at this point, so we can eliminate (zero out and then remove) any row indexed by nonmultiliear monomials. To be compatible with Definition 5.1 we choose not to further reduce / remove rows despite the factB might have zero rows, for example, rows with ψ(mf ) = ψ(m ′ f ′ ), but f ̸ = f ′ .
As in the general case letB = [M − ⃗ b] define the Boolean Macaulay linear system as M ⃗ y = ⃗ b, where the entries of ⃗ y are labeled by the nontrivial multilinear monomials and ⃗ b = 1 0 .
Recall that a matrix is s-sparse if it has at most s entries in any row or column.
Lemma 5.4. The Boolean Macaulay matrixB of total degree d of F 1 is an O (m · #F 1 )-sparse matrix.
Proof. The Boolean Macaulay matrixB is constructed by placing ψ(tf ) in a row for a multilinear monomial t and f ∈ F 1 , so the support of each row has size at most O (#F 1 ).
For the column sparsity first consider the Boolean Maculay matrix of {t}, which has a 1 matrix element at column m ′ and row (m ′′ , t) if and only if m ′ = ψ(m ′′ · t). This can only happen if t divides m ′ , so we can definet := m ′ /t. It is easy to see that m ′ = ψ(m ′′ · t) if and only if m ′′ =t · m d for some monomial m d that divides t. This implies that the column sparsity of the Boolean Maculay matrix of {t} equals the number of divisors of t which is at most 4 if the multilinear monomial t has (total) degree at most 2. Now consider the Boolean Maculay matrix of {f } for some (at most) quadratic polynomial f ∈ F 1 . Observe that f is a linear combination of at most #F 1 monomials of degree at most 2 and the Boolean Maculay matrix of {f } is likewise the linear combination of the Boolean Maculay matrix of these (at most) quadratic monomials. So the column sparsity of the Boolean Maculay matrix of {f } is at most 4 · #F 1 . Finally, the entire Boolean Macaulay matrix of F 1 is simply given by stacking the Boolean Maculay matrices of {f } for f ∈ F 1 , so the total column sparsity is at most 4m · #F 1 .
Note that this also implies that M , a submatrix ofB, is also sparse. Moreover, the location and value of the nonzero entries of each column/row ofB can be efficiently computed. Now we show that the Boolean Macaulay linear system is equivalent to the Macaulay linear system. It follows that solving the Boolean Macaulay linear system returns a correct solution of the Boolean polynomial system. 12 Lemma 5.5. Let M 1 ⃗ y 1 = ⃗ b 1 be the Macaulay linear system of a polynomial system F = F 1 ∪ F 2 and let M 2 ⃗ y 2 = ⃗ b 2 be the corresponding Boolean Therefore 0B Since performing row operations on the augmented matrix of a linear system does not change the set of solutions, solving the Macaulay linear system M 1 ⃗ y 1 = ⃗ b 1 is equivalent to solving the lin- , where the entries of ⃗ y 2 and ⃗ z 1 are indexed by nontrivial multilinear monomials and nonmultilinear monomials respectively. For the linear system we have M 2 ⃗ y 2 = ⃗ b 2 , which is the Boolean Macaulay linear system, and ⃗ z 1 + B ′ 2 ⃗ y 2 = 0. Ifŷ 2 is a solution of the Boolean Macaulay linear system M 2 ⃗ y 2 = ⃗ b 2 , setẑ 1 to be −B ′ 2ŷ 2 , then ẑ 1 y 2 is a solution of the linear system . Because the Macaulay linear system is equivalent to the linear system 0 M 2 I 2 B ′ 2 ⃗ z 1 ⃗ y 2 = ⃗ b 2 0 , therefore, a solutionŷ 2 of the Boolean Macaulay linear system M 2 ⃗ y 2 = ⃗ b 2 corresponds to a solutionŷ 1 of the Macaulay linear system M 1 ⃗ y 1 = ⃗ b 1 .
As M is a O (m · #F)-sparse row / column computable matrix and we can efficiently prepare the sparse vector ⃗ b as quantum state |b⟩, we can apply a QLS algorithm to "solve" the Boolean The key parameter in the running time is the condition number of the matrix M . Next we will provide a lower bound of the tQLScn of M and thus also a lower bound on known QLS algorithms.

Lower bound on the tQLScn κ ⃗ b (M )
Suppose a 1 , a 2 , · · · , a t ∈ {0, 1} n are the t solu- . , x 2 n − x n }, and let h be the minimum Hamming weight of the t solutions a 1 , a 2 , · · · , a t . Let ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t be the corresponding solution vectors of the Boolean Macaulay linear system M ⃗ y = ⃗ b under the assignments a 1 , a 2 , · · · , a t respectively.
In this case, we have ∥M ∥ ≥ 1/2 as M has at least one matrix element which has an absolute value at least 1/2 13 . Analogously to Theorem 4.5 we get: be the Boolean Macaulay matrix of F with columns labeled by multilinear monomials. Let h be the minimum Hamming weight of the t solutions a 1 , a 2 , · · · , a t . If all the t solutions a 1 , a 2 , · · · , a t have the same Hamming weight h or the minimum ℓ 2 -norm solution vector ⃗ y = M + ⃗ b is in the convex hull of ⃗ y 1 , ⃗ y 2 , · · · , ⃗ y t , then the tQLScn 13 After applying Red2 there is at least one polynomial f with a constant term of magnitude 1. If f does not have a degree-1 monomial xi, then xi · f has a magnitude 1 degree-1 monomial xi, so the row (xi, f ) will have a matrix element of magnitude 1. Otherwise, suppose the coefficient of x1 in f is c1, then the rows (1, f ) and (x1, f ) will have a matrix element of magnitude c1 and c1 − 1 respectively. Therefore, at least one of them has a magnitude of at least 1/2.
For h = Θ(log n), this lower bound does not rule out the possibility that the Macaulay matrix has a polynomial condition number, which would result in the quantum algorithm beating the brute-force classical algorithm that runs in time O n log n .

Details comparing running times
As we have discussed in Section 4.3 the classical brute-force algorithm tries all n j choices for the locations of the 1's in the solution a for each j ≤ h, and its running time can be bounded Comparing the above expression with our tQLScn lower bound, we saw that the Gorver-  = (a 1 , a 2 , · · · , a n ) ∈ {0, 1} n , then for some d less than or equal to n, the corresponding Boolean Macaulay linear system M ⃗ y = ⃗ b of total degree d has a unique solution ⃗ y = M + ⃗ b, where the entries of ⃗ y are indexed by multilinear monomials in x 1 , . . . , x n with total degree at most d. Let U = {x 1 , x 2 , . . . , x n }. There is a oneto-one correspondence between the subsets of U of size at most d and multilinear monomials in x 1 , . . . , x n with total degree at most d. Let S be the largest subset of U such that all the variables x k ∈ S have assignment a k = 1 and S d be the set containing all nonempty subsets of S that have size at most d. There is a one-to-one correspondence between the elements of the set S d and the nonzero entries of ⃗ y. Given implicit access to matrix M and sparse vector b, the QLS algorithm outputs the solution vector ⃗ y as the quantum state |⃗ y⟩, which encodes the nonzero entries of ⃗ y. Because the unique solution ⃗ y = M + ⃗ b of the Boolean Macaulay linear system is a 0/1 vector, the quantum state |⃗ y⟩ can be represented by If we measure the quantum state |⃗ y⟩, we will get a uniformly random subset R ∈ S d , where all the variables x k ∈ R have assignment a k = 1. Given copies of the quantum state |⃗ y⟩, the goal is to compute S. Next, we will reformulate this problem as a variant of the quantum coupon collector problem.
which is a superposition of subsets of S of size at most d. The goal is to compute S.
Specially, when d equals 1, this is the quantum coupon collector problem defined in [ABC + 20]. They proved that Θ(|S| log(min{|S|, n − |S|})) copies of the states |⃗ y⟩ are necessary to compute S.
Without loss of generality, we can assume d is at most |S| because when d is greater than |S|, the quantum state |⃗ y⟩ is the same as the case of d equals |S|. Then, we have the following result of Problem 6.1.

Theorem 6.2. Let r = O ((|S|/d) log(|S|/ε)).
Measuring r copies of the quantum superposition state in Problem 6.1, the set S can be computed with probability at least 1 − ε.
Since the only quantum operation is a measurement in the computational basis, this is essentially a classical coupon collector problem, where we can sample a uniformly random subset.
Proof. For any x ∈ S, the number of sets Hence, when 0 ≤ d ≤ ⌊ |S| 3 ⌋, the probability of not seeing x after r tries is at most Since , the when |S| = 3d + 2 the probability of seeing x is greater than d |S| · |S|−2d+1 |S|−d+1 ≥ 1 6 . Hence, the probability of not finding x after r tries is at most (1 − 1 6 ) r . Let r = O ((|S|/d) log(|S|/ε)), and by the union bound, the probability of not collecting all the elements x in S is at most ε. That is, if r = O ((|S|/d) log(|S|/ε)), the entire set S can be recovered with probability at least 1 − ε.

The algorithm
When a set of polynomials F has a unique solution, Algorithm 1 finds the solution. If a set of polynomials has more than one solution, we apply the Valiant-Vazirani reduction Red1 to get a set of polynomials F that have a unique solution.
Step 1: Apply a quantum linear system algorithm to the Boolean Macaulay linear system M ⃗ y = ⃗ b of total degree n and get the solution ⃗ y in quantum state Step 2: Perform measurement on the quantum state |⃗ y⟩ and get outcome |R⟩, then let all the variables in the set R equal 1.
Step 3: Repeat Step 1 and Step 2 O(log n) times, and then set all left remaining variables a j = 0.
Compared with Chen and Gao's [CG21] algorithm, there are two differences with Algorithm 1. First, the size of the Boolean Macaulay matrix in Algorithm 1 is m2 n × 2 n , which leads to a smaller lower bound of the tQLScn and leaves a possibility of superpolynomial speedup using Algorithm 1. By contrast, the size of the Macaulay matrix in Chen and Gao's algorithm is (m + n)(3n + 1) n × (3n + 1) n 15 , which leads to a larger lower bound of the tQLScn that prohibits a potential quantum speedup. Second, in Algorithm 1, the polynomial system have a unique solution, so in contrast to [CG21] the Boolean Macaulay linear system stays the same for every iteration and the number of iterations (measurements) required to obtain the solution of the polynomial system is O (log n). However, the Valiant-Vazirani reduction needs O (n) iterations to generate a polynomial system that has a unique solution with high probability. This amounts to O (n log n) iterations in total to find a solution.
On the other hand, in Chen and Gao's algorithm, the polynomial system could have any finite number of solutions, so the Macaulay linear system needs to be updated after each iteration (measurement) and the number of iterations is O (n).

Discussion
The Boolean Macaulay linear system approach is an interesting framework to study giving insights to the limitations and capabilities of quantum computation. On one hand, a lot of problems such as Factoring, Graph isomorphism, and Learning with binary errors can be put into this single framework. On the other hand, the QLS algorithm used for the (Boolean) Macaulay linear system is BQP-complete and the Factoring problem is known to be inside BQP by Shor's algorithm, therefore, if we can find an approach to get around the curse of the condition number of the Boolean Macaulay linear system derived from the Factoring problem, then we might be able to extend the result to other problems, such as Graph Isomorphism and Learning with binary errors, revealing new capabilities of quantum computation.
Our analytical lower bound on the condition number decreases when there are multiple solutions of the polynomial systems, but the polyno-mial systems used for cryptography purpose usually have one or few solutions [CG17], so our result gives strong evidence that the QLS algorithm cannot be used for attacking cryptosystems via the Macaulay matrix approach. Also, we suspect that having many solutions will not make the QLS algorithm to work substantially better. For example consider adding l new field equations y 2 i − y i = 0, then the number of solution of the new polynomial system will increase by a factor of 2 l , however the length of the shortest vector stays the same -indeed one can see that the shortest vector is an affine combination of solutions where all the new variables are set to 0. 16 Given an ill-conditioned QLSP, two main approaches have been proposed, one is the truncated QLS algorithm, the other one is the preconditioned QLS algorithm [HHL09]. Our lower bound on the tQLScn prohibits further speedup by the truncated QLS algorithm, however further investigation is needed regarding the possibility of using preconditioned QLS algorithms, such as parallel sparse approximate inverse preconditioner [CJS13], circulant preconditioner [SX18], fast inversion [TAWL20], or develop new preconditioned QLS algorithms for the Boolean Macaulay linear system. A promising feature of the Boolean Macaulay linear system is that it cannot be de-quantized by known classical techniques [CGL + 20] since the Boolean Macaulay linear system has high (full) column rank.
The introduced variant of the quantum coupon collector problem provides an example of how to extract the solution efficiently if the values in the solution vector of a QLSP are correlated in a nice pattern. Since applications of the QLS algorithm usually gain restricted access to the solution vector, a generalization of our extraction method, utilizing more general correlated patterns, could have interesting applications.
When the Hamming weight of the solution of a Boolean polynomial system is logarithmic in the number of variables, our lower bound does not rule out a superpolynomial speedup over exhaustive search. Finding a real application that exhibits such a superpolynomial speedup would be very interesting. However, for applications like polynomial systems over a finite field, the employed reductions usually increase the number of 16 To see this consider the coordinates corresponding to monomials not including the new variables. variables in the corresponding polynomial system significantly -likely compromising the potential speedup.
[AFI + 04] [Src22] The Mathematica source code is also available at the Woflram Notebook Archive: