diff --git a/groth16/bn128/curves.nim b/groth16/bn128/curves.nim
index b4e0d1d..ac44cc9 100644
--- a/groth16/bn128/curves.nim
+++ b/groth16/bn128/curves.nim
@@ -321,6 +321,15 @@ func addG1*(p,q: G1): G1 =
   prj.affine(s, r)
   return s
 
+func subG1*(p,q: G1): G1 =
+  var r, x, y : ProjG1
+  prj.fromAffine(x, p)
+  prj.fromAffine(y, q)
+  prj.diff(r, x, y)
+  var s : G1
+  prj.affine(s, r)
+  return s
+
 #---------------------------------------
 
 func addG2*(p,q: G2): G2 =
@@ -332,6 +341,15 @@ func addG2*(p,q: G2): G2 =
   prj.affine(s, r)
   return s
 
+func subG2*(p,q: G2): G2 =
+  var r, x, y : ProjG2
+  prj.fromAffine(x, p)
+  prj.fromAffine(y, q)
+  prj.diff(r, x, y)
+  var s : G2
+  prj.affine(s, r)
+  return s
+
 func negG1*(p: G1): G1 =
   var r : G1 = p
   neg(r)
@@ -347,11 +365,14 @@ func negG2*(p: G2): G2 =
 func `+`*(p,q: G1): G1 = addG1(p,q)
 func `+`*(p,q: G2): G2 = addG2(p,q)
 
+func `-`*(p,q: G1): G1 = subG1(p,q)
+func `-`*(p,q: G2): G2 = subG2(p,q)
+
 func `+=`*(p: var G1, q: G1) =    p = addG1(p,q)
 func `+=`*(p: var G2, q: G2) =    p = addG2(p,q)
 
-func `-=`*(p: var G1, q: G1) =    p = addG1(p,negG1(q))
-func `-=`*(p: var G2, q: G2) =    p = addG2(p,negG2(q))
+func `-=`*(p: var G1, q: G1) =    p = subG1(p,q)
+func `-=`*(p: var G2, q: G2) =    p = subG2(p,q)
 
 #-------------------------------------------------------------------------------
 #
