To make the transformation from R1CS to QAP less abstract, let’s use a real example.
Let’s say we are encoding
out = x⁴ - 5y²x²
This will break down as
v1 = x * x
v2 = v1 * v1 # x^4
v3 = -5y * y
-v2 + out = v3*v1 # -5y^2 * x^2We need to pick a prime field we will do this over. When we later combine this with elliptic curves, the order of our prime field needs to equal the order of the elliptic curve. (This is a very common mistake).
But for now, we will pick a small number to make this manageable. We will pick the prime number 79.
First, we define our matrices R, L, and O
import numpy as np
# 1, out, x, y, v1, v2, v3
L = np.array([
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, -5, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1],
])
R = np.array([
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
])
O = np.array([
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, -1, 0],
])
To verify we constructed the R1CS correctly (it’s very easy to mess up when doing manually!) we create a valid witness and do the matrix multiplication
x = 4
y = -2
v1 = x * x
v2 = v1 * v1 # x^4
v3 = -5*y * y
out = v3*v1 + v2 # -5y^2 * x^2
witness = np.array([1, out, x, y, v1, v2, v3])
assert all(np.equal(np.matmul(L, witness) * np.matmul(R, witness), np.matmul(O, witness))), "not equal"Finite Field Arithmetic in Python
The next step is to convert this to a field array. Doing modular arithmetic in numpy will get very messy, but thankfully Python has libraries for this. We will use the python galois library.
Here is a quick overview of how the library works
import galois
GF = galois.GF(79)
a = GF(70)
b = GF(10)
print(a + b)
# prints 1We cannot give it negative values such as GF(-1) or it will throw an exception. For our purposes, it is okay to re-write the values manually, but we leave it as an exercise for the reader to come up with a solution that is more general.
Our new matrices are
L = np.array([
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 74, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1],
])
R = np.array([
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
])
O = np.array([
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 78, 0],
])
We can convert them to field arrays simply by wrapping them with GF now. We will also need to recompute our witness, because it contains negative values.
L_galois = GF(L)
R_galois = GF(R)
O_galois = GF(O)
x = GF(4)
y = GF(79-2) # we are using 79 as the field size, so 79 - 2 is -2
v1 = x * x
v2 = v1 * v1 # x^4
v3 = GF(79-5)*y * y
out = v3*v1 + v2 # -5y^2 * x^2
witness = GF(np.array([1, out, x, y, v1, v2, v3]))
assert all(np.equal(np.matmul(L_galois, witness) * np.matmul(R_galois, witness), np.matmul(O_galois, witness))), "not equal"Polynomial interpolation in finite fields
Now, we need to turn each of the columns of the matrices into a list of galois polynomials that interpolate the columns. The points we will interpolate are x = [1,2,3,4], since we have 4 rows.
def interpolate_column(col):
xs = GF(np.array([1,2,3,4]))
return galois.lagrange_poly(xs, col)
# axis 0 is the columns. apply_along_axis is the same as doing a for loop over the columns and collecting the results in an array
U_polys = np.apply_along_axis(interpolate_column, 0, L_galois)
V_polys = np.apply_along_axis(interpolate_column, 0, R_galois)
W_polys = np.apply_along_axis(interpolate_column, 0, O_galois)
If we look again at the contents of our matrices, we expect the first two polynomials of U and V to be zero, and the first column of W to be zero also.
L = np.array([
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 74, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1],
])
R = np.array([
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
])
O = np.array([
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 78, 0],
])We can sanity check this:
print(U_polys[:2])
print(V_polys[:2])
print(W_polys[:1])
# [Poly(0, GF(79)) Poly(0, GF(79))]# [Poly(0, GF(79)) Poly(0, GF(79))]# [Poly(0, GF(79))]The term Poly(0, GF(79)) is simply a polynomial where all the coefficients are zero.
The reader is encouraged to evaluate the polynomials at the values in the R1CS to see they interpolate the matrix values correctly.
Computing h(x)
We already know t(x) will be (x - 1)(x - 2)(x - 3)(x - 4) since there are four rows.
By way of reminder, this is the formula for a Quadratic Arithmetic Program. The variable a is the witness (1, a₁, …, aₘ).
Each of the terms is taking the inner product of the witness with the column-interpolating polynomials. That is, each of summation terms are effectively the inner product between <1, a₁, …, aₘ> and <u₀(x), u₁(x), ..., uₘ(x)>
def inner_product_polynomials_with_witness(polys, witness):
mul_ = lambda x, y: x * y
sum_ = lambda x, y: x + y
return reduce(sum_, map(mul_, polys, witness))
term_1 = inner_product_polynomials_with_witness(U_polys, witness)
term_2 = inner_product_polynomials_with_witness(V_polys, witness)
term_3 = inner_product_polynomials_with_witness(W_polys, witness)To compute h(x), we simply solve for it. Note that we cannot compute h(x) unless we have a valid witness, as we have combined the witness into our polynomials above.
# t = (x - 1)(x - 2)(x - 3)(x - 4)
t = galois.Poly([1, 78], field = GF) * galois.Poly([1, 77], field = GF) * galois.Poly([1, 76], field = GF) * galois.Poly([1, 75], field = GF)
h = (term_1 * term_2 - term_3) // t
Unlike poly1d from numpy, the galois library won’t indicate to us if there is a remainder, so we need to check if the QAP formula is still true.
assert term_1 * term_2 == term_3 + h * t, "division has a remainder"The scheme above will not work when we evaluate the polynomials on a hidden point from a trusted setup. However, the computer doing the trusted setup will still have to execute many of the computations above.
Learn more with RareSkills
This material is from our Zero Knowledge Course.