diff --git a/benchmarks/bench_pairing_bls12_381.nim b/benchmarks/bench_pairing_bls12_381.nim
new file mode 100644
index 0000000..0518612
--- /dev/null
+++ b/benchmarks/bench_pairing_bls12_381.nim
@@ -0,0 +1,51 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  # Internals
+  ../constantine/config/curves,
+  ../constantine/arithmetic,
+  ../constantine/towers,
+  # Helpers
+  ../helpers/static_for,
+  ./bench_pairing_template,
+  # Standard library
+  std/strutils
+
+# ############################################################
+#
+#               Benchmark of pairings
+#                   for BLS12-381
+#
+# ############################################################
+
+
+const Iters = 50
+const AvailableCurves = [
+  BLS12_381,
+]
+
+proc main() =
+  separator()
+  staticFor i, 0, AvailableCurves.len:
+    const curve = AvailableCurves[i]
+    lineDoubleBench(curve, Iters)
+    lineAddBench(curve, Iters)
+    mulFp12byLineBench(curve, Iters)
+    separator()
+    finalExpEasyBench(curve, Iters)
+    finalExpHardBLS12Bench(curve, Iters)
+    separator()
+    millerLoopBLS12Bench(curve, Iters)
+    finalExpBLS12Bench(curve, Iters)
+    separator()
+    pairingBLS12Bench(curve, Iters)
+    separator()
+
+main()
+notes()
diff --git a/benchmarks/bench_pairing_template.nim b/benchmarks/bench_pairing_template.nim
new file mode 100644
index 0000000..33be3fe
--- /dev/null
+++ b/benchmarks/bench_pairing_template.nim
@@ -0,0 +1,206 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+# ############################################################
+#
+#             Benchmark of pairings
+#
+# ############################################################
+
+import
+  # Internals
+  ../constantine/config/[curves, common],
+  ../constantine/arithmetic,
+  ../constantine/io/io_bigints,
+  ../constantine/towers,
+  ../constantine/elliptic/[ec_weierstrass_projective, ec_weierstrass_affine],
+  ../constantine/hash_to_curve/cofactors,
+  ../constantine/pairing/[
+    cyclotomic_fp12,
+    lines_projective,
+    mul_fp12_by_lines,
+    pairing_bls12
+  ],
+  # Helpers
+  ../helpers/[prng_unsafe, static_for],
+  ./platforms,
+  # Standard library
+  std/[monotimes, times, strformat, strutils, macros]
+
+var rng: RngState
+let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+rng.seed(seed)
+echo "bench xoshiro512** seed: ", seed
+
+# warmup
+proc warmup*() =
+  # Warmup - make sure cpu is on max perf
+  let start = cpuTime()
+  var foo = 123
+  for i in 0 ..< 300_000_000:
+    foo += i*i mod 456
+    foo = foo mod 789
+
+  # Compiler shouldn't optimize away the results as cpuTime rely on sideeffects
+  let stop = cpuTime()
+  echo &"Warmup: {stop - start:>4.4f} s, result {foo} (displayed to avoid compiler optimizing warmup away)\n"
+
+warmup()
+
+when defined(gcc):
+  echo "\nCompiled with GCC"
+elif defined(clang):
+  echo "\nCompiled with Clang"
+elif defined(vcc):
+  echo "\nCompiled with MSVC"
+elif defined(icc):
+  echo "\nCompiled with ICC"
+else:
+  echo "\nCompiled with an unknown compiler"
+
+echo "Optimization level => "
+echo "  no optimization: ", not defined(release)
+echo "  release: ", defined(release)
+echo "  danger: ", defined(danger)
+echo "  inline assembly: ", UseASM_X86_64
+
+when (sizeof(int) == 4) or defined(Constantine32):
+  echo "⚠️ Warning: using Constantine with 32-bit limbs"
+else:
+  echo "Using Constantine with 64-bit limbs"
+
+when SupportsCPUName:
+  echo "Running on ", cpuName(), ""
+
+when SupportsGetTicks:
+  echo "\n⚠️ Cycles measurements are approximate and use the CPU nominal clock: Turbo-Boost and overclocking will skew them."
+  echo "i.e. a 20% overclock will be about 20% off (assuming no dynamic frequency scaling)"
+
+echo "\n=================================================================================================================\n"
+
+proc separator*() =
+  echo "-".repeat(177)
+
+proc report(op, curve: string, start, stop: MonoTime, startClk, stopClk: int64, iters: int) =
+  let ns = inNanoseconds((stop-start) div iters)
+  let throughput = 1e9 / float64(ns)
+  when SupportsGetTicks:
+    echo &"{op:<60} {curve:<15} {throughput:>15.3f} ops/s     {ns:>9} ns/op     {(stopClk - startClk) div iters:>9} CPU cycles (approx)"
+  else:
+    echo &"{op:<60} {curve:<15} {throughput:>15.3f} ops/s     {ns:>9} ns/op"
+
+proc notes*() =
+  echo "Notes:"
+  echo "  - Compilers:"
+  echo "    Compilers are severely limited on multiprecision arithmetic."
+  echo "    Constantine compile-time assembler is used by default (nimble bench_fp)."
+  echo "    GCC is significantly slower than Clang on multiprecision arithmetic due to catastrophic handling of carries."
+  echo "    GCC also seems to have issues with large temporaries and register spilling."
+  echo "    This is somewhat alleviated by Constantine compile-time assembler."
+  echo "    Bench on specific compiler with assembler: \"nimble bench_ec_g1_gcc\" or \"nimble bench_ec_g1_clang\"."
+  echo "    Bench on specific compiler with assembler: \"nimble bench_ec_g1_gcc_noasm\" or \"nimble bench_ec_g1_clang_noasm\"."
+  echo "  - The simplest operations might be optimized away by the compiler."
+  echo "  - Fast Squaring and Fast Multiplication are possible if there are spare bits in the prime representation (i.e. the prime uses 254 bits out of 256 bits)"
+
+template bench(op: string, C: static Curve, iters: int, body: untyped): untyped =
+  let start = getMonotime()
+  when SupportsGetTicks:
+    let startClk = getTicks()
+  for _ in 0 ..< iters:
+    body
+  when SupportsGetTicks:
+    let stopClk = getTicks()
+  let stop = getMonotime()
+
+  when not SupportsGetTicks:
+    let startClk = -1'i64
+    let stopClk = -1'i64
+
+  report(op, $C, start, stop, startClk, stopClk, iters)
+
+func random_point*(rng: var RngState, EC: typedesc): EC {.noInit.} =
+  result = rng.random_unsafe(EC)
+  result.clearCofactorReference()
+
+proc lineDoubleBench*(C: static Curve, iters: int) =
+  var line: Line[Fp2[C], C.getSexticTwist()]
+  var T = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  let P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+  var Paff: ECP_SWei_Aff[Fp[C]]
+  Paff.affineFromProjective(P)
+  bench("Line double", C, iters):
+    line.line_double(T, Paff)
+
+proc lineAddBench*(C: static Curve, iters: int) =
+  var line: Line[Fp2[C], C.getSexticTwist()]
+  var T = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  let
+    P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+    Q = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  var
+    Paff: ECP_SWei_Aff[Fp[C]]
+    Qaff: ECP_SWei_Aff[Fp2[C]]
+  Paff.affineFromProjective(P)
+  Qaff.affineFromProjective(Q)
+  bench("Line add", C, iters):
+    line.line_add(T, Qaff, Paff)
+
+proc mulFp12byLineBench*(C: static Curve, iters: int) =
+  var line: Line[Fp2[C], C.getSexticTwist()]
+  var T = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  let P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+  var Paff: ECP_SWei_Aff[Fp[C]]
+  Paff.affineFromProjective(P)
+
+  line.line_double(T, Paff)
+  var f = rng.random_unsafe(Fp12[C])
+
+  bench("Mul 𝔽p12 by line xy000z", C, iters):
+    f.mul_sparse_by_line_xy000z(line)
+
+proc millerLoopBLS12Bench*(C: static Curve, iters: int) =
+  let
+    P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+    Q = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  var
+    Paff: ECP_SWei_Aff[Fp[C]]
+    Qaff: ECP_SWei_Aff[Fp2[C]]
+  Paff.affineFromProjective(P)
+  Qaff.affineFromProjective(Q)
+
+  var f: Fp12[C]
+
+  bench("Miller Loop BLS12", C, iters):
+    f.millerLoopGenericBLS12(Paff, Qaff)
+
+proc finalExpEasyBench*(C: static Curve, iters: int) =
+  var r = rng.random_unsafe(Fp12[C])
+  bench("Final Exponentiation Easy", C, iters):
+    r.finalExpEasy()
+
+proc finalExpHardBLS12Bench*(C: static Curve, iters: int) =
+  var r = rng.random_unsafe(Fp12[C])
+  r.finalExpEasy()
+  bench("Final Exponentiation Hard BLS12", C, iters):
+    r.finalExpHard_BLS12()
+
+proc finalExpBLS12Bench*(C: static Curve, iters: int) =
+  var r = rng.random_unsafe(Fp12[C])
+  bench("Final Exponentiation BLS12", C, iters):
+    r.finalExpEasy()
+    r.finalExpHard_BLS12()
+
+proc pairingBLS12Bench*(C: static Curve, iters: int) =
+  let
+    P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+    Q = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+
+  var f: Fp12[C]
+
+  bench("Pairing BLS12", C, iters):
+    f.pairing_bls12(P, Q)
diff --git a/constantine.nimble b/constantine.nimble
index 4c59b62..44b2686 100644
--- a/constantine.nimble
+++ b/constantine.nimble
@@ -71,7 +71,8 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # Edge cases highlighted by past bugs
   ("tests/t_ec_wstrass_prj_edge_cases.nim", false),
   # Pairing
