pedersen-libsnark/ecgadget.hpp

#include "libsnark_headers.hpp"

using namespace libsnark;

// There are two types of values:
//   _constants_ are values known at circuit generation time; they
//       are global constants known to everyone
//   _variables_ are values that change in each use of the circuit;
//       they have two subtypes:
//
//   _public variables_ are values known to both the prover
//       and verifier but change in each use of the circuit
//   _private variables_ are values known only to the prover
//       and change in each use of the circuit

// Double a constant EC point (inx,iny) to yield (outx,outy).  The input
// point must not be the point at infinity.
template<typename FieldT>
static void ec_double_point(FieldT &outx, FieldT &outy,
    const FieldT &inx, const FieldT &iny)
{
    FieldT xsq = inx.squared();
    FieldT lambda = (xsq * 3 - 3) * (iny * 2).inverse();
    outx = lambda.squared() - inx * 2;
    outy = lambda * (inx - outx) - iny;
}

// Add constant EC points (inx, iny) and (addx, addy) to yield (outx, outy).
// inx and addx must not be equal.
template<typename FieldT>
static void ec_add_points(FieldT &outx, FieldT &outy,
    const FieldT &inx, const FieldT &iny,
    const FieldT &addx, const FieldT &addy)
{
    FieldT lambda = (addy - iny) * (addx - inx).inverse();
    outx = lambda.squared() - (addx + inx);
    outy = lambda * (inx - outx) - iny;
}

// Double the variable EC point (inx,iny) to yield (outx,outy).  The
// input point must not be the point at infinity.
template<typename FieldT>
class ec_double_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> lambda, inxsq;
public:
  const pb_variable<FieldT> outx, outy, inx, iny;

  ec_double_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_linear_combination<FieldT> &inx,
              const pb_linear_combination<FieldT> &iny) :
    gadget<FieldT>(pb, "ec_double_gadget"), outx(outx), outy(outy),
        inx(inx), iny(iny)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    lambda.allocate(this->pb, "lambda");
    inxsq.allocate(this->pb, "inxsq");
  }

  void generate_r1cs_constraints()
  {
    // inxsq = inx * inx
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(inx, inx, inxsq));

    // 2 * iny * lambda = 3 * inxsq - 3  (a = -3 on our curve)
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(2 * iny, lambda, 3 * inxsq - 3));

    // outx = lambda^2 - 2 * inx
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, lambda, outx + 2 * inx));

    // outy = lambda * (inx - outx) - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, inx - outx, outy + iny));

  }

  void generate_r1cs_witness()
  {
    this->pb.val(inxsq) = this->pb.lc_val(inx) * this->pb.lc_val(inx);
    this->pb.val(lambda) = (this->pb.val(inxsq) * 3 - 3) * (this->pb.lc_val(iny) * 2).inverse();
    this->pb.val(outx) = this->pb.val(lambda).squared() - this->pb.lc_val(inx) * 2;
    this->pb.val(outy) = this->pb.val(lambda) * (this->pb.lc_val(inx) - this->pb.val(outx)) - this->pb.lc_val(iny);
  }
};

// Add the variable EC point (addx,addy) to the variable EC point
// (inx,iny) to yield (outx,outy).  The input point must not be the
// point at infinity.
template<typename FieldT>
class ec_add_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> lambda;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_linear_combination<FieldT> inx, iny, addx, addy;

  ec_add_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_linear_combination<FieldT> &inx,
              const pb_linear_combination<FieldT> &iny,
              const pb_linear_combination<FieldT> &addx,
              const pb_linear_combination<FieldT> &addy) :
    gadget<FieldT>(pb, "ec_add_gadget"),
    outx(outx), outy(outy), inx(inx), iny(iny), addx(addx), addy(addy)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    lambda.allocate(this->pb, "lambda");
  }

  void generate_r1cs_constraints()
  {
    // (addx - inx) * lambda = addy - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(addx - inx, lambda, addy - iny));

    // outx = lambda^2 - (addx + inx)
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, lambda, outx + addx + inx));

    // outy = lambda * (inx - outx) - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, inx - outx, outy + iny));

  }

  void generate_r1cs_witness()
  {
    this->pb.val(lambda) = (this->pb.lc_val(addy) - this->pb.lc_val(iny)) * (this->pb.lc_val(addx) - this->pb.lc_val(inx)).inverse();
    this->pb.val(outx) = this->pb.val(lambda).squared() - (this->pb.lc_val(addx) + this->pb.lc_val(inx));
    this->pb.val(outy) = this->pb.val(lambda) * (this->pb.lc_val(inx) - this->pb.val(outx)) - this->pb.lc_val(iny);
  }
};

