diff --git a/groth16/dynamic/finish.nim b/groth16/dynamic/finish.nim
index 8247581..eb72eb8 100644
--- a/groth16/dynamic/finish.nim
+++ b/groth16/dynamic/finish.nim
@@ -28,6 +28,7 @@ import groth16/files/witness
 import groth16/math/domain
 import groth16/math/ntt
 import groth16/math/poly
+import groth16/math/convolution
 
 import groth16/prover/types
 import groth16/prover/shared
@@ -111,6 +112,9 @@ proc finishDynaProofWithMaskV1*( zkey: ZKey, wtns: Witness, dynaPreProof: DynaPr
 
   let partialMask = dynaPreProof.partialProof.partial_mask
 
+  let wvec    = dynaPreProof.dynaSetup.weightVec 
+  let wvecRev = fftReverseVec( wvec )
+
   var deltaAB: OnlyAB
   var witnessDelta: seq[Option[F]] = newSeq[Option[F]]( M )
   withMeasureTime(printTimings,"building deltaAz, deltaBz"):
@@ -125,18 +129,67 @@ proc finishDynaProofWithMaskV1*( zkey: ZKey, wtns: Witness, dynaPreProof: DynaPr
   withMeasureTime(printTimings,"finishing the linear terms"):
     proof = finishPartialProofWithMaskGeneric( false, zkey, wtns, dynaPreProof.partialProof, mask, pool, false )
 
-  var nonlin: G1
+  var nonlin: G1 = infG1
+
+#[ 
+
+  --- this only makes sense if we can restrict to a subgroup ---
+  
   var cs: seq[F]
   withMeasureTime(printTimings,"computing the nonlinear \"cross\" term"):
     withMeasureTime(printTimings," - computing the cross coeffs took"):
       cs = crossTermCoeffs(D , deltaAB.valuesAz , deltaAB.valuesBz )
     withMeasureTime(printTimings," - computing the cross MSM took"):
       nonlin = msmMultiThreadedG1( cs , dynaPreProof.dynaSetup.pointsDeltaLZ , pool )
-
+  
   echo " - nonzero coefficients in deltaA = " & $countNonZerosFr(deltaAB.valuesAz)
   echo " - nonzero coefficients in deltaB = " & $countNonZerosFr(deltaAB.valuesBz)
   echo " - nonzero coefficients in the cross-term = " & $countNonZerosFr(cs)
 
+]#
+
+  # following the Dynark paper
+  withMeasureTime(printTimings,"computing the nonlinear \"cross\" term (Dynark-style)"):
+
+    let idxs = trueIndices( dynaPreProof.deltaImages.imageAB ) 
+    let K = idxs.len
+
+    withMeasureTime(printTimings," - computing the diagonal part of the cross-term"):
+      var cs: seq[F]  = newSeq[F ]( K )
+      var ps: seq[G1] = newSeq[G1]( K )
+      for (k,i) in idxs.pairs:
+        cs[k] = deltaAB.valuesAz[i] * deltaAB.valuesBz[i] 
+        ps[k] = dynaPreProof.dynaSetup.diagPhiPoints[i]
+      nonlin += msmMultiThreadedG1( cs , ps , pool ) 
+
+    withMeasureTime(printTimings," - computing the off-diagonal part of the cross-term"):
+      let As = deltaAB.valuesAz 
+      let Bs = deltaAB.valuesBz 
+      let AstarW = fieldConvolution( As , wvecRev )
+      let BstarW = fieldConvolution( Bs , wvecRev )     # TODO: make this sparse somehow??!
+
+      var ds: seq[F] = newSeq[F]( K )
+      for (k,i) in idxs.pairs:
+        ds[k] = As[i] * BstarW[i] + Bs[i] * AstarW[i]    
+
+      let ls = selectTrues( dynaPreProof.deltaImages.imageAB , dynaPreProof.dynaSetup.pointsDeltaLZ )
+      nonlin += msmMultiThreadedG1( ds , ls , pool ) 
+
+#[
+      var ds: seq[F]  = newSeq[F ]( K )
+      var qs: seq[G1] = newSeq[G1]( K )
+      for (k,i) in idxs.pairs:
+        var x: F = zeroFr
+        for j in idxs: 
+          if (i != j):
+            let ab = deltaAB.valuesAz[i] * deltaAB.valuesBz[j]  +
+                     deltaAB.valuesAz[j] * deltaAB.valuesBz[i]
+            x += ab * wvec[ safeMod( j - i , N ) ]
+        ds[k] = x
+        qs[k] = dynaPreProof.dynaSetup.pointsDeltaLZ[i]
+      nonlin += msmMultiThreadedG1( ds , qs , pool ) 
+ ]#
+
   withMeasureTime(printTimings,"computing the nonlinear \"projection\" terms"):
     let comprAz = selectTrues( dynaPreProof.deltaImages.imageA , deltaAB.valuesAz )
     let comprBz = selectTrues( dynaPreProof.deltaImages.imageB , deltaAB.valuesBz )
