Naive pairings + Naive cofactor clearing (#82)

* Pairing - initial commit - line functions - sparse Fp12 functions * Small fixes: - Line parametrized by twist for generic algorithm - Add a conjugate operator for quadratic extensions - Have frobenius use it - Create an Affine coordinate type for elliptic curve * Implement (failing) pairing test * Stash pairing debug session, temp switch Fp12 over Fp4 * Proper naive pairing on BLS12-381 * Frobenius map * Implement naive pairing for BN curves * Add pairing tests to CI + reduce time spent on lower-level tests * Test without assembler in Github Actions + less base layers test iterations
2020-09-21 23:24:00 +02:00 · 2020-09-21 23:24:00 +02:00 · d84edcd217
parent 28e83e7b49
commit d84edcd217
67 changed files with 2346 additions and 255 deletions
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@ -55,6 +55,11 @@ jobs:
    name: '${{ matrix.target.os }}-${{ matrix.target.cpu }}-${{ matrix.target.TEST_LANG }} (${{ matrix.branch }})'
    runs-on: ${{ matrix.builder }}
    steps:
+      - name: Cancel Previous Runs
+        uses: styfle/cancel-workflow-action@0.5.0
+        with:
+          access_token: ${{ github.token }}
+
      - name: Checkout constantine
        uses: actions/checkout@v2
        with:
@ -197,13 +202,20 @@ jobs:
        run: |
          nimble refresh
          nimble install -y gmp stew
-      - name: Run Constantine tests (with GMP)
+      - name: Run Constantine tests (with Assembler & with GMP)
        if: (runner.os == 'Linux' || runner.os == 'macOS') && matrix.target.cpu == 'amd64'
        shell: bash
        run: |
          export UCPU="$cpu"
          cd constantine
          nimble test_parallel
+      - name: Run Constantine tests (no Assembler & with GMP)
+        if: (runner.os == 'Linux' || runner.os == 'macOS') && matrix.target.cpu == 'amd64'
+        shell: bash
+        run: |
+          export UCPU="$cpu"
+          cd constantine
+          nimble test_parallel_no_assembler
      - name: Run Constantine tests (without GMP)
        if: runner.os == 'Linux' && matrix.target.cpu == 'i386'
        shell: bash
@ -211,3 +223,10 @@ jobs:
          export UCPU="$cpu"
          cd constantine
          nimble test_parallel_no_gmp
+      - name: Run Constantine tests (without Assembler or GMP)
+        if: runner.os == 'Linux' && matrix.target.cpu == 'i386'
+        shell: bash
+        run: |
+          export UCPU="$cpu"
+          cd constantine
+          nimble test_parallel_no_gmp_no_assembler
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@ -240,5 +240,5 @@ steps:
      echo "PATH=${PATH}"
      export ucpu=${UCPU}
      nimble test_no_gmp
-    displayName: 'Testing the package (without Assembler or GMP)'
+    displayName: 'Testing the package (without GMP and ASM on Windows)'
    condition: eq(variables['Agent.OS'], 'Windows_NT')
--- a/benchmarks/bench_ec_g1.nim
+++ b/benchmarks/bench_ec_g1.nim
@ -26,8 +26,8 @@ import
 # ############################################################


-const Iters = 1_000_000
-const MulIters = 1000
+const Iters = 10_000
+const MulIters = 100
 const AvailableCurves = [
  # P224,
  # BN254_Nogami,
--- a/benchmarks/bench_ec_g2.nim
+++ b/benchmarks/bench_ec_g2.nim
@ -27,7 +27,7 @@ import
 # ############################################################


-const Iters = 500_000
+const Iters = 10_000
 const MulIters = 500
 const AvailableCurves = [
  # P224,
--- a/benchmarks/bench_fp6.nim
+++ b/benchmarks/bench_fp6.nim
@ -23,7 +23,7 @@ import
 # ############################################################


-const Iters = 1_000_000
+const Iters = 100_000
 const InvIters = 1000
 const AvailableCurves = [
  # Pairing-Friendly curves
--- a/constantine.nimble
+++ b/constantine.nimble
@ -46,6 +46,9 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
  ("tests/t_fp12_bn254_snarks.nim", false),
  ("tests/t_fp12_bls12_377.nim", false),
  ("tests/t_fp12_bls12_381.nim", false),
+  ("tests/t_fp12_exponentiation.nim", false),
+
+  ("tests/t_fp4_frobenius.nim", false),
  # Elliptic curve arithmetic G1
  ("tests/t_ec_wstrass_prj_g1_add_double.nim", false),
  ("tests/t_ec_wstrass_prj_g1_mul_sanity.nim", false),
@ -66,7 +69,12 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
  ("tests/t_ec_sage_bn254.nim", false),
  ("tests/t_ec_sage_bls12_381.nim", false),
  # Edge cases highlighted by past bugs
-  ("tests/t_ec_wstrass_prj_edge_cases.nim", false)
+  ("tests/t_ec_wstrass_prj_edge_cases.nim", false),
+  # Pairing
+  ("tests/t_pairing_fp12_sparse.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)
 ]

 # For temporary (hopefully) investigation that can only be reproduced in CI
@ -139,8 +147,6 @@ task test, "Run all tests":
      else:
        test "-d:Constantine32", td.path

-  # Benchmarks compile and run
-  # ignore Windows 32-bit for the moment
  # Ensure benchmarks stay relevant. Ignore Windows 32-bit at the moment
  if not defined(windows) or not (existsEnv"UCPU" or getEnv"UCPU" == "i686"):
    runBench("bench_fp")
@ -167,9 +173,6 @@ task test_no_gmp, "Run tests that don't require GMP":
        else:
          test "-d:Constantine32", td.path

-
-  # Benchmarks compile and run
-  # ignore Windows 32-bit for the moment
  # Ensure benchmarks stay relevant. Ignore Windows 32-bit at the moment
  if not defined(windows) or not (existsEnv"UCPU" or getEnv"UCPU" == "i686"):
    runBench("bench_fp")
@ -298,6 +301,47 @@ task test_parallel_no_gmp, "Run all tests in parallel (via GNU parallel)":
    runBench("bench_ec_g1")
    runBench("bench_ec_g2")

+task test_parallel_no_gmp_no_assembler, "Run all tests in parallel (via GNU parallel)":
+  # -d:testingCurves is configured in a *.nim.cfg for convenience
+  let cmdFile = true # open(buildParallel, mode = fmWrite) # Nimscript doesn't support IO :/
+  exec "> " & buildParallel
+
+  for td in testDesc:
+    if not td.useGMP:
+      if td.path in useDebug:
+        test "-d:debugConstantine -d:ConstantineASM=false", td.path, cmdFile
+      else:
+        test "-d:ConstantineASM=false", td.path, cmdFile
+
+  # cmdFile.close()
+  # Execute everything in parallel with GNU parallel
+  exec "parallel --keep-order --group < " & buildParallel
+
+  exec "> " & buildParallel
+  if sizeof(int) == 8: # 32-bit tests on 64-bit arch
+    for td in testDesc:
+      if not td.useGMP:
+        if td.path in useDebug:
+          test "-d:Constantine32 -d:debugConstantine", td.path, cmdFile
+        else:
+          test "-d:Constantine32", td.path, cmdFile
+    # cmdFile.close()
+    # Execute everything in parallel with GNU parallel
+    exec "parallel --keep-order --group < " & buildParallel
+
+  # Now run the benchmarks
+  #
+  # Benchmarks compile and run
+  # ignore Windows 32-bit for the moment
+  # Ensure benchmarks stay relevant. Ignore Windows 32-bit at the moment
+  if not defined(windows) or not (existsEnv"UCPU" or getEnv"UCPU" == "i686"):
+    runBench("bench_fp")
+    runBench("bench_fp2")
+    runBench("bench_fp6")
+    runBench("bench_fp12")
+    runBench("bench_ec_g1")
+    runBench("bench_ec_g2")
+
 task bench_fp, "Run benchmark 𝔽p with your default compiler":
  runBench("bench_fp")

--- a/constantine/arithmetic/bigints.nim
+++ b/constantine/arithmetic/bigints.nim
@ -318,6 +318,87 @@ func bit0*(a: BigInt): Ct[uint8] =
  ## Access the least significant bit
  ct(a.limbs[0] and SecretWord(1), uint8)

+# Multiplication by small cosntants
+# ------------------------------------------------------------
+
+func `*=`*(a: var BigInt, b: static int) {.inline.} =
+  ## Multiplication by a small integer known at compile-time
+  # Implementation:
+  #
+  # we hardcode addition chains for small integer
+  const negate = b < 0
+  const b = if negate: -b
+            else: b
+  when negate:
+    a.neg(a)
+  when b == 0:
+    a.setZero()
+  elif b == 1:
+    return
+  elif b == 2:
+    discard a.double()
+  elif b == 3:
+    let t1 = a
+    discard a.double()
+    a += t1
+  elif b == 4:
+    discard a.double()
+    discard a.double()
+  elif b == 5:
+    let t1 = a
+    discard a.double()
+    discard a.double()
+    a += t1
+  elif b == 6:
+    discard a.double()
+    let t2 = a
+    discard a.double() # 4
+    a += t2
+  elif b == 7:
+    let t1 = a
+    discard a.double()
+    let t2 = a
+    discard a.double() # 4
+    a += t2
+    a += t1
+  elif b == 8:
+    discard a.double()
+    discard a.double()
+    discard a.double()
+  elif b == 9:
+    let t1 = a
+    discard a.double()
+    discard a.double()
+    discard a.double() # 8
+    a += t1
+  elif b == 10:
+    discard a.double()
+    let t2 = a
+    discard a.double()
+    discard a.double() # 8
+    a += t2
+  elif b == 11:
+    let t1 = a
+    discard a.double()
+    let t2 = a
+    discard a.double()
+    discard a.double() # 8
+    a += t2
+    a += t1
+  elif b == 12:
+    discard a.double()
+    discard a.double() # 4
+    let t4 = a
+    discard a.double() # 8
+    a += t4
+  else:
+    {.error: "Multiplication by this small int not implemented".}
+
+func `*`*(b: static int, a: BigInt): BigInt {.noinit, inline.} =
+  ## Multiplication by a small integer known at compile-time
+  result = a
+  result *= b
+
 # ############################################################
 #
 #                   Modular BigInt
@ -417,58 +498,6 @@ func montySquare*(r: var BigInt, a, M: BigInt, negInvModWord: static BaseType, c
  ## to avoid duplicating with Nim zero-init policy
  montySquare(r.limbs, a.limbs, M.limbs, negInvModWord, canUseNoCarryMontyMul)

-func montyPow*[mBits, eBits: static int](
-       a: var BigInt[mBits], exponent: BigInt[eBits],
-       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
-       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
-      ) =
-  ## Compute a <- a^exponent (mod M)
-  ## ``a`` in the Montgomery domain
-  ## ``exponent`` is any BigInt, in the canonical domain
-  ##
-  ## This uses fixed window optimization
-  ## A window size in the range [1, 5] must be chosen
-  ##
-  ## This is constant-time: the window optimization does
-  ## not reveal the exponent bits or hamming weight
-  mixin exportRawUint # exported in io_bigints which depends on this module ...
-
-  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
-  expBE.exportRawUint(exponent, bigEndian)
-
-  const scratchLen = if windowSize == 1: 2
-                     else: (1 shl windowSize) + 1
-  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
-  montyPow(a.limbs, expBE, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
-
-func montyPowUnsafeExponent*[mBits, eBits: static int](
-       a: var BigInt[mBits], exponent: BigInt[eBits],
-       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
-       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
-      ) =
-  ## Compute a <- a^exponent (mod M)
-  ## ``a`` in the Montgomery domain
-  ## ``exponent`` is any BigInt, in the canonical domain
-  ##
-  ## Warning ⚠️ :
-  ## This is an optimization for public exponent
-  ## Otherwise bits of the exponent can be retrieved with:
-  ## - memory access analysis
-  ## - power analysis
-  ## - timing analysis
-  ##
-  ## This uses fixed window optimization
-  ## A window size in the range [1, 5] must be chosen
-  mixin exportRawUint # exported in io_bigints which depends on this module ...
-
-  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
-  expBE.exportRawUint(exponent, bigEndian)
-
-  const scratchLen = if windowSize == 1: 2
-                     else: (1 shl windowSize) + 1
-  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
-  montyPowUnsafeExponent(a.limbs, expBE, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
-
 func montyPow*[mBits: static int](
       a: var BigInt[mBits], exponent: openarray[byte],
       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
@ -517,5 +546,51 @@ func montyPowUnsafeExponent*[mBits: static int](
  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
  montyPowUnsafeExponent(a.limbs, exponent, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)

+func montyPow*[mBits, eBits: static int](
+       a: var BigInt[mBits], exponent: BigInt[eBits],
+       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
+       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
+      ) =
+  ## Compute a <- a^exponent (mod M)
+  ## ``a`` in the Montgomery domain
+  ## ``exponent`` is any BigInt, in the canonical domain
+  ##
+  ## This uses fixed window optimization
+  ## A window size in the range [1, 5] must be chosen
+  ##
+  ## This is constant-time: the window optimization does
+  ## not reveal the exponent bits or hamming weight
+  mixin exportRawUint # exported in io_bigints which depends on this module ...
+
+  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
+  expBE.exportRawUint(exponent, bigEndian)
+
+  montyPow(a, expBE, M, one, negInvModWord, windowSize, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
+
+func montyPowUnsafeExponent*[mBits, eBits: static int](
+       a: var BigInt[mBits], exponent: BigInt[eBits],
+       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
+       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
+      ) =
+  ## Compute a <- a^exponent (mod M)
+  ## ``a`` in the Montgomery domain
+  ## ``exponent`` is any BigInt, in the canonical domain
+  ##
+  ## Warning ⚠️ :
+  ## This is an optimization for public exponent
+  ## Otherwise bits of the exponent can be retrieved with:
+  ## - memory access analysis
+  ## - power analysis
+  ## - timing analysis
+  ##
+  ## This uses fixed window optimization
+  ## A window size in the range [1, 5] must be chosen
+  mixin exportRawUint # exported in io_bigints which depends on this module ...
+
+  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
+  expBE.exportRawUint(exponent, bigEndian)
+
+  montyPowUnsafeExponent(a, expBE, M, one, negInvModWord, windowSize, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
+
 {.pop.} # inline
 {.pop.} # raises no exceptions
--- a/constantine/arithmetic/finite_fields.nim
+++ b/constantine/arithmetic/finite_fields.nim
@ -225,6 +225,10 @@ func neg*(r: var Fp, a: Fp) {.inline.} =
  else:
    discard r.mres.diff(Fp.C.Mod, a.mres)

+func neg*(a: var Fp) {.inline.} =
+  ## Negate modulo p
+  a.neg(a)
+
 func div2*(a: var Fp) {.inline.} =
  ## Modular division by 2
  a.mres.div2_modular(Fp.C.getPrimePlus1div2())
--- a/constantine/config/curves.nim
+++ b/constantine/config/curves.nim
@ -94,6 +94,7 @@ macro get_CNR_Fp2*(C: static Curve): untyped =
  result = bindSym($C & "_nonresidue_cube_fp2")

 macro getSexticTwist*(C: static Curve): untyped =
+  ## Returns if D-Twist or M-Twist
  result = bindSym($C & "_sexticTwist")

 macro get_SNR_Fp2*(C: static Curve): untyped =
--- a/constantine/config/curves_declaration.nim
+++ b/constantine/config/curves_declaration.nim
@ -85,6 +85,19 @@ declareCurves:
    family: BarretoNaehrig
    # Equation: Y^2 = X^3 + 2
    # u: -(2^62 + 2^55 + 1)
+
+    order: "0x2523648240000001ba344d8000000007ff9f800000000010a10000000000000d"
+    orderBitwidth: 254
+    cofactor: 1
+    eq_form: ShortWeierstrass
+    coef_a: 0
+    coef_b: 2
+    nonresidue_quad_fp: -1       #      -1   is not a square in 𝔽p
+    nonresidue_cube_fp2: (1, 1)  # 1+𝑖   1+𝑖  is not a cube in 𝔽p²
+
+    sexticTwist: D_Twist
+    sexticNonResidue_fp2: (1, 1) # 1+𝑖
+
  curve BN254_Snarks: # Zero-Knowledge proofs curve (SNARKS, STARKS, Ethereum)
    bitwidth: 254
    modulus: "0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47"
--- a/constantine/elliptic/ec_weierstrass_affine.nim
+++ b/constantine/elliptic/ec_weierstrass_affine.nim
@ -13,6 +13,20 @@ import
  ../towers,
  ../io/io_bigints

+# ############################################################
+#
+#             Elliptic Curve in Short Weierstrass form
+#                 with Projective Coordinates
+#
+# ############################################################
+
+type ECP_SWei_Aff*[F] = object
+  ## Elliptic curve point for a curve in Short Weierstrass form
+  ##   y² = x³ + a x + b
+  ##
+  ## over a field F
+  x*, y*: F
+
 func curve_eq_rhs*[F](y2: var F, x: F) =
  ## Compute the curve equation right-hand-side from field element `x`
  ## i.e.  `y²` in `y² = x³ + a x + b`
@ -65,3 +79,35 @@ func isOnCurve*[F](x, y: F): SecretBool =
  rhs.curve_eq_rhs(x)

  return y2 == rhs
+
+func trySetFromCoordX*[F](P: var ECP_SWei_Aff[F], x: F): SecretBool =
+  ## Try to create a point the elliptic curve
+  ## y² = x³ + a x + b     (affine coordinate)
+  ##
+  ## The `Z` coordinates is set to 1
+  ##
+  ## return true and update `P` if `x` leads to a valid point
+  ## return false otherwise, in that case `P` is undefined.
+  ##
+  ## Note: Dedicated robust procedures for hashing-to-curve
+  ##       will be provided, this is intended for testing purposes.
+  P.y.curve_eq_rhs(x)
+  # TODO: supports non p ≡ 3 (mod 4) modulus like BLS12-377
+  result = sqrt_if_square(P.y)
+
+func neg*(P: var ECP_SWei_Aff, Q: ECP_SWei_Aff) =
+  ## Negate ``P``
+  P.x = Q.x
+  P.y.neg(Q.y)
+
+func neg*(P: var ECP_SWei_Aff) =
+  ## Negate ``P``
+  P.y.neg()
+
+func cneg*(P: var ECP_SWei_Aff, ctl: CTBool) =
+  ## Conditional negation.
+  ## Negate if ``ctl`` is true
+  var Q{.noInit.}: typeof(P)
+  Q.x = P.x
+  Q.y.neg(P.y)
+  P.ccopy(Q, ctl)
--- a/constantine/elliptic/ec_weierstrass_projective.nim
+++ b/constantine/elliptic/ec_weierstrass_projective.nim
@ -15,7 +15,7 @@ import

 # ############################################################
 #
-#             Elliptic Curve in Weierstrass form
+#             Elliptic Curve in Short Weierstrass form
 #                 with Projective Coordinates
 #
 # ############################################################
@ -104,9 +104,15 @@ func trySetFromCoordX*[F](P: var ECP_SWei_Proj[F], x: F): SecretBool =
  P.x = x
  P.z.setOne()

+func neg*(P: var ECP_SWei_Proj, Q: ECP_SWei_Proj) =
+  ## Negate ``P``
+  P.x = Q.x
+  P.y.neg(Q.y)
+  P.z = Q.z
+
 func neg*(P: var ECP_SWei_Proj) =
  ## Negate ``P``
-  P.y.neg(P.y)
+  P.y.neg()

 func cneg*(P: var ECP_SWei_Proj, ctl: CTBool) =
  ## Conditional negation.
@ -309,13 +315,17 @@ func diff*[F](r: var ECP_SWei_Proj[F],
  nQ.neg()
  r.sum(P, nQ)

-func affineFromProjective*[F](aff: var ECP_SWei_Proj[F], proj: ECP_SWei_Proj) =
+func affineFromProjective*[F](aff: var ECP_SWei_Aff[F], proj: ECP_SWei_Proj) =
  # TODO: for now we reuse projective coordinate backend
  # TODO: simultaneous inversion

  var invZ {.noInit.}: F
  invZ.inv(proj.z)

-  aff.z.setOne()
  aff.x.prod(proj.x, invZ)
  aff.y.prod(proj.y, invZ)
+
+func projectiveFromAffine*[F](proj: var ECP_SWei_Proj, aff: ECP_SWei_Aff[F]) {.inline.} =
+  proj.x = aff.x
+  proj.y = aff.y
+  proj.z.setOne()
--- a/constantine/hash_to_curve/cofactors.nim
+++ b/constantine/hash_to_curve/cofactors.nim
@ -0,0 +1,69 @@
+# 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
+  ../config/common,
+  ../arithmetic,
+  ../primitives,
+  ../towers,
+  ../config/curves,
+  ../io/io_bigints,
+  ../elliptic/[ec_weierstrass_projective, ec_scalar_mul]
+
+# ############################################################
+#
+#                Clear Cofactor - Naive
+#
+# ############################################################
+
+const Cofactor_Eff_BN254_Nogami_G1 = BigInt[1].fromHex"0x1"
+const Cofactor_Eff_BN254_Nogami_G2 = BigInt[254].fromHex"0x2523648240000001ba344d8000000008c2a2800000000016ad00000000000019"
+  ## G2.order // r
+
+const Cofactor_Eff_BN254_Snarks_G1 = BigInt[1].fromHex"0x1"
+const Cofactor_Eff_BN254_Snarks_G2 = BigInt[254].fromHex"0x30644e72e131a029b85045b68181585e06ceecda572a2489345f2299c0f9fa8d"
+  ## G2.order // r
+
+# https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-09#section-8.8
+const Cofactor_Eff_BLS12_381_G1 = BigInt[64].fromHex"0xd201000000010001"
+  ## P -> (1 - x) P
+const Cofactor_Eff_BLS12_381_G2 = BigInt[636].fromHex"0xbc69f08f2ee75b3584c6a0ea91b352888e2a8e9145ad7689986ff031508ffe1329c2f178731db956d82bf015d1212b02ec0ec69d7477c1ae954cbc06689f6a359894c0adebbf6b4e8020005aaa95551"
+  ## P -> (x^2 - x - 1) P + (x - 1) psi(P) + psi(psi(2P))
+
+func clearCofactorReference*(P: var ECP_SWei_Proj[Fp[BN254_Nogami]]) {.inline.} =
+  ## Clear the cofactor of BN254_Nogami G1
+  ## BN curve have a G1 cofactor of 1 so this is a no-op
+  discard
+
+func clearCofactorReference*(P: var ECP_SWei_Proj[Fp2[BN254_Nogami]]) {.inline.} =
+  ## Clear the cofactor of BN254_Snarks G2
+  # Endomorphism acceleration cannot be used if cofactor is not cleared
+  P.scalarMulGeneric(Cofactor_Eff_BN254_Nogami_G2)
+
+func clearCofactorReference*(P: var ECP_SWei_Proj[Fp[BN254_Snarks]]) {.inline.} =
+  ## Clear the cofactor of BN254_Snarks G1
+  ## BN curve have a G1 cofactor of 1 so this is a no-op
+  discard
+
+func clearCofactorReference*(P: var ECP_SWei_Proj[Fp2[BN254_Snarks]]) {.inline.} =
+  ## Clear the cofactor of BN254_Snarks G2
+  # Endomorphism acceleration cannot be used if cofactor is not cleared
+  P.scalarMulGeneric(Cofactor_Eff_BN254_Snarks_G2)
+
+func clearCofactorReference*(P: var ECP_SWei_Proj[Fp[BLS12_381]]) {.inline.} =
+  ## Clear the cofactor of BLS12_381 G1
+  ## BN curve have a G1 cofactor of 1 so this is a no-op
+  P.scalarMulGeneric(Cofactor_Eff_BLS12_381_G1)
+
+func clearCofactorReference*(P: var ECP_SWei_Proj[Fp2[BLS12_381]]) {.inline.} =
+  ## Clear the cofactor of BLS12_381 G2
+  # Endomorphism acceleration cannot be used if cofactor is not cleared
+  P.scalarMulGeneric(Cofactor_Eff_BLS12_381_G2)
--- a/constantine/io/io_ec.nim
+++ b/constantine/io/io_ec.nim
@ -40,7 +40,7 @@ func toHex*(P: ECP_SWei_Proj): string =
  ##
  ## This proc output may change format in the future

-  var aff {.noInit.}: typeof(P)
+  var aff {.noInit.}: ECP_SWei_Aff[ECP_SWei_Proj.F]
  aff.affineFromProjective(P)

  result = "ECP[" & $aff.F & "](\n  x: "
@ -70,3 +70,23 @@ func fromHex*(dst: var ECP_SWei_Proj, x0, x1, y0, y1: string): bool {.raises: [V
  dst.y.fromHex(y0, y1)
  dst.z.setOne()
  return bool(isOnCurve(dst.x, dst.y))
+
+func fromHex*(dst: var ECP_SWei_Aff, x, y: string): bool {.raises: [ValueError].}=
+  ## Convert hex strings to a G1 curve point
+  ## Returns `false`
+  ## if there is no point with coordinates (`x`, `y`) on the curve
+  ## In that case, `dst` content is undefined.
+  static: doAssert dst.F is Fp, "dst must be on G1, an elliptic curve over 𝔽p"
+  dst.x.fromHex(x)
+  dst.y.fromHex(y)
+  return bool(isOnCurve(dst.x, dst.y))
+
+func fromHex*(dst: var ECP_SWei_Aff, x0, x1, y0, y1: string): bool {.raises: [ValueError].}=
+  ## Convert hex strings to a G2 curve point
+  ## Returns `false`
+  ## if there is no point with coordinates (`x`, `y`) on the curve
+  ## In that case, `dst` content is undefined.
+  static: doAssert dst.F is Fp2, "dst must be on G2, an elliptic curve over 𝔽p2"
+  dst.x.fromHex(x0, x1)
+  dst.y.fromHex(y0, y1)
+  return bool(isOnCurve(dst.x, dst.y))
--- a/constantine/io/io_towers.nim
+++ b/constantine/io/io_towers.nim
@ -26,7 +26,7 @@ import
 #
 # ############################################################

-func appendHex*(accum: var string, f: Fp2 or Fp6 or Fp12, order: static Endianness = bigEndian) =
+func appendHex*(accum: var string, f: Fp2 or Fp4 or Fp6 or Fp12, order: static Endianness = bigEndian) =
  ## Hex accumulator
  accum.add static($f.typeof.genericHead() & '(')
  for fieldName, fieldValue in fieldPairs(f):
@ -36,7 +36,7 @@ func appendHex*(accum: var string, f: Fp2 or Fp6 or Fp12, order: static Endianne
    accum.appendHex(fieldValue, order)
  accum.add ")"

-func toHex*(f: Fp2 or Fp6 or Fp12, order: static Endianness = bigEndian): string =
+func toHex*(f: Fp2 or Fp4 or Fp6 or Fp12, order: static Endianness = bigEndian): string =
  ## Stringify a tower field element to hex.
  ## Note. Leading zeros are not removed.
  ## Result is prefixed with 0x
@ -59,3 +59,49 @@ func fromHex*(T: typedesc[Fp2], c0, c1: string): T {.raises: [ValueError].}=
  ## with dst = c0 + β * c1
  ## β is the quadratic non-residue chosen to construct 𝔽p2
  result.fromHex(c0, c1)
+
+func fromHex*(dst: var Fp4,
+              c0, c1, c2, c3: string) {.raises: [ValueError].}=
+  ## Convert 4 coordinates to an element of 𝔽p4
+  dst.c0.fromHex(c0, c1)
+  dst.c1.fromHex(c2, c3)
+
+func fromHex*(T: typedesc[Fp4],
+              c0, c1, c2: string,
+              c3, c4, c5: string): T {.raises: [ValueError].}=
+  ## Convert 4 coordinates to an element of 𝔽p4
+  result.fromHex(c0, c1, c2, c3)
+
+func fromHex*(dst: var Fp6,
+              c0, c1, c2: string,
+              c3, c4, c5: string) {.raises: [ValueError].}=
+  ## Convert 6 coordinates to an element of 𝔽p6
+  dst.c0.fromHex(c0, c1)
+  dst.c1.fromHex(c2, c3)
+  dst.c2.fromHex(c4, c5)
+
+func fromHex*(T: typedesc[Fp6],
+              c0, c1, c2: string,
+              c3, c4, c5: string): T {.raises: [ValueError].}=
+  ## Convert 6 coordinates to an element of 𝔽p6
+  result.fromHex(c0, c1, c2, c3, c4, c5)
+
+func fromHex*(dst: var Fp12,
+              c0, c1, c2, c3: string,
+              c4, c5, c6, c7: string,
+              c8, c9, c10, c11: string) {.raises: [ValueError].}=
+  ## Convert 12 coordinates to an element of 𝔽p12
+  when Fp12.c0 is Fp6:
+    dst.c0.fromHex(c0, c1, c2, c3, c4, c5)
+    dst.c1.fromHex(c6, c7, c8, c9, c10, c11)
+  else:
+    dst.c0.fromHex(c0, c1, c2, c3)
+    dst.c1.fromHex(c4, c5, c6, c7)
+    dst.c2.fromHex(c8, c9, c10, c11)
+
+func fromHex*(T: typedesc[Fp12],
+              c0, c1, c2, c3: string,
+              c4, c5, c6, c7: string,
+              c8, c9, c10, c11: string): T {.raises: [ValueError].}=
+  ## Convert 12 coordinates to an element of 𝔽p12
+  result.fromHex(c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11)
--- a/constantine/isogeny/frobenius.nim
+++ b/constantine/isogeny/frobenius.nim
@ -12,7 +12,7 @@ import
  ../io/io_towers,
  ../towers, ../arithmetic

-# Frobenius automorphism
+# Frobenius Map
 # ------------------------------------------------------------
 #
 # https://en.wikipedia.org/wiki/Frobenius_endomorphism
@ -35,91 +35,13 @@ import
 # whether u = √-1 = i
 # or          √-2 or √-5

-func frobenius*(r: var Fp2, a: Fp2, k: static int = 1) {.inline.} =
+func frobenius_map*(r: var Fp2, a: Fp2, k: static int = 1) {.inline.} =
  ## Computes a^(p^k)
  ## The p-power frobenius automorphism on 𝔽p2
-  r.c0 = a.c0
  when (k and 1) == 1:
-    r.c1.neg(a.c1)
+    r.conj(a)
  else:
-    r.c1 = a.c1
-
-# Frobenius endomorphism
-# ------------------------------------------------------------
-# TODO: generate those constants via Sage in a Json file
-#       and parse at compile-time
-
-# Constants:
-#   Assuming embedding degree of 12 and a sextic twist
-#   with SNR the sextic non-residue
-#
-#   BN254_Snarks is a D-Twist: SNR^((p-1)/6)
-const FrobConst_BN254_Snarks_psi1_coef1 = Fp2[BN254_Snarks].fromHex(
-  "0x1284b71c2865a7dfe8b99fdd76e68b605c521e08292f2176d60b35dadcc9e470",
-  "0x246996f3b4fae7e6a6327cfe12150b8e747992778eeec7e5ca5cf05f80f362ac"
-)
-#  SNR^((p-1)/3)
-const FrobConst_BN254_Snarks_psi1_coef2 = Fp2[BN254_Snarks].fromHex(
-  "0x2fb347984f7911f74c0bec3cf559b143b78cc310c2c3330c99e39557176f553d",
-  "0x16c9e55061ebae204ba4cc8bd75a079432ae2a1d0b7c9dce1665d51c640fcba2"
-)
-#  SNR^((p-1)/2)
-const FrobConst_BN254_Snarks_psi1_coef3 = Fp2[BN254_Snarks].fromHex(
-  "0x63cf305489af5dcdc5ec698b6e2f9b9dbaae0eda9c95998dc54014671a0135a",
-  "0x7c03cbcac41049a0704b5a7ec796f2b21807dc98fa25bd282d37f632623b0e3"
-)
-#  norm(SNR)^((p-1)/3)
-const FrobConst_BN254_Snarks_psi2_coef2 = Fp2[BN254_Snarks].fromHex(
-  "0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd48",
-  "0x0"
-)
-
-#   BLS12_377 is a D-Twist: SNR^((p-1)/6)
-const FrobConst_BLS12_377_psi1_coef1 = Fp2[BLS12_377].fromHex(
-  "0x9a9975399c019633c1e30682567f915c8a45e0f94ebc8ec681bf34a3aa559db57668e558eb0188e938a9d1104f2031",
-  "0x0"
-)
-#  SNR^((p-1)/3)
-const FrobConst_BLS12_377_psi1_coef2 = Fp2[BLS12_377].fromHex(
-  "0x9b3af05dd14f6ec619aaf7d34594aabc5ed1347970dec00452217cc900000008508c00000000002",
-  "0x0"
-)
-#  SNR^((p-1)/2)
-const FrobConst_BLS12_377_psi1_coef3 = Fp2[BLS12_377].fromHex(
-  "0x1680a40796537cac0c534db1a79beb1400398f50ad1dec1bce649cf436b0f6299588459bff27d8e6e76d5ecf1391c63",
-  "0x0"
-)
-#  norm(SNR)^((p-1)/3)
-const FrobConst_BLS12_377_psi2_coef2 = Fp2[BLS12_377].fromHex(
-  "0x9b3af05dd14f6ec619aaf7d34594aabc5ed1347970dec00452217cc900000008508c00000000001",
-  "0x0"
-)
-
-#   BLS12_381 is a M-twist: (1/SNR)^((p-1)/6)
-const FrobConst_BLS12_381_psi1_coef1 = Fp2[BLS12_381].fromHex(
-  "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116",
-  "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116"
-)
-#  (1/SNR)^((p-1)/3)
-const FrobConst_BLS12_381_psi1_coef2 = Fp2[BLS12_381].fromHex(
-  "0x0",
-  "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad"
-)
-#  (1/SNR)^((p-1)/2)
-const FrobConst_BLS12_381_psi1_coef3 = Fp2[BLS12_381].fromHex(
-  "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2",
-  "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09"
-)
-#  norm(SNR)^((p-1)/3)
-const FrobConst_BLS12_381_psi2_coef2 = Fp2[BLS12_381].fromHex(
-  "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac",
-  "0x0"
-)
-
-{.experimental: "dynamicBindSym".}
-
-macro frobConst(C: static Curve, psipow, coefpow: static int): untyped =
-  return bindSym("FrobConst_" & $C & "_psi" & $psipow & "_coef" & $coefpow)
+    r = a

 template mulCheckSparse[Fp2](a: var Fp2, b: Fp2) =
  when b.c0.isZero().bool:
@ -129,11 +51,129 @@ template mulCheckSparse[Fp2](a: var Fp2, b: Fp2) =
  else:
    a *= b

+# Frobenius map - on extension fields
+# -----------------------------------------------------------------
+
+{.experimental: "dynamicBindSym".}
+
+macro frobMapConst(C: static Curve, coefpow, frobpow: static int): untyped =
+  return bindSym("FrobMapConst_" & $C & "_coef" & $coefpow & "_pow" & $frobpow)
+
+# SNR^((p-1)/6)^3
+const FrobMapConst_BLS12_381_coef3_pow1 = Fp2[BLS12_381].fromHex(
+    "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09",
+    "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09"
+  )
+
+func frobenius_map*(r: var Fp4, a: Fp4, k: static int = 1) {.inline.} =
+  ## Computes a^(p^k)
+  ## 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(Fp4.C, 3, k)
+
+# ψ (Psi) - Untwist-Frobenius-Twist Endomorphisms on twisted curves
+# -----------------------------------------------------------------
+# TODO: generate those constants via Sage in a Json file
+#       and parse at compile-time
+
+# Constants:
+#   Assuming embedding degree of 12 and a sextic twist
+#   with SNR the sextic non-residue
+#
+#   BN254_Snarks is a D-Twist: SNR^((p-1)/6)
+const FrobPsiConst_BN254_Snarks_psi1_coef1 = Fp2[BN254_Snarks].fromHex(
+  "0x1284b71c2865a7dfe8b99fdd76e68b605c521e08292f2176d60b35dadcc9e470",
+  "0x246996f3b4fae7e6a6327cfe12150b8e747992778eeec7e5ca5cf05f80f362ac"
+)
+#  SNR^((p-1)/3)
+const FrobPsiConst_BN254_Snarks_psi1_coef2 = Fp2[BN254_Snarks].fromHex(
+  "0x2fb347984f7911f74c0bec3cf559b143b78cc310c2c3330c99e39557176f553d",
+  "0x16c9e55061ebae204ba4cc8bd75a079432ae2a1d0b7c9dce1665d51c640fcba2"
+)
+#  SNR^((p-1)/2)
+const FrobPsiConst_BN254_Snarks_psi1_coef3 = Fp2[BN254_Snarks].fromHex(
+  "0x63cf305489af5dcdc5ec698b6e2f9b9dbaae0eda9c95998dc54014671a0135a",
+  "0x7c03cbcac41049a0704b5a7ec796f2b21807dc98fa25bd282d37f632623b0e3"
+)
+#  norm(SNR)^((p-1)/3)
+const FrobPsiConst_BN254_Snarks_psi2_coef2 = Fp2[BN254_Snarks].fromHex(
+  "0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd48",
+  "0x0"
+)
+
+#   BN254_Nogami is a D-Twist: SNR^((p-1)/6)
+const FrobPsiConst_BN254_Nogami_psi1_coef1 = Fp2[BN254_Nogami].fromHex(
+  "0x1b377619212e7c8cb6499b50a846953f850974924d3f77c2e17de6c06f2a6de9",
+  "0x9ebee691ed1837503eab22f57b96ac8dc178b6db2c08850c582193f90d5922a"
+)
+#  SNR^((p-1)/3)
+const FrobPsiConst_BN254_Nogami_psi1_coef2 = Fp2[BN254_Nogami].fromHex(
+  "0x0",
+  "0x25236482400000017080eb4000000006181800000000000cd98000000000000b"
+)
+#  SNR^((p-1)/2)
+const FrobPsiConst_BN254_Nogami_psi1_coef3 = Fp2[BN254_Nogami].fromHex(
+  "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5",
+  "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5"
+)
+#  norm(SNR)^((p-1)/3)
+const FrobPsiConst_BN254_Nogami_psi2_coef2 = Fp2[BN254_Nogami].fromHex(
+  "0x49b36240000000024909000000000006cd80000000000007",
+  "0x0"
+)
+
+#   BLS12_377 is a D-Twist: SNR^((p-1)/6)
+const FrobPsiConst_BLS12_377_psi1_coef1 = Fp2[BLS12_377].fromHex(
+  "0x9a9975399c019633c1e30682567f915c8a45e0f94ebc8ec681bf34a3aa559db57668e558eb0188e938a9d1104f2031",
+  "0x0"
+)
+#  SNR^((p-1)/3)
+const FrobPsiConst_BLS12_377_psi1_coef2 = Fp2[BLS12_377].fromHex(
+  "0x9b3af05dd14f6ec619aaf7d34594aabc5ed1347970dec00452217cc900000008508c00000000002",
+  "0x0"
+)
+#  SNR^((p-1)/2)
+const FrobPsiConst_BLS12_377_psi1_coef3 = Fp2[BLS12_377].fromHex(
+  "0x1680a40796537cac0c534db1a79beb1400398f50ad1dec1bce649cf436b0f6299588459bff27d8e6e76d5ecf1391c63",
+  "0x0"
+)
+#  norm(SNR)^((p-1)/3)
+const FrobPsiConst_BLS12_377_psi2_coef2 = Fp2[BLS12_377].fromHex(
+  "0x9b3af05dd14f6ec619aaf7d34594aabc5ed1347970dec00452217cc900000008508c00000000001",
+  "0x0"
+)
+
+#   BLS12_381 is a M-twist: (1/SNR)^((p-1)/6)
+const FrobPsiConst_BLS12_381_psi1_coef1 = Fp2[BLS12_381].fromHex(
+  "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116",
+  "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116"
+)
+#  (1/SNR)^((p-1)/3)
+const FrobPsiConst_BLS12_381_psi1_coef2 = Fp2[BLS12_381].fromHex(
+  "0x0",
+  "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad"
+)
+#  (1/SNR)^((p-1)/2)
+const FrobPsiConst_BLS12_381_psi1_coef3 = Fp2[BLS12_381].fromHex(
+  "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2",
+  "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09"
+)
+#  norm(SNR)^((p-1)/3)
+const FrobPsiConst_BLS12_381_psi2_coef2 = Fp2[BLS12_381].fromHex(
+  "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac",
+  "0x0"
+)
+
+macro frobPsiConst(C: static Curve, psipow, coefpow: static int): untyped =
+  return bindSym("FrobPsiConst_" & $C & "_psi" & $psipow & "_coef" & $coefpow)
+
 func frobenius_psi*[PointG2](r: var PointG2, P: PointG2) =
  ## "Untwist-Frobenius-Twist" endomorphism
  ## r = ψ(P)
  for coordR, coordP in fields(r, P):
-    coordR.frobenius(coordP, 1)
+    coordR.frobenius_map(coordP, 1)

  # With ξ (xi) the sextic non-residue
  # c = ξ^((p-1)/6) for D-Twist
@ -142,14 +182,14 @@ func frobenius_psi*[PointG2](r: var PointG2, P: PointG2) =
  # c1_2 = c²
  # c1_3 = c³

-  r.x.mulCheckSparse frobConst(PointG2.F.C, psipow=1, coefpow=2)
-  r.y.mulCheckSparse frobConst(PointG2.F.C, psipow=1, coefpow=3)
+  r.x.mulCheckSparse frobPsiConst(PointG2.F.C, psipow=1, coefpow=2)
+  r.y.mulCheckSparse frobPsiConst(PointG2.F.C, psipow=1, coefpow=3)

 func frobenius_psi2*[PointG2](r: var PointG2, P: PointG2) =
  ## "Untwist-Frobenius-Twist" endomorphism applied twice
  ## r = ψ(ψ(P))
  for coordR, coordP in fields(r, P):
-    coordR.frobenius(coordP, 2)
+    coordR.frobenius_map(coordP, 2)

  # With ξ (xi) the sextic non-residue
  # c = ξ for D-Twist
@ -181,5 +221,5 @@ func frobenius_psi2*[PointG2](r: var PointG2, P: PointG2) =
  # c2_3 ≡ -1 (mod p²)
  # QED

-  r.x.mulCheckSparse frobConst(PointG2.F.C, psipow=2, coefpow=2)
+  r.x.mulCheckSparse frobPsiConst(PointG2.F.C, psipow=2, coefpow=2)
  r.y.neg(r.y)
--- a/constantine/pairing/README.md
+++ b/constantine/pairing/README.md
@ -15,10 +15,28 @@
  Ben Lynn, 2007\
  https://crypto.stanford.edu/pbc/thesis.pdf

+- Faster Pairing Computations on Curves with High-Degree Twists
+  Craig Costello, Tanja Lange, and Michael Naehrig, 2009
+  https://eprint.iacr.org/2009/615.pdf
+
+- High-Speed Software Implementation of the Optimal Ate Pairing over Barreto-Naehrig Curves\
+  Jean-Luc Beuchat and Jorge Enrique González Díaz and Shigeo Mitsunari and Eiji Okamoto and Francisco Rodríguez-Henríquez and Tadanori Teruya, 2010\
+  https://eprint.iacr.org/2010/354
+
 - An Analysis of Affine Coordinates for Pairing Computation\
  Kristin Lauter, Peter L. Montgomery, and Michael Naehrig, 2010\
  https://eprint.iacr.org/2010/363.pdf

+- Faster Explicit Formulas for Computing Pairings over Ordinary Curves\
+  Diego F. Aranha and Koray Karabina and Patrick Longa and Catherine H. Gebotys and Julio López, 2010\
+  https://eprint.iacr.org/2010/526.pdf\
+  https://www.iacr.org/archive/eurocrypt2011/66320047/66320047.pdf
+
+- Avoiding Full Extension Field Arithmetic in Pairing Computations
+  Craig Costello, Colin Boyd,
+  Juan Manuel Gonzalez Nieto, and Kenneth Koon-Ho Wong, 2010
+  https://eprint.iacr.org/2010/104.pdf
+
 - Pairings for beginners\
  Craig Costello, 2012 (?)\
  http://www.craigcostello.com.au/pairings/PairingsForBeginners.pdf
@ -28,6 +46,16 @@
  Craig Costello, 2012\
  https://eprints.qut.edu.au/61037/1/Craig_Costello_Thesis.pdf

+- Efficient Implementation of Bilinear Pairings on ARM Processors
+  Gurleen Grewal, Reza Azarderakhsh,
+  Patrick Longa, Shi Hu, and David Jao, 2012
+  https://eprint.iacr.org/2012/408.pdf
+
+- The Realm of the Pairings\
+  Diego F. Aranha and Paulo S. L. M. Barreto and Patrick Longa and Jefferson E. Ricardini\
+  https://eprint.iacr.org/2013/722.pdf\
+  http://sac2013.irmacs.sfu.ca/slides/s1.pdf
+
 - Efficient Implementations of Pairing-Based Cryptography on Embedded Systems\
  Master Thesis\
  Rajeev Verma, 2015\
@ -38,10 +66,21 @@
  Razvan Barbulescu, Nadia El Mrabet, and Loubna Ghammam, 2019\
  https://hal.archives-ouvertes.fr/hal-02129868/file/2019-485.pdf

- A short-list of pairing-friendly curves resistantto Special TNFS at the 128-bit security level\
-  Aurore Guillevic\
+- Pairing Implementation Revisited\
+  Mike Scott, 2019\
+  https://eprint.iacr.org/2019/077.pdf
+
+- A short-list of pairing-friendly curves resistant to Special TNFS at the 128-bit security level\
+  Aurore Guillevic, 2019\
  https://eprint.iacr.org/2019/1371.pdf

+- 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
+
 ### Presentations

 - Introduction to pairings\
--- a/constantine/pairing/gt_fp12.nim
+++ b/constantine/pairing/gt_fp12.nim
@ -0,0 +1,168 @@
+# 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,
+  ./lines_projective
+
+# ############################################################
+#
+#                 Sparse Multiplication
+#                and cyclotomic squaring
+#              for elements of Gₜ = E(Fp¹²)
+#
+# ############################################################
+
+# - Pairing Implementation Revisited
+#   Michael Scott, 2019
+#   https://eprint.iacr.org/2019/077
+#
+# - Efficient Implementation of Bilinear Pairings on ARM Processors
+#   Gurleen Grewal, Reza Azarderakhsh,
+#   Patrick Longa, Shi Hu, and David Jao, 2012
+#   https://eprint.iacr.org/2012/408.pdf
+#
+# - High-Speed Software Implementation of the Optimal Ate Pairing over Barreto-Naehrig Curves\
+#   Jean-Luc Beuchat and Jorge Enrique González Díaz and Shigeo Mitsunari and Eiji Okamoto and Francisco Rodríguez-Henríquez and Tadanori Teruya, 2010\
+#   https://eprint.iacr.org/2010/354
+#
+# - Faster Pairing Computations on Curves with High-Degree Twists
+#   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
+# ----------------------------------------------------------------
+
+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)
+  # TODO: make generic and move to tower_field_extensions
+
+  # v0 = a0 b0 = 0
+  # v1 = a1 b1
+  # v2 = a2 b2 = 0
+  #
+  # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
+  #    = ξ (a1 b1 + a2 b1 - v1)
+  #    = ξ a2 b1
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + β v2
+  #    = a0 b1 + a1 b1 - v1
+  #    = a0 b1
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
+  #    = v1
+  #    = a1 b1
+
+  r.c0.prod(a.c2, b)
+  r.c0 *= ξ
+  r.c1.prod(a.c0, b)
+  r.c2.prod(a.c1, b)
+
+func mul_by_line_xy0*[C: static Curve, twist: static SexticTwist](
+       r: var Fp6[C],
+       a: Fp6[C],
+       b: Line[Fp2[C], twist]) =
+  ## Sparse multiplication of an 𝔽p6
+  ## with coordinates (a₀, a₁, a₂) by a line (x, y, 0)
+  ## The z coordinates in the line will be ignored.
+  ## `r` and `a` must not alias
+  var
+    v0 {.noInit.}: Fp2[C]
+    v1 {.noInit.}: Fp2[C]
+
+  v0.prod(a.c0, b.x)
+  v1.prod(a.c1, b.y)
+  r.c0.prod(a.c2, b.y)
+  r.c0 *= ξ
+  r.c0 += v0
+
+  r.c1.sum(a.c0, a.c1) # Error when r and a alias as r.c0 was updated
+  r.c2.sum(b.x, b.y)
+  r.c1 *= r.c2
+  r.c1 -= v0
+  r.c1 -= v1
+
+  r.c2.prod(a.c2, b.x)
+  r.c2 += v1
+
+func mul_sparse_by_line_xy00z0*[C: static Curve, Tw: static SexticTwist](
+      f: var Fp12[C], l: Line[Fp2[C], Tw]) =
+  ## Sparse multiplication of an 𝔽p12 element
+  ## by a sparse 𝔽p12 element coming from an D-Twist line function.
+  ## The sparse element is represented by a packed Line type
+  ## with coordinate (x,y,z) matching 𝔽p12 coordinates xy00z0 (TODO: verify this)
+
+  static: doAssert f.c0.typeof is Fp6, "This assumes 𝔽p12 as a quadratic extension of 𝔽p6"
+
+  var
+    v0 {.noInit.}: Fp6[C]
+    v1 {.noInit.}: Fp6[C]
+    v2 {.noInit.}: Line[Fp2[C], Tw]
+    v3 {.noInit.}: Fp6[C]
+
+  v0.mul_by_line_xy0(f.c0, l)
+  v1.mul_sparse_by_0y0(f.c1, l.z)
+
+  v2.x = l.x
+  v2.y.sum(l.y, l.z)
+  f.c1 += f.c0
+  v3.mul_by_line_xy0(f.c1, v2)
+  v3 -= v0
+  v3 -= v1
+  f.c1 = v3
+
+  v3.c0 = ξ * v1.c2
+  v3.c0 += v0.c0
+  v3.c1.sum(v0.c1, v1.c0)
+  v3.c2.sum(v0.c2, v1.c1)
+  f.c0 = v3
+
+func mul_sparse_by_line_xyz000*[C: static Curve, Tw: static SexticTwist](
+       f: var Fp12[C], l: Line[Fp2[C], Tw]) =
+  ## Sparse multiplication of an 𝔽p12 element
+  ## by a sparse 𝔽p12 element coming from an D-Twist line function.
+  ## The sparse element is represented by a packed Line type
+  ## with coordinates (x,y,z) matching 𝔽p12 coordinates xyz000
+
+  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.c1.c0 = l.z
+
+  f *= v
+
+func mul_sparse_by_line_xy000z*[C: static Curve, Tw: static SexticTwist](
+       f: var Fp12[C], l: Line[Fp2[C], Tw]) =
+
+  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
+
+  f *= v
+
+# 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
--- a/constantine/pairing/lines_projective.nim
+++ b/constantine/pairing/lines_projective.nim
@ -0,0 +1,234 @@
+# 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
+  std/typetraits,
+  ../primitives,
+  ../config/[common, curves],
+  ../arithmetic,
+  ../towers,
+  ../elliptic/[
+    ec_weierstrass_affine,
+    ec_weierstrass_projective
+  ],
+  ../io/io_towers
+
+# ############################################################
+#
+#            Miller Loop's Line Evaluation
+#             with projective coordinates
+#
+# ############################################################
+#
+# - Pairing Implementation Revisited
+#   Michael Scott, 2019
+#   https://eprint.iacr.org/2019/077
+#
+# - The Realm of the Pairings
+#   Diego F. Aranha and Paulo S. L. M. Barreto
+#   and Patrick Longa and Jefferson E. Ricardini, 2013
+#   https://eprint.iacr.org/2013/722.pdf
+#   http://sac2013.irmacs.sfu.ca/slides/s1.pdf
+
+#
+# TODO: Implement fused line doubling and addition
+#       from Costello2009 or Aranha2010
+#       We don't need the complete formulae in the Miller Loop
+
+type
+  Line*[F; twist: static SexticTwist] = object
+    ## Packed line representation over a E'(Fp^k/d)
+    ## with k the embedding degree and d the twist degree
+    ## i.e. for a curve with embedding degree 12 and sextic twist
+    ## F is Fp2
+    ##
+    ## Assuming a Sextic Twist
+    ##
+    ## Out of 6 Fp2 coordinates, 3 are 0 and
+    ## the non-zero coordinates depend on the twist kind.
+    ##
+    ## For a D-twist,
+    ##   (x, y, z) corresponds to an sparse element of Fp12
+    ##   with Fp2 coordinates: xy00z0
+    ## For a M-Twist
+    ##   (x, y, z) corresponds to an sparse element of Fp12
+    ##   with Fp2 coordinates: xyz000
+    x*, y*, z*: F
+
+func toHex*(line: Line, order: static Endianness = bigEndian): string =
+  result = static($line.typeof.genericHead() & '(')
+  for fieldName, fieldValue in fieldPairs(line):
+    when fieldName != "x":
+      result.add ", "
+    result.add fieldName & ": "
+    result.appendHex(fieldValue, order)
+  result.add ")"
+
+# Line evaluation
+# --------------------------------------------------
+
+func `*=`(a: var Fp2, b: Fp) =
+  ## Multiply an element of Fp2 by an element of Fp
+  # TODO: make generic and move to tower_field_extensions
+  a.c0 *= b
+  a.c1 *= b
+
+func line_update(line: var Line, P: ECP_SWei_Aff) =
+  ## Update the line evaluation with P
+  ## after addition or doubling
+  ## P in G1
+  line.x *= P.y
+  line.z *= P.x
+
+func line_eval_double*(line: var Line, T: ECP_SWei_Proj) =
+  ## Evaluate the line function for doubling
+  ## i.e. the tangent at T
+  ##
+  ## With T in homogenous projective coordinates (X, Y, Z)
+  ## And ξ the sextic non residue to construct 𝔽p4 / 𝔽p6 / 𝔽p12
+  ##
+  ## M-Twist:
+  ##   A = -2ξ Y.Z
+  ##   B = 3bξ Z² - Y²
+  ##   C = 3 X²
+  ##
+  ## D-Twist are scaled by ξ to avoid dividing by ξ:
+  ##   A = -2ξ Y.Z
+  ##   B = 3b Z² - ξY²
+  ##   C = 3ξ X²
+  ##
+  ## Instead of
+  ##   - equation 10 from The Real of pairing, Aranha et al, 2013
+  ##   - or chapter 3 from pairing Implementation Revisited, Scott 2019
+  ##   A = -2 Y.Z
+  ##   B = 3b/ξ Z² - Y²
+  ##   C = 3 X²
+  ##
+  ## A constant factor will be wiped by the final exponentiation
+  ## as for all non-zero α ∈ GF(pᵐ)
+  ## with
+  ## - p odd prime
+  ## - and gcd(α,pᵐ) = 1 (i.e. the extension field pᵐ is using irreducible polynomials)
+  ##
+  ## Little Fermat holds and we have
+  ## α^(pᵐ - 1) ≡ 1 (mod pᵐ)
+  ##
+  ## The final exponent is of the form
+  ## (pᵏ-1)/r
+  ##
+  ## A constant factor on twisted coordinates pᵏᐟᵈ
+  ## is a constant factor on pᵏ with d the twisting degree
+  ## and so will be elminated. QED.
+  var v {.noInit.}: Line.F
+  const b3 = 3 * ECP_SWei_Proj.F.C.getCoefB()
+
+  template A: untyped = line.x
+  template B: untyped = line.y
+  template C: untyped = line.z
+
+  A = T.y               # A = Y
+  v = T.y               # v = Y
+  B = T.z               # B = Z
+  C = T.x               # C = X
+
+  A *= B                # A = Y.Z
+  C.square()            # C = X²
+  v.square()            # v = Y²
+  B.square()            # B = Z²
+
+  A.double()            # A =  2 Y.Z
+  A.neg()               # A = -2 Y.Z
+  A *= SexticNonResidue # A = -2 ξ Y.Z
+
+  B *= b3               # B = 3b Z²
+  C *= 3                # C = 3X²
+  when ECP_SWei_Proj.F.C.getSexticTwist() == M_Twist:
+    B *= SexticNonResidue # B = 3b' Z² = 3bξ Z²
+  elif ECP_SWei_Proj.F.C.getSexticTwist() == D_Twist:
+    v *= SexticNonResidue # v =  ξ Y²
+    C *= SexticNonResidue # C = 3ξ X²
+  else:
+    {.error: "unreachable".}
+
+  B -= v                # B = 3bξ Z² - Y²  (M-twist)
+                        # B = 3b Z² - ξ Y² (D-twist)
+
+func line_eval_add*(line: var Line, T: ECP_SWei_Proj, Q: ECP_SWei_Aff) =
+  ## Evaluate the line function for addition
+  ## i.e. the line between T and Q
+  ##
+  ## With T in homogenous projective coordinates (X, Y, Z)
+  ## And ξ the sextic non residue to construct 𝔽p4 / 𝔽p6 / 𝔽p12
+  ##
+  ## M-Twist:
+  ##   A = ξ (X₁ - Z₁X₂)
+  ##   B = (Y₁ - Z₁Y₂) X₂ - (X₁ - Z₁X₂) Y₂
+  ##   C = - (Y₁ - Z₁Y₂)
+  ##
+  ## D-Twist:
+  ##   A = X₁ - Z₁X₂
+  ##   B = (Y₁ - Z₁Y₂) X₂ - (X₁ - Z₁X₂) Y₂
+  ##   C = - (Y₁ - Z₁Y₂)
+  ##
+  ## Note: There is no need for complete formula as
+  ## we have T ∉ [Q, -Q] in the Miller loop doubling-and-add
+  ## i.e. the line cannot be vertical
+  var v {.noInit.}: Line.F
+
+  template A: untyped = line.x
+  template B: untyped = line.y
+  template C: untyped = line.z
+
+  A = T.x     # A = X₁
+  v = T.z     # v = Z₁
+  B = T.z     # B = Z₁
+  C = T.y     # C = Y₁
+
+  v *= Q.y    # v = Z₁Y₂
+  B *= Q.x    # B = Z₁X₂
+
+  A -= B      # A = X₁-Z₁X₂
+  C -= v      # C = Y₁-Z₁Y₂
+
+  v = A       # v = X₁-Z₁X₂
+  when ECP_SWei_Proj.F.C.getSexticTwist() == M_Twist:
+    A *= SexticNonResidue # A = ξ (X₁ - Z₁X₂)
+
+  v *= Q.y    # v = (X₁-Z₁X₂) Y₂
+  B = C       # B = Y₁-Z₁Y₂
+  B *= Q.x    # B = (Y₁-Z₁Y₂) X₂
+  B -= v      # B = (Y₁-Z₁Y₂) X₂ - (X₁-Z₁X₂) Y₂
+
+  C.neg()     # C = -(Y₁-Z₁Y₂)
+
+func line_double*(line: var Line, T: var ECP_SWei_Proj, P: ECP_SWei_Aff) =
+  ## Doubling step of the Miller loop
+  ## T in G2, P in G1
+  ##
+  ## Compute lt,t(P)
+  ##
+  # TODO fused line doubling from Costello 2009, Grewal 2012, Aranha 2013
+  line_eval_double(line, T)
+  line.line_update(P)
+  T.double()
+
+func line_add*[C](
+       line: var Line,
+       T: var ECP_SWei_Proj[Fp2[C]],
+       Q: ECP_SWei_Aff[Fp2[C]], P: ECP_SWei_Aff[Fp[C]]) =
+  ## Addition step of the Miller loop
+  ## T and Q in G2, P in G1
+  ##
+  ## Compute lt,q(P)
+  # TODO fused line addition from Costello 2009, Grewal 2012, Aranha 2013
+  line_eval_add(line, T, Q)
+  line.line_update(P)
+  # TODO: mixed addition
+  var QProj {.noInit.}: ECP_SWei_Proj[Fp2[C]]
+  QProj.projectiveFromAffine(Q)
+  T += QProj
--- a/constantine/pairing/pairing_bls12.nim
+++ b/constantine/pairing/pairing_bls12.nim
@ -0,0 +1,152 @@
+# 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,
+  ../io/io_bigints,
+  ../elliptic/[
+    ec_weierstrass_affine,
+    ec_weierstrass_projective
+  ],
+  ./lines_projective,
+  ./gt_fp12
+
+# ############################################################
+#
+#                 Optimal ATE pairing for
+#                      BLS12 curves
+#
+# ############################################################
+
+# - 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
+#
+# - Faster Pairing Computations on Curves with High-Degree Twists
+#   Craig Costello, Tanja Lange, and Michael Naehrig, 2009
+#   https://eprint.iacr.org/2009/615.pdf
+
+# TODO: implement quadruple-and-add and octuple-and-add
+#       from Costello2009 to trade multiplications in Fpᵏ
+#       for multiplications in Fp
+
+# TODO: should be part of curve parameters
+const BLS12_381_param = block:
+  # BLS Miller loop is parametrized by u
+  BigInt[64+2].fromHex("0xd201000000010000") # +2 so that we can take *3 and NAF encode it
+
+const BLS12_381_param_isNeg = true
+
+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](
+       f: var Fp12[C],
+       P: ECP_SWei_Aff[Fp[C]],
+       Q: ECP_SWei_Aff[Fp2[C]]
+     ) =
+  ## Generic Miller Loop for BLS12 curve
+  ## Computes f{u,Q}(P) with u the BLS curve parameter
+  # TODO: retrieve the curve parameter from the curve declaration
+
+  # Boundary cases
+  #   Loop start
+  #     The litterature starts from both L-1 or L-2:
+  #     L-1:
+  #     - Scott2019, Pairing Implementation Revisited, Algorithm 1
+  #     - Aranha2010, Faster Explicit Formulas ..., Algorithm 1
+  #     L-2
+  #     - Beuchat2010, High-Speed Software Implementation ..., Algorithm 1
+  #     - Aranha2013, The Realm of The Pairings, Algorithm 1
+  #     - Costello, Thesis, Algorithm 2.1
+  #     - Costello2012, Pairings for Beginners, Algorithm 5.1
+  #
+  #     Even the guide to pairing based cryptography has both
+  #     Chapter 3: L-1 (Algorithm 3.1)
+  #     Chapter 11: L-2 (Algorithm 11.1) but it explains why L-2 (unrolling)
+  #  Loop end
+  #    - Some implementation, for example Beuchat2010 or the Guide to Pairing-Based Cryptography
+  #      have extra line additions after the main loop,
+  #      this is needed for BN curves.
+  #    - With r the order of G1 / G2 / GT,
+  #      we have [r]T = Inf
+  #      Hence, [r-1]T = -T
+  #      so either we use complete addition
+  #      or we special case line addition of T and -T (it's a vertical line)
+  #      or we ensure the loop is done for a number of iterations strictly less
+  #      than the curve order which is the case for BLS12 curves
+
+  static: doAssert C == BLS12_381, "Only BLS12-381 is supported at the moment"
+
+  var
+    T {.noInit.}: ECP_SWei_Proj[Fp2[C]]
+    line {.noInit.}: Line[Fp2[C], C.getSexticTwist()]
+    nQ{.noInit.}: typeof(Q)
+
+  T.projectiveFromAffine(Q)
+  nQ.neg(Q)
+  f.setOne()
+
+  template u: untyped = BLS12_381_param
+  let u3 = 3*BLS12_381_param
+  for i in countdown(u3.bits - 2, 1):
+    f.square()
+    line.line_double(T, P)
+    f.mul_sparse_by_line_xy000z(line)
+
+    let naf = u3.bit(i).int8 - u.bit(i).int8 # This can throw exception
+    if naf == 1:
+      line.line_add(T, Q, P)
+      f.mul_sparse_by_line_xy000z(line)
+    elif naf == -1:
+      line.line_add(T, nQ, P)
+      f.mul_sparse_by_line_xy000z(line)
+
+  when BLS12_381_param_isNeg: # TODO generic
+    # In GT, x^-1 == conjugate(x)
+    # Remark 7.1, chapter 7.1.1 of Guide to Pairing-Based Cryptography, El Mrabet, 2017
+    f.conj()
+
+func finalExpGeneric[C: static Curve](f: var Fp12[C]) =
+  ## A generic and slow implementation of final exponentiation
+  ## for sanity checks purposes.
+  static: doAssert C == BLS12_381, "Only BLS12-381 is supported at the moment"
+  f.powUnsafeExponent(BLS12_381_finalexponent, window = 3)
+
+func pairing_bls12_reference*[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
+  ##
+  ## 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.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)
+  ##
+  ## With an embedding degree of 12, the easy part of final exponentiation is
+  ##
+  ##  f^(p⁶−1)(p²+1)
+  discard
--- a/constantine/pairing/pairing_bn.nim
+++ b/constantine/pairing/pairing_bn.nim
@ -0,0 +1,178 @@
+# 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
+  std/macros,
+  ../primitives,
+  ../config/[common, curves],
+  ../arithmetic,
+  ../towers,
+  ../io/io_bigints,
+  ../elliptic/[
+    ec_weierstrass_affine,
+    ec_weierstrass_projective
+  ],
+  ./lines_projective,
+  ./gt_fp12,
+  ../isogeny/frobenius
+
+# ############################################################
+#
+#                 Optimal ATE pairing for
+#                      BN curves
+#
+# ############################################################
+
+# - 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
+#
+# - Faster Pairing Computations on Curves with High-Degree Twists
+#   Craig Costello, Tanja Lange, and Michael Naehrig, 2009
+#   https://eprint.iacr.org/2009/615.pdf
+
+# TODO: implement quadruple-and-add and octuple-and-add
+#       from Costello2009 to trade multiplications in Fpᵏ
+#       for multiplications in Fp
+
+# TODO: should be part of curve parameters
+const BN254_Snarks_ate_param = block:
+  # BN Miller loop is parametrized by 6u+2
+  BigInt[67].fromHex"0x19d797039be763ba8"
+
+const BN254_Snarks_ate_param_isNeg = false
+
+const BN254_Snarks_finalexponent = block:
+  # (p^12 - 1) / r
+  BigInt[2790].fromHex"0x2f4b6dc97020fddadf107d20bc842d43bf6369b1ff6a1c71015f3f7be2e1e30a73bb94fec0daf15466b2383a5d3ec3d15ad524d8f70c54efee1bd8c3b21377e563a09a1b705887e72eceaddea3790364a61f676baaf977870e88d5c6c8fef0781361e443ae77f5b63a2a2264487f2940a8b1ddb3d15062cd0fb2015dfc6668449aed3cc48a82d0d602d268c7daab6a41294c0cc4ebe5664568dfc50e1648a45a4a1e3a5195846a3ed011a337a02088ec80e0ebae8755cfe107acf3aafb40494e406f804216bb10cf430b0f37856b42db8dc5514724ee93dfb10826f0dd4a0364b9580291d2cd65664814fde37ca80bb4ea44eacc5e641bbadf423f9a2cbf813b8d145da90029baee7ddadda71c7f3811c4105262945bba1668c3be69a3c230974d83561841d766f9c9d570bb7fbe04c7e8a6c3c760c0de81def35692da361102b6b9b2b918837fa97896e84abb40a4efb7e54523a486964b64ca86f120"
+
+const BN254_Nogami_ate_param = block:
+  # BN Miller loop is parametrized by 6u+2
+  BigInt[67].fromHex"0x18300000000000004" # 65+2 bit for NAF x3 encoding
+
+const BN254_Nogami_ate_param_isNeg = true
+
+const BN254_Nogami_finalexponent = block:
+  # (p^12 - 1) / r
+  BigInt[2786].fromHex"0x2928fbb36b391596ee3fe4cbe857330da83e46fedf04d235a4a8daf5ff9f6eabcb4e3f20aa06f0a0d96b24f9af0cbbce750d61627dcbf5fec9139b8f1c46c86b49b4f8a202af26e4504f2c0f56570e9bd5b94c403f385d1908556486e24b396ddc2cdf13d06542f84fe8e82ccbad7b7423fc1ef4e8cc73d605e3e867c0a75f45ea7f6356d9846ce35d5a34f30396938818ad41914b97b99c289a7259b5d2e09477a77bd3c409b19f19e893f8ade90b0aed1b5fc8a07a3cebb41d4e9eee96b21a832ddb1e93e113edfb704fa532848c18593cd0ee90444a1b3499a800177ea38bdec62ec5191f2b6bbee449722f98d2173ad33077545c2ad10347e125a56fb40f086e9a4e62ad336a72c8b202ac3c1473d73b93d93dc0795ca0ca39226e7b4c1bb92f99248ec0806e0ad70744e9f2238736790f5185ea4c70808442a7d530c6ccd56b55a6973867ec6c73599bbd020bbe105da9c6b5c009ad8946cd6f0"
+
+{.experimental: "dynamicBindSym".}
+
+macro get(C: static Curve, value: untyped): untyped =
+  return bindSym($C & "_" & $value)
+
+func millerLoopGenericBN[C: static Curve](
+       f: var Fp12[C],
+       P: ECP_SWei_Aff[Fp[C]],
+       Q: ECP_SWei_Aff[Fp2[C]]
+     ) =
+  ## Generic Miller Loop for BN curves
+  ## Computes f{6u+2,Q}(P) with u the BN curve parameter
+  # TODO: retrieve the curve parameter from the curve declaration
+
+  # TODO - boundary cases
+  #   Loop start
+  #     The literatture starts from both L-1 or L-2:
+  #     L-1:
+  #     - Scott2019, Pairing Implementation Revisited, Algorithm 1
+  #     - Aranha2010, Faster Explicit Formulas ..., Algorithm 1
+  #     L-2
+  #     - Beuchat2010, High-Speed Software Implementation ..., Algorithm 1
+  #     - Aranha2013, The Realm of The Pairings, Algorithm 1
+  #     - Costello, Thesis, Algorithm 2.1
+  #     - Costello2012, Pairings for Beginners, Algorithm 5.1
+  #
+  #     Even the guide to pairing based cryptography has both
+  #     Chapter 3: L-1 (Algorithm 3.1)
+  #     Chapter 11: L-2 (Algorithm 11.1) but it explains why L-2 (unrolling)
+  #  Loop end
+  #    - Some implementation, for example Beuchat2010 or the Guide to Pairing-Based Cryptography
+  #      have an extra line addition after the main loop, this seems related to
+  #      the NAF recoding and not Miller Loop
+  #    - With r the order of G1 / G2 / GT,
+  #      we have [r]T = Inf
+  #      Hence, [r-1]T = -T
+  #      so either we use complete addition
+  #      or we special case line addition of T and -T (it's a vertical line)
+  #      or we ensure the loop is done for a number of iterations strictly less
+  #      than the curve order which is the case for BN curves
+
+  var
+    T {.noInit.}: ECP_SWei_Proj[Fp2[C]]
+    line {.noInit.}: Line[Fp2[C], C.getSexticTwist()]
+    nQ{.noInit.}: typeof(Q)
+
+  T.projectiveFromAffine(Q)
+  nQ.neg(Q)
+  f.setOne()
+
+  template u: untyped = C.get(ate_param)
+  let u3 = 3*C.get(ate_param)
+  for i in countdown(u3.bits - 2, 1):
+    f.square()
+    line.line_double(T, P)
+    f.mul_sparse_by_line_xyz000(line)
+
+    let naf = u3.bit(i).int8 - u.bit(i).int8 # This can throw exception
+    if naf == 1:
+      line.line_add(T, Q, P)
+      f.mul_sparse_by_line_xyz000(line)
+    elif naf == -1:
+      line.line_add(T, nQ, P)
+      f.mul_sparse_by_line_xyz000(line)
+
+  when C.get(ate_param_isNeg): # TODO generic
+    # In GT, x^-1 == conjugate(x)
+    # Remark 7.1, chapter 7.1.1 of Guide to Pairing-Based Cryptography, El Mrabet, 2017
+    f.conj()
+
+  # Ate pairing for BN curves need adjustment after Miller loop
+  when C.get(ate_param_isNeg):
+    T.neg()
+  var V {.noInit.}: typeof(Q)
+
+  V.frobenius_psi(Q)
+  line.line_add(T, V, P)
+  f.mul_sparse_by_line_xyz000(line)
+
+  V.frobenius_psi2(Q)
+  V.neg()
+  line.line_add(T, V, P)
+  f.mul_sparse_by_line_xyz000(line)
+
+func finalExpGeneric[C: static Curve](f: var Fp12[C]) =
+  ## A generic and slow implementation of final exponentiation
+  ## for sanity checks purposes.
+  f.powUnsafeExponent(C.get(finalexponent), window = 3)
+
+func pairing_bn_reference*[C](gt: var Fp12[C], P: ECP_SWei_Proj[Fp[C]], Q: ECP_SWei_Proj[Fp2[C]]) =
+  ## Compute the optimal Ate Pairing for BN curves
+  ## Input: P ∈ G1, Q ∈ G2
+  ## Output: e(P, Q) ∈ Gt
+  ##
+  ## 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.millerLoopGenericBN(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)
+  ##
+  ## With an embedding degree of 12, the easy part of final exponentiation is
+  ##
+  ##  f^(p⁶−1)(p²+1)
+  discard
--- a/constantine/tower_field_extensions/cubic_extensions.nim
+++ b/constantine/tower_field_extensions/cubic_extensions.nim
@ -190,9 +190,31 @@ func inv*(r: var CubicExt, a: CubicExt) =
  r.c1.prod(v1, v3)
  r.c2.prod(v2, v3)

+func inv*(a: var CubicExt) =
+  ## Compute the multiplicative inverse of ``a``
+  ##
+  ## The inverse of 0 is 0.
+  ## Incidentally this avoids extra check
+  ## to convert Jacobian and Projective coordinates
+  ## to affine for elliptic curve
+  let t = a
+  a.inv(t)
+
 func `*=`*(a: var CubicExt, b: CubicExt) {.inline.} =
  ## In-place multiplication
  # On higher extension field like 𝔽p6,
  # if `prod` is called on shared in and out buffer, the result is wrong
  let t = a
  a.prod(t, b)
+
+func conj*(a: var CubicExt) {.inline.} =
+  ## Computes the conjugate in-place
+  mixin conj, conjneg
+  a.c0.conj()
+  a.c1.conjneg()
+  a.c2.conj()
+
+func square*(a: var CubicExt) {.inline.} =
+  ## In-place squaring
+  let t = a
+  a.square(t)
--- a/constantine/tower_field_extensions/exponentiations.nim
+++ b/constantine/tower_field_extensions/exponentiations.nim
@ -11,6 +11,7 @@ import
  ../arithmetic,
  ../config/[common, curves],
  ../primitives,
+  ../io/io_bigints,
  ./tower_common,
  ./quadratic_extensions,
  ./cubic_extensions
@ -60,7 +61,7 @@ func powPrologue[F](a: var F, scratchspace: var openarray[F]): uint =
  else:
    scratchspace[2] = a
    for k in 2 ..< 1 shl result:
-      scratchspace[k+1]
+      scratchspace[k+1].prod(scratchspace[k], a)
  a.setOne()

 func powSquarings[F](
@ -69,7 +70,7 @@ func powSquarings[F](
       tmp: var F,
       window: uint,
       acc, acc_len: var uint,
-       e: var uint
+       e: var int
     ): tuple[k, bits: uint] {.inline.}=
  ## Squaring step of exponentiation by squaring
  ## Get the next k bits in range [1, window)
@ -109,10 +110,10 @@ func powSquarings[F](

  return (k, bits)

-func powUnsafeExponent(
-       a: var ExtensionField,
+func powUnsafeExponent[F](
+       a: var F,
       exponent: openArray[byte],
-       scratchspace: var openArray[byte]
+       scratchspace: var openArray[F]
     ) =
  ## Extension field exponentiation r = a^exponent (mod p^m)
  ##
@ -145,6 +146,55 @@ func powUnsafeExponent(
        scratchspace[0].prod(a, scratchspace[1])
      a = scratchspace[0]

+func powUnsafeExponent*[F](
+       a: var F,
+       exponent: openArray[byte],
+       window: static int
+     ) =
+  ## Extension field exponentiation r = a^exponent (mod p^m)
+  ## exponent is an big integer in canonical octet-string format
+  ##
+  ## Window is used for window optimization.
+  ## 2^window field elements are allocated for scratchspace.
+  ##
+  ## - On Fp2, with a 256-bit base field, a window of size 5 requires
+  ##   2*256*2^5 = 16KiB
+  ##
+  ## Warning ⚠️ :
+  ## This is an optimization for public exponent
+  ## Otherwise bits of the exponent can be retrieved with:
+  ## - memory access analysis
+  ## - power analysis
+  ## - timing analysis
+  const scratchLen = if window == 1: 2
+                     else: (1 shl window) + 1
+  var scratchSpace {.noInit.}: array[scratchLen, typeof(a)]
+  a.powUnsafeExponent(exponent, scratchspace)
+
+func powUnsafeExponent*[F; bits: static int](
+       a: var F,
+       exponent: BigInt[bits],
+       window: static int
+     ) =
+  ## Extension field exponentiation r = a^exponent (mod p^m)
+  ## exponent is an big integer in canonical octet-string format
+  ##
+  ## Window is used for window optimization.
+  ## 2^window field elements are allocated for scratchspace.
+  ##
+  ## - On Fp2, with a 256-bit base field, a window of size 5 requires
+  ##   2*256*2^5 = 16KiB
+  ##
+  ## Warning ⚠️ :
+  ## This is an optimization for public exponent
+  ## Otherwise bits of the exponent can be retrieved with:
+  ## - memory access analysis
+  ## - power analysis
+  ## - timing analysis
+  var expBE {.noInit.}: array[(bits + 7) div 8, byte]
+  expBE.exportRawUint(exponent, bigEndian)
+  a.powUnsafeExponent(expBE, window)
+
 # Square root
 # -----------------------------------------------------------
 #
--- a/constantine/tower_field_extensions/quadratic_extensions.nim
+++ b/constantine/tower_field_extensions/quadratic_extensions.nim
@ -246,6 +246,19 @@ func mul_sparse_generic_by_x0(r: var QuadraticExt, a, sparseB: QuadraticExt) =
 # Exported symbols
 # -------------------------------------------------------------------

+func conj*(a: var QuadraticExt) {.inline.} =
+  ## Computes the conjugate in-place
+  a.c1.neg()
+
+func conj*(r: var QuadraticExt, a: QuadraticExt) {.inline.} =
+  ## Computes the conjugate out-of-place
+  r.c0 = a.c0
+  r.c1.neg(a.c1)
+
+func conjneg*(a: var QuadraticExt) {.inline.} =
+  ## Computes the negated conjugate in-place
+  a.c0.neg()
+
 func square*(r: var QuadraticExt, a: QuadraticExt) {.inline.} =
  mixin fromComplexExtension
  when r.fromComplexExtension():
@ -294,6 +307,15 @@ func inv*(r: var QuadraticExt, a: QuadraticExt) =
  v0.neg(v1)              # v0 = -1 / (a0² - β a1²)
  r.c1.prod(a.c1, v0)     # r1 = -a1 / (a0² - β a1²)

+func inv*(a: var QuadraticExt) =
+  ## Compute the multiplicative inverse of ``a``
+  ##
+  ## The inverse of 0 is 0.
+  ## Incidentally this avoids extra check
+  ## to convert Jacobian and Projective coordinates
+  ## to affine for elliptic curve
+  a.inv(a)
+
 func `*=`*(a: var QuadraticExt, b: QuadraticExt) {.inline.} =
  ## In-place multiplication
  # On higher extension field like 𝔽p12,
--- a/constantine/tower_field_extensions/tower_common.nim
+++ b/constantine/tower_field_extensions/tower_common.nim
@ -107,6 +107,11 @@ func neg*(r: var ExtensionField, a: ExtensionField) =
  for fR, fA in fields(r, a):
    fR.neg(fA)

+func neg*(a: var ExtensionField) =
+  ## Field out-of-place negation
+  for fA in fields(a):
+    fA.neg()
+
 func `+=`*(a: var ExtensionField, b: ExtensionField) =
  ## Addition in the extension field
  for fA, fB in fields(a, b):
--- a/constantine/towers.nim
+++ b/constantine/towers.nim
@ -121,18 +121,22 @@ func `/=`*(a: var Fp2, _: typedesc[SexticNonResidue]) {.inline.} =
 # ----------------------------------------------------------------

 type
+  Fp4*[C: static Curve] = object
+    c0*, c1*: Fp2[C]
+
  Fp6*[C: static Curve] = object
    c0*, c1*, c2*: Fp2[C]

-  ξ = NonResidue
-    # We call the non-residue ξ on 𝔽p6 to avoid confusion between non-residue
+  ξ* = NonResidue
+    # We call the non-residue ξ on 𝔽p4/𝔽p6 to avoid confusion
+    # between non-residue
    # of different tower level

 func `*`*(_: typedesc[ξ], a: Fp2): Fp2 {.inline, noInit.} =
-  ## Multiply an element of 𝔽p2 by the cubic non-residue
-  ## chosen to construct 𝔽p6
+  ## Multiply an element of 𝔽p2 by the quadratic and cubic non-residue
+  ## chosen to construct 𝔽p4/𝔽p6
  # Yet another const tuple unpacking bug
-  const u = Fp2.C.get_CNR_Fp2()[0] # Cubic non-residue to construct 𝔽p6
+  const u = Fp2.C.get_CNR_Fp2()[0] # Quadratic & Cubic non-residue to construct 𝔽p4/𝔽p6
  const v = Fp2.C.get_CNR_Fp2()[1]
  const Beta = Fp2.C.get_QNR_Fp()  # Quadratic non-residue to construct 𝔽p2
  # ξ = u + v x
@ -176,21 +180,20 @@ func `*=`*(a: var Fp2, _: typedesc[ξ]) {.inline.} =

 type
  Fp12*[C: static Curve] = object
-    c0*, c1*: Fp6[C]
+    c0*, c1*, c2*: Fp4[C]

  γ = NonResidue
    # We call the non-residue γ (Gamma) on 𝔽p6 to avoid confusion between non-residue
    # of different tower level

-func `*`*(_: typedesc[γ], a: Fp6): Fp6 {.noInit, inline.} =
+func `*`*(_: typedesc[γ], a: Fp4): Fp4 {.noInit, inline.} =
  ## Multiply an element of 𝔽p6 by the cubic non-residue
  ## chosen to construct 𝔽p12
  ## For all curves γ = v with v the factor for 𝔽p6 coordinate
  ## and v³ = ξ
  ## (c0 + c1 v + c2 v²) v => ξ c2 + c0 v + c1 v²
-  result.c0 = ξ * a.c2
+  result.c0 = ξ * a.c1
  result.c1 = a.c0
-  result.c2 = a.c1

-func `*=`*(a: var Fp6, _: typedesc[γ]) {.inline.} =
+func `*=`*(a: var Fp4, _: typedesc[γ]) {.inline.} =
  a = γ * a
--- a/helpers/prng_unsafe.nim
+++ b/helpers/prng_unsafe.nim
@ -226,7 +226,7 @@ func random_long01Seq[T](rng: var RngState, a: var T, C: static Curve) =
 # Elliptic curves
 # ------------------------------------------------------------

-func random_unsafe[F](rng: var RngState, a: var ECP_SWei_Proj[F]) =
+func random_unsafe[F](rng: var RngState, a: var (ECP_SWei_Proj[F] or ECP_SWei_Aff[F])) =
  ## Initialize a random curve point with Z coordinate == 1
  ## Unsafe: for testing and benchmarking purposes only
  var fieldElem {.noInit.}: F
@ -251,7 +251,7 @@ func random_unsafe_with_randZ[F](rng: var RngState, a: var ECP_SWei_Proj[F]) =
    rng.random_unsafe(fieldElem, F.C)
    success = trySetFromCoordsXandZ(a, fieldElem, Z)

-func random_highHammingWeight[F](rng: var RngState, a: var ECP_SWei_Proj[F]) =
+func random_highHammingWeight[F](rng: var RngState, a: var (ECP_SWei_Proj[F] or ECP_SWei_Aff[F])) =
  ## Initialize a random curve point with Z coordinate == 1
  ## This will be generated with a biaised RNG with high Hamming Weight
  ## to trigger carry bugs
@ -264,7 +264,7 @@ func random_highHammingWeight[F](rng: var RngState, a: var ECP_SWei_Proj[F]) =
    rng.random_highHammingWeight(fieldElem, F.C)
    success = trySetFromCoordX(a, fieldElem)

-func random_highHammingWeight_with_randZ[F](rng: var RngState, a: var ECP_SWei_Proj[F]) =
+func random_highHammingWeight_with_randZ[F](rng: var RngState, a: var (ECP_SWei_Proj[F] or ECP_SWei_Aff[F])) =
  ## Initialize a random curve point with Z coordinate == 1
  ## This will be generated with a biaised RNG with high Hamming Weight
  ## to trigger carry bugs
@ -278,7 +278,7 @@ func random_highHammingWeight_with_randZ[F](rng: var RngState, a: var ECP_SWei_P
    rng.random_highHammingWeight(fieldElem, F.C)
    success = trySetFromCoordsXandZ(a, fieldElem, Z)

-func random_long01Seq[F](rng: var RngState, a: var ECP_SWei_Proj[F]) =
+func random_long01Seq[F](rng: var RngState, a: var (ECP_SWei_Proj[F] or ECP_SWei_Aff[F])) =
  ## Initialize a random curve point with Z coordinate == 1
  ## This will be generated with a biaised RNG
  ## that produces long bitstrings of 0 and 1
@ -313,7 +313,7 @@ func random_long01Seq_with_randZ[F](rng: var RngState, a: var ECP_SWei_Proj[F])
 func random_unsafe*(rng: var RngState, T: typedesc): T =
  ## Create a random Field or Extension Field or Curve Element
  ## Unsafe: for testing and benchmarking purposes only
-  when T is ECP_SWei_Proj:
+  when T is (ECP_SWei_Proj or ECP_SWei_Aff):
    rng.random_unsafe(result)
  elif T is SomeNumber:
    cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
@ -330,7 +330,7 @@ func random_unsafe_with_randZ*(rng: var RngState, T: typedesc[ECP_SWei_Proj]): T
 func random_highHammingWeight*(rng: var RngState, T: typedesc): T =
  ## Create a random Field or Extension Field or Curve Element
  ## Skewed towards high Hamming Weight
-  when T is ECP_SWei_Proj:
+  when T is (ECP_SWei_Proj or ECP_SWei_Aff):
    rng.random_highHammingWeight(result)
  elif T is SomeNumber:
    cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
@ -347,7 +347,7 @@ func random_highHammingWeight_with_randZ*(rng: var RngState, T: typedesc[ECP_SWe
 func random_long01Seq*(rng: var RngState, T: typedesc): T =
  ## Create a random Field or Extension Field or Curve Element
  ## Skewed towards long bitstrings of 0 or 1
-  when T is ECP_SWei_Proj:
+  when T is (ECP_SWei_Proj or ECP_SWei_Aff):
    rng.random_long01Seq(result)
  elif T is SomeNumber:
    cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
--- a/helpers/static_for.nim
+++ b/helpers/static_for.nim
@ -34,6 +34,14 @@ macro staticFor*(ident: untyped{nkIdent}, choices: typed, body: untyped): untype
  ##   staticFor(curve, TestCurves):
  ##     body
  ## and unroll the body for each curve in TestCurves
+
+  let choices = if choices.kind == nnkSym:
+                  # Unpack symbol
+                  choices.getImpl()
+                else:
+                  choices.expectKind(nnkBracket)
+                  choices
+
  result = newStmtList()
  for choice in choices:
    result.add nnkBlockStmt.newTree(
--- a/sage/frobenius_bls12_381.sage
+++ b/sage/frobenius_bls12_381.sage
@ -50,7 +50,15 @@ def fp2_to_hex(a):
    v = vector(a)
    return Integer(v[0]).hex() + ' + β * ' + Integer(v[1]).hex()

-# Frobenius constants (D type: use SNR, M type use 1/SNR)
+# Frobenius map constants (D type: use SNR, M type use 1/SNR)
+print('\nFrobenius extension field constants')
+FrobConst_map = SNR^((p-1)/6)
+print('FrobConst_map     : ' + fp2_to_hex(FrobConst_map))
+FrobConst_map_fp4 = FrobConst_map^(12//4)
+print('FrobConst_map_fp4_y : ' + fp2_to_hex(FrobConst_map_fp4))
+
+# Frobenius psi constants (D type: use SNR, M type use 1/SNR)
+print('\nψ (Psi) - Untwist-Frobenius-Twist constants')
 print('1/sextic_non_residue: ' + fp2_to_hex(1/SNR))
 FrobConst_psi = (1/SNR)^((p-1)/6)
 FrobConst_psi_2 = FrobConst_psi * FrobConst_psi
--- a/sage/frobenius_bn254_nogami.sage
+++ b/sage/frobenius_bn254_nogami.sage
@ -8,7 +8,7 @@

 # ############################################################
 #
-#                         BLS12-381
+#                         BN254-Nogami
 #                   Frobenius Endomorphism
 #              Untwist-Frobenius-Twist isogeny
 #
--- a/sage/frobenius_bn254_snarks.sage
+++ b/sage/frobenius_bn254_snarks.sage
@ -8,7 +8,7 @@

 # ############################################################
 #
-#                         BLS12-381
+#                         BN254-Snarks
 #                   Frobenius Endomorphism
 #              Untwist-Frobenius-Twist isogeny
 #
--- a/tests/t_bigints_mod_vs_gmp.nim
+++ b/tests/t_bigints_mod_vs_gmp.nim
@ -92,7 +92,7 @@ proc main() =
  mpz_init(m)
  mpz_init(r)

-  testRandomModSizes(24, aBits, mBits):
+  testRandomModSizes(12, aBits, mBits):
    # echo "--------------------------------------------------------------------------------"
    echo "Testing: random dividend (" & align($aBits, 4) & "-bit) -- random modulus (" & align($mBits, 4) & "-bit)"

--- a/tests/t_bigints_mul_high_words_vs_gmp.nim
+++ b/tests/t_bigints_mul_high_words_vs_gmp.nim
@ -67,7 +67,7 @@ proc main() =
  mpz_init(a)
  mpz_init(b)

-  testRandomModSizes(24, rBits, aBits, bBits, wordsStartIndex):
+  testRandomModSizes(12, rBits, aBits, bBits, wordsStartIndex):
    # echo "--------------------------------------------------------------------------------"
    echo "Testing: random mul_high_words  r (", align($rBits, 4),
      "-bit, keeping from ", wordsStartIndex,
--- a/tests/t_bigints_mul_vs_gmp.nim
+++ b/tests/t_bigints_mul_vs_gmp.nim
@ -64,7 +64,7 @@ proc main() =
  mpz_init(a)
  mpz_init(b)

-  testRandomModSizes(24, rBits, aBits, bBits):
+  testRandomModSizes(12, rBits, aBits, bBits):
    # echo "--------------------------------------------------------------------------------"
    echo "Testing: random mul  r (", align($rBits, 4), "-bit) <- a (", align($aBits, 4), "-bit) * b (", align($bBits, 4), "-bit) (full mul bits: ", align($(aBits+bBits), 4), "), r large enough? ", rBits >= aBits+bBits

--- a/tests/t_ec_template.nim
+++ b/tests/t_ec_template.nim
@ -31,21 +31,21 @@ type
    HighHammingWeight
    Long01Sequence

-func random_point*(rng: var RngState, F: typedesc, randZ: static bool, gen: static RandomGen): F {.inline, noInit.} =
-  when not randZ:
-    when gen == Uniform:
-      result = rng.random_unsafe(F)
+func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen): EC {.noInit.} =
+  if not randZ:
+    if gen == Uniform:
+      result = rng.random_unsafe(EC)
    elif gen == HighHammingWeight:
-      result = rng.random_highHammingWeight(F)
+      result = rng.random_highHammingWeight(EC)
    else:
-      result = rng.random_long01Seq(F)
+      result = rng.random_long01Seq(EC)
  else:
-    when gen == Uniform:
-      result = rng.random_unsafe_with_randZ(F)
+    if gen == Uniform:
+      result = rng.random_unsafe_with_randZ(EC)
    elif gen == HighHammingWeight:
-      result = rng.random_highHammingWeight_with_randZ(F)
+      result = rng.random_highHammingWeight_with_randZ(EC)
    else:
-      result = rng.random_long01Seq_with_randZ(F)
+      result = rng.random_long01Seq_with_randZ(EC)

 proc run_EC_addition_tests*(
       ec: typedesc,
@ -69,7 +69,7 @@ proc run_EC_addition_tests*(

  suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
    test "The infinity point is the neutral element w.r.t. to EC " & G1_or_G2 & " addition":
-      proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, randZ: bool, gen: RandomGen) =
        var inf {.noInit.}: EC
        inf.setInf()
        check: bool inf.isInf()
@ -92,7 +92,7 @@ proc run_EC_addition_tests*(
      test(ec, randZ = true, gen = Long01Sequence)

    test "Adding opposites gives an infinity point":
-      proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< Iters:
          var r{.noInit.}: EC
          let P = rng.random_point(EC, randZ, gen)
@ -113,7 +113,7 @@ proc run_EC_addition_tests*(
      test(ec, randZ = true, gen = Long01Sequence)

    test "EC " & G1_or_G2 & " add is commutative":
-      proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< Iters:
          var r0{.noInit.}, r1{.noInit.}: EC
          let P = rng.random_point(EC, randZ, gen)
@ -131,7 +131,7 @@ proc run_EC_addition_tests*(
      test(ec, randZ = true, gen = Long01Sequence)

    test "EC " & G1_or_G2 & " add is associative":
-      proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< Iters:
          let a = rng.random_point(EC, randZ, gen)
          let b = rng.random_point(EC, randZ, gen)
@ -180,7 +180,7 @@ proc run_EC_addition_tests*(
      test(ec, randZ = true, gen = Long01Sequence)

    test "EC " & G1_or_G2 & " double and EC " & G1_or_G2 & " add are consistent":
-      proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< Iters:
          let a = rng.random_point(EC, randZ, gen)

@ -220,7 +220,7 @@ proc run_EC_mul_sanity_tests*(

  suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
    test "EC " & G1_or_G2 & " mul [0]P == Inf":
-      proc test(EC: typedesc, bits: static int, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, bits: static int, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< ItersMul:
          let a = rng.random_point(EC, randZ, gen)

@ -243,7 +243,7 @@ proc run_EC_mul_sanity_tests*(
      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = Long01Sequence)

    test "EC " & G1_or_G2 & " mul [1]P == P":
-      proc test(EC: typedesc, bits: static int, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, bits: static int, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< ItersMul:
          let a = rng.random_point(EC, randZ, gen)

@ -269,7 +269,7 @@ proc run_EC_mul_sanity_tests*(
      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = Long01Sequence)

    test "EC " & G1_or_G2 & " mul [2]P == P.double()":
-      proc test(EC: typedesc, bits: static int, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, bits: static int, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< ItersMul:
          let a = rng.random_point(EC, randZ, gen)

@ -319,7 +319,7 @@ proc run_EC_mul_distributive_tests*(
  suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":

    test "EC " & G1_or_G2 & " mul is distributive over EC add":
-      proc test(EC: typedesc, bits: static int, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, bits: static int, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< ItersMul:
          let a = rng.random_point(EC, randZ, gen)
          let b = rng.random_point(EC, randZ, gen)
@ -388,7 +388,7 @@ proc run_EC_mul_vs_ref_impl*(

  suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
    test "EC " & G1_or_G2 & " mul constant-time is equivalent to a simple double-and-add algorithm":
-      proc test(EC: typedesc, bits: static int, randZ: static bool, gen: static RandomGen) =
+      proc test(EC: typedesc, bits: static int, randZ: bool, gen: RandomGen) =
        for _ in 0 ..< ItersMul:
          let a = rng.random_point(EC, randZ, gen)

--- a/tests/t_ec_wstrass_prj_g1_add_double.nim
+++ b/tests/t_ec_wstrass_prj_g1_add_double.nim
@ -19,7 +19,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 8

 run_EC_addition_tests(
    ec = ECP_SWei_Proj[Fp[BN254_Snarks]],
--- a/tests/t_ec_wstrass_prj_g1_mul_distri.nim
+++ b/tests/t_ec_wstrass_prj_g1_mul_distri.nim
@ -20,7 +20,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_distributive_tests(
--- a/tests/t_ec_wstrass_prj_g1_mul_sanity.nim
+++ b/tests/t_ec_wstrass_prj_g1_mul_sanity.nim
@ -20,7 +20,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_sanity_tests(
--- a/tests/t_ec_wstrass_prj_g1_mul_vs_ref.nim
+++ b/tests/t_ec_wstrass_prj_g1_mul_vs_ref.nim
@ -20,7 +20,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_vs_ref_impl(
--- a/tests/t_ec_wstrass_prj_g2_add_double_bls12_381.nim
+++ b/tests/t_ec_wstrass_prj_g2_add_double_bls12_381.nim
@ -20,7 +20,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 8

 run_EC_addition_tests(
    ec = ECP_SWei_Proj[Fp2[BLS12_381]],
--- a/tests/t_ec_wstrass_prj_g2_add_double_bn254_snarks.nim
+++ b/tests/t_ec_wstrass_prj_g2_add_double_bn254_snarks.nim
@ -20,7 +20,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 8

 run_EC_addition_tests(
    ec = ECP_SWei_Proj[Fp2[BN254_Snarks]],
--- a/tests/t_ec_wstrass_prj_g2_mul_distri_bls12_381.nim
+++ b/tests/t_ec_wstrass_prj_g2_mul_distri_bls12_381.nim
@ -21,7 +21,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_distributive_tests(
--- a/tests/t_ec_wstrass_prj_g2_mul_distri_bn254_snarks.nim
+++ b/tests/t_ec_wstrass_prj_g2_mul_distri_bn254_snarks.nim
@ -21,7 +21,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_distributive_tests(
--- a/tests/t_ec_wstrass_prj_g2_mul_sanity_bls12_381.nim
+++ b/tests/t_ec_wstrass_prj_g2_mul_sanity_bls12_381.nim
@ -21,7 +21,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_sanity_tests(
--- a/tests/t_ec_wstrass_prj_g2_mul_sanity_bn254_snarks.nim
+++ b/tests/t_ec_wstrass_prj_g2_mul_sanity_bn254_snarks.nim
@ -21,7 +21,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_sanity_tests(
--- a/tests/t_ec_wstrass_prj_g2_mul_vs_ref_bls12_381.nim
+++ b/tests/t_ec_wstrass_prj_g2_mul_vs_ref_bls12_381.nim
@ -21,7 +21,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_vs_ref_impl(
--- a/tests/t_ec_wstrass_prj_g2_mul_vs_ref_bn254_snarks.nim
+++ b/tests/t_ec_wstrass_prj_g2_mul_vs_ref_bn254_snarks.nim
@ -21,7 +21,7 @@ import
  ./t_ec_template

 const
-  Iters = 24
+  Iters = 12
  ItersMul = Iters div 4

 run_EC_mul_vs_ref_impl(
--- a/tests/t_finite_fields_powinv.nim
+++ b/tests/t_finite_fields_powinv.nim
@ -20,7 +20,7 @@ import

 static: doAssert defined(testingCurves), "This modules requires the -d:testingCurves compile option"

-const Iters = 24
+const Iters = 8

 var rng: RngState
 let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
--- a/tests/t_finite_fields_sqrt.nim
+++ b/tests/t_finite_fields_sqrt.nim
@ -17,7 +17,7 @@ import
  ../helpers/prng_unsafe


-const Iters = 24
+const Iters = 8

 var rng: RngState
 let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
--- a/tests/t_fp12_bls12_377.nim
+++ b/tests/t_fp12_bls12_377.nim
@ -19,7 +19,7 @@ const TestCurves = [

 runTowerTests(
  ExtDegree = 12,
-  Iters = 24,
+  Iters = 12,
  TestCurves = TestCurves,
  moduleName = "test_fp12_" & $BLS12_377,
  testSuiteDesc = "𝔽p12 = 𝔽p6[w] " & $BLS12_377
--- a/tests/t_fp12_bls12_381.nim
+++ b/tests/t_fp12_bls12_381.nim
@ -19,7 +19,7 @@ const TestCurves = [

 runTowerTests(
  ExtDegree = 12,
-  Iters = 24,
+  Iters = 12,
  TestCurves = TestCurves,
  moduleName = "test_fp12_" & $BLS12_381,
  testSuiteDesc = "𝔽p12 = 𝔽p6[w] " & $BLS12_381
--- a/tests/t_fp12_bn254_snarks.nim
+++ b/tests/t_fp12_bn254_snarks.nim
@ -19,7 +19,7 @@ const TestCurves = [

 runTowerTests(
  ExtDegree = 12,
-  Iters = 24,
+  Iters = 12,
  TestCurves = TestCurves,
  moduleName = "test_fp12_" & $BN254_Snarks,
  testSuiteDesc = "𝔽p12 = 𝔽p6[w] " & $BN254_Snarks
--- a/tests/t_fp12_exponentiation.nim
+++ b/tests/t_fp12_exponentiation.nim
@ -0,0 +1,219 @@
+# 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,
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+const
+  Iters = 2
+  TestCurves = [
+    BN254_Snarks,
+    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_fp12_exponentiation 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)
+
+proc test_sameBaseProduct(C: static Curve, gen: RandomGen) =
+  ## xᴬ xᴮ = xᴬ⁺ᴮ - product of power
+  let x = rng.random_elem(Fp12[C], gen)
+
+  var a = rng.random_elem(BigInt[128], gen)
+  var b = rng.random_elem(BigInt[128], gen)
+  # div by 2 to ensure their sum doesn't overflow
+  # 128 bits
+  a.div2()
+  b.div2()
+
+  var xa = x
+  xa.powUnsafeExponent(a, window = 3)
+
+  var xb = x
+  xb.powUnsafeExponent(b, window = 3)
+
+  var xapb = x
+  var apb: BigInt[128]
+  discard apb.sum(a, b)
+  xapb.powUnsafeExponent(apb, window = 3)
+
+  xa *= xb
+  check: bool(xa == xapb)
+
+proc test_powpow(C: static Curve, gen: RandomGen) =
+  ## (xᴬ)ᴮ = xᴬᴮ - power of power
+  var x = rng.random_elem(Fp12[C], gen)
+
+  var a = rng.random_elem(BigInt[128], gen)
+  var b = rng.random_elem(BigInt[128], gen)
+
+  var ab: BigInt[256]
+  ab.prod(a, b)
+
+  var y = x
+
+  x.powUnsafeExponent(a, window = 3)
+  x.powUnsafeExponent(b, window = 3)
+
+  y.powUnsafeExponent(ab, window = 3)
+  check: bool(x == y)
+
+proc test_powprod(C: static Curve, gen: RandomGen) =
+  ## (xy)ᴬ = xᴬyᴬ - power of product
+  var x = rng.random_elem(Fp12[C], gen)
+  var y = rng.random_elem(Fp12[C], gen)
+
+  let a = rng.random_elem(BigInt[128], gen)
+
+  var xy{.noInit.}: Fp12[C]
+  xy.prod(x, y)
+
+  xy.powUnsafeExponent(a, window=3)
+
+  x.powUnsafeExponent(a, window=3)
+  y.powUnsafeExponent(a, window=3)
+
+  x *= y
+
+  check: bool(x == xy)
+
+proc test_pow0(C: static Curve, gen: RandomGen) =
+  ## x⁰ = 1
+  var x = rng.random_elem(Fp12[C], gen)
+  var a: BigInt[128] # 0-init
+
+  x.powUnsafeExponent(a, window=3)
+  check: bool x.isOne()
+
+proc test_0pow0(C: static Curve, gen: RandomGen) =
+  ## 0⁰ = 1
+  var x: Fp12[C] # 0-init
+  var a: BigInt[128] # 0-init
+
+  x.powUnsafeExponent(a, window=3)
+  check: bool x.isOne()
+
+proc test_powinv(C: static Curve, gen: RandomGen) =
+  ## xᴬ / xᴮ = xᴬ⁻ᴮ - quotient of power
+  let x = rng.random_elem(Fp12[C], gen)
+
+  var a = rng.random_elem(BigInt[128], gen)
+  var b = rng.random_elem(BigInt[128], gen)
+  # div by 2 to ensure their sum doesn't overflow
+  # 128 bits
+  a.div2()
+  b.div2()
+  # Ensure a > b
+  cswap(a, b, a < b)
+
+  var xa = x
+  xa.powUnsafeExponent(a, window = 3)
+
+  var xb = x
+  xb.powUnsafeExponent(b, window = 3)
+
+  xb.inv()
+  xa *= xb
+
+  var xamb = x
+  var amb: BigInt[128]
+  discard amb.diff(a, b)
+  xamb.powUnsafeExponent(amb, window = 3)
+
+  check: bool(xa == xamb)
+
+proc test_invpow(C: static Curve, gen: RandomGen) =
+  ## (x / y)ᴬ = xᴬ / yᴬ - power of quotient
+  let x = rng.random_elem(Fp12[C], gen)
+  let y = rng.random_elem(Fp12[C], gen)
+
+  var a = rng.random_elem(BigInt[128], gen)
+
+  var xa = x
+  xa.powUnsafeExponent(a, window = 3)
+
+  var ya = y
+  ya.powUnsafeExponent(a, window = 3)
+  ya.inv()
+  xa *= ya
+
+  var xqya = x
+  var invy = y
+  invy.inv()
+  xqya *= invy
+  xqya.powUnsafeExponent(a, window = 3)
+
+  check: bool(xa == xqya)
+
+suite "Exponentiation in 𝔽p12" & " [" & $WordBitwidth & "-bit mode]":
+  staticFor(curve, TestCurves):
+    test "xᴬ xᴮ = xᴬ⁺ᴮ on " & $curve:
+      test_sameBaseProduct(curve, gen = Uniform)
+      test_sameBaseProduct(curve, gen = HighHammingWeight)
+      test_sameBaseProduct(curve, gen = Long01Sequence)
+
+  staticFor(curve, TestCurves):
+    test "(xᴬ)ᴮ = xᴬᴮ on " & $curve:
+      test_powpow(curve, gen = Uniform)
+      test_powpow(curve, gen = HighHammingWeight)
+      test_powpow(curve, gen = Long01Sequence)
+
+  staticFor(curve, TestCurves):
+    test "(xy)ᴬ = xᴬyᴬ on " & $curve:
+      test_powprod(curve, gen = Uniform)
+      test_powprod(curve, gen = HighHammingWeight)
+      test_powprod(curve, gen = Long01Sequence)
+
+  staticFor(curve, TestCurves):
+    test "x⁰ = 1 on " & $curve:
+      test_pow0(curve, gen = Uniform)
+      test_pow0(curve, gen = HighHammingWeight)
+      test_pow0(curve, gen = Long01Sequence)
+
+  staticFor(curve, TestCurves):
+    test "0⁰ = 1 on " & $curve:
+      test_0pow0(curve, gen = Uniform)
+      test_0pow0(curve, gen = HighHammingWeight)
+      test_0pow0(curve, gen = Long01Sequence)
+
+  staticFor(curve, TestCurves):
+    test "xᴬ / xᴮ = xᴬ⁻ᴮ on " & $curve:
+      test_powinv(curve, gen = Uniform)
+      test_powinv(curve, gen = HighHammingWeight)
+      test_powinv(curve, gen = Long01Sequence)
+
+  staticFor(curve, TestCurves):
+    test "(x / y)ᴬ = xᴬ / yᴬ on " & $curve:
+      test_invpow(curve, gen = Uniform)
+      test_invpow(curve, gen = HighHammingWeight)
+      test_invpow(curve, gen = Long01Sequence)
--- a/tests/t_fp2_frobenius.nim
+++ b/tests/t_fp2_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 = 2,
+  Iters = 8,
+  TestCurves = TestCurves,
+  moduleName = "test_fp2_frobenius",
+  testSuiteDesc = "𝔽p2 Frobenius map: Frobenius(a, k) = a^(p^k) (mod p²)"
+)
--- a/tests/t_fp2_sqrt.nim
+++ b/tests/t_fp2_sqrt.nim
@ -16,9 +16,20 @@ import
  ../constantine/config/curves,
  ../constantine/io/io_towers,
  # Test utilities
-  ../helpers/prng_unsafe
+  ../helpers/[prng_unsafe, static_for]

-const Iters = 24
+const
+  Iters = 8
+  TestCurves = [
+    BN254_Snarks,
+    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
@ -26,33 +37,43 @@ rng.seed(seed)
 echo "\n------------------------------------------------------\n"
 echo "test_fp2_sqrt xoshiro512** seed: ", seed

-proc randomSqrtCheck_p3mod4(C: static Curve) =
-  test "[𝔽p2] Random square root check for p ≡ 3 (mod 4) on " & $Curve(C):
-    for _ in 0 ..< Iters:
-      let a = rng.random_unsafe(Fp2[C])
-      var na{.noInit.}: Fp2[C]
-      na.neg(a)
+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)

-      var a2 = a
-      var na2 = na
-      a2.square()
-      na2.square()
-      check:
-        bool a2 == na2
-        bool a2.isSquare()
+proc randomSqrtCheck_p3mod4(C: static Curve, gen: RandomGen) =
+  for _ in 0 ..< Iters:
+    let a = rng.random_elem(Fp2[C], gen)
+    var na{.noInit.}: Fp2[C]
+    na.neg(a)

-      var r, s = a2
-      # r.sqrt()
-      let ok = s.sqrt_if_square()
-      check:
-        bool ok
-        # bool(r == s)
-        bool(s == a or s == na)
+    var a2 = a
+    var na2 = na
+    a2.square()
+    na2.square()
+    check:
+      bool a2 == na2
+      bool a2.isSquare()
+
+    var r, s = a2
+    # r.sqrt()
+    let ok = s.sqrt_if_square()
+    check:
+      bool ok
+      # bool(r == s)
+      bool(s == a or s == na)

 proc main() =
  suite "Modular square root" & " [" & $WordBitwidth & "-bit mode]":
-    randomSqrtCheck_p3mod4 BN254_Snarks
-    randomSqrtCheck_p3mod4 BLS12_381
+    staticFor(curve, TestCurves):
+      test "[𝔽p2] Random square root check for p ≡ 3 (mod 4) on " & $curve:
+        randomSqrtCheck_p3mod4(curve, gen = Uniform)
+        randomSqrtCheck_p3mod4(curve, gen = HighHammingWeight)
+        randomSqrtCheck_p3mod4(curve, gen = Long01Sequence)

  suite "Modular square root - 32-bit bugs highlighted by property-based testing " & " [" & $WordBitwidth & "-bit mode]":
    test "sqrt_if_square invalid square BLS12_381 - #64":
--- a/tests/t_fp4_frobenius.nim
+++ b/tests/t_fp4_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 = 4,
+  Iters = 8,
+  TestCurves = TestCurves,
+  moduleName = "test_fp4_frobenius",
+  testSuiteDesc = "𝔽p4 Frobenius map: Frobenius(a, k) = a^(p^k) (mod p⁴)"
+)
--- a/tests/t_fp6_bls12_377.nim
+++ b/tests/t_fp6_bls12_377.nim
@ -19,7 +19,7 @@ const TestCurves = [

 runTowerTests(
  ExtDegree = 6,
-  Iters = 24,
+  Iters = 12,
  TestCurves = TestCurves,
  moduleName = "test_fp6_" & $BLS12_377,
  testSuiteDesc = "𝔽p6 = 𝔽p2[v] " & $BLS12_377
--- a/tests/t_fp6_bls12_381.nim
+++ b/tests/t_fp6_bls12_381.nim
@ -19,7 +19,7 @@ const TestCurves = [

 runTowerTests(
  ExtDegree = 6,
-  Iters = 24,
+  Iters = 12,
  TestCurves = TestCurves,
  moduleName = "test_fp6_" & $BLS12_381,
  testSuiteDesc = "𝔽p6 = 𝔽p2[v] " & $BLS12_381
--- a/tests/t_fp6_bn254_snarks.nim
+++ b/tests/t_fp6_bn254_snarks.nim
@ -19,7 +19,7 @@ const TestCurves = [

 runTowerTests(
  ExtDegree = 6,
-  Iters = 128,
+  Iters = 12,
  TestCurves = TestCurves,
  moduleName = "test_fp6_" & $BN254_Snarks,
  testSuiteDesc = "𝔽p6 = 𝔽p2[w] " & $BN254_Snarks
--- a/tests/t_fp_tower_frobenius_template.nim
+++ b/tests/t_fp_tower_frobenius_template.nim
@ -0,0 +1,115 @@
+# 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.
+
+# ############################################################
+#
+#         Template tests for towered extension fields
+#
+# ############################################################
+
+
+import
+  # Standard library
+  std/[unittest, times],
+  # Internals
+  ../constantine/towers,
+  ../constantine/config/[common, curves],
+  ../constantine/arithmetic,
+  ../constantine/isogeny/frobenius,
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+echo "\n------------------------------------------------------\n"
+
+template ExtField(degree: static int, curve: static Curve): untyped =
+  when degree == 2:
+    Fp2[curve]
+  elif degree == 4:
+    Fp4[curve]
+  elif degree == 6:
+    Fp6[curve]
+  elif degree == 12:
+    Fp12[curve]
+  else:
+    {.error: "Unconfigured extension degree".}
+
+type
+  RandomGen = enum
+    Uniform
+    HighHammingWeight
+    Long01Sequence
+
+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)
+
+proc runFrobeniusTowerTests*[N](
+      ExtDegree: static int,
+      Iters: static int,
+      TestCurves: static array[N, Curve],
+      moduleName: string,
+      testSuiteDesc: string
+    ) =
+  # Random seed for reproducibility
+  var rng: RngState
+  let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+  rng.seed(seed)
+  echo moduleName, " xoshiro512** seed: ", seed
+
+  suite testSuiteDesc & " [" & $WordBitwidth & "-bit mode]":
+    test "Frobenius(a) = a^p (mod p^" & $ExtDegree & ")":
+      proc test(Field: typedesc, Iters: static int, gen: RandomGen) =
+        for _ in 0 ..< Iters:
+          var a = rng.random_elem(Field, gen)
+          var fa {.noInit.}: typeof(a)
+          fa.frobenius_map(a, k = 1)
+
+          a.powUnsafeExponent(Field.C.Mod, window = 3)
+          check: bool(a == fa)
+
+      staticFor(curve, TestCurves):
+        test(ExtField(ExtDegree, curve), Iters, gen = Uniform)
+        test(ExtField(ExtDegree, curve), Iters, gen = HighHammingWeight)
+        test(ExtField(ExtDegree, curve), Iters, gen = Long01Sequence)
+
+    # test "Frobenius(a, 2) = a^(p^2) (mod p^" & $ExtDegree & ")":
+    #   proc test(Field: typedesc, Iters: static int, gen: RandomGen) =
+    #     for _ in 0 ..< Iters:
+    #       var a = rng.random_elem(Field, gen)
+    #       var fa {.noInit.}: typeof(a)
+    #       fa.frobenius_map(a, k = 2)
+
+    #       a.powUnsafeExponent(Field.C.Mod, window = 3)
+    #       a.powUnsafeExponent(Field.C.Mod, window = 3)
+    #       check: bool(a == fa)
+
+    #   staticFor(curve, TestCurves):
+    #     test(ExtField(ExtDegree, curve), Iters, gen = Uniform)
+    #     test(ExtField(ExtDegree, curve), Iters, gen = HighHammingWeight)
+    #     test(ExtField(ExtDegree, curve), Iters, gen = Long01Sequence)
+
+    # test "Frobenius(a, 3) = a^(p^3) (mod p^" & $ExtDegree & ")":
+    #   proc test(Field: typedesc, Iters: static int, gen: RandomGen) =
+    #     for _ in 0 ..< Iters:
+    #       var a = rng.random_elem(Field, gen)
+    #       var fa {.noInit.}: typeof(a)
+    #       fa.frobenius_map(a, k = 3)
+
+    #       a.powUnsafeExponent(Field.C.Mod, window = 3)
+    #       a.powUnsafeExponent(Field.C.Mod, window = 3)
+    #       a.powUnsafeExponent(Field.C.Mod, window = 3)
+    #       check: bool(a == fa)
+
+    #   staticFor(curve, TestCurves):
+    #     test(ExtField(ExtDegree, curve), Iters, gen = Uniform)
+    #     test(ExtField(ExtDegree, curve), Iters, gen = HighHammingWeight)
+    #     test(ExtField(ExtDegree, curve), Iters, gen = Long01Sequence)
--- a/tests/t_fp_tower_template.nim
+++ b/tests/t_fp_tower_template.nim
@ -41,8 +41,8 @@ type
    HighHammingWeight
    Long01Sequence

-func random_elem(rng: var RngState, F: typedesc, gen: static RandomGen): F {.inline, noInit.} =
-  when gen == Uniform:
+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)
@ -77,7 +77,7 @@ proc runTowerTests*[N](
        test(ExtField(ExtDegree, curve))

    test "Addition, substraction negation are consistent":
-      proc test(Field: typedesc, Iters: static int, gen: static RandomGen) =
+      proc test(Field: typedesc, Iters: static int, gen: RandomGen) =
        # Try to exercise all code paths for in-place/out-of-place add/sum/sub/diff/double/neg
        # (1 - (-a) - b + (-a) - 2a) + (2a + 2b + (-b))  == 1
        var accum {.noInit.}, One {.noInit.}, a{.noInit.}, na{.noInit.}, b{.noInit.}, nb{.noInit.}, a2 {.noInit.}, b2 {.noInit.}: Field
--- a/tests/t_pairing_bls12_381_line_functions.nim
+++ b/tests/t_pairing_bls12_381_line_functions.nim
@ -0,0 +1,114 @@
+# 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/elliptic/[
+    ec_weierstrass_affine,
+    ec_weierstrass_projective,
+    ec_scalar_mul],
+  ../constantine/pairing/lines_projective,
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+const
+  Iters = 4
+  TestCurves = [
+    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_bls12_381_line_functions xoshiro512** seed: ", seed
+
+func random_point*(rng: var RngState, EC: typedesc, gen: RandomGen): EC {.noInit.} =
+  if gen == Uniform:
+    result = rng.random_unsafe(EC)
+  elif gen == HighHammingWeight:
+    result = rng.random_highHammingWeight(EC)
+  else:
+    result = rng.random_long01Seq(EC)
+
+func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen): EC {.noInit.} =
+  if not randZ:
+    if gen == Uniform:
+      result = rng.random_unsafe(EC)
+    elif gen == HighHammingWeight:
+      result = rng.random_highHammingWeight(EC)
+    else:
+      result = rng.random_long01Seq(EC)
+  else:
+    if gen == Uniform:
+      result = rng.random_unsafe_with_randZ(EC)
+    elif gen == HighHammingWeight:
+      result = rng.random_highHammingWeight_with_randZ(EC)
+    else:
+      result = rng.random_long01Seq_with_randZ(EC)
+
+suite "Pairing - Line Functions on BLS12-381":
+  test "Line double - lt,t(P)":
+    proc test_line_double(C: static Curve, randZ: bool, gen: RandomGen) =
+      for _ in 0 ..< Iters:
+        let P = rng.random_point(ECP_SWei_Aff[Fp[C]], gen)
+        var T = rng.random_point(ECP_SWei_Proj[Fp2[C]], randZ, gen)
+        let Q = rng.random_point(ECP_SWei_Proj[Fp2[C]], randZ, gen)
+        var l: Line[Fp2[C], C.getSexticTwist()]
+
+        var T2: typeof(Q)
+        T2.double(T)
+        l.line_double(T, P)
+
+        doAssert: bool(T == T2)
+
+    staticFor(curve, TestCurves):
+      test_line_double(curve, randZ = false, gen = Uniform)
+      test_line_double(curve, randZ = true, gen = Uniform)
+      test_line_double(curve, randZ = false, gen = HighHammingWeight)
+      test_line_double(curve, randZ = true, gen = HighHammingWeight)
+      test_line_double(curve, randZ = false, gen = Long01Sequence)
+      test_line_double(curve, randZ = true, gen = Long01Sequence)
+
+  test "Line add - lt,q(P)":
+    proc test_line_add(C: static Curve, randZ: bool, gen: RandomGen) =
+      for _ in 0 ..< Iters:
+        let P = rng.random_point(ECP_SWei_Aff[Fp[C]], gen)
+        let Q = rng.random_point(ECP_SWei_Proj[Fp2[C]], randZ, gen)
+        var T = rng.random_point(ECP_SWei_Proj[Fp2[C]], randZ, gen)
+        var l: Line[Fp2[C], C.getSexticTwist()]
+
+        var TQ{.noInit.}: typeof(T)
+        TQ.sum(T, Q)
+
+        var Qaff{.noInit.}: ECP_SWei_Aff[Fp2[C]]
+        Qaff.affineFromProjective(Q)
+        l.line_add(T, Qaff, P)
+
+        doAssert: bool(T == TQ)
+
+    staticFor(curve, TestCurves):
+      test_line_add(curve, randZ = false, gen = Uniform)
+      test_line_add(curve, randZ = true, gen = Uniform)
+      test_line_add(curve, randZ = false, gen = HighHammingWeight)
+      test_line_add(curve, randZ = true, gen = HighHammingWeight)
+      test_line_add(curve, randZ = false, gen = Long01Sequence)
+      test_line_add(curve, randZ = true, gen = Long01Sequence)
--- a/tests/t_pairing_bls12_381_optate.nim
+++ b/tests/t_pairing_bls12_381_optate.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
+  ../constantine/config/curves,
+  ../constantine/pairing/pairing_bls12,
+  # Test utilities
+  ./t_pairing_template
+
+runPairingTests(4, BLS12_381, pairing_bls12_reference)
--- a/tests/t_pairing_bn254_nogami_optate.nim
+++ b/tests/t_pairing_bn254_nogami_optate.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
+  ../constantine/config/curves,
+  ../constantine/pairing/pairing_bn,
+  # Test utilities
+  ./t_pairing_template
+
+runPairingTests(4, BN254_Nogami, pairing_bn_reference)
--- a/tests/t_pairing_bn254_snarks_optate.nim
+++ b/tests/t_pairing_bn254_snarks_optate.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
+  ../constantine/config/curves,
+  ../constantine/pairing/pairing_bn,
+  # Test utilities
+  ./t_pairing_template
+
+runPairingTests(4, BN254_Snarks, pairing_bn_reference)
--- a/tests/t_pairing_fp12_sparse.nim
+++ b/tests/t_pairing_fp12_sparse.nim
@ -0,0 +1,143 @@
+# 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/[
+    lines_projective,
+    gt_fp12
+  ],
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+const
+  Iters = 8
+  TestCurves = [
+    # BN254_Snarks,
+    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 - Sparse 𝔽p12 multiplication by line function is consistent with dense 𝔽p12 mul":
+  test "Dense 𝔽p6 by Sparse 0y0":
+    proc test_fp6_0y0(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        let a = rng.random_elem(Fp6[C], gen)
+        let y = rng.random_elem(Fp2[C], gen)
+        let b = Fp6[C](c1: y)
+
+        var r {.noInit.}, r2 {.noInit.}: Fp6[C]
+
+        r.prod(a, b)
+        r2.mul_sparse_by_0y0(a, y)
+
+        check: bool(r == r2)
+
+    staticFor(curve, TestCurves):
+      test_fp6_0y0(curve, gen = Uniform)
+      test_fp6_0y0(curve, gen = HighHammingWeight)
+      test_fp6_0y0(curve, gen = Long01Sequence)
+
+  test "Dense 𝔽p6 by Sparse xy0":
+    proc test_fp6_0y0(C: static Curve, gen: static RandomGen) =
+      for _ in 0 ..< Iters:
+        let a = rng.random_elem(Fp6[C], gen)
+        let x = rng.random_elem(Fp2[C], gen)
+        let y = rng.random_elem(Fp2[C], gen)
+        let b = Fp6[C](c0: x, c1: y)
+        let line = Line[Fp2[C], M_twist](x: x, y: y)
+
+        var r {.noInit.}, r2 {.noInit.}: Fp6[C]
+
+        r.prod(a, b)
+        r2.mul_by_line_xy0(a, line)
+
+        check: bool(r == r2)
+
+    staticFor(curve, TestCurves):
+      test_fp6_0y0(curve, gen = Uniform)
+      test_fp6_0y0(curve, gen = HighHammingWeight)
+      test_fp6_0y0(curve, gen = Long01Sequence)
+
+  when Fp12[BN254_Snarks]().c0.typeof is Fp6:
+    test "Sparse 𝔽p12 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)
+          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: Fp6[C](c0: x, c1: y),
+            c1: Fp6[C](c1: z)
+          )
+
+          a *= b
+          a2.mul_sparse_by_line_xy00z0(line)
+
+          check: bool(a == a2)
+
+      staticFor(curve, TestCurves):
+        test_fp12_xy00z0(curve, gen = Uniform)
+        test_fp12_xy00z0(curve, gen = HighHammingWeight)
+        test_fp12_xy00z0(curve, gen = Long01Sequence)
+
+    test "Sparse 𝔽p12 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)
+          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: Fp6[C](c0: x, c1: y, c2: z)
+          )
+
+          a *= b
+          a2.mul_sparse_by_line_xyz000(line)
+
+          check: bool(a == a2)
+
+      staticFor(curve, TestCurves):
+        test_fp12_xyz000(curve, gen = Uniform)
+        test_fp12_xyz000(curve, gen = HighHammingWeight)
+        test_fp12_xyz000(curve, gen = Long01Sequence)
--- a/tests/t_pairing_template.nim
+++ b/tests/t_pairing_template.nim
@ -0,0 +1,92 @@
+# 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/unittest, times,
+  # Internals
+  ../constantine/[arithmetic, primitives],
+  ../constantine/towers,
+  ../constantine/config/curves,
+  ../constantine/elliptic/ec_weierstrass_projective,
+  ../constantine/hash_to_curve/cofactors,
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+export
+  prng_unsafe, times, unittest,
+  ec_weierstrass_projective, arithmetic, towers,
+  primitives
+
+type
+  RandomGen* = enum
+    Uniform
+    HighHammingWeight
+    Long01Sequence
+
+func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen): EC {.noInit.} =
+  if not randZ:
+    if gen == Uniform:
+      result = rng.random_unsafe(EC)
+      result.clearCofactorReference()
+    elif gen == HighHammingWeight:
+      result = rng.random_highHammingWeight(EC)
+      result.clearCofactorReference()
+    else:
+      result = rng.random_long01Seq(EC)
+      result.clearCofactorReference()
+  else:
+    if gen == Uniform:
+      result = rng.random_unsafe_with_randZ(EC)
+      result.clearCofactorReference()
+    elif gen == HighHammingWeight:
+      result = rng.random_highHammingWeight_with_randZ(EC)
+      result.clearCofactorReference()
+    else:
+      result = rng.random_long01Seq_with_randZ(EC)
+      result.clearCofactorReference()
+
+template runPairingTests*(Iters: static int, C: static Curve, pairing_fn: untyped): 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_",$C,"_optate xoshiro512** seed: ", timeseed
+
+  proc test_bilinearity_double_impl(randZ: bool, gen: RandomGen) =
+    for _ in 0 ..< Iters:
+      let P = rng.random_point(ECP_SWei_Proj[Fp[C]], randZ, gen)
+      let Q = rng.random_point(ECP_SWei_Proj[Fp2[C]], randZ, gen)
+      var P2: typeof(P)
+      var Q2: typeof(Q)
+
+      var r {.noInit.}, r2 {.noInit.}, r3 {.noInit.}: Fp12[C]
+
+      P2.double(P)
+      Q2.double(Q)
+
+      r.pairing_fn(P, Q)
+      r.square()
+      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)
+
+  suite "Pairing - Optimal Ate on " & $C:
+    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)
+      test_bilinearity_double_impl(randZ = false, gen = HighHammingWeight)
+      test_bilinearity_double_impl(randZ = true, gen = HighHammingWeight)
+      test_bilinearity_double_impl(randZ = false, gen = Long01Sequence)
+      test_bilinearity_double_impl(randZ = true, gen = Long01Sequence)