// Add the variable EC point P to the constant EC point (inx,iny) to
// yield (outx,outy).  The input point must not be the point at
// infinity.
template<typename FieldT>
class ec_constant_add_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> lambda;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_linear_combination<FieldT> inx, iny;
  const FieldT Px, Py;

  ec_constant_add_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_linear_combination<FieldT> &inx,
              const pb_linear_combination<FieldT> &iny,
              const FieldT &Px, const FieldT &Py) :
    gadget<FieldT>(pb, "ec_constant_add_gadget"),
    outx(outx), outy(outy), inx(inx), iny(iny), Px(Px), Py(Py)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    lambda.allocate(this->pb, "lambda");
  }

  void generate_r1cs_constraints()
  {
    // (Px - inx) * lambda = Py - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(Px - inx, lambda, Py - iny));

    // outx = lambda^2 - (Px + inx)
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, lambda, outx + Px + inx));

    // outy = lambda * (inx - outx) - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, inx - outx, outy + iny));

  }

  void generate_r1cs_witness()
  {
    this->pb.val(lambda) = (Py - this->pb.lc_val(iny)) * (Px - this->pb.lc_val(inx)).inverse();
    this->pb.val(outx) = this->pb.val(lambda).squared() - (Px + this->pb.lc_val(inx));
    this->pb.val(outy) = this->pb.val(lambda) * (this->pb.lc_val(inx) - this->pb.val(outx)) - this->pb.lc_val(iny);
  }
};

// Add the constant EC point P0 or the constant EC point P1 to the
// variable EC point (inx,iny) to yield (outx,outy).  The input point
// must not be the point at infinity.  The input bit which controls
// which addition is done.
template<typename FieldT>
class ec_2_constant_add_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> sumx, sumy;
  pb_linear_combination<FieldT> addx, addy;
  std::vector<ec_add_gadget<FieldT> > adder;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_linear_combination<FieldT> inx, iny;
  const pb_variable<FieldT> which;
  const FieldT P0x, P0y, P1x, P1y;

  ec_2_constant_add_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_linear_combination<FieldT> &inx,
              const pb_linear_combination<FieldT> &iny,
              const pb_variable<FieldT> &which,
              const FieldT &P0x, const FieldT &P0y,
              const FieldT &P1x, const FieldT &P1y) :
    gadget<FieldT>(pb, "ec_2_constant_add_gadget"),
    outx(outx), outy(outy), inx(inx), iny(iny), which(which),
    P0x(P0x), P0y(P0y), P1x(P1x), P1y(P1y)
  {
    // Allocate variables to protoboard

    addx.assign(pb, which * (P1x-P0x) + P0x);
    addy.assign(pb, which * (P1y-P0y) + P0y);
    adder.emplace_back(this->pb, outx, outy, inx, iny, addx, addy);
  }

  void generate_r1cs_constraints()
  {
    adder[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    addx.evaluate(this->pb);
    addy.evaluate(this->pb);
    adder[0].generate_r1cs_witness();
  }
};