diff --git a/groth16/dynamic/setup.nim b/groth16/dynamic/setup.nim
index f0b2739..1f3aa6f 100644
--- a/groth16/dynamic/setup.nim
+++ b/groth16/dynamic/setup.nim
@@ -3,7 +3,10 @@
 
 # import constantine/named/properties_fields
 
+import taskpools
+
 import groth16/bn128
+import groth16/bn128/arrays
 import groth16/misc
 
 import groth16/math/domain
@@ -26,15 +29,54 @@ import groth16/dynamic/shared
 
 #-------------------------------------------------------------------------------
 
+# the points `delta^{-1} * phi_ii(tau) * (tau^N - 1) * g1` as in the Dynark paper
+#
+# where
+#
+# > phi_ii = - sum_j phi_ij = - sum_j ( W[j-i]*L[i] + W[i-j]*L[j] )
+# > phi_ij = W[j-i]*L[i] + W[i-j]*L[j] 
+#
+func calculateDiagPhiFFT*( wvec: seq[F], deltaLZ: seq[G1] ): seq[G1] =
+  let N = wvec.len
+  assert( N == deltaLZ.len )
+
+  let sumW = sumSeqFr( wvec )
+
+  var hs: seq[G1] = groupConvolution( wvec , deltaLZ )
+  for i in 0..<N:
+    hs[i] += (sumW ** deltaLZ[i])
+    hs[i] =  negG1(hs[i])
+
+  return hs
+
+# NOTE: this is EXTREMELY SLOW
+proc calculateDiagPhiNaive*( wvec: seq[F], deltaLZ: seq[G1], pool: Taskpool ): seq[G1] =
+  let N = wvec.len
+  assert( N == deltaLZ.len )
+  
+  let sumW = sumSeqFr( wvec )
+
+  var diagPhi: seq[G1] = newSeq[G1]( N )  
+  for i in 0..<N:
+    var ws: seq[F] = newSeq[F]( N )
+    for j in 0..<N:
+      if i != j:
+        ws[j] = wvec[ safeMod(i-j , N) ] 
+      else:
+        ws[j] = sumW
+    diagPhi[i] = negG1( msmMultiThreadedG1( ws, deltaLZ, pool ) )
+
+  return diagPhi
+
+#-------------------------------------------------------------------------------
+
 # does the setup from the ZKey (prover key)
-proc dynaSetupV1FromZKey*(zkey: Zkey): DynaSetupV1 = 
+proc dynaSetupV1FromZKey*(zkey: Zkey, pool: Taskpool ): DynaSetupV1 = 
 
   let N = zkey.header.domainSize
   let D = createDomain(N)
 
