Browse Source

Initial commit of pedersen-libsnark code

Using libsnark, create a zkSNARK for knowledge of a preimage
of a Pedersen commitment.
Ian Goldberg 4 years ago
commit
7ea1a3d913
4 changed files with 480 additions and 0 deletions
  1. 339 0
      ecgadget.hpp
  2. 86 0
      pedersen.cpp
  3. 31 0
      sage/findcurve
  4. 24 0
      sage/testcurve

+ 339 - 0
ecgadget.hpp

@@ -0,0 +1,339 @@
+#include "libsnark/gadgetlib1/gadgets/basic_gadgets.hpp"
+
+using namespace libsnark;
+
+// Double the 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_variable<FieldT> &inx,
+              const pb_variable<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.val(inx) * this->pb.val(inx);
+    this->pb.val(lambda) = (this->pb.val(inxsq) * 3 - 3) * (this->pb.val(iny) * 2).inverse();
+    this->pb.val(outx) = this->pb.val(lambda).squared() - this->pb.val(inx) * 2;
+    this->pb.val(outy) = this->pb.val(lambda) * (this->pb.val(inx) - this->pb.val(outx)) - this->pb.val(iny);
+  }
+};
+
+// Add nothing, G, H, or G+H 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> {
+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;
+
+  ec_add_GH_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"),
+    Gx(0),
+    Gy("11977228949870389393715360594190192321220966033310912010610740966317727761886"),
+    Hx(1),
+    Hy("21803877843449984883423225223478944275188924769286999517937427649571474907279"),
+    GHx("2864090850787705444524344020850508438903451433901276387624248428140647539638"),
+    GHy("3350168998338968221269367365107720885864670493693161027931048546881356285970"),
+    outx(outx), outy(outy), inx(inx), iny(iny), addG(addG), addH(addH)
+  {
+    // 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 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)
+
+    // In particular:
+    // G = (0, 11977228949870389393715360594190192321220966033310912010610740966317727761886)
+    // H = (1, 21803877843449984883423225223478944275188924769286999517937427649571474907279)
+    // G+H = (2864090850787705444524344020850508438903451433901276387624248428140647539638, 3350168998338968221269367365107720885864670493693161027931048546881356285970)
+    // so the point to add is ( (GHx - Hx) * addG + (GHx - Gx) * addH + (Gx + Hx - GHx), (GHy - Hy) * addG + (GHy - Gy) * addH + (Gy + Hy - GHy))
+    // = (2864090850787705444524344020850508438903451433901276387624248428140647539637 * addG + 2864090850787705444524344020850508438903451433901276387624248428140647539638 * addH - 2864090850787705444524344020850508438903451433901276387624248428140647539637, -18453708845111016662153857858371223389324254275593838490006379102690118621309 * addG - 8627059951531421172445993229082471435356295539617750982679692419436371475916 * addH + 30430937794981406055869218452561415710545220308904750500617120069007846383195)
+
+    // (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));
+
+    // 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));
+
+    // 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
+
+    // so we compute move = addG || addH, 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>(sumx - inx, move, outx - inx));
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(sumy - iny, move, outy - iny));
+
+  }
+
+  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);
+  }
+};
+
+// Compute a*G + b*H as (outx, outy), given a and b as bit vectors.
+// a and b must be of the same size.  The _caller_ is responsible for
+// proving that the elements of avec and bvec are bits.
+template<typename FieldT>
+class ec_pedersen_vec_gadget : public gadget<FieldT> {
+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;
+  const FieldT Cx, Cy, CGx, CGy, CHx, CHy, CGHx, CGHy;
+  FieldT m2nCx, m2nCy;
+
+  // Compute m2nC = -2^n * C.  We can precomute the answers for values
+  // of n we expect to get.
+  void compute_m2nC(FieldT &m2nCx, FieldT &m2nCy, size_t n)
+  {
+    if (n == 253) {
+        m2nCx = FieldT("2630025903576807331238993847875694711243784786568881628418508626984487096258");
+        m2nCy = FieldT("17628834417659968531880949658739649785248429713924280788649629869316127047701");
+    } else {
+        // Invariant: m2iC = -2^i * C
+        FieldT m2iCx = Cx;
+        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;
+            m2iCx = m2iCxo;
+            m2iCy = m2iCyo;
+            ++i;
+        }
+        m2nCx = m2iCx;
+        m2nCy = m2iCy;
+    }
+  }
+
+public:
+  const pb_variable<FieldT> outx, outy;
+  const pb_variable_array<FieldT> avec, bvec;
+
+  ec_pedersen_vec_gadget(protoboard<FieldT> &pb,
+              const pb_variable<FieldT> &outx,
+              const pb_variable<FieldT> &outy,
+              const pb_variable_array<FieldT> &avec,
+              const pb_variable_array<FieldT> &bvec) :
+    gadget<FieldT>(pb, "ec_pedersen_vec_gadget"),
+    Cx(2),
+    Cy("4950745124018817972378217179409499695353526031437053848725554590521829916331"),
+    CGx("4998993376791159436553350546778310121346937620672073819457843493128326049156"),
+    CGy("11119675827304465476900978353730540420130346377889406728458325551400357147144"),
+    CHx("19614539896004018833724771305328960655474424364705508053472946746883341111010"),
+    CHy("9853241351900213537247225242092949438866383394579783148395572971112906592855"),
+    CGHx("10755582242294898568680134375159803731902153202607833320871336755950640390928"),
+    CGHy("3110667473759844579409644567672992116704859238881299917617768683686288881761"),
+    outx(outx), outy(outy), avec(avec), bvec(bvec)
+  {
+    // Allocate variables to protoboard
+    // The strings (like "x") are only for debugging purposes
+	  
+    size_t numbits = avec.size();
+    accumx.allocate(this->pb, numbits, "accumx");
+    accumy.allocate(this->pb, numbits, "accumy");
+    daccumx.allocate(this->pb, numbits-1, "daccumx");
+    daccumy.allocate(this->pb, numbits-1, "daccumy");
+    lambda.allocate(this->pb, "lambda");
+
+    for (size_t i = 0; i < numbits-1; ++i) {
+        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]);
+    }
+
+    compute_m2nC(m2nCx, m2nCy, numbits-1);
+  }
+
+  void generate_r1cs_constraints()
+  {
+    // Strategy: We do a basic double-and-add, using the variant of a
+    // single double and an add of one of (O, G, H, G+H) at each step.
+    // *However*, there's a twist.  Our doubling and adding routines
+    // don't handle the point at infinity O as the point to add to or to
+    // double.  (Adding O as one of the four options above is fine.)
+    // Normally, the double-and-add algorithm starts with an accumulator
+    // of O, and that won't work for us.  So instead, we start the
+    // accumulator at a different base point C, whose discrete log
+    // with respect to the (G,H) basis is unknown.  Then we'll end up
+    // with an extra 2^n * C in the accumulator (where n is the number
+    // of doublings we do), so at the end, we'll add the (constant!)
+    // point -2^n * C to get the final result.  That the discrete log of
+    // C is unknown means we won't encounter O along the way, either (if
+    // we did, we could compute the DL of C in the (G,H) basis).
+
+    // For the first bit, we just precompute C, C+G, C+H, C+G+H and use
+    // the top bit of a and b to choose which one to start with.
+
+    // accumx[0] = Cx + (CGx - Cx) * avec[numbits-1] + (CHx - Cx) *
+    // bvec[numbits-1] + (CGHx - CGx - CHx + Cx) *
+    // avec[numbits-1]*bvec[numbits-1] (and similarly for y)
+
+    // We could possibly optimize this later by computing the a*b
+    // product once, but then we'd have to pass a large linear
+    // combination to a gadget, which it probably doesn't like?  It
+    // would save just one constraint, so probably not so important?
+    size_t numbits = avec.size();
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>((CGHx - CGx - CHx + Cx) * avec[numbits-1], bvec[numbits-1], accumx[0] - (Cx + (CGx - Cx) * avec[numbits-1] + (CHx - Cx) * bvec[numbits-1])));
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>((CGHy - CGy - CHy + Cy) * avec[numbits-1], bvec[numbits-1], accumy[0] - (Cy + (CGy - Cy) * avec[numbits-1] + (CHy - Cy) * bvec[numbits-1])));
+
+    // After that, for each remaining bit of a and b, we use the
+    // ec_double_gadget and the ec_add_GH_gadget to accumulate
+    // the answer.
+    for (size_t i = 0; i < numbits-1; ++i) {
+        doublers[i].generate_r1cs_constraints();
+        adders[i].generate_r1cs_constraints();
+    }
+
+    // Finally, we add the constant point -2^n * C to the result, where
+    // n = numbits-1
+    // (m2nCx - accumx[numbits-1]) * lambda = m2nCy - accumy[numbits-1]
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(m2nCx - accumx[numbits-1], lambda, m2nCy - accumy[numbits-1]));
+    // outx = lambda^2 - (m2nCx + accumx[numbits-1])
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, lambda, outx + m2nCx + accumx[numbits-1]));
+    // outy = lambda * (accumx[numbits-1] - outx) - accumy[numbits-1]
+    this->pb.add_r1cs_constraint(r1cs_constraint<FieldT>(lambda, accumx[numbits-1] - outx, outy + accumy[numbits-1]));
+
+  }
+
+  void generate_r1cs_witness()
+  {
+    size_t numbits = avec.size();
+    this->pb.val(accumx[0]) = Cx + (CGx - Cx) * this->pb.val(avec[numbits-1]) + (CHx - Cx) * this->pb.val(bvec[numbits-1]) + (CGHx - CGx - CHx + Cx) * this->pb.val(avec[numbits-1])*this->pb.val(bvec[numbits-1]);
+    this->pb.val(accumy[0]) = Cy + (CGy - Cy) * this->pb.val(avec[numbits-1]) + (CHy - Cy) * this->pb.val(bvec[numbits-1]) + (CGHy - CGy - CHy + Cy) * this->pb.val(avec[numbits-1])*this->pb.val(bvec[numbits-1]);
+
+    for (size_t i = 0; i < numbits-1; ++i) {
+        doublers[i].generate_r1cs_witness();
+        adders[i].generate_r1cs_witness();
+    }
+
+    this->pb.val(lambda) = (m2nCy - this->pb.val(accumy[numbits-1])) *
+                           (m2nCx - this->pb.val(accumx[numbits-1])).inverse();
+    this->pb.val(outx) = this->pb.val(lambda).squared() -
+                         (m2nCx + this->pb.val(accumx[numbits-1]));
+    this->pb.val(outy) = this->pb.val(lambda) *
+                         (this->pb.val(accumx[numbits-1]) - this->pb.val(outx))
+                         - this->pb.val(accumy[numbits-1]);
+    
+  }
+};
+
+// 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_array<FieldT> avec, bvec;
+  std::vector<packing_gadget<FieldT> > packers;
+  std::vector<ec_pedersen_vec_gadget<FieldT> > vecgadget;
+
+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)
+  {
+    // Allocate variables to protoboard
+    // The strings (like "x") are only for debugging purposes
+	  
+    FieldT minus1(-1);
+    size_t numbits = FieldT::num_bits;
+    avec.allocate(this->pb, numbits, "a");
+    bvec.allocate(this->pb, numbits, "b");
+    packers.emplace_back(this->pb, avec, a);
+    packers.emplace_back(this->pb, bvec, b);
+    vecgadget.emplace_back(this->pb, outx, outy, avec, bvec);
+  }
+
+  void generate_r1cs_constraints()
+  {
+    packers[0].generate_r1cs_constraints(true);
+    packers[1].generate_r1cs_constraints(true);
+    vecgadget[0].generate_r1cs_constraints();
+  }
+
+  void generate_r1cs_witness()
+  {
+    packers[0].generate_r1cs_witness_from_packed();
+    packers[1].generate_r1cs_witness_from_packed();
+    vecgadget[0].generate_r1cs_witness();
+  }
+};