// Add the variable EC point P0 or the variable EC point P1 to the
// variable EC point (inx,iny) to yield (outx,outy).  The input point
// must not be the point at infinity.  The input bit which controls
// which addition is done.
template<typename FieldT>
class ec_2_add_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> sumx, sumy;
  pb_variable<FieldT> addx, addy;
  std::vector<ec_add_gadget<FieldT> > adder;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_linear_combination<FieldT> inx, iny;
  const pb_variable<FieldT> which;
  const pb_variable<FieldT> P0x, P0y, P1x, P1y;

  ec_2_add_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_linear_combination<FieldT> &inx,
              const pb_linear_combination<FieldT> &iny,
              const pb_variable<FieldT> &which,
              const pb_variable<FieldT> &P0x,
              const pb_variable<FieldT> &P0y,
              const pb_variable<FieldT> &P1x,
              const pb_variable<FieldT> &P1y) :
    gadget<FieldT>(pb, "ec_2_add_gadget"),
    outx(outx), outy(outy), inx(inx), iny(iny), which(which),
    P0x(P0x), P0y(P0y), P1x(P1x), P1y(P1y)
  {
    // Allocate variables to protoboard

    addx.allocate(this->pb, "addx");
    addy.allocate(this->pb, "addy");
    adder.emplace_back(this->pb, outx, outy, inx, iny, addx, addy);
  }

  void generate_r1cs_constraints()
  {
    // Set (addx,addy) = which ? (P0x, P0y) : (P1x, P1y)
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(P1x - P0x, which, addx - P0x));
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(P1y - P0y, which, addy - P0y));
    adder[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    this->pb.val(addx) = this->pb.val(which) ? this->pb.val(P1x) : this->pb.val(P0x);
    this->pb.val(addy) = this->pb.val(which) ? this->pb.val(P1y) : this->pb.val(P0y);
    adder[0].generate_r1cs_witness();
  }
};

#if 0
// Add nothing, or one of the constant EC points P1, P2, or P3 to the EC
// point (inx,iny) to yield (outx,outy).  The input point must not be
// the point at infinity.  The two input bits add1 and add2 control what
// is added.  Typically, P3 will equal P1+P2, in which case this gadget
// does two conditional constant adds simultaneously in just 6 constraints.
template<typename FieldT>
class ec_add_P123_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> lambda, sumx, sumy, move;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_linear_combination<FieldT> inx, iny;
  const pb_variable<FieldT> add1, add2;
  const FieldT P1x, P1y, P2x, P2y, P3x, P3y;

  ec_add_P123_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_linear_combination<FieldT> &inx,
              const pb_linear_combination<FieldT> &iny,
              const pb_variable<FieldT> &add1,
              const pb_variable<FieldT> &add2,
              const FieldT &P1x, const FieldT &P1y,
              const FieldT &P2x, const FieldT &P2y,
              const FieldT &P3x, const FieldT &P3y) : 
    gadget<FieldT>(pb, "ec_add_P123_gadget"),
    outx(outx), outy(outy), inx(inx), iny(iny), add1(add1), add2(add2),
    P1x(P1x), P1y(P1y), P2x(P2x), P2y(P2y), P3x(P3x), P3y(P3y)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    lambda.allocate(this->pb, "lambda");
    sumx.allocate(this->pb, "sumx");
    sumy.allocate(this->pb, "sumy");
    move.allocate(this->pb, "move");
  }

  void generate_r1cs_constraints()
  {
    // Strategy: if add1 = add2 = 0, we compute some nonsense but throw
    // it away later.  Otherwise, the coordinates of the point to add
    // are a _linear_ function of add1 and add2 (since P1, P2, and P3
    // are public constants)

    // In particular, the point to add is ( (P3x - P2x) * add1 + (P3x -
    // P1x) * add2 + (P1x + P2x - P3x), (P3y - P2y) * add1 + (P3y - P1y) *
    // add2 + (P1y + P2y - P3y))

    // (addx - inx) * lambda = addy - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>((P3x - P2x) * add1 + (P3x - P1x) * add2 + (P1x + P2x - P3x) - inx, lambda, (P3y - P2y) * add1 + (P3y - P1y) * add2 + (P1y + P2y - P3y) - iny));

    // sumx = lambda^2 - (addx + inx)
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, lambda, sumx + (P3x - P2x) * add1 + (P3x - P1x) * add2 + (P1x + P2x - P3x) + inx));

    // sumy = lambda * (inx - sumx) - iny
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, inx - sumx, sumy + iny));

    // Now we want to conditionally move the sum.  We want that
    // outx = (add1 || add2) ? sumx : inx
    // outy = (add1 || add2) ? sumy : iny

    // so we compute move = add1 || add2, and then
    // outx = inx + (sumx - inx) * move
    // outy = iny + (sumy - iny) * move

    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(1 - add1, 1 - add2, 1 - move));
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(sumx - inx, move, outx - inx));
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(sumy - iny, move, outy - iny));

  }

  void generate_r1cs_witness()
  {
    FieldT addxval = (P3x - P2x) * this->pb.val(add1) + (P3x - P1x) * this->pb.val(add2) + (P1x + P2x - P3x);
    FieldT addyval = (P3y - P2y) * this->pb.val(add1) + (P3y - P1y) * this->pb.val(add2) + (P1y + P2y - P3y);
    this->pb.val(lambda) = (addyval - this->pb.lc_val(iny)) * (addxval - this->pb.lc_val(inx)).inverse();
    this->pb.val(sumx) = this->pb.val(lambda).squared() - (addxval + this->pb.lc_val(inx));
    this->pb.val(sumy) = this->pb.val(lambda) * (this->pb.lc_val(inx) - this->pb.val(sumx)) - this->pb.lc_val(iny);

    bool a1 = this->pb.val(add1) != FieldT(0);
    bool a2 = this->pb.val(add2) != FieldT(0);
    this->pb.val(move) = a1 || a2;
    this->pb.val(outx) = (a1 || a2) ? this->pb.val(sumx) : this->pb.lc_val(inx);
    this->pb.val(outy) = (a1 || a2) ? this->pb.val(sumy) : this->pb.lc_val(iny);
  }
};
#endif

