Fix libsnark dependency build.

This changes libsnark to build in-place, instead of copying first to
a build directory. Previously, modifications made to the original
sources wouldn't get rebuilt without a 'make clean' because users
would be pointing to the copies.

This closes #2689.

src/snark/libsnark/relations/arithmetic_programs/qap/qap.hpp

@@ -0,0 +1,193 @@
+/** @file
+ *****************************************************************************
+
+ Declaration of interfaces for a QAP ("Quadratic Arithmetic Program").
+
+ QAPs are defined in \[GGPR13].
+
+ References:
+
+ \[GGPR13]:
+ "Quadratic span programs and succinct NIZKs without PCPs",
+ Rosario Gennaro, Craig Gentry, Bryan Parno, Mariana Raykova,
+ EUROCRYPT 2013,
+ <http://eprint.iacr.org/2012/215>
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef QAP_HPP_
+#define QAP_HPP_
+
+#include "algebra/evaluation_domain/evaluation_domain.hpp"
+
+namespace libsnark {
+
+/* forward declaration */
+template<typename FieldT>
+class qap_witness;
+
+/**
+ * A QAP instance.
+ *
+ * Specifically, the datastructure stores:
+ * - a choice of domain (corresponding to a certain subset of the field);
+ * - the number of variables, the degree, and the number of inputs; and
+ * - coefficients of the A,B,C polynomials in the Lagrange basis.
+ *
+ * There is no need to store the Z polynomial because it is uniquely
+ * determined by the domain (as Z is its vanishing polynomial).
+ */
+template<typename FieldT>
+class qap_instance {
+private:
+    size_t num_variables_;
+    size_t degree_;
+    size_t num_inputs_;
+
+public:
+    std::shared_ptr<evaluation_domain<FieldT> > domain;
+
+    std::vector<std::map<size_t, FieldT> > A_in_Lagrange_basis;
+    std::vector<std::map<size_t, FieldT> > B_in_Lagrange_basis;
+    std::vector<std::map<size_t, FieldT> > C_in_Lagrange_basis;
+
+    qap_instance(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                 const size_t num_variables,
+                 const size_t degree,
+                 const size_t num_inputs,
+                 const std::vector<std::map<size_t, FieldT> > &A_in_Lagrange_basis,
+                 const std::vector<std::map<size_t, FieldT> > &B_in_Lagrange_basis,
+                 const std::vector<std::map<size_t, FieldT> > &C_in_Lagrange_basis);
+
+    qap_instance(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                 const size_t num_variables,
+                 const size_t degree,
+                 const size_t num_inputs,
+                 std::vector<std::map<size_t, FieldT> > &&A_in_Lagrange_basis,
+                 std::vector<std::map<size_t, FieldT> > &&B_in_Lagrange_basis,
+                 std::vector<std::map<size_t, FieldT> > &&C_in_Lagrange_basis);
+
+    qap_instance(const qap_instance<FieldT> &other) = default;
+    qap_instance(qap_instance<FieldT> &&other) = default;
+    qap_instance& operator=(const qap_instance<FieldT> &other) = default;
+    qap_instance& operator=(qap_instance<FieldT> &&other) = default;
+
+    size_t num_variables() const;
+    size_t degree() const;
+    size_t num_inputs() const;
+
+    bool is_satisfied(const qap_witness<FieldT> &witness) const;
+};
+
+/**
+ * A QAP instance evaluation is a QAP instance that is evaluated at a field element t.
+ *
+ * Specifically, the datastructure stores:
+ * - a choice of domain (corresponding to a certain subset of the field);
+ * - the number of variables, the degree, and the number of inputs;
+ * - a field element t;
+ * - evaluations of the A,B,C (and Z) polynomials at t;
+ * - evaluations of all monomials of t;
+ * - counts about how many of the above evaluations are in fact non-zero.
+ */
+template<typename FieldT>
+class qap_instance_evaluation {
+private:
+    size_t num_variables_;
+    size_t degree_;
+    size_t num_inputs_;
+public:
+    std::shared_ptr<evaluation_domain<FieldT> > domain;
+
+    FieldT t;
+
+    std::vector<FieldT> At, Bt, Ct, Ht;
+
+    FieldT Zt;
+
+    qap_instance_evaluation(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                            const size_t num_variables,
+                            const size_t degree,
+                            const size_t num_inputs,
+                            const FieldT &t,
+                            const std::vector<FieldT> &At,
+                            const std::vector<FieldT> &Bt,
+                            const std::vector<FieldT> &Ct,
+                            const std::vector<FieldT> &Ht,
+                            const FieldT &Zt);
+    qap_instance_evaluation(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                            const size_t num_variables,
+                            const size_t degree,
+                            const size_t num_inputs,
+                            const FieldT &t,
+                            std::vector<FieldT> &&At,
+                            std::vector<FieldT> &&Bt,
+                            std::vector<FieldT> &&Ct,
+                            std::vector<FieldT> &&Ht,
+                            const FieldT &Zt);
+
+    qap_instance_evaluation(const qap_instance_evaluation<FieldT> &other) = default;
+    qap_instance_evaluation(qap_instance_evaluation<FieldT> &&other) = default;
+    qap_instance_evaluation& operator=(const qap_instance_evaluation<FieldT> &other) = default;
+    qap_instance_evaluation& operator=(qap_instance_evaluation<FieldT> &&other) = default;
+
+    size_t num_variables() const;
+    size_t degree() const;
+    size_t num_inputs() const;
+
+    bool is_satisfied(const qap_witness<FieldT> &witness) const;
+};
+
+/**
+ * A QAP witness.
+ */
+template<typename FieldT>
+class qap_witness {
+private:
+    size_t num_variables_;
+    size_t degree_;
+    size_t num_inputs_;
+
+public:
+    FieldT d1, d2, d3;
+
+    std::vector<FieldT> coefficients_for_ABCs;
+    std::vector<FieldT> coefficients_for_H;
+
+    qap_witness(const size_t num_variables,
+                const size_t degree,
+                const size_t num_inputs,
+                const FieldT &d1,
+                const FieldT &d2,
+                const FieldT &d3,
+                const std::vector<FieldT> &coefficients_for_ABCs,
+                const std::vector<FieldT> &coefficients_for_H);
+
+    qap_witness(const size_t num_variables,
+                const size_t degree,
+                const size_t num_inputs,
+                const FieldT &d1,
+                const FieldT &d2,
+                const FieldT &d3,
+                const std::vector<FieldT> &coefficients_for_ABCs,
+                std::vector<FieldT> &&coefficients_for_H);
+
+    qap_witness(const qap_witness<FieldT> &other) = default;
+    qap_witness(qap_witness<FieldT> &&other) = default;
+    qap_witness& operator=(const qap_witness<FieldT> &other) = default;
+    qap_witness& operator=(qap_witness<FieldT> &&other) = default;
+
+    size_t num_variables() const;
+    size_t degree() const;
+    size_t num_inputs() const;
+};
+
+} // libsnark
+
+#include "relations/arithmetic_programs/qap/qap.tcc"
+
+#endif // QAP_HPP_

src/snark/libsnark/relations/arithmetic_programs/qap/qap.tcc

