diff --git a/research/kzg_poly_commit/fft_g1.nim b/research/kzg_poly_commit/fft_g1.nim
index 9f7c499..f256e72 100644
--- a/research/kzg_poly_commit/fft_g1.nim
+++ b/research/kzg_poly_commit/fft_g1.nim
@@ -10,7 +10,7 @@ import
   ../../constantine/config/curves,
   ../../constantine/[arithmetic, primitives],
   ../../constantine/elliptic/[
-    ec_endomorphism_accel,
+    ec_scalar_mul,
     ec_shortweierstrass_affine,
     ec_shortweierstrass_projective,
     ec_shortweierstrass_jacobian,
@@ -117,10 +117,10 @@ func simpleFT[EC; bits: static int](
 
   for i in 0 ..< L:
     last = vals[0]
-    last.scalarMulGLV_m2w2(rootsOfUnity[0])
+    last.scalarMul(rootsOfUnity[0])
     for j in 1 ..< L:
       v = vals[j]
-      v.scalarMulGLV_m2w2(rootsOfUnity[(i*j) mod L])
+      v.scalarMul(rootsOfUnity[(i*j) mod L])
       last += v
     output[i] = last
 
@@ -147,7 +147,7 @@ func fft_internal[EC; bits: static int](
   for i in 0 ..< half:
     # FFT Butterfly
     y_times_root = output[i+half]
-    y_times_root   .scalarMulGLV_m2w2(rootsOfUnity[i])
+    y_times_root   .scalarMul(rootsOfUnity[i])
     output[i+half] .diff(output[i], y_times_root)
     output[i]      += y_times_root
 
@@ -192,7 +192,7 @@ func ifft*[EC](
   let inv = invLen.toBig()
 
   for i in 0..< output.len:
-    output[i].scalarMulGLV_m2w2(inv)
+    output[i].scalarMul(inv)
 
   return FFTS_Success
 
@@ -276,7 +276,7 @@ when isMainModule:
 
     warmup()
 
-    for scale in 4 ..< 16:
+    for scale in 4 ..< 10:
       # Setup
 
       let desc = FFTDescriptor[G1].init(uint8 scale)
diff --git a/research/kzg_poly_commit/polynomials.nim b/research/kzg_poly_commit/polynomials.nim
new file mode 100644
index 0000000..8bab5fc
--- /dev/null
+++ b/research/kzg_poly_commit/polynomials.nim
@@ -0,0 +1,67 @@
+import
+  ../../constantine/config/curves,
+  ../../constantine/[arithmetic, primitives],
+  ../../constantine/elliptic/[
+    ec_scalar_mul,
+    ec_shortweierstrass_projective,
+  ],
+  ../../constantine/io/[io_fields, io_ec],
+  ../../constantine/pairings/[
+    pairings_bls12,
+    miller_loops
+  ]
+
+type
+  G1 = ECP_ShortW_Prj[Fp[BLS12_381], NotOnTwist]
+  G2 = ECP_ShortW_Prj[Fp2[BLS12_381], OnTwist]
+  G1aff = ECP_ShortW_Aff[Fp[BLS12_381], NotOnTwist]
+  G2aff = ECP_ShortW_Aff[Fp2[BLS12_381], OnTwist]
+  GT = Fp12[BLS12_381]
+
+func linear_combination*(
+       r: var ,
+       points: openarray[G1],
+       coefs: openarray[Fr[BLS12_381]]
+     ) =
+  ## Polynomial evaluation
+  ## TODO: multi scalar mul
+  doAssert points.len == coefs.len
+
+  r.setInf()
+  for i in 0 ..< points.len:
+    var tmp = points[i]
+    tmp.scalarMul(coefs[i].toBig())
+    r += tmp
+
+func pair_verify*(
+       P1: G1,
+       Q1: G2,
+       P2: G1,
+       Q2: G2,
+     ): bool =
+  ## TODO, multi-pairings.
+
+  ## Affine
+  var P1a, P2a: G1aff
+  var Q1a, Q2a: G2aff
+
+  P1a.affineFromProjective(P1)
+  Q1a.affineFromProjective(Q1)
+  P2a.affineFromProjective(P2)
+  Q2a.affineFromProjective(Q2)
+
+  # To verify if e(P1, Q1) == e(P2, Q2)
+  # we can do e(P1, Q1) / e(P2, Q2) == 1
+  # <=> e(P1, Q1) . e(P2, Q2)^-1
+  # <=> e(P1, Q1) . e(-P2, Q2) due to pairings bilinearity
+  # we can negate any of the points but it's cheaper to use a G1
+  P1a.neg()
+
+  # Merge 2 miller loops.
+  var gt1, gt2: GT
+  gt1.millerLoopAddchain(Q1a, P1a)
+  gt2.millerLoopAddchain(Q2a, P2a)
+
+  gt1 *= gt2
+  gt.finalExpEasy()
+  gt.finalExpHard_BLS12()