[WIP] Pairings for bw6 761 (#108)

* Prepare BW6-761 pairing constants

* Extract the basic miller loop from pairings

* template and method call syntax issue

* Layout pairing for BW6-761

* Fix rebasing woes

* Try to match the paper (still buggy)

* Stash BW6-761

constantine/curves/bw6_761_pairing.nim

@@ -0,0 +1,57 @@
+# 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
+  ../config/[curves, type_bigint],
+  ../io/io_bigints
+
+# Slow generic implementation
+# ------------------------------------------------------------
+
+# 1st part: f_{u+1,Q}(P)
+const BW6_761_pairing_ate_param_1_unopt* = block:
+  # BW6-761 unoptimized Miller loop first part is parametrized by u+1
+  # +1 to bitlength so that we can mul by 3 for NAF encoding
+  BigInt[64+1].fromHex"0x8508c00000000002"
+
+const BW6_761_pairing_ate_param_1_unopt_isNeg* = false
+
+
+# 2nd part: f_{u*(u²-u-1),Q}(P) followed by Frobenius application
+const BW6_761_pairing_ate_param_2_unopt* = block:
+  # BW6 unoptimized Miller loop second part is parametrized by u*(u²-u-1)
+  # +1 to bitlength so that we can mul by 3 for NAF encoding
+  BigInt[190+1].fromHex"0x23ed1347970dec008a442f991fffffffffffffffffffffff"
+
+const BW6_761_pairing_ate_param_2_unopt_isNeg* = false
+
+
+# 1st part: f_{u,Q}(P)
+const BW6_761_pairing_ate_param_1_opt* = block:
+  # BW6 Miller loop first part is parametrized by u
+  # no NAF for the optimized first Miller loop
+  BigInt[64].fromHex"0x8508c00000000001"
+
+const BW6_761_pairing_ate_param_1_opt_isNeg* = false
+
+
+# 2nd part: f_{u²-u-1,Q}(P) followed by Frobenius application
+const BW6_761_pairing_ate_param_opt_2* = block:
+  # BW6 Miller loop second part is parametrized by u²-u-1
+  # +1 to bitlength so that we can mul by 3 for NAF encoding
+  BigInt[127+1].fromHex"0x452217cc900000008508bfffffffffff"
+
+const BW6_761_pairing_ate_param_2_opt_isNeg* = false
+
+
+const BW6_761_pairing_finalexponent* = block:
+  # (p^6 - 1) / r * 3
+  BigInt[4186].fromHex"0x3d7fafd4d00189a67bdf3e3e099095571b3671b450e1430228baeca99efec770d2499a6732e8891ede83d26c08c7afdcb004a074ccea612933db92ba5b26a6683f2b782d91befd4170c3203b47ecb246847cd292b51591c00f608b6bd51942243a3042325356d537c26dc5cbe2c64656bf2aed4b94c66bf8629eb027698ebde2b14cbeda063db5d74b44c16ffd421206094832fe5b7ec54d68e312f5bfa26f87ea2c85578de4a05d1283d040a9a13ee0c9b4dfaf4116599b14ffbde13fb06415e28945def8dc5ada9692d40c49b675718ca8865551b0cca4c87bbb2becd0a90db08638c5bd777015ae4f34d19c66bb5de3e9929deb7de11789fb4100a0d1bbd75cabb2d52979693cd2f2c7bbb77016161f43722b3b1a32f3cf150df07f282193c7bd573c046e3b17775c3f007b2ba146b8fd2434604c0f29fb56edf981d37ad4c312c3daa27314b14db0d0c4d030a5dd7641899e685efcb9d41791a84ed44ef6b8f6f86522ef26e63f53693df95706fc1264a062f93d499cfdc033465a582b86fe0329b011a4536505fcd30aa0e09dfc3c57fc4a9e95246d4d4519160cb6088828f5082ebc1775012c6868441ee831d897fabe8de92fde56533968e8bd25fd04cccb2d932f768350e8b0eaebbcab3649380640e01daba898ae6c5085a149cf14bb0b2f465391d8393298b2c3caf1c30a8496035a5c00c8327c30f7d1d5c24f02a65f7d3d0b413ade8564b78"
+
+# Addition chain
+# ------------------------------------------------------------

constantine/curves/zoo_pairings.nim

@@ -12,7 +12,8 @@ import
   ./bls12_377_pairing,
   ./bls12_381_pairing,
   ./bn254_nogami_pairing,
-  ./bn254_snarks_pairing
+  ./bn254_snarks_pairing,
+  ./bw6_761_pairing
 
 {.experimental: "dynamicBindSym".}
 

constantine/hash_to_curve/cofactors.nim

@@ -34,13 +34,21 @@ const Cofactor_Eff_BN254_Snarks_G2 = BigInt[254].fromHex"0x30644e72e131a029b8504
 const Cofactor_Eff_BLS12_377_G1 = BigInt[125].fromHex"0x170b5d44300000000000000000000000"
   ## P -> (1 - x) P
 const Cofactor_Eff_BLS12_377_G2 = BigInt[502].fromHex"0x26ba558ae9562addd88d99a6f6a829fbb36b00e1dcc40c8c505634fae2e189d693e8c36676bd09a0f3622fba094800452217cc900000000000000000000001"
-  ## P -> (x^2 - x - 1) P + (x - 1) psi(P) + psi(psi(2P))
+  ## P -> (x^2 - x - 1) P + (x - 1) ψ(P) + ψ(ψ(2P))
 
 # 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))
+  ## P -> (x^2 - x - 1) P + (x - 1) ψ(P) + ψ(ψ(2P))
+
+# TODO https://eprint.iacr.org/2020/351.pdf p12
+const Cofactor_Eff_BW6_761_G1 = BigInt[384].fromHex"0xad1972339049ce762c77d5ac34cb12efc856a0853c9db94cc61c554757551c0c832ba4061000003b3de580000000007c"
+  ## P -> 103([u³]P)− 83([u²]P)−40([u]P)+136P + φ(7([u²]P)+89([u]P)+130P)
+
+# TODO https://eprint.iacr.org/2020/351.pdf p13
+const Cofactor_Eff_BW6_761_G2 = BigInt[384].fromHex"0xad1972339049ce762c77d5ac34cb12efc856a0853c9db94cc61c554757551c0c832ba4061000003b3de580000000007c"
+  ## P -> (103([u³]P) − 83([u²]P) − 143([u]P) + 27P) + ψ(7([u²]P) − 117([u]P) − 109P)
 
 func clearCofactorReference*(P: var ECP_ShortW_Prj[Fp[BN254_Nogami], NotOnTwist]) {.inline.} =
   ## Clear the cofactor of BN254_Nogami G1
@@ -79,3 +87,12 @@ func clearCofactorReference*(P: var ECP_ShortW_Prj[Fp2[BLS12_381], OnTwist]) {.i
   ## Clear the cofactor of BLS12_381 G2
   # Endomorphism acceleration cannot be used if cofactor is not cleared
   P.scalarMulGeneric(Cofactor_Eff_BLS12_381_G2)
+
+func clearCofactorReference*(P: var ECP_ShortW_Prj[Fp[BW6_761], NotOnTwist]) {.inline.} =
+  ## Clear the cofactor of BW6_761 G1
+  P.scalarMulGeneric(Cofactor_Eff_BW6_761_G1)
+
+func clearCofactorReference*(P: var ECP_ShortW_Prj[Fp[BW6_761], OnTwist]) {.inline.} =
+  ## Clear the cofactor of BW6_761 G2
+  # Endomorphism acceleration cannot be used if cofactor is not cleared
+  P.scalarMulGeneric(Cofactor_Eff_BW6_761_G2)

constantine/pairing/miller_loops.nim

@@ -0,0 +1,59 @@
+# 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
+  ../elliptic/[
+    ec_shortweierstrass_affine,
+    ec_shortweierstrass_projective
+  ],
+  ./lines_projective,
+  ./mul_fp6_by_lines, ./mul_fp12_by_lines,
+  ../curves/zoo_pairings
+
+# ############################################################
+#
+#                 Basic Miller Loop
+#
+# ############################################################
+
+template basicMillerLoop*[FT, F1, F2](
+       f: var FT,
+       T: var ECP_ShortW_Prj[F2, OnTwist],
+       line: var Line[F2],
+       P: ECP_ShortW_Aff[F1, NotOnTwist],
+       Q, nQ: ECP_ShortW_Aff[F2, OnTwist],
+       ate_param: untyped,
+       ate_param_isNeg: untyped
+    ) =
+  ## Basic Miller loop iterations
+  static:
+    doAssert FT.C == F1.C
+    doAssert FT.C == F2.C
+
+  f.setOne()
+
+  template u: untyped = pairing(C, ate_param)
+  var u3 = pairing(C, ate_param)
+  u3 *= 3
+  for i in countdown(u3.bits - 2, 1):
+    square(f)
+    line_double(line, T, P)
+    mul(f, line)
+
+    let naf = bit(u3, i).int8 - bit(u, i).int8 # This can throw exception
+    if naf == 1:
+      line_add(line, T, Q, P)
+      mul(f, line)
+    elif naf == -1:
+      line_add(line, T, nQ, P)
+      mul(f, line)
+
+  when pairing(C, ate_param_isNeg):
+    # In GT, x^-1 == conjugate(x)
+    # Remark 7.1, chapter 7.1.1 of Guide to Pairing-Based Cryptography, El Mrabet, 2017
+    conj(f)

constantine/pairing/mul_fp6_by_lines.nim

@@ -0,0 +1,140 @@
+# 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/curves,
+  ../arithmetic,
+  ../towers,
+  ./lines_projective
+
+
+# ############################################################
+#
+#                 Sparse Multiplication
+#                        by lines
+#
+# ############################################################
+
+# 𝔽p6 by line - Sparse functions
+# ----------------------------------------------------------------
+
+func mul_sparse_by_line_xyz000*[C: static Curve](
+       f: var Fp6[C], l: Line[Fp[C]]) =
+  ## 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 C.getSexticTwist() == D_Twist
+    doAssert f.c0.typeof is Fp2, "This assumes 𝔽p6 as a cubic extension of 𝔽p6"
+
+  # 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
+
+  var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fp2[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](
+       f: var Fp6[C], l: Line[Fp[C]]) =
+
+  static:
+    doAssert C.getSexticTwist() == M_Twist
+    doAssert f.c0.typeof is Fp2, "This assumes 𝔽p6 as a cubic extension of 𝔽p2"
+
+  # In the following equations (taken from cubic extension implementation)
+  # a = f
+  # b0 = (x, y)
+  # b1 = (0, 0)
+  # b2 = (0, z)
+  #
+  # v0 = a0 b0 = (f00, f01).(x, y)
+  # v1 = a1 b1 = (f10, f11).(0, 0)
+  # v2 = a2 b2 = (f20, f21).(0, z)
+  #
+  # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
+  #    = ξ (a1 b2 + a2 b2 - v2) + v0
+  #    = ξ a1 b2 + v0
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2
+  #    = a0 b0 + a1 b0 - v0 + ξ v2
+  #    = a1 b0 + ξ v2
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
+  #    = (a0 + a2) * (b0 + b2) - v0 - v2
+
+  var b0 {.noInit.}, v0{.noInit.}, v2{.noInit.}, t{.noInit.}: Fp2[C]
+
+  b0.c0 = l.x
+  b0.c1 = l.y
+
+  v0.prod(f.c0, b0)
+  v2.mul_sparse_by_0y(f.c2, l.z)
+
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2
+  f.c2 += f.c0 # r2 = a0 + a2
+  t = b0
+  t.c1 += l.z  # t = b0 + b2
+  f.c2 *= t    # r2 = (a0 + a2)(b0 + b2)
+  f.c2 -= v0
+  f.c2 -= v2   # r2 = (a0 + a2)(b0 + b2) - v0 - v2
+
+  # r0 = ξ a1 b2 + v0
+  f.c0.mul_sparse_by_0y(f.c1, l.z)
+  f.c0 *= NonResidue
+  f.c0 += v0
+
+  # r1 = a1 b0 + ξ v2
+  f.c1 *= b0
+  v2 *= NonResidue
+  f.c1 += v2
+
+func mul*[C](f: var Fp6[C], line: Line[Fp[C]]) {.inline.} =
+  when C.getSexticTwist() == D_Twist:
+    f.mul_sparse_by_line_xyz000(line)
+  elif C.getSexticTwist() == M_Twist:
+    f.mul_sparse_by_line_xy000z(line)
+  else:
+    {.error: "A line function assumes that the curve has a twist".}

constantine/pairing/pairing_bls12.nim

@@ -14,10 +14,10 @@ import
     ec_shortweierstrass_projective
   ],
   ../isogeny/frobenius,
-  ./lines_projective,
+  ../curves/zoo_pairings,
-  ./mul_fp12_by_lines,
   ./cyclotomic_fp12,
-  ../curves/zoo_pairings
+  ./lines_common,
+  ./miller_loops
 
 # ############################################################
 #
@@ -53,32 +53,6 @@ func millerLoopGenericBLS12*[C](
   ## Generic Miller Loop for BLS12 curve
   ## Computes f{u,Q}(P) with u the BLS curve parameter
 
-  # 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
   var
     T {.noInit.}: ECP_ShortW_Prj[Fp2[C], OnTwist]
     line {.noInit.}: Line[Fp2[C]]
@@ -86,29 +60,12 @@ func millerLoopGenericBLS12*[C](
 
   T.projectiveFromAffine(Q)
   nQ.neg(Q)
-  f.setOne()
 
-  template u: untyped = C.pairing(ate_param)
+  basicMillerLoop(
-  var u3 = C.pairing(ate_param)
+    f, T, line,
-  u3 *= 3
+    P, Q, nQ,
-  for i in countdown(u3.bits - 2, 1):
+    ate_param, ate_param_isNeg
-    f.square()
+  )
-    line.line_double(T, P)
-
-    f.mul(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(line)
-    elif naf == -1:
-      line.line_add(T, nQ, P)
-      f.mul(line)
-
-  when C.pairing(ate_param_isNeg):
-    # 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

constantine/pairing/pairing_bn.nim

@@ -13,11 +13,12 @@ import
     ec_shortweierstrass_affine,
     ec_shortweierstrass_projective
   ],
+  ../isogeny/frobenius,
+  ../curves/zoo_pairings,
   ./lines_projective,
   ./mul_fp12_by_lines,
   ./cyclotomic_fp12,
-  ../isogeny/frobenius,
+  ./miller_loops
-  ../curves/zoo_pairings
 
 # ############################################################
 #
@@ -50,33 +51,6 @@ func millerLoopGenericBN*[C](
   ## Generic Miller Loop for BN curves
   ## Computes f{6u+2,Q}(P) with u the BN curve parameter
 
-  # 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_ShortW_Prj[Fp2[C], OnTwist]
     line {.noInit.}: Line[Fp2[C]]
@@ -84,30 +58,14 @@ func millerLoopGenericBN*[C](
 
   T.projectiveFromAffine(Q)
   nQ.neg(Q)
-  f.setOne()
 
-  template u: untyped = C.pairing(ate_param)
+  basicMillerLoop(
-  var u3 = C.pairing(ate_param)
+    f, T, line,
-  u3 *= 3
+    P, Q, nQ,
-  for i in countdown(u3.bits - 2, 1):
+    ate_param, ate_param_isNeg
-    f.square()
+  )
-    line.line_double(T, P)
-    f.mul(line)
 
-    let naf = u3.bit(i).int8 - u.bit(i).int8 # This can throw exception
+  # Ate pairing for BN curves need adjustment after basic Miller loop
-    if naf == 1:
-      line.line_add(T, Q, P)
-      f.mul(line)
-    elif naf == -1:
-      line.line_add(T, nQ, P)
-      f.mul(line)
-
-  when C.pairing(ate_param_isNeg):
-    # 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.pairing(ate_param_isNeg):
     T.neg()
   var V {.noInit.}: typeof(Q)

constantine/pairing/pairing_bw6_761.nim

@@ -0,0 +1,165 @@
+# 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
+  ../config/[curves, type_ff],
+  ../arithmetic,
+  ../towers,
+  ../elliptic/[
+    ec_shortweierstrass_affine,
+    ec_shortweierstrass_projective
+  ],
+  ../isogeny/frobenius,
+  ../curves/zoo_pairings,
+  ./lines_projective,
+  ./mul_fp6_by_lines,
+  ./miller_loops
+
+# ############################################################
+#
+#                 Optimal ATE pairing for
+#                      BW6-761 curve
+#
+# ############################################################
+
+# Generic pairing implementation
+# ----------------------------------------------------------------
+# TODO: debug this
+
+func millerLoopBW6_761_naive[C](
+       f: var Fp6[C],
+       P: ECP_ShortW_Aff[Fp[C], NotOnTwist],
+       Q: ECP_ShortW_Aff[Fp[C], OnTwist]
+     ) =
+  ## Miller Loop for BW6_761 curve
+  ## Computes f_{u+1,Q}(P)*Frobenius(f_{u*(u^2-u-1),Q}(P))
+
+  var
+    T {.noInit.}: ECP_ShortW_Prj[Fp[C], OnTwist]
+    line {.noInit.}: Line[Fp[C]]
+    nQ{.noInit.}: typeof(Q)
+
+  T.projectiveFromAffine(Q)
+  nQ.neg(Q)
+
+  basicMillerLoop(
+    f, T, line,
+    P, Q, nQ,
+    ate_param_1_unopt, ate_param_1_unopt_isNeg
+  )
+
+  var f2 {.noInit.}: typeof(f)
+  T.projectiveFromAffine(Q)
+
+  basicMillerLoop(
+    f2, T, line,
+    P, Q, nQ,
+    ate_param_1_unopt, ate_param_1_unopt_isNeg
+  )
+
+  let t = f2
+  f2.frobenius_map(t)
+  f *= f2
+
+func finalExpGeneric[C: static Curve](f: var Fp6[C]) =
+  ## A generic and slow implementation of final exponentiation
+  ## for sanity checks purposes.
+  f.powUnsafeExponent(C.pairing(finalexponent), window = 3)
+
+# Optimized pairing implementation
+# ----------------------------------------------------------------
+
+func millerLoopBW6_761_opt_to_debug[C](
+       f: var Fp6[C],
+       P: ECP_ShortW_Aff[Fp[C], NotOnTwist],
+       Q: ECP_ShortW_Aff[Fp[C], OnTwist]
+     ) {.used.} =
+  ## Miller Loop Otpimized for BW6_761 curve
+
+  # 1st part: f_{u,Q}(P)
+  # ------------------------------
+  var
+    T {.noInit.}: ECP_ShortW_Prj[Fp[C], OnTwist]
+    line {.noInit.}: Line[Fp[C]]
+
+  T.projectiveFromAffine(Q)
+  f.setOne()
+
+  template u: untyped = pairing(C, ate_param_1_opt)
+  for i in countdown(u.bits - 2, 1):
+    square(f)
+    line_double(line, T, P)
+    mul(f, line)
+
+    let bit = u.bit(i).int8
+    if bit == 1:
+      line_add(line, T, Q, P)
+      mul(f, line)
+
+  # Fixup
+  # ------------------------------
+  var minvu {.noInit.}, mu {.noInit.}, muplusone: typeof(f)
+  var Qu {.noInit.}, nQu {.noInit.}: typeof(Q)
+
+  mu = f
+  minvu.inv(f)
+  Qu.affineFromProjective(T)
+  nQu.neg(Qu)
+
+  # Drop the vertical line
+  line.line_add(T, Q, P) # TODO: eval without updating T
+  muplusone = mu
+  muplusone.mul(line)
+
+  # 2nd part: f_{u²-u-1,Q}(P)
+  # ------------------------------
+  # We restart from `f` and `T`
+  T.projectiveFromAffine(Qu)
+
+  template u: untyped = pairing(C, ate_param_2_opt)
+  var u3 = pairing(C, ate_param_2_opt)
+  u3 *= 3
+  for i in countdown(u3.bits - 2, 1):
+    square(f)
+    line_double(line, T, P)
+    mul(f, line)
+
+    let naf = bit(u3, i).int8 - bit(u, i).int8 # This can throw exception
+    if naf == 1:
+      line_add(line, T, Qu, P)
+      mul(f, line)
+      f *= mu
+    elif naf == -1:
+      line_add(line, T, nQu, P)
+      mul(f, line)
+      f *= minvu
+
+  # Final
+  # ------------------------------
+  let t = f
+  f.frobenius_map(t)
+  f *= muplusone
+
+# Public
+# ----------------------------------------------------------------
+
+func pairing_bw6_761_reference*[C](
+       gt: var Fp6[C],
+       P: ECP_ShortW_Prj[Fp[C], NotOnTwist],
+       Q: ECP_ShortW_Prj[Fp[C], OnTwist]) =
+  ## Compute the optimal Ate Pairing for BW6 curves
+  ## Input: P ∈ G1, Q ∈ G2
+  ## Output: e(P, Q) ∈ Gt
+  ##
+  ## Reference implementation
+  var Paff {.noInit.}: ECP_ShortW_Aff[Fp[C], NotOnTwist]
+  var Qaff {.noInit.}: ECP_ShortW_Aff[Fp[C], OnTwist]
+  Paff.affineFromProjective(P)
+  Qaff.affineFromProjective(Q)
+  gt.millerLoopBW6_761_naive(Paff, Qaff)
+  gt.finalExpGeneric()

sage/derive_pairing.sage

@@ -50,6 +50,11 @@ def genAteParam(curve_name, curve_config):
   elif family == 'BN':
     ate_param = 6*u+2
     ate_comment = '  # BN Miller loop is parametrized by 6u+2\n'
+  elif family == 'BW6':
+    result = genAteParam_BW6_unoptimized(curve_name, curve_config)
+    result += '\n\n'
+    result += genAteParam_BW6_opt(curve_name, curve_config)
+    return result
   else:
     raise ValueError(f'family: {family} is not implemented')
 
@@ -67,18 +72,130 @@ def genAteParam(curve_name, curve_config):
 
   return buf
 
+def genAteParam_BW6_unoptimized(curve_name, curve_config):
+  u = curve_config[curve_name]['field']['param']
+  family = curve_config[curve_name]['field']['family']
+  assert family == 'BW6'
+
+  # Algorithm 5 - https://eprint.iacr.org/2020/351.pdf
+  ate_param = u+1
+  ate_param_2 = u*(u^2 - u - 1)
+
+  ate_comment = '  # BW6-761 unoptimized Miller loop first part is parametrized by u+1\n'
+  ate_comment_2 = '  # BW6 unoptimized Miller loop second part is parametrized by u*(u²-u-1)\n'
+
+  # Note we can use the fact that
+  #  f_{u+1,Q}(P) = f_{u,Q}(P) . l_{[u]Q,Q}(P)
+  #  f_{u³-u²-u,Q}(P) = f_{u (u²-u-1),Q}(P)
+  #                   = (f_{u,Q}(P))^(u²-u-1) * f_{v,[u]Q}(P)
+  #
+  #  to have a common computation f_{u,Q}(P)
+  # but this require a scalar mul [u]Q
+  # and then its inversion to plug it back in the second Miller loop
+
+  # f_{u+1,Q}(P)
+  # ---------------------------------------------------------
+  buf = '# 1st part: f_{u+1,Q}(P)\n'
+  buf += f'const {curve_name}_pairing_ate_param_1_unopt* = block:\n'
+  buf += ate_comment
+
+  ate_bits = int(ate_param).bit_length()
+  naf_bits = int(3*ate_param).bit_length() - ate_bits
+
+  buf += f'  # +{naf_bits} to bitlength so that we can mul by 3 for NAF encoding\n'
+  buf += f'  BigInt[{ate_bits}+{naf_bits}].fromHex"0x{Integer(abs(ate_param)).hex()}"\n\n'
+
+  buf += f'const {curve_name}_pairing_ate_param_1_unopt_isNeg* = {"true" if ate_param < 0 else "false"}'
+
+  # frobenius(f_{u*(u²-u-1),Q}(P))
+  # ---------------------------------------------------------
+
+  buf += '\n\n\n'
+  buf += '# 2nd part: f_{u*(u²-u-1),Q}(P) followed by Frobenius application\n'
+  buf += f'const {curve_name}_pairing_ate_param_2_unopt* = block:\n'
+  buf += ate_comment_2
+
+  ate_2_bits = int(ate_param_2).bit_length()
+  naf_2_bits = int(3*ate_param_2).bit_length() - ate_2_bits
+
+  buf += f'  # +{naf_2_bits} to bitlength so that we can mul by 3 for NAF encoding\n'
+  buf += f'  BigInt[{ate_2_bits}+{naf_2_bits}].fromHex"0x{Integer(abs(ate_param_2)).hex()}"\n\n'
+
+  buf += f'const {curve_name}_pairing_ate_param_2_unopt_isNeg* = {"true" if ate_param_2 < 0 else "false"}'
+
+  buf += '\n'
+  return buf
+
+def genAteParam_BW6_opt(curve_name, curve_config):
+  u = curve_config[curve_name]['field']['param']
+  family = curve_config[curve_name]['field']['family']
+  assert family == 'BW6'
+
+  # Algorithm 5 - https://eprint.iacr.org/2020/351.pdf
+  ate_param = u
+  ate_param_2 = u^2 - u - 1
+
+  ate_comment = '  # BW6 Miller loop first part is parametrized by u\n'
+  ate_comment_2 = '  # BW6 Miller loop second part is parametrized by u²-u-1\n'
+
+  # Note we can use the fact that
+  #  f_{u+1,Q}(P) = f_{u,Q}(P) . l_{[u]Q,Q}(P)
+  #  f_{u³-u²-u,Q}(P) = f_{u (u²-u-1),Q}(P)
+  #                   = (f_{u,Q}(P))^(u²-u-1) * f_{v,[u]Q}(P)
+  #
+  #  to have a common computation f_{u,Q}(P)
+  # but this require a scalar mul [u]Q
+  # and then its inversion to plug it back in the second Miller loop
+
+  # f_{u,Q}(P)
+  # ---------------------------------------------------------
+  buf = '# 1st part: f_{u,Q}(P)\n'
+  buf += f'const {curve_name}_pairing_ate_param_1_opt* = block:\n'
+  buf += ate_comment
+
+  ate_bits = int(ate_param).bit_length()
+  naf_bits = 0 # int(3*ate_param).bit_length() - ate_bits
+
+  buf += f'  # no NAF for the optimized first Miller loop\n'
+  buf += f'  BigInt[{ate_bits}].fromHex"0x{Integer(abs(ate_param)).hex()}"\n\n'
+
+  buf += f'const {curve_name}_pairing_ate_param_1_opt_isNeg* = {"true" if ate_param < 0 else "false"}'
+
+  # frobenius(f_{u²-u-1,Q}(P))
+  # ---------------------------------------------------------
+
+  buf += '\n\n\n'
+  buf += '# 2nd part: f_{u²-u-1,Q}(P) followed by Frobenius application\n'
+  buf += f'const {curve_name}_pairing_ate_param_opt_2* = block:\n'
+  buf += ate_comment_2
+
+  ate_2_bits = int(ate_param_2).bit_length()
+  naf_2_bits = int(3*ate_param_2).bit_length() - ate_2_bits
+
+  buf += f'  # +{naf_2_bits} to bitlength so that we can mul by 3 for NAF encoding\n'
+  buf += f'  BigInt[{ate_2_bits}+{naf_2_bits}].fromHex"0x{Integer(abs(ate_param_2)).hex()}"\n\n'
+
+  buf += f'const {curve_name}_pairing_ate_param_2_opt_isNeg* = {"true" if ate_param_2 < 0 else "false"}'
+
+  buf += '\n'
+  return buf
+
 def genFinalExp(curve_name, curve_config):
   p = curve_config[curve_name]['field']['modulus']
   r = curve_config[curve_name]['field']['order']
   k = curve_config[curve_name]['tower']['embedding_degree']
   family = curve_config[curve_name]['field']['family']
 
+  # For BLS12 and BW6, 3*hard part has a better expression
+  # in the q basis with LLL algorithm
+  fexpMul3 = family == 'BLS12' or family == 'BW6'
+
   fexp = (p^k - 1)//r
-  if family == 'BLS12':
+  if fexpMul3:
     fexp *= 3
 
   buf = f'const {curve_name}_pairing_finalexponent* = block:\n'
-  buf += f'  # (p^{k} - 1) / r' + (' * 3' if family == 'BLS12' else '')
+  buf += f'  # (p^{k} - 1) / r' + (' * 3' if fexpMul3 else '')
   buf += '\n'
   buf += f'  BigInt[{int(fexp).bit_length()}].fromHex"0x{Integer(fexp).hex()}"'
 

tests/t_pairing_bls12_377_optate.nim

@@ -13,4 +13,9 @@ import
   # Test utilities
   ./t_pairing_template
 
-runPairingTests(4, BLS12_377, pairing_bls12)
+runPairingTests(
+  4, BLS12_377,
+  G1 = ECP_ShortW_Prj[Fp[BLS12_377], NotOnTwist],
+  G2 = ECP_ShortW_Prj[Fp2[BLS12_377], OnTwist],
+  GT = Fp12[BLS12_377],
+  pairing_bls12)

tests/t_pairing_bls12_381_optate.nim

@@ -13,5 +13,9 @@ import
   # Test utilities
   ./t_pairing_template
 
-# runPairingTests(4, BLS12_381, pairing_bls12_reference)
+runPairingTests(
-runPairingTests(4, BLS12_381, pairing_bls12)
+  4, BLS12_381,
+  G1 = ECP_ShortW_Prj[Fp[BLS12_381], NotOnTwist],
+  G2 = ECP_ShortW_Prj[Fp2[BLS12_381], OnTwist],
+  GT = Fp12[BLS12_381],
+  pairing_bls12)

tests/t_pairing_bn254_nogami_optate.nim

@@ -13,5 +13,9 @@ import
   # Test utilities
   ./t_pairing_template
 
-# runPairingTests(4, BN254_Nogami, pairing_bn_reference)
+runPairingTests(
-runPairingTests(4, BN254_Nogami, pairing_bn)
+  4, BN254_Nogami,
+  G1 = ECP_ShortW_Prj[Fp[BN254_Nogami], NotOnTwist],
+  G2 = ECP_ShortW_Prj[Fp2[BN254_Nogami], OnTwist],
+  GT = Fp12[BN254_Nogami],
+  pairing_bn)

tests/t_pairing_bn254_snarks_optate.nim

@@ -13,5 +13,9 @@ import
   # Test utilities
   ./t_pairing_template
 
-# runPairingTests(4, BN254_Snarks, pairing_bn_reference)
+runPairingTests(
-runPairingTests(4, BN254_Snarks, pairing_bn)
+  4, BN254_Snarks,
+  G1 = ECP_ShortW_Prj[Fp[BN254_Snarks], NotOnTwist],
+  G2 = ECP_ShortW_Prj[Fp2[BN254_Snarks], OnTwist],
+  GT = Fp12[BN254_Snarks],
+  pairing_bn)

tests/t_pairing_bw6_761_optate.nim

@@ -0,0 +1,21 @@
+# 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/common,
+  ../constantine/config/curves,
+  ../constantine/pairing/pairing_bw6_761,
+  # Test utilities
+  ./t_pairing_template
+
+runPairingTests(
+  4, BW6_761,
+  G1 = ECP_ShortW_Prj[Fp[BW6_761], NotOnTwist],
+  G2 = ECP_ShortW_Prj[Fp[BW6_761], OnTwist],
+  GT = Fp6[BW6_761],
+  pairing_bw6_761_reference)

tests/t_pairing_template.nim

@@ -52,7 +52,7 @@ func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen)
       result = rng.random_long01Seq_with_randZ(EC)
       result.clearCofactorReference()
 
-template runPairingTests*(Iters: static int, C: static Curve, pairing_fn: untyped): untyped {.dirty.}=
+template runPairingTests*(Iters: static int, C: static Curve, G1, G2, GT: typedesc, 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)
@@ -61,12 +61,12 @@ template runPairingTests*(Iters: static int, C: static Curve, pairing_fn: untype
 
   proc test_bilinearity_double_impl(randZ: bool, gen: RandomGen) =
     for _ in 0 ..< Iters:
-      let P = rng.random_point(ECP_ShortW_Prj[Fp[C], NotOnTwist], randZ, gen)
+      let P = rng.random_point(G1, randZ, gen)
-      let Q = rng.random_point(ECP_ShortW_Prj[Fp2[C], OnTwist], randZ, gen)
+      let Q = rng.random_point(G2, randZ, gen)
       var P2: typeof(P)
       var Q2: typeof(Q)
 
-      var r {.noInit.}, r2 {.noInit.}, r3 {.noInit.}: Fp12[C]
+      var r {.noInit.}, r2 {.noInit.}, r3 {.noInit.}: GT
 
       P2.double(P)
       Q2.double(Q)