@@ -0,0 +1,324 @@
+/** @file
+*****************************************************************************
+
+Implementation of interfaces for a QAP ("Quadratic Arithmetic Program").
+
+See qap.hpp .
+
+*****************************************************************************
+* @author     This file is part of libsnark, developed by SCIPR Lab
+*             and contributors (see AUTHORS).
+* @copyright  MIT license (see LICENSE file)
+*****************************************************************************/
+
+#ifndef QAP_TCC_
+#define QAP_TCC_
+
+#include "common/profiling.hpp"
+#include "common/utils.hpp"
+#include "algebra/evaluation_domain/evaluation_domain.hpp"
+#include "algebra/scalar_multiplication/multiexp.hpp"
+
+namespace libsnark {
+
+template<typename FieldT>
+qap_instance<FieldT>::qap_instance(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                                   const size_t num_variables,
+                                   const size_t degree,
+                                   const size_t num_inputs,
+                                   const std::vector<std::map<size_t, FieldT> > &A_in_Lagrange_basis,
+                                   const std::vector<std::map<size_t, FieldT> > &B_in_Lagrange_basis,
+                                   const std::vector<std::map<size_t, FieldT> > &C_in_Lagrange_basis) :
+    num_variables_(num_variables),
+    degree_(degree),
+    num_inputs_(num_inputs),
+    domain(domain),
+    A_in_Lagrange_basis(A_in_Lagrange_basis),
+    B_in_Lagrange_basis(B_in_Lagrange_basis),
+    C_in_Lagrange_basis(C_in_Lagrange_basis)
+{
+}
+
+template<typename FieldT>
+qap_instance<FieldT>::qap_instance(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                                   const size_t num_variables,
+                                   const size_t degree,
+                                   const size_t num_inputs,
+                                   std::vector<std::map<size_t, FieldT> > &&A_in_Lagrange_basis,
+                                   std::vector<std::map<size_t, FieldT> > &&B_in_Lagrange_basis,
+                                   std::vector<std::map<size_t, FieldT> > &&C_in_Lagrange_basis) :
+    num_variables_(num_variables),
+    degree_(degree),
+    num_inputs_(num_inputs),
+    domain(domain),
+    A_in_Lagrange_basis(std::move(A_in_Lagrange_basis)),
+    B_in_Lagrange_basis(std::move(B_in_Lagrange_basis)),
+    C_in_Lagrange_basis(std::move(C_in_Lagrange_basis))
+{
+}
+
+template<typename FieldT>
+size_t qap_instance<FieldT>::num_variables() const
+{
+    return num_variables_;
+}
+
+template<typename FieldT>
+size_t qap_instance<FieldT>::degree() const
+{
+    return degree_;
+}
+
+template<typename FieldT>
+size_t qap_instance<FieldT>::num_inputs() const
+{
+    return num_inputs_;
+}
+
+template<typename FieldT>
+bool qap_instance<FieldT>::is_satisfied(const qap_witness<FieldT> &witness) const
+{
+    const FieldT t = FieldT::random_element();
+
+    std::vector<FieldT> At(this->num_variables()+1, FieldT::zero());
+    std::vector<FieldT> Bt(this->num_variables()+1, FieldT::zero());
+    std::vector<FieldT> Ct(this->num_variables()+1, FieldT::zero());
+    std::vector<FieldT> Ht(this->degree()+1);
+
+    const FieldT Zt = this->domain->compute_Z(t);
+
+    const std::vector<FieldT> u = this->domain->lagrange_coeffs(t);
+
+    for (size_t i = 0; i < this->num_variables()+1; ++i)
+    {
+        for (auto &el : A_in_Lagrange_basis[i])
+        {
+            At[i] += u[el.first] * el.second;
+        }
+
+        for (auto &el : B_in_Lagrange_basis[i])
+        {
+            Bt[i] += u[el.first] * el.second;
+        }
+
+        for (auto &el : C_in_Lagrange_basis[i])
+        {
+            Ct[i] += u[el.first] * el.second;
+        }
+    }
+
+    FieldT ti = FieldT::one();
+    for (size_t i = 0; i < this->degree()+1; ++i)
+    {
+        Ht[i] = ti;
+        ti *= t;
+    }
+
+    const qap_instance_evaluation<FieldT> eval_qap_inst(this->domain,
+                                                        this->num_variables(),
+                                                        this->degree(),
+                                                        this->num_inputs(),
+                                                        t,
+                                                        std::move(At),
+                                                        std::move(Bt),
+                                                        std::move(Ct),
+                                                        std::move(Ht),
+                                                        Zt);
+    return eval_qap_inst.is_satisfied(witness);
+}
+
+template<typename FieldT>
+qap_instance_evaluation<FieldT>::qap_instance_evaluation(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                                                         const size_t num_variables,
+                                                         const size_t degree,
+                                                         const size_t num_inputs,
+                                                         const FieldT &t,
+                                                         const std::vector<FieldT> &At,
+                                                         const std::vector<FieldT> &Bt,
+                                                         const std::vector<FieldT> &Ct,
+                                                         const std::vector<FieldT> &Ht,
+                                                         const FieldT &Zt) :
+    num_variables_(num_variables),
+    degree_(degree),
+    num_inputs_(num_inputs),
+    domain(domain),
+    t(t),
+    At(At),
+    Bt(Bt),
+    Ct(Ct),
+    Ht(Ht),
+    Zt(Zt)
+{
+}
+
+template<typename FieldT>
+qap_instance_evaluation<FieldT>::qap_instance_evaluation(const std::shared_ptr<evaluation_domain<FieldT> > &domain,
+                                                         const size_t num_variables,
+                                                         const size_t degree,
+                                                         const size_t num_inputs,
+                                                         const FieldT &t,
+                                                         std::vector<FieldT> &&At,
+                                                         std::vector<FieldT> &&Bt,
+                                                         std::vector<FieldT> &&Ct,
+                                                         std::vector<FieldT> &&Ht,
+                                                         const FieldT &Zt) :
+    num_variables_(num_variables),
+    degree_(degree),
+    num_inputs_(num_inputs),
+    domain(domain),
+    t(t),
+    At(std::move(At)),
+    Bt(std::move(Bt)),
+    Ct(std::move(Ct)),
+    Ht(std::move(Ht)),
+    Zt(Zt)
+{
+}
+
+template<typename FieldT>
+size_t qap_instance_evaluation<FieldT>::num_variables() const
+{
+    return num_variables_;
+}
+
+template<typename FieldT>
+size_t qap_instance_evaluation<FieldT>::degree() const
+{
+    return degree_;
+}
+
+template<typename FieldT>
+size_t qap_instance_evaluation<FieldT>::num_inputs() const
+{
+    return num_inputs_;
+}
+
+template<typename FieldT>
+bool qap_instance_evaluation<FieldT>::is_satisfied(const qap_witness<FieldT> &witness) const
+{
+
+    if (this->num_variables() != witness.num_variables())
+    {
+        return false;
+    }
+
+    if (this->degree() != witness.degree())
+    {
+        return false;
+    }
+
+    if (this->num_inputs() != witness.num_inputs())
+    {
+        return false;
+    }
+
+    if (this->num_variables() != witness.coefficients_for_ABCs.size())
+    {
+        return false;
+    }
+
+    if (this->degree()+1 != witness.coefficients_for_H.size())
+    {
+        return false;
+    }
+
+    if (this->At.size() != this->num_variables()+1 || this->Bt.size() != this->num_variables()+1 || this->Ct.size() != this->num_variables()+1)
+    {
+        return false;
+    }
+
+    if (this->Ht.size() != this->degree()+1)
+    {
+        return false;
+    }
+
+    if (this->Zt != this->domain->compute_Z(this->t))
+    {
+        return false;
+    }
+
+    FieldT ans_A = this->At[0] + witness.d1*this->Zt;
+    FieldT ans_B = this->Bt[0] + witness.d2*this->Zt;
+    FieldT ans_C = this->Ct[0] + witness.d3*this->Zt;
+    FieldT ans_H = FieldT::zero();
+
+    ans_A = ans_A + naive_plain_exp<FieldT, FieldT>(this->At.begin()+1, this->At.begin()+1+this->num_variables(),
+                                                    witness.coefficients_for_ABCs.begin(), witness.coefficients_for_ABCs.begin()+this->num_variables());
+    ans_B = ans_B + naive_plain_exp<FieldT, FieldT>(this->Bt.begin()+1, this->Bt.begin()+1+this->num_variables(),
+                                                    witness.coefficients_for_ABCs.begin(), witness.coefficients_for_ABCs.begin()+this->num_variables());
+    ans_C = ans_C + naive_plain_exp<FieldT, FieldT>(this->Ct.begin()+1, this->Ct.begin()+1+this->num_variables(),
+                                                    witness.coefficients_for_ABCs.begin(), witness.coefficients_for_ABCs.begin()+this->num_variables());
+    ans_H = ans_H + naive_plain_exp<FieldT, FieldT>(this->Ht.begin(), this->Ht.begin()+this->degree()+1,
+                                                    witness.coefficients_for_H.begin(), witness.coefficients_for_H.begin()+this->degree()+1);
+
+    if (ans_A * ans_B - ans_C != ans_H * this->Zt)
+    {
+        return false;
+    }
+
+    return true;
+}
+
+template<typename FieldT>
+qap_witness<FieldT>::qap_witness(const size_t num_variables,
+                                 const size_t degree,
+                                 const size_t num_inputs,
+                                 const FieldT &d1,
+                                 const FieldT &d2,
+                                 const FieldT &d3,
+                                 const std::vector<FieldT> &coefficients_for_ABCs,
+                                 const std::vector<FieldT> &coefficients_for_H) :
+    num_variables_(num_variables),
+    degree_(degree),
+    num_inputs_(num_inputs),
+    d1(d1),
+    d2(d2),
+    d3(d3),
+    coefficients_for_ABCs(coefficients_for_ABCs),
+    coefficients_for_H(coefficients_for_H)
+{
+}
+
+template<typename FieldT>
+qap_witness<FieldT>::qap_witness(const size_t num_variables,
+                                 const size_t degree,
+                                 const size_t num_inputs,
+                                 const FieldT &d1,
+                                 const FieldT &d2,
+                                 const FieldT &d3,
+                                 const std::vector<FieldT> &coefficients_for_ABCs,
+                                 std::vector<FieldT> &&coefficients_for_H) :
+    num_variables_(num_variables),
+    degree_(degree),
+    num_inputs_(num_inputs),
+    d1(d1),
+    d2(d2),
+    d3(d3),
+    coefficients_for_ABCs(coefficients_for_ABCs),
+    coefficients_for_H(std::move(coefficients_for_H))
+{
+}
+
+
+template<typename FieldT>
+size_t qap_witness<FieldT>::num_variables() const
+{
+    return num_variables_;
+}
+
+template<typename FieldT>
+size_t qap_witness<FieldT>::degree() const
+{
+    return degree_;
+}
+
+template<typename FieldT>
+size_t qap_witness<FieldT>::num_inputs() const
+{
+    return num_inputs_;
+}
+
+
+} // libsnark
+
+#endif // QAP_TCC_