diff --git a/groth16/dynamic/README.md b/groth16/dynamic/README.md
new file mode 100644
index 0000000..1fd7181
--- /dev/null
+++ b/groth16/dynamic/README.md
@@ -0,0 +1,14 @@
+
+"Semi-dynamic" Groth16 proofs
+-----------------------------
+
+This is loosely based on the [Dynark paper](https://eprint.iacr.org/2025/1897):
+
+- _"Dynark: Making Groth16 Dynamic"_ by Tianyu Zhang, Yupeng Ouyang and Yupeng Zhang
+
+See also [this write-up](https://hackmd.io/@bkomuves/HyBL5V5xMe) (mostly following the paper).
+
+Here some details are different from the paper though, and our implementation 
+is noticeably more efficient in practice.
+
+
diff --git a/groth16/dynamic/preprocess.nim b/groth16/dynamic/preprocess.nim
new file mode 100644
index 0000000..1591f4e
--- /dev/null
+++ b/groth16/dynamic/preprocess.nim
@@ -0,0 +1,174 @@
+
+#
+# the goal of Dynark-style preprocessing is to compute the `2d` group elements
+#
+#  U_k := L_k(\tau) * A0(tau) * g1 
+#  V_k := L_k(\tau) * B0(tau) * g1 
+#
+# where A0(x), B0(x) are the polynomials corresponding to the partial witness
+#
+# this is V1, the version closest to the Dynark paper
+#
+
+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/convolution
+import groth16/math/poly
+
+import groth16/dynamic/types
+import groth16/dynamic/setup
+
+#-------------------------------------------------------------------------------
+
+# computes the vectors A*z, B*z (but skips C*z)
+func buildPartialAB*( zkey: ZKey, pwitness: seq[Option[Fr[BN254_Snarks]]] ): PartialAB =
+  let hdr: GrothHeader = zkey.header
+  let domSize = hdr.domainSize
+
+  var valuesAz = newSeq[Fr[BN254_Snarks]](domSize)
+  var valuesBz = newSeq[Fr[BN254_Snarks]](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 )
+
+#-------------------------------------------------------------------------------
+
+func projectionElementsV1*(setup: DynaSetupV1, D: Domain, As: seq[F], complImage: seq[bool]): seq[G1] =
+  let N = D.domainSize
+  var Us: seq[G1] = newSeq[G1]( N )
+
+  # reverse indexed wvec: wvecBar[i] = wvec[-i]
+  let wvecBar: seq[F] = fftReverseVec( setup.weightVec )
+
+  let fldN    : F = intToFr( N ) 
+  let negSumW : F = (fldN - oneFr) / (fldN + fldN) 
+
+  let WBarStarA : seq[F]  = fieldConvolution( As , wvecBar ) 
+  let ALstarW   : seq[G1] = groupConvolution( setup.weightVec , pointwiseScaleG1( As , setup.pointsDeltaLZ ) )
+
+  # 2*d scalar multiplications + the group convolution above
+  for k in 0..<N: 
+ 
+    # we only compute for the _image of_ the complementer of the partial witness
+    if complImage[k]:         
+      let cf : F = WBarStarA[k] - As[k] * negSumW
+      Us[k] = ALstarW[k] - (As[k] ** setup.wConvDeltaLZ[k]) + (cf ** setup.pointsDeltaLZ[k])
+
+  return Us
+
+#---------------------------------------
+
+# for testing purposes
+func simulateProjectionElementsV1*( D: Domain, tau: F, delta: F, As: seq[F], complImage: seq[bool] ): seq[G1] =
+
+  let N = D.domainSize
+ 
+  let setup = simulateDynaSetupV1( D, tau, delta )
+
+  let deltaInv : F = invFr(delta)
+
+  # reverse indexed wvec: wvecBar[i] = wvec[-i]
+  let wvecBar: seq[F] = fftReverseVec( setup.weightVec )
+
+  # sum_i A[i]*L_i(tau)
+  var asLTau: F = zeroFr
+  for i in 0..<N:
+    asLTau += As[i] * evalLagrangePolyAt( D, i, tau )
+
+  # delta^-1 * sum A[i] * L_i(tau)
+  var deltaAsLTau: F = deltaInv * asLTau 
+ 
+  #
+  # a note about the correction term
+  #
+  # so in the actual protocol, we simply calculate things modulo `(x^N - 1)`. 
+  # This very useful because otherwise we would have to change the trusted setup ceremony...
+  # however, to make this simulation compatible, we have to do the same here!
+  #
+  # fortunately, here we know explicitly that the remainder modulo `(x^N - 1)` 
+  # in `U[k]` is `delta^-1 * As[k] * L_k(x)`
+  #
+
+  var Us: seq[G1] = newSeq[G1]( N )
+  for k in 0..<N:
+    if complImage[k]:
+      let Lk_tau : F = evalLagrangePolyAt(D, k, tau)
+      let y      : F = deltaAsLTau * Lk_tau
+      let corr   : F = deltaInv * As[k] * Lk_tau
+      Us[k] = (y - corr) ** gen1
+ 
+  return Us
+
+#---------------------------------------
+
+proc testProjectionElementsV1*(N: int, tau: F, delta: F ): bool =
+
+  let D  = createDomain( N )
+  let As = randFrSeq( N )
+
+  # image of the changes mask, let's just compute everything, easier for testing
+  var trues: seq[bool] = newSeq[bool](N) ; for i in 0..<N: trues[i] = true
+
+  let setup = simulateDynaSetupV1( D, tau, delta )
+
+  let simulated = simulateProjectionElementsV1(D, tau, delta, As, trues)
+  let computed  = projectionElementsV1( setup, D, As, trues )
+
+  return isEqualG1Seq( simulated , computed )
+
+#---------------------------------------
+
+func dynaPreprocessV1*(zkey: ZKey, setup: DynaSetupV1, partialWitness: PartialWitness): DynaPreprocess =
+  let N = zkey.header.domainSize
+  let D = createDomain(N)
+
+  # let wtnsMask = partialWitnessMask(partialWitness)
+
+  let partialAB = buildPartialAB( zkey, partialWitness.values )
+
+  let projA = projectionElementsV1( setup , D , partialAB.valuesAz , partialAB.complImageA )
+  let projB = projectionElementsV1( setup , D , partialAB.valuesBz , partialAB.complImageB )
+  
+  return DynaPreprocess( projA0: projA, projB0: projB )
+
+#-------------------------------------------------------------------------------
+
+
+
+
+
diff --git a/groth16/dynamic/setup.nim b/groth16/dynamic/setup.nim
new file mode 100644
index 0000000..8ee1322
--- /dev/null
+++ b/groth16/dynamic/setup.nim
@@ -0,0 +1,85 @@
+
+{.push raises:[].}
+
+# import constantine/named/properties_fields
+
+import groth16/bn128
+
+import groth16/math/domain
+import groth16/math/group_fft
+import groth16/math/poly
+import groth16/math/convolution
+# import groth16/math/convert
+# 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
+
+#-------------------------------------------------------------------------------
+
+# does the setup from the ZKey (prover key)
+func dynaSetupV1FromZKey*(zkey: Zkey): DynaSetupV1 = 
+
+  assert( zkey.header.flavour == JensGroth , "DynaSetupV1 requires classic (quotient) flavour, not Jordi's one!" )
+
+  let N = zkey.header.domainSize
+  let D = createDomain(N)
+  
+  let deltaLZ = inverseGroupFFT( zkey.pPoints.pointsH1 , D )
+  let wvec    = calculateWVec( D )
+  let conv    = groupConvolution( wvec , deltaLZ )
+
+  return DynaSetupV1( pointsDeltaLZ : deltaLZ , 
+                      weightVec     : wvec    ,
+                      wConvDeltaLZ  : conv    )
+
+#---------------------------------------
+
+# simulates a setup (from the "toxic waste" values `tau` and `delta`), so that 
+# we can test components of the system
+func simulateDynaSetupV1*( D: Domain, tau: F, delta: F ): DynaSetupV1 = 
+
+  let N = D.domainSize
+  let ztau:      F = smallPowFr(tau,N) - oneFr
+  let deltaZTau: F = ztau / delta
+
+  # compute `delta^-1 * L_i(tau) * (tau^N - 1) ** g1 
+  var deltaLZ: seq[G1] = newSeq[G1]( N ) 
+  for i in 0..<N: 
+    let y: F = deltaZTau * evalLagrangePolyAt( D, i, tau )    
+    deltaLZ[i] = y ** gen1
+
+  let wvec   = calculateWVec( D )
+  let conv   = groupConvolution( wvec , deltaLZ )
+
+  return DynaSetupV1( weightVec     : wvec    ,
+                      pointsDeltaLZ : deltaLZ , 
+                      wConvDeltaLZ  : conv    )
+
+#-------------------------------------------------------------------------------
diff --git a/groth16/dynamic/types.nim b/groth16/dynamic/types.nim
index ff2759e..94c922f 100644
--- a/groth16/dynamic/types.nim
+++ b/groth16/dynamic/types.nim
@@ -31,9 +31,9 @@ type
 
   # things we can compute at circuit setup time
   DynaSetupV1* = object
-    pointsDeltaL* : seq[G1]      # the points `delta^-1 * L_i(tau) * Z(tau) * g1`
-    weightVec*    : seq[F]       # the weights `W_k = 1/N/(omega^-k - 1)`
-    wConvDeltaL*  : seq[G1]      # the convolution of `W` and `deltaL`
+    weightVec*     : seq[F]       # the weights `W_k = 1/N/(omega^-k - 1)`
+    pointsDeltaLZ* : seq[G1]      # the points `delta^-1 * L_i(tau) * Z(tau) * g1` where `Z(x) = x^N-1`
+    wConvDeltaLZ*  : seq[G1]      # the convolution of `W` and `pointsDeltaLZ`
 
   # things we can compute from the partial witness
   DynaPreprocess* = object
diff --git a/tests/groth16/testCurve.nim b/tests/groth16/testCurve.nim
index 9c1667d..d70c35a 100644
--- a/tests/groth16/testCurve.nim
+++ b/tests/groth16/testCurve.nim
@@ -9,8 +9,32 @@ import constantine/math/io/io_fields
 import constantine/named/properties_fields 
 import constantine/math/extension_fields/towers 
 
+import groth16/bn128
 import groth16/bn128/fields
 import groth16/bn128/curves
+import groth16/bn128/rnd
+
+#-------------------------------------------------------------------------------
+
+suite "curve arithmetic":
+
+  test "G1 add/sub":
+    let g = randG1()
+    let h = randG1()
+    let lhs = g - h
+    let rhs = g + negG1(h)
+    var x = g
+    x -= h
+    check ((lhs === rhs) and (x === lhs))
+
+  test "G2 add/sub":
+    let g = randG2()
+    let h = randG2()
+    let lhs = g - h
+    let rhs = g + negG2(h)
+    var x = g
+    x -= h
+    check ((lhs === rhs) and (x === lhs))
 
 #-------------------------------------------------------------------------------
 
diff --git a/tests/groth16/testDyna.nim b/tests/groth16/testDyna.nim
new file mode 100644
index 0000000..9f1c8d9
--- /dev/null
+++ b/tests/groth16/testDyna.nim
@@ -0,0 +1,29 @@
+
+{.used.} 
+
+import std/unittest
+
+# import constantine/math/arithmetic
+# import constantine/named/properties_fields
+
+import groth16/bn128
+import groth16/bn128/rnd
+
+import groth16/dynamic/types
+import groth16/dynamic/preprocess
+#import groth16/dynamic/setup
+
+#-------------------------------------------------------------------------------
+
+suite "dynamic proof tests":
+
+  let N : int = 128
+
+  let tau   : F = randFr()
+  let delta : F = randFr()
+
+  test "projection term V1 calculation":
+
+    check testProjectionElementsV1( N , tau , delta )
+
+#-------------------------------------------------------------------------------
diff --git a/tests/test.nim b/tests/test.nim
index 4ff0c94..f46e1a6 100644
--- a/tests/test.nim
+++ b/tests/test.nim
@@ -4,6 +4,7 @@ import ./groth16/testCurve
 import ./groth16/testPoly
 import ./groth16/testFFT
 import ./groth16/testLagrangeTau
+import ./groth16/testDyna
 import ./groth16/testPtCompression
 import ./groth16/testProver
 import ./groth16/testMultithreading