-  ("tests/t_pairing_fp12_sparse.nim", false),
+  ("tests/t_pairing_mul_fp12_by_lines.nim", false),
+  ("tests/t_pairing_cyclotomic_fp12.nim", false),
   ("tests/t_pairing_bn254_nogami_optate.nim", false),
   ("tests/t_pairing_bn254_snarks_optate.nim", false),
   ("tests/t_pairing_bls12_381_optate.nim", false)
@@ -417,7 +418,7 @@ task bench_fp12_gcc_noasm, "Run benchmark 𝔽p12 with gcc - no Assembly":
 task bench_fp12_clang_noasm, "Run benchmark 𝔽p12 with clang - no Assembly":
   runBench("bench_fp12", "clang", useAsm = false)
 
-task bench_ec_g1, "Run benchmark on Elliptic Curve group 𝔾1 - Short Weierstrass with Projective Coordinates - GCC":
+task bench_ec_g1, "Run benchmark on Elliptic Curve group 𝔾1 - Short Weierstrass with Projective Coordinates - Default compiler":
   runBench("bench_ec_g1")
 
 task bench_ec_g1_gcc, "Run benchmark on Elliptic Curve group 𝔾1 - Short Weierstrass with Projective Coordinates - GCC":
@@ -432,7 +433,7 @@ task bench_ec_g1_gcc_noasm, "Run benchmark on Elliptic Curve group 𝔾1 - Short
 task bench_ec_g1_clang_noasm, "Run benchmark on Elliptic Curve group 𝔾1 - Short Weierstrass with Projective Coordinates - Clang no Assembly":
   runBench("bench_ec_g1", "clang", useAsm = false)
 
-task bench_ec_g2, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - GCC":
+task bench_ec_g2, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - Default compiler":
   runBench("bench_ec_g2")
 
 task bench_ec_g2_gcc, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - GCC":
@@ -446,3 +447,18 @@ task bench_ec_g2_gcc_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Short
 
 task bench_ec_g2_clang_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - Clang no Assembly":
   runBench("bench_ec_g2", "clang", useAsm = false)
+
+task bench_pairing_bls12_381, "Run pairings benchmarks for BLS12-381 - Default compiler":
+  runBench("bench_pairing_bls12_381")
+
+task bench_pairing_bls12_381_gcc, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - GCC":
+  runBench("bench_pairing_bls12_381", "gcc")
+
+task bench_pairing_bls12_381_clang, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - Clang":
+  runBench("bench_pairing_bls12_381", "clang")
+
+task bench_pairing_bls12_381_gcc_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - GCC no Assembly":
+  runBench("bench_pairing_bls12_381", "gcc", useAsm = false)
+
+task bench_pairing_bls12_381_clang_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - Clang no Assembly":
+  runBench("bench_pairing_bls12_381", "clang", useAsm = false)
diff --git a/constantine/isogeny/frobenius.nim b/constantine/isogeny/frobenius.nim
index 9036273..2f32102 100644
--- a/constantine/isogeny/frobenius.nim
+++ b/constantine/isogeny/frobenius.nim
@@ -79,28 +79,28 @@ const FrobMapConst_BLS12_381 = [
     "0x1904d3bf02bb0667c231beb4202c0d1f0fd603fd3cbd5f4f7b2443d784bab9c4f67ea53d63e7813d8d0775ed92235fb8",
     "0xfc3e2b36c4e03288e9e902231f9fb854a14787b6c7b36fec0c8ec971f63c5f282d5ac14d6c7ec22cf78a126ddc4af3"
   ),