src/snark/libsnark/relations/arithmetic_programs/qap/tests/test_qap.cpp

@@ -0,0 +1,104 @@
+/**
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstring>
+#include <vector>
+
+#include "algebra/curves/alt_bn128/alt_bn128_pp.hpp"
+#include "algebra/fields/field_utils.hpp"
+#include "common/profiling.hpp"
+#include "common/utils.hpp"
+#include "reductions/r1cs_to_qap/r1cs_to_qap.hpp"
+#include "relations/constraint_satisfaction_problems/r1cs/examples/r1cs_examples.hpp"
+
+#include <gtest/gtest.h>
+
+using namespace libsnark;
+
+template<typename FieldT>
+void test_qap(const size_t qap_degree, const size_t num_inputs, const bool binary_input)
+{
+    /*
+      We construct an instance where the QAP degree is qap_degree.
+      So we generate an instance of R1CS where the number of constraints qap_degree - num_inputs - 1.
+      See the transformation from R1CS to QAP for why this is the case.
+      So we need that qap_degree >= num_inputs + 1.
+    */
+    ASSERT_LE(num_inputs + 1, qap_degree);
+    enter_block("Call to test_qap");
+
+    const size_t num_constraints = qap_degree - num_inputs - 1;
+
+    print_indent(); printf("* QAP degree: %zu\n", qap_degree);
+    print_indent(); printf("* Number of inputs: %zu\n", num_inputs);
+    print_indent(); printf("* Number of R1CS constraints: %zu\n", num_constraints);
+    print_indent(); printf("* Input type: %s\n", binary_input ? "binary" : "field");
+
+    enter_block("Generate constraint system and assignment");
+    r1cs_example<FieldT> example;
+    if (binary_input)
+    {
+        example = generate_r1cs_example_with_binary_input<FieldT>(num_constraints, num_inputs);
+    }
+    else
+    {
+        example = generate_r1cs_example_with_field_input<FieldT>(num_constraints, num_inputs);
+    }
+    leave_block("Generate constraint system and assignment");
+
+    enter_block("Check satisfiability of constraint system");
+    EXPECT_TRUE(example.constraint_system.is_satisfied(example.primary_input, example.auxiliary_input));
+    leave_block("Check satisfiability of constraint system");
+
+    const FieldT t = FieldT::random_element(),
+    d1 = FieldT::random_element(),
+    d2 = FieldT::random_element(),
+    d3 = FieldT::random_element();
+
+    enter_block("Compute QAP instance 1");
+    qap_instance<FieldT> qap_inst_1 = r1cs_to_qap_instance_map(example.constraint_system);
+    leave_block("Compute QAP instance 1");
+
+    enter_block("Compute QAP instance 2");
+    qap_instance_evaluation<FieldT> qap_inst_2 = r1cs_to_qap_instance_map_with_evaluation(example.constraint_system, t);
+    leave_block("Compute QAP instance 2");
+
+    enter_block("Compute QAP witness");
+    qap_witness<FieldT> qap_wit = r1cs_to_qap_witness_map(example.constraint_system, example.primary_input, example.auxiliary_input, d1, d2, d3);
+    leave_block("Compute QAP witness");
+
+    enter_block("Check satisfiability of QAP instance 1");
+    EXPECT_TRUE(qap_inst_1.is_satisfied(qap_wit));
+    leave_block("Check satisfiability of QAP instance 1");
+
+    enter_block("Check satisfiability of QAP instance 2");
+    EXPECT_TRUE(qap_inst_2.is_satisfied(qap_wit));
+    leave_block("Check satisfiability of QAP instance 2");
+
+    leave_block("Call to test_qap");
+}
+
+TEST(relations, qap)
+{
+    start_profiling();
+
+    const size_t num_inputs = 10;
+
+    enter_block("Test QAP with binary input");
+
+    test_qap<Fr<alt_bn128_pp> >(1ul << 21, num_inputs, true);
+
+    leave_block("Test QAP with binary input");
+
+    enter_block("Test QAP with field input");
+
+    test_qap<Fr<alt_bn128_pp> >(1ul << 21, num_inputs, false);
+
+    leave_block("Test QAP with field input");
+}

src/snark/libsnark/relations/constraint_satisfaction_problems/r1cs/examples/r1cs_examples.hpp

@@ -0,0 +1,73 @@
+/** @file
+ *****************************************************************************
+
+ Declaration of interfaces for a R1CS example, as well as functions to sample
+ R1CS examples with prescribed parameters (according to some distribution).
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef R1CS_EXAMPLES_HPP_
+#define R1CS_EXAMPLES_HPP_
+
+#include "relations/constraint_satisfaction_problems/r1cs/r1cs.hpp"
+
+namespace libsnark {
+
+/**
+ * A R1CS example comprises a R1CS constraint system, R1CS input, and R1CS witness.
+ */
+template<typename FieldT>
+struct r1cs_example {
+    r1cs_constraint_system<FieldT> constraint_system;
+    r1cs_primary_input<FieldT> primary_input;
+    r1cs_auxiliary_input<FieldT> auxiliary_input;
+
+    r1cs_example<FieldT>() = default;
+    r1cs_example<FieldT>(const r1cs_example<FieldT> &other) = default;
+    r1cs_example<FieldT>(const r1cs_constraint_system<FieldT> &constraint_system,
+                         const r1cs_primary_input<FieldT> &primary_input,
+                         const r1cs_auxiliary_input<FieldT> &auxiliary_input) :
+        constraint_system(constraint_system),
+        primary_input(primary_input),
+        auxiliary_input(auxiliary_input)
+    {};
+    r1cs_example<FieldT>(r1cs_constraint_system<FieldT> &&constraint_system,
+                         r1cs_primary_input<FieldT> &&primary_input,
+                         r1cs_auxiliary_input<FieldT> &&auxiliary_input) :
+        constraint_system(std::move(constraint_system)),
+        primary_input(std::move(primary_input)),
+        auxiliary_input(std::move(auxiliary_input))
+    {};
+};
+
+/**
+ * Generate a R1CS example such that:
+ * - the number of constraints of the R1CS constraint system is num_constraints;
+ * - the number of variables of the R1CS constraint system is (approximately) num_constraints;
+ * - the number of inputs of the R1CS constraint system is num_inputs;
+ * - the R1CS input consists of ``full'' field elements (typically require the whole log|Field| bits to represent).
+ */
+template<typename FieldT>
+r1cs_example<FieldT> generate_r1cs_example_with_field_input(const size_t num_constraints,
+                                                            const size_t num_inputs);
+
+/**
+ * Generate a R1CS example such that:
+ * - the number of constraints of the R1CS constraint system is num_constraints;
+ * - the number of variables of the R1CS constraint system is (approximately) num_constraints;
+ * - the number of inputs of the R1CS constraint system is num_inputs;
+ * - the R1CS input consists of binary values (as opposed to ``full'' field elements).
+ */
+template<typename FieldT>
+r1cs_example<FieldT> generate_r1cs_example_with_binary_input(const size_t num_constraints,
+                                                             const size_t num_inputs);
+
+} // libsnark
+
+#include "relations/constraint_satisfaction_problems/r1cs/examples/r1cs_examples.tcc"
+
+#endif // R1CS_EXAMPLES_HPP_

src/snark/libsnark/relations/constraint_satisfaction_problems/r1cs/examples/r1cs_examples.tcc

