Pairings for BN254-Nogami and BN254-Snarks (#86)

* Implement optimized final exponentiation for BN254-Nogami * And BN254 Snarks support * Optimize D-Twist sparse Fp12 x line multiplication * Move quadruple/octuple and add to Github issues: https://github.com/mratsim/constantine/issues/88 [skip ci]
2020-09-25 21:58:20 +02:00 · 2020-09-25 21:58:20 +02:00 · 03ecb31c57
parent f78ed23dad
commit 03ecb31c57
22 changed files with 708 additions and 75 deletions
--- a/benchmarks/bench_pairing_bls12_381.nim
+++ b/benchmarks/bench_pairing_bls12_381.nim
@ -36,7 +36,7 @@ proc main() =
    const curve = AvailableCurves[i]
    lineDoubleBench(curve, Iters)
    lineAddBench(curve, Iters)
-    mulFp12byLineBench(curve, Iters)
+    mulFp12byLine_xy000z_Bench(curve, Iters)
    separator()
    finalExpEasyBench(curve, Iters)
    finalExpHardBLS12Bench(curve, Iters)
--- a/benchmarks/bench_pairing_bn254_nogami.nim
+++ b/benchmarks/bench_pairing_bn254_nogami.nim
@ -0,0 +1,51 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  # Internals
+  ../constantine/config/curves,
+  ../constantine/arithmetic,
+  ../constantine/towers,
+  # Helpers
+  ../helpers/static_for,
+  ./bench_pairing_template,
+  # Standard library
+  std/strutils
+
+# ############################################################
+#
+#               Benchmark of pairings
+#                 for BN254-Nogami
+#
+# ############################################################
+
+
+const Iters = 50
+const AvailableCurves = [
+  BN254_Nogami,
+]
+
+proc main() =
+  separator()
+  staticFor i, 0, AvailableCurves.len:
+    const curve = AvailableCurves[i]
+    lineDoubleBench(curve, Iters)
+    lineAddBench(curve, Iters)
+    mulFp12byLine_xyz000_Bench(curve, Iters)
+    separator()
+    finalExpEasyBench(curve, Iters)
+    finalExpHardBNBench(curve, Iters)
+    separator()
+    millerLoopBNBench(curve, Iters)
+    finalExpBNBench(curve, Iters)
+    separator()
+    pairingBNBench(curve, Iters)
+    separator()
+
+main()
+notes()
--- a/benchmarks/bench_pairing_bn254_snarks.nim
+++ b/benchmarks/bench_pairing_bn254_snarks.nim
@ -0,0 +1,51 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  # Internals
+  ../constantine/config/curves,
+  ../constantine/arithmetic,
+  ../constantine/towers,
+  # Helpers
+  ../helpers/static_for,
+  ./bench_pairing_template,
+  # Standard library
+  std/strutils
+
+# ############################################################
+#
+#               Benchmark of pairings
+#                 for BN254-Snarks
+#
+# ############################################################
+
+
+const Iters = 50
+const AvailableCurves = [
+  BN254_Snarks,
+]
+
+proc main() =
+  separator()
+  staticFor i, 0, AvailableCurves.len:
+    const curve = AvailableCurves[i]
+    lineDoubleBench(curve, Iters)
+    lineAddBench(curve, Iters)
+    mulFp12byLine_xyz000_Bench(curve, Iters)
+    separator()
+    finalExpEasyBench(curve, Iters)
+    finalExpHardBNBench(curve, Iters)
+    separator()
+    millerLoopBNBench(curve, Iters)
+    finalExpBNBench(curve, Iters)
+    separator()
+    pairingBNBench(curve, Iters)
+    separator()
+
+main()
+notes()
--- a/benchmarks/bench_pairing_template.nim
+++ b/benchmarks/bench_pairing_template.nim
@ -24,7 +24,8 @@ import
    cyclotomic_fp12,
    lines_projective,
    mul_fp12_by_lines,
-    pairing_bls12
+    pairing_bls12,
+    pairing_bn
  ],
  # Helpers
  ../helpers/[prng_unsafe, static_for],
@ -150,7 +151,20 @@ proc lineAddBench*(C: static Curve, iters: int) =
  bench("Line add", C, iters):
    line.line_add(T, Qaff, Paff)

-proc mulFp12byLineBench*(C: static Curve, iters: int) =
+proc mulFp12byLine_xyz000_Bench*(C: static Curve, iters: int) =
+  var line: Line[Fp2[C], C.getSexticTwist()]
+  var T = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  let P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+  var Paff: ECP_SWei_Aff[Fp[C]]
+  Paff.affineFromProjective(P)
+
+  line.line_double(T, Paff)
+  var f = rng.random_unsafe(Fp12[C])
+
+  bench("Mul 𝔽p12 by line xyz000", C, iters):
+    f.mul_sparse_by_line_xyz000(line)
+
+proc mulFp12byLine_xy000z_Bench*(C: static Curve, iters: int) =
  var line: Line[Fp2[C], C.getSexticTwist()]
  var T = rng.random_point(ECP_SWei_Proj[Fp2[C]])
  let P = rng.random_point(ECP_SWei_Proj[Fp[C]])
@ -178,6 +192,21 @@ proc millerLoopBLS12Bench*(C: static Curve, iters: int) =
  bench("Miller Loop BLS12", C, iters):
    f.millerLoopGenericBLS12(Paff, Qaff)

+proc millerLoopBNBench*(C: static Curve, iters: int) =
+  let
+    P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+    Q = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+  var
+    Paff: ECP_SWei_Aff[Fp[C]]
+    Qaff: ECP_SWei_Aff[Fp2[C]]
+  Paff.affineFromProjective(P)
+  Qaff.affineFromProjective(Q)
+
+  var f: Fp12[C]
+
+  bench("Miller Loop BN", C, iters):
+    f.millerLoopGenericBN(Paff, Qaff)
+
 proc finalExpEasyBench*(C: static Curve, iters: int) =
  var r = rng.random_unsafe(Fp12[C])
  bench("Final Exponentiation Easy", C, iters):
