Understand how to optimize a circuit implementation of Groth16 verification

[daira]

This ticket will be a discussion of optimizations for implementing recursive SNARKs.

For Groth16 verification we need pairings, and for pairings we need extension field arithmetic. Assume that the base field arithmetic is efficient in the circuit — which in practice requires it to be the same as the R1CS field. Call that field F_p. We'll use C to refer to the cost of an R1CS constraint.

Table 1 of Beuchat, González-Díaz, Mitsunari, Okamoto, Rodríguez-Henríquez, and Teruya gives costs for operations in F_p², F_p⁶, and F_p¹², using Karatsuba multiplication and squaring. (This part, and in fact most of the pairing algorithm, doesn't differ between BN and BLS curves.)

Let (a, m, s, d), (a~, m~, s~, d~), and (A, M, S, D) denote the cost of field addition, multiplication, squaring, and division in F_p, F_p², and F_p⁶, respectively. This matches the paper except that we use "d" for division instead of "i" for inversion, since a division is its own operation which need not be decomposed into a multiplication and an inversion. In the circuit we have a = 0 and m = s = d = 1C.

Karatsuba is certainly optimal for multiplication and squaring in F_p², so that gives m~ = 3C and s~ = 2C. [Edit: I was thinking of Complex squaring and calling it Karatsuba here.] Division (or inversion) can be implemented via multiplication: d~ = 3C.

It isn't clear to me that Karatsuba is optimal for the extension of F_p² to F_p⁶. That gives 18C for multiplication, but I believe we can do it more efficiently by applying Toom-3 to get M = 5m~ = 15C. Note that the reduction from 5 F_p² terms to 3 terms is linear, and so free. On the other hand, that appears to beat Montgomery's 17C here, so maybe I made a mistake. (Major impostor syndrome here :-p )

[daira]

F_p⁶ squaring can be done by Toom-3 in 5s~ = 10C, which beats "Squaring Method 3" (1m~ + 4s~) in Chung and Hasan by an m~ - s~ tradeoff. Imposter syndrome intensifies! (Granted, they were considering the case where linear combinations are not free.)

Edit: actually Chung and Hasan do consider this variant; it's the first row of Table 2.

[Pratyush]

Just as a friendly heads up, Bls12's Fq is not highly 2-adic (q -1 is not a multiple of a large power of 2), and so using it in a SNARK is not so efficient (one would need to use different, more inefficient evaluation domains than the efficient powers-of-2 stuff that most SNARK impls use). For comparison, Bls12-381's Fr is highly 2-adic.

(We ran into the same issues in our project, and had to sample a new BLS curve that was highly 2-adic on both sides)

[Pratyush]

Also, there's a good comparison of the number of multiplications (and hence constraints) required to implement various finite-field extension arithmetic here: https://eprint.iacr.org/2006/471.pdf

[daira]

@Pratyush Yes, I knew that but thanks for clarifying. Do you mind giving the parameters of your curve?

The comments on F_p² and F_p⁶ arithmetic above are not specific to the curve, only (somewhat) to the extension degree. Thanks for the reference.

Edit: the best (circuit) algorithm for multiplication in F_p⁶ seems to be Toom-3 over Karatsuba, which costs 15m = 15C as given in table 9 of https://eprint.iacr.org/2006/471.pdf . So we have M = D = 15C and S = 10C.

(A = a~ = a = 0C, and reductions are also free.)

[daira]

The next step is curve arithmetic on E/F_p and E'/F_p². For short Weierstrass curves (e.g. any BN or BLS curve), it is most efficient to use affine coordinates; the formulae are the same as for Montgomery curves except for the numerator of λ in doubling, which does not affect the constraint count.

So, for E/F_p, incomplete addition costs 1s + 1m + 1d = 3C, and doubling costs 2s + 1m + 1d = 4C.

For E'/F_p², incomplete addition costs 1s~ + 1m~ + 1d~, and doubling costs 2s~ + 1m~ + 1d~.
With F_p² implemented using Karatsuba multiplication (m~ = d~ = 3m) and Complex squaring (s~ = 2m), that is (2 + 3 + 3)m = 8C for incomplete addition, and (4 + 3 + 3)m = 10C for doubling.

[Edit: use d~ for division instead of i~ for inversion, since this more accurately represents the circuit operation.]

[daira]

Next is the final exponentiation, for which we need arithmetic in F_p¹² and its cyclotomic subgroup G_Φ6(F_p²).

For multiplication in F_p¹², we can use Karatsuba over F_p⁶ which gives 3M = 45C. For squaring we can use Complex squaring over F_p⁶ which gives 2M = 30C.

(Those costs are for the tower F_p² → F_p⁶ → F_p¹². If we instead used F_p² → F_p⁴ → F_p¹² with Toom-3 over Karatsuba over Karatsuba, multiplication would require 15m~ = 45C and squaring 10m~ = 30C, so there would be no difference.)

For squaring in G_Φ6, we represent elements in the compressed representation (g₂, g₃, g₄, g₅) of Karabina. We can use the main algorithm in Karabina's paper which requires 4m~ = 12C, or the "SQR₂₃₄₅" one from that paper (with 2 squarings in F_p⁴) which is also 12C, or the algorithm in section 5.2 of Aranha, Karabina, Longa, Gebotys and López which requires 6s~ which again is 12C. So we can see that this saves significantly over the 30C needed to do squarings directly in F_p¹².

For multiplication in G_Φ6, we cannot use the compressed representation. In the case needed for the final exponentiation where the exponent is sparse, it is most efficient to compute all of the squarings in compressed format and then decompress the ones we need. A decompression uses the formulae:

    g₀ = (2.g₁² + g₂.g₅ - 3.g₃.g₄).ξ + 1
    g₁ = { (g₅².ξ + 3.g₄² - 2.g₃) / 4.g₂, if g₂ ≠ 0
            { (2.g₄.g₅) / g₃, if g₂ = 0

The formula for g₀ can be computed in terms of g₁ without any conditionals. We require 3C for a "boolean scaling" which computes a boolean b which is 0 when g₂ = 0, otherwise 1:

    (1 - b) × (b) = 0
    (g₂) × (λ) = (b)
    (μ) × (b) = (g₂)

and then 4C to select the numerator and denominator of the division (since each element of F_p² is 2 field elements). Note that computing two divisions and then paying 2C for a single selection is less efficient; trying to combine the boolean scaling with the selection (which would work if we only had a single selection) is also less efficient.

We then need 2s~ + m~ to compute g₅², g₄², and g₄.g₅, and d~ for the division. We need s~ + 2m~ to compute g₀. Overall this costs 3C + 4C + 2s~ + m~ + d~ + s~ + 2m~ = (3 + 4 + 4 + 3 + 3 + 2 + 6)C = 25C for each decompression.

So when u is sparse with k bits set and is ℓ bits in length, exponentiation by u (which is a component of the "hard part" of the final exponentiation) costs ℓ * 12C + (k-1) * (25 + 45)C = ℓ * 12C + (k-1) * 70C. The Karabina squaring and decompression pays off whenever (30 - 12)*ℓ > 25*(k-1), i.e. whenever k < 0.72*ℓ + 1, which is true for typical pairing-friendly curves.

We also need p-power and p²-power Frobenius operations, i.e. g ↦ g^p and g ↦ g^p2, for which we can use the formulae in section 3.2 of Beuchat, González-Díaz, Mitsunari, Okamoto, Rodríguez-Henríquez, and Teruya. These formulae only involve constant multiplications and conjugations so they are free in the circuit.

For the curve used in Aranha, Karabina, Longa, Gebotys and López:

The final exponentiation executes in total 1 inversion, 4 conjugations, 15 multiplications, 3 u-th powers, 4 cyclotomic squarings, 5 p-power Frobenius, 3 p²-power Frobenius.

So that is:

• 45C for the inversion in F_p¹²;
• the conjugations are free;
• 15 * 45C for the multiplications in F_p¹²;
• 3 * (62 * 12C + 2 * 70C) = 3 * 884C for the 3 u-th powers;
• 4 * 12C for the cyclotomic squarings;
• the Frobenius computations are free

giving a total of 3420C for the final exponentiation.

[Edit: corrected arithmetic error.]

[Edit: we can combine the inversion and one of the multiplications into a division in F_p¹², saving 45C. This gives a total of 3375C.]

[daira]

The remaining part of a pairing computation is the Miller loop. The basic operations in the body of the loop are "point doubling with line evaluation" and "point addition with line evaluation".

Doubling with line evaluation

A doubling on E'/F_p² requires 10C as computed above. Note that this already computes x₁². So the line evaluation (based on formula (2) of Aranha et al):

    l = (−2.y₁.y_P).v.w + (3.x₁².x_P).v² + ξ.(3b' - y₁²)

can be computed in an additional 2m~ + s~ = 8C, for a total of 18C. There are alternative ways of combining doubling and line evaluation, but I wasn't able to find any further improvement in the constraint cost.

TODO: check whether (x_P, y_P) is sometimes constant (this depends on how the pairing is used).

Addition with line evaluation

An incomplete addition on E'/F_p² costs 8C, also as computed above. Incomplete addition is sufficient provided that the input points P and Q are of large order. Aranha et al don't give formulae for the line evaluation used in the addition step, but their work is based on section 5 of Costello, Lange, and Naehrig (in turn based on section 4 of Costello, Hısıl, Boyd, Nieto, and Wong) which has very similar formulae for the doubling case. Switching to affine coordinates:

    l = (y₁ - y_P).(x_P - x_Q) + (x₁ - x_P).(y_Q - y_P)

or something like that. (This formula is probably slightly wrong because it is the version that eliminates the twisting of the point P, which is an optimization from Costello, Lange, and Naehrig that is handled differently in Aranha et al; but I don't think that affects the operation count.) So that is an additional 2m~, for a total of 14C.

Miller loop

Aranha et al section 6 says:

For the selected parameters and with the presented improvements, the Miller Loop in Algorithm 6 executes 64 point doublings with line evaluations, 6 point additions with line evaluations (4 inside Miller Loop and 2 more at the final steps), 1 negation in F_p² to precompute y_P, 1 p-power Frobenius, 1 p²-power Frobenius and 2 negations in E'(F_p²); and 1 conjugation, 1 multiplication, 66 sparse multiplications, 2 sparser multiplications and 63 squarings in F_p¹².

For now we won't bother to optimize the sparse[r] multiplications –partly because the paper doesn't describe that– so we will count those as full multiplications requiring 45C each. (This adds up to 3060C so is an obvious area for improvement.)

The resulting total cost of the rest of the pairing, i.e. everything but the final exponentiation, is 64 * 18C + 6 * 14C + (1+66+2) * 45C + 63 * 30C = 6231C.

[daira]

Grewal, Azarderakhsh, Longa, Hu, and Jao section 4.1 gives doubling/addition-with-line-evaluation formulae for affine coordinates on y² = x³ + b curves. In both cases they require only an extra 2m~ for the line evaluation, i.e.

• 2s~ + 1m~ + 1d~ + 2m~ = 16C for doubling with line evaluation (improving on 18C above);
• 1s~ + 1m~ + 1d~ + 2m~ = 14C for incomplete addition with line evaluation (same as above).

The sparse multiplication is described in "Pseudo 8-Sparse Multiplication for Efficient Ate-Based Pairing on Barreto-Naerig Curve", by Yuki Mori, Shoichi Akagi, Yasuyuki Nogani, and Masaaki Shirase, DOI: 10.1007/978-3-319-04873-4_11. This paper does not seem to be online, but the journal it is published in is available from Library Genesis. The relevant part is this:

This has 10 general multiplications (m~_u in the right-hand column), so the cost is 10m~ = 30C. We'll treat the "sparser" multiplications as just sparse for now.

So we can improve the cost excluding the final exponentiation to 64 * 16C + 6 * 14C + 45C + (66+2) * 30C + 63 * 30C = 5083C.

The total cost of a single pairing is (3375 + 5083)C = 8458C (less than a tenth of the cost of the Sapling Spend circuit).

[daira]

We may be able to improve the cost of the F_p¹² arithmetic by implementing F_p⁴ multiplication directly using Toom-4, and using that in an F_p⁴ → F_p¹² tower. Toom-4 multiplication costs 7m, instead of 3m~ = 9m for Karatsuba-over-Karatsuba. So Toom-3-over-Toom-4 costs 5 * 7m = 35C rather than 45C. Note that conversions between F_p⁴ → F_p¹² and F_p² → F_p⁶ → F_p¹² representations are free.

This helps most for the final exponentiation, where it saves 10C for the F_p¹² inversion and each of 21 F_p¹² multiplications (for the Aranha et al curve). Note that 6 of the multiplications are in the u-th power computations. There's also a single F_p¹² multiplication in the Miller loop (I'm conservatively assuming we can't apply Toom-4 to the sparse multiplications). So overall we save 230C giving a total cost of 8228C; around a 2.7% improvement.

[daira]

It seems that we can also just use Toom-12 directly to implement F_p¹² multiplication and squaring. This costs 23m for multiplication and 23s for squaring, i.e. 23C in both cases (this is better than the sparse multiplication algorithm; in this approach sparse multiplications have fewer terms but no fewer constraints).

It turns out that the cyclotomic squaring method still pays off whenever k < 0.44*ℓ + 1. For the Aranha et al curve, the final exponentiation takes (excluding free operations):

• 23C for the inversion in F_p¹²;
• 15 * 23C for the multiplications in F_p¹²;
• 3 * (62 * 12C + 2 * (25+23)C) = 3 * 840C for the 3 u-th powers;
• 4 * 12C for the cyclotomic squarings

giving a total of 2629C, a saving of 22.1%.

For the same curve, the Miller loop takes (again excluding free operations):

• 64 * 16C for the doublings with line evaluation;
• 6 * 14C for the incomplete additions with line evaluation;
• (1+66+2) * 23C for the multiplications in F_p¹²;
• 63 * 23C for the squarings in F_p¹²

giving a total of 4144C, a saving of 18.5%.

The overall cost of a pairing on the Aranha et al curve is (2629 + 4144)C = 6773C, which is a 19.9% improvement.

[daira]

From Appendix B.2 of the protocol spec, adapting notation to be consistent with above:

Neglecting multiplications in the r-order subgroup of F_p¹² and F_p, and other trivial operations, the cost of batched verification [for N proofs] is therefore

• for each proof: the cost of decoding the proof representation to the form Groth16.Proof, which requires three point decompressions and three prime-order subgroup checks (two for G₁ and one for G₂);
• for each successfully decoded proof: a Miller loop; and a 128-bit scalar multiplication by z_j;
• for each verification key: two Miller loops; an exponentiation in F_p¹²; a multiscalar multiplication with N 128-bit terms to compute Accum_∆; and a multiscalar multiplication with ℓ + 1 255-bit terms to compute ∑_i = 0..ℓ [Accum_Γ,i] Ψ_i;
• one final exponentiation.

A point decompression for Aranha-curve G₁ costs 291C for the range check + 3C for the on-curve check. (The Aranha BN curve has p = 36.u⁴ + 36.u³ + 24.u² + 6.u + 1 where u = -(2⁶² + 2⁵⁵ + 1). The comparison in the range check is with (p-1)/2 because we're using SORT compression.)

A point decompression for Aranha-curve G₂ costs 291C + 5C + 3C (where the 5C is for a trick to avoid having to do two comparisons, by selecting the element of the y polynomial to compare based on whether the most significant coefficient of the polynomial is 0).

The Aranha G₁ is prime-order so no subgroup checks are needed for π_A and π_C (this might be different for BLS and other curves). A subgroup check is needed for π_B. This can be done using the QED-it trick of witnessing π_B/h and then checking that [h] (π_B/h) = π_B ≠ 𝓞, where h is the cofactor in G₂. Since h is fixed, the best way to do this is to find a short addition-subtraction chain; unfortunately I don't know the order of E'/F_p² (the curve of which G₂ is a subgroup) so I can't compute h exactly. However we know that by the Hasse bound, h must be around lg(p² / r) ~= 254 bits, so we can implement the cofactor multiplication by h with roughly 253 doublings (10C each) and 84 incomplete additions or subtractions (8C each); probably better. That costs 253 * 10C + 84 * 8C = 3202C (quite expensive compared to the pairing, so worth optimizing further).

We will also need 128-bit variable-base scalar multiplications in G₁. Using "Montgomery without clamping", each of these costs (128 * 9 - 1)C = 1151C. (This includes conversion of π_j,A and π_j,C from Edwards to Montgomery.)

Let n = ceiling(lg(N)), so that Accum_Y has 128+n bits. (This is always less than r in practice, so reduction mod r is not useful.) The exponentiation Y^Accum_Y in F_p¹², using a simple square-and-conditional-multiply algorithm, costs (127+n) * 23C + (128+n) * (23 + 12)C = (7401 + 58 * n)C.

The computation of [Accum_Γ,i] Ψ_i is tricky. It's possible to do it by splitting each a_j,i into three 85-bit limbs. (One limb could be 84 bits but let's ignore that for simplicity.) Then we precompute bases [2^85.k] Ψ_i for k = 0..2, and do 3.(ℓ + 1) scalar multiplications each with a (128+85+n)-bit scalar. If we use the "fixed-base mostly Montgomery" approach, we only need a single adjustment at the end to correct for the digit offsets. If we assume for simplicity that n is a multiple of 3, then we have 71 + n/3 windows for each multiplication, so the cost per multiplication is (71 + n/3)*3 for the window lookups, (70 + n/3)*3 for the incomplete additions, and 4C to add the result of this multiplication to the partial result for [Accum_Γ,i] Ψ_i (with error checking for equal x-coordinates). So the cost of this part is 3 * (ℓ + 1) * (426 + n * 2 + 4)C = (ℓ + 1) * (1290 + n * 6)C. This is expensive so I'm sure there must be a better way; I'll see what I can come up with later.

Assume a single verification key, with ℓ field elements for the public input. In that case we need, for a batch of N Groth16 verifications on the Aranha curve:

• N * 294C for the G₁ decompression;
• N * 299C for the G₂ decompression;
• N * 3202C for the G₂ subgroup check;
• (N+2) * 4144C for N+2 Miller loops;
• (N+1) * 1151C for N+1 128-bit scalar multiplications by z_j in G₁;
• N * ℓ * 3 for computation of the three limbs of Accum_Γ,i = ∑_j = 0..N-1 z_j · a_j,i for each i ∈ {1..ℓ};
• (ℓ + 1) * (1290 + n * 6)C for the computation of [Accum_Γ,i] Ψ_i;
• (7401 + n * 58)C for an exponentiation in F_p¹²;
• 2629C for the final exponentiation;
• 3 * 23C for three multiplications in F_p¹²

for a total of (N * 9090 + 20828 + ℓ * (1290 + N * 3 + n * 6) + n * 64)C.

[Edit: the cost of the exponentiation was too low because it didn't take into account that the multiplies are conditional. The fix added an extra (128+n) * 12C.]

[Edit: fixed an arithmetic error in the cost of the F_p¹² exponentiation, adding (n * 23)C.]

[Edit: I had forgotten to include the final exponentiation. Oops!]