@@ -0,0 +1,164 @@
+/** @file
+ *****************************************************************************
+
+ Implementation of functions to sample R1CS examples with prescribed parameters
+ (according to some distribution).
+
+ See r1cs_examples.hpp .
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef R1CS_EXAMPLES_TCC_
+#define R1CS_EXAMPLES_TCC_
+
+#include <cassert>
+
+#include "common/utils.hpp"
+
+namespace libsnark {
+
+template<typename FieldT>
+r1cs_example<FieldT> generate_r1cs_example_with_field_input(const size_t num_constraints,
+                                                            const size_t num_inputs)
+{
+    enter_block("Call to generate_r1cs_example_with_field_input");
+
+    assert(num_inputs <= num_constraints + 2);
+
+    r1cs_constraint_system<FieldT> cs;
+    cs.primary_input_size = num_inputs;
+    cs.auxiliary_input_size = 2 + num_constraints - num_inputs; // TODO: explain this
+
+    r1cs_variable_assignment<FieldT> full_variable_assignment;
+    FieldT a = FieldT::random_element();
+    FieldT b = FieldT::random_element();
+    full_variable_assignment.push_back(a);
+    full_variable_assignment.push_back(b);
+
+    for (size_t i = 0; i < num_constraints-1; ++i)
+    {
+        linear_combination<FieldT> A, B, C;
+
+        if (i % 2)
+        {
+            // a * b = c
+            A.add_term(i+1, 1);
+            B.add_term(i+2, 1);
+            C.add_term(i+3, 1);
+            FieldT tmp = a*b;
+            full_variable_assignment.push_back(tmp);
+            a = b; b = tmp;
+        }
+        else
+        {
+            // a + b = c
+            B.add_term(0, 1);
+            A.add_term(i+1, 1);
+            A.add_term(i+2, 1);
+            C.add_term(i+3, 1);
+            FieldT tmp = a+b;
+            full_variable_assignment.push_back(tmp);
+            a = b; b = tmp;
+        }
+
+        cs.add_constraint(r1cs_constraint<FieldT>(A, B, C));
+    }
+
+    linear_combination<FieldT> A, B, C;
+    FieldT fin = FieldT::zero();
+    for (size_t i = 1; i < cs.num_variables(); ++i)
+    {
+        A.add_term(i, 1);
+        B.add_term(i, 1);
+        fin = fin + full_variable_assignment[i-1];
+    }
+    C.add_term(cs.num_variables(), 1);
+    cs.add_constraint(r1cs_constraint<FieldT>(A, B, C));
+    full_variable_assignment.push_back(fin.squared());
+
+    /* split variable assignment */
+    r1cs_primary_input<FieldT> primary_input(full_variable_assignment.begin(), full_variable_assignment.begin() + num_inputs);
+    r1cs_primary_input<FieldT> auxiliary_input(full_variable_assignment.begin() + num_inputs, full_variable_assignment.end());
+
+    /* sanity checks */
+    assert(cs.num_variables() == full_variable_assignment.size());
+    assert(cs.num_variables() >= num_inputs);
+    assert(cs.num_inputs() == num_inputs);
+    assert(cs.num_constraints() == num_constraints);
+    assert(cs.is_satisfied(primary_input, auxiliary_input));
+
+    leave_block("Call to generate_r1cs_example_with_field_input");
+
+    return r1cs_example<FieldT>(std::move(cs), std::move(primary_input), std::move(auxiliary_input));
+}
+
+template<typename FieldT>
+r1cs_example<FieldT> generate_r1cs_example_with_binary_input(const size_t num_constraints,
+                                                             const size_t num_inputs)
+{
+    enter_block("Call to generate_r1cs_example_with_binary_input");
+
+    assert(num_inputs >= 1);
+
+    r1cs_constraint_system<FieldT> cs;
+    cs.primary_input_size = num_inputs;
+    cs.auxiliary_input_size = num_constraints; /* we will add one auxiliary variable per constraint */
+
+    r1cs_variable_assignment<FieldT> full_variable_assignment;
+    for (size_t i = 0; i < num_inputs; ++i)
+    {
+        full_variable_assignment.push_back(FieldT(std::rand() % 2));
+    }
+
+    size_t lastvar = num_inputs-1;
+    for (size_t i = 0; i < num_constraints; ++i)
+    {
+        ++lastvar;
+        const size_t u = (i == 0 ? std::rand() % num_inputs : std::rand() % i);
+        const size_t v = (i == 0 ? std::rand() % num_inputs : std::rand() % i);
+
+        /* chose two random bits and XOR them together:
+           res = u + v - 2 * u * v
+           2 * u * v = u + v - res
+        */
+        linear_combination<FieldT> A, B, C;
+        A.add_term(u+1, 2);
+        B.add_term(v+1, 1);
+        if (u == v)
+        {
+            C.add_term(u+1, 2);
+        }
+        else
+        {
+            C.add_term(u+1, 1);
+            C.add_term(v+1, 1);
+        }
+        C.add_term(lastvar+1, -FieldT::one());
+
+        cs.add_constraint(r1cs_constraint<FieldT>(A, B, C));
+        full_variable_assignment.push_back(full_variable_assignment[u] + full_variable_assignment[v] - full_variable_assignment[u] * full_variable_assignment[v] - full_variable_assignment[u] * full_variable_assignment[v]);
+    }
+
+    /* split variable assignment */
+    r1cs_primary_input<FieldT> primary_input(full_variable_assignment.begin(), full_variable_assignment.begin() + num_inputs);
+    r1cs_primary_input<FieldT> auxiliary_input(full_variable_assignment.begin() + num_inputs, full_variable_assignment.end());
+
+    /* sanity checks */
+    assert(cs.num_variables() == full_variable_assignment.size());
+    assert(cs.num_variables() >= num_inputs);
+    assert(cs.num_inputs() == num_inputs);
+    assert(cs.num_constraints() == num_constraints);
+    assert(cs.is_satisfied(primary_input, auxiliary_input));
+
+    leave_block("Call to generate_r1cs_example_with_binary_input");
+
+    return r1cs_example<FieldT>(std::move(cs), std::move(primary_input), std::move(auxiliary_input));
+}
+
+} // libsnark
+
+#endif // R1CS_EXAMPLES_TCC

src/snark/libsnark/relations/constraint_satisfaction_problems/r1cs/r1cs.hpp

