Improved algorithm for scalar multiplication

The old algorithm was a single multiexponentiation, using 3045
constraints for a 254-bit modulus.  The new algorithm implements
scalar multiplication of a public base in 4 constraints per bit of
modulus (including 1 for the bit decomposition of the exponent), plus a
5-constraint constant overhead, for 1021 constraints per multiplcation,
and 2045 for the whole Pedersen commitment.

Makefile

@@ -1,4 +1,4 @@
-all: pedersen
+all: pedersen scalarmul
 
 LIBSNARK=libsnark
 
@@ -8,3 +8,6 @@ LDFLAGS=-L$(LIBSNARK)/build -L$(LIBSNARK)/build/libsnark -L$(LIBSNARK)/build/dep
 
 pedersen: pedersen.cpp ecgadget.hpp
 	g++ $(CXXFLAGS) -o pedersen pedersen.cpp $(LDFLAGS)
+
+scalarmul: scalarmul.cpp ecgadget.hpp
+	g++ $(CXXFLAGS) -o scalarmul scalarmul.cpp $(LDFLAGS)

README.md

@@ -2,9 +2,9 @@
 
 *Ian Goldberg (iang@uwaterloo.ca), August 2019*
 
-I spent a day learning how to use [libsnark](https://github.com/scipr-lab/libsnark), and thought an interesting first project would be to create a zkSNARK for knowledge of a preimage for a Pedersen commitment.
+I spent a day learning how to use [libsnark](https://github.com/scipr-lab/libsnark), and thought an interesting first project would be to create a zkSNARK for knowledge of a preimage for a Pedersen commitment.  I spent another day reimplementing it with a better scalar multiplication algorithm.
 
-The algorithm is a pretty naive square-and-multiply; there are surely better ones.  This circuit ends up with 3045 constraints for a Pedersen commitment over a 254-bit elliptic curve.
+The algorithm is still a pretty naive square-and-multiply; there are surely better ones.  This circuit ends up with 2045 constraints for a Pedersen commitment over a 254-bit elliptic curve.
 
 It uses libsnark's BN128 implementation, which has an order (not modulus) of r=21888242871839275222246405745257275088548364400416034343698204186575808495617.  Then using [Sage](http://www.sagemath.org/), the [findcurve](sage/findcurve) script in this repo searches for an elliptic curve with _modulus_ r, and with both prime order and whose twist has prime order.  (You do not need to run the findcurve script yourself.)  The resulting curve (over F_r) is E: y^2 = x^3 - 3*x + b, where b=7950939520449436327800262930799465135910802758673292356620796789196167463969.  The order of this curve is the prime 21888242871839275222246405745257275088760161411100494528458776273921456643749.
 

ecgadget.hpp

@@ -2,6 +2,30 @@
 
 using namespace libsnark;
 
+// Double a public 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 public 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 EC point (inx,iny) to yield (outx,outy).  The input point
 // must not be the point at infinity.
 template<typename FieldT>
@@ -51,33 +75,192 @@ public:
   }
 };
 
-// Add nothing, G, H, or G+H to the EC point (inx,iny) to yield
+// Add the EC point (addx,addy) to the EC point (inx,iny) to yield
 // (outx,outy).  The input point must not be the point at infinity.
-// The two input bits addG and addH control what is added.
 template<typename FieldT>
