diff --git a/benchmarks/bench_fp6.nim b/benchmarks/bench_fp6.nim
index 53f4995..93a3c8a 100644
--- a/benchmarks/bench_fp6.nim
+++ b/benchmarks/bench_fp6.nim
@@ -40,10 +40,10 @@ proc main() =
   echo "-".repeat(80)
   staticFor i, 0, AvailableCurves.len:
     const curve = AvailableCurves[i]
-    # addBench(Fp6[curve], Iters)
-    # subBench(Fp6[curve], Iters)
-    # negBench(Fp6[curve], Iters)
-    # mulBench(Fp6[curve], Iters)
+    addBench(Fp6[curve], Iters)
+    subBench(Fp6[curve], Iters)
+    negBench(Fp6[curve], Iters)
+    mulBench(Fp6[curve], Iters)
     sqrBench(Fp6[curve], Iters)
     # invBench(Fp6[curve], InvIters)
     echo "-".repeat(80)
diff --git a/constantine.nimble b/constantine.nimble
index f4f35e5..9278846 100644
--- a/constantine.nimble
+++ b/constantine.nimble
@@ -48,8 +48,9 @@ task test, "Run all tests":
 
   test "", "tests/test_finite_fields_vs_gmp.nim"
 
-  # 𝔽p2
+  # Towers of extension fields
   test "", "tests/test_fp2.nim"
+  test "", "tests/test_fp6.nim"
 
   if sizeof(int) == 8: # 32-bit tests
     # Primitives
@@ -70,8 +71,9 @@ task test, "Run all tests":
 
     test "-d:Constantine32", "tests/test_finite_fields_vs_gmp.nim"
 
-    # 𝔽p2
+    # Towers of extension fields
     test "", "tests/test_fp2.nim"
+    test "", "tests/test_fp6.nim"
 
 task test_no_gmp, "Run tests that don't require GMP":
   # -d:testingCurves is configured in a *.nim.cfg for convenience
@@ -90,8 +92,9 @@ task test_no_gmp, "Run tests that don't require GMP":
   test "", "tests/test_finite_fields_mulsquare.nim"
   test "", "tests/test_finite_fields_powinv.nim"
 
-  # 𝔽p2
+  # Towers of extension fields
   test "", "tests/test_fp2.nim"
+  test "", "tests/test_fp6.nim"
 
   if sizeof(int) == 8: # 32-bit tests
     # Primitives
@@ -108,8 +111,9 @@ task test_no_gmp, "Run tests that don't require GMP":
     test "-d:Constantine32", "tests/test_finite_fields_mulsquare.nim"
     test "-d:Constantine32", "tests/test_finite_fields_powinv.nim"
 
-    # 𝔽p2
+    # Towers of extension fields
     test "", "tests/test_fp2.nim"
+    test "", "tests/test_fp6.nim"
 
 proc runBench(benchName: string, compiler = "") =
   if not dirExists "build":
diff --git a/constantine/tower_field_extensions/abelian_groups.nim b/constantine/tower_field_extensions/abelian_groups.nim
index c3d5bfe..65b46c0 100644
--- a/constantine/tower_field_extensions/abelian_groups.nim
+++ b/constantine/tower_field_extensions/abelian_groups.nim
@@ -99,6 +99,19 @@ func diff*(r: var CubicExtAddGroup, a, b: CubicExtAddGroup) =
   r.c1.diff(a.c1, b.c1)
   r.c2.diff(a.c2, b.c2)
 
+func neg*(r: var CubicExtAddGroup, a: CubicExtAddGroup) =
+  ## Negate ``a`` and store the result into r
+  r.c0.neg(a.c0)
+  r.c1.neg(a.c1)
+  r.c2.neg(a.c2)
+
+# TODO: The type system is lost with out-of-place result concepts
+# func `+`*(a, b: CubicExtAddGroup): CubicExtAddGroup =
+#   result.sum(a, b)
+#
+# func `-`*(a, b: CubicExtAddGroup): CubicExtAddGroup =
+#   result.sum(a, b)
+
 # ############################################################
 #
 #                 Quadratic Extension fields
@@ -165,3 +178,10 @@ func neg*(r: var QuadExtAddGroup, a: QuadExtAddGroup) =
   ## Negate ``a`` into ``r``
   r.c0.neg(a.c0)
   r.c1.neg(a.c1)
