From 258e7e516f3df9edd00eee8dc25dadea7abdd06e Mon Sep 17 00:00:00 2001 From: Mamy Ratsimbazafy Date: Sun, 7 Feb 2021 09:46:41 +0100 Subject: [PATCH] [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 | 57 ++++++++ constantine/curves/zoo_pairings.nim | 3 +- constantine/hash_to_curve/cofactors.nim | 21 ++- constantine/pairing/miller_loops.nim | 59 ++++++++ constantine/pairing/mul_fp6_by_lines.nim | 140 +++++++++++++++++++ constantine/pairing/pairing_bls12.nim | 59 ++------ constantine/pairing/pairing_bn.nim | 60 ++------- constantine/pairing/pairing_bw6_761.nim | 165 +++++++++++++++++++++++ sage/derive_pairing.sage | 121 ++++++++++++++++- tests/t_pairing_bls12_377_optate.nim | 7 +- tests/t_pairing_bls12_381_optate.nim | 8 +- tests/t_pairing_bn254_nogami_optate.nim | 8 +- tests/t_pairing_bn254_snarks_optate.nim | 8 +- tests/t_pairing_bw6_761_optate.nim | 21 +++ tests/t_pairing_template.nim | 8 +- 15 files changed, 627 insertions(+), 118 deletions(-) create mode 100644 constantine/curves/bw6_761_pairing.nim create mode 100644 constantine/pairing/miller_loops.nim create mode 100644 constantine/pairing/mul_fp6_by_lines.nim create mode 100644 constantine/pairing/pairing_bw6_761.nim create mode 100644 tests/t_pairing_bw6_761_optate.nim diff --git a/constantine/curves/bw6_761_pairing.nim b/constantine/curves/bw6_761_pairing.nim new file mode 100644 index 0000000..1516ae8 --- /dev/null +++ b/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 +# ------------------------------------------------------------ \ No newline at end of file diff --git a/constantine/curves/zoo_pairings.nim b/constantine/curves/zoo_pairings.nim index c25ac21..0f56c75 100644 --- a/constantine/curves/zoo_pairings.nim +++ b/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".} diff --git a/constantine/hash_to_curve/cofactors.nim b/constantine/hash_to_curve/cofactors.nim index cdf5965..eb5cce1 100644 --- a/constantine/hash_to_curve/cofactors.nim +++ b/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) diff --git a/constantine/pairing/miller_loops.nim b/constantine/pairing/miller_loops.nim new file mode 100644 index 0000000..0029734 --- /dev/null +++ b/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) diff --git a/constantine/pairing/mul_fp6_by_lines.nim b/constantine/pairing/mul_fp6_by_lines.nim new file mode 100644 index 0000000..5d012a8 --- /dev/null +++ b/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".} diff --git a/constantine/pairing/pairing_bls12.nim b/constantine/pairing/pairing_bls12.nim index 079243a..72dac8f 100644 --- a/constantine/pairing/pairing_bls12.nim +++ b/constantine/pairing/pairing_bls12.nim @@ -14,10 +14,10 @@ import ec_shortweierstrass_projective ], ../isogeny/frobenius, - ./lines_projective, - ./mul_fp12_by_lines, + ../curves/zoo_pairings, ./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) - var u3 = C.pairing(ate_param) - u3 *= 3 - for i in countdown(u3.bits - 2, 1): - 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() + basicMillerLoop( + f, T, line, + P, Q, nQ, + ate_param, ate_param_isNeg + ) func finalExpGeneric[C: static Curve](f: var Fp12[C]) = ## A generic and slow implementation of final exponentiation diff --git a/constantine/pairing/pairing_bn.nim b/constantine/pairing/pairing_bn.nim index 49bef54..f711e6f 100644 --- a/constantine/pairing/pairing_bn.nim +++ b/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, - ../curves/zoo_pairings + ./miller_loops # ############################################################ # @@ -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) - var u3 = C.pairing(ate_param) - u3 *= 3 - for i in countdown(u3.bits - 2, 1): - f.square() - line.line_double(T, P) - f.mul(line) + basicMillerLoop( + f, T, line, + P, Q, nQ, + ate_param, ate_param_isNeg + ) - 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() - - # Ate pairing for BN curves need adjustment after Miller loop + # Ate pairing for BN curves need adjustment after basic Miller loop when C.pairing(ate_param_isNeg): T.neg() var V {.noInit.}: typeof(Q) diff --git a/constantine/pairing/pairing_bw6_761.nim b/constantine/pairing/pairing_bw6_761.nim new file mode 100644 index 0000000..c0d2625 --- /dev/null +++ b/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() diff --git a/sage/derive_pairing.sage b/sage/derive_pairing.sage index e86149b..adf9c3d 100644 --- a/sage/derive_pairing.sage +++ b/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()}"' diff --git a/tests/t_pairing_bls12_377_optate.nim b/tests/t_pairing_bls12_377_optate.nim index f963b79..4f367e5 100644 --- a/tests/t_pairing_bls12_377_optate.nim +++ b/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) diff --git a/tests/t_pairing_bls12_381_optate.nim b/tests/t_pairing_bls12_381_optate.nim index 522432e..5c117b5 100644 --- a/tests/t_pairing_bls12_381_optate.nim +++ b/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(4, BLS12_381, pairing_bls12) +runPairingTests( + 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) diff --git a/tests/t_pairing_bn254_nogami_optate.nim b/tests/t_pairing_bn254_nogami_optate.nim index 20f59be..4be9a24 100644 --- a/tests/t_pairing_bn254_nogami_optate.nim +++ b/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(4, BN254_Nogami, pairing_bn) +runPairingTests( + 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) diff --git a/tests/t_pairing_bn254_snarks_optate.nim b/tests/t_pairing_bn254_snarks_optate.nim index 3554b03..46a5233 100644 --- a/tests/t_pairing_bn254_snarks_optate.nim +++ b/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(4, BN254_Snarks, pairing_bn) +runPairingTests( + 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) diff --git a/tests/t_pairing_bw6_761_optate.nim b/tests/t_pairing_bw6_761_optate.nim new file mode 100644 index 0000000..954de38 --- /dev/null +++ b/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) diff --git a/tests/t_pairing_template.nim b/tests/t_pairing_template.nim index 0dc5c0d..e4a1c51 100644 --- a/tests/t_pairing_template.nim +++ b/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 Q = rng.random_point(ECP_ShortW_Prj[Fp2[C], OnTwist], randZ, gen) + let P = rng.random_point(G1, 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)