-  Fp2[BLS12_381].fromHex(
+  Fp2[BLS12_381].fromHex(  # SNR^((p-1)/6)^2 = SNR^((p-1)/3)
     "0x0",
     "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac"
   ),
-  Fp2[BLS12_381].fromHex(
+  Fp2[BLS12_381].fromHex( # SNR^((p-1)/6)^3 = SNR^((p-1)/2)
     "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09",
     "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09"
   ),
-  Fp2[BLS12_381].fromHex(
+  Fp2[BLS12_381].fromHex( # SNR^((p-1)/6)^4 = SNR^(2(p-1)/3)
     "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad",
     "0x0"
   ),
-  Fp2[BLS12_381].fromHex(
+  Fp2[BLS12_381].fromHex( # SNR^((p-1)/6)^5
     "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116",
     "0x144e4211384586c16bd3ad4afa99cc9170df3560e77982d0db45f3536814f0bd5871c1908bd478cd1ee605167ff82995"
   )],
   # frobenius(2)
-  [Fp2[BLS12_381].fromHex(
+  [Fp2[BLS12_381].fromHex( # norm(SNR)^((p-1)/6)^1
     "0x1",
     "0x0"
   ),
-  Fp2[BLS12_381].fromHex(
+  Fp2[BLS12_381].fromHex( # norm(SNR)^((p-1)/6)^2
     "0x5f19672fdf76ce51ba69c6076a0f77eaddb3a93be6f89688de17d813620a00022e01fffffffeffff",
     "0x0"
   ),
@@ -151,15 +151,24 @@ func frobenius_map*(r: var Fp4, a: Fp4, k: static int = 1) {.inline.} =
   ## The p-power frobenius automorphism on 𝔽p4
   r.c0.frobenius_map(a.c0, k)
   r.c1.frobenius_map(a.c1, k)
-  r.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][4-1]
+  r.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][3]
+
+func frobenius_map*(r: var Fp6, a: Fp6, k: static int = 1) {.inline.} =
+  ## Computes a^(p^k)
+  ## The p-power frobenius automorphism on 𝔽p6
+  r.c0.frobenius_map(a.c0, k)
+  r.c1.frobenius_map(a.c1, k)
+  r.c2.frobenius_map(a.c2, k)
+  r.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][2]
+  r.c2.mulCheckSparse FrobMapConst_BLS12_381[k-1][4]
 
 func frobenius_map*(r: var Fp12, a: Fp12, k: static int = 1) {.inline.} =
   ## Computes a^(p^k)
-  ## The p-power frobenius automorphism on 𝔽p4
+  ## The p-power frobenius automorphism on 𝔽p12
   static: doAssert r.c0 is Fp4
   for r_fp4, a_fp4 in fields(r, a):
     for r_fp2, a_fp2 in fields(r_fp4, a_fp4):
-      r_fp2.frobenius_map(a_fp2)
+      r_fp2.frobenius_map(a_fp2, k)
 
   r.c0.c0.mulCheckSparse FrobMapConst_BLS12_381[k-1][0]
   r.c0.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][3]
diff --git a/constantine/pairing/README.md b/constantine/pairing/README.md
index 23e34a3..9a07fa8 100644
--- a/constantine/pairing/README.md
+++ b/constantine/pairing/README.md
@@ -24,7 +24,7 @@
   Scott, Benger, Charlemagne, Perez, Kachisa, 2008\
   https://eprint.iacr.org/2008/490.pdf
 
-- Faster Squaring in the Cyclotomic Subgroup ofSixth Degree Extensions\
+- Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions\
   Granger, Scott, 2009\
   https://eprint.iacr.org/2009/565.pdf
 
@@ -87,6 +87,11 @@
   Aurore Guillevic, 2019\
   https://eprint.iacr.org/2019/1371.pdf
 
+- Improving the computation of the optimal ate pairing
+  for a high security level.
+  Loubna Ghammam, Emmanuel Fouotsa
+  J. Appl. Math. Comput.59, 21–36 (2019)
+
 - Efficient Final Exponentiation
   via Cyclotomic Structure for Pairings
   over Families of Elliptic Curves
diff --git a/constantine/pairing/cyclotomic_fp12.nim b/constantine/pairing/cyclotomic_fp12.nim
new file mode 100644
index 0000000..b41d479
--- /dev/null
+++ b/constantine/pairing/cyclotomic_fp12.nim
@@ -0,0 +1,226 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  ../primitives,
+  ../config/[common, curves],
+  ../arithmetic,
+  ../towers,
+  ../isogeny/frobenius
+
+# ############################################################
+#
+#               Gϕ₁₂, Cyclotomic subgroup of Fp12
+#         with GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) ≡ 1 (mod pⁿ)}
+#
+# ############################################################
+
+# - Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions
+#   Granger, Scott, 2009
+#   https://eprint.iacr.org/2009/565.pdf
+#
+# - On the final exponentiation for calculating
+#   pairings on ordinary elliptic curves
+#   Scott, Benger, Charlemagne, Perez, Kachisa, 2008
+#   https://eprint.iacr.org/2008/490.pdf
+
+# 𝔽p12 -> Gϕ₁₂ - Mapping to Cyclotomic group
+# ----------------------------------------------------------------
+func finalExpEasy*[C: static Curve](f: var Fp12[C]) =
+  ## Easy part of the final exponentiation
+  ##
+  ## This maps the result of the Miller loop into the cyclotomic subgroup Gϕ₁₂
+  ##
+  ## We need to clear the Gₜ cofactor to obtain
+  ## an unique Gₜ representation
+  ## (reminder, Gₜ is a multiplicative group hence we exponentiate by the cofactor)
+  ##
+  ## i.e. Fp^12 --> (fexp easy) --> Gϕ₁₂ --> (fexp hard) --> Gₜ
+  ##
+  ## The final exponentiation is fexp = f^((p^12 - 1) / r)
+  ## It is separated into:
+  ## f^((p^12 - 1) / r) = (p^12 - 1) / ϕ₁₂(p)  * ϕ₁₂(p) / r
+  ##
+  ## with the cyclotomic polynomial ϕ₁₂(p) = (p⁴-p²+1)
+  ##
+  ## With an embedding degree of 12, the easy part of final exponentiation is
+  ##
+  ##  f^(p⁶−1)(p²+1)
+  ##
+  ## And properties are
+  ## 0. f^(p⁶) ≡ conj(f) (mod p^12) for all f in Fp12
+  ##
+  ## After g = f^(p⁶−1) the result g is on the cyclotomic subgroup
+  ## 1. g^(-1) ≡ g^(p⁶) (mod p^12)
+  ## 2. Inversion can be done with conjugate
+  ## 3. g is unitary, its norm |g| (the product of conjugates) is 1
+  ## 4. Squaring has a fast compressed variant.
+  #
+  # Proofs:
+  #
+  # Fp12 can be defined as a quadratic extension over Fp⁶
+  # with g = g₀ + x g₁ with x a quadratic non-residue
+  #
+  # with q = p⁶
+  # The frobenius map f^q ≡ (f₀ + x f₁)^q (mod q²)
+  #                       ≡ f₀^q + x^q f₁^q (mod q²)
+  #                       ≡ f₀ + x^q f₁ (mod q²)
+  #                       ≡ f₀ - x f₁ (mod q²)
+  # hence
+  # f^p⁶ ≡ conj(f) (mod p^12)
+  # Q.E.D. of (0)
+  #
+  # ----------------
+  #
+  # p^12 - 1 = (p⁶−1)(p⁶+1) = (p⁶−1)(p²+1)(p⁴-p²+1)
+  # by Fermat's little theorem we have
+  # f^(p^12 - 1) ≡ 1 (mod p^12)
+  #
+  # Hence f^(p⁶−1)(p⁶+1) ≡ 1 (mod p^12)
+  #
+  # We call g = f^(p⁶−1) we have
+  # g^(p⁶+1) ≡ 1 (mod p^12) <=> g^(p⁶) * g ≡ 1 (mod p^12)
+  # hence g^(-1) ≡ g^(p⁶) (mod p^12)
+  # Q.E.D. of (1)
+  #
+  # --
+  #
+  # From (1) g^(-1) ≡ g^(p⁶) (mod p^12) for g = f^(p⁶−1)
+  # and  (0) f^p⁶ ≡ conj(f) (mod p^12)  for all f in fp12
+  #
+  # so g^(-1) ≡ conj(g) (mod p^12) for g = f^(p⁶−1)
+  # Q.E.D. of (2)
+  #
+  # --
+  #
+  # f^(p^12 - 1) ≡ 1 (mod p^12) by Fermat's Little Theorem
+  # f^(p⁶−1)(p⁶+1) ≡ 1 (mod p^12)
+  # g^(p⁶+1) ≡ 1 (mod p^12)
+  # g * g^p⁶ ≡ 1 (mod p^12)
+  # g * conj(g) ≡ 1 (mod p^12)
+  # Q.E.D. of (3)
+  var g {.noinit.}: typeof(f)
+  g.inv(f)              # g = f^-1
+  conj(f)               # f = f^p⁶
+  g *= f                # g = f^(p⁶-1)
+  f.frobenius_map(g, 2) # f = f^((p⁶-1) p²)
+  f *= g                # f = f^((p⁶-1) (p²+1))
+
+# Gϕ₁₂ - Cyclotomic functions
+# ----------------------------------------------------------------
+# A cyclotomic group is a subgroup of Fp^n defined by
+#
+# GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) = 1}
+#
+# The result of any pairing is in a cyclotomic subgroup
+
+func cyclotomic_inv*(a: var Fp12) =
+  ## Fast inverse for a
+  ## `a` MUST be in the cyclotomic subgroup
+  ## consequently `a` MUST be unitary
+  a.conj()
+
+func cyclotomic_inv*(r: var Fp12, a: Fp12) =
+  ## Fast inverse for a
+  ## `a` MUST be in the cyclotomic subgroup
+  ## consequently `a` MUST be unitary
+  r.conj(a)
+
+func cyclotomic_square*[C](r: var Fp12[C], a: Fp12[C]) =
+  ## Square `a` into `r`
+  ## `a` MUST be in the cyclotomic subgroup
+  ## consequently `a` MUST be unitary
+  #
+  # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions
+  # Granger, Scott, 2009
+  # https://eprint.iacr.org/2009/565.pdf
+
+  when a.c0 is Fp4:
+    # Cubic over quadratic
+    # A = 3a² − 2 ̄a
+    # B = 3 √i c² + 2 ̄b
+    # C = 3b² − 2 ̄c
+    var A{.noinit.}, B{.noinit.}, C{.noinit.}, D{.noinit.}: Fp4[C]
+
+    A = a.c0
+
+    r.c0.square(a.c0)  # r0 = a²
+    D.double(r.c0)     # D  = 2a²
+    r.c0 += D          # r0 = 3a²
+
+    A.conjneg()        # A = − ̄a
+    A.double()         # A = − 2 ̄a
+    r.c0 += A          # r0 = 3a² − 2 ̄a
+
+    B.square(a.c2)     # B = c²
+    B *= NonResidue    # B = √i c²
+    D.double(B)        # B = 2 √i c²
+    B += D             # B = 3 √i c²
+
+    r.c1.conj(a.c1)    # r1 = ̄b
+    r.c1.double()      # r1 = 2 ̄b
+    r.c1 += B          # r1 = 3 √i c² + 2 ̄b
+
+    C.square(a.c1)     # C = b²
+    D.double(C)        # D = 2b²
+    C += D             # C = 3b²
+
+    r.c2.conjneg(a.c2) # r2 = - ̄c
+    r.c2.double()      # r2 = - 2 ̄c
+    r.c2 += C          # r2 = 3b² - 2 ̄c
+
+  else:
+    {.error: "Not implemented".}
+
+func cyclotomic_square*[C](a: var Fp12[C]) =
+  ## Square `a` into `r`
+  ## `a` MUST be in the cyclotomic subgroup
+  ## consequently `a` MUST be unitary
+  #
+  # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions
+  # Granger, Scott, 2009
+  # https://eprint.iacr.org/2009/565.pdf
+
+  when a.c0 is Fp4:
+    # Cubic over quadratic
+    # A = 3a² − 2 ̄a
+    # B = 3 √i c² + 2 ̄b
+    # C = 3b² − 2 ̄c
+    var A{.noinit.}, B{.noinit.}, C{.noinit.}, D{.noinit.}: Fp4[C]
+
+    A = a.c0
+
+    a.c0.square()      # r0 = a²
+    D.double(a.c0)     # D  = 2a²
+    a.c0 += D          # r0 = 3a²
+
+    A.conjneg()        # A = − ̄a
+    A.double()         # A = − 2 ̄a
+    a.c0 += A          # r0 = 3a² − 2 ̄a
+
+    B.square(a.c2)     # B = c²
+    B *= NonResidue    # B = √i c²
+    D.double(B)        # B = 2 √i c²
+    B += D             # B = 3 √i c²
+
+    A = a.c1
+
+    a.c1.conj()        # r1 = ̄b
+    a.c1.double()      # r1 = 2 ̄b
+    a.c1 += B          # r1 = 3 √i c² + 2 ̄b
+
+    C.square(A)        # C = b²
+    D.double(C)        # D = 2b²
+    C += D             # C = 3b²
+
+    a.c2.conjneg()     # r2 = - ̄c
+    a.c2.double()      # r2 = - 2 ̄c
+    a.c2 += C          # r2 = 3b² - 2 ̄c
+
+  else:
+    {.error: "Not implemented".}
diff --git a/constantine/pairing/gt_fp12.nim b/constantine/pairing/mul_fp12_by_lines.nim
similarity index 75%
rename from constantine/pairing/gt_fp12.nim
rename to constantine/pairing/mul_fp12_by_lines.nim
index 1fb9c97..7305248 100644
--- a/constantine/pairing/gt_fp12.nim
+++ b/constantine/pairing/mul_fp12_by_lines.nim
@@ -13,11 +13,11 @@ import
   ../towers,
   ./lines_projective
 
+
 # ############################################################
 #
 #                 Sparse Multiplication
-#                and cyclotomic squaring
-#              for elements of Gₜ = E(Fp¹²)
+#                        by lines
 #
 # ############################################################
 
@@ -38,16 +38,16 @@ import
 #   Craig Costello, Tanja Lange, and Michael Naehrig, 2009
 #   https://eprint.iacr.org/2009/615.pdf
 
-# TODO: we assume an embedding degree k of 12 and a sextic twist.
-#       -> Generalize to KSS (k=18), BLS24 and BLS48 curves
-#
-# TODO: we assume a 2-3-2 towering scheme
-#
-# TODO: merge that in the quadratic/cubic files
-
-# 𝔽p12 - Sparse functions
+# 𝔽p12 by line - Sparse functions
 # ----------------------------------------------------------------
 
+func mul_sparse_by_0y*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
+  ## Sparse multiplication of an Fp4 element
+  ## with coordinates (a₀, a₁) by (0, b₁, 0)
+  r.c0.prod(a.c1, b)
+  r.c0 *= NonResidue
+  r.c1.prod(a.c0, b)
+
 func mul_sparse_by_0y0*[C: static Curve](r: var Fp6[C], a: Fp6[C], b: Fp2[C]) =
   ## Sparse multiplication of an Fp6 element
   ## with coordinates (a₀, a₁, a₂) by (0, b₁, 0)
@@ -60,7 +60,7 @@ func mul_sparse_by_0y0*[C: static Curve](r: var Fp6[C], a: Fp6[C], b: Fp2[C]) =
   # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
   #    = ξ (a1 b1 + a2 b1 - v1)
   #    = ξ a2 b1
-  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + β v2
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2
   #    = a0 b1 + a1 b1 - v1
   #    = a0 b1
   # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
@@ -152,17 +152,47 @@ func mul_sparse_by_line_xy000z*[C: static Curve, Tw: static SexticTwist](
 
   static: doAssert f.c0.typeof is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4"
 
-  var v: Fp12[C]
-  v.c0.c0 = l.x
-  v.c0.c1 = l.y
-  v.c2.c1 = l.z
+  # In the following equations (taken from cubic extension implementation)
+  # a = f
+  # b0 = (x, y)
+  # b1 = (0, 0)
+  # b2 = (0, z)
+  #
+  # v0 = a0 b0 = (f00, f01).(x, y)
+  # v1 = a1 b1 = (f10, f11).(0, 0)
+  # v2 = a2 b2 = (f20, f21).(0, z)
+  #
+  # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
+  #    = ξ (a1 b2 + a2 b2 - v2) + v0
+  #    = ξ a1 b2 + v0
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2
+  #    = a0 b0 + a1 b0 - v0 + ξ v2
+  #    = a1 b0 + ξ v2
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
+  #    = (a0 + a2) * (b0 + b2) - v0 - v2
 
-  f *= v
+  var b0 {.noInit.}, v0{.noInit.}, v2{.noInit.}, t{.noInit.}: Fp4[C]
 
-# Gₜ = 𝔽p12 - Cyclotomic functions
-# ----------------------------------------------------------------
-# A cyclotomic group is a subgroup of Fp^n defined by
-#
-# GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) = 1}
-#
-# The result of any pairing is in a cyclotomic subgroup
+  b0.c0 = l.x
+  b0.c1 = l.y
+
+  v0.prod(f.c0, b0)
+  v2.mul_sparse_by_0y(f.c2, l.z)
+
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2
+  f.c2 += f.c0 # r2 = a0 + a2
+  t = b0
+  t.c1 += l.z  # t = b0 + b2
+  f.c2 *= t    # r2 = (a0 + a2)(b0 + b2)
+  f.c2 -= v0
+  f.c2 -= v2   # r2 = (a0 + a2)(b0 + b2) - v0 - v2
+
+  # r0 = ξ a1 b2 + v0
+  f.c0.mul_sparse_by_0y(f.c1, l.z)
+  f.c0 *= NonResidue
+  f.c0 += v0
+
+  # r1 = a1 b0 + ξ v2
+  f.c1 *= b0
+  v2 *= NonResidue
+  f.c1 += v2
diff --git a/constantine/pairing/pairing_bls12.nim b/constantine/pairing/pairing_bls12.nim
index 12a8079..b0885c7 100644
--- a/constantine/pairing/pairing_bls12.nim
+++ b/constantine/pairing/pairing_bls12.nim
@@ -16,8 +16,10 @@ import
     ec_weierstrass_affine,
     ec_weierstrass_projective
   ],
+  ../isogeny/frobenius,
   ./lines_projective,
-  ./gt_fp12
+  ./mul_fp12_by_lines,
+  ./cyclotomic_fp12
 
 # ############################################################
 #
@@ -33,6 +35,11 @@ import
 #   and Tadanori Teruya, 2020
 #   https://eprint.iacr.org/2020/875.pdf
 #
+# - Improving the computation of the optimal ate pairing
+#   for a high security level.
+#   Loubna Ghammam, Emmanuel Fouotsa
+#   J. Appl. Math. Comput.59, 21–36 (2019)
+#
 # - Faster Pairing Computations on Curves with High-Degree Twists
 #   Craig Costello, Tanja Lange, and Michael Naehrig, 2009
 #   https://eprint.iacr.org/2009/615.pdf
@@ -48,13 +55,16 @@ const BLS12_381_param = block:
 
 const BLS12_381_param_isNeg = true
 
+# Generic slow pairing implementation
+# ----------------------------------------------------------------
+
 const BLS12_381_finalexponent = block:
   # (p^12 - 1) / r
   # BigInt[4314].fromHex"0x2ee1db5dcc825b7e1bda9c0496a1c0a89ee0193d4977b3f7d4507d07363baa13f8d14a917848517badc3a43d1073776ab353f2c30698e8cc7deada9c0aadff5e9cfee9a074e43b9a660835cc872ee83ff3a0f0f1c0ad0d6106feaf4e347aa68ad49466fa927e7bb9375331807a0dce2630d9aa4b113f414386b0e8819328148978e2b0dd39099b86e1ab656d2670d93e4d7acdd350da5359bc73ab61a0c5bf24c374693c49f570bcd2b01f3077ffb10bf24dde41064837f27611212596bc293c8d4c01f25118790f4684d0b9c40a68eb74bb22a40ee7169cdc1041296532fef459f12438dfc8e2886ef965e61a474c5c85b0129127a1b5ad0463434724538411d1676a53b5a62eb34c05739334f46c02c3f0bd0c55d3109cd15948d0a1fad20044ce6ad4c6bec3ec03ef19592004cedd556952c6d8823b19dadd7c2498345c6e5308f1c511291097db60b1749bf9b71a9f9e0100418a3ef0bc627751bbd81367066bca6a4c1b6dcfc5cceb73fc56947a403577dfa9e13c24ea820b09c1d9f7c31759c3635de3f7a3639991708e88adce88177456c49637fd7961be1a4c7e79fb02faa732e2f3ec2bea83d196283313492caa9d4aff1c910e9622d2a73f62537f2701aaef6539314043f7bbce5b78c7869aeb2181a67e49eeed2161daf3f881bd88592d767f67c4717489119226c2f011d4cab803e9d71650a6f80698e2f8491d12191a04406fbc8fbd5f48925f98630e68bfb24c0bcb9b55df57510"
   # (p^12 - 1) / r * 3
   BigInt[4316].fromHex"0x8ca592196587127a538fd40dc3e541f9dca04bb7dc671be77cf17715a2b2fe3bea73dfb468d8f473094aecb7315a664019fbd84913caba6579c08fd42009fe1bd6fcbce15eacb2cf3218a165958cb8bfdae2d2d54207282314fc0dea9d6ff3a07dbd34efb77b732ba5f994816e296a72928cfee133bdc3ca9412b984b9783d9c6aa81297ab1cd294a502304773528bbae8706979f28efa0d355b0224e2513d6e4a5d3bb4dde0523678105d9167ff1323d6e99ac312d8a7d762336370c4347bb5a7e405d6f3496b2dd38e722d4c1f3ac25e3167ec2cb543d69430c37c2f98fcdd0dd36caa9f5aa7994cec31b24ed5e515911037b376e521070d29c9d56cfa8c3574363efb20f28c19e4105ab99edd44084bd23725017931d6740bda71e5f07600ce6b407e543c4bc40bcd4c0b600e6c98003bf8548986b14d9098746dc89d154af91ad54f337b31c79222145dd3ed254fdeda0300c49ebcd2352765f533883a3513435f3ee452496f5166c25bf503bd6ec0a0679efda3b46ebf86211d458de749460d4a2a19abe6ea2accb451ab9a096b98465d044dc2a7f86c253a4ee57b6df108eff598a8dbc483bf8b74c2789939db85ffd7e0fd55b32bc26877f5be26fa7d750500ce2fab93c0cbe7336b126a5693d0c16484f37addccc7642590dbe98538990b88637e374d545d9b34b67448d0357e60280bbd8542f1f4e813caa8e8db57364b4e0cc14f35af381dd9b71ec9292b3a3f16e42362d2019e05f30"
 
-func millerLoopGenericBLS12[C: static Curve](
+func millerLoopGenericBLS12*[C: static Curve](
        f: var Fp12[C],
        P: ECP_SWei_Aff[Fp[C]],
        Q: ECP_SWei_Aff[Fp2[C]]
@@ -140,13 +150,104 @@ func pairing_bls12_reference*[C](gt: var Fp12[C], P: ECP_SWei_Proj[Fp[C]], Q: EC
   gt.millerLoopGenericBLS12(Paff, Qaff)
   gt.finalExpGeneric()
 
-func finalExpEasy[C: static Curve](f: var Fp12[C]) =
-  ## Easy part of the final exponentiation
-  ## We need to clear the GT cofactor to obtain
-  ## an unique GT representation
-  ## (reminder, GT is a multiplicative group hence we exponentiate by the cofactor)
+# Optimized pairing implementation
+# ----------------------------------------------------------------
+
+func cycl_sqr_repeated(f: var Fp12, num: int) =
+  ## Repeated cyclotomic squarings
+  for _ in 0 ..< num:
+    f.cyclotomic_square()
+
+func pow_xdiv2(r: var Fp12[BLS12_381], a: Fp12[BLS12_381], invert = BLS12_381_param_isNeg) =
+  ## f^(x/2) with x the curve parameter
+  ## For BLS12_381 f^-0xd201000000010000
+
+  r.cyclotomic_square(a)
+  r *= a
+  r.cycl_sqr_repeated(2)
+  r *= a
+  r.cycl_sqr_repeated(3)
+  r *= a
+  r.cycl_sqr_repeated(9)
+  r *= a
+  r.cycl_sqr_repeated(32)   # TODO: use Karabina?
+  r *= a
+  r.cycl_sqr_repeated(16-1) # Don't do the last iteration
+
+  if invert:
+    r.cyclotomic_inv()
+
+func pow_x(r: var Fp12[BLS12_381], a: Fp12[BLS12_381], invert = BLS12_381_param_isNeg) =
+  ## f^x with x the curve parameter
+  ## For BLS12_381 f^-0xd201000000010000
+  r.pow_xdiv2(a, invert)
+  r.cyclotomic_square()
+
+func finalExpHard_BLS12*[C: static Curve](f: var Fp12[C]) =
+  ## Hard part of the final exponentiation
+  ## Specialized for BLS12 curves
   ##
-  ## With an embedding degree of 12, the easy part of final exponentiation is
+  # - Efficient Final Exponentiation
+  #   via Cyclotomic Structure for Pairings
+  #   over Families of Elliptic Curves
+  #   Daiki Hayashida and Kenichiro Hayasaka
+  #   and Tadanori Teruya, 2020
+  #   https://eprint.iacr.org/2020/875.pdf
+  #
+  # p14: 3 Φ₁₂(p(x))/r(x) = (x−1)² (x+p) (x²+p²−1) + 3
+  #
+  # TODO: paper costs are 4Eₓ+Eₓ/₂+7M₁₂+S₁₂+F₁+F₂
+  #       so we have an extra cyclotomic squaring since we use 5Eₓ
+  #
+  # with
+  # - Eₓ being f^x
+  # - Eₓ/₂ being f^(x/2)
+  # - M₁₂ being mul in Fp12
+  # - S₁₂ being cyclotomic squaring
+  # - Fₙ being n Frobenius applications
+
+  var v0 {.noInit.}, v1 {.noInit.}, v2 {.noInit.}: Fp12[C]
+
+  # Save for f³ and (x−1)²
+  v2.cyclotomic_square(f)      # v2 = f²
+
+  # (x−1)²
+  v0.pow_xdiv2(v2)             # v0 = (f²)^(x/2) = f^x
+  v1.cyclotomic_inv(f)         # v1 = f^-1
+  v0 *= v1                     # v0 = f^(x-1)
+  v1.pow_x(v0)                 # v1 = (f^(x-1))^x
+  v0.cyclotomic_inv()          # v0 = (f^(x-1))^-1
+  v0 *= v1                     # v0 = (f^(x-1))^(x-1) = f^((x-1)*(x-1)) = f^((x-1)²)
+
+  # (x+p)
+  v1.pow_x(v0)                 # v1 = f^((x-1)².x)
+  v0.frobenius_map(v0)         # v0 = f^((x-1)².p)
+  v0 *= v1                     # v0 = f^((x-1)².(x+p))
+
+  # + 3
+  f *= v2                      # f = f³
+
+  # (x²+p²−1)
+  v2.pow_x(v0, invert = false)
+  v1.pow_x(v2, invert = false) # v1 = f^((x-1)².(x+p).x²)
+  v2.frobenius_map(v0, 2)      # v2 = f^((x-1)².(x+p).p²)
+  v0.cyclotomic_inv()          # v0 = f^((x-1)².(x+p).-1)
+  v0 *= v1                     # v0 = f^((x-1)².(x+p).(x²-1))
+  v0 *= v2                     # v0 = f^((x-1)².(x+p).(x²+p²-1))
+
+  # (x−1)².(x+p).(x²+p²−1) + 3
+  f *= v0
+
+func pairing_bls12*[C](gt: var Fp12[C], P: ECP_SWei_Proj[Fp[C]], Q: ECP_SWei_Proj[Fp2[C]]) =
+  ## Compute the optimal Ate Pairing for BLS12 curves
+  ## Input: P ∈ G1, Q ∈ G2
+  ## Output: e(P, Q) ∈ Gt
   ##
-  ##  f^(p⁶−1)(p²+1)
-  discard
+  ## Reference implementation
+  var Paff {.noInit.}: ECP_SWei_Aff[Fp[C]]
+  var Qaff {.noInit.}: ECP_SWei_Aff[Fp2[C]]
+  Paff.affineFromProjective(P)
+  Qaff.affineFromProjective(Q)
+  gt.millerLoopGenericBLS12(Paff, Qaff)
+  gt.finalExpEasy()
+  gt.finalExpHard_BLS12()
diff --git a/constantine/pairing/pairing_bn.nim b/constantine/pairing/pairing_bn.nim
index 5e928a2..989af37 100644
--- a/constantine/pairing/pairing_bn.nim
+++ b/constantine/pairing/pairing_bn.nim
@@ -18,7 +18,8 @@ import
     ec_weierstrass_projective
   ],
   ./lines_projective,
-  ./gt_fp12,
+  ./mul_fp12_by_lines,
+  ./cyclotomic_fp12,
   ../isogeny/frobenius
 
 # ############################################################
diff --git a/constantine/tower_field_extensions/cubic_extensions.nim b/constantine/tower_field_extensions/cubic_extensions.nim
index c61bda1..409bf6f 100644
--- a/constantine/tower_field_extensions/cubic_extensions.nim
+++ b/constantine/tower_field_extensions/cubic_extensions.nim
@@ -214,6 +214,13 @@ func conj*(a: var CubicExt) {.inline.} =
   a.c1.conjneg()
   a.c2.conj()
 
+func conj*(r: var CubicExt, a: CubicExt) {.inline.} =
+  ## Computes the conjugate out-of-place
+  mixin conj, conjneg
+  r.c0.conj(a.c0)
+  r.c1.conjneg(a.c1)
+  r.c2.conj(a.c2)
+
 func square*(a: var CubicExt) {.inline.} =
   ## In-place squaring
   let t = a
diff --git a/constantine/tower_field_extensions/quadratic_extensions.nim b/constantine/tower_field_extensions/quadratic_extensions.nim
index d004ba4..d6abcf5 100644
--- a/constantine/tower_field_extensions/quadratic_extensions.nim
+++ b/constantine/tower_field_extensions/quadratic_extensions.nim
@@ -259,6 +259,11 @@ func conjneg*(a: var QuadraticExt) {.inline.} =
   ## Computes the negated conjugate in-place
   a.c0.neg()
 
+func conjneg*(r: var QuadraticExt, a: QuadraticExt) {.inline.} =
+  ## Computes the negated conjugate out-of-place
+  r.c0.neg(a.c0)
+  r.c1 = a.c1
+
 func square*(r: var QuadraticExt, a: QuadraticExt) {.inline.} =
   mixin fromComplexExtension
   when r.fromComplexExtension():
diff --git a/tests/t_fp6_frobenius.nim b/tests/t_fp6_frobenius.nim
new file mode 100644
index 0000000..68efef1
--- /dev/null
+++ b/tests/t_fp6_frobenius.nim
@@ -0,0 +1,33 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  # Internals
+  ../constantine/towers,
+  ../constantine/config/curves,
+  # Test utilities
+  ./t_fp_tower_frobenius_template
+
+const TestCurves = [
+    # BN254_Nogami
+    # BN254_Snarks,
+    # BLS12_377,
+    BLS12_381,
+    # BN446
+    # FKM12_447
+    # BLS12_461
+    # BN462
+  ]
+
+runFrobeniusTowerTests(
+  ExtDegree =6,
+  Iters = 8,
+  TestCurves = TestCurves,
+  moduleName = "test_fp6_frobenius",
+  testSuiteDesc = "𝔽p6 Frobenius map: Frobenius(a, k) = a^(p^k) (mod p⁶)"
+)
diff --git a/tests/t_pairing_bls12_381_line_functions.nim b/tests/t_pairing_bls12_381_line_functions.nim
index bf4edc0..ee465c0 100644
--- a/tests/t_pairing_bls12_381_line_functions.nim
+++ b/tests/t_pairing_bls12_381_line_functions.nim
@@ -65,7 +65,7 @@ func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen)
     else:
       result = rng.random_long01Seq_with_randZ(EC)
 
-suite "Pairing - Line Functions on BLS12-381":
+suite "Pairing - Line Functions on BLS12-381" & " [" & $WordBitwidth & "-bit mode]":
   test "Line double - lt,t(P)":
     proc test_line_double(C: static Curve, randZ: bool, gen: RandomGen) =
       for _ in 0 ..< Iters:
diff --git a/tests/t_pairing_bls12_381_optate.nim b/tests/t_pairing_bls12_381_optate.nim
index 6b824f8..522432e 100644
--- a/tests/t_pairing_bls12_381_optate.nim
+++ b/tests/t_pairing_bls12_381_optate.nim
@@ -7,9 +7,11 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.
 
 import
+  ../constantine/config/common,
   ../constantine/config/curves,
   ../constantine/pairing/pairing_bls12,
   # Test utilities
   ./t_pairing_template
 
-runPairingTests(4, BLS12_381, pairing_bls12_reference)
+# runPairingTests(4, BLS12_381, pairing_bls12_reference)
+runPairingTests(4, BLS12_381, pairing_bls12)
diff --git a/tests/t_pairing_bn254_nogami_optate.nim b/tests/t_pairing_bn254_nogami_optate.nim
index a32b939..207c9ed 100644
--- a/tests/t_pairing_bn254_nogami_optate.nim
+++ b/tests/t_pairing_bn254_nogami_optate.nim
@@ -7,6 +7,7 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.
 
 import
+  ../constantine/config/common,
   ../constantine/config/curves,
   ../constantine/pairing/pairing_bn,
   # Test utilities
diff --git a/tests/t_pairing_bn254_snarks_optate.nim b/tests/t_pairing_bn254_snarks_optate.nim
index 07946d0..b42bb6f 100644
--- a/tests/t_pairing_bn254_snarks_optate.nim
+++ b/tests/t_pairing_bn254_snarks_optate.nim
@@ -7,6 +7,7 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.
 
 import
+  ../constantine/config/common,
   ../constantine/config/curves,
   ../constantine/pairing/pairing_bn,
   # Test utilities
diff --git a/tests/t_pairing_cyclotomic_fp12.nim b/tests/t_pairing_cyclotomic_fp12.nim
new file mode 100644
index 0000000..cf1c9fb
--- /dev/null
+++ b/tests/t_pairing_cyclotomic_fp12.nim
@@ -0,0 +1,127 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  # Standard library
+  std/[tables, unittest, times],
+  # Internals
+  ../constantine/config/common,
+  ../constantine/[arithmetic, primitives],
+  ../constantine/towers,
+  ../constantine/config/curves,
+  ../constantine/io/io_towers,
+  ../constantine/pairing/cyclotomic_fp12,
+  ../constantine/isogeny/frobenius,
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+const
+  Iters = 4
+  TestCurves = [
+    # BN254_Nogami,
+    # BN254_Snarks,
+    # BLS12_377,
+    BLS12_381
+  ]
+
+type
+  RandomGen = enum
+    Uniform
+    HighHammingWeight
+    Long01Sequence
+
+var rng: RngState
+let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+rng.seed(seed)
+echo "\n------------------------------------------------------\n"
+echo "test_pairing_fp12_sparse xoshiro512** seed: ", seed
+
+func random_elem(rng: var RngState, F: typedesc, gen: RandomGen): F {.inline, noInit.} =
+  if gen == Uniform:
+    result = rng.random_unsafe(F)
+  elif gen == HighHammingWeight:
+    result = rng.random_highHammingWeight(F)
+  else:
+    result = rng.random_long01Seq(F)
+
+suite "Pairing - Cyclotomic subgroup - GΦ₁₂(p) = {α ∈ Fp¹² : α^Φ₁₂(p) ≡ 1 (mod p¹²)}" & " [" & $WordBitwidth & "-bit mode]":
+  test "Easy part of the final exponentiation maps to the cyclotomic subgroup":
+    proc test_final_exp_easy_cycl(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        var f = rng.random_elem(Fp12[C], gen)
+
+        f.finalExpEasy()
+
+        var f4, minus_f2: typeof(f)
+        minus_f2.frobenius_map(f, 2)  # f^p²
+        f4.frobenius_map(minus_f2, 2) # f^p⁴
+        minus_f2.conj()               # f^⁻²p
+
+        f *= f4
+        f *= minus_f2                 # f^(p⁴-p²+1) = f^Φ₁₂(p)
+
+        check: bool(f.isOne())
+
+    staticFor(curve, TestCurves):
+      test_final_exp_easy_cycl(curve, gen = Uniform)
+      test_final_exp_easy_cycl(curve, gen = HighHammingWeight)
+      test_final_exp_easy_cycl(curve, gen = Long01Sequence)
+
+  test "Cyclotomic inverse":
+    proc test_cycl_inverse(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        var f = rng.random_elem(Fp12[C], gen)
+
+        f.finalExpEasy()
+        var g = f
+
+        f.cyclotomic_inv()
+        f *= g
+
+        check: bool(f.isOne())
+
+    staticFor(curve, TestCurves):
+      test_cycl_inverse(curve, gen = Uniform)
+      test_cycl_inverse(curve, gen = HighHammingWeight)
+      test_cycl_inverse(curve, gen = Long01Sequence)
+
+  test "Cyclotomic squaring":
+    proc test_cycl_squaring_in_place(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        var f = rng.random_elem(Fp12[C], gen)
+
+        f.finalExpEasy()
+        var g = f
+
+        f.square()
+        g.cyclotomic_square()
+
+        check: bool(f == g)
+
+    staticFor(curve, TestCurves):
+      test_cycl_squaring_in_place(curve, gen = Uniform)
+      test_cycl_squaring_in_place(curve, gen = HighHammingWeight)
+      test_cycl_squaring_in_place(curve, gen = Long01Sequence)
+
+    proc test_cycl_squaring_out_place(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        var f = rng.random_elem(Fp12[C], gen)
+
+        f.finalExpEasy()
+        var g = f
+        var r: typeof(f)
+
+        f.square()
+        r.cyclotomic_square(g)
+
+        check: bool(f == r)
+
+    staticFor(curve, TestCurves):
+      test_cycl_squaring_out_place(curve, gen = Uniform)
+      test_cycl_squaring_out_place(curve, gen = HighHammingWeight)
+      test_cycl_squaring_out_place(curve, gen = Long01Sequence)
diff --git a/tests/t_pairing_fp12_sparse.nim b/tests/t_pairing_mul_fp12_by_lines.nim
similarity index 72%
rename from tests/t_pairing_fp12_sparse.nim
rename to tests/t_pairing_mul_fp12_by_lines.nim
index 671d050..c7a731d 100644
--- a/tests/t_pairing_fp12_sparse.nim
+++ b/tests/t_pairing_mul_fp12_by_lines.nim
@@ -17,7 +17,7 @@ import
   ../constantine/io/io_towers,
   ../constantine/pairing/[
     lines_projective,
-    gt_fp12
+    mul_fp12_by_lines
   ],
   # Test utilities
   ../helpers/[prng_unsafe, static_for]
@@ -50,6 +50,25 @@ func random_elem(rng: var RngState, F: typedesc, gen: RandomGen): F {.inline, no
     result = rng.random_long01Seq(F)
 
 suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent with dense 𝔽p12 mul":
+  test "Dense 𝔽p4 by Sparse 0y":
+    proc test_fp4_0y(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        let a = rng.random_elem(Fp4[C], gen)
+        let y = rng.random_elem(Fp2[C], gen)
+        let b = Fp4[C](c1: y)
+
+        var r {.noInit.}, r2 {.noInit.}: Fp4[C]
+
+        r.prod(a, b)
+        r2.mul_sparse_by_0y(a, y)
+
+        check: bool(r == r2)
+
+    staticFor(curve, TestCurves):
+      test_fp4_0y(curve, gen = Uniform)
+      test_fp4_0y(curve, gen = HighHammingWeight)
+      test_fp4_0y(curve, gen = Long01Sequence)
+
   test "Dense 𝔽p6 by Sparse 0y0":
     proc test_fp6_0y0(C: static Curve, gen: static RandomGen) =
       for _ in 0 ..< Iters:
@@ -91,7 +110,7 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
       test_fp6_0y0(curve, gen = Long01Sequence)
 
   when Fp12[BN254_Snarks]().c0.typeof is Fp6:
-    test "Sparse 𝔽p12 resulting from xy00z0 line function":
+    test "Sparse 𝔽p12/𝔽p6 resulting from xy00z0 line function":
       proc test_fp12_xy00z0(C: static Curve, gen: static RandomGen) =
         for _ in 0 ..< Iters:
           var a = rng.random_elem(Fp12[C], gen)
@@ -117,7 +136,7 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
         test_fp12_xy00z0(curve, gen = HighHammingWeight)
         test_fp12_xy00z0(curve, gen = Long01Sequence)
 
-    test "Sparse 𝔽p12 resulting from xyz000 line function":
+    test "Sparse 𝔽p12/𝔽p6 resulting from xyz000 line function":
       proc test_fp12_xyz000(C: static Curve, gen: static RandomGen) =
         for _ in 0 ..< Iters:
           var a = rng.random_elem(Fp12[C], gen)
@@ -141,3 +160,32 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
         test_fp12_xyz000(curve, gen = Uniform)
         test_fp12_xyz000(curve, gen = HighHammingWeight)
         test_fp12_xyz000(curve, gen = Long01Sequence)
+  else:
+    static: doAssert Fp12[BN254_Snarks]().c0.typeof is Fp4
+
+    test "Sparse 𝔽p12/𝔽p4 resulting from xy000z line function":
+      proc test_fp12_xy000z(C: static Curve, gen: static RandomGen) =
+        for _ in 0 ..< Iters:
+          var a = rng.random_elem(Fp12[C], gen)
+          var a2 = a
+
+          var x = rng.random_elem(Fp2[C], gen)
+          var y = rng.random_elem(Fp2[C], gen)
+          var z = rng.random_elem(Fp2[C], gen)
+
+          let line = Line[Fp2[C], Mtwist](x: x, y: y, z: z)
+          let b = Fp12[C](
+            c0: Fp4[C](c0: x, c1: y),
+            # c1
+            c2: Fp4[C](       c1: z),
+          )
+
+          a *= b
+          a2.mul_sparse_by_line_xy000z(line)
+
+          check: bool(a == a2)
+
+      staticFor(curve, TestCurves):
+        test_fp12_xy000z(curve, gen = Uniform)
+        test_fp12_xy000z(curve, gen = HighHammingWeight)
+        test_fp12_xy000z(curve, gen = Long01Sequence)
diff --git a/tests/t_pairing_template.nim b/tests/t_pairing_template.nim
index 39155d0..d1b2785 100644
--- a/tests/t_pairing_template.nim
+++ b/tests/t_pairing_template.nim
@@ -10,6 +10,7 @@ import
   # Standard library
   std/unittest, times,
   # Internals
+  ../constantine/config/common,
   ../constantine/[arithmetic, primitives],
   ../constantine/towers,
   ../constantine/config/curves,
@@ -75,14 +76,13 @@ template runPairingTests*(Iters: static int, C: static Curve, pairing_fn: untype
       r2.pairing_fn(P2, Q)
       r3.pairing_fn(P, Q2)
 
-      check:
-        bool(not r.isZero())
-        bool(not r.isOne())
-        bool(r == r2)
-        bool(r == r3)
-        bool(r2 == r3)
+      doAssert bool(not r.isZero())
+      doAssert bool(not r.isOne())
+      doAssert bool(r == r2)
+      doAssert bool(r == r3)
+      doAssert bool(r2 == r3)
 
-  suite "Pairing - Optimal Ate on " & $C:
+  suite "Pairing - Optimal Ate on " & $C & " [" & $WordBitwidth & "-bit mode]":
     test "Bilinearity e([2]P, Q) = e(P, [2]Q) = e(P, Q)^2":
       test_bilinearity_double_impl(randZ = false, gen = Uniform)
       test_bilinearity_double_impl(randZ = true, gen = Uniform)