-class ec_add_GH_gadget : public gadget<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 public EC point P to the EC point (inx,iny) to yield
+// (outx,outy).  The input point must not be the point at infinity.
+template<typename FieldT>
+class ec_public_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_public_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_public_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 nothing or the public EC point P to the EC point (inx,iny) to
+// yield (outx,outy).  The input point must not be the point at
+// infinity.  The input bit do_add controls whether the addition is
+// done.
+template<typename FieldT>
+class ec_conditional_add_gadget : public gadget<FieldT> {
+private:
+  pb_variable<FieldT> sumx, sumy;
+  std::vector<ec_public_add_gadget<FieldT> > adder;
+public:
+  const pb_variable<FieldT> outx, outy;
+  const pb_linear_combination<FieldT> inx, iny;
+  const pb_variable<FieldT> do_add;
+  const FieldT Px, Py;
+
+  ec_conditional_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> &do_add,
+              const FieldT &Px, const FieldT &Py) :
+    gadget<FieldT>(pb, "ec_conditional_add_gadget"),
+    outx(outx), outy(outy), inx(inx), iny(iny), do_add(do_add),
+    Px(Px), Py(Py)
+  {
+    // Allocate variables to protoboard
+    // The strings (like "x") are only for debugging purposes
+	  
+    sumx.allocate(this->pb, "sumx");
+    sumy.allocate(this->pb, "sumy");
+    adder.emplace_back(this->pb, sumx, sumy, inx, iny, Px, Py);
+  }
+
+  void generate_r1cs_constraints()
+  {
+    // Strategy: we always do the addition, but if do_add = 0, we throw
+    // it away later.
+
+    adder[0].generate_r1cs_constraints();
+
+    // Now we want to conditionally move the sum.  We want that
+    // outx = do_add ? sumx : inx
+    // outy = do_add ? sumy : iny
+
+    // so we compute
+    // outx = inx + (sumx - inx) * do_add
+    // outy = iny + (sumy - iny) * do_add
+
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(sumx - inx, do_add, outx - inx));
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(sumy - iny, do_add, outy - iny));
+
+  }
+
+  void generate_r1cs_witness()
+  {
+    adder[0].generate_r1cs_witness();
+
+    bool move = this->pb.val(do_add) != FieldT(0);
+    this->pb.val(outx) = move ? this->pb.val(sumx) : this->pb.lc_val(inx);
+    this->pb.val(outy) = move ? this->pb.val(sumy) : this->pb.lc_val(iny);
+  }
+};
+
+// Add nothing, or one of the public 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 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;
-  const FieldT Gx, Gy, Hx, Hy, GHx, GHy;
 public:
