From d8721a1a131cdad2af452dae6c97b05eec64c814 Mon Sep 17 00:00:00 2001
From: Mamy Ratsimbazafy <mamy_github@numforge.co>
Date: Sun, 14 Jul 2024 23:31:47 +0000
Subject: [PATCH] =?UTF-8?q?=F0=9D=94=BE=E2=82=9C=20exponentiation,=20with?=
 =?UTF-8?q?=20endomorphism=20acceleration=20(#429)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* feat(𝔾ₜ exponentiation): Add square-multiply and signed recoding and wNAF 𝔾ₜ exponentiation

* feat(𝔾ₜ exponentiation): Add andomorphism acceleration + wNAF

* feat(𝔾ₜ exponentiation): Add benchmarks

* fix(bench): forgot declaring var in refactoring
---
 README-PERFORMANCE.md                         |   4 +-
 README.md                                     |   6 +-
 benchmarks/bench_fields_template.nim          |   8 +-
 benchmarks/bench_fp.nim                       |   2 +-
 benchmarks/bench_gt.nim                       |  64 ++++
 benchmarks/bench_gt_template.nim              | 167 ++++++++++
 constantine.nimble                            |  28 +-
 .../math/elliptic/ec_endomorphism_accel.nim   |  49 +--
 .../math/elliptic/ec_multi_scalar_mul.nim     |   4 +-
 .../elliptic/ec_multi_scalar_mul_parallel.nim |   4 +-
 .../math/elliptic/ec_scalar_mul_vartime.nim   |  10 +-
 constantine/math/extension_fields/towers.nim  |   3 -
 .../pairings/gt_exponentiations_vartime.nim   | 312 ++++++++++++++++++
 constantine/named/zoo_endomorphisms.nim       |  19 +-
 sage/derive_endomorphisms.sage                |   6 +-
 sage/derive_frobenius.sage                    |   4 +-
 .../t_pairing_bls12_381_gt_exp.nim            |  15 +
 .../t_pairing_bn254_snarks_gt_exp.nim         |  15 +
 tests/math_pairings/t_pairing_template.nim    |  80 ++++-
 19 files changed, 732 insertions(+), 68 deletions(-)
 create mode 100644 benchmarks/bench_gt.nim
 create mode 100644 benchmarks/bench_gt_template.nim
 create mode 100644 constantine/math/pairings/gt_exponentiations_vartime.nim
 create mode 100644 tests/math_pairings/t_pairing_bls12_381_gt_exp.nim
 create mode 100644 tests/math_pairings/t_pairing_bn254_snarks_gt_exp.nim

diff --git a/README-PERFORMANCE.md b/README-PERFORMANCE.md
index 2e6eb9d7..85f3c443 100644
--- a/README-PERFORMANCE.md
+++ b/README-PERFORMANCE.md
@@ -58,7 +58,7 @@ The full list of benchmarks is available in the [`benchmarks`](./benchmarks) fol
 
 As mentioned in the [Compiler caveats](#compiler-caveats) section, GCC is up to 2x slower than Clang due to mishandling of carries and register usage.
 
-#### Ethereum BLS signatures (over BLS12-381 G2)
+#### Ethereum BLS signatures (over BLS12-381 𝔾₂)
 
 ![Bench Ethereum BLS signature](./media/ethereum_bls_signatures.png)
 
@@ -192,4 +192,4 @@ Constantine does not use heap allocation.
 At the moment Constantine is optimized for 32-bit and 64-bit CPUs.
 
 When performance and code size conflicts, a careful and informed default is chosen.
-In the future, a compile-time flag that goes beyond the compiler `-Os` might be provided.
\ No newline at end of file
+In the future, a compile-time flag that goes beyond the compiler `-Os` might be provided.
diff --git a/README.md b/README.md
index d5ea9dcf..8b731eb3 100644
--- a/README.md
+++ b/README.md
@@ -90,8 +90,8 @@ For all elliptic curves, the following arithmetic is supported
     - on Fr (i.e. modulo the 255-bit curve order)
     - on Fp (i.e. modulo the 381-bit prime modulus)
   - elliptic curve arithmetic:
-    - on elliptic curve over Fp (EC G1) with affine, jacobian and homogenous projective coordinates
-    - on elliptic curve over Fp2 (EC G2) with affine, jacobian and homogenous projective coordinates
+    - on elliptic curve over Fp (EC 𝔾₁) with affine, jacobian and homogenous projective coordinates
+    - on elliptic curve over Fp2 (EC 𝔾₂) with affine, jacobian and homogenous projective coordinates
     - including scalar multiplication, multi-scalar-multiplication (MSM) and parallel MSM
 
 _All operations are constant-time unless explicitly mentioned_ vartime.
@@ -222,7 +222,7 @@ and modify Constantine's [`build.rs`](https://github.com/mratsim/constantine/blo
     ```
     > [!IMPORTANT]
     > Constantine uses a separate modfile for tests.<br />It has no dependencies (key to avoid supply chain attacks) except for testing.
-    
+
 ### From C
 
 1. Install a C compiler, `clang` is recommended, for example:
diff --git a/benchmarks/bench_fields_template.nim b/benchmarks/bench_fields_template.nim
index ca9c414b..5ad318e2 100644
--- a/benchmarks/bench_fields_template.nim
+++ b/benchmarks/bench_fields_template.nim
@@ -229,13 +229,13 @@ proc sqrtRatioVartimeBench*(T: typedesc, iters: int) =
 proc powBench*(T: typedesc, iters: int) =
   let x = rng.random_unsafe(T)
   let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
+  var r = x
   bench("Exp curve order (constant-time) - " & $exponent.bits & "-bit", T, iters):
-    var r = x
     r.pow(exponent)
 
-proc powUnsafeBench*(T: typedesc, iters: int) =
+proc powVartimeBench*(T: typedesc, iters: int) =
   let x = rng.random_unsafe(T)
   let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
-  bench("Exp curve order (Leak exponent bits) - " & $exponent.bits & "-bit", T, iters):
-    var r = x
+  var r = x
+  bench("Exp by curve order (vartime) - " & $exponent.bits & "-bit", T, iters):
     r.pow_vartime(exponent)
diff --git a/benchmarks/bench_fp.nim b/benchmarks/bench_fp.nim
index 9867e3cb..7678340b 100644
--- a/benchmarks/bench_fp.nim
+++ b/benchmarks/bench_fp.nim
@@ -70,7 +70,7 @@ proc main() =
       sqrtRatioVartimeBench(Fp[curve], ExponentIters)
     # Exponentiation by a "secret" of size ~the curve order
     powBench(Fp[curve], ExponentIters)
-    powUnsafeBench(Fp[curve], ExponentIters)
+    powVartimeBench(Fp[curve], ExponentIters)
     separator()
 
 main()
diff --git a/benchmarks/bench_gt.nim b/benchmarks/bench_gt.nim
new file mode 100644
index 00000000..aeb4d93d
--- /dev/null
+++ b/benchmarks/bench_gt.nim
@@ -0,0 +1,64 @@
+# 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/named/algebras,
+  constantine/math/extension_fields,
+  # Helpers
+  ./bench_gt_template
+
+# ############################################################
+#
+#               Benchmark of the 𝔾ₜ group of
+#                  Pairing Friendly curves
+#
+# ############################################################
+
+const Iters = 10000
+const ExpIters = 1000
+const AvailableCurves = [
+  # BN254_Nogami,
+  BN254_Snarks,
+  # BLS12_377,
+  BLS12_381,
+]
+
+proc main() =
+  separator()
+  staticFor i, 0, AvailableCurves.len:
+    const curve = AvailableCurves[i]
+    const bits = Fr[curve].bits()
+    separator()
+    mulBench(Fp12[curve], Iters)
+    sqrBench(Fp12[curve], Iters)
+    invBench(Fp12[curve], Iters)
+    separator()
+    cyclotomicSquare_Bench(Fp12[curve], Iters)
+    cyclotomicInv_Bench(Fp12[curve], Iters)
+    cyclotomicSquareCompressed_Bench(Fp12[curve], Iters)
+    cyclotomicDecompression_Bench(Fp12[curve], Iters)
+    separator()
+    powVartimeBench(Fp12[curve], window = 2, ExpIters)
+    powVartimeBench(Fp12[curve], window = 3, ExpIters)
+    powVartimeBench(Fp12[curve], window = 4, ExpIters)
+    separator()
+    gtExp_sqrmul_vartimeBench(Fp12[curve], ExpIters)
+    gtExp_minHammingWeight_vartimeBench(Fp12[curve], ExpIters)
+    separator()
+    gtExp_wNAF_vartimeBench(Fp12[curve], window = 2, ExpIters)
+    gtExp_wNAF_vartimeBench(Fp12[curve], window = 3, ExpIters)
+    gtExp_wNAF_vartimeBench(Fp12[curve], window = 4, ExpIters)
+    separator()
+    gtExp_endo_wNAF_vartimeBench(Fp12[curve], window = 2, ExpIters)
+    gtExp_endo_wNAF_vartimeBench(Fp12[curve], window = 3, ExpIters)
+    gtExp_endo_wNAF_vartimeBench(Fp12[curve], window = 4, ExpIters)
+    separator()
+
+main()
+notes()
diff --git a/benchmarks/bench_gt_template.nim b/benchmarks/bench_gt_template.nim
new file mode 100644
index 00000000..f27420a9
--- /dev/null
+++ b/benchmarks/bench_gt_template.nim
@@ -0,0 +1,167 @@
+# 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.
+
+# ############################################################
+#
+#             Summary of the performance of a curve
+#
+# ############################################################
+
+import
+  # Internals
+  constantine/platforms/abstractions,
+  constantine/named/algebras,
+  constantine/math/[arithmetic, extension_fields],
+  constantine/math/pairings/[
+    pairings_generic,
+    cyclotomic_subgroups,
+    gt_exponentiations_vartime
+  ],
+  # Helpers
+  helpers/prng_unsafe,
+  ./bench_blueprint
+
+
+export notes
+export abstractions
+proc separator*() = separator(168)
+
+proc report(op, domain: 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:<68} {domain:<20} {throughput:>15.3f} ops/s     {ns:>9} ns/op     {(stopClk - startClk) div iters:>9} CPU cycles (approx)"
+  else:
+    echo &"{op:<68} {domain:<20} {throughput:>15.3f} ops/s     {ns:>9} ns/op"
+
+macro fixFieldDisplay(T: typedesc): untyped =
+  # At compile-time, enums are integers and their display is buggy
+  # we get the Curve ID instead of the curve name.
+  let instantiated = T.getTypeInst()
+  var name = $instantiated[1][0] # Fp
+  name.add "[" & $Algebra(instantiated[1][1].intVal) & "]"
+  result = newLit name
+
+func fixDisplay(T: typedesc): string =
+  when T is (Fp or Fp2 or Fp4 or Fp6 or Fp12):
+    fixFieldDisplay(T)
+  else:
+    $T
+
+func fixDisplay(T: Algebra): string =
+  $T
+
+template bench(op: string, T: typed, iters: int, body: untyped): untyped =
+  measure(iters, startTime, stopTime, startClk, stopClk, body)
+  report(op, fixDisplay(T), startTime, stopTime, startClk, stopClk, iters)
+
+func random_gt*(rng: var RngState, F: typedesc): F {.inline, noInit.} =
+  result = rng.random_unsafe(F)
+  result.finalExp()
+
+proc mulBench*(T: typedesc, iters: int) =
+  var r: T
+  let x = rng.random_gt(T)
+  let y = rng.random_gt(T)
+  preventOptimAway(r)
+  bench("Multiplication", T, iters):
+    r.prod(x, y)
+
+proc sqrBench*(T: typedesc, iters: int) =
+  var r: T
+  let x = rng.random_gt(T)
+  preventOptimAway(r)
+  bench("Squaring", T, iters):
+    r.square(x)
+
+proc invBench*(T: typedesc, iters: int) =
+  var r: T
+  let x = rng.random_gt(T)
+  preventOptimAway(r)
+  bench("Inversion", T, iters):
+    r.inv(x)
+
+proc cyclotomicSquare_Bench*(T: typedesc, iters: int) =
+  var f = rng.random_gt(T)
+
+  bench("Squaring in cyclotomic subgroup", T, iters):
+    f.cyclotomic_square()
+
+proc cyclotomicInv_Bench*(T: typedesc, iters: int) =
+  var f = rng.random_gt(T)
+
+  bench("Inversion in cyclotomic subgroup", T, iters):
+    f.cyclotomic_inv()
+
+proc cyclotomicSquareCompressed_Bench*(T: typedesc, iters: int) =
+  var f = rng.random_gt(T)
+
+  when T is Fp12:
+    type F = Fp2[T.Name]
+  else:
+    {.error: "Only compression of Fp12 extension is configured".}
+
+  var g: G2345[F]
+  g.fromFpk(f)
+
+  bench("Cyclotomic Compressed Squaring", T, iters):
+    g.cyclotomic_square_compressed()
+
+proc cyclotomicDecompression_Bench*(T: typedesc, iters: int) =
+  var f = rng.random_gt(T)
+
+  when T is Fp12:
+    type F = Fp2[T.Name]
+  else:
+    {.error: "Only compression of Fp12 extension is configured".}
+
+  var gs: array[1, G2345[F]]
+  gs[0].fromFpk(f)
+
+  var g1s_ratio: array[1, tuple[g1_num, g1_den: F]]
+  var g0s, g1s: array[1, F]
+
+  bench("Cyclotomic Decompression", T, iters):
+    recover_g1(g1s_ratio[0].g1_num, g1s_ratio[0].g1_den, gs[0])
+    g1s.batch_ratio_g1s(g1s_ratio)
+    g0s[0].recover_g0(g1s[0], gs[0])
+
+proc powVartimeBench*(T: typedesc, window: static int, iters: int) =
+  let x = rng.random_gt(T)
+  let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
+  var r = x
+  bench("Field Exponentiation " & $exponent.bits & "-bit (window-" & $window & ", vartime)", T, iters):
+    r.pow_vartime(exponent, window)
+
+proc gtExp_sqrmul_vartimeBench*(T: typedesc, iters: int) =
+  let x = rng.random_gt(T)
+  let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
+  var r {.noInit.}: T
+  bench("𝔾ₜ Exponentiation " & $exponent.bits & "-bit (cyclotomic square-multiply, vartime)", T, iters):
+    r.gtExp_sqrmul_vartime(x, exponent)
+
+proc gtExp_minHammingWeight_vartimeBench*(T: typedesc, iters: int) =
+  let x = rng.random_gt(T)
+  let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
+  var r {.noInit.}: T
+  bench("𝔾ₜ Exponentiation " & $exponent.bits & "-bit (signed recoding, vartime)", T, iters):
+    r.gtExp_minHammingWeight_vartime(x, exponent)
+
+proc gtExp_wNAF_vartimeBench*(T: typedesc, window: static int, iters: int) =
+  let x = rng.random_gt(T)
+  let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
+  var r {.noInit.}: T
+  bench("𝔾ₜ Exponentiation " & $exponent.bits & "-bit (wNAF-" & $window & ", vartime)", T, iters):
+    r.gtExp_minHammingWeight_windowed_vartime(x, exponent, window)
+
+proc gtExp_endo_wNAF_vartimeBench*(T: typedesc, window: static int, iters: int) =
+  let x = rng.random_gt(T)
+  let exponent = rng.random_unsafe(BigInt[Fr[T.Name].bits()])
+  var r {.noInit.}: T
+  bench("𝔾ₜ Exponentiation " & $exponent.bits & "-bit (endomorphism, wNAF-" & $window & ", vartime)", T, iters):
+    r.gtExpEndo_minHammingWeight_windowed_vartime(x, exponent, window)
diff --git a/constantine.nimble b/constantine.nimble
index dbdccb3f..4ad7e001 100644
--- a/constantine.nimble
+++ b/constantine.nimble
@@ -430,7 +430,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # ----------------------------------------------------------
   ("tests/math_elliptic_curves/t_ec_conversion.nim", false),
 
-  # Elliptic curve arithmetic G1
+  # Elliptic curve arithmetic 𝔾₁
   # ----------------------------------------------------------
   ("tests/math_elliptic_curves/t_ec_shortw_prj_g1_add_double.nim", false),
   # ("tests/math_elliptic_curves/t_ec_shortw_prj_g1_mul_sanity.nim", false),
@@ -458,7 +458,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/math_elliptic_curves/t_ec_twedwards_mul_endomorphism_bandersnatch", false),
 
 
-  # Elliptic curve arithmetic G2
+  # Elliptic curve arithmetic 𝔾₂
   # ----------------------------------------------------------
   # ("tests/math_elliptic_curves/t_ec_shortw_prj_g2_add_double_bn254_snarks.nim", false),
   # ("tests/math_elliptic_curves/t_ec_shortw_prj_g2_mul_sanity_bn254_snarks.nim", false),
@@ -551,7 +551,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # ("tests/math_pairings/t_pairing_bls12_381_gt_subgroup.nim", false),
   # ("tests/math_pairings/t_pairing_bw6_761_gt_subgroup.nim", false),
 
-  # Pairing
+  # Pairing &
   # ----------------------------------------------------------
   # ("tests/math_pairings/t_pairing_bls12_377_line_functions.nim", false),
   # ("tests/math_pairings/t_pairing_bls12_381_line_functions.nim", false),
@@ -562,6 +562,9 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # ("tests/math_pairings/t_pairing_bls12_377_optate.nim", false),
   # ("tests/math_pairings/t_pairing_bls12_381_optate.nim", false),
 
+  ("tests/math_pairings/t_pairing_bn254_snarks_gt_exp.nim", false),
+  ("tests/math_pairings/t_pairing_bls12_381_gt_exp.nim", false),
+
   # Multi-Pairing
   # ----------------------------------------------------------
   ("tests/math_pairings/t_pairing_bn254_nogami_multi.nim", false),
@@ -650,6 +653,7 @@ const benchDesc = [
   "bench_pairing_bls12_381",
   "bench_pairing_bn254_nogami",
   "bench_pairing_bn254_snarks",
+  "bench_gt",
   "bench_summary_bls12_377",
   "bench_summary_bls12_381",
   "bench_summary_bn254_nogami",
@@ -977,25 +981,25 @@ task bench_fp6, "Run benchmark 𝔽p6 with your CC compiler":
 task bench_fp12, "Run benchmark 𝔽p12 with your CC compiler":
   runBench("bench_fp12")
 
-# Elliptic curve G1
+# Elliptic curve 𝔾₁
 # ------------------------------------------
 
 task bench_ec_g1, "Run benchmark on Elliptic Curve group 𝔾1 - CC compiler":
   runBench("bench_ec_g1")
 
-# Elliptic curve G1 - batch operations
+# Elliptic curve 𝔾₁ - batch operations
 # ------------------------------------------
 
 task bench_ec_g1_batch, "Run benchmark on Elliptic Curve group 𝔾1 (batch ops) - CC compiler":
   runBench("bench_ec_g1_batch")
 
-# Elliptic curve G1 - scalar multiplication
+# Elliptic curve 𝔾₁ - scalar multiplication
 # ------------------------------------------
 
 task bench_ec_g1_scalar_mul, "Run benchmark on Elliptic Curve group 𝔾1 (Scalar Multiplication) - CC compiler":
   runBench("bench_ec_g1_scalar_mul")
 
-# Elliptic curve G1 - Multi-scalar-mul
+# Elliptic curve 𝔾₁ - Multi-scalar-mul
 # ------------------------------------------
 
 task bench_ec_msm_pasta, "Run benchmark: Multi-Scalar-Mul for Pasta curves - CC compiler":
@@ -1014,18 +1018,24 @@ task bench_ec_msm_bandersnatch, "Run benchmark: Multi-Scalar-Mul for Bandersnatc
   runBench("bench_ec_msm_bandersnatch")
 
 
-# Elliptic curve G2
+# Elliptic curve 𝔾₂
 # ------------------------------------------
 
 task bench_ec_g2, "Run benchmark on Elliptic Curve group 𝔾2 - CC compiler":
   runBench("bench_ec_g2")
 
-# Elliptic curve G2 - scalar multiplication
+# Elliptic curve 𝔾₂ - scalar multiplication
 # ------------------------------------------
 
 task bench_ec_g2_scalar_mul, "Run benchmark on Elliptic Curve group 𝔾2 (Multi-Scalar-Mul) - CC compiler":
   runBench("bench_ec_g2_scalar_mul")
 
+# 𝔾ₜ
+# ------------------------------------------
+
+task bench_gt, "Run 𝔾ₜ benchmarks - CC compiler":
+  runBench("bench_gt")
+
 # Pairings
 # ------------------------------------------
 
diff --git a/constantine/math/elliptic/ec_endomorphism_accel.nim b/constantine/math/elliptic/ec_endomorphism_accel.nim
index bb518a88..d4be13b5 100644
--- a/constantine/math/elliptic/ec_endomorphism_accel.nim
+++ b/constantine/math/elliptic/ec_endomorphism_accel.nim
@@ -44,11 +44,12 @@ type
 
 template decomposeEndoImpl[scalBits: static int](
        scalar: BigInt[scalBits],
-       F: typedesc[Fp or Fp2],
+       frBits: static int,
+       Name: static Algebra,
+       G: static Subgroup,
        copyMiniScalarsResult: untyped) =
   static: doAssert scalBits >= L, "Cannot decompose a scalar smaller than a mini-scalar or the decomposition coefficient"
   # Equal when no window or no negative handling, greater otherwise
-  const frBits = Fr[F.Name].bits()
   static: doAssert frBits >= scalBits
   static: doAssert L >= ceilDiv_vartime(frBits, M) + 1
   const w = frBits.wordsRequired()
@@ -59,24 +60,24 @@ template decomposeEndoImpl[scalBits: static int](
 
   when M == 2:
     var alphas{.noInit, inject.}: (
-      BigInt[frBits + babai(F)[0][0].bits],
-      BigInt[frBits + babai(F)[1][0].bits]
+      BigInt[frBits + babai(Name, G)[0][0].bits],
+      BigInt[frBits + babai(Name, G)[1][0].bits]
     )
   elif M == 4:
     var alphas{.noInit, inject.}: (
-      BigInt[frBits + babai(F)[0][0].bits],
-      BigInt[frBits + babai(F)[1][0].bits],
-      BigInt[frBits + babai(F)[2][0].bits],
-      BigInt[frBits + babai(F)[3][0].bits]
+      BigInt[frBits + babai(Name, G)[0][0].bits],
+      BigInt[frBits + babai(Name, G)[1][0].bits],
+      BigInt[frBits + babai(Name, G)[2][0].bits],
+      BigInt[frBits + babai(Name, G)[3][0].bits]
     )
   else:
     {.error: "The decomposition degree " & $M & " is not configured".}
 
   staticFor i, 0, M:
-    when bool babai(F)[i][0].isZero():
+    when bool babai(Name, G)[i][0].isZero():
       alphas[i].setZero()
     else:
-      alphas[i].prod_high_words(babai(F)[i][0], scalar, w)
+      alphas[i].prod_high_words(babai(Name, G)[i][0], scalar, w)
 
   # We have k0 = s - 𝛼0 b00 - 𝛼1 b10 ... - 𝛼m bm0
   # and     kj = 0 - 𝛼j b0j - 𝛼1 b1j ... - 𝛼m bmj
@@ -86,13 +87,13 @@ template decomposeEndoImpl[scalBits: static int](
   k[0].copyTruncatedFrom(scalar)
   staticFor miniScalarIdx, 0, M:
     staticFor basisIdx, 0, M:
-      when not bool lattice(F)[basisIdx][miniScalarIdx][0].isZero():
-        when bool lattice(F)[basisIdx][miniScalarIdx][0].isOne():
+      when not bool lattice(Name, G)[basisIdx][miniScalarIdx][0].isZero():
+        when bool lattice(Name, G)[basisIdx][miniScalarIdx][0].isOne():
           alphaB.copyTruncatedFrom(alphas[basisIdx])
         else:
-          alphaB.prod(alphas[basisIdx], lattice(F)[basisIdx][miniScalarIdx][0])
+          alphaB.prod(alphas[basisIdx], lattice(Name, G)[basisIdx][miniScalarIdx][0])
 
-        when lattice(F)[basisIdx][miniScalarIdx][1] xor babai(F)[basisIdx][1]:
+        when lattice(Name, G)[basisIdx][miniScalarIdx][1] xor babai(Name, G)[basisIdx][1]:
           k[miniScalarIdx] += alphaB
         else:
           k[miniScalarIdx] -= alphaB
@@ -103,7 +104,9 @@ func decomposeEndo*[M, scalBits, L: static int](
        miniScalars: var MultiScalar[M, L],
        negatePoints: var array[M, SecretBool],
        scalar: BigInt[scalBits],
-       F: typedesc[Fp or Fp2]) =
+       frBits: static int,
+       Name: static Algebra,
+       G: static Subgroup) =
   ## Decompose a secret scalar into M mini-scalars
   ## using a curve endomorphism(s) characteristics.
   ##
@@ -119,7 +122,7 @@ func decomposeEndo*[M, scalBits, L: static int](
   ##     and negate it as well.
   ##
   ## This implements solution 1.
-  decomposeEndoImpl(scalar, F):
+  decomposeEndoImpl(scalar, frBits, Name, G):
     # Negative miniscalars are turned positive
     # Caller should negate the corresponding Elliptic Curve points
     let isNeg = k[miniScalarIdx].isMsbSet()
@@ -130,7 +133,9 @@ func decomposeEndo*[M, scalBits, L: static int](
 func decomposeEndo*[M, scalBits, L: static int](
        miniScalars: var MultiScalar[M, L],
        scalar: BigInt[scalBits],
-       F: typedesc[Fp or Fp2]) =
+       frBits: static int,
+       Name: static Algebra,
+       G: static Subgroup) =
   ## Decompose a secret scalar into M mini-scalars
   ## using a curve endomorphism(s) characteristics.
   ##
@@ -150,7 +155,7 @@ func decomposeEndo*[M, scalBits, L: static int](
   ## Also for partitioned GLV-SAC (with 8-way decomposition) it is necessary.
   ##
   ## This implements solution 2.
-  decomposeEndoImpl(scalar, F):
+  decomposeEndoImpl(scalar, frBits, Name, G):
     miniScalars[miniScalarIdx].copyTruncatedFrom(k[miniScalarIdx])
 
 # Secret scalar + dynamic point
@@ -339,6 +344,8 @@ func scalarMulEndo*[scalBits; EC](
   const M = when P.F is Fp:  2
             elif P.F is Fp2: 4
             else: {.error: "Unconfigured".}
+  const G = when EC isnot EC_ShortW_Aff|EC_ShortW_Jac|EC_ShortW_Prj: G1
+            else: EC.G
 
   var endos {.noInit.}: array[M-1, EC]
   endos.computeEndomorphisms(P)
@@ -347,7 +354,7 @@ func scalarMulEndo*[scalBits; EC](
   const L = EC.getScalarField().bits().ceilDiv_vartime(M) + 1
   var miniScalars {.noInit.}: array[M, BigInt[L]]
   var negatePoints {.noInit.}: array[M, SecretBool]
-  miniScalars.decomposeEndo(negatePoints, scalar, P.F)
+  miniScalars.decomposeEndo(negatePoints, scalar, EC.getScalarField().bits(), EC.getName(), G)
 
   # 3. Handle negative mini-scalars
   # A scalar decomposition might lead to negative miniscalar.
@@ -510,6 +517,8 @@ func scalarMulGLV_m2w2*[scalBits; EC](P0: var EC, scalar: BigInt[scalBits]) {.me
   mixin affine
   const C = P0.F.Name # curve
   static: doAssert: scalBits <= EC.getScalarField().bits()
+  const G = when EC isnot EC_ShortW_Aff|EC_ShortW_Jac|EC_ShortW_Prj: G1
+            else: EC.G
 
   # 1. Compute endomorphisms
   var P1 {.noInit.}: EC
@@ -519,7 +528,7 @@ func scalarMulGLV_m2w2*[scalBits; EC](P0: var EC, scalar: BigInt[scalBits]) {.me
   const L = computeRecodedLength(EC.getScalarField().bits(), 2)
   var miniScalars {.noInit.}: array[2, BigInt[L]]
   var negatePoints {.noInit.}: array[2, SecretBool]
-  miniScalars.decomposeEndo(negatePoints, scalar, P0.F)
+  miniScalars.decomposeEndo(negatePoints, scalar, EC.getScalarField().bits(), EC.getName(), G)
 
   # 3. Handle negative mini-scalars
   #    Either negate the associated base and the scalar (in the `endomorphisms` array)
diff --git a/constantine/math/elliptic/ec_multi_scalar_mul.nim b/constantine/math/elliptic/ec_multi_scalar_mul.nim
index a9541601..cee871f4 100644
--- a/constantine/math/elliptic/ec_multi_scalar_mul.nim
+++ b/constantine/math/elliptic/ec_multi_scalar_mul.nim
@@ -400,6 +400,8 @@ proc applyEndomorphism[bits: static int, ECaff](
   const M = when ECaff.F is Fp:  2
             elif ECaff.F is Fp2: 4
             else: {.error: "Unconfigured".}
+  const G = when ECaff isnot EC_ShortW_Aff: G1
+            else: ECaff.G
 
   const L = ECaff.getScalarField().bits().ceilDiv_vartime(M) + 1
   let splitCoefs   = allocHeapArray(array[M, BigInt[L]], N)
@@ -407,7 +409,7 @@ proc applyEndomorphism[bits: static int, ECaff](
 
   for i in 0 ..< N:
     var negatePoints {.noinit.}: array[M, SecretBool]
-    splitCoefs[i].decomposeEndo(negatePoints, coefs[i], ECaff.F)
+    splitCoefs[i].decomposeEndo(negatePoints, coefs[i], ECaff.getScalarField().bits(), ECaff.getName(), G)
     if negatePoints[0].bool:
       endoBasis[i][0].neg(points[i])
     else:
diff --git a/constantine/math/elliptic/ec_multi_scalar_mul_parallel.nim b/constantine/math/elliptic/ec_multi_scalar_mul_parallel.nim
index a132caff..4dcbaea4 100644
--- a/constantine/math/elliptic/ec_multi_scalar_mul_parallel.nim
+++ b/constantine/math/elliptic/ec_multi_scalar_mul_parallel.nim
@@ -462,6 +462,8 @@ proc applyEndomorphism_parallel[bits: static int, ECaff](
   const M = when ECaff.F is Fp:  2
             elif ECaff.F is Fp2: 4
             else: {.error: "Unconfigured".}
+  const G = when ECaff isnot EC_ShortW_Aff: G1
+            else: ECaff.G
 
   const L = ECaff.getScalarField().bits().ceilDiv_vartime(M) + 1
   let splitCoefs   = allocHeapArray(array[M, BigInt[L]], N)
@@ -472,7 +474,7 @@ proc applyEndomorphism_parallel[bits: static int, ECaff](
       captures: {coefs, points, splitCoefs, endoBasis}
 
       var negatePoints {.noinit.}: array[M, SecretBool]
-      splitCoefs[i].decomposeEndo(negatePoints, coefs[i], ECaff.F)
+      splitCoefs[i].decomposeEndo(negatePoints, coefs[i], ECaff.getScalarField().bits(), ECaff.getName(), G)
       if negatePoints[0].bool:
         endoBasis[i][0].neg(points[i])
       else:
diff --git a/constantine/math/elliptic/ec_scalar_mul_vartime.nim b/constantine/math/elliptic/ec_scalar_mul_vartime.nim
index 825a4394..7486d39f 100644
--- a/constantine/math/elliptic/ec_scalar_mul_vartime.nim
+++ b/constantine/math/elliptic/ec_scalar_mul_vartime.nim
@@ -180,7 +180,7 @@ func initNAF[precompSize, NafMax: static int, EC, ECaff](
     P.fromAffine(tab[digit shr 1])
     return true
   elif digit < 0:
-    P.fromAffine(tab[digit shr 1])
+    P.fromAffine(tab[-digit shr 1])
     P.neg()
     return true
   else:
@@ -213,7 +213,7 @@ func scalarMul_minHammingWeight_windowed_vartime*[EC](P: var EC, scalar: BigInt,
   # Odd-only divides precomputation table size by another 2
 
   const precompSize = 1 shl (window - 2)
-  static: doAssert window < 8, "Window is too large and precomputation would use " & $(precompSize * sizeof(EC)) & " stack space."
+  static: doAssert window < 8, "Window of size " & $window & " is too large and precomputation would use " & $(precompSize * sizeof(EC)) & " stack space."
 
   var tabEC {.noinit.}: array[precompSize, EC]
   var P2{.noInit.}: EC
@@ -252,12 +252,14 @@ func scalarMulEndo_minHammingWeight_windowed_vartime*[scalBits: static int; EC](
   # Signed digits divides precomputation table size by 2
   # Odd-only divides precomputation table size by another 2
   const precompSize = 1 shl (window - 2)
-  static: doAssert window < 8, "Window is too large and precomputation would use " & $(precompSize * sizeof(EC)) & " stack space."
+  static: doAssert window < 8, "Window of size " & $window & " is too large and precomputation would use " & $(precompSize * sizeof(EC)) & " stack space."
 
   # 1. Compute endomorphisms
   const M = when P.F is Fp:  2
             elif P.F is Fp2: 4
             else: {.error: "Unconfigured".}
+  const G = when EC isnot EC_ShortW_Aff|EC_ShortW_Jac|EC_ShortW_Prj: G1
+            else: EC.G
 
   var endos {.noInit.}: array[M-1, EC]
   endos.computeEndomorphisms(P)
@@ -266,7 +268,7 @@ func scalarMulEndo_minHammingWeight_windowed_vartime*[scalBits: static int; EC](
   const L = EC.getScalarField().bits().ceilDiv_vartime(M) + 1
   var miniScalars {.noInit.}: array[M, BigInt[L]]
   var negatePoints {.noInit.}: array[M, SecretBool]
-  miniScalars.decomposeEndo(negatePoints, scalar, EC.F)
+  miniScalars.decomposeEndo(negatePoints, scalar, EC.getScalarField().bits(), EC.getName(), G)
 
   # 3. Handle negative mini-scalars
   if negatePoints[0].bool:
diff --git a/constantine/math/extension_fields/towers.nim b/constantine/math/extension_fields/towers.nim
index b70fa255..1d76c379 100644
--- a/constantine/math/extension_fields/towers.nim
+++ b/constantine/math/extension_fields/towers.nim
@@ -329,9 +329,6 @@ func has2extraBits*(E: type ExtensionField): bool =
   ## We construct extensions only on Fp (and not Fr)
   getSpareBits(Fp[E.F.Name]) >= 2
 
-template A(E: type ExtensionField2x): Algebra =
-  E.F.Name
-
 template c0*(a: ExtensionField2x): auto =
   a.coords[0]
 template c1*(a: ExtensionField2x): auto =
diff --git a/constantine/math/pairings/gt_exponentiations_vartime.nim b/constantine/math/pairings/gt_exponentiations_vartime.nim
new file mode 100644
index 00000000..cbc9fb31
--- /dev/null
+++ b/constantine/math/pairings/gt_exponentiations_vartime.nim
@@ -0,0 +1,312 @@
+# 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/math/elliptic/ec_endomorphism_accel,
+  constantine/math/arithmetic,
+  constantine/math/extension_fields,
+  constantine/math/io/io_bigints,
+  constantine/platforms/abstractions,
+  constantine/math_arbitrary_precision/arithmetic/limbs_views,
+  constantine/named/zoo_endomorphisms,
+  constantine/named/algebras,
+  ./cyclotomic_subgroups
+
+from constantine/math/elliptic/ec_shortweierstrass_affine import G2
+
+{.push raises: [].} # No exceptions allowed in core cryptographic operations
+{.push checks: off.} # No defects due to array bound checking or signed integer overflow allowed
+
+iterator unpackBE(scalarByte: byte): bool =
+  for i in countdown(7, 0):
+    yield bool((scalarByte shr i) and 1)
+
+func gtExp_sqrmul_vartime*[Gt: ExtensionField](r: var Gt, a: Gt, scalar: BigInt) {.tags:[VarTime], meter.} =
+  ## **Variable-time** Exponentiation in 𝔾ₜ
+  ##
+  ##   r <- aᵏ
+  ##
+  ## This uses the square-and-multiply algorithm
+  ## This MUST NOT be used with secret data.
+  ##
+  ## This is highly VULNERABLE to timing attacks and power analysis attacks.
+  var scalarCanonical: array[scalar.bits.ceilDiv_vartime(8), byte]
+  scalarCanonical.marshal(scalar, bigEndian)
+
+  let a {.noInit.} = a # Avoid aliasing issues
+
+  r.setOne()
+  var isNeutral = true
+
+  for scalarByte in scalarCanonical:
+    for bit in unpackBE(scalarByte):
+      if not isNeutral:
+        r.cyclotomic_square()
+      if bit:
+        if isNeutral:
+          r = a
+          isNeutral = false
+        else:
+          r *= a
+
+func gtExp_addchain_4bit_vartime[Gt: ExtensionField](r: var Gt, a: Gt, scalar: BigInt) {.tags:[VarTime], meter.} =
+  ## **Variable-time** Exponentiation in 𝔾ₜ
+  ## This can only handle for small scalars up to 2⁴ = 16 excluded
+  let s = uint scalar.limbs[0]
+
+  case s
+  of 0:
+    r.setNeutral()
+  of 1:
+    discard
+  of 2:
+    r.cyclotomic_square(a)
+  of 3:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    r.prod(a, t)
+  of 4:
+    r.cyclotomic_square(a)
+    r.cyclotomic_square()
+  of 5:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    t.cyclotomic_square()
+    r.prod(a, t)
+  of 6:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    r.prod(a, t)
+    r.cyclotomic_square()
+  of 7:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    t.cyclotomic_square()
+    t.cyclotomic_square()
+    r.cyclotomic_inv(a)
+    r *= t
+  of 8:
+    r.cyclotomic_square(a)
+    r.cyclotomic_square()
+    r.cyclotomic_square()
+  of 9:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    t.cyclotomic_square()
+    t.cyclotomic_square()
+    r.prod(a, t)
+  of 10:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    t.cyclotomic_square()
+    r.prod(a, t)
+    r.cyclotomic_square()
+  of 11:
+    var t1 {.noInit.}, t2 {.noInit.}: Gt
+    t1.cyclotomic_square(a)  # [2]P
+    t2.cyclotomic_square(t1)
+    t2.cyclotomic_square()   # [8]P
+    t1 *= t2
+    r.prod(a, t1)
+  of 12:
+    var t1 {.noInit.}, t2 {.noInit.}: Gt
+    t1.cyclotomic_square(a)
+    t1.cyclotomic_square()   # [4]P
+    t2.cyclotomic_square(t1) # [8]P
+    r.prod(t1, t2)
+  of 13:
+    var t1 {.noInit.}, t2 {.noInit.}: Gt
+    t1.cyclotomic_square(a)
+    t1.cyclotomic_square()   # [4]P
+    t2.cyclotomic_square(t1) # [8]P
+    t1 *= t2
+    r.prod(a, t1)
+  of 14:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    t.cyclotomic_square()
+    t.cyclotomic_square()
+    r.cyclotomic_inv(a)
+    t *= r  # [7]P
+    r.cyclotomic_square(t)
+  of 15:
+    var t {.noInit.}: Gt
+    t.cyclotomic_square(a)
+    t.cyclotomic_square()
+    t.cyclotomic_square()
+    t.cyclotomic_square()
+    r.cyclotomic_inv(a)
+    r *= t
+  else:
+    unreachable()
+
+func gtExp_minHammingWeight_vartime*[Gt: ExtensionField](r: var Gt, a: Gt, scalar: BigInt) {.tags:[VarTime].}  =
+  ## **Variable-time** Exponentiation in 𝔾ₜ
+  ##
+  ##   r <- aᵏ
+  ##
+  ## This uses an online recoding with minimum Hamming Weight
+  ## (which is not NAF, NAF is least-significant bit to most)
+  ## This MUST NOT be used with secret data.
+  ##
+  ## This is highly VULNERABLE to timing attacks and power analysis attacks
+  let a {.noInit.} = a # Avoid aliasing issues
+  var na {.noInit.}: Gt
+  na.cyclotomic_inv(a)
+
+  r.setOne()
+  for bit in recoding_l2r_signed_vartime(scalar):
+    r.cyclotomic_square()
+    if bit == 1:
+      r *= a
+    elif bit == -1:
+      r *= na
+
+func initNAF[precompSize, NafMax: static int, Gt: ExtensionField](
+       acc: var Gt,
+       tab: array[precompSize, Gt],
+       naf: array[NafMax, int8], nafLen: int,
+       nafIteratorIdx: int): bool {.inline.} =
+
+  let digit = naf[nafLen-1-nafIteratorIdx]
+  if digit > 0:
+    acc = tab[digit shr 1]
+    return true
+  elif digit < 0:
+    acc.cyclotomic_inv(tab[digit shr 1])
+    return true
+  else:
+    acc.setOne()
+    return false
+
+func accumNAF[precompSize, NafMax: static int, Gt: ExtensionField](
+       acc: var Gt,
+       tab: array[precompSize, Gt],
+       naf: array[NafMax, int8], nafLen: int,
+       nafIteratorIdx: int) {.inline.} =
+
+    let digit = naf[nafLen-1-nafIteratorIdx]
+    if digit > 0:
+      acc *= tab[digit shr 1]
+    elif digit < 0:
+      var neg {.noInit.}: Gt
+      neg.cyclotomic_inv(tab[-digit shr 1])
+      acc *= neg
+
+func gtExp_minHammingWeight_windowed_vartime*[Gt: ExtensionField](
+        r: var Gt, a: Gt, scalar: BigInt, window: static int) {.tags:[VarTime], meter.} =
+  ## **Variable-time** Exponentiation in 𝔾ₜ
+  ##
+  ##   r <- aᵏ
+  ##
+  ## This uses windowed-NAF (wNAF)
+  ## This MUST NOT be used with secret data.
+  ##
+  ## This is highly VULNERABLE to timing attacks and power analysis attacks
+
+  # Signed digits divides precomputation table size by 2
+  # Odd-only divides precomputation table size by another 2
+  const precompSize = 1 shl (window - 2)
+  static: doAssert window <= 4, "Window of size " & $window & " is too large and precomputation would use " & $(precompSize * sizeof(Gt)) & " stack space."
+
+  var tab {.noinit.}: array[precompSize, Gt]
+  var a2{.noInit.}: Gt
+  tab[0] = a
+  a2.cyclotomic_square(a)
+  for i in 1 ..< tab.len:
+    tab[i].prod(tab[i-1], a2)
+
+  var naf {.noInit.}: array[BigInt.bits+1, int8]
+  let nafLen = naf.recode_r2l_signed_window_vartime(scalar, window)
+
+  var isInit = false
+  for i in 0 ..< nafLen:
+    if isInit:
+      r.cyclotomic_square()
+      r.accumNAF(tab, naf, nafLen, i)
+    else:
+      isInit = r.initNAF(tab, naf, nafLen, i)
+
+func gtExpEndo_minHammingWeight_windowed_vartime*[Gt: ExtensionField, scalBits: static int](
+        r: var Gt, a: Gt, scalar: BigInt[scalBits], window: static int) {.tags:[VarTime], meter.} =
+  ## Endomorphism accelerated **Variable-time** Exponentiation in 𝔾ₜ
+  ##
+  ##   r <- aᵏ
+  ##
+  ## This uses windowed-NAF (wNAF)
+  ## This MUST NOT be used with secret data.
+  ##
+  ## This is highly VULNERABLE to timing attacks and power analysis attacks
+
+  # Signed digits divides precomputation table size by 2
+  # Odd-only divides precomputation table size by another 2
+  const precompSize = 1 shl (window - 2)
+  static: doAssert window <= 4, "Window of size " & $window & " is too large and precomputation would use " & $(precompSize * sizeof(Gt)) & " stack space."
+
+  # 1. Compute endomorphisms
+  const M = when Gt is Fp6:  2
+            elif Gt is Fp12: 4
+            else: {.error: "Unconfigured".}
+
+  var endos {.noInit.}: array[M-1, Gt]
+  endos.computeEndomorphisms(a)
+
+  # 2. Decompose scalar into mini-scalars
+  const L = Fr[Gt.Name].bits().ceilDiv_vartime(M) + 1
+  var miniScalars {.noInit.}: array[M, BigInt[L]]
+  var negateElems {.noInit.}: array[M, SecretBool]
+  miniScalars.decomposeEndo(negateElems, scalar, Fr[Gt.Name].bits(), Gt.Name, G2) # 𝔾ₜ has same decomposition as 𝔾₂
+
+  # 3. Handle negative mini-scalars
+  if negateElems[0].bool:
+    r.cyclotomic_inv(a)
+  else:
+    r = a
+  for m in 1 ..< M:
+    if negateElems[m].bool:
+      endos[m-1].cyclotomic_inv()
+
+  # It's OK if r aliases a, we don't need a anymore
+
+  # 4. Precomputed table
+  var tab {.noinit.}: array[M, array[precompSize, Gt]]
+  for m in 0 ..< M:
+    var a2{.noInit.}: Gt
+    if m == 0:
+      tab[0][0] = r
+      a2.cyclotomic_square(r)
+    else:
+      tab[m][0] = endos[m-1]
+      a2.cyclotomic_square(endos[m-1])
+    for i in 1 ..< tab[m].len:
+      tab[m][i].prod(tab[m][i-1], a2)
+
+  # 5. wNAF precomputed tables
+  const NafLen = L+1
+  var tabNaf {.noinit.}: array[M, array[NafLen, int8]]
+
+  for m in 0 ..< M:
+    # tabNaf returns NAF from least-significant to most significant bits
+    let miniScalarLen = tabNaf[m].recode_r2l_signed_window_vartime(miniScalars[m], window)
+    # We compute from most significant to least significant
+    # so we pad with 0
+    for i in miniScalarLen ..< NafLen:
+      tabNaf[m][i] = 0
+
+  # 6. Compute
+  var isInit = false
+
+  for i in 0 ..< NafLen:
+    if isInit:
+      r.cyclotomic_square()
+    for m in 0 ..< M:
+      if isInit:
+        r.accumNAF(tab[m], tabNaf[m], NafLen, i)
+      else:
+        isInit = r.initNAF(tab[m], tabNaf[m], NafLen, i)
diff --git a/constantine/named/zoo_endomorphisms.nim b/constantine/named/zoo_endomorphisms.nim
index 0561b985..f49fd365 100644
--- a/constantine/named/zoo_endomorphisms.nim
+++ b/constantine/named/zoo_endomorphisms.nim
@@ -8,6 +8,7 @@
 
 import
   std/macros,
+  constantine/platforms/abstractions,
   constantine/math/extension_fields,
   constantine/math/isogenies/frobenius,
   constantine/math/elliptic/[
@@ -34,17 +35,13 @@ import
 macro dispatch(Name: static Algebra, tag: static string, G: static string): untyped =
   result = bindSym($Name & "_" & tag & "_" & G)
 
-template babai*(F: typedesc[Fp or Fp2]): untyped =
+template babai*(Name: static Algebra, G: static Subgroup): untyped =
   ## Return the GLV Babai roundings vector
-  const G = if F is Fp: "G1"
-            else: "G2"
-  dispatch(F.Name, "Babai", G)
+  dispatch(Name, "Babai", $G)
 
-template lattice*(F: typedesc[Fp or Fp2]): untyped =
+template lattice*(Name: static Algebra, G: static Subgroup): untyped =
   ## Returns the GLV Decomposition Lattice
-  const G = if F is Fp: "G1"
-            else: "G2"
-  dispatch(F.Name, "Lattice", G)
+  dispatch(Name, "Lattice", $G)
 
 macro getCubicRootOfUnity_mod_p*(Name: static Algebra): untyped =
   ## Get a non-trivial cubic root of unity (mod p) with p the prime field
@@ -90,7 +87,7 @@ func computeEndomorphism*[EC](endo: var EC, P: EC) =
   else: # For BW6-761, both G1 and G2 are on Fp
     endo.frobenius_psi(P, 2)
 
-func computeEndomorphisms*[EC; M: static int](endos: var array[M-1, EC], P: EC) =
+func computeEndomorphisms*[EC: not ExtensionField; M: static int](endos: var array[M-1, EC], P: EC) =
   ## An endomorphism decomposes M-way.
   when P.F is Fp:
     static: doAssert M == 2
@@ -103,6 +100,10 @@ func computeEndomorphisms*[EC; M: static int](endos: var array[M-1, EC], P: EC)
   else:
     {.error: "Unconfigured".}
 
+func computeEndomorphisms*[Gt: ExtensionField; M: static int](endos: var array[M-1, Gt], a: Gt) =
+  staticFor i, 0, M-1:
+    endos[i].frobenius_map(a, i+1)
+
 func hasEndomorphismAcceleration*(Name: static Algebra): bool {.compileTime.} =
   Name in {
     Bandersnatch,
diff --git a/sage/derive_endomorphisms.sage b/sage/derive_endomorphisms.sage
index 0f17413c..c17ecace 100644
--- a/sage/derive_endomorphisms.sage
+++ b/sage/derive_endomorphisms.sage
@@ -171,7 +171,7 @@ def genCubicRootEndo(curve_name, curve_config):
 
   return phi1, lattice, babai
 
-# G2 Endomorphism
+# 𝔾₂ Endomorphism
 # ---------------------------------------------------------
 
 def genPsiEndo(curve_name, curve_config):
@@ -258,7 +258,7 @@ if __name__ == "__main__":
     print('\nPrecomputing G1 - 𝜑 (phi) cubic root endomorphism')
     print('----------------------------------------------------\n')
     cubeRootModP, g1lat, g1babai = genCubicRootEndo(curve, Curves)
-    print('\n\nPrecomputing G2 - ψ (Psi) - untwist-Frobenius-twist endomorphism')
+    print('\n\nPrecomputing 𝔾₂ - ψ (Psi) - untwist-Frobenius-twist endomorphism')
     print('----------------------------------------------------\n')
     g2lat, g2babai = genPsiEndo(curve, Curves)
 
@@ -290,7 +290,7 @@ if __name__ == "__main__":
       ))
       f.write('\n\n')
       f.write(inspect.cleandoc(f"""
-        # {curve} G2
+        # {curve} 𝔾₂
         # ------------------------------------------------------------
       """))
       f.write('\n\n')
diff --git a/sage/derive_frobenius.sage b/sage/derive_frobenius.sage
index 12932644..8746437c 100644
--- a/sage/derive_frobenius.sage
+++ b/sage/derive_frobenius.sage
@@ -216,7 +216,7 @@ def genFrobeniusPsiConstants(curve_name, curve_config):
     # psi_2 ≡ ξ^((p-1)/6)^2 ≡ ξ^((p-1)/3)
     # psi_3 ≡ psi_2 * ξ^((p-1)/6) ≡ ξ^((p-1)/3) * ξ^((p-1)/6) ≡ ξ^((p-1)/2)
     #
-    # In Fp² (i.e. embedding degree of 12, G2 on Fp2)
+    # In Fp² (i.e. embedding degree of 12, 𝔾₂ on Fp2)
     # - quadratic non-residues respect the equation a^((p²-1)/2) ≡ -1 (mod p²) by the Legendre symbol
     # - sextic non-residues are also quadratic non-residues so ξ^((p²-1)/2) ≡ -1 (mod p²)
     # - QRT(1/a) = QRT(a) with QRT the quadratic residuosity test
@@ -229,7 +229,7 @@ def genFrobeniusPsiConstants(curve_name, curve_config):
     # So psi2_3 ≡ -1 (mod p²)
     #
     #
-    # In Fp (i.e. embedding degree of 6, G2 on Fp)
+    # In Fp (i.e. embedding degree of 6, 𝔾₂ on Fp)
     # - Fermat's Little Theorem gives us a^(p-1) ≡ 1 (mod p)
     #
     # psi2_3 ≡ ξ^((p-1)(p+1)/2) (mod p)
diff --git a/tests/math_pairings/t_pairing_bls12_381_gt_exp.nim b/tests/math_pairings/t_pairing_bls12_381_gt_exp.nim
new file mode 100644
index 00000000..2309c628
--- /dev/null
+++ b/tests/math_pairings/t_pairing_bls12_381_gt_exp.nim
@@ -0,0 +1,15 @@
+# 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
+  # Test utilities
+  ./t_pairing_template
+
+runGTexponentiationTests(
+  Iters = 4,
+  GT = Fp12[BLS12_381])
diff --git a/tests/math_pairings/t_pairing_bn254_snarks_gt_exp.nim b/tests/math_pairings/t_pairing_bn254_snarks_gt_exp.nim
new file mode 100644
index 00000000..dea4efd8
--- /dev/null
+++ b/tests/math_pairings/t_pairing_bn254_snarks_gt_exp.nim
@@ -0,0 +1,15 @@
+# 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
+  # Test utilities
+  ./t_pairing_template
+
+runGTexponentiationTests(
+  Iters = 4,
+  GT = Fp12[BN254_Nogami])
diff --git a/tests/math_pairings/t_pairing_template.nim b/tests/math_pairings/t_pairing_template.nim
index c7378d1f..fb3d3de7 100644
--- a/tests/math_pairings/t_pairing_template.nim
+++ b/tests/math_pairings/t_pairing_template.nim
@@ -16,7 +16,7 @@ import
   constantine/named/algebras,
   constantine/named/zoo_subgroups,
   constantine/math/elliptic/[ec_shortweierstrass_affine, ec_shortweierstrass_projective],
-  constantine/math/pairings/cyclotomic_subgroups,
+  constantine/math/pairings/[cyclotomic_subgroups, gt_exponentiations_vartime, pairings_generic],
   constantine/math/io/io_extfields,
 
   # Test utilities
@@ -27,7 +27,7 @@ export
   ec_shortweierstrass_affine, ec_shortweierstrass_projective,
   arithmetic, extension_fields,
   io_extfields,
-  cyclotomic_subgroups,
+  cyclotomic_subgroups, gt_exponentiations_vartime,
   abstractions, algebras
 
 type
@@ -105,7 +105,7 @@ template runPairingTests*(Iters: static int, Name: static Algebra, G1, G2, GT: t
       test_bilinearity_double_impl(randZ = false, gen = HighHammingWeight)
       test_bilinearity_double_impl(randZ = false, gen = Long01Sequence)
 
-func random_elem*(rng: var RngState, F: typedesc, gen: RandomGen): F {.inline, noInit.} =
+func random_elem(rng: var RngState, F: typedesc, gen: RandomGen): F {.inline, noInit.} =
   if gen == Uniform:
     result = rng.random_unsafe(F)
   elif gen == HighHammingWeight:
@@ -114,7 +114,7 @@ func random_elem*(rng: var RngState, F: typedesc, gen: RandomGen): F {.inline, n
     result = rng.random_long01Seq(F)
 
 template runGTsubgroupTests*(Iters: static int, GT: typedesc, finalExpHard_fn: untyped): untyped {.dirty.}=
-  bind affineType
+  bind affineType, random_elem
 
   var rng: RngState
   let timeseed = uint32(toUnix(getTime()) and (1'i64 shl 32 - 1)) # unixTime mod 2^32
@@ -139,8 +139,76 @@ template runGTsubgroupTests*(Iters: static int, GT: typedesc, finalExpHard_fn: u
 
     stdout.write '\n'
 
-  suite "Pairing - GT subgroup " & $GT.Name & " [" & $WordBitWidth & "-bit words]":
-    test "Final Exponentiation and GT-subgroup membership":
+  suite "Pairing - 𝔾ₜ subgroup " & $GT.Name & " [" & $WordBitWidth & "-bit words]":
+    test "Final Exponentiation and 𝔾ₜ-subgroup membership":
       test_gt_impl(gen = Uniform)
       test_gt_impl(gen = HighHammingWeight)
       test_gt_impl(gen = Long01Sequence)
+
+func random_gt*(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)
+
+  result.finalExp()
+
+template runGTexponentiationTests*(Iters: static int, GT: typedesc): untyped {.dirty.} =
+  var rng: RngState
+  let timeseed = uint32(toUnix(getTime()) and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+  seed(rng, timeseed)
+  echo "\n------------------------------------------------------\n"
+  echo "test_pairing_",$GT.Name,"_gt_exponentiation xoshiro512** seed: ", timeseed
+
+  proc test_gt_exponentiation_impl(gen: RandomGen) =
+    stdout.write "    "
+    for _ in 0 ..< Iters:
+      let a = rng.random_gt(GT, gen)
+      let kUnred = rng.random_long01seq(GT.Name.getBigInt(kScalarField))
+      var k {.noInit.}: GT.Name.getBigInt(kScalarField)
+      discard k.reduce_vartime(kUnred, GT.Name.scalarFieldModulus())
+
+      # Reference impl using exponentiation with tables on any field/extension field
+      var r_ref = a
+      r_ref.pow_vartime(k, window = 3)
+
+      # Square-and-multiply
+      var r_sqrmul {.noInit.}: GT
+      r_sqrmul.gtExp_sqrmul_vartime(a, k)
+      doAssert bool(r_ref == r_sqrmul)
+
+      # MSB->LSB min Hamming Weight signed recoding
+      var r_l2r_recoding {.noInit.}: GT
+      r_l2r_recoding.gtExp_minHammingWeight_vartime(a, k)
+      doAssert bool(r_ref == r_l2r_recoding)
+
+      # Windowed NAF
+      var r_wNAF {.noInit.}: GT
+      r_wNAF.gtExp_minHammingWeight_windowed_vartime(a, k, window = 2)
+      doAssert bool(r_ref == r_wNAF)
+      r_wNAF.gtExp_minHammingWeight_windowed_vartime(a, k, window = 3)
+      doAssert bool(r_ref == r_wNAF)
+      r_wNAF.gtExp_minHammingWeight_windowed_vartime(a, k, window = 4)
+      doAssert bool(r_ref == r_wNAF)
+
+      # Windowed NAF + endomorphism acceleration
+      var r_endoWNAF {.noInit.}: GT
+      r_endoWNAF.gtExpEndo_minHammingWeight_windowed_vartime(a, k, window = 2)
+      doAssert bool(r_ref == r_endoWNAF)
+      r_endoWNAF.gtExpEndo_minHammingWeight_windowed_vartime(a, k, window = 3)
+      doAssert bool(r_ref == r_endoWNAF)
+      r_endoWNAF.gtExpEndo_minHammingWeight_windowed_vartime(a, k, window = 4)
+      doAssert bool(r_ref == r_endoWNAF)
+
+      stdout.write '.'
+
+    stdout.write '\n'
+
+
+  suite "Pairing - Exponentiation for 𝔾ₜ " & $GT.Name & " [" & $WordBitWidth & "-bit words]":
+    test "𝔾ₜ exponentiation consistency":
+      test_gt_exponentiation_impl(gen = Uniform)
+      test_gt_exponentiation_impl(gen = HighHammingWeight)
+      test_gt_exponentiation_impl(gen = Long01Sequence)