// Compute A + s*P as (outx, outy) for an accumulator A, a given
// constant point P, and s given as a bit vector.  The _caller_ is
// responsible for proving that the elements of svec are bits.  The
// (constant) accumulator excess (AXS) will be updated; when all the
// computations are complete, AXS should be subtracted from the
// accumulator A.
template<typename FieldT>
class ec_constant_scalarmul_vec_accum_gadget : public gadget<FieldT> {
private:
  FieldT Cx, Cy;
  pb_variable_array<FieldT> accumx, accumy;
  std::vector<ec_2_constant_add_gadget<FieldT> > twoadders;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable<FieldT> Ax, Ay;
  const pb_variable_array<FieldT> svec;
  const FieldT Px, Py;

  // Strategy: We compute (as compile-time constants) (powers of 2)
  // times P.  Based on each bit of s, we add one of the constant points
  // C or (2^i * P) + C to the accumulator, and regardless of s, add C
  // to the excess.

  ec_constant_scalarmul_vec_accum_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &Ax,
              const pb_variable<FieldT> &Ay,
              const pb_variable_array<FieldT> &svec,
              const FieldT &Px, const FieldT &Py,
              FieldT &AXSx, FieldT &AXSy) :
    gadget<FieldT>(pb, "ec_constant_scalarmul_vec_accum_gadget"),
    // Precomputed coordinates of C
    Cx(2),
    Cy("4950745124018817972378217179409499695353526031437053848725554590521829916331"),
    outx(outx), outy(outy), Ax(Ax), Ay(Ay), svec(svec), Px(Px), Py(Py)
  {
    size_t numbits = svec.size();
    accumx.allocate(this->pb, numbits-1, "accumx");
    accumy.allocate(this->pb, numbits-1, "accumy");

    FieldT twoiPx = Px, twoiPy = Py;
    size_t i = 0;

    while(i < numbits) {
        // Invariant: twoiP = 2^i * P
        FieldT twoiPCx, twoiPCy;
        ec_add_points(twoiPCx, twoiPCy, twoiPx, twoiPy, Cx, Cy);

        twoadders.emplace_back(this->pb,
            (i == numbits-1 ? outx : accumx[i]),
            (i == numbits-1 ? outy : accumy[i]),
            (i == 0 ? Ax : accumx[i-1]),
            (i == 0 ? Ay : accumy[i-1]),
            svec[i], Cx, Cy, twoiPCx, twoiPCy);

        FieldT newtwoiPx, newtwoiPy, newAXSx, newAXSy;
        ec_double_point(newtwoiPx, newtwoiPy, twoiPx, twoiPy);
        twoiPx = newtwoiPx;
        twoiPy = newtwoiPy;
        i += 1;
        ec_add_points(newAXSx, newAXSy, AXSx, AXSy, Cx, Cy);
        AXSx = newAXSx;
        AXSy = newAXSy;
    }
  }

  void generate_r1cs_constraints()
  {
    for (auto&& gadget : twoadders) {
        gadget.generate_r1cs_constraints();
    }
  }

  void generate_r1cs_witness()
  {
    for (auto&& gadget : twoadders) {
        gadget.generate_r1cs_witness();
    }
  }
};

