diff --git a/benchmarks/bench_fp6.nim b/benchmarks/bench_fp6.nim
index 93a3c8a..2acede8 100644
--- a/benchmarks/bench_fp6.nim
+++ b/benchmarks/bench_fp6.nim
@@ -45,7 +45,7 @@ proc main() =
     negBench(Fp6[curve], Iters)
     mulBench(Fp6[curve], Iters)
     sqrBench(Fp6[curve], Iters)
-    # invBench(Fp6[curve], InvIters)
+    invBench(Fp6[curve], InvIters)
     echo "-".repeat(80)
 
 main()
diff --git a/constantine/arithmetic/limbs.nim b/constantine/arithmetic/limbs.nim
index 3b18839..cc764fa 100644
--- a/constantine/arithmetic/limbs.nim
+++ b/constantine/arithmetic/limbs.nim
@@ -48,10 +48,9 @@ debug:
 
   func toString*(a: Limbs): string =
     result = "["
-    result.add $BaseType(a[0]) & " (0x" & toHex(BaseType(a[0])) & ')'
+    result.add " 0x" & toHex(BaseType(a[0]))
     for i in 1 ..< a.len:
-      result.add ", "
-      result.add $BaseType(a[i]) & " (0x" & toHex(BaseType(a[i])) & ')'
+      result.add ", 0x" & toHex(BaseType(a[i]))
     result.add "])"
 
   func toHex*(a: Limbs): string =
diff --git a/constantine/tower_field_extensions/fp2_complex.nim b/constantine/tower_field_extensions/fp2_complex.nim
index 6921799..eb29b6c 100644
--- a/constantine/tower_field_extensions/fp2_complex.nim
+++ b/constantine/tower_field_extensions/fp2_complex.nim
@@ -134,7 +134,7 @@ func prod*(r: var Fp2, a, b: Fp2) =
   r.c1 -= a1b1          # r1 = (b0 + b1)(a0 + a1) - a0b0 - a1b1 # [3 Mul, 2 Add, 3 Sub]
 
 func inv*(r: var Fp2, a: Fp2) =
-  ## Compute the modular multiplicative inverse of ``a``
+  ## Compute the multiplicative inverse of ``a``
   ## in 𝔽p2 = 𝔽p[𝑖]
   #
   # Algorithm: (the inverse exist if a != 0 which might cause constant-time issue)
diff --git a/constantine/tower_field_extensions/fp6_1_plus_i.nim b/constantine/tower_field_extensions/fp6_1_plus_i.nim
index e31c607..af5cc35 100644
--- a/constantine/tower_field_extensions/fp6_1_plus_i.nim
+++ b/constantine/tower_field_extensions/fp6_1_plus_i.nim
@@ -46,7 +46,7 @@ type
   Xi = object
     ## ξ (Xi) the cubic non-residue
 
-func `*`(_: typedesc[Xi], a: Fp2): Fp2 =
+func `*`(_: typedesc[Xi], a: Fp2): Fp2 {.inline.}=
   ## Multiply an element of 𝔽p2 by 𝔽p6 cubic non-residue 1 + 𝑖
   ## (c0 + c1 𝑖) (1 + 𝑖) => c0 + (c0 + c1)𝑖 + c1 𝑖²
   ##                     => c0 - c1 + (c0 + c1) 𝑖
@@ -56,7 +56,7 @@ func `*`(_: typedesc[Xi], a: Fp2): Fp2 =
 template `*`(a: Fp2, _: typedesc[Xi]): Fp2 =
   Xi * a
 
-func `*=`(a: var Fp2, _: typedesc[Xi]) =
+func `*=`(a: var Fp2, _: typedesc[Xi]) {.inline.}=
   ## Inplace multiply an element of 𝔽p2 by 𝔽p6 cubic non-residue 1 + 𝑖
   let t = a.c0
   a.c0 -= a.c1