@ -189,12 +218,24 @@ proc finalExpHardBLS12Bench*(C: static Curve, iters: int) =
  bench("Final Exponentiation Hard BLS12", C, iters):
    r.finalExpHard_BLS12()

+proc finalExpHardBNBench*(C: static Curve, iters: int) =
+  var r = rng.random_unsafe(Fp12[C])
+  r.finalExpEasy()
+  bench("Final Exponentiation Hard BN", C, iters):
+    r.finalExpHard_BN()
+
 proc finalExpBLS12Bench*(C: static Curve, iters: int) =
  var r = rng.random_unsafe(Fp12[C])
  bench("Final Exponentiation BLS12", C, iters):
    r.finalExpEasy()
    r.finalExpHard_BLS12()

+proc finalExpBNBench*(C: static Curve, iters: int) =
+  var r = rng.random_unsafe(Fp12[C])
+  bench("Final Exponentiation BN", C, iters):
+    r.finalExpEasy()
+    r.finalExpHard_BN()
+
 proc pairingBLS12Bench*(C: static Curve, iters: int) =
  let
    P = rng.random_point(ECP_SWei_Proj[Fp[C]])
@ -204,3 +245,13 @@ proc pairingBLS12Bench*(C: static Curve, iters: int) =

  bench("Pairing BLS12", C, iters):
    f.pairing_bls12(P, Q)
+
+proc pairingBNBench*(C: static Curve, iters: int) =
+  let
+    P = rng.random_point(ECP_SWei_Proj[Fp[C]])
+    Q = rng.random_point(ECP_SWei_Proj[Fp2[C]])
+
+  var f: Fp12[C]
+
+  bench("Pairing BN", C, iters):
+    f.pairing_bn(P, Q)
--- a/constantine.nimble
+++ b/constantine.nimble
@ -451,14 +451,44 @@ task bench_ec_g2_clang_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Sho
 task bench_pairing_bls12_381, "Run pairings benchmarks for BLS12-381 - Default compiler":
  runBench("bench_pairing_bls12_381")

-task bench_pairing_bls12_381_gcc, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - GCC":
+task bench_pairing_bls12_381_gcc, "Run pairings benchmarks for BLS12-381 - GCC":
  runBench("bench_pairing_bls12_381", "gcc")

-task bench_pairing_bls12_381_clang, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - Clang":
+task bench_pairing_bls12_381_clang, "Run pairings benchmarks for BLS12-381 - Clang":
  runBench("bench_pairing_bls12_381", "clang")

-task bench_pairing_bls12_381_gcc_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - GCC no Assembly":
+task bench_pairing_bls12_381_gcc_noasm, "Run pairings benchmarks for BLS12-381 - GCC no Assembly":
  runBench("bench_pairing_bls12_381", "gcc", useAsm = false)

-task bench_pairing_bls12_381_clang_noasm, "Run benchmark on Elliptic Curve group 𝔾2 - Short Weierstrass with Projective Coordinates - Clang no Assembly":
+task bench_pairing_bls12_381_clang_noasm, "Run pairings benchmarks for BLS12-381 - Clang no Assembly":
  runBench("bench_pairing_bls12_381", "clang", useAsm = false)
+
+task bench_pairing_bn254_nogami, "Run pairings benchmarks for BN254-Nogami - Default compiler":
+  runBench("bench_pairing_bn254_nogami")
+
+task bench_pairing_bn254_nogami_gcc, "Run pairings benchmarks for BN254-Nogami - GCC":
+  runBench("bench_pairing_bn254_nogami", "gcc")
+
+task bench_pairing_bn254_nogami_clang, "Run pairings benchmarks for BN254-Nogami - Clang":
+  runBench("bench_pairing_bn254_nogami", "clang")
+
+task bench_pairing_bn254_nogami_gcc_noasm, "Run pairings benchmarks for BN254-Nogami - GCC no Assembly":
+  runBench("bench_pairing_bn254_nogami", "gcc", useAsm = false)
+
+task bench_pairing_bn254_nogami_clang_noasm, "Run pairings benchmarks for BN254-Nogami - Clang no Assembly":
+  runBench("bench_pairing_bn254_nogami", "clang", useAsm = false)
+
+task bench_pairing_bn254_snarks, "Run pairings benchmarks for BN254-Snarks - Default compiler":
+  runBench("bench_pairing_bn254_snarks")
+
+task bench_pairing_bn254_snarks_gcc, "Run pairings benchmarks for BN254-Snarks - GCC":
+  runBench("bench_pairing_bn254_snarks", "gcc")
+
+task bench_pairing_bn254_snarks_clang, "Run pairings benchmarks for BN254-Snarks - Clang":
+  runBench("bench_pairing_bn254_snarks", "clang")
+
+task bench_pairing_bn254_snarks_gcc_noasm, "Run pairings benchmarks for BN254-Snarks - GCC no Assembly":
+  runBench("bench_pairing_bn254_snarks", "gcc", useAsm = false)
+
+task bench_pairing_bn254_snarks_clang_noasm, "Run pairings benchmarks for BN254-Snarks - Clang no Assembly":
+  runBench("bench_pairing_bn254_snarks", "clang", useAsm = false)
--- a/constantine/config/curves_declaration.nim
+++ b/constantine/config/curves_declaration.nim
@ -103,7 +103,7 @@ declareCurves:
    modulus: "0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47"
    family: BarretoNaehrig
    bn_u_bitwidth: 63
-    bn_u: "0x44E992B44A6909F1" # u: 4965661367192848881
+    bn_u: "0x44e992b44a6909f1" # u: 4965661367192848881
    cubicRootOfUnity_modP: "0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd48"
    # For sanity checks
    cubicRootOfUnity_modR: "0x30644e72e131a029048b6e193fd84104cc37a73fec2bc5e9b8ca0b2d36636f23"
--- a/constantine/isogeny/frobenius.nim
+++ b/constantine/isogeny/frobenius.nim
@ -68,7 +68,7 @@ template mulCheckSparse[Fp2](a: var Fp2, b: Fp2) =

 # c = (SNR^((p-1)/6)^coef).
 # Then for frobenius(2): c * conjugate(c)
