ZKP — Deep into Bellman Library

Bellman is a zk-SNARK software library developed by the Zcash team in Rust language, which implements the Groth16 algorithm.
https://github.com/zcash/librustzcash/tree/master/bellman

The overall logic of Groth16 proof scheme is illustrated. And the data structures in bellman library and proof creation and verification logic are explained in detail.

1. Overall process

The overall process can be roughly divided into the following steps:
1. Flatten the problem polynomial and construct the corresponding circuit. This step is made by the upper-level application.
2. Generate R1CS (Rank 1 Constraint System) according to the circuit
3. Convert R1CS to QAP (Quadratic Arithmetic Program). The traditional method is to use Lagrange interpolation, but in order to reduce the computational complexity, it can be implemented by Fast Fourier Transform.
4. Setup parameters of the QAP problem, that is, CRS (Common Reference Strings).
5. Create proof based on CRS and input.
6. Verify proof.

The details of each step are described below in turn.

2. Setup Stage

The main work of the Setup stage is to generate CRS data. The related formulas are:

Note:

1. x in the formula corresponds to the variable tau in the code, and the corresponding x^i corresponds to powers_of_tau.

which is calculated by `powers_of_tau.z (& tau)` in the code.

2.1 Data Structure of Parameters

pub struct VerifyingKey<E: Engine> {
    pub alpha_g1: E::G1Affine,
    pub beta_g1: E::G1Affine,
    pub beta_g2: E::G2Affine,
    pub gamma_g2: E::G2Affine,
    pub delta_g1: E::G1Affine,
    pub delta_g2: E::G2Affine,
    pub ic: Vec<E::G1Affine>
}

where

pub struct Parameters<E: Engine> {
    pub vk: VerifyingKey<E>,
    pub h: Arc<Vec<E::G1Affine>>,
    pub l: Arc<Vec<E::G1Affine>>,
    pub a: Arc<Vec<E::G1Affine>>,
    pub b_g1: Arc<Vec<E::G1Affine>>,
    pub b_g2: Arc<Vec<E::G2Affine>>
}

where

The last three parameters a, b_g1, b_g2 do not appear in the formula. In fact, they are the values of the QAP polynomial calculated from x^i.

2.2 Variable

The Variable type represents each value in the input data, and is divided into public statement data and private witness data:

• Input type: statement data
• Aux type: witness data

pub enum Index {
    Input(usize),
    Aux(usize)
}
pub struct Variable(Index);

2.3 ConstraintSystem

ConstraintSystem is an interface that defines the following functions for generating different types of variables:

• one (): produces a variable of type Input with index 0
• alloc (): produces Aux type variables, the index is incremented
• alloc_input (): produces a variable of type Input, the index is incremented

For example:

let a = cs.alloc(...)
let b = cs.alloc(...)
let c = cs.alloc_input(...)
cs.enforce(
    || "a*b=c",
    |lc| lc + a,
    |lc| lc + b,
    |lc| lc + c
);

In the above example, c is a statement, and a and b are witnesses. The sample circuit is used to verify a * b = c.
If you want to verify a + b = c, you need to write the circuit like this:

cs.enforce(
    || "a*b=c",
    |lc| lc + a + b,
    |lc| lc + CS::one(),
    |lc| lc + c
);

2.4 Building R1CS

Circuit’s syntheticsize () will call enforce () of ConstraintSystem to build R1CS.
KeypairAssembly is an implementation of ConstraintSystem. The parameters of R1CS are stored in its member variables:

struct KeypairAssembly<E: Engine> {
    num_inputs: usize,
    num_aux: usize,
    num_constraints: usize,
    at_inputs: Vec<Vec<(E::Fr, usize)>>,
    bt_inputs: Vec<Vec<(E::Fr, usize)>>,
    ct_inputs: Vec<Vec<(E::Fr, usize)>>,
    at_aux: Vec<Vec<(E::Fr, usize)>>,
    bt_aux: Vec<Vec<(E::Fr, usize)>>,
    ct_aux: Vec<Vec<(E::Fr, usize)>>
}

Take the above a * b = c as an example:
num_inputs = 1
num_aux = 2
num_constraints = 1
The next six fields correspond to the A, B, and C matrices in R1CS:

2.5 Building QAP

The next step is to complete the conversion from R1CS to QAP. One step is to implement Lagrange interpolation using inverse discrete fast Fourier transform:

powers_of_tau.ifft(&worker);
let powers_of_tau = powers_of_tau.into_coeffs();

Here is the concept of a root of unity:
If ωⁿ = 1 , then omega is called the unit root. Taking the complex plane as an example, ωⁱ actually divides the unit circle into n parts.

Now we are under the condition of a finite cyclic group, according to Fermat’s Little Theorem: if p is a prime number and ω and p are coprime, then ω ᵖ⁻¹ mod p= 1 . Therefore, primes that are coprime with p can be used as unit roots.
At the same time, we found that any element can be represented by a linear combination of ω ⁻¹, ω⁻² , …, ω ᵖ⁻¹ , that is, they can be used as one of the cyclic groups. Group base. We use this set of bases to perform inverse discrete Fourier transform, we can get Lagrange interpolation coefficients, and convert R1CS to QAP.

2.6 Prepare verification parameters

Finally, the parameters with brackets in the following formula are needed to verify proof:

For quick verification, the values of these parameters are calculated in advance.

pub fn prepare_verifying_key<E: Engine>(
    vk: &VerifyingKey<E>
) -> PreparedVerifyingKey<E>
{
    let mut gamma = vk.gamma_g2;
    gamma.negate();
    let mut delta = vk.delta_g2;
    delta.negate();
PreparedVerifyingKey {
        alpha_g1_beta_g2: E::pairing(vk.alpha_g1, vk.beta_g2),
        neg_gamma_g2: gamma.prepare(),
        neg_delta_g2: delta.prepare(),
        ic: vk.ic.clone()
    }
}