@@ -122,3 +122,56 @@ func prod*[C](r: var Fp6[C], a, b: Fp6[C]) =
   r.c2 -= v0
   r.c2 -= v2
   r.c2 += v1
+
+func inv*[C](r: var Fp6[C], a: Fp6[C]) =
+  ## Compute the multiplicative inverse of ``a``
+  ## in 𝔽p6 = 𝔽p2[∛(1 + 𝑖)]
+  #
+  # Algorithm 5.23
+  #
+  # Arithmetic of Finite Fields
+  # Chapter 5 of Guide to Pairing-Based Cryptography
+  # Jean Luc Beuchat, Luis J. Dominguez Perez, Sylvain Duquesne, Nadia El Mrabet, Laura Fuentes-Castañeda, Francisco Rodríguez-Henríquez, 2017\
+  # https://www.researchgate.net/publication/319538235_Arithmetic_of_Finite_Fields
+  #
+  # We optimize for stack usage and use 4 temporaries (+r as temporary)
+  # instead of 9, because 5 * 2 (𝔽p2) * Bitsize would be:
+  # - ~2540 bits for BN254
+  # - ~3810 bits for BLS12-381
+  var
+    v1 {.noInit.}, v2 {.noInit.}, v3 {.noInit.}: Fp2[C]
+
+  # A in r0
+  # A <- a0² - ξ(a1 a2)
+  r.c0.square(a.c0)
+  v1.prod(a.c1, a.c2)
+  v1 *= Xi
+  r.c0 -= v1
+
+  # B in v1
+  # B <- ξ a2² - a0 a1
+  v1.square(a.c2)
+  v1 *= Xi
+  v2.prod(a.c0, a.c1)
+  v1 -= v2
+
+  # C in v2
+  # C <- a1² - a0 a2
+  v2.square(a.c1)
+  v3.prod(a.c0, a.c2)
+  v2 -= v3
+
+  # F in v3
+  # F <- ξ a1 C + a0 A + ξ a2 B
+  r.c1.prod(v1, Xi * a.c2)
+  r.c2.prod(v2, Xi * a.c1)
+  v3.prod(r.c0, a.c0)
+  v3 += r.c1
+  v3 += r.c2
+
+  v3.inv(v3)
+
+  # (a0 + a1 ξ + a2 ξ²)^-1 = (A + B ξ + C ξ²) / F
+  r.c0 *= v3
+  r.c1.prod(v1, v3)
+  r.c2.prod(v2, v3)
diff --git a/tests/test_fp6.nim b/tests/test_fp6.nim
index 63c27fe..39fdfd3 100644
--- a/tests/test_fp6.nim
+++ b/tests/test_fp6.nim
@@ -183,7 +183,7 @@ suite "𝔽p6 = 𝔽p2[∛(1+𝑖)] (irreducible polynomial x³ - (1+𝑖))":
             O.setOne()
             O
 
-          for _ in 0 ..< 1: # Iters:
+          for _ in 0 ..< Iters:
             let x {.inject.} = rng.random(Fp6[C])
             var r{.noinit, inject.}: Fp6[C]
             body
@@ -226,3 +226,135 @@ suite "𝔽p6 = 𝔽p2[∛(1+𝑖)] (irreducible polynomial x³ - (1+𝑖))":
     test(BN462):
       r.prod(One, x)
       check: bool(r == x)