+ 86 - 0
pedersen.cpp

@@ -0,0 +1,86 @@
+#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");
+  b.allocate(pb, "b");
+
+  // 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_pedersen_gadget<FieldT> ped(pb, outx, outy, a, b);
+  ped.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();
+  pb.val(b) = FieldT::random_element();
+  cout << "Computing " << pb.val(a) << "*G + " << pb.val(b) << "*H" << endl;
+
+  ped.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_pedersen");
+  pkfile << keypair.pk;
+  pkfile.close();
+  ofstream vkfile("vk_pedersen");
+  vkfile << keypair.vk;
+  vkfile.close();
+  ofstream pffile("proof_pedersen");
+  pffile << proof;
+  pffile.close();
+
+  return 0;
+}

+ 31 - 0
sage/findcurve

@@ -0,0 +1,31 @@
+#!/usr/bin/env sage
+
+import sys
+from sage.all import *
+
+# The prime order of libsnark's BN128 curve
+r = 21888242871839275222246405745257275088548364400416034343698204186575808495617
+F = GF(r)
+
+count = 0
+primecount = 0
+done = False
+while not done:
+    b = F.random_element()
+    try:
+	E = EllipticCurve(F, [-3, b])
+	Ec = E.cardinality()
+	count += 1
+	print "b =", b, "Ec =", Ec, "c =", count, "pc =", primecount
+	Etc = 2 * p + 2 - Ec
+	if Ec in Primes():
+	    primecount += 1
+	    print "Ec is prime; primecount =", primecount
+	    if Etc in Primes():
+		print "Solution found:", b
+		done = True
+    except KeyboardInterrupt:
+	print "Terminating on keyboard interrupt"
+	sys.exit(1)
+    except:
+    	pass

+ 24 - 0
sage/testcurve

@@ -0,0 +1,24 @@
+#!/usr/bin/env sage
+
+import sys
+from sage.all import *
+
+r = 21888242871839275222246405745257275088548364400416034343698204186575808495617
+F = GF(r)
+b = 7950939520449436327800262930799465135910802758673292356620796789196167463969
+btwist = F(b)*125
+print "b      =", b
+print "btwist =", btwist
+
+E = EllipticCurve(F, [-3, b])
+print E
+Ec = E.cardinality()
+print "Ec  =", Ec
+assert(Ec in Primes())
+
+Et = EllipticCurve(F, [-75, btwist])
+print Et
+Etc = Et.cardinality()
+print "Etc =", Etc
+assert(Etc in Primes())
+assert(Ec + Etc == 2*r+2)