+
+# TODO: The type system is lost with out-of-place result concepts
+# func `+`*(a, b: CubicExtAddGroup): CubicExtAddGroup =
+#   result.sum(a, b)
+#
+# func `-`*(a, b: CubicExtAddGroup): CubicExtAddGroup =
+#   result.sum(a, b)
diff --git a/constantine/tower_field_extensions/fp2_complex.nim b/constantine/tower_field_extensions/fp2_complex.nim
index 25ea7fc..6921799 100644
--- a/constantine/tower_field_extensions/fp2_complex.nim
+++ b/constantine/tower_field_extensions/fp2_complex.nim
@@ -157,3 +157,9 @@ func inv*(r: var Fp2, a: Fp2) =
   r.c0.prod(a.c0, t0)  # r0 = a0 / (a0² + a1²)
   t1.neg(t0)           # t0 = -1 / (a0² + a1²)
   r.c1.prod(a.c1, t1)  # r1 = -a1 / (a0² + a1²)
+
+func `*=`*(a: var Fp2, b: Fp2) {.inline.} =
+  a.prod(a, b)
+
+func `*`*(a, b: Fp2): Fp2 {.inline.} =
+  result.prod(a, b)
diff --git a/constantine/tower_field_extensions/fp6_1_plus_i.nim b/constantine/tower_field_extensions/fp6_1_plus_i.nim
index a223c62..e31c607 100644
--- a/constantine/tower_field_extensions/fp6_1_plus_i.nim
+++ b/constantine/tower_field_extensions/fp6_1_plus_i.nim
@@ -50,17 +50,25 @@ func `*`(_: typedesc[Xi], a: Fp2): Fp2 =
   ## Multiply an element of 𝔽p2 by 𝔽p6 cubic non-residue 1 + 𝑖
   ## (c0 + c1 𝑖) (1 + 𝑖) => c0 + (c0 + c1)𝑖 + c1 𝑖²
   ##                     => c0 - c1 + (c0 + c1) 𝑖
-  result.c0 = a.c0 - a.c1
-  result.c1 = a.c0 + a.c1
+  result.c0.diff(a.c0, a.c1)
+  result.c1.sum(a.c0, a.c1)
 
 template `*`(a: Fp2, _: typedesc[Xi]): Fp2 =
   Xi * a
 
+func `*=`(a: var Fp2, _: typedesc[Xi]) =
+  ## Inplace multiply an element of 𝔽p2 by 𝔽p6 cubic non-residue 1 + 𝑖
+  let t = a.c0
+  a.c0 -= a.c1
+  a.c1 += t
+
 func square*[C](r: var Fp6[C], a: Fp6[C]) =
-  ## Return a²
-  ##
-  # Algorithm is Chung-Hasan Squaring SQR3
+  ## Returns r = a²
+  # Algorithm is Chung-Hasan Squaring SQR2
   # http://cacr.uwaterloo.ca/techreports/2006/cacr2006-24.pdf
+  #
+  # TODO: change to SQR1 or SQR3 (requires div2)
+  #       which are faster for the sizes we are interested in.
   var v2{.noInit.}, v3{.noInit.}, v4{.noInit.}, v5{.noInit.}: Fp2[C]
 
   v4.prod(a.c0, a.c1)
@@ -80,3 +88,37 @@ func square*[C](r: var Fp6[C], a: Fp6[C]) =
   r.c2.sum(v2, v4)
   r.c2 += v5
   r.c2 -= v3
+
+func prod*[C](r: var Fp6[C], a, b: Fp6[C]) =
+  ## Returns r = a * b
+  # Algorithm is Karatsuba
+  var v0{.noInit.}, v1{.noInit.}, v2{.noInit.}, t{.noInit.}: Fp2[C]
+
+  v0.prod(a.c0, b.c0)
+  v1.prod(a.c1, b.c1)
+  v2.prod(a.c2, b.c2)
+
+  # r.c0 = ((a.c1 + a.c2) * (b.c1 + b.c2) - v1 - v2) * Xi + v0
+  r.c0.sum(a.c1, a.c2)
+  t.sum(b.c1, b.c2)
+  r.c0 *= t
+  r.c0 -= v1
+  r.c0 -= v2
+  r.c0 *= Xi
+  r.c0 += v0
+
+  # r.c1 = (a.c0 + a.c1) * (b.c0 + b.c1) - v0 - v1 + Xi * v2
+  r.c1.sum(a.c0, a.c1)
+  t.sum(b.c0, b.c1)
+  r.c1 *= t
+  r.c1 -= v0
+  r.c1 -= v1
+  r.c1 += Xi * v2
+
+  # r.c2 = (a.c0 + a.c2) * (b.c0 + b.c2) - v0 - v2 + v1
+  r.c2.sum(a.c0, a.c2)
+  t.sum(b.c0, b.c2)
+  r.c2 *= t
+  r.c2 -= v0
+  r.c2 -= v2
+  r.c2 += v1
diff --git a/tests/test_fp6.nim b/tests/test_fp6.nim
index 5d4f925..63c27fe 100644
--- a/tests/test_fp6.nim
+++ b/tests/test_fp6.nim
@@ -43,16 +43,14 @@ suite "𝔽p6 = 𝔽p2[∛(1+𝑖)] (irreducible polynomial x³ - (1+𝑖))":
             var r{.noinit.}: Fp6[C]
             r.square(One)
             check: bool(r == One)