+
+  test "𝔽p6 = 𝔽p2[∛(1+𝑖)] addition is associative and commutative":
+    proc abelianGroup(curve: static Curve) =
+      for _ in 0 ..< Iters:
+        let a = rng.random(Fp6[curve])
+        let b = rng.random(Fp6[curve])
+        let c = rng.random(Fp6[curve])
+
+        var tmp1{.noInit.}, tmp2{.noInit.}: Fp6[curve]
+
+        # r0 = (a + b) + c
+        tmp1.sum(a, b)
+        tmp2.sum(tmp1, c)
+        let r0 = tmp2
+
+        # r1 = a + (b + c)
+        tmp1.sum(b, c)
+        tmp2.sum(a, tmp1)
+        let r1 = tmp2
+
+        # r2 = (a + c) + b
+        tmp1.sum(a, c)
+        tmp2.sum(tmp1, b)
+        let r2 = tmp2
+
+        # r3 = a + (c + b)
+        tmp1.sum(c, b)
+        tmp2.sum(a, tmp1)
+        let r3 = tmp2
+
+        # r4 = (c + a) + b
+        tmp1.sum(c, a)
+        tmp2.sum(tmp1, b)
+        let r4 = tmp2
+
+        # ...
+
+        check:
+          bool(r0 == r1)
+          bool(r0 == r2)
+          bool(r0 == r3)
+          bool(r0 == r4)
+
+    abelianGroup(BN254)
+    abelianGroup(P256)
+    abelianGroup(Secp256k1)
+    abelianGroup(BLS12_377)
+    abelianGroup(BLS12_381)
+    abelianGroup(BN446)
+    abelianGroup(FKM12_447)
+    abelianGroup(BLS12_461)
+    abelianGroup(BN462)
+
+  test "𝔽p6 = 𝔽p2[∛(1+𝑖)] multiplication is associative and commutative":
+    proc commutativeRing(curve: static Curve) =
+      for _ in 0 ..< Iters:
+        let a = rng.random(Fp6[curve])
+        let b = rng.random(Fp6[curve])
+        let c = rng.random(Fp6[curve])
+
+        var tmp1{.noInit.}, tmp2{.noInit.}: Fp6[curve]
+
+        # r0 = (a * b) * c
+        tmp1.prod(a, b)
+        tmp2.prod(tmp1, c)
+        let r0 = tmp2
+
+        # r1 = a * (b * c)
+        tmp1.prod(b, c)
+        tmp2.prod(a, tmp1)
+        let r1 = tmp2
+
+        # r2 = (a * c) * b
+        tmp1.prod(a, c)
+        tmp2.prod(tmp1, b)
+        let r2 = tmp2
+
+        # r3 = a * (c * b)
+        tmp1.prod(c, b)
+        tmp2.prod(a, tmp1)
+        let r3 = tmp2
+
+        # r4 = (c * a) * b
+        tmp1.prod(c, a)
+        tmp2.prod(tmp1, b)
+        let r4 = tmp2
+
+        # ...
+
+        check:
+          bool(r0 == r1)
+          bool(r0 == r2)
+          bool(r0 == r3)
+          bool(r0 == r4)
+
+    commutativeRing(BN254)
+    commutativeRing(BLS12_377)
+    commutativeRing(BLS12_381)
+    commutativeRing(BN446)
+    commutativeRing(FKM12_447)
+    commutativeRing(BLS12_461)
+    commutativeRing(BN462)
+
+  test "𝔽p6 = 𝔽p2[∛(1+𝑖)] extension field multiplicative inverse":
+    proc mulInvOne(curve: static Curve) =
+      var one: Fp6[curve]
+      one.setOne()
+
+      block: # Inverse of 1 is 1
+        var r {.noInit.}: Fp6[curve]
+        r.inv(one)
+        check: bool(r == one)
+
+      var aInv, r{.noInit.}: Fp6[curve]
+
+      for _ in 0 ..< 1: # Iters:
+        let a = rng.random(Fp6[curve])
+
+        aInv.inv(a)
+        r.prod(a, aInv)
+        check: bool(r == one)
+
+        r.prod(aInv, a)
+        check: bool(r == one)
+
+    mulInvOne(BN254)
+    mulInvOne(BLS12_377)
+    mulInvOne(BLS12_381)
+    mulInvOne(BN446)
+    mulInvOne(FKM12_447)
+    mulInvOne(BLS12_461)
+    mulInvOne(BN462)