@@ -0,0 +1,153 @@
+/** @file
+ *****************************************************************************
+
+ Declaration of interfaces for:
+ - a R1CS constraint,
+ - a R1CS variable assignment, and
+ - a R1CS constraint system.
+
+ Above, R1CS stands for "Rank-1 Constraint System".
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef R1CS_HPP_
+#define R1CS_HPP_
+
+#include <cstdlib>
+#include <iostream>
+#include <map>
+#include <string>
+#include <vector>
+
+#include "relations/variable.hpp"
+
+namespace libsnark {
+
+/************************* R1CS constraint ***********************************/
+
+template<typename FieldT>
+class r1cs_constraint;
+
+template<typename FieldT>
+std::ostream& operator<<(std::ostream &out, const r1cs_constraint<FieldT> &c);
+
+template<typename FieldT>
+std::istream& operator>>(std::istream &in, r1cs_constraint<FieldT> &c);
+
+/**
+ * A R1CS constraint is a formal expression of the form
+ *
+ *                < A , X > * < B , X > = < C , X > ,
+ *
+ * where X = (x_0,x_1,...,x_m) is a vector of formal variables and A,B,C each
+ * consist of 1+m elements in <FieldT>.
+ *
+ * A R1CS constraint is used to construct a R1CS constraint system (see below).
+ */
+template<typename FieldT>
+class r1cs_constraint {
+public:
+
+    linear_combination<FieldT> a, b, c;
+
+    r1cs_constraint() {};
+    r1cs_constraint(const linear_combination<FieldT> &a,
+                    const linear_combination<FieldT> &b,
+                    const linear_combination<FieldT> &c);
+
+    r1cs_constraint(const std::initializer_list<linear_combination<FieldT> > &A,
+                    const std::initializer_list<linear_combination<FieldT> > &B,
+                    const std::initializer_list<linear_combination<FieldT> > &C);
+
+    bool operator==(const r1cs_constraint<FieldT> &other) const;
+
+    friend std::ostream& operator<< <FieldT>(std::ostream &out, const r1cs_constraint<FieldT> &c);
+    friend std::istream& operator>> <FieldT>(std::istream &in, r1cs_constraint<FieldT> &c);
+};
+
+/************************* R1CS variable assignment **************************/
+
+/**
+ * A R1CS variable assignment is a vector of <FieldT> elements that represents
+ * a candidate solution to a R1CS constraint system (see below).
+ */
+
+/* TODO: specify that it does *NOT* include the constant 1 */
+template<typename FieldT>
+using r1cs_primary_input = std::vector<FieldT>;
+
+template<typename FieldT>
+using r1cs_auxiliary_input = std::vector<FieldT>;
+
+template<typename FieldT>
+using r1cs_variable_assignment = std::vector<FieldT>; /* note the changed name! (TODO: remove this comment after primary_input transition is complete) */
+
+/************************* R1CS constraint system ****************************/
+
+template<typename FieldT>
+class r1cs_constraint_system;
+
+template<typename FieldT>
+std::ostream& operator<<(std::ostream &out, const r1cs_constraint_system<FieldT> &cs);
+
+template<typename FieldT>
+std::istream& operator>>(std::istream &in, r1cs_constraint_system<FieldT> &cs);
+
+/**
+ * A system of R1CS constraints looks like
+ *
+ *     { < A_k , X > * < B_k , X > = < C_k , X > }_{k=1}^{n}  .
+ *
+ * In other words, the system is satisfied if and only if there exist a
+ * USCS variable assignment for which each R1CS constraint is satisfied.
+ *
+ * NOTE:
+ * The 0-th variable (i.e., "x_{0}") always represents the constant 1.
+ * Thus, the 0-th variable is not included in num_variables.
+ */
+template<typename FieldT>
+class r1cs_constraint_system {
+public:
+    size_t primary_input_size;
+    size_t auxiliary_input_size;
+
+    std::vector<r1cs_constraint<FieldT> > constraints;
+
+    r1cs_constraint_system() : primary_input_size(0), auxiliary_input_size(0) {}
+
+    size_t num_inputs() const;
+    size_t num_variables() const;
+    size_t num_constraints() const;
+
+#ifdef DEBUG
+    std::map<size_t, std::string> constraint_annotations;
+    std::map<size_t, std::string> variable_annotations;
+#endif
+
+    bool is_valid() const;
+    bool is_satisfied(const r1cs_primary_input<FieldT> &primary_input,
+                      const r1cs_auxiliary_input<FieldT> &auxiliary_input) const;
+
+    void add_constraint(const r1cs_constraint<FieldT> &c);
+    void add_constraint(const r1cs_constraint<FieldT> &c, const std::string &annotation);
+
+    void swap_AB_if_beneficial();
+
+    bool operator==(const r1cs_constraint_system<FieldT> &other) const;
+
+    friend std::ostream& operator<< <FieldT>(std::ostream &out, const r1cs_constraint_system<FieldT> &cs);
+    friend std::istream& operator>> <FieldT>(std::istream &in, r1cs_constraint_system<FieldT> &cs);
+
+    void report_linear_constraint_statistics() const;
+};
+
+
+} // libsnark
+
+#include "relations/constraint_satisfaction_problems/r1cs/r1cs.tcc"
+
+#endif // R1CS_HPP_

src/snark/libsnark/relations/constraint_satisfaction_problems/r1cs/r1cs.tcc

