more WIP dynark stuff; the cross-term coefficients seem to work

README.md

@@ -2,7 +2,7 @@
 Groth16 prover written in Nim
 -----------------------------
 
-This is Groth16 prover implementation in Nim, using the 
+This is a Groth16 prover implementation in Nim, using the 
 [`constantine`](https://github.com/mratsim/constantine)
 library as an arithmetic / curve backend.
 
@@ -10,6 +10,36 @@ The implementation is compatible with the `circom` + `snarkjs` ecosystem.
 
 At the moment only the `BN254` (aka. `alt-bn128`) curve is supported.
 
+### Partial proofs
+
+Groth16 has an interesting internal structure in that most of the proof
+computation is _linear_ in the witness.
+
+This implies that if you have to generate a lot of proofs where a large part
+of the witness is unchanged, then you can precalculate quite a lot, making
+subsequent proofs faster. A situation where this happens is when you have
+a long-term identity, which however you don't want to reveal (eg. anonymous 
+credentials).
+
+As splitting the proof into two has negligible overhead in practice, this
+is an essentially free win in such situations. However the speedup potential
+is limited, because the single nonlinear term depends on the full witness.
+In practice we've seen 2x-3x speedups when 10% of the witness changes.
+
+This is implemented under `groth16/partial/`.
+
+### Semi-dynamic proofs
+
+In a recent paper ["Dynark: Making Groth16 Dynamic"](https://eprint.iacr.org/2025/1897) 
+(2025), a solution was proposed for also speeding up the computation of the nonlinear 
+term. This results in an asymptotic speedup, however, unlike the previous 
+(essentially trivial) solution, this one is not _for free_: There is a somewhat 
+heavy precalculation when the full witness changes.
+
+Whether that's worth it depends on the particular use-case.
+
+A version of this idea is implemented under `groth16/dynamic/`.
+
 ### License
 
 Licensed and distributed under either of the
@@ -19,16 +49,15 @@ at your choice.
 
 ### TODO
 
-- [x] find and fix the _second_ totally surreal bug
 - [ ] clean up the code
-- [x] make it compatible with the latest constantine and also Nim 2.0.x
-- [x] make it a nimble package
+- [ ] more comprehensive testing
 - [ ] compare `.r1cs` to the "coeffs" section of `.zkey`
 - [x] generate fake circuit-specific setup ourselves
 - [x] make a CLI interface
 - [x] multithreading support (MSM, and possibly also FFT)
-- [ ] add Groth16 notes
+- [ ] add Groth16 explanatory notes
 - [ ] document the `snarkjs` circuit-specific setup `H` points convention
 - [x] precalculate stuff for "partial" proofs
 - [ ] implement Dynark style "dynamic proofs" too
+- [ ] benchmarks
 - [ ] make it work for different curves

groth16/bn128/arrays.nim

@@ -6,6 +6,29 @@ import groth16/bn128/fields
 import groth16/bn128/curves
 import groth16/bn128/rnd
 
+#-------------------------------------------------------------------------------
+
+# proper mathematical modulo operation.
+# (all mainstream programming languages are fucked up... Except Haskell, Haskell gets it right)
+#
+# TODO: maybe move this somewhere else?
+func safeMod*(x: int, N: int): int =
+  if x >= 0:
+    return (x mod N)
+  else:
+    let d = ((-x) div N)
+    let y = x + (d+1)*N
+    return (y mod N)
+
+#-------------------------------------------------------------------------------
+# generic arrays
+
+func constSeq*[T]( N: int , y: T ): seq[T] =
+  var xs : seq[T] = newSeq[T]( N )
+  for i in 0..<N:
+    xs[i] = y
+  return xs
+
 #-------------------------------------------------------------------------------
 # Fr arrays
 
@@ -48,6 +71,13 @@ func dotProdFr*(xs, ys: seq[Fr[BN254_Snarks]]): Fr[BN254_Snarks] =
     s += xs[i] * ys[i]
   return s
 
+func sumSeqFr*(xs : seq[Fr[BN254_Snarks]]): Fr[BN254_Snarks] =
+  let n = xs.len
+  var s : Fr[BN254_Snarks] = zeroFr
+  for i in 0..<n:
+    s += xs[i]
+  return s
+
 #-------------------------------------------------------------------------------
 # G1 arrays
 

groth16/bn128/fields.nim

@@ -135,6 +135,10 @@ func `/`*(x, y: Fp ): Fp  =  ( var z : Fp  = x ; z *= invFr(y) ; return z )
 func `/`*(x, y: Fp2): Fp2 =  ( var z : Fp2 = x ; z *= invFr(y) ; return z )
 func `/`*(x, y: Fr ): Fr  =  ( var z : Fr  = x ; z *= invFr(y) ; return z )
 
+func divBy2Fp* (x: Fp ): Fp = ( var z : Fp  = x ; div2(z) ; return z )
+func divBy2Fp2*(x: Fp2): Fp2= ( var z : Fp2 = x ; div2(z) ; return z )
+func divBy2Fr* (x: Fr ): Fr = ( var z : Fr  = x ; div2(z) ; return z )
+
 # template/generic instantiation of `pow_vartime` from here
 # /Users/bkomuves/.nimble/pkgs/constantine-0.0.1/constantine/math/arithmetic/finite_fields.nim(389, 7) template/generic instantiation of `fieldMod` from here
 # /Users/bkomuves/.nimble/pkgs/constantine-0.0.1/constantine/math/config/curves_prop_field_derived.nim(67, 5) Error: undeclared identifier: 'getCurveOrder'

groth16/dynamic/finish.nim

@@ -0,0 +1,104 @@
+
+#
+# Finishing a Dynark-style proof
+#
+# We are implementing the cross-term differently from the paper:
+# 
+# - in the paper there are only two convolutions (4 field FFT-s), but 2 MSM-s
+# - we have three convolutions (6 field FFT-s), but only 1 MSM. 
+#
+# In practice our version is significantly faster, as MSM-s are much more expensive
+# than field FFTs (at least for practical sizes)
+#
+
+import std/options
+
+import constantine/math/arithmetic
+import constantine/named/properties_fields
+
+import groth16/bn128
+import groth16/bn128/arrays
+
+import groth16/zkey_types
+
+import groth16/math/domain
+import groth16/math/ntt
+import groth16/math/poly
+
+import groth16/dynamic/types
+import groth16/dynamic/setup
+import groth16/dynamic/shared
+
+#-------------------------------------------------------------------------------
+
+#
+# Computes the expansion `f(x)` where `f(x) = A(x)B(x) mod (x^N - 1)` in
+# terms of `(x^N-1)*L_i(x)`
+#
+# The inputs `As` and `Bs` are the Lagrange-basis coefficients of A(x) and B(x).
+# 
+# Note: the remainder modulo `Z(x) := x^N-1` is `sum_k a[k]*b[k]*L_k(x)`
+# 
+func crossTermCoeffs*(D: Domain, As: seq[F], Bs: seq[F]) : seq[F] =
+
+  let N = D.domainSize
+  assert( N == As.len )
+  assert( N == Bs.len )
+
+  let ABs = pointwiseProdFr( As, Bs )
+
+  var Ahat  = forwardNTT( As  , D )
+  var Bhat  = forwardNTT( Bs  , D )
+  var ABhat = forwardNTT( ABs , D )
+
+  inplaceMulByFFTofWVecBar( Ahat  )
+  inplaceMulByFFTofWVecBar( Bhat  )
+  inplaceMulByFFTofWVecBar( ABhat )
+
+  let Aconv  = inverseNTT( Ahat  , D)
+  let Bconv  = inverseNTT( Bhat  , D)
+  let ABconv = inverseNTT( ABhat , D)
+
+  let sumW = sumOfWVec( N )
+
+  var output: seq[F] = newSeq[F]( N )
+ 
+  for k in 0..<N:
+    output[k] = As[k]*Bconv[k] + Bs[k]*Aconv[k] - ABconv[k] - ABs[k]*sumW
+
+  return output
+
+#-------------------------------------------------------------------------------
+
+proc testCrossTermCoeffs*(N : int): bool = 
+
+  let D = createDomain( N )
+
+  let As = randFrSeq(N)
+  let Bs = randFrSeq(N)
+
+  let tau  = randFr()
+  let ztau = smallPowFr(tau,N) - oneFr                      # Z(tau) = tau^N - 1
+
+  # reference `A(tau)*B(tau)`
+  var Atau: F = zeroFr
+  var Btau: F = zeroFr
+  for i in 0..<N:
+    Atau += As[i] * evalLagrangePolyAt( D , i , tau )       # A(x) = sum_i A_i * L_i(x)
+    Btau += Bs[i] * evalLagrangePolyAt( D , i , tau )
+  var reference = Atau * Btau
+
+  # correction term (because of the modulo Z(x) behaviour of what we test)
+  for k in 0..<N:
+    reference -= As[k] * Bs[k] * evalLagrangePolyAt( D , k , tau )
+
+  # the thing we want to test
+  let coeffs = crossTermCoeffs( D , As , Bs )
+  var smart: F = zeroFr
+  for i in 0..<N:
+    let lztau = zTau * evalLagrangePolyAt( D , i , tau ) 
+    smart += coeffs[i] * lztau
+
+  return (smart === reference)
+
+#-------------------------------------------------------------------------------

groth16/dynamic/preprocess.nim

@@ -26,6 +26,7 @@ import groth16/math/poly
 
 import groth16/dynamic/types
 import groth16/dynamic/setup
+import groth16/dynamic/shared
 
 #-------------------------------------------------------------------------------
 

groth16/dynamic/setup.nim

@@ -13,33 +13,9 @@ import groth16/math/convolution
 # import groth16/math/ntt
 
 import groth16/zkey_types
+
 import groth16/dynamic/types
-
-#-------------------------------------------------------------------------------
-
-# the "weight vector" from the Dynark paper
-#
-# these weights appear in the expansion of the
-# product of Lagrange polynomials `L_i(x)L_k(x)`)
-#
-func calculateWVec*( D: Domain ): seq[F] = 
-  let N = D.domainSize
-  var wvec : seq[F] = newSeq[F]( N )
-  let invN     : F = invFr( intToFr(N)  )
-  let invOmega : F = invFr( D.domainGen )
-  wvec[0] = zeroFr
-  for i in 1..<N:
-    wvec[i] = invN / ( smallPowFr(invOmega,i) - oneFr )
-  return wvec
-
-# reverse indexing vecBar[i] = vec[-i] 
-func fftReverseVec*[T]( vec: seq[T] ): seq[T] =
-  let N = vec.len
-  var vecBar: seq[T] = newSeq[T]( N )
-  vecBar[0] = vec[0]
-  for i in 1..<N:
-    vecBar[N-i] = vec[i]
-  return vecBar
+import groth16/dynamic/shared
 
 #-------------------------------------------------------------------------------
 

groth16/dynamic/shared.nim

@@ -0,0 +1,176 @@
+
+import std/options
+
+import constantine/math/arithmetic
+import constantine/named/properties_fields
+
+import groth16/bn128
+import groth16/bn128/arrays
+
+import groth16/math/domain
+import groth16/math/poly
+#import groth16/math/convolution
+
+import groth16/zkey_types
+
+import groth16/dynamic/types
+
+#-------------------------------------------------------------------------------
+
+# reverse indexing vecBar[i] = vec[-i] 
+func fftReverseVec*[T]( vec: seq[T] ): seq[T] =
+  let N = vec.len
+  var vecBar: seq[T] = newSeq[T]( N )
+  vecBar[0] = vec[0]
+  for i in 1..<N:
+    vecBar[N-i] = vec[i]
+  return vecBar
+
+#-------------------------------------------------------------------------------
+
+# the k-th element `W_k` of the "weight vector"
+func Wcoeff*( D: Domain, k0: int ): F = 
+  let N = D.domainSize
+  let k = safeMod( k0 , N )
+  let invOmega : F = invFr( D.domainGen )
+  return invFr( intToFr(N) * ( smallPowFr(invOmega,k) - oneFr ) )
+
+# the "weight vector" from the Dynark paper
+#
+# these weights appear in the expansion of the
+# product of Lagrange polynomials `L_i(x)L_k(x)`)
+#
+func calculateWVec*( D: Domain ): seq[F] = 
+  let N = D.domainSize
+  var wvec : seq[F] = newSeq[F]( N )
+  let invN     : F = invFr( intToFr(N)  )
+  let invOmega : F = invFr( D.domainGen )
+  wvec[0] = zeroFr
+  for i in 1..<N:
+    wvec[i] = invN / ( smallPowFr(invOmega,i) - oneFr )
+  return wvec
+
+func calculateWVecBar*( D: Domain ): seq[F] = 
+  fftReverseVec( calculateWVec(D) )
+
+func sumOfWVec*( N: int ): F = 
+  let fN : F = intToFr( N ) 
+  return (oneFr - fN) / (fN + fN)
+
+#-------------------------------------------------------------------------------
+
+# pointwise multiply by `FFT[W]_k = (k + (1-N)/2) / N`
+proc inplaceMulByFFTofWVec*( xs: var seq[F] ) = 
+  let N    = xs.len
+  let fN   = intToFr( N )
+  let invN = invFr(  fN )
+  let c    = divBy2Fr(oneFr - fN) * invN
+  for k in 0..<N:
+    let u = c + intToFr(k) * invN
+    xs[k] *= u
+
+# pointwise multiply by `FFT[Wbar]_k = Bar[FFT[W]]_k`
+proc inplaceMulByFFTofWVecBar*( xs: var seq[F] ) = 
+  let N    = xs.len
+  let fN   = intToFr( N )
+  let invN = invFr(  fN )
+  let c    = divBy2Fr(oneFr - fN) * invN
+  for k in 0..<N:
+    let i = (if k==0: 0 else: N-k)
+    let u = c + intToFr(i) * invN
+    xs[k] *= u
+
+#-------------------------------------------------------------------------------
+
+# computes the vectors A*z, B*z (but skips C*z)
+func buildOnlyAB*( zkey: ZKey, pwitness: seq[Option[F]] ): OnlyAB =
+  let hdr: GrothHeader = zkey.header
+  let domSize = hdr.domainSize
+
+  var valuesAz = newSeq[F](domSize)
+  var valuesBz = newSeq[F](domSize)
+ 
+  for entry in zkey.coeffs:
+    case entry.matrix 
+
+      of MatrixA: 
+        if isSome(pwitness[entry.col]):
+          valuesAz[entry.row] += entry.coeff * pwitness[entry.col].unsafeGet()
+
+      of MatrixB: 
+        if isSome(pwitness[entry.col]):
+          valuesBz[entry.row] += entry.coeff * pwitness[entry.col].unsafeGet()
+
+      else: raise newException(AssertionDefect, "fatal error")
+
+  return OnlyAB( valuesAz:valuesAz ,
+                 valuesBz:valuesBz )
+
+#---------------------------------------
+
+# computes the vectors A*z, B*z (but skips C*z), and also some image masks under A and B
+func buildPartialAB*( zkey: ZKey, pwitness: seq[Option[F]] ): PartialAB =
+  let hdr: GrothHeader = zkey.header
+  let domSize = hdr.domainSize
+
+  var valuesAz = newSeq[F](domSize)
+  var valuesBz = newSeq[F](domSize)
+ 
+  # we also compute the image of the complement of the partial witness under A and B
+  var complImageA = newSeq[bool](domSize) 
+  var complImageB = newSeq[bool](domSize) 
+  for i in 0..<domSize: 
+    complImageA[i] = false
+    complImageB[i] = false
+
+  for entry in zkey.coeffs:
+    case entry.matrix 
+
+      of MatrixA: 
+        if isSome(pwitness[entry.col]):
+          valuesAz[entry.row] += entry.coeff * pwitness[entry.col].unsafeGet()
+        else:
+          complImageA[entry.row] = true
+
+      of MatrixB: 
+        if isSome(pwitness[entry.col]):
+          valuesBz[entry.row] += entry.coeff * pwitness[entry.col].unsafeGet()
+        else:
+          complImageB[entry.row] = true
+
+      else: raise newException(AssertionDefect, "fatal error")
+
+  return PartialAB( valuesAz:valuesAz,
+                    valuesBz:valuesBz, 
+                    complImageA:complImageA,
+                    complImageB:complImageB )
+
+#-------------------------------------------------------------------------------
+# the phi(x) polynomials and Lagrange product decomposition (for testing purposes)
+
+# evaluates the `phi_{ik}(x)` polynomial at a point `tau`
+# we assume that `i != k`
+func evalPhiAt*(D: Domain, i: int, k: int, tau: F): F =
+  assert( not (i == k) )
+  return ( Wcoeff(D , k - i) * evalLagrangePolyAt(D , i , tau) +
+           Wcoeff(D , i - k) * evalLagrangePolyAt(D , k , tau) ) 
+           
+# evaluates the `phi_{ii}(x)` polynomial (as defined in the Dynark parper) at a point `tau`
+func evalPhiDiagonalAt*(D: Domain, i:int, tau: F): F =
+  let N = D.domainSize
+  var s: F = zeroFr
+  for k in 0..<N:
+    if not (k == i):
+      s -= evalPhiAt(D, i, k, tau)
+  return s
+
+# the product L_i(tau)L_k(tau) via the key lemma
+func evalLagrangeProductViaLemma*(D: Domain, i: int, k: int, tau: F): F =
+  let N = D.domainSize
+  let ztau = smallPowFr(tau, N) - oneFr
+  if i == k:
+    return ztau * evalPhiDiagonalAt(D, i, tau) + evalLagrangePolyAt(D, i, tau) 
+  else:
+    return ztau * evalPhiAt(D, i, k, tau)
+
+#-------------------------------------------------------------------------------

groth16/dynamic/types.nim

@@ -44,6 +44,10 @@ type
     partialProof*   : PartialProof
     dyanPreprocess* : DynaPreprocess
 
+  OnlyAB* = object
+    valuesAz*    : seq[F]
+    valuesBz*    : seq[F]
+
   PartialAB* = object
     valuesAz*    : seq[F]
     valuesBz*    : seq[F]

groth16/math/group_fft.nim

@@ -2,6 +2,8 @@
 #
 # FFT for groups elements
 #
+# TODO: flip the order of arguments, so that the domain is first...
+#
 
 #-------------------------------------------------------------------------------
 

groth16/partial/finish.nim

@@ -12,7 +12,7 @@ import constantine/math/arithmetic
 import constantine/named/properties_fields
 
 import groth16/bn128
-#import groth16/math/domain
+import groth16/math/domain
 import groth16/math/poly
 import groth16/zkey_types
 import groth16/files/witness

tests/groth16/testDyna.nim

@@ -8,22 +8,74 @@ import std/unittest
 
 import groth16/bn128
 import groth16/bn128/rnd
+import groth16/bn128/arrays
+
+import groth16/math/domain
+import groth16/math/ntt
+import groth16/math/poly
 
 import groth16/dynamic/types
+import groth16/dynamic/shared
+import groth16/dynamic/setup
 import groth16/dynamic/preprocess
-#import groth16/dynamic/setup
+import groth16/dynamic/finish
 
 #-------------------------------------------------------------------------------
 
-suite "dynamic proof tests":
+func ones(N : int): seq[F] = constSeq(N, oneFr)
+
+#-------------------------------------------------------------------------------
+
+suite "dynamic proof (Dynark) tests":
 
   let N : int = 128
-
+  let D = createDomain( N )
   let tau   : F = randFr()
   let delta : F = randFr()
 
-  test "projection term V1 calculation":
+  test "Wvec vs Wcoeff":
+    let xs = calculateWVec(D) 
+    var ys: seq[F] = newSeq[F]( N )
+    for i in 0..<N:
+      ys[i] = Wcoeff( D , i )
+    check isEqualFrSeq( xs , ys )
 
+  test "WvecBar vs Wcoeff":
+    let xs = calculateWVecBar(D) 
+    var ys: seq[F] = newSeq[F]( N )
+    for i in 0..<N:
+      ys[i] = Wcoeff( D , -i )
+    check isEqualFrSeq( xs , ys )
+
+  test "sum of Wvec":
+    check sumOfWVec(N) === sumSeqFr( calculateWVec(D) )
+
+  test "FFT of Wvec":
+    var xs = ones(N)
+    inplaceMulByFFTofWVec( xs ) 
+    let ys = forwardNTT( calculateWVec(D) , D ) 
+    check isEqualFrSeq( xs , ys )
+
+  test "FFT of WvecBar":
+    var xs = ones(N)
+    inplaceMulByFFTofWVecBar( xs ) 
+    let ys = forwardNTT( calculateWVecBar(D) , D ) 
+    check isEqualFrSeq( xs , ys)
+
+  test "expansion of the product of Lagrange polynomials":
+    var ok = true
+    for i in 0..<N:
+      for k in 0..<N:
+        let lhs = evalLagrangeProductViaLemma(D, i, k, tau)
+        let rhs = evalLagrangePolyAt(D, i, tau) * evalLagrangePolyAt(D, k, tau)
+        let local_ok = (lhs === rhs) 
+        ok = ok and local_ok
+    check ok
+
+  test "cross-term coefficients":
+    check testCrossTermCoeffs(N)
+
+  test "projection term V1 calculation":
     check testProjectionElementsV1( N , tau , delta )
 
 #-------------------------------------------------------------------------------

tests/groth16/testFFT.nim

@@ -26,17 +26,17 @@ suite "field NTT checks":
   let N = D.domainSize
   let xs : seq[Fr[BN254_Snarks]] = randFrSeq(N)
 
-  test "INTT(NTT(xs) == xs":
+  test "INTT(NTT(xs)) == xs":
     let ys = forwardNTT(xs, D)
     let zs = inverseNTT(ys, D) 
     check isEqualFrSeq( xs , zs )
   
-  test "NTT(INTT(xs) == xs":
+  test "NTT(INTT(xs)) == xs":
     let ys = inverseNTT(xs, D)
     let zs = forwardNTT(ys, D) 
     check isEqualFrSeq( xs , zs )
 
-  test "unscaledINTT(NTT(xs) == N * xs":
+  test "unscaledINTT(NTT(xs)) == N * xs":
     let ys = forwardNTT(xs, D)
     let zs = unscaledInverseNTT(ys, D) 
     var ws = xs
@@ -51,17 +51,17 @@ suite "group FFT checks (for the group G1)":
   let N = D.domainSize
   let gs : seq[G1] = randG1Seq(N)
 
-  test "IFFT(FFT(gs) == gs":
+  test "IFFT(FFT(gs)) == gs":
     let hs = forwardGroupFFT(gs, D)
     let rs = inverseGroupFFT(hs, D) 
     check isEqualG1Seq( gs , rs )
   
-  test "FFT(IFFT(xs) == xs":
+  test "FFT(IFFT(xs)) == xs":
     let hs = inverseGroupFFT(gs, D)
     let rs = forwardGroupFFT(hs, D) 
     check isEqualG1Seq( gs , rs )
 
-  test "FFT(unscaledIFFT(gs) == N * gs":
+  test "FFT(unscaledIFFT(gs)) == N * gs":
     let hs = unscaledInverseGroupFFT(gs, D)
     let rs = forwardGroupFFT(hs, D) 
     var qs = gs