The importance of zkSNARKs [1] is ever-increasing in the cryptocurrency and smart contracts ecosystem. Due to the significant threat of financial loss a bug represents in such applications, there is also a surge of interest in the formal verification of zkSNARK programs. These programs are ultimately encoded as a system of polynomial equations over large finite prime fields, and to prove statements about such a program is to prove statements about its system of equations. In this ongoing work, we investigate algebraic techniques with the goal of writing a mixed algebraic-SMT decision procedure to compute satisfiability of a new theory of polynomials over large prime fields. Ideally, the new theory and decision procedure could be implemented in existing SMT solvers as well as a standalone tool, in order to perform verification tasks over real world applications of zkSNARKs.
Zero-Knowledge Succinct Non-Interactive Argument of Knowledge (zkSNARKs) [
The need for secrecy over publicly available data makes blockchains a fertile ground for
zkSNARKs applications. This technology is fundamental, for example, for the ZCash [
r e t u r n ( y ∗ y ) == x
A zkSNARK program is often written in a high level DSL, such as ZoKrates [
The most complex zkSNARK applications have systems on the order of thousands of variables and millions of polynomial constraints, at which point the execution cost for generating a zkSNARK proofs becomes a limiting factor. Naturally, verifying properties of such systems is not any easier. Some verification tasks that are of particular interest to zkSNARKs application developers are: • given a constraint system, check whether it represents a function, i.e. whether the output variables are uniquely determined by the input; • given two constraint systems, check whether they are equisatisfiable; • more generically, check whether a constraint system conforms to a specification.
In this extended abstract, we present our work in progress trying to build a decision procedure that can be used to perform automated verification of such systems and their properties in practice, mixing algebraic and SMT approaches. We survey the existing solutions considering the challenges that our context enforces, and how they may be handled.
As an example, take the ZoKrates program in Figure 1, and consider there are two interacting persons: the verifier and the prover 3.
The verifier publishes a certain unsigned 8-bit integer , and is interested in whether someone
can prove that they know an unsigned 8-bit integer such that is the square root of . ZoKrates
compiles this code into the R1CS [
After the system is published by the verifier, the prover privately generates a witness, that is, an assignment for all the variables in the constraint system, including the private in our 2The exact value is 21888242871839275222246405745257275088548364400416034343698204186575808495617. 3The word “prover” here designates a party in the cryptosystem. Not to be confused with an automated reasoning engine. 1 . b1 ∗ b1 == b1 2 . b2 ∗ b2 == b2 3 . b3 ∗ b3 == b3 4 . b4 ∗ b4 == b4 5 . b5 ∗ b5 == b5 6 . b6 ∗ b6 == b6 7 . b7 ∗ b7 == b7 8 . b8 ∗ b8 == b8 9 . ( 1 2 8 ∗ b1 + 64 ∗ b2 + 32 ∗ b3 + 16 ∗ b4 +
8 ∗ b5 + 4 ∗ b6 + 2 ∗ b7 + 1 ∗ b8 ) == x example. Using the witness, the constraint system, and the proving key published by the verifier, the prover generates a zero knowledge proof that they know and that they have executed the circuit computations specified by the constraint system. Finally, the prover can submit the proof to a smart contract deployed by the verifier for an on-chain check that the proof is valid, without revealing any private data.
Besides the expected constraints for multiplication and equality, when compiling the ZoKrates high level code into circuit-like polynomial constraints the compiler also conveniently adds constraints that ensure that and are indeed 8-bit numbers, represented in Figure 2 for (analogous for ).
Since is a prime field element, the compiler decomposes as a binary number of 8 bits in constraint 9. However, the variables 1, . . . , 8 that represent those bits are also prime field elements. Thus, the compiler also adds constraints to ensure that each of those variables can only be 0 or 1.
Using current SMT solvers, we can represent the decomposition above as a predicate bits(, 1, 2, 3, 4, 5, 6, 7, 8) in the theory of Nonlinear Integer Arithmetic (NIA) applying (mod ) to both sides of each equality, where is the order of the large prime field employed by zkSNARK proof system, is the number to be decomposed, and if the predicate is true then 1, . . . , 8 are the bits that compose . Then we can ask whether this decomposition is deterministic via the assertion in Figure 3.
The full smtlib2 file can be found in [
) )
SMT solvers can often handle large software verification problems, but there is no theory
specifically written for large prime fields, and applying NIA solvers [
Ecne [
There are several algebraic approaches that can solve the systems of polynomials we propose in theory, however the complexity of the presented algorithms is prohibitively high to be used in practice. We discuss these works in detail in Sec. 3.
Given a system of polynomial statements, such as the following: 43 + 222 + 12 − = 5 3 − 32 ̸= 73 + 1 ̸= , our goal is to find if there is an assignment of the variables ( and in this case) in a field modulo a certain prime where all the statements are true.
It seems that inequations of the kind 2 ≤ do not make sense in a prime field because the values wrap around, so a total ordering among the elements is not quite meaningfully defined for most applications of a satisfiability solver. Alternatively, ranges of the kind ≤ 2 − < look much more useful, but they are not being investigated in the scope of this work.
The problem above can be reduced to the problem of finding the zeros of a set of polynomials, i.e. the value of x such that 1(x) = 2(x) = ... = (x) = 0, where x is the vector (1, 2, ..., ). For equalities this is simple, as one can just rearrange all terms to the lefthand side, leaving 0 on the right-hand side. The equality in the example above becomes the polynomial
1 = 43 + 222 + 12 − − 5.
For disequalities, we can use the “Rabinowitsch trick” [
In the above example, the disequalities produce the polynomial
2 = (3 − 32 − 73 − 1)( − ) − 1.
Thus, we ultimately obtain a set of polynomials over a finite field, and the problem is reduced to finding their common roots in that same field. Like the set of real numbers R, a finite field is not algebraically closed, and like 2 + 1 = 0 over R have roots and − in C (which is an extension field of R), a polynomial over a finite field might have roots that only exists in an extension field. We are not interested in such roots.
Every polynomial over a field has a unique factorization of the kind = ℎ1ℎ22ℎ33 . . . ℎ, where the ℎ are square-free (not divisible by any non-constant polynomial square) and pairwise coprime (have no common non-constant polynomial divisors). This is called square-free decomposition and is cheap to compute. It is easy to see that the polynomial ℎ1ℎ2ℎ3 . . . ℎ (the square-free component of ) has the same roots of .
It might prove beneficial to replace the polynomials with their square-free components, (and 12... with its square-free component, for the disequality polynomial). This is because the square-free decomposition may reduce the degree and size of the polynomial, thus simplifying the solving procedure. The solution would be identical, except for the introduced variable , which is not important to the problem.
Once we reduce our problem to a system of polynomial equations, we have to solve it. We
now describe the approaches considered, and what we have already implemented. In this
discussion, we assume the degrees of the input polynomials are small (not much larger than
10) and the prime field is large (around 2254). For a definition of the basic algebra concepts
used in this section (rings, ideals and irreducibility), check [
Before we try any expensive solving algorithm, let us consider some special cases that are easily handled. It is convenient to first try to reduce the set of polynomials against itself (reduction is a step of any Gröbner Basis algorithm, see Sec. 3.3), so that the system is sanitized (all trivial zeros and obviously equivalent polynomials are removed), and some easy to spot relations might be exposed, thus reducing the degree of some polynomials.
Special case #1: if the set of polynomials is empty, we consider it vacuously satisfiable and ifnish. This is consistent with the fact that the 0 polynomial is satisfiable, and that the initial reduction step removes all 0 polynomials from the set, so an empty set is equivalent to a set with only zero polynomials.
Special case #2: non-zero constant polynomial case: if the self-reduction returns a set with a non-zero constant polynomial (the reduction definition guarantees it will be the only polynomial left), it obviously cannot be zero, so the system is unsatisafible. In the current implementation of reduction, if a non-zero constant polynomial is ever encountered, the algorithm is cut short and returns immediately.
We then iterate through all the polynomials and count 3 things: the maximum degree of the system , the number of variables , and the total number of zero degree terms .
Special case #3: if = 0, there are no zero degree terms. Thus (0, 0, ..., 0) is a solution of the system, so we can stop.
Special case #4: if = 1 (it cannot be 0 because that was checked in special cases #1 and #2) then the system is linear. Reduction is the generalization of Gaussian elimination, so whatever remained of the system is necessarily solvable (otherwise it would have finished in #2) and already in row echelon form — from which it is trivial to compute the unique solution — and we can stop.
Special case #5: if = 1, reduction guarantees there will be only one polynomial left. Since
it is monic, we can extract a root by using the Cantor-Zassenhaus probabilistic algorithm [
Special case #6: this is not yet implemented, and we are not sure if it is worth it. In case of
a single polynomial with > 1, [
We now discuss the alternatives being explored for the general case.
If the set of polynomials has some particular properties, namely, it defines a prime ideal (an ideal
that is both primary and radical), and it defines an absolutely irreducible algebraic set (a variety),
then the algorithm in [
The problem is, of course, that this method is incomplete, as not all problems will obey these properties. It is unclear how many problems this technique would be able to solve on its own. It is also unclear what happens if the algorithm is executed on a polynomial set without those properties. Moreover, it is also unknown how to test whether an ideal is a radical without a full Gröbner Basis, and how to test whether an algebraic set is a variety. Nothing on this front has been implemented.
There is another class of algorithms based on algebraic geometry (original probabilistic version
given in [
This algorithm’s complexity is double exponential on the number of variables, and only works if > () , where is the order of the prime field, is the maximum degree of the polynomials and is the number of variables. Using our assumptions of ≈ 2254 and ≈ 10, this algorithm does not work in general for much more than 3 variables. It is unclear for how many instances of the problem this algorithm might work in practice, but at least for the probabilistic version, this algorithm does not look practical at all.
Nothing of the algebraic set decomposition algorithm has been implemented. From here, we
only implemented a test that tells whether a single multivariate polynomial that happens to be
absolutely irreducible has any roots, originally from [
Finding the reduced Gröbner Basis [
At the moment, we have implemented only a naive version of Buchberger’s algorithm, which is much worse in practice than Faugère’s F4 and F5 (but has the same worst case complexity). For a very small system it is capable of computing the Gröbner Basis in grevlex ordering, but we do not believe it would perform well for lexicographical ordering.
For this approach, we need the Gröbner Basis in lexicographical order, so that the system is in a
triangular format where we can eliminate one variable at a time. However Gröbner Basis is
much faster to compute on grevlex, so much that in practice lexicographical order is obtained
by transforming from grevlex via Gröbner Walk [
This yields a new set of polynomials with the same roots, but in a triangular format. Once the free variables are assigned, the system can be solved one variable at a time, so that the first polynomial has only one variable to be solved. This leaves the next polynomial with only one variable after replacing the variables already known, and so forth until the system is completely solved.
For instance, in 3 variable systems, if we define > > > 1 in our ordering (1 must always be the smallest), then the first variable that can be solved is , then , then , assuming the system is fully determined. If the system is not fully determined, and there are free variables, then the lowest variables are the ones to be assigned.
Having an algorithm to find the roots one variable at a time, we still have a branching problem, where each variable assignment leads to a new set of possible assignments for the next variable, and so on. Thus, the naive search expands into a typical NP problem.
Finding roots of a univariate polynomial in a prime field is relatively cheap, with the best known algorithm being factorizing the polynomial into unique factors (− 1)(− 2)...(− ) where are the roots of the polynomial. The most eficient algorithm for factorization in a prime field is called the Cantor-Zassenhaus algorithm, with complexity (3 log ) where is the degree of the polynomial and is the field’s characteristic prime.
Besides the Gröbner Basis in grevlex, we have not yet implemented the techniques mentioned above, Gröbner Walk (or some alternative) and Cantor-Zassenhaus being the most notable absences.
The unsolved problems of this approach include how to choose the variable ordering optimally and a search for variable assignments in the triangular system that is better than exhaustive trial and error. It is unclear how to perform a local search in the triangular system, much less how to prove unsatisfiability.
Given polynomials in the ring F[1, 2, ..., ], every polynomial set defines an algebraic set, which is the set of points in F where all polynomials evaluate to zero. This is the set from where we seek solutions: by finding one point in the algebraic set, we have proven the problem to be satisfiable, and by proving the algebraic set to be empty, we have proven the problem to be unsatisfiable.
The set of polynomials we are trying to solve defines an ideal, and by Lasker-Noether theorem we know that such polynomial ideal is the intersection of a finite set of primary ideals (called the primary decomposition of the ideal), which implies any algebraic set is a union of irreducible algebraic sets, named varieties.
So the idea of this approach is to find the radical of the ideal defined by our polynomial set,
compute its primary decomposition, and for each prime ideal found (which is primary and
radical), try to find a zero in its algebraic set. This seems to go back to the initial problem:
we have a set of polynomials for which we must find a common zero point. The diference is
that since we know we are dealing with an algebraic variety, we may be able to handle it with
algebraic geometry methods such as the one described in [
This can also be the first approach of a hybrid method, where the bigger problem is reduced to smaller sub-problems that can then be handled in parallel. Each subproblem is attempted to be solved first by the Fast Algebraic Geometry method, and if that fails, we can apply diferent approaches such as the triangular system method, local search, CDCL or others.
Computing the Gröbner Basis is the first step of the primary decomposition and radical ideal
algorithms, such as in [
The biggest issue of this method is that it is not complete: it decides for satisfiable if it is easy to extract a rational point from one of the algebraic sets (i.e. it is a variety) and it decides for unsatisfiable if the algebraic set is empty on the extension field (Gröbner Basis turns out to be a single non-zero constant). However, besides partitioning the problem, it does not help in case all the algebraic sets are reducible in the extension field.
The satisfiability problem can be solved by finding the Gröbner Basis of the system extended with the equations
− = 0 for each variable in the system (as per Fermat’s little theorem, = mod ). This removes all roots lying in an extension field but not in the base field.
For instance, for the system: 1 = 43 + 222 + 12 − − 5 2 = (3 − 32 − 73 − 1)( − ) − 1 the new polynomials would be: 3 = − 4 = − 5 = − where is the field’s characteristic prime.
The Gröbner Basis of the extended set contains a non-zero constant polynomial if the system is unsatisfiable. Such Gröbner Basis can be calculated using any monomial ordering, such as degree reverse lexicographical order, which is cheaper to compute than the lexicographical order.
It is unclear how to turn this into a practical algorithm for the case where ≫ 2, as the procedure to compute the Gröbner Basis is prohibitively expensive due to the time complexity being (2) on the reduction step of − = 0. We suspect it can be done in ( log ) by using exponentiation by squaring, but even linear complexity on is impractical, and would not help with the huge size that each reduced polynomial would have (on the order of ()).
This is more a complementary approach than an alternative, as it cannot prove unsatisfiability, and it may find a solution quickly if they are abundant.
A local search is the standard way of approximating floating point solutions for most equation systems found in engineering and simulations. The problem is treated as a black box function (called objective function), and the input is varied (the solution you want) according to how it afects the output, usually following the gradient towards a local minimum. A solution is found when zero is reached.
The traditional way of building an objective function for systems of equations is to take some norm (usually 2 or ∞) of the residual vector. In case of polynomials, the elements of the residual vector are the evaluation of each polynomial at the given input. In case of prime fields, it is important to consider that the distance to 0 is at most /2 when calculating the norm.
In SAT solving, discrete probabilistic methods (such as simulated annealing [
That said, the problem is much harder on the discrete case compared to continuous, because
the search might narrow down to oscillate around zero, but never reach it because it is in
between the minimum discrete steps. Such problem does not happen over reals. A metaheuristic
such as Tabu Search [
Some ideas that could be considered are: it might be beneficial to just use the equalities for the local search, removing the disequalities polynomial and considering each individual disequality as a constraint. If the disequalities happen to be useful in the search, it might be better to apply the Rabinowitsch trick individually for each one, because adding new variables is not a problem for this method, and this makes each restriction more explicit for the search algorithm.
The work in [
When a conflict is detected, the set of conflicting polynomial constraints is identified, and Cylindrical Algebraic Decomposition (CAD) is used to isolate a whole region in the search space where this conflict can happen. This region is then described as a set of polynomial constraints that are added to the problem as learned clauses.
The dificulty of generalizing this technique for prime fields is that we need a way of generalizing a conflict in one particular variable evaluation to a whole region of the search space (the role of CAD in the original algorithm). Moreover, we also need to be able to describe it in a meaningful way for the deduction system, i.e. define the new restrictions only in terms of = and ̸=, not with >, ≥ , < and ≤ .
What also complicates the matter is that we are only interested in solutions that happen to
be on the base field, not in the whole extension field, a restriction that [
The work in [
The problem of roots outside the base field is solved by including the polynomials − for
every variable , which is most likely the main cause for the poor performance displayed in
practice by the algorithms in [
We suspect it might be possible to improve the approach in [
In linear systems, the search can be reduced by looking at the bounds of the variables with the SIMPLEX algorithm. It may be that there is a SIMPLEX generalization that works here to speed up the search. However, that is unlikely since SIMPLEX works by iteratively restricting the convex search space limited by straight lines/hyperplanes, but the curves designated by non-linear polynomials are not planes, so the limited space will not be convex.
This possibility is entertained for the sake of completeness and has not been considered in practice.
An algorithm named Mutant Zhuang-Zi is given in [
In this work in progress, we investigate how to decide satisfiability of a system of polynomials over a large prime field. We are also currently implementing a tool that will hopefully be used for the formal verification of zkSNARKs constraint systems using the investigated techniques.
It seems that the bane of large prime fields are the field equations − = 0. When used either with Gröbner Basis or other algorithms, these equations makes the running time (and space) of the algorithms a superlinear function of the prime number, which is so large that even () is impractical.
The role of field equations in the existing algorithms is to restrict the roots found on the
prime field itself (i.e. filter out solutions found only on the extension field). Embedding them
explicitly on the system to be solved might only work for small prime fields, but with creative
indirect use, there are probabilistic algorithms such as Cantor-Zassenhaus and [
The authors thank Barry Whitehat, Mahmut Levent Doğan, Nikolaj Bjørner and Thibaut Schaefer for the valuable discussions and insights on Zero Knowledge Proofs, ZoKrates, polynomial constraint systems, Gröbner Basis and algebraic geometry.