A small detail: how to calculate division in a finite field?
Obviously, a⋅b = a⋅b⁻¹, how to calculate the multiplicative inverse b⁻¹?
According to Fermat’s little theorem: if p is prime and a and p are prime, then aᵖ⁻¹(mod p) = 1.
Therefore, in a finite field, b⁻¹ = bᵖ⁻².

fn inverse(&self) -> Option<Self> {
    if <Fr as Field>::is_zero(self) {
        None
    } else {
        Some(self.pow(&[(MODULUS_R.0 as u64) - 2]))
    }
}

3 Create Proof Stage

Groth16 algorithm generates Proof’s formula as follows:

Randomly select r and s and calculate the three values of [A]₁, [B]₂, [C]₁.

pub struct Proof<E: Engine> {
    pub a: E::G1Affine,
    pub b: E::G2Affine,
    pub c: E::G1Affine
}

The other parameters in the formula are known. The main difficulty is to calculate h (x), which will be described in detail below.

3.1 ​Calculate s⋅A(x), s⋅B(x), s⋅C(x)

Looking back at QAP, we need to transform R1CS into a series of polynomials through Lagrange interpolation:

And verifcation is needed for:

For that reason, s⋅A(x), s⋅B(x) and s⋅C(x) should be calculated. One method is to perform polynomial operations directly, and the other is to calculate their values at ω⁰, ω¹, ω²…, ωⁿ⁻¹, Polynomial coefficients are then obtained by iFFT (Inverse Fast Fourier Transform).
It is worth noting that when x takes ω⁰, ω¹, ω²…, ωⁿ⁻¹, A (x), B (x) and C (x) is actually the result of R1CS. Therefore, s⋅A (x) is actually the left input, s⋅B (x) is the right input, and s⋅C (x) is the output of the circuit. We can reuse Circuit’s syntheticsize () function for calculation.
When generating the proof, another implementation of ConstraintSystem, ProvingAssignment, is used to call Circuit’s syntheticsize () function.

struct ProvingAssignment<E: Engine> {
    ...
    a: Vec<Scalar<E>>,
    b: Vec<Scalar<E>>,
    c: Vec<Scalar<E>>,
input_assignment: Vec<E::Fr>,
    aux_assignment: Vec<E::Fr>
}

Basically similar to the Setup phase:

-alloc (): send witness data to aux_assignment
-alloc_input (): send statement data to input_assignment
-enforce () / eval (): Calculate s⋅A (x), s⋅B (x) and s⋅C (x), put them into 3 member variables a, b, and c, respectively in

Then the corresponding polynomial coefficients can be obtained by iFFT.

3.2 Calculate H(x)

With sA(x), sB(x) and sC(x), H(x) can be calculated as:

1. Define ω = σ⁽ʳ⁻¹⁾/ⁿ, then σⁿ= 1 (Fermat ’s Little Theorem)
2. First calculate the polynomial coefficients of a, b, and c by iFFT 3 times
3. Calculate the value of the polynomial on σω⁰, σω¹, σω²…, σωⁿ⁻¹.
4. After 3 times of FFT, get the a’, b’, c’
5. Since t (σωⁿⁱ) = (σωⁱ)ⁿ-1= σⁿ-1, calculate the point h (x) on the offset set

6. Do an iFFT once to calculate the polynomial coefficients with σ offset

7. According to the scaling theorem:

, divide the result obtained in the previous step by σ, and you get h (x) polynomial coefficients.

Based on the above analysis, a total of 3 FFTs and 4 iFFTs need to be performed.let’s check the according source code in detail:

First, get the calculation results from the previous section:

let mut a = EvaluationDomain::from_coeffs(prover.a)?;
let mut b = EvaluationDomain::from_coeffs(prover.b)?;
let mut c = EvaluationDomain::from_coeffs(prover.c)?;

Then, obtain the polynomial coefficients through iFFT, calculate the value on coset through FFT:

a.ifft(&worker);
a.coset_fft(&worker);
b.ifft(&worker);
b.coset_fft(&worker);
c.ifft(&worker);
c.coset_fft(&worker);

Next to calculate:

a.mul_assign(&worker, &b);
drop(b);
a.sub_assign(&worker, &c);
drop(c);
a.divide_by_z_on_coset(&worker);

Finally, get the polynomial coefficients by iFFT and divide by σ:

a.icoset_fft(&worker);

In addition, the next step is to calculate

as part of [C]₁ :

multiexp(&worker, params.get_h(a.len())?, FullDensity, a)

4. Verification Stage

The proof verification formula of Groth16 algorithm is as follows:

Check the following code for summation (using ic in CRS) to generate the intermediate variable acc:

let mut acc = pvk.ic[0].into_projective();
for (i, b) in public_inputs.iter().zip(pvk.ic.iter().skip(1)) {
        acc.add_assign(&b.mul(i.into_repr()));
    }

Then deform the formula slightly, and turn the problem into verifying the following equation:

The corresponding code is as follows:

Ok(E::final_exponentiation(
        &E::miller_loop([
            (&proof.a.prepare(), &proof.b.prepare()),
            (&acc.into_affine().prepare(), &pvk.neg_gamma_g2),
            (&proof.c.prepare(), &pvk.neg_delta_g2)
        ].into_iter())
    ).unwrap() == pvk.alpha_g1_beta_g2)

Now the basic process of the Bellman source code is analyzed.

Welcome to discuss Bellman source code and optimization with me :)