// Compute A + s*P as (outx, outy) for an accumulator A, a given
// constant point P, and s given as a field element.  The (constant)
// accumulator excess (AXS) will be updated; when all the computations
// are complete, AXS should be subtracted from the accumulator A.
template<typename FieldT>
class ec_constant_scalarmul_accum_gadget : public gadget<FieldT> {
private:
  pb_variable_array<FieldT> svec;
  std::vector<packing_gadget<FieldT> > packers;
  std::vector<ec_constant_scalarmul_vec_accum_gadget<FieldT> > vecgadget;

public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable<FieldT> Ax, Ay;
  const pb_variable<FieldT> s;
  const FieldT Px, Py;

  ec_constant_scalarmul_accum_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &Ax,
              const pb_variable<FieldT> &Ay,
              const pb_variable<FieldT> &s,
              const FieldT &Px, const FieldT &Py,
              FieldT &AXSx, FieldT &AXSy) :
    gadget<FieldT>(pb, "ec_constant_scalarmul_accum_gadget"),
    outx(outx), outy(outy), Ax(Ax), Ay(Ay), s(s), Px(Px), Py(Py)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    size_t numbits = FieldT::num_bits;
    svec.allocate(this->pb, numbits, "svec");
    packers.emplace_back(this->pb, svec, s);
    vecgadget.emplace_back(this->pb, outx, outy, Ax, Ay, svec, Px, Py, AXSx, AXSy);
  }

  void generate_r1cs_constraints()
  {
    packers[0].generate_r1cs_constraints(true);
    vecgadget[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    packers[0].generate_r1cs_witness_from_packed();
    vecgadget[0].generate_r1cs_witness();
  }
};

// Compute s*P as (outx, outy) for a given constant point P, and s given
// as a bit vector.  The _caller_ is responsible for proving that the
// elements of svec are bits.
template<typename FieldT>
class ec_constant_scalarmul_vec_gadget : public gadget<FieldT> {
private:
  FieldT Cx, Cy, Ax, Ay, AXSx, AXSy;
  pb_variable<FieldT> accinx, acciny, accoutx, accouty;
  std::vector<ec_constant_scalarmul_vec_accum_gadget<FieldT> > scalarmuls;
  std::vector<ec_constant_add_gadget<FieldT> > adders;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable_array<FieldT> svec;
  const FieldT Px, Py;

  ec_constant_scalarmul_vec_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable_array<FieldT> &svec,
              const FieldT &Px, const FieldT &Py) :
    gadget<FieldT>(pb, "ec_constant_scalarmul_vec_gadget"),
    // Precomputed coordinates of C and A
    Cx(2),
    Cy("4950745124018817972378217179409499695353526031437053848725554590521829916331"),
    Ax("7536839002660211356286040193441766649532044555061394833845553337792579131020"),
    Ay("11391058648720923807988142436733355540810929560298907319389650598553246451302"),
    outx(outx), outy(outy), svec(svec), Px(Px), Py(Py)
  {
    AXSx = Ax;
    AXSy = Ay;
    accinx.allocate(this->pb, "accinx");
    acciny.allocate(this->pb, "acciny");
    accoutx.allocate(this->pb, "accoutx");
    accouty.allocate(this->pb, "accouty");

    scalarmuls.emplace_back(pb, accoutx, accouty, accinx, acciny, svec, Px, Py, AXSx, AXSy);
    adders.emplace_back(pb, outx, outy, accoutx, accouty, AXSx, -AXSy);
  }

  void generate_r1cs_constraints()
  {
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(accinx, 1, Ax));
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(acciny, 1, Ay));
    scalarmuls[0].generate_r1cs_constraints();
    adders[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    this->pb.val(accinx) = Ax;
    this->pb.val(acciny) = Ay;
    scalarmuls[0].generate_r1cs_witness();
    adders[0].generate_r1cs_witness();
  }
};

