Jeffrey Scholz

# R1CS to Quadratic Arithmetic Program over a Finite Field in Python

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^2
```

We 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 1
```

We 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__.