-# And for frobenius(2): c² * conjugate(c)
+# And for frobenius(3): c² * conjugate(c)
 const FrobMapConst_BLS12_381 = [
  # frobenius(1)
  [Fp2[BLS12_381].fromHex( # SNR^((p-1)/6)^0
@ -146,23 +146,182 @@ const FrobMapConst_BLS12_381 = [
    "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2"
  )]]

-func frobenius_map*(r: var Fp4, a: Fp4, k: static int = 1) {.inline.} =
+const FrobMapConst_BN254_Nogami = [
+  # frobenius(1)
+  [Fp2[BN254_Nogami].fromHex( # SNR^((p-1)/6)^0
+    "0x1",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(  # SNR^((p-1)/6)^1
+    "0x1b377619212e7c8cb6499b50a846953f850974924d3f77c2e17de6c06f2a6de9",
+    "0x9ebee691ed1837503eab22f57b96ac8dc178b6db2c08850c582193f90d5922a"
+  ),
+  Fp2[BN254_Nogami].fromHex(  # SNR^((p-1)/6)^2 = SNR^((p-1)/3)
+    "0x0",
+    "0x25236482400000017080eb4000000006181800000000000cd98000000000000b"
+  ),
+  Fp2[BN254_Nogami].fromHex( # SNR^((p-1)/6)^3 = SNR^((p-1)/2)
+    "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5",
+    "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5"
+  ),
+  Fp2[BN254_Nogami].fromHex( # SNR^((p-1)/6)^4 = SNR^(2(p-1)/3)
+    "0x25236482400000017080eb4000000006181800000000000cd98000000000000c",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex( # SNR^((p-1)/6)^5
+    "0x19f3db6884cdca43c2b0d5792cd135accb1baea0b017046e859975ab54b5ef9b",
+    "0xb2f8919bb3235bdf7837806d32eca5b9605515f4fe8fba521668a54ab4a1078"
+  )],
+  # frobenius(2)
+  [Fp2[BN254_Nogami].fromHex( # norm(SNR)^((p-1)/6)^1
+    "0x1",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex( # norm(SNR)^((p-1)/6)^2
+    "0x49b36240000000024909000000000006cd80000000000008",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x49b36240000000024909000000000006cd80000000000007",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x2523648240000001ba344d80000000086121000000000013a700000000000012",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x25236482400000017080eb4000000006181800000000000cd98000000000000b",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x25236482400000017080eb4000000006181800000000000cd98000000000000c",
+    "0x0"
+  )],
+  # frobenius(3)
+  [Fp2[BN254_Nogami].fromHex(
+    "0x1",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x1439ab09c60b248f398c5d77b755f92b9edc5f19d2873545be471151a747e4e",
+    "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x0",
+    "0x1"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x1439ab09c60b248f398c5d77b755f92b9edc5f19d2873545be471151a747e4e",
+    "0x1439ab09c60b248f398c5d77b755f92b9edc5f19d2873545be471151a747e4e"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x2523648240000001ba344d80000000086121000000000013a700000000000012",
+    "0x0"
+  ),
+  Fp2[BN254_Nogami].fromHex(
+    "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5",
+    "0x1439ab09c60b248f398c5d77b755f92b9edc5f19d2873545be471151a747e4e"
+  )]]
+
+const FrobMapConst_BN254_Snarks = [
+  # frobenius(1)
+  [Fp2[BN254_Snarks].fromHex( # SNR^((p-1)/6)^0
+    "0x1",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex(  # SNR^((p-1)/6)^1
+    "0x1284b71c2865a7dfe8b99fdd76e68b605c521e08292f2176d60b35dadcc9e470",
+    "0x246996f3b4fae7e6a6327cfe12150b8e747992778eeec7e5ca5cf05f80f362ac"
+  ),
+  Fp2[BN254_Snarks].fromHex(  # SNR^((p-1)/6)^2 = SNR^((p-1)/3)
+    "0x2fb347984f7911f74c0bec3cf559b143b78cc310c2c3330c99e39557176f553d",
+    "0x16c9e55061ebae204ba4cc8bd75a079432ae2a1d0b7c9dce1665d51c640fcba2"
+  ),
+  Fp2[BN254_Snarks].fromHex( # SNR^((p-1)/6)^3 = SNR^((p-1)/2)
+    "0x63cf305489af5dcdc5ec698b6e2f9b9dbaae0eda9c95998dc54014671a0135a",
+    "0x7c03cbcac41049a0704b5a7ec796f2b21807dc98fa25bd282d37f632623b0e3"
+  ),
+  Fp2[BN254_Snarks].fromHex( # SNR^((p-1)/6)^4 = SNR^(2(p-1)/3)
+    "0x5b54f5e64eea80180f3c0b75a181e84d33365f7be94ec72848a1f55921ea762",
+    "0x2c145edbe7fd8aee9f3a80b03b0b1c923685d2ea1bdec763c13b4711cd2b8126"
+  ),
+  Fp2[BN254_Snarks].fromHex( # SNR^((p-1)/6)^5
+    "0x183c1e74f798649e93a3661a4353ff4425c459b55aa1bd32ea2c810eab7692f",
+    "0x12acf2ca76fd0675a27fb246c7729f7db080cb99678e2ac024c6b8ee6e0c2c4b"
+  )],
+  # frobenius(2)
+  [Fp2[BN254_Snarks].fromHex( # norm(SNR)^((p-1)/6)^1
+    "0x1",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex( # norm(SNR)^((p-1)/6)^2
+    "0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd49",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd48",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd46",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x59e26bcea0d48bacd4f263f1acdb5c4f5763473177fffffe",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x59e26bcea0d48bacd4f263f1acdb5c4f5763473177ffffff",
+    "0x0"
+  )],
+  # frobenius(3)
+  [Fp2[BN254_Snarks].fromHex(
+    "0x1",
+    "0x0"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x19dc81cfcc82e4bbefe9608cd0acaa90894cb38dbe55d24ae86f7d391ed4a67f",
+    "0xabf8b60be77d7306cbeee33576139d7f03a5e397d439ec7694aa2bf4c0c101"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x856e078b755ef0abaff1c77959f25ac805ffd3d5d6942d37b746ee87bdcfb6d",
+    "0x4f1de41b3d1766fa9f30e6dec26094f0fdf31bf98ff2631380cab2baaa586de"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x2a275b6d9896aa4cdbf17f1dca9e5ea3bbd689a3bea870f45fcc8ad066dce9ed",
+    "0x28a411b634f09b8fb14b900e9507e9327600ecc7d8cf6ebab94d0cb3b2594c64"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0xbc58c6611c08dab19bee0f7b5b2444ee633094575b06bcb0e1a92bc3ccbf066",
+    "0x23d5e999e1910a12feb0f6ef0cd21d04a44a9e08737f96e55fe3ed9d730c239f"
+  ),
+  Fp2[BN254_Snarks].fromHex(
+    "0x13c49044952c0905711699fa3b4d3f692ed68098967c84a5ebde847076261b43",
+    "0x16db366a59b1dd0b9fb1b2282a48633d3e2ddaea200280211f25041384282499"
+  )]]
+
+{.experimental: "dynamicBindSym".}
+
+macro frobMapConst(C: static Curve): untyped =
+  return bindSym("FrobMapConst_" & $C)
+
+func frobenius_map*[C](r: var Fp4[C], a: Fp4[C], 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_BLS12_381[k-1][3]
+  r.c1.mulCheckSparse frobMapConst(C)[k-1][3]

-func frobenius_map*(r: var Fp6, a: Fp6, k: static int = 1) {.inline.} =
+func frobenius_map*[C](r: var Fp6[C], a: Fp6[C], k: static int = 1) {.inline.} =
  ## Computes a^(p^k)
  ## The p-power frobenius automorphism on 𝔽p6
  r.c0.frobenius_map(a.c0, k)
  r.c1.frobenius_map(a.c1, k)
  r.c2.frobenius_map(a.c2, k)
-  r.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][2]
-  r.c2.mulCheckSparse FrobMapConst_BLS12_381[k-1][4]
+  r.c1.mulCheckSparse frobMapConst(C)[k-1][2]
+  r.c2.mulCheckSparse frobMapConst(C)[k-1][4]

-func frobenius_map*(r: var Fp12, a: Fp12, k: static int = 1) {.inline.} =
+func frobenius_map*[C](r: var Fp12[C], a: Fp12[C], k: static int = 1) {.inline.} =
  ## Computes a^(p^k)
  ## The p-power frobenius automorphism on 𝔽p12
  static: doAssert r.c0 is Fp4
@ -170,12 +329,12 @@ func frobenius_map*(r: var Fp12, a: Fp12, k: static int = 1) {.inline.} =
    for r_fp2, a_fp2 in fields(r_fp4, a_fp4):
      r_fp2.frobenius_map(a_fp2, k)

-  r.c0.c0.mulCheckSparse FrobMapConst_BLS12_381[k-1][0]
-  r.c0.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][3]
-  r.c1.c0.mulCheckSparse FrobMapConst_BLS12_381[k-1][1]
-  r.c1.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][4]
-  r.c2.c0.mulCheckSparse FrobMapConst_BLS12_381[k-1][2]
-  r.c2.c1.mulCheckSparse FrobMapConst_BLS12_381[k-1][5]
+  r.c0.c0.mulCheckSparse frobMapConst(C)[k-1][0]
+  r.c0.c1.mulCheckSparse frobMapConst(C)[k-1][3]
+  r.c1.c0.mulCheckSparse frobMapConst(C)[k-1][1]
+  r.c1.c1.mulCheckSparse frobMapConst(C)[k-1][4]
+  r.c2.c0.mulCheckSparse frobMapConst(C)[k-1][2]
+  r.c2.c1.mulCheckSparse frobMapConst(C)[k-1][5]

 # ψ (Psi) - Untwist-Frobenius-Twist Endomorphisms on twisted curves
 # -----------------------------------------------------------------
@ -270,8 +429,6 @@ const FrobPsiConst_BLS12_381_psi2_coef2 = Fp2[BLS12_381].fromHex(
  "0x0"
 )

-{.experimental: "dynamicBindSym".}
-
 macro frobPsiConst(C: static Curve, psipow, coefpow: static int): untyped =
  return bindSym("FrobPsiConst_" & $C & "_psi" & $psipow & "_coef" & $coefpow)

--- a/constantine/pairing/README.md
+++ b/constantine/pairing/README.md
@ -50,6 +50,11 @@
  Juan Manuel Gonzalez Nieto, and Kenneth Koon-Ho Wong, 2010
  https://eprint.iacr.org/2010/104.pdf

+- Faster hashing to G2\
+  Laura Fuentes-Castañeda, Edward Knapp,\
+  Francisco Jose Rodríguez-Henríquez, 2011\
+  https://link.springer.com/content/pdf/10.1007%2F978-3-642-28496-0_25.pdf
+
 - Pairings for beginners\
  Craig Costello, 2012 (?)\
  http://www.craigcostello.com.au/pairings/PairingsForBeginners.pdf
@ -65,7 +70,8 @@
  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\
+  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

@ -75,6 +81,10 @@
  https://scholarworks.rit.edu/cgi/viewcontent.cgi?referer=&httpsredir=1&article=10083&context=theses
  https://github.com/rajeevakarv/FiniteFieldComputations

+- Memory-saving computation of the pairing final exponentiation on BN curves\
+  Sylvain Duquesne and Loubna Ghammam, 2015\
+  https://eprint.iacr.org/2015/192
+
 - A taxonomy of pairings, their security, their complexity\
  Razvan Barbulescu, Nadia El Mrabet, and Loubna Ghammam, 2019\
  https://hal.archives-ouvertes.fr/hal-02129868/file/2019-485.pdf
--- a/constantine/pairing/cyclotomic_fp12.nim
+++ b/constantine/pairing/cyclotomic_fp12.nim
@ -224,3 +224,26 @@ func cyclotomic_square*[C](a: var Fp12[C]) =

  else:
    {.error: "Not implemented".}
+
+iterator unpack(scalarByte: byte): bool =
+  yield bool((scalarByte and 0b10000000) shr 7)
+  yield bool((scalarByte and 0b01000000) shr 6)
+  yield bool((scalarByte and 0b00100000) shr 5)
+  yield bool((scalarByte and 0b00010000) shr 4)
+  yield bool((scalarByte and 0b00001000) shr 3)
+  yield bool((scalarByte and 0b00000100) shr 2)
+  yield bool((scalarByte and 0b00000010) shr 1)
+  yield bool( scalarByte and 0b00000001)
+
+func cyclotomic_exp*[C](r: var Fp12[C], a: Fp12[C], exponent: BigInt, invert: bool) =
+    var eBytes: array[(exponent.bits+7) div 8, byte]
+    eBytes.exportRawUint(exponent, bigEndian)
+
+    r.setOne()
+    for b in eBytes:
+      for bit in unpack(b):
+        r.cyclotomic_square()
+        if bit:
+          r *= a
+    if invert:
+      r.cyclotomic_inv()
--- a/constantine/pairing/mul_fp12_by_lines.nim
+++ b/constantine/pairing/mul_fp12_by_lines.nim
@ -41,9 +41,15 @@ import
 # 𝔽p12 by line - Sparse functions
 # ----------------------------------------------------------------

+func mul_sparse_by_y0*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
+  ## Sparse multiplication of an Fp4 element
+  ## with coordinates (a₀, a₁) by (b₀, 0)
+  r.c0.prod(a.c0, b)
+  r.c1.prod(a.c1, b)
+
 func mul_sparse_by_0y*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
  ## Sparse multiplication of an Fp4 element
-  ## with coordinates (a₀, a₁) by (0, b₁, 0)
+  ## with coordinates (a₀, a₁) by (0, b₁)
  r.c0.prod(a.c1, b)
  r.c0 *= NonResidue
  r.c1.prod(a.c0, b)
@ -140,12 +146,49 @@ func mul_sparse_by_line_xyz000*[C: static Curve, Tw: static SexticTwist](

  static: doAssert f.c0.typeof is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4"

-  var v: Fp12[C]
-  v.c0.c0 = l.x
-  v.c0.c1 = l.y
-  v.c1.c0 = l.z
+  # In the following equations (taken from cubic extension implementation)
+  # a = f
+  # b0 = (x, y)
+  # b1 = (z, 0)
+  # b2 = (0, 0)
+  #
+  # v0 = a0 b0 = (f00, f01).(x, y)
+  # v1 = a1 b1 = (f10, f11).(z, 0)
+  # v2 = a2 b2 = (f20, f21).(0, 0)
+  #
+  # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
+  #    = ξ (a1 b1 + a2 b1 - v1) + v0
+  #    = ξ a2 b1 + v0
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2
+  #    = (a0 + a1) * (b0 + b1) - v0 - v1
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
+  #    = a0 b0 + a2 b0 - v0 + v1
+  #    = a2 b0 + v1

-  f *= v
+  var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fp4[C]
+
+  b0.c0 = l.x
+  b0.c1 = l.y
+
+  v0.prod(f.c0, b0)
+  v1.mul_sparse_by_y0(f.c1, l.z)
+
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1
+  f.c1 += f.c0 # r1 = a0 + a1
+  t = b0
+  t.c0 += l.z  # t = b0 + b1
+  f.c1 *= t    # r2 = (a0 + a1)(b0 + b1)
+  f.c1 -= v0
+  f.c1 -= v1   # r2 = (a0 + a1)(b0 + b1) - v0 - v1
+
+  # r0 = ξ a2 b1 + v0
+  f.c0.mul_sparse_by_y0(f.c2, l.z)
+  f.c0 *= NonResidue
+  f.c0 += v0
+
+  # r2 = a2 b0 + v1
+  f.c2 *= b0
+  f.c2 += v1

 func mul_sparse_by_line_xy000z*[C: static Curve, Tw: static SexticTwist](
       f: var Fp12[C], l: Line[Fp2[C], Tw]) =
--- a/constantine/pairing/pairing_bls12.nim
+++ b/constantine/pairing/pairing_bls12.nim
@ -44,10 +44,6 @@ import
 #   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
@ -242,8 +238,6 @@ func pairing_bls12*[C](gt: var Fp12[C], P: ECP_SWei_Proj[Fp[C]], Q: ECP_SWei_Pro
  ## 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)
--- a/constantine/pairing/pairing_bn.nim
+++ b/constantine/pairing/pairing_bn.nim
@ -29,21 +29,19 @@ import
 #
 # ############################################################

-# - 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
+# - Memory-saving computation of the pairing final exponentiation on BN curves
+#   Sylvain Duquesne and Loubna Ghammam, 2015
+#   https://eprint.iacr.org/2015/192
+#
+# - Faster hashing to G2
+#   Laura Fuentes-Castañeda, Edward Knapp,
+#   Francisco Jose Rodríguez-Henríquez, 2011
+#   https://link.springer.com/content/pdf/10.1007%2F978-3-642-28496-0_25.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
@ -51,16 +49,19 @@ const BN254_Snarks_ate_param = block:

 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

+# Generic slow pairing implementation
+# ----------------------------------------------------------------
+
+const BN254_Snarks_finalexponent = block:
+  # (p^12 - 1) / r
+  BigInt[2790].fromHex"0x2f4b6dc97020fddadf107d20bc842d43bf6369b1ff6a1c71015f3f7be2e1e30a73bb94fec0daf15466b2383a5d3ec3d15ad524d8f70c54efee1bd8c3b21377e563a09a1b705887e72eceaddea3790364a61f676baaf977870e88d5c6c8fef0781361e443ae77f5b63a2a2264487f2940a8b1ddb3d15062cd0fb2015dfc6668449aed3cc48a82d0d602d268c7daab6a41294c0cc4ebe5664568dfc50e1648a45a4a1e3a5195846a3ed011a337a02088ec80e0ebae8755cfe107acf3aafb40494e406f804216bb10cf430b0f37856b42db8dc5514724ee93dfb10826f0dd4a0364b9580291d2cd65664814fde37ca80bb4ea44eacc5e641bbadf423f9a2cbf813b8d145da90029baee7ddadda71c7f3811c4105262945bba1668c3be69a3c230974d83561841d766f9c9d570bb7fbe04c7e8a6c3c760c0de81def35692da361102b6b9b2b918837fa97896e84abb40a4efb7e54523a486964b64ca86f120"
+
 const BN254_Nogami_finalexponent = block:
  # (p^12 - 1) / r
  BigInt[2786].fromHex"0x2928fbb36b391596ee3fe4cbe857330da83e46fedf04d235a4a8daf5ff9f6eabcb4e3f20aa06f0a0d96b24f9af0cbbce750d61627dcbf5fec9139b8f1c46c86b49b4f8a202af26e4504f2c0f56570e9bd5b94c403f385d1908556486e24b396ddc2cdf13d06542f84fe8e82ccbad7b7423fc1ef4e8cc73d605e3e867c0a75f45ea7f6356d9846ce35d5a34f30396938818ad41914b97b99c289a7259b5d2e09477a77bd3c409b19f19e893f8ade90b0aed1b5fc8a07a3cebb41d4e9eee96b21a832ddb1e93e113edfb704fa532848c18593cd0ee90444a1b3499a800177ea38bdec62ec5191f2b6bbee449722f98d2173ad33077545c2ad10347e125a56fb40f086e9a4e62ad336a72c8b202ac3c1473d73b93d93dc0795ca0ca39226e7b4c1bb92f99248ec0806e0ad70744e9f2238736790f5185ea4c70808442a7d530c6ccd56b55a6973867ec6c73599bbd020bbe105da9c6b5c009ad8946cd6f0"
@ -70,7 +71,7 @@ const BN254_Nogami_finalexponent = block:
 macro get(C: static Curve, value: untyped): untyped =
  return bindSym($C & "_" & $value)

-func millerLoopGenericBN[C: static Curve](
+func millerLoopGenericBN*[C: static Curve](
       f: var Fp12[C],
       P: ECP_SWei_Aff[Fp[C]],
       Q: ECP_SWei_Aff[Fp2[C]]
@ -167,13 +168,169 @@ func pairing_bn_reference*[C](gt: var Fp12[C], P: ECP_SWei_Proj[Fp[C]], Q: ECP_S
  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)
+# Optimized pairing implementation
+# ----------------------------------------------------------------
+
+func cycl_sqr_repeated(f: var Fp12, num: int) =
+  ## Repeated cyclotomic squarings
+  for _ in 0 ..< num:
+    f.cyclotomic_square()
+
+func pow_u(r: var Fp12[BN254_Nogami], a: Fp12[BN254_Nogami], invert = BN254_Nogami_ate_param_isNeg) =
+  ## f^u with u the curve parameter
+  ## For BN254_Nogami f^-0x4080000000000001
+  r = a
+  r.cycl_sqr_repeated(7)
+  r *= a
+  r.cycl_sqr_repeated(55)
+  r *= a
+
+  if invert:
+    r.cyclotomic_inv()
+
+func pow_u(r: var Fp12[BN254_Snarks], a: Fp12[BN254_Snarks], invert = BN254_Snarks_ate_param_isNeg) =
+  ## f^u with u the curve parameter
+  ## For BN254_Snarks f^0x44e992b44a6909f1
+  when false:
+    cyclotomic_exp(
+      r, a,
+      BigInt[63].fromHex("0x44e992b44a6909f1"),
+      invert
+    )
+  else:
+    var # Hopefully the compiler optimizes away unused Fp12
+        # because those are huge
+      x10       {.noInit.}: Fp12[BN254_Snarks]
+      x11       {.noInit.}: Fp12[BN254_Snarks]
+      x100      {.noInit.}: Fp12[BN254_Snarks]
+      x110      {.noInit.}: Fp12[BN254_Snarks]
+      x1100     {.noInit.}: Fp12[BN254_Snarks]
+      x1111     {.noInit.}: Fp12[BN254_Snarks]
+      x10010    {.noInit.}: Fp12[BN254_Snarks]
+      x10110    {.noInit.}: Fp12[BN254_Snarks]
+      x11100    {.noInit.}: Fp12[BN254_Snarks]
+      x101110   {.noInit.}: Fp12[BN254_Snarks]
+      x1001010  {.noInit.}: Fp12[BN254_Snarks]
+      x1111000  {.noInit.}: Fp12[BN254_Snarks]
+      x10001110 {.noInit.}: Fp12[BN254_Snarks]
+
+    x10       .cyclotomic_square(a)
+    x11       .prod(x10, a)
+    x100      .prod(x11, a)
+    x110      .prod(x10, x100)
+    x1100     .cyclotomic_square(x110)
+    x1111     .prod(x11, x1100)
+    x10010    .prod(x11, x1111)
+    x10110    .prod(x100, x10010)
+    x11100    .prod(x110, x10110)
+    x101110   .prod(x10010, x11100)
+    x1001010  .prod(x11100, x101110)
+    x1111000  .prod(x101110, x1001010)
+    x10001110 .prod(x10110, x1111000)
+
+    var
+      i15 {.noInit.}: Fp12[BN254_Snarks]
+      i16 {.noInit.}: Fp12[BN254_Snarks]
+      i17 {.noInit.}: Fp12[BN254_Snarks]
+      i18 {.noInit.}: Fp12[BN254_Snarks]
+      i20 {.noInit.}: Fp12[BN254_Snarks]
+      i21 {.noInit.}: Fp12[BN254_Snarks]
+      i22 {.noInit.}: Fp12[BN254_Snarks]
+      i26 {.noInit.}: Fp12[BN254_Snarks]
+      i27 {.noInit.}: Fp12[BN254_Snarks]
+      i61 {.noInit.}: Fp12[BN254_Snarks]
+
+    i15.cyclotomic_square(x10001110)
+    i15 *= x1001010
+    i16.prod(x10001110, i15)
+    i17.prod(x1111, i16)
+    i18.prod(i16, i17)
+
+    i20.cyclotomic_square(i18)
+    i20 *= i17
+    i21.prod(x1111000, i20)
+    i22.prod(i15, i21)
+
+    i26.cyclotomic_square(i22)
+    i26.cyclotomic_square()
+    i26 *= i22
+    i26 *= i18
+
+    i27.prod(i22, i26)
+
+    i61.prod(i26, i27)
+    i61.cycl_sqr_repeated(17)
+    i61 *= i27
+    i61.cycl_sqr_repeated(14)
+    i61 *= i21
+
+    r = i61
+    r.cycl_sqr_repeated(16)
+    r *= i20
+
+    if invert:
+      r.cyclotomic_inv()
+
+func finalExpHard_BN*[C: static Curve](f: var Fp12[C]) =
+  ## Hard part of the final exponentiation
+  ## Specialized for BN curves
  ##
-  ## With an embedding degree of 12, the easy part of final exponentiation is
-  ##
-  ##  f^(p⁶−1)(p²+1)
-  discard
+  # - Memory-saving computation of the pairing final exponentiation on BN curves
+  #   Sylvain Duquesne and Loubna Ghammam, 2015
+  #   https://eprint.iacr.org/2015/192
+  # - Faster hashing to G2
+  #   Laura Fuentes-Castañeda, Edward Knapp,
+  #   Francisco Jose Rodríguez-Henríquez, 2011
+  #   https://link.springer.com/content/pdf/10.1007%2F978-3-642-28496-0_25.pdf
+  #
+  # We use the Fuentes-Castañeda et al algorithm without
+  # memory saving optimization
+  # as that variant has an exponentiation by -2u-1
+  # that requires another addition chain
+  var t0 {.noInit.}, t1 {.noinit.}, t2 {.noinit.}, t3 {.noinit.}, t4 {.noinit.}: Fp12[C]
+
+  t0.pow_u(f, invert = false)  # t0 = f^|u|
+  t0.cyclotomic_square()       # t0 = f^2|u|
+  t1.cyclotomic_square(t0)     # t1 = f^4|u|
+  t1 *= t0                     # t1 = f^6|u|
+  t2.pow_u(t1, invert = false) # t2 = f^6u²
+
+  if C.get(ate_param_is_Neg):
+    t3.cyclotomic_inv(t1)      # t3 = f^6u
+  else:
+    t3 = t1                    # t3 = f^6u
+  t1.prod(t2, t3)              # t1 = f^6u.f^6u²
+  t3.cyclotomic_square(t2)     # t3 = f^12u²
+  t4.pow_u(t3)                 # t4 = f^12u³
+  t4 *= t1                     # t4 = f^(6u + 6u² + 12u³) = f^λ₂
+
+  if not C.get(ate_param_is_Neg):
+    t0.cyclotomic_inv()        # t0 = f^-2u
+  t3.prod(t4, t0)              # t3 = f^(4u + 6u² + 12u³)
+
+  t0.prod(t2, t4)              # t0 = f^6u.f^12u².f^12u³
+  t0 *= f                      # t0 = f^(1 + 6u + 12u² + 12u³) = f^λ₀
+
+  t2.frobenius_map(t3)         # t2 = f^(4u + 6u² + 12u³)p = f^λ₁p
+  t0 *= t2                     # t0 = f^(λ₀+λ₁p)
+
+  t2.frobenius_map(t4, 2)      # t2 = f^λ₂p²
+  t0 *= t2                     # t0 = f^(λ₀ + λ₁p + λ₂p²)
+
+  t2.cyclotomic_inv(f)         # t2 = f⁻¹
+  t2 *= t3                     # t3 = f^(-1 + 4u + 6u² + 12u³) = f^λ₃
+
+  f.frobenius_map(t2, 3)       # r = f^λ₃p³
+  f *= t0                      # r = f^(λ₀ + λ₁p + λ₂p² + λ₃p³) = f^((p⁴-p²+1)/r)
+
+func pairing_bn*[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
+  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.finalExpEasy()
+  gt.finalExpHard_BN()
--- a/sage/frobenius_bls12_381.sage
+++ b/sage/frobenius_bls12_381.sage
@ -50,7 +50,7 @@ def fp2_to_hex(a):
    v = vector(a)
    return '0x' + Integer(v[0]).hex() + ' + β * ' + '0x' + Integer(v[1]).hex()

-# Frobenius map constants (D type: use SNR, M type use 1/SNR)
+# Frobenius map constants
 print('\nFrobenius extension field constants')
 FrobConst_map = SNR^((p-1)/6)
 FrobConst_map_list = []
--- a/sage/frobenius_bn254_nogami.sage
+++ b/sage/frobenius_bn254_nogami.sage
@ -49,9 +49,27 @@ print('')
 # Utilities
 def fp2_to_hex(a):
    v = vector(a)
-    return Integer(v[0]).hex() + ' + β * ' + Integer(v[1]).hex()
+    return '0x' + Integer(v[0]).hex() + ' + β * ' + '0x' + Integer(v[1]).hex()

-# Frobenius constants (D type: use SNR, M type use 1/SNR)
+# Frobenius map constants
+print('\nFrobenius extension field constants')
+FrobConst_map = SNR^((p-1)/6)
+FrobConst_map_list = []
+cur = Fp2([1, 0])
+
+for i in range(6):
+    FrobConst_map_list.append(cur)
+    print(f'FrobConst_map_{i}     : {fp2_to_hex(cur)}')
+    cur *= FrobConst_map
+print('')
+for i in range(6):
+    print(f'FrobConst_map_{i}_pow2     : {fp2_to_hex(FrobConst_map_list[i]*conjugate(FrobConst_map_list[i]))}')
+print('')
+for i in range(6):
+    print(f'FrobConst_map_{i}_pow3     : {fp2_to_hex(FrobConst_map_list[i]**2 * conjugate(FrobConst_map_list[i]))}')
+
+# Frobenius psi constants (D type: use SNR, M type use 1/SNR)
+print('\nψ (Psi) - Untwist-Frobenius-Twist constants')
 FrobConst_psi = SNR^((p-1)/6)
 FrobConst_psi_2 = FrobConst_psi * FrobConst_psi
 FrobConst_psi_3 = FrobConst_psi_2 * FrobConst_psi
--- a/sage/frobenius_bn254_snarks.sage
+++ b/sage/frobenius_bn254_snarks.sage
@ -49,9 +49,27 @@ print('')
 # Utilities
 def fp2_to_hex(a):
    v = vector(a)
-    return Integer(v[0]).hex() + ' + β * ' + Integer(v[1]).hex()
+    return '0x' + Integer(v[0]).hex() + ' + β * ' + '0x' + Integer(v[1]).hex()

-# Frobenius constants (D type: use SNR, M type use 1/SNR)
+# Frobenius map constants
+print('\nFrobenius extension field constants')
+FrobConst_map = SNR^((p-1)/6)
+FrobConst_map_list = []
+cur = Fp2([1, 0])
+
+for i in range(6):
+    FrobConst_map_list.append(cur)
+    print(f'FrobConst_map_{i}     : {fp2_to_hex(cur)}')
+    cur *= FrobConst_map
+print('')
+for i in range(6):
+    print(f'FrobConst_map_{i}_pow2     : {fp2_to_hex(FrobConst_map_list[i]*conjugate(FrobConst_map_list[i]))}')
+print('')
+for i in range(6):
+    print(f'FrobConst_map_{i}_pow3     : {fp2_to_hex(FrobConst_map_list[i]**2 * conjugate(FrobConst_map_list[i]))}')
+
+# Frobenius psi constants (D type: use SNR, M type use 1/SNR)
+print('\nψ (Psi) - Untwist-Frobenius-Twist constants')
 FrobConst_psi = SNR^((p-1)/6)
 FrobConst_psi_2 = FrobConst_psi * FrobConst_psi
 FrobConst_psi_3 = FrobConst_psi_2 * FrobConst_psi
--- a/tests/t_fp12_frobenius.nim
+++ b/tests/t_fp12_frobenius.nim
@ -14,8 +14,8 @@ import
  ./t_fp_tower_frobenius_template

 const TestCurves = [
-    # BN254_Nogami
-    # BN254_Snarks,
+    BN254_Nogami,
+    BN254_Snarks,
    # BLS12_377,
    BLS12_381,
    # BN446
--- a/tests/t_fp2_frobenius.nim
+++ b/tests/t_fp2_frobenius.nim
@ -14,7 +14,7 @@ import
  ./t_fp_tower_frobenius_template

 const TestCurves = [
-    # BN254_Nogami
+    BN254_Nogami,
    BN254_Snarks,
    BLS12_377,
    BLS12_381,
--- a/tests/t_fp4_frobenius.nim
+++ b/tests/t_fp4_frobenius.nim
@ -14,8 +14,8 @@ import
  ./t_fp_tower_frobenius_template

 const TestCurves = [
-    # BN254_Nogami
-    # BN254_Snarks,
+    BN254_Nogami,
+    BN254_Snarks,
    # BLS12_377,
    BLS12_381,
    # BN446
--- a/tests/t_fp6_frobenius.nim
+++ b/tests/t_fp6_frobenius.nim
@ -14,8 +14,8 @@ import
  ./t_fp_tower_frobenius_template

 const TestCurves = [
-    # BN254_Nogami
-    # BN254_Snarks,
+    BN254_Nogami,
+    BN254_Snarks,
    # BLS12_377,
    BLS12_381,
    # BN446
--- a/tests/t_pairing_bn254_nogami_optate.nim
+++ b/tests/t_pairing_bn254_nogami_optate.nim
@ -13,4 +13,5 @@ import
  # Test utilities
  ./t_pairing_template

-runPairingTests(4, BN254_Nogami, pairing_bn_reference)
+# runPairingTests(4, BN254_Nogami, pairing_bn_reference)
+runPairingTests(4, BN254_Nogami, pairing_bn)
--- a/tests/t_pairing_bn254_snarks_optate.nim
+++ b/tests/t_pairing_bn254_snarks_optate.nim
@ -13,4 +13,5 @@ import
  # Test utilities
  ./t_pairing_template

-runPairingTests(4, BN254_Snarks, pairing_bn_reference)
+# runPairingTests(4, BN254_Snarks, pairing_bn_reference)
+runPairingTests(4, BN254_Snarks, pairing_bn)
--- a/tests/t_pairing_mul_fp12_by_lines.nim
+++ b/tests/t_pairing_mul_fp12_by_lines.nim
@ -25,7 +25,8 @@ import
 const
  Iters = 8
  TestCurves = [
-    # BN254_Snarks,
+    BN254_Nogami,
+    BN254_Snarks,
    BLS12_381
  ]

@ -189,3 +190,30 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
        test_fp12_xy000z(curve, gen = Uniform)
        test_fp12_xy000z(curve, gen = HighHammingWeight)
        test_fp12_xy000z(curve, gen = Long01Sequence)
+
+    test "Sparse 𝔽p12/𝔽p4 resulting from xyz000 line function":
+      proc test_fp12_xy000z(C: static Curve, gen: static RandomGen) =
+        for _ in 0 ..< Iters:
+          var a = rng.random_elem(Fp12[C], gen)
+          var a2 = a
+
+          var x = rng.random_elem(Fp2[C], gen)
+          var y = rng.random_elem(Fp2[C], gen)
+          var z = rng.random_elem(Fp2[C], gen)
+
+          let line = Line[Fp2[C], Dtwist](x: x, y: y, z: z)
+          let b = Fp12[C](
+            c0: Fp4[C](c0: x, c1: y),
+            c1: Fp4[C](c0: z       ),
+            # c2:
+          )
+
+          a *= b
+          a2.mul_sparse_by_line_xyz000(line)
+
+          check: bool(a == a2)
+
+      staticFor(curve, TestCurves):
+        test_fp12_xy000z(curve, gen = Uniform)
+        test_fp12_xy000z(curve, gen = HighHammingWeight)
+        test_fp12_xy000z(curve, gen = Long01Sequence)