// Compute s*P as (outx, outy) for a given constant point P, and s given
// as a field element.
template<typename FieldT>
class ec_constant_scalarmul_gadget : public gadget<FieldT> {
private:
  pb_variable_array<FieldT> svec;
  std::vector<packing_gadget<FieldT> > packers;
  std::vector<ec_constant_scalarmul_vec_gadget<FieldT> > vecgadget;

public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable<FieldT> s;
  const FieldT Px, Py;

  ec_constant_scalarmul_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &s,
              const FieldT &Px, const FieldT &Py) :
    gadget<FieldT>(pb, "ec_constant_scalarmul_gadget"),
    outx(outx), outy(outy), s(s), Px(Px), Py(Py)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    size_t numbits = FieldT::num_bits;
    svec.allocate(this->pb, numbits, "svec");
    packers.emplace_back(this->pb, svec, s);
    vecgadget.emplace_back(this->pb, outx, outy, svec, Px, Py);
  }

  void generate_r1cs_constraints()
  {
    packers[0].generate_r1cs_constraints(true);
    vecgadget[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    packers[0].generate_r1cs_witness_from_packed();
    vecgadget[0].generate_r1cs_witness();
  }
};

// Compute A + s*P as (outx, outy) for an accumulator A, a precomputed
// addition table Ptable for a variable point P, and s given as a bit
// vector.  The _caller_ is responsible for proving that the elements of
// svec are bits.  The (constant) accumulator excess (AXS) will be
// updated; when all the computations are complete, AXS should be
// subtracted from the accumulator A.  The addition table is a variable
// array of length 2*numbits (where numbits is the length of svec) such
// that Ptable[2*i] and Ptable[2*i+1] are the (x,y) coordinates of
// 2^i * P + C.
template<typename FieldT>
class ec_scalarmul_vec_accum_gadget : public gadget<FieldT> {
private:
  FieldT Cx, Cy;
  pb_variable_array<FieldT> accumx, accumy;
  std::vector<ec_2_add_gadget<FieldT> > twoadders;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable<FieldT> Ax, Ay;
  const pb_variable_array<FieldT> svec, Ptable;

  ec_scalarmul_vec_accum_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &Ax,
              const pb_variable<FieldT> &Ay,
              const pb_variable_array<FieldT> &svec,
              const pb_variable_array<FieldT> &Ptable,
              FieldT &AXSx, FieldT &AXSy) :
    gadget<FieldT>(pb, "ec_scalarmul_vec_accum_gadget"),
    // Precomputed coordinates of C
    Cx(2),
    Cy("4950745124018817972378217179409499695353526031437053848725554590521829916331"),
    outx(outx), outy(outy), Ax(Ax), Ay(Ay), svec(svec), Ptable(Ptable)
  {
    size_t numbits = svec.size();
    assert(Ptable.size() == 2*numbits);
    accumx.allocate(this->pb, numbits-1, "accumx");
    accumy.allocate(this->pb, numbits-1, "accumy");

    for (size_t i = 0; i < numbits; ++i) {
        twoadders.emplace_back(this->pb,
            (i == numbits-1 ? outx : accumx[i]),
            (i == numbits-1 ? outy : accumy[i]),
            (i == 0 ? Ax : accumx[i-1]),
            (i == 0 ? Ay : accumy[i-1]),
            svec[i], Cx, Cy, Ptable[2*i], Ptable[2*i+1]);

        FieldT newAXSx, newAXSy;
        ec_add_points(newAXSx, newAXSy, AXSx, AXSy, Cx, Cy);
        AXSx = newAXSx;
        AXSy = newAXSy;
    }
  }

  void generate_r1cs_constraints()
  {
    for (auto&& gadget : twoadders) {
        gadget.generate_r1cs_constraints();
    }
  }

  void generate_r1cs_witness()
  {
    for (auto&& gadget : twoadders) {
        gadget.generate_r1cs_witness();
    }
  }
};

