[MPC] Optimizations for linear combination

src/snarks/groth16/mpc_utils.hpp

@@ -57,6 +57,7 @@ template<typename ppT> class srs_mpc_layer_L1
 template<typename ppT>
 srs_mpc_layer_L1<ppT> mpc_compute_linearcombination(
     const srs_powersoftau &pot,
+    const srs_lagrange_evaluations &lagrange,
     const libsnark::qap_instance<libff::Fr<ppT>> &qap);
 
 /// Given the output from the first layer of the MPC, perform the 2nd

src/snarks/groth16/mpc_utils.tcc

@@ -5,11 +5,30 @@
 #include "mpc_utils.hpp"
 #include "multi_exp.hpp"
 
+#include <algorithm>
+#include <exception>
 #include <libff/algebra/scalar_multiplication/multiexp.hpp>
+#include <libfqfft/evaluation_domain/domains/basic_radix2_domain_aux.tcc>
 
 namespace libzeth
 {
 
+template<typename T>
+static void fill_vector_from_map(
+    std::vector<T> &out_vector,
+    const std::map<size_t, T> index_map,
+    const size_t index_bound)
+{
+    out_vector.resize(index_bound);
+    for (auto &it : index_map) {
+        const size_t out_idx = it.first;
+        const T &value = it.second;
+        if (!value.is_zero()) {
+            out_vector[out_idx] = value;
+        }
+    }
+}
+
 template<typename ppT>
 srs_mpc_layer_L1<ppT>::srs_mpc_layer_L1(
     libff::G1_vector<ppT> &&T_tau_powers_g1,
@@ -28,13 +47,13 @@ srs_mpc_layer_L1<ppT>::srs_mpc_layer_L1(
 template<typename ppT>
 srs_mpc_layer_L1<ppT> mpc_compute_linearcombination(
     const srs_powersoftau &pot,
+    const srs_lagrange_evaluations &lagrange,
     const libsnark::qap_instance<libff::Fr<ppT>> &qap)
 {
     using Fr = libff::Fr<ppT>;
     using G1 = libff::G1<ppT>;
     using G2 = libff::G2<ppT>;
-
-    libfqfft::evaluation_domain<Fr> &domain = *qap.domain;
+    libff::enter_block("Call to mpc_compute_linearcombination");
 
     // n = number of constraints in r1cs, or equivalently, n = deg(t(x))
     // t(x) being the target polynomial of the QAP
@@ -44,87 +63,112 @@ srs_mpc_layer_L1<ppT> mpc_compute_linearcombination(
     const size_t n = qap.degree();
     const size_t num_variables = qap.num_variables();
 
+    if (n != 1ull << libff::log2(n)) {
+        throw std::invalid_argument("non-pow-2 domain");
+    }
+    if (n != lagrange.degree) {
+        throw std::invalid_argument(
+            "domain size differs from Lagrange evaluation");
+    }
+
+    libff::print_indent();
+    printf("n=%zu\n", n);
+
     // The QAP polynomials A, B, C are of degree (n-1) as we know they
     // are created by interpolation of an r1cs of n constraints.
     // As a consequence, the polynomial (A.B - C) is of degree 2n-2,
     // while the target polynomial t is of degree n.
     // Thus, we need to have access (in the SRS) to powers up to 2n-2.
     // To represent such polynomials we need {x^i} for in {0, ... n-2}
-    // hence why we check below that we have at least n-1 elements
+    // hence we check below that we have at least n-1 elements
     // in the set of powers of tau
     assert(pot.tau_powers_g1.size() >= 2 * n - 1);
 
-    // n+1 coefficients of t
-    std::vector<Fr> t_coefficients(n + 1, Fr::zero());
-    qap.domain->add_poly_Z(Fr::one(), t_coefficients);
-
-    // Compute [ t(x) . x^i ]_1 for i = 0 .. n-2
-    libff::G1_vector<ppT> t_x_pow_i(n - 1);
+    // Domain uses n-roots of unity, so
+    //      t(x)       = x^n - 1
+    //  =>  t(x) . x^i = x^(n+i) - x^i
+    libff::G1_vector<ppT> t_x_pow_i(n - 1, G1::zero());
+    libff::enter_block("computing [t(x) . x^i]_1");
     for (size_t i = 0; i < n - 1; ++i) {
-        // Use { [x^i] , ... , [x^(i+n)] } with coefficients of t() to
-        // compute t(x).x^i.
-        t_x_pow_i[i] = multi_exp<ppT, G1>(
-            pot.tau_powers_g1.begin() + i,
-            pot.tau_powers_g1.begin() + i + n + 1,
-            t_coefficients.begin(),
-            t_coefficients.end());
+        t_x_pow_i[i] = pot.tau_powers_g1[n + i] - pot.tau_powers_g1[i];
     }
-
-    // Compute [ beta.A_i(x) + alpha.B_i(x) + C_i(x) ]_1
-    //
-    // For each i, get the Lagrange factors of A_i, B_i, C_i.  For
-    // each j, if A, B or C has a non-zero factor, grab the Lagrange
-    // coefficients, evaluate at [t]_1, multiply by the factor and
-    // accumulate.
-    libff::G1_vector<ppT> A_i_g1(num_variables + 1);
-    libff::G1_vector<ppT> B_i_g1(num_variables + 1);
-    libff::G2_vector<ppT> B_i_g2(num_variables + 1);
-    libff::G1_vector<ppT> ABC_i_g1(num_variables + 1);
-
-    evaluator_from_lagrange<ppT, G1> tau_eval_g1(pot.tau_powers_g1, domain);
-    evaluator_from_lagrange<ppT, G2> tau_eval_g2(pot.tau_powers_g2, domain);
-    evaluator_from_lagrange<ppT, G1> alpha_tau_eval(
-        pot.alpha_tau_powers_g1, domain);
-    evaluator_from_lagrange<ppT, G1> beta_tau_eval(
-        pot.beta_tau_powers_g1, domain);
-
-    for (size_t i = 0; i < num_variables + 1; ++i) {
-        // Compute [beta.A_i(x)], [alpha.B_i(x)] . [C_i(x)]
-        const std::map<size_t, Fr> &A_i_in_lagrange =
-            qap.A_in_Lagrange_basis[i];
-        const std::map<size_t, Fr> &B_i_in_lagrange =
-            qap.B_in_Lagrange_basis[i];
-        const std::map<size_t, Fr> &C_i_in_lagrange =
-            qap.C_in_Lagrange_basis[i];
-
-        A_i_g1[i] = tau_eval_g1.evaluate_from_lagrange_factors(A_i_in_lagrange);
-        B_i_g1[i] = tau_eval_g1.evaluate_from_lagrange_factors(B_i_in_lagrange);
-        B_i_g2[i] = tau_eval_g2.evaluate_from_lagrange_factors(B_i_in_lagrange);
-
-        G1 beta_A_at_t =
-            beta_tau_eval.evaluate_from_lagrange_factors(A_i_in_lagrange);
-        G1 alpha_B_at_t =
-            alpha_tau_eval.evaluate_from_lagrange_factors(B_i_in_lagrange);
-        G1 C_at_t = tau_eval_g1.evaluate_from_lagrange_factors(C_i_in_lagrange);
-
-        ABC_i_g1[i] = beta_A_at_t + alpha_B_at_t + C_at_t;
+    libff::leave_block("computing [t(x) . x^i]_1");
+
+    libff::enter_block("computing A_i, B_i, C_i, ABC_i at x");
+    libff::G1_vector<ppT> As_g1(num_variables + 1);
+    libff::G1_vector<ppT> Bs_g1(num_variables + 1);
+    libff::G2_vector<ppT> Bs_g2(num_variables + 1);
+    libff::G1_vector<ppT> Cs_g1(num_variables + 1);
+    libff::G1_vector<ppT> ABCs_g1(num_variables + 1);
+#ifdef MULTICORE
+#pragma omp parallel for
+#endif
+    for (size_t j = 0; j < num_variables + 1; ++j) {
+        G1 ABC_j_at_x = G1::zero();
+
+        {
+            G1 A_j_at_x = G1::zero();
+            const std::map<size_t, Fr> &A_j_lagrange =
+                qap.A_in_Lagrange_basis[j];
+            for (const auto &entry : A_j_lagrange) {
+                A_j_at_x = A_j_at_x +
+                           (entry.second * lagrange.lagrange_g1[entry.first]);
+                ABC_j_at_x =
+                    ABC_j_at_x +
+                    (entry.second * lagrange.beta_lagrange_g1[entry.first]);
+            }
+
+            As_g1[j] = A_j_at_x;
+        }
+
+        {
+            G1 B_j_at_x_g1 = G1::zero();
+            G2 B_j_at_x_g2 = G2::zero();
+
+            const std::map<size_t, Fr> &B_j_lagrange =
+                qap.B_in_Lagrange_basis[j];
+            for (const auto &entry : B_j_lagrange) {
+                B_j_at_x_g1 = B_j_at_x_g1 + (entry.second *
+                                             lagrange.lagrange_g1[entry.first]);
+                B_j_at_x_g2 = B_j_at_x_g2 + (entry.second *
+                                             lagrange.lagrange_g2[entry.first]);
+                ABC_j_at_x =
+                    ABC_j_at_x +
+                    (entry.second * lagrange.alpha_lagrange_g1[entry.first]);
+            }
+
+            Bs_g1[j] = B_j_at_x_g1;
+            Bs_g2[j] = B_j_at_x_g2;
+        }
+
+        {
+            G1 C_j_at_x = G1::zero();
+            const std::map<size_t, Fr> &C_j_lagrange =
+                qap.C_in_Lagrange_basis[j];
+            for (const auto &entry : C_j_lagrange) {
+                C_j_at_x =
+                    C_j_at_x + entry.second * lagrange.lagrange_g1[entry.first];
+            }
+
+            Cs_g1[j] = C_j_at_x;
+            ABC_j_at_x = ABC_j_at_x + C_j_at_x;
+        }
+
+        ABCs_g1[j] = ABC_j_at_x;
     }
-
-    assert(num_variables + 1 == A_i_g1.size());
-    assert(num_variables + 1 == B_i_g1.size());
-    assert(num_variables + 1 == B_i_g2.size());
-    assert(num_variables + 1 == ABC_i_g1.size());
+    libff::leave_block("computing A_i, B_i, C_i, ABC_i at x");
 
     // TODO: Consider dropping those entries we know will not be used
     // by this circuit and using sparse vectors where it makes sense
     // (as is done for B_i's in r1cs_gg_ppzksnark_proving_key).
+    libff::leave_block("Call to mpc_compute_linearcombination");
 
     return srs_mpc_layer_L1<ppT>(
         std::move(t_x_pow_i),
-        std::move(A_i_g1),
-        std::move(B_i_g1),
-        std::move(B_i_g2),
-        std::move(ABC_i_g1));
+        std::move(As_g1),
+        std::move(Bs_g1),
+        std::move(Bs_g2),
+        std::move(ABCs_g1));
 }
 
 template<typename ppT>