Add faster FpMul implementation

circuits/bigint.circom

@@ -464,12 +464,10 @@ template BigSubModP(n, k){
     component add = BigAdd(n, k);
     for (var i = 0; i < k; i++){
         add.a[i] <== sub.out[i];
-        add.b[i] <== p[i];
+        add.b[i] <== flag * p[i];
     }
-    signal tmp[k];
     for (var i = 0; i < k; i++){
-        tmp[i] <== (1 - flag) * sub.out[i];
-        out[i] <== tmp[i] + flag * add.out[i];
+        out[i] <== add.out[i];
     }
 }
 
@@ -540,6 +538,8 @@ template CheckCarryToZero(n, m, k) {
 
     var EPSILON = 3;
 
+    assert(m + EPSILON <= 253);
+
     signal input in[k];
 
     signal carry[k];

circuits/bigint_func.circom

@@ -34,25 +34,6 @@ function SplitThreeFn(in, n, m, k) {
     return [in % (1 << n), (in \ (1 << n)) % (1 << m), (in \ (1 << n + m)) % (1 << k)];
 }
 
-// in is an m bit number
-// split into ceil(m/n) n-bit registers
-function splitOverflowedRegister(m, n, in) {
-    var out[100];
-
-    for (var i = 0; i < 100; i++) {
-        out[i] = 0;
-    }
-
-    var nRegisters = div_ceil(m, n);
-    var running = in;
-    for (var i = 0; i < nRegisters; i++) {
-        out[i] = running % (1<<n);
-        running>>=n;
-    }
-
-    return out;
-}
-
 // m bits per overflowed register (values are potentially negative)
 // n bits per properly-sized register
 // in has k registers
@@ -61,60 +42,73 @@ function splitOverflowedRegister(m, n, in) {
 // - 1 since the last register is included in the last ceil(m/n) array
 // + 1 since the carries from previous registers could push you over
 function getProperRepresentation(m, n, k, in) {
-    var ceilMN = 0; // ceil(m/n)
-    if (m % n == 0) {
-        ceilMN = m \ n;
-    } else {
-        ceilMN = m \ n + 1;
-    }
+    var ceilMN = div_ceil(m, n);
 
-    var pieces[100][100]; // should be pieces[k][ceilMN]
+    var out[100]; // should be out[k + ceilMN]
+    assert(k + ceilMN < 100);
     for (var i = 0; i < k; i++) {
-        for (var j = 0; j < 100; j++) {
-            pieces[i][j] = 0;
-        }
-        if (isNegative(in[i]) == 1) {
-            var negPieces[100] = splitOverflowedRegister(m, n, -1 * in[i]);
-            for (var j = 0; j < ceilMN; j++) {
-                pieces[i][j] = -1 * negPieces[j];
-            }
-        } else {
-            pieces[i] = splitOverflowedRegister(m, n, in[i]);
-        }
+        out[i] = in[i];
+    }
+    for (var i = k; i < 100; i++) {
+        out[i] = 0;
+    }
+    assert(n <= m);
+    for (var i = 0; i+1 < k + ceilMN; i++) {
+        assert((1 << m) >= out[i] && out[i] >= -(1 << m));
+        var shifted_val = out[i] + (1 << m);
+        assert(0 <= shifted_val && shifted_val <= (1 << (m+1)));
+        out[i] = shifted_val & ((1 << n) - 1);
+        out[i+1] += (shifted_val >> n) - (1 << (m - n));
     }
 
-    var out[100]; // should be out[k + ceilMN]
-    var carries[100]; // should be carries[k + ceilMN]
-    for (var i = 0; i < 100; i++) {
+    return out;
+}
+
+// Evaluate polynomial a at point x
+function poly_eval(len, a, x) {
+    var v = 0;
+    for (var i = 0; i < len; i++) {
+        v += a[i] * (x ** i);
+    }
+    return v;
+}
+
+// Interpolate a degree len-1 polynomial given its evaluations at 0..len-1
+function poly_interp(len, v) {
+    assert(len <= 200);
+    var out[200];
+    for (var i = 0; i < len; i++) {
         out[i] = 0;
-        carries[i] = 0;
     }
-    for (var registerIdx = 0; registerIdx < k + ceilMN; registerIdx++) {
-        var thisRegisterValue = 0;
-        if (registerIdx > 0) {
-            thisRegisterValue = carries[registerIdx - 1];
-        }
 
-        var start = 0;
-        if (registerIdx >= ceilMN) {
-            start = registerIdx - ceilMN + 1;
+    // Product_{i=0..len-1} (x-i)
+    var full_poly[201];
+    full_poly[0] = 1;
+    for (var i = 0; i < len; i++) {
+        full_poly[i+1] = 0;
+        for (var j = i; j >= 0; j--) {
+            full_poly[j+1] += full_poly[j];
+            full_poly[j] *= -i;
         }
+    }
 
-        // go from start to min(registerIdx, len(pieces)-1)
-        for (var i = start; i <= registerIdx; i++) {
-            if (i < k) {
-                thisRegisterValue += pieces[i][registerIdx - i];
+    for (var i = 0; i < len; i++) {
+        var cur_v = 1;
+        for (var j = 0; j < len; j++) {
+            if (i == j) {
+                // do nothing
+            } else {
+                cur_v *= i-j;
             }
         }
+        cur_v = v[i] / cur_v;
 
-        if (isNegative(thisRegisterValue) == 1) {
-            var thisRegisterAbs = -1 * thisRegisterValue;
-            out[registerIdx] = (1<<n) - (thisRegisterAbs % (1<<n));
-            carries[registerIdx] = -1 * (thisRegisterAbs >> n) - 1;
-        } else {
-            out[registerIdx] = thisRegisterValue % (1<<n);
-            carries[registerIdx] = thisRegisterValue >> n;
+        var cur_rem = full_poly[len];
+        for (var j = len-1; j >= 0; j--) {
+            out[j] += cur_v * cur_rem;
+            cur_rem = full_poly[j] + i * cur_rem;
         }
+        assert(cur_rem == 0);
     }
 
     return out;
@@ -184,9 +178,15 @@ function long_scalar_mult(n, k, a, b) {
 // out[0] has length m + 1 -- quotient
 // out[1] has length k -- remainder
 // implements algorithm of https://people.eecs.berkeley.edu/~fateman/282/F%20Wright%20notes/week4.pdf
-// b[k-1] must be nonzero!
 function long_div(n, k, m, a, b){
     var out[2][100];
+    m += k;
+    while (b[k-1] == 0) {
+        out[1][k] = 0;
+        k--;
+        assert(k > 0);
+    }
+    m -= k;
 
     var remainder[200];
     for (var i = 0; i < m + k; i++) {

circuits/fp.circom

@@ -0,0 +1,84 @@
+pragma circom 2.0.3;
+
+include "../node_modules/circomlib/circuits/bitify.circom";
+
+include "./bigint.circom";
+include "./bigint_func.circom";
+
+// These functions operate over values in Z/Zp for some integer p (typically,
+// but not necessarily prime). Values are stored as standard bignums with k
+// chunks of n bits, but intermediate values often have "overflow" bits inside
+// various chunks.
+//
+// These Fp functions will always correctly generate witnesses mod p, but they
+// do not *check* that values are normalized to < p; they only check that
+// values are correct mod p. This is to save the comparison circuit.
+// They *will* always check for intended results mod p (soundness), but it may
+// not have a unique intermediate signal.
+//
+// Conversely, some templates may not be satisfiable if the input witnesses are
+// not < p. This does not break completeness, as honest provers will always
+// generate witnesses which are canonical (between 0 and p).
+
+// a * b = r mod p
+// a * b - p * q - r for some q
+template FpMul(n, k) {
+    assert(n + n + log_ceil(k) + 2 <= 252);
+    signal input a[k];
+    signal input b[k];
+    signal input p[k];
+
+    signal output out[k];
+
+    signal v_ab[2*k-1];
+    for (var x = 0; x < 2*k-1; x++) {
+        var v_a = poly_eval(k, a, x);
+        var v_b = poly_eval(k, b, x);
+        v_ab[x] <== v_a * v_b;
+    }
+
+    var ab[200] = poly_interp(2*k-1, v_ab);
+    // ab_proper has length 2*k
+    var ab_proper[200] = getProperRepresentation(n + n + log_ceil(k), n, 2*k-1, ab);
+
+    var long_div_out[2][100] = long_div(n, k, k, ab_proper, p);
+
+    // Since we're only computing a*b, we know that q < p will suffice, so we
+    // know it fits into k chunks and can do size n range checks.
+    signal q[k];
+    component q_range_check[k];
+    signal r[k];
+    component r_range_check[k];
+    for (var i = 0; i < k; i++) {
+        q[i] <-- long_div_out[0][i];
+        q_range_check[i] = Num2Bits(n);
+        q_range_check[i].in <== q[i];
+
+        r[i] <-- long_div_out[1][i];
+        r_range_check[i] = Num2Bits(n);
+        r_range_check[i].in <== r[i];
+    }
+
+    signal v_pq_r[2*k-1];
+    for (var x = 0; x < 2*k-1; x++) {
+        var v_p = poly_eval(k, p, x);
+        var v_q = poly_eval(k, q, x);
+        var v_r = poly_eval(k, r, x);
+        v_pq_r[x] <== v_p * v_q + v_r;
+    }
+
+    signal v_t[2*k-1];
+    for (var x = 0; x < 2*k-1; x++) {
+        v_t[x] <== v_ab[x] - v_pq_r[x];
+    }
+
+    var t[200] = poly_interp(2*k-1, v_t);
+    component tCheck = CheckCarryToZero(n, n + n + log_ceil(k) + 2, 2*k-1);
+    for (var i = 0; i < 2*k-1; i++) {
+        tCheck.in[i] <== t[i];
+    }
+
+    for (var i = 0; i < k; i++) {
+        out[i] <== r[i];
+    }
+}

package.json

@@ -17,7 +17,8 @@
   "dependencies": {
     "circomlib": "^2.0.2",
     "ethers": "^5.5.2",
-    "snarkjs": "^0.4.10"
+    "ffjavascript": "^0.2.51",
+    "snarkjs": "^0.4.15"
   },
   "devDependencies": {
     "@noble/secp256k1": "1.3.4",