// Compute A + s*P as (outx, outy) for an accumulator A, a precomputed
// addition table Ptable for a variable point P, and s given as a field
// element.  The _caller_ is responsible for proving that the elements
// of svec are bits.  The (constant) accumulator excess (AXS) will be
// updated; when all the computations are complete, AXS should be
// subtracted from the accumulator A.  The addition table is a variable
// array of length 2*numbits (where numbits is the length of the FieldT
// size) such that Ptable[2*i] and Ptable[2*i+1] are the (x,y)
// coordinates of 2^i * P + C.
template<typename FieldT>
class ec_scalarmul_accum_gadget : public gadget<FieldT> {
private:
  pb_variable_array<FieldT> svec;
  std::vector<packing_gadget<FieldT> > packers;
  std::vector<ec_scalarmul_vec_accum_gadget<FieldT> > vecgadget;

public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable<FieldT> Ax, Ay;
  const pb_variable<FieldT> s;
  const pb_variable_array<FieldT> Ptable;

  ec_scalarmul_accum_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &Ax,
              const pb_variable<FieldT> &Ay,
              const pb_variable<FieldT> &s,
              const pb_variable_array<FieldT> &Ptable,
              FieldT &AXSx, FieldT &AXSy) :
    gadget<FieldT>(pb, "ec_scalarmul_accum_gadget"),
    outx(outx), outy(outy), Ax(Ax), Ay(Ay), s(s), Ptable(Ptable)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    size_t numbits = FieldT::num_bits;
    svec.allocate(this->pb, numbits, "svec");
    packers.emplace_back(this->pb, svec, s);
    vecgadget.emplace_back(this->pb, outx, outy, Ax, Ay, svec, Ptable, AXSx, AXSy);
  }

  void generate_r1cs_constraints()
  {
    packers[0].generate_r1cs_constraints(true);
    vecgadget[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    packers[0].generate_r1cs_witness_from_packed();
    vecgadget[0].generate_r1cs_witness();
  }
};

// Compute s*P as (outx, outy) for a precomputed addition table Ptable
// for a variable point P, and s given as a bit vector.  The _caller_ is
// responsible for proving that the elements of svec are bits.
// The addition table is a variable array of length 2*numbits (where
// numbits is the length of svec) such that Ptable[2*i] and
// Ptable[2*i+1] are the (x,y) coordinates of 2^i * P + C.
template<typename FieldT>
class ec_scalarmul_vec_gadget : public gadget<FieldT> {
private:
  FieldT Cx, Cy, Ax, Ay, AXSx, AXSy;
  pb_variable<FieldT> accinx, acciny, accoutx, accouty;
  std::vector<ec_scalarmul_vec_accum_gadget<FieldT> > scalarmuls;
  std::vector<ec_constant_add_gadget<FieldT> > adders;
public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable_array<FieldT> svec;
  const pb_variable_array<FieldT> Ptable;

  ec_scalarmul_vec_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable_array<FieldT> &svec,
              const pb_variable_array<FieldT> &Ptable) :
    gadget<FieldT>(pb, "ec_scalarmul_vec_gadget"),
    // Precomputed coordinates of C and A
    Cx(2),
    Cy("4950745124018817972378217179409499695353526031437053848725554590521829916331"),
    Ax("7536839002660211356286040193441766649532044555061394833845553337792579131020"),
    Ay("11391058648720923807988142436733355540810929560298907319389650598553246451302"),
    outx(outx), outy(outy), svec(svec), Ptable(Ptable)
  {
    AXSx = Ax;
    AXSy = Ay;
    accinx.allocate(this->pb, "accinx");
    acciny.allocate(this->pb, "acciny");
    accoutx.allocate(this->pb, "accoutx");
    accouty.allocate(this->pb, "accouty");

    scalarmuls.emplace_back(pb, accoutx, accouty, accinx, acciny, svec, Ptable, AXSx, AXSy);
    adders.emplace_back(pb, outx, outy, accoutx, accouty, AXSx, -AXSy);
  }

  void generate_r1cs_constraints()
  {
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(accinx, 1, Ax));
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(acciny, 1, Ay));
    scalarmuls[0].generate_r1cs_constraints();
    adders[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    this->pb.val(accinx) = Ax;
    this->pb.val(acciny) = Ay;
    scalarmuls[0].generate_r1cs_witness();
    adders[0].generate_r1cs_witness();
  }
};