-  # assert( zkey.header.flavour == JensGroth , "DynaSetupV1 requires classic (quotient) flavour, not Jordi's one!" )
-
-  var deltaZTau : seq[G1]                     # the points `delta^-1 * (tau^N-1) * tau^i * g1`
+  var deltaZTau : seq[G1]      # the points `delta^-1 * (tau^N-1) * tau^i * g1`
   case zkey.header.flavour     
 
     of JensGroth:
@@ -45,12 +87,22 @@ proc dynaSetupV1FromZKey*(zkey: Zkey): DynaSetupV1 =
       withMeasureTime(true,"Jordi-to-Jens conversion"):
         deltaZTau = convertPointsFromJordi(D , zkey.pPoints.pointsH1)
   
+  let sumW    = sumOfWVec( N )
   let wvec    = calculateWVec( D )
   let deltaLZ = inverseGroupFFT( deltaZTau , D )
   let conv    = groupConvolution( wvec , deltaLZ )
 
-  return DynaSetupV1( pointsDeltaLZ : deltaLZ , 
-                      wConvDeltaLZ  : conv    )
+  var diagPhi  : seq[G1]
+  var diagPhi2 : seq[G1]
+  withMeasureTime(true,"phi diagonal took"):
+    diagPhi = calculateDiagPhiFFT( wvec , deltaLZ )
+
+  echo "agree = " & $isEqualG1Seq( diagPhi , diagPhi2 ) 
+
+  return DynaSetupV1( weightVec     : wvec    ,
+                      pointsDeltaLZ : deltaLZ , 
+                      wConvDeltaLZ  : conv    ,
+                      diagPhiPoints : diagPhi )
 
 #---------------------------------------
 
@@ -68,10 +120,13 @@ func simulateDynaSetupV1*( D: Domain, tau: F, delta: F ): DynaSetupV1 =
     let y: F = deltaZTau * evalLagrangePolyAt( D, i, tau )    
     deltaLZ[i] = y ** gen1
 
-  let wvec   = calculateWVec( D )
-  let conv   = groupConvolution( wvec , deltaLZ )
+  let wvec    = calculateWVec( D )
+  let conv    = groupConvolution(    wvec , deltaLZ )
+  let diagPhi = calculateDiagPhiFFT( wvec , deltaLZ )
 
-  return DynaSetupV1( pointsDeltaLZ : deltaLZ , 
-                      wConvDeltaLZ  : conv    )
+  return DynaSetupV1( weightVec     : wvec    ,
+                      pointsDeltaLZ : deltaLZ , 
+                      wConvDeltaLZ  : conv    ,
+                      diagPhiPoints : diagPhi )
 
 #-------------------------------------------------------------------------------
diff --git a/groth16/dynamic/types.nim b/groth16/dynamic/types.nim
index d3952fe..964d0f7 100644
--- a/groth16/dynamic/types.nim
+++ b/groth16/dynamic/types.nim
@@ -36,17 +36,20 @@ type
 
   # things we can compute at circuit setup time
   DynaSetupV1* = object
+    weightVec*     : seq[F]       # the values `W_k` (it's surprisingly slow to re-compute in the finish part!)
     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`
+    diagPhiPoints* : seq[G1]      # the points `delta^-1 * phi_ii(tau) * Z(tau) * g1` from the Dynark paper
 
   # things we can compute from the partial witness
   DynaPreprocessV1* = object
     projA0*       : seq[G1]      # the points `delta^-1 * A0(tau) * L_i(tau) * g1` (well, "modulo Z(tau)")
     projB0*       : seq[G1]      # the same for B0
 
+  # "pre-proof" of the Dynark paper, which we can finish 
   DynaPreProofV1* = object
     partialProof*   : PartialProof       # linear part of the partial proof
-    deltaImages*    : DeltaImages       # image of the witness delta under A and B (where can update happen)
+    deltaImages*    : DeltaImages        # image of the witness delta under A and B (where can update happen)
     dynaSetup*      : DynaSetupV1        # circuit-setup time calculations
     dynaPreprocess* : DynaPreprocessV1   # "Dynark"-style preprocessing