@@ -0,0 +1,310 @@
+/** @file
+ *****************************************************************************
+
+ Declaration of interfaces for:
+ - a R1CS constraint,
+ - a R1CS variable assignment, and
+ - a R1CS constraint system.
+
+ See r1cs.hpp .
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef R1CS_TCC_
+#define R1CS_TCC_
+
+#include <algorithm>
+#include <cassert>
+#include <set>
+#include "common/utils.hpp"
+#include "common/profiling.hpp"
+#include "algebra/fields/bigint.hpp"
+
+namespace libsnark {
+
+template<typename FieldT>
+r1cs_constraint<FieldT>::r1cs_constraint(const linear_combination<FieldT> &a,
+                                         const linear_combination<FieldT> &b,
+                                         const linear_combination<FieldT> &c) :
+    a(a), b(b), c(c)
+{
+}
+
+template<typename FieldT>
+r1cs_constraint<FieldT>::r1cs_constraint(const std::initializer_list<linear_combination<FieldT> > &A,
+                                         const std::initializer_list<linear_combination<FieldT> > &B,
+                                         const std::initializer_list<linear_combination<FieldT> > &C)
+{
+    for (auto lc_A : A)
+    {
+        a.terms.insert(a.terms.end(), lc_A.terms.begin(), lc_A.terms.end());
+    }
+    for (auto lc_B : B)
+    {
+        b.terms.insert(b.terms.end(), lc_B.terms.begin(), lc_B.terms.end());
+    }
+    for (auto lc_C : C)
+    {
+        c.terms.insert(c.terms.end(), lc_C.terms.begin(), lc_C.terms.end());
+    }
+}
+
+template<typename FieldT>
+bool r1cs_constraint<FieldT>::operator==(const r1cs_constraint<FieldT> &other) const
+{
+    return (this->a == other.a &&
+            this->b == other.b &&
+            this->c == other.c);
+}
+
+template<typename FieldT>
+std::ostream& operator<<(std::ostream &out, const r1cs_constraint<FieldT> &c)
+{
+    out << c.a;
+    out << c.b;
+    out << c.c;
+
+    return out;
+}
+
+template<typename FieldT>
+std::istream& operator>>(std::istream &in, r1cs_constraint<FieldT> &c)
+{
+    in >> c.a;
+    in >> c.b;
+    in >> c.c;
+
+    return in;
+}
+
+template<typename FieldT>
+size_t r1cs_constraint_system<FieldT>::num_inputs() const
+{
+    return primary_input_size;
+}
+
+template<typename FieldT>
+size_t r1cs_constraint_system<FieldT>::num_variables() const
+{
+    return primary_input_size + auxiliary_input_size;
+}
+
+
+template<typename FieldT>
+size_t r1cs_constraint_system<FieldT>::num_constraints() const
+{
+    return constraints.size();
+}
+
+template<typename FieldT>
+bool r1cs_constraint_system<FieldT>::is_valid() const
+{
+    if (this->num_inputs() > this->num_variables()) return false;
+
+    for (size_t c = 0; c < constraints.size(); ++c)
+    {
+        if (!(constraints[c].a.is_valid(this->num_variables()) &&
+              constraints[c].b.is_valid(this->num_variables()) &&
+              constraints[c].c.is_valid(this->num_variables())))
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+template<typename FieldT>
+void dump_r1cs_constraint(const r1cs_constraint<FieldT> &constraint,
+                          const r1cs_variable_assignment<FieldT> &full_variable_assignment,
+                          const std::map<size_t, std::string> &variable_annotations)
+{
+    printf("terms for a:\n"); constraint.a.print_with_assignment(full_variable_assignment, variable_annotations);
+    printf("terms for b:\n"); constraint.b.print_with_assignment(full_variable_assignment, variable_annotations);
+    printf("terms for c:\n"); constraint.c.print_with_assignment(full_variable_assignment, variable_annotations);
+}
+
+template<typename FieldT>
+bool r1cs_constraint_system<FieldT>::is_satisfied(const r1cs_primary_input<FieldT> &primary_input,
+                                                  const r1cs_auxiliary_input<FieldT> &auxiliary_input) const
+{
+    assert(primary_input.size() == num_inputs());
+    assert(primary_input.size() + auxiliary_input.size() == num_variables());
+
+    r1cs_variable_assignment<FieldT> full_variable_assignment = primary_input;
+    full_variable_assignment.insert(full_variable_assignment.end(), auxiliary_input.begin(), auxiliary_input.end());
+
+    for (size_t c = 0; c < constraints.size(); ++c)
+    {
+        const FieldT ares = constraints[c].a.evaluate(full_variable_assignment);
+        const FieldT bres = constraints[c].b.evaluate(full_variable_assignment);
+        const FieldT cres = constraints[c].c.evaluate(full_variable_assignment);
+
+        if (!(ares*bres == cres))
+        {
+#ifdef DEBUG
+            auto it = constraint_annotations.find(c);
+            printf("constraint %zu (%s) unsatisfied\n", c, (it == constraint_annotations.end() ? "no annotation" : it->second.c_str()));
+            printf("<a,(1,x)> = "); ares.print();
+            printf("<b,(1,x)> = "); bres.print();
+            printf("<c,(1,x)> = "); cres.print();
+            printf("constraint was:\n");
+            dump_r1cs_constraint(constraints[c], full_variable_assignment, variable_annotations);
+#endif // DEBUG
+            return false;
+        }
+    }
+
+    return true;
+}
+
+template<typename FieldT>
+void r1cs_constraint_system<FieldT>::add_constraint(const r1cs_constraint<FieldT> &c)
+{
+    constraints.emplace_back(c);
+}
+
+template<typename FieldT>
+void r1cs_constraint_system<FieldT>::add_constraint(const r1cs_constraint<FieldT> &c, const std::string &annotation)
+{
+#ifdef DEBUG
+    constraint_annotations[constraints.size()] = annotation;
+#endif
+    constraints.emplace_back(c);
+}
+
+template<typename FieldT>
+void r1cs_constraint_system<FieldT>::swap_AB_if_beneficial()
+{
+    enter_block("Call to r1cs_constraint_system::swap_AB_if_beneficial");
+
+    enter_block("Estimate densities");
+    bit_vector touched_by_A(this->num_variables() + 1, false), touched_by_B(this->num_variables() + 1, false);
+
+    for (size_t i = 0; i < this->constraints.size(); ++i)
+    {
+        for (size_t j = 0; j < this->constraints[i].a.terms.size(); ++j)
+        {
+            touched_by_A[this->constraints[i].a.terms[j].index] = true;
+        }
+
+        for (size_t j = 0; j < this->constraints[i].b.terms.size(); ++j)
+        {
+            touched_by_B[this->constraints[i].b.terms[j].index] = true;
+        }
+    }
+
+    size_t non_zero_A_count = 0, non_zero_B_count = 0;
+    for (size_t i = 0; i < this->num_variables() + 1; ++i)
+    {
+        non_zero_A_count += touched_by_A[i] ? 1 : 0;
+        non_zero_B_count += touched_by_B[i] ? 1 : 0;
+    }
+
+    if (!inhibit_profiling_info)
+    {
+        print_indent(); printf("* Non-zero A-count (estimate): %zu\n", non_zero_A_count);
+        print_indent(); printf("* Non-zero B-count (estimate): %zu\n", non_zero_B_count);
+    }
+    leave_block("Estimate densities");
+
+    if (non_zero_B_count > non_zero_A_count)
+    {
+        enter_block("Perform the swap");
+        for (size_t i = 0; i < this->constraints.size(); ++i)
+        {
+            std::swap(this->constraints[i].a, this->constraints[i].b);
+        }
+        leave_block("Perform the swap");
+    }
+    else
+    {
+        print_indent(); printf("Swap is not beneficial, not performing\n");
+    }
+
+    leave_block("Call to r1cs_constraint_system::swap_AB_if_beneficial");
+}
+
+template<typename FieldT>
+bool r1cs_constraint_system<FieldT>::operator==(const r1cs_constraint_system<FieldT> &other) const
+{
+    return (this->constraints == other.constraints &&
+            this->primary_input_size == other.primary_input_size &&
+            this->auxiliary_input_size == other.auxiliary_input_size);
+}
+
+template<typename FieldT>
+std::ostream& operator<<(std::ostream &out, const r1cs_constraint_system<FieldT> &cs)
+{
+    out << cs.primary_input_size << "\n";
+    out << cs.auxiliary_input_size << "\n";
+
+    out << cs.num_constraints() << "\n";
+    for (const r1cs_constraint<FieldT>& c : cs.constraints)
+    {
+        out << c;
+    }
+
+    return out;
+}
+
+template<typename FieldT>
+std::istream& operator>>(std::istream &in, r1cs_constraint_system<FieldT> &cs)
+{
+    in >> cs.primary_input_size;
+    in >> cs.auxiliary_input_size;
+
+    cs.constraints.clear();
+
+    size_t s;
+    in >> s;
+
+    char b;
+    in.read(&b, 1);
+
+    cs.constraints.reserve(s);
+
+    for (size_t i = 0; i < s; ++i)
+    {
+        r1cs_constraint<FieldT> c;
+        in >> c;
+        cs.constraints.emplace_back(c);
+    }
+
+    return in;
+}
+
+template<typename FieldT>
+void r1cs_constraint_system<FieldT>::report_linear_constraint_statistics() const
+{
+#ifdef DEBUG
+    for (size_t i = 0; i < constraints.size(); ++i)
+    {
+        auto &constr = constraints[i];
+        bool a_is_const = true;
+        for (auto &t : constr.a.terms)
+        {
+            a_is_const = a_is_const && (t.index == 0);
+        }
+
+        bool b_is_const = true;
+        for (auto &t : constr.b.terms)
+        {
+            b_is_const = b_is_const && (t.index == 0);
+        }
+
+        if (a_is_const || b_is_const)
+        {
+            auto it = constraint_annotations.find(i);
+            printf("%s\n", (it == constraint_annotations.end() ? FORMAT("", "constraint_%zu", i) : it->second).c_str());
+        }
+    }
+#endif
+}
+
+} // libsnark
+#endif // R1CS_TCC_

src/snark/libsnark/relations/variable.hpp

@@ -0,0 +1,213 @@
+/** @file
+ *****************************************************************************
+
+ Declaration of interfaces for:
+ - a variable (i.e., x_i),
+ - a linear term (i.e., a_i * x_i), and
+ - a linear combination (i.e., sum_i a_i * x_i).
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef VARIABLE_HPP_
+#define VARIABLE_HPP_
+
+#include <cstddef>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace libsnark {
+
+/**
+ * Mnemonic typedefs.
+ */
+typedef size_t var_index_t;
+typedef long integer_coeff_t;
+
+/**
+ * Forward declaration.
+ */
+template<typename FieldT>
+class linear_term;
+
+/**
+ * Forward declaration.
+ */
+template<typename FieldT>
+class linear_combination;
+
+/********************************* Variable **********************************/
+
+/**
+ * A variable represents a formal expresison of the form "x_{index}".
+ */
+template<typename FieldT>
+class variable {
+public:
+
+    var_index_t index;
+
+    variable(const var_index_t index = 0) : index(index) {};
+
+    linear_term<FieldT> operator*(const integer_coeff_t int_coeff) const;
+    linear_term<FieldT> operator*(const FieldT &field_coeff) const;
+
+    linear_combination<FieldT> operator+(const linear_combination<FieldT> &other) const;
+    linear_combination<FieldT> operator-(const linear_combination<FieldT> &other) const;
+
+    linear_term<FieldT> operator-() const;
+
+    bool operator==(const variable<FieldT> &other) const;
+};
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const integer_coeff_t int_coeff, const variable<FieldT> &var);
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const FieldT &field_coeff, const variable<FieldT> &var);
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const integer_coeff_t int_coeff, const variable<FieldT> &var);
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const FieldT &field_coeff, const variable<FieldT> &var);
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const integer_coeff_t int_coeff, const variable<FieldT> &var);
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const FieldT &field_coeff, const variable<FieldT> &var);
+
+
+/****************************** Linear term **********************************/
+
+/**
+ * A linear term represents a formal expression of the form "coeff * x_{index}".
+ */
+template<typename FieldT>
+class linear_term {
+public:
+
+    var_index_t index = 0;
+    FieldT coeff;
+
+    linear_term() {};
+    linear_term(const variable<FieldT> &var);
+    linear_term(const variable<FieldT> &var, const integer_coeff_t int_coeff);
+    linear_term(const variable<FieldT> &var, const FieldT &field_coeff);
+
+    linear_term<FieldT> operator*(const integer_coeff_t int_coeff) const;
+    linear_term<FieldT> operator*(const FieldT &field_coeff) const;
+
+    linear_combination<FieldT> operator+(const linear_combination<FieldT> &other) const;
+    linear_combination<FieldT> operator-(const linear_combination<FieldT> &other) const;
+
+    linear_term<FieldT> operator-() const;
+
+    bool operator==(const linear_term<FieldT> &other) const;
+};
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const integer_coeff_t int_coeff, const linear_term<FieldT> &lt);
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const FieldT &field_coeff, const linear_term<FieldT> &lt);
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const integer_coeff_t int_coeff, const linear_term<FieldT> &lt);
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const FieldT &field_coeff, const linear_term<FieldT> &lt);
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const integer_coeff_t int_coeff, const linear_term<FieldT> &lt);
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const FieldT &field_coeff, const linear_term<FieldT> &lt);
+
+
+/***************************** Linear combination ****************************/
+
+template<typename FieldT>
+class linear_combination;
+
+template<typename FieldT>
+std::ostream& operator<<(std::ostream &out, const linear_combination<FieldT> &lc);
+
+template<typename FieldT>
+std::istream& operator>>(std::istream &in, linear_combination<FieldT> &lc);
+
+/**
+ * A linear combination represents a formal expression of the form "sum_i coeff_i * x_{index_i}".
+ */
+template<typename FieldT>
+class linear_combination {
+public:
+
+    std::vector<linear_term<FieldT> > terms;
+
+    linear_combination() {};
+    linear_combination(const integer_coeff_t int_coeff);
+    linear_combination(const FieldT &field_coeff);
+    linear_combination(const variable<FieldT> &var);
+    linear_combination(const linear_term<FieldT> &lt);
+    linear_combination(const std::vector<linear_term<FieldT> > &all_terms);
+
+    /* for supporting range-based for loops over linear_combination */
+    typename std::vector<linear_term<FieldT> >::const_iterator begin() const;
+    typename std::vector<linear_term<FieldT> >::const_iterator end() const;
+
+    void add_term(const variable<FieldT> &var);
+    void add_term(const variable<FieldT> &var, const integer_coeff_t int_coeff);
+    void add_term(const variable<FieldT> &var, const FieldT &field_coeff);
+
+    void add_term(const linear_term<FieldT> &lt);
+
+    FieldT evaluate(const std::vector<FieldT> &assignment) const;
+
+    linear_combination<FieldT> operator*(const integer_coeff_t int_coeff) const;
+    linear_combination<FieldT> operator*(const FieldT &field_coeff) const;
+
+    linear_combination<FieldT> operator+(const linear_combination<FieldT> &other) const;
+
+    linear_combination<FieldT> operator-(const linear_combination<FieldT> &other) const;
+    linear_combination<FieldT> operator-() const;
+
+    bool operator==(const linear_combination<FieldT> &other) const;
+
+    bool is_valid(const size_t num_variables) const;
+
+    void print(const std::map<size_t, std::string> &variable_annotations = std::map<size_t, std::string>()) const;
+    void print_with_assignment(const std::vector<FieldT> &full_assignment, const std::map<size_t, std::string> &variable_annotations = std::map<size_t, std::string>()) const;
+
+    friend std::ostream& operator<< <FieldT>(std::ostream &out, const linear_combination<FieldT> &lc);
+    friend std::istream& operator>> <FieldT>(std::istream &in, linear_combination<FieldT> &lc);
+};
+
+template<typename FieldT>
+linear_combination<FieldT> operator*(const integer_coeff_t int_coeff, const linear_combination<FieldT> &lc);
+
+template<typename FieldT>
+linear_combination<FieldT> operator*(const FieldT &field_coeff, const linear_combination<FieldT> &lc);
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const integer_coeff_t int_coeff, const linear_combination<FieldT> &lc);
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const FieldT &field_coeff, const linear_combination<FieldT> &lc);
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const integer_coeff_t int_coeff, const linear_combination<FieldT> &lc);
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const FieldT &field_coeff, const linear_combination<FieldT> &lc);
+
+} // libsnark
+
+#include "relations/variable.tcc"
+
+#endif // VARIABLE_HPP_

src/snark/libsnark/relations/variable.tcc

@@ -0,0 +1,512 @@
+/** @file
+ *****************************************************************************
+
+ Implementation of interfaces for:
+ - a variable (i.e., x_i),
+ - a linear term (i.e., a_i * x_i), and
+ - a linear combination (i.e., sum_i a_i * x_i).
+
+ See variabe.hpp .
+
+ *****************************************************************************
+ * @author     This file is part of libsnark, developed by SCIPR Lab
+ *             and contributors (see AUTHORS).
+ * @copyright  MIT license (see LICENSE file)
+ *****************************************************************************/
+
+#ifndef VARIABLE_TCC_
+#define VARIABLE_TCC_
+
+#include <algorithm>
+#include <cassert>
+
+#include "algebra/fields/bigint.hpp"
+
+namespace libsnark {
+
+template<typename FieldT>
+linear_term<FieldT> variable<FieldT>::operator*(const integer_coeff_t int_coeff) const
+{
+    return linear_term<FieldT>(*this, int_coeff);
+}
+
+template<typename FieldT>
+linear_term<FieldT> variable<FieldT>::operator*(const FieldT &field_coeff) const
+{
+    return linear_term<FieldT>(*this, field_coeff);
+}
+
+template<typename FieldT>
+linear_combination<FieldT> variable<FieldT>::operator+(const linear_combination<FieldT> &other) const
+{
+    linear_combination<FieldT> result;
+
+    result.add_term(*this);
+    result.terms.insert(result.terms.begin(), other.terms.begin(), other.terms.end());
+
+    return result;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> variable<FieldT>::operator-(const linear_combination<FieldT> &other) const
+{
+    return (*this) + (-other);
+}
+
+template<typename FieldT>
+linear_term<FieldT> variable<FieldT>::operator-() const
+{
+    return linear_term<FieldT>(*this, -FieldT::one());
+}
+
+template<typename FieldT>
+bool variable<FieldT>::operator==(const variable<FieldT> &other) const
+{
+    return (this->index == other.index);
+}
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const integer_coeff_t int_coeff, const variable<FieldT> &var)
+{
+    return linear_term<FieldT>(var, int_coeff);
+}
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const FieldT &field_coeff, const variable<FieldT> &var)
+{
+    return linear_term<FieldT>(var, field_coeff);
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const integer_coeff_t int_coeff, const variable<FieldT> &var)
+{
+    return linear_combination<FieldT>(int_coeff) + var;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const FieldT &field_coeff, const variable<FieldT> &var)
+{
+    return linear_combination<FieldT>(field_coeff) + var;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const integer_coeff_t int_coeff, const variable<FieldT> &var)
+{
+    return linear_combination<FieldT>(int_coeff) - var;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const FieldT &field_coeff, const variable<FieldT> &var)
+{
+    return linear_combination<FieldT>(field_coeff) - var;
+}
+
+template<typename FieldT>
+linear_term<FieldT>::linear_term(const variable<FieldT> &var) :
+    index(var.index), coeff(FieldT::one())
+{
+}
+
+template<typename FieldT>
+linear_term<FieldT>::linear_term(const variable<FieldT> &var, const integer_coeff_t int_coeff) :
+    index(var.index), coeff(FieldT(int_coeff))
+{
+}
+
+template<typename FieldT>
+linear_term<FieldT>::linear_term(const variable<FieldT> &var, const FieldT &coeff) :
+    index(var.index), coeff(coeff)
+{
+}
+
+template<typename FieldT>
+linear_term<FieldT> linear_term<FieldT>::operator*(const integer_coeff_t int_coeff) const
+{
+    return (this->operator*(FieldT(int_coeff)));
+}
+
+template<typename FieldT>
+linear_term<FieldT> linear_term<FieldT>::operator*(const FieldT &field_coeff) const
+{
+    return linear_term<FieldT>(this->index, field_coeff * this->coeff);
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const integer_coeff_t int_coeff, const linear_term<FieldT> &lt)
+{
+    return linear_combination<FieldT>(int_coeff) + lt;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const FieldT &field_coeff, const linear_term<FieldT> &lt)
+{
+    return linear_combination<FieldT>(field_coeff) + lt;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const integer_coeff_t int_coeff, const linear_term<FieldT> &lt)
+{
+    return linear_combination<FieldT>(int_coeff) - lt;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const FieldT &field_coeff, const linear_term<FieldT> &lt)
+{
+    return linear_combination<FieldT>(field_coeff) - lt;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_term<FieldT>::operator+(const linear_combination<FieldT> &other) const
+{
+    return linear_combination<FieldT>(*this) + other;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_term<FieldT>::operator-(const linear_combination<FieldT> &other) const
+{
+    return (*this) + (-other);
+}
+
+template<typename FieldT>
+linear_term<FieldT> linear_term<FieldT>::operator-() const
+{
+    return linear_term<FieldT>(this->index, -this->coeff);
+}
+
+template<typename FieldT>
+bool linear_term<FieldT>::operator==(const linear_term<FieldT> &other) const
+{
+    return (this->index == other.index &&
+            this->coeff == other.coeff);
+}
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const integer_coeff_t int_coeff, const linear_term<FieldT> &lt)
+{
+    return FieldT(int_coeff) * lt;
+}
+
+template<typename FieldT>
+linear_term<FieldT> operator*(const FieldT &field_coeff, const linear_term<FieldT> &lt)
+{
+    return linear_term<FieldT>(lt.index, field_coeff * lt.coeff);
+}
+
+template<typename FieldT>
+linear_combination<FieldT>::linear_combination(const integer_coeff_t int_coeff)
+{
+    this->add_term(linear_term<FieldT>(0, int_coeff));
+}
+
+template<typename FieldT>
+linear_combination<FieldT>::linear_combination(const FieldT &field_coeff)
+{
+    this->add_term(linear_term<FieldT>(0, field_coeff));
+}
+
+template<typename FieldT>
+linear_combination<FieldT>::linear_combination(const variable<FieldT> &var)
+{
+    this->add_term(var);
+}
+
+template<typename FieldT>
+linear_combination<FieldT>::linear_combination(const linear_term<FieldT> &lt)
+{
+    this->add_term(lt);
+}
+
+template<typename FieldT>
+typename std::vector<linear_term<FieldT> >::const_iterator linear_combination<FieldT>::begin() const
+{
+    return terms.begin();
+}
+
+template<typename FieldT>
+typename std::vector<linear_term<FieldT> >::const_iterator linear_combination<FieldT>::end() const
+{
+    return terms.end();
+}
+
+template<typename FieldT>
+void linear_combination<FieldT>::add_term(const variable<FieldT> &var)
+{
+    this->terms.emplace_back(linear_term<FieldT>(var.index, FieldT::one()));
+}
+
+template<typename FieldT>
+void linear_combination<FieldT>::add_term(const variable<FieldT> &var, const integer_coeff_t int_coeff)
+{
+    this->terms.emplace_back(linear_term<FieldT>(var.index, int_coeff));
+}
+
+template<typename FieldT>
+void linear_combination<FieldT>::add_term(const variable<FieldT> &var, const FieldT &coeff)
+{
+    this->terms.emplace_back(linear_term<FieldT>(var.index, coeff));
+}
+
+template<typename FieldT>
+void linear_combination<FieldT>::add_term(const linear_term<FieldT> &other)
+{
+    this->terms.emplace_back(other);
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_combination<FieldT>::operator*(const integer_coeff_t int_coeff) const
+{
+    return (*this) * FieldT(int_coeff);
+}
+
+template<typename FieldT>
+FieldT linear_combination<FieldT>::evaluate(const std::vector<FieldT> &assignment) const
+{
+    FieldT acc = FieldT::zero();
+    for (auto &lt : terms)
+    {
+        acc += (lt.index == 0 ? FieldT::one() : assignment[lt.index-1]) * lt.coeff;
+    }
+    return acc;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_combination<FieldT>::operator*(const FieldT &field_coeff) const
+{
+    linear_combination<FieldT> result;
+    result.terms.reserve(this->terms.size());
+    for (const linear_term<FieldT> &lt : this->terms)
+    {
+        result.terms.emplace_back(lt * field_coeff);
+    }
+    return result;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_combination<FieldT>::operator+(const linear_combination<FieldT> &other) const
+{
+    linear_combination<FieldT> result;
+
+    auto it1 = this->terms.begin();
+    auto it2 = other.terms.begin();
+
+    /* invariant: it1 and it2 always point to unprocessed items in the corresponding linear combinations */
+    while (it1 != this->terms.end() && it2 != other.terms.end())
+    {
+        if (it1->index < it2->index)
+        {
+            result.terms.emplace_back(*it1);
+            ++it1;
+        }
+        else if (it1->index > it2->index)
+        {
+            result.terms.emplace_back(*it2);
+            ++it2;
+        }
+        else
+        {
+            /* it1->index == it2->index */
+            result.terms.emplace_back(linear_term<FieldT>(variable<FieldT>(it1->index), it1->coeff + it2->coeff));
+            ++it1;
+            ++it2;
+        }
+    }
+
+    if (it1 != this->terms.end())
+    {
+        result.terms.insert(result.terms.end(), it1, this->terms.end());
+    }
+    else
+    {
+        result.terms.insert(result.terms.end(), it2, other.terms.end());
+    }
+
+    return result;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_combination<FieldT>::operator-(const linear_combination<FieldT> &other) const
+{
+    return (*this) + (-other);
+}
+
+template<typename FieldT>
+linear_combination<FieldT> linear_combination<FieldT>::operator-() const
+{
+    return (*this) * (-FieldT::one());
+}
+
+template<typename FieldT>
+bool linear_combination<FieldT>::operator==(const linear_combination<FieldT> &other) const
+{
+    return (this->terms == other.terms);
+}
+
+template<typename FieldT>
+bool linear_combination<FieldT>::is_valid(const size_t num_variables) const
+{
+    /* check that all terms in linear combination are sorted */
+    for (size_t i = 1; i < terms.size(); ++i)
+    {
+        if (terms[i-1].index >= terms[i].index)
+        {
+            return false;
+        }
+    }
+
+    /* check that the variables are in proper range. as the variables
+       are sorted, it suffices to check the last term */
+    if ((--terms.end())->index >= num_variables)
+    {
+        return false;
+    }
+
+    return true;
+}
+
+template<typename FieldT>
+void linear_combination<FieldT>::print(const std::map<size_t, std::string> &variable_annotations) const
+{
+    for (auto &lt : terms)
+    {
+        if (lt.index == 0)
+        {
+            printf("    1 * ");
+            lt.coeff.print();
+        }
+        else
+        {
+            auto it = variable_annotations.find(lt.index);
+            printf("    x_%zu (%s) * ", lt.index, (it == variable_annotations.end() ? "no annotation" : it->second.c_str()));
+            lt.coeff.print();
+        }
+    }
+}
+
+template<typename FieldT>
+void linear_combination<FieldT>::print_with_assignment(const std::vector<FieldT> &full_assignment, const std::map<size_t, std::string> &variable_annotations) const
+{
+    for (auto &lt : terms)
+    {
+        if (lt.index == 0)
+        {
+            printf("    1 * ");
+            lt.coeff.print();
+        }
+        else
+        {
+            printf("    x_%zu * ", lt.index);
+            lt.coeff.print();
+
+            auto it = variable_annotations.find(lt.index);
+            printf("    where x_%zu (%s) was assigned value ", lt.index,
+                   (it == variable_annotations.end() ? "no annotation" : it->second.c_str()));
+            full_assignment[lt.index-1].print();
+            printf("      i.e. negative of ");
+            (-full_assignment[lt.index-1]).print();
+        }
+    }
+}
+
+template<typename FieldT>
+std::ostream& operator<<(std::ostream &out, const linear_combination<FieldT> &lc)
+{
+    out << lc.terms.size() << "\n";
+    for (const linear_term<FieldT>& lt : lc.terms)
+    {
+        out << lt.index << "\n";
+        out << lt.coeff << OUTPUT_NEWLINE;
+    }
+
+    return out;
+}
+
+template<typename FieldT>
+std::istream& operator>>(std::istream &in, linear_combination<FieldT> &lc)
+{
+    lc.terms.clear();
+
+    size_t s;
+    in >> s;
+
+    consume_newline(in);
+
+    lc.terms.reserve(s);
+
+    for (size_t i = 0; i < s; ++i)
+    {
+        linear_term<FieldT> lt;
+        in >> lt.index;
+        consume_newline(in);
+        in >> lt.coeff;
+        consume_OUTPUT_NEWLINE(in);
+        lc.terms.emplace_back(lt);
+    }
+
+    return in;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator*(const integer_coeff_t int_coeff, const linear_combination<FieldT> &lc)
+{
+    return lc * int_coeff;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator*(const FieldT &field_coeff, const linear_combination<FieldT> &lc)
+{
+    return lc * field_coeff;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const integer_coeff_t int_coeff, const linear_combination<FieldT> &lc)
+{
+    return linear_combination<FieldT>(int_coeff) + lc;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator+(const FieldT &field_coeff, const linear_combination<FieldT> &lc)
+{
+    return linear_combination<FieldT>(field_coeff) + lc;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const integer_coeff_t int_coeff, const linear_combination<FieldT> &lc)
+{
+    return linear_combination<FieldT>(int_coeff) - lc;
+}
+
+template<typename FieldT>
+linear_combination<FieldT> operator-(const FieldT &field_coeff, const linear_combination<FieldT> &lc)
+{
+    return linear_combination<FieldT>(field_coeff) - lc;
+}
+
+template<typename FieldT>
+linear_combination<FieldT>::linear_combination(const std::vector<linear_term<FieldT> > &all_terms)
+{
+    if (all_terms.empty())
+    {
+        return;
+    }
+
+    terms = all_terms;
+    std::sort(terms.begin(), terms.end(), [](linear_term<FieldT> a, linear_term<FieldT> b) { return a.index < b.index; });
+
+    auto result_it = terms.begin();
+    for (auto it = ++terms.begin(); it != terms.end(); ++it)
+    {
+        if (it->index == result_it->index)
+        {
+            result_it->coeff += it->coeff;
+        }
+        else
+        {
+            *(++result_it) = *it;
+        }
+    }
+    terms.resize((result_it - terms.begin()) + 1);
+}
+
+} // libsnark
+
+#endif // VARIABLE_TCC