// Compute s*P as (outx, outy) for a precomputed addition table Ptable
// for a variable point P, and s given as a field element.  The addition
// table is a variable array of length 2*numbits (where numbits is the
// length of the FieldT size) such that Ptable[2*i] and Ptable[2*i+1]
// are the (x,y) coordinates of 2^i * P + C.
template<typename FieldT>
class ec_scalarmul_gadget : public gadget<FieldT> {
private:
  pb_variable_array<FieldT> svec;
  std::vector<packing_gadget<FieldT> > packers;
  std::vector<ec_constant_scalarmul_vec_gadget<FieldT> > vecgadget;

public:
  const pb_variable<FieldT> outx, outy;
  const pb_variable<FieldT> s;
  const pb_variable_array<FieldT> Ptable;

  ec_scalarmul_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &s,
              const pb_variable_array<FieldT> &Ptable) :
    gadget<FieldT>(pb, "ec_scalarmul_gadget"),
    outx(outx), outy(outy), s(s), Ptable(Ptable)
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    size_t numbits = FieldT::num_bits;
    svec.allocate(this->pb, numbits, "svec");
    packers.emplace_back(this->pb, svec, s);
    vecgadget.emplace_back(this->pb, outx, outy, svec, Ptable);
  }

  void generate_r1cs_constraints()
  {
    packers[0].generate_r1cs_constraints(true);
    vecgadget[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    packers[0].generate_r1cs_witness_from_packed();
    vecgadget[0].generate_r1cs_witness();
  }
};

// Compute a*G + b*H as (outx, outy), given a and b as field elements.
template<typename FieldT>
class ec_pedersen_gadget : public gadget<FieldT> {
private:
  pb_variable<FieldT> accinx, acciny, accmidx, accmidy, accoutx, accouty;
  std::vector<ec_constant_scalarmul_accum_gadget<FieldT> > mulgadgets;
  std::vector<ec_constant_add_gadget<FieldT> > addgadget;
  const FieldT Gx, Gy, Hx, Hy, Ax, Ay;

public:
  const pb_variable<FieldT> outx, outy, a, b;

  ec_pedersen_gadget(protoboard<FieldT> &pb,
              const pb_variable<FieldT> &outx,
              const pb_variable<FieldT> &outy,
              const pb_variable<FieldT> &a,
              const pb_variable<FieldT> &b) :
    gadget<FieldT>(pb, "ec_pedersen_gadget"),
    outx(outx), outy(outy), a(a), b(b),
  // Precomputed coordinates of G, H, and A
  Gx(0),
  Gy("11977228949870389393715360594190192321220966033310912010610740966317727761886"),
  Hx(1),
  Hy("21803877843449984883423225223478944275188924769286999517937427649571474907279"),
  Ax("7536839002660211356286040193441766649532044555061394833845553337792579131020"),
  Ay("11391058648720923807988142436733355540810929560298907319389650598553246451302")
  {
    // Allocate variables to protoboard
    // The strings (like "x") are only for debugging purposes

    accinx.allocate(this->pb, "accinx");
    acciny.allocate(this->pb, "acciny");
    accmidx.allocate(this->pb, "accmidx");
    accmidy.allocate(this->pb, "accmidy");
    accoutx.allocate(this->pb, "accoutx");
    accouty.allocate(this->pb, "accouty");

    // Initialize the accumulator
    FieldT AXSx = Ax;
    FieldT AXSy = Ay;
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(accinx, 1, Ax));
    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(acciny, 1, Ay));

    // Initialize the gadgets
    mulgadgets.emplace_back(this->pb, accmidx, accmidy, accinx, acciny, a, Gx, Gy, AXSx, AXSy);
    mulgadgets.emplace_back(this->pb, accoutx, accouty, accmidx, accmidy, b, Hx, Hy, AXSx, AXSy);
    // Subtract the accumulator excess to get the result
    addgadget.emplace_back(this->pb, outx, outy, accoutx, accouty, AXSx, -AXSy);
  }

  void generate_r1cs_constraints()
  {
    this->pb.val(accinx) = Ax;
    this->pb.val(acciny) = Ay;
    mulgadgets[0].generate_r1cs_constraints();
    mulgadgets[1].generate_r1cs_constraints();
    addgadget[0].generate_r1cs_constraints();
  }

  void generate_r1cs_witness()
  {
    mulgadgets[0].generate_r1cs_witness();
    mulgadgets[1].generate_r1cs_witness();
    addgadget[0].generate_r1cs_witness();
  }
};