-          # block:
-          #   var r{.noinit.}: Fp6[C]
-          #   r.prod(One, One)
-          #   check: bool(r == One)
+          block:
+            var r{.noinit.}: Fp6[C]
+            r.prod(One, One)
+            check: bool(r == One)
 
         testInstance()
 
     test(BN254)
-    test(P256)
-    test(Secp256k1)
     test(BLS12_377)
     test(BLS12_381)
     test(BN446)
@@ -75,16 +73,20 @@ suite "𝔽p6 = 𝔽p2[∛(1+𝑖)] (irreducible polynomial x³ - (1+𝑖))":
           var Four: Fp6[C]
           Four.double(Two)
 
-          var r: Fp6[C]
-          r.square(Two)
+          block:
+            var r: Fp6[C]
+            r.square(Two)
 
-          check: bool(r == Four)
+            check: bool(r == Four)
+          block:
+            var r: Fp6[C]
+            r.prod(Two, Two)
+
+            check: bool(r == Four)
 
         testInstance()
 
     test(BN254)
-    test(P256)
-    test(Secp256k1)
     test(BLS12_377)
     test(BLS12_381)
     test(BN446)
@@ -109,16 +111,20 @@ suite "𝔽p6 = 𝔽p2[∛(1+𝑖)] (irreducible polynomial x³ - (1+𝑖))":
           for _ in 0 ..< 9:
             Nine += One
 
-          var u: Fp6[C]
-          u.square(Three)
+          block:
+            var u: Fp6[C]
+            u.square(Three)
 
-          check: bool(u == Nine)
+            check: bool(u == Nine)
+          block:
+            var u: Fp6[C]
+            u.prod(Three, Three)
+
+            check: bool(u == Nine)
 
         testInstance()
 
     test(BN254)
-    test(P256)
-    test(Secp256k1)
     test(BLS12_377)
     test(BLS12_381)
     test(BN446)
@@ -143,19 +149,80 @@ suite "𝔽p6 = 𝔽p2[∛(1+𝑖)] (irreducible polynomial x³ - (1+𝑖))":
           for _ in 0 ..< 9:
             Nine += One
 
-          var u: Fp6[C]
-          u.square(MinusThree)
+          block:
+            var u: Fp6[C]
+            u.square(MinusThree)
 
-          check: bool(u == Nine)
+            check: bool(u == Nine)
+          block:
+            var u: Fp6[C]
+            u.prod(MinusThree, MinusThree)
+
+            check: bool(u == Nine)
 
         testInstance()
 
     test(BN254)
-    test(P256)
-    test(Secp256k1)
     test(BLS12_377)
     test(BLS12_381)
     test(BN446)
     test(FKM12_447)
     test(BLS12_461)
     test(BN462)
+
+  test "Multiplication by 0 and 1":
+    template test(C: static Curve, body: untyped) =
+      block:
+        proc testInstance() =
+          let Zero {.inject.} = block:
+            var Z{.noInit.}: Fp6[C]
+            Z.setZero()
+            Z
+          let One {.inject.} = block:
+            var O{.noInit.}: Fp6[C]
+            O.setOne()
+            O
+
+          for _ in 0 ..< 1: # Iters:
+            let x {.inject.} = rng.random(Fp6[C])
+            var r{.noinit, inject.}: Fp6[C]
+            body
+
+        testInstance()
+
+    test(BN254):
+      r.prod(x, Zero)
+      check: bool(r == Zero)
+    test(BN254):
+      r.prod(Zero, x)
+      check: bool(r == Zero)
+    test(BN254):
+      r.prod(x, One)
+      check: bool(r == x)
+    test(BN254):
+      r.prod(One, x)
+      check: bool(r == x)
+    test(BLS12_381):
+      r.prod(x, Zero)
+      check: bool(r == Zero)
+    test(BLS12_381):
+      r.prod(Zero, x)
+      check: bool(r == Zero)
+    test(BLS12_381):
+      r.prod(x, One)
+      check: bool(r == x)
+    test(BLS12_381):
+      r.prod(One, x)
+      check: bool(r == x)
+    test(BN462):
+      r.prod(x, Zero)
+      check: bool(r == Zero)
+    test(BN462):
+      r.prod(Zero, x)
+      check: bool(r == Zero)
+    test(BN462):
+      r.prod(x, One)
+      check: bool(r == x)
+    test(BN462):
+      r.prod(One, x)
+      check: bool(r == x)