-  const pb_variable<FieldT> outx, outy, inx, iny, addG, addH;
+  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_GH_gadget(protoboard<FieldT> &pb,
+  ec_add_P123_gadget(protoboard<FieldT> &pb,
               const pb_variable<FieldT> &outx,
               const pb_variable<FieldT> &outy,
-              const pb_variable<FieldT> &inx,
-              const pb_variable<FieldT> &iny,
-              const pb_variable<FieldT> &addG,
-              const pb_variable<FieldT> &addH) : 
-    gadget<FieldT>(pb, "ec_add_GH_gadget"),
-    // The coordinates of G, H, and G+H
-    Gx(0),
-    Gy("11977228949870389393715360594190192321220966033310912010610740966317727761886"),
-    Hx(1),
-    Hy("21803877843449984883423225223478944275188924769286999517937427649571474907279"),
-    GHx("2864090850787705444524344020850508438903451433901276387624248428140647539638"),
-    GHy("3350168998338968221269367365107720885864670493693161027931048546881356285970"),
-    outx(outx), outy(outy), inx(inx), iny(iny), addG(addG), addH(addH)
+              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
@@ -90,33 +273,33 @@ public:
 
   void generate_r1cs_constraints()
   {
-    // Strategy: if addG = addH = 0, we compute some nonsense but throw
-    // it away later.  Otherwise, the coordinates of addG * G + addH * H
-    // are a _linear_ function of addG and addH (since G and H are global
-    // constants)
+    // 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 ( (GHx - Hx) * addG + (GHx -
-    // Gx) * addH + (Gx + Hx - GHx), (GHy - Hy) * addG + (GHy - Gy) *
-    // addH + (Gy + Hy - GHy))
+    // 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>((GHx - Hx) * addG + (GHx - Gx) * addH + (Gx + Hx - GHx) - inx, lambda, (GHy - Hy) * addG + (GHy - Gy) * addH + (Gy + Hy - GHy) - 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 + (GHx - Hx) * addG + (GHx - Gx) * addH + (Gx + Hx - GHx) + 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 = (addG || addH) ? sumx : inx
-    // outy = (addG || addH) ? sumy : iny
+    // outx = (add1 || add2) ? sumx : inx
+    // outy = (add1 || add2) ? sumy : iny
 
-    // so we compute move = addG || addH, and then
+    // 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 - addG, 1 - addH, 1 - 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));
 
@@ -124,17 +307,179 @@ public:
 
   void generate_r1cs_witness()
   {
-    FieldT addxval = (GHx - Hx) * this->pb.val(addG) + (GHx - Gx) * this->pb.val(addH) + (Gx + Hx - GHx);
-    FieldT addyval = (GHy - Hy) * this->pb.val(addG) + (GHy - Gy) * this->pb.val(addH) + (Gy + Hy - GHy);
-    this->pb.val(lambda) = (addyval - this->pb.val(iny)) * (addxval - this->pb.val(inx)).inverse();
-    this->pb.val(sumx) = this->pb.val(lambda).squared() - (addxval + this->pb.val(inx));
-    this->pb.val(sumy) = this->pb.val(lambda) * (this->pb.val(inx) - this->pb.val(sumx)) - this->pb.val(iny);
-
-    bool aG = this->pb.val(addG) != FieldT(0);
-    bool aH = this->pb.val(addH) != FieldT(0);
-    this->pb.val(move) = aG || aH;
-    this->pb.val(outx) = this->pb.val(aG || aH ? sumx : inx);
-    this->pb.val(outy) = this->pb.val(aG || aH ? sumy : iny);
+    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);
+  }
+};
+
+// Compute a*P as (outx, outy) for a given public point P, given a
+// as a bit vector.  The _caller_ is responsible for proving that the
+// elements of avec are bits.
+template<typename FieldT>
+class ec_scalarmul_vec_gadget : public gadget<FieldT> {
+private:
+  FieldT Cx, Cy, CPx, CPy;
+  pb_variable_array<FieldT> accumx, accumy;
+  std::vector<ec_add_P123_gadget<FieldT> > conddoubleadders;
+  std::vector<ec_conditional_add_gadget<FieldT> > condsingleadders;
+  std::vector<ec_public_add_gadget<FieldT> > singleadders;
+public:
+  const pb_variable<FieldT> outx, outy;
+  const pb_variable_array<FieldT> avec;
+  const FieldT Px, Py;
+
+  // Strategy: We compute (in public) (powers of 2) times P, and then
+  // conditionally add them into an accumulator.  Because our adder cannot
+  // handle the point at infinity O, we start the accumulator with a value
+  // of C, whose discrete log with respect to P should be unknown, so that
+  // we won't encounter O along the way.  (Also, a should not be 0 or
+  // the group order.)  We actually start the accumulator with either
+  // C or C+P depending on avec[0], so we get the first conditional add
+  // "for free".  Then we use the ec_add_P123_gadget to do the conditional
+  // adds two at a time (at a cost of 6 constraints per pair, as opposed to
+  // 5 for a single conditional add).  If the length of avec is even, then
+  // there will be one left over, and we do a single conditional add for
+  // that one.  Finally, we add the public point -C.
+
+  ec_scalarmul_vec_gadget(protoboard<FieldT> &pb,
+              const pb_variable<FieldT> &outx,
+              const pb_variable<FieldT> &outy,
+              const pb_variable_array<FieldT> &avec,
+              const FieldT &Px, const FieldT &Py) :
+    gadget<FieldT>(pb, "ec_pedersen_vec_gadget"),
+    // Precomputed coordinates of C
+    Cx(2),
+    Cy("4950745124018817972378217179409499695353526031437053848725554590521829916331"),
+    outx(outx), outy(outy), avec(avec), Px(Px), Py(Py)
+  {
+    size_t numbits = avec.size();
+    accumx.allocate(this->pb, numbits/2+1, "accumx");
+    accumy.allocate(this->pb, numbits/2+1, "accumy");
+
+    ec_add_points(CPx, CPy, Cx, Cy, Px, Py);
+
+    FieldT twoiPx, twoiPy, twoi1Px, twoi1Py, twoi3Px, twoi3Py;
+    size_t i = 1;
+
+    ec_double_point(twoiPx, twoiPy, Px, Py);
+
+    while(i < numbits) {
+        // Invariants: i is odd, and twoiP = 2^i * P
+        // Compute twoi1P = 2^{i+1} * P = 2 * twoiP and
+        //         twoi3P = 2^i * 3 * P = 3 * twoiP
+        ec_double_point(twoi1Px, twoi1Py, twoiPx, twoiPy);
+        ec_add_points(twoi3Px, twoi3Py, twoi1Px, twoi1Py, twoiPx, twoiPy);
+
+        if (i == numbits-1) {
+            // There's only one bit of avec left; use a single conditional
+            // add.
+            condsingleadders.emplace_back(this->pb,
+                accumx[(i+1)/2], accumy[(i+1)/2],
+                accumx[(i-1)/2], accumy[(i-1)/2],
+                avec[i],
+                twoiPx, twoiPy);
+        } else {
+            conddoubleadders.emplace_back(this->pb,
+                accumx[(i+1)/2], accumy[(i+1)/2],
+                accumx[(i-1)/2], accumy[(i-1)/2],
+                avec[i], avec[i+1],
+                twoiPx, twoiPy, twoi1Px, twoi1Py, twoi3Px, twoi3Py);
+        }
+
+        ec_double_point(twoiPx, twoiPy, twoi1Px, twoi1Py);
+        i += 2;
+    }
+
+    // If numbits is even, the output so far is in accum[(numbits)/2].
+    // If numbits is odd, it is in accum[(numbits-1)/2].  So in either
+    // case, it is in accum[numbits/2].
+    singleadders.emplace_back(this->pb,
+        outx, outy, accumx[numbits/2], accumy[numbits/2],
+        Cx, -Cy);
+  }
+
+  void generate_r1cs_constraints()
+  {
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(Cx + (CPx-Cx) * avec[0], 1, accumx[0]));
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(Cy + (CPy-Cy) * avec[0], 1, accumy[0]));
+
+    for (auto&& gadget : conddoubleadders) {
+        gadget.generate_r1cs_constraints();
+    }
+    for (auto&& gadget : condsingleadders) {
+        gadget.generate_r1cs_constraints();
+    }
+    for (auto&& gadget : singleadders) {
+        gadget.generate_r1cs_constraints();
+    }
+  }
+
+  void generate_r1cs_witness()
+  {
+    this->pb.val(accumx[0]) = Cx + (CPx-Cx) * this->pb.val(avec[0]);
+    this->pb.val(accumy[0]) = Cy + (CPy-Cy) * this->pb.val(avec[0]);
+
+    for (auto&& gadget : conddoubleadders) {
+        gadget.generate_r1cs_witness();
+    }
+    for (auto&& gadget : condsingleadders) {
+        gadget.generate_r1cs_witness();
+    }
+    for (auto&& gadget : singleadders) {
+        gadget.generate_r1cs_witness();
+    }
+  }
+};
+
+// Compute a*P as (outx, outy) for a given public point P, given a
+// as a field element.
+template<typename FieldT>
+class ec_scalarmul_gadget : public gadget<FieldT> {
+private:
+  pb_variable_array<FieldT> avec;
+  std::vector<packing_gadget<FieldT> > packers;
+  std::vector<ec_scalarmul_vec_gadget<FieldT> > vecgadget;
+
+public:
+  const pb_variable<FieldT> outx, outy, a;
+  const FieldT Px, Py;
+
+  ec_scalarmul_gadget(protoboard<FieldT> &pb,
+              const pb_variable<FieldT> &outx,
+              const pb_variable<FieldT> &outy,
+              const pb_variable<FieldT> &a,
+              const FieldT &Px, const FieldT &Py) :
+    gadget<FieldT>(pb, "ec_scalarmul_gadget"),
+    outx(outx), outy(outy), a(a), Px(Px), Py(Py)
+  {
+    // Allocate variables to protoboard
+    // The strings (like "x") are only for debugging purposes
+	  
+    size_t numbits = FieldT::num_bits;
+    avec.allocate(this->pb, numbits, "avec");
+    packers.emplace_back(this->pb, avec, a);
+    vecgadget.emplace_back(this->pb, outx, outy, avec, 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();
   }
 };
 
@@ -147,7 +492,7 @@ private:
   pb_variable_array<FieldT> accumx, accumy, daccumx, daccumy;
   pb_variable<FieldT> lambda;
   std::vector<ec_double_gadget<FieldT> > doublers;
-  std::vector<ec_add_GH_gadget<FieldT> > adders;
+  std::vector<ec_add_P123_gadget<FieldT> > adders;
   const FieldT Cx, Cy, CGx, CGy, CHx, CHy, CGHx, CGHy;
   FieldT m2nCx, m2nCy;
 
@@ -165,10 +510,8 @@ private:
         FieldT m2iCy = -Cy;
         size_t i = 0;
         while (i < n) {
-            FieldT xsq = m2iCx.squared();
-            FieldT lambda = (xsq * 3 - 3) * (m2iCy * 2).inverse();
-            FieldT m2iCxo = lambda.squared() - m2iCx * 2;
-            FieldT m2iCyo = lambda * (m2iCx - m2iCxo) - m2iCy;
+            FieldT m2iCxo, m2iCyo;
+            ec_double_point(m2iCxo, m2iCyo, m2iCx, m2iCy);
             m2iCx = m2iCxo;
             m2iCy = m2iCyo;
             ++i;
@@ -213,7 +556,12 @@ public:
         doublers.emplace_back(this->pb, daccumx[i], daccumy[i],
             accumx[i], accumy[i]);
         adders.emplace_back(this->pb, accumx[i+1], accumy[i+1],
-            daccumx[i], daccumy[i], avec[numbits-2-i], bvec[numbits-2-i]);
+            daccumx[i], daccumy[i], avec[numbits-2-i], bvec[numbits-2-i],
+            // The coordinates of G, H, and G+H
+            FieldT(0), FieldT("11977228949870389393715360594190192321220966033310912010610740966317727761886"),
+            FieldT(1), FieldT("21803877843449984883423225223478944275188924769286999517937427649571474907279"),
+            FieldT("2864090850787705444524344020850508438903451433901276387624248428140647539638"),
+            FieldT("3350168998338968221269367365107720885864670493693161027931048546881356285970"));
     }
 
     compute_m2nC(m2nCx, m2nCy, numbits-1);
@@ -295,7 +643,7 @@ public:
 
 // 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> {
+class ec_old_pedersen_gadget : public gadget<FieldT> {
 private:
   pb_variable_array<FieldT> avec, bvec;
   std::vector<packing_gadget<FieldT> > packers;
@@ -304,20 +652,20 @@ private:
 public:
   const pb_variable<FieldT> outx, outy, a, b;
 
-  ec_pedersen_gadget(protoboard<FieldT> &pb,
+  ec_old_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"),
+    gadget<FieldT>(pb, "ec_old_pedersen_gadget"),
     outx(outx), outy(outy), a(a), b(b)
   {
     // Allocate variables to protoboard
     // The strings (like "x") are only for debugging purposes
 	  
     size_t numbits = FieldT::num_bits;
-    avec.allocate(this->pb, numbits, "a");
-    bvec.allocate(this->pb, numbits, "b");
+    avec.allocate(this->pb, numbits, "avec");
+    bvec.allocate(this->pb, numbits, "bvec");
     packers.emplace_back(this->pb, avec, a);
     packers.emplace_back(this->pb, bvec, b);
     vecgadget.emplace_back(this->pb, outx, outy, avec, bvec);
@@ -337,3 +685,55 @@ public:
     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> aoutx, aouty, boutx, bouty;
+  std::vector<ec_scalarmul_gadget<FieldT> > mulgadgets;
+  std::vector<ec_add_gadget<FieldT> > addgadget;
+  const FieldT Gx, Gy, Hx, Hy;
+
+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 and H
+  Gx(0),
+  Gy("11977228949870389393715360594190192321220966033310912010610740966317727761886"),
+  Hx(1),
+  Hy("21803877843449984883423225223478944275188924769286999517937427649571474907279")
+  {
+    // Allocate variables to protoboard
+    // The strings (like "x") are only for debugging purposes
+	  
+    aoutx.allocate(this->pb, "aoutx");
+    aouty.allocate(this->pb, "aouty");
+    boutx.allocate(this->pb, "boutx");
+    bouty.allocate(this->pb, "bouty");
+    mulgadgets.emplace_back(this->pb, aoutx, aouty, a, Gx, Gy);
+    mulgadgets.emplace_back(this->pb, boutx, bouty, b, Hx, Hy);
+    addgadget.emplace_back(this->pb, outx, outy, aoutx, aouty, boutx, bouty);
+  }
+
+  void generate_r1cs_constraints()
+  {
+    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();
+  }
+};

scalarmul.cpp

@@ -0,0 +1,84 @@
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+
+#include "libff/algebra/fields/field_utils.hpp"
+#include "libsnark/zk_proof_systems/ppzksnark/r1cs_ppzksnark/r1cs_ppzksnark.hpp"
+#include "libsnark/common/default_types/r1cs_ppzksnark_pp.hpp"
+#include "libsnark/gadgetlib1/pb_variable.hpp"
+
+#include "ecgadget.hpp"
+
+using namespace libsnark;
+using namespace std;
+
+int main()
+{
+  // Initialize the curve parameters
+
+  default_r1cs_ppzksnark_pp::init_public_params();
+
+  typedef libff::Fr<default_r1cs_ppzksnark_pp> FieldT;
+  
+  // Create protoboard
+
+  libff::start_profiling();
+
+  cout << "Keypair" << endl;
+
+  protoboard<FieldT> pb;
+  pb_variable<FieldT> outx, outy;
+  pb_variable<FieldT> a, b;
+
+  // Allocate variables
+
+  outx.allocate(pb, "outx");
+  outy.allocate(pb, "outy");
+  a.allocate(pb, "a");
+
+  // This sets up the protoboard variables so that the first n of them
+  // represent the public input and the rest is private input
+
+  pb.set_input_sizes(2);
+
+  // Initialize gadget
+
+  ec_scalarmul_gadget<FieldT> sm(pb, outx, outy, a, FieldT(0), FieldT("11977228949870389393715360594190192321220966033310912010610740966317727761886"));
+  sm.generate_r1cs_constraints();
+  
+  const r1cs_constraint_system<FieldT> constraint_system = pb.get_constraint_system();
+
+  const r1cs_ppzksnark_keypair<default_r1cs_ppzksnark_pp> keypair = r1cs_ppzksnark_generator<default_r1cs_ppzksnark_pp>(constraint_system);
+
+  // Add witness values
+
+  cout << "Prover" << endl;
+  
+  pb.val(a) = FieldT::random_element();
+  cout << "Computing " << pb.val(a) << "*G" << endl;
+
+  sm.generate_r1cs_witness();
+
+  const r1cs_ppzksnark_proof<default_r1cs_ppzksnark_pp> proof = r1cs_ppzksnark_prover<default_r1cs_ppzksnark_pp>(keypair.pk, pb.primary_input(), pb.auxiliary_input());
+
+  cout << "Verifier" << endl;
+
+  bool verified = r1cs_ppzksnark_verifier_strong_IC<default_r1cs_ppzksnark_pp>(keypair.vk, pb.primary_input(), proof);
+
+  cout << "Number of R1CS constraints: " << constraint_system.num_constraints() << endl;
+  cout << "Primary (public) input: " << pb.primary_input() << endl;
+  cout << "Auxiliary (private) input: " << pb.auxiliary_input() << endl;
+  cout << "Verification status: " << verified << endl;
+
+  ofstream pkfile("pk_scalarmul");
+  pkfile << keypair.pk;
+  pkfile.close();
+  ofstream vkfile("vk_scalarmul");
+  vkfile << keypair.vk;
+  vkfile.close();
+  ofstream pffile("proof_scalarmul");
+  pffile << proof;
+  pffile.close();
+
+  return 0;
+}