From a022db1c08a435fc386d0a505202e1c66ed14c86 Mon Sep 17 00:00:00 2001 From: Mamy Ratsimbazafy Date: Sun, 13 Dec 2020 17:08:32 +0100 Subject: [PATCH] Sqrt fp2 acceleration (#109) * Use sqrt.square() == a instead of sqrt * invsqrt = -1 (Euler criterion) for sqrt existence. * Accelerate sqrt_fp2 by 33% --- .../arithmetic/finite_fields_square_root.nim | 6 +- ...377_square_root.nim => bls12_377_sqrt.nim} | 0 constantine/curves/bls12_381_sqrt_fp2.nim | 26 ++ constantine/curves/bn254_nogami_sqrt_fp2.nim | 26 ++ constantine/curves/bn254_snarks_sqrt_fp2.nim | 26 ++ constantine/curves/zoo_square_roots.nim | 2 +- constantine/curves/zoo_square_roots_fp2.nim | 19 ++ constantine/io/io_towers.nim | 4 +- .../exponentiations.nim | 85 ------ .../square_root_fp2.nim | 205 ++++++++++++++ .../tower_instantiation.nim | 259 ++++++++++++++++++ constantine/towers.nim | 250 +---------------- sage/derive_square_root.sage | 151 ++++++++++ 13 files changed, 723 insertions(+), 336 deletions(-) rename constantine/curves/{bls12_377_square_root.nim => bls12_377_sqrt.nim} (100%) create mode 100644 constantine/curves/bls12_381_sqrt_fp2.nim create mode 100644 constantine/curves/bn254_nogami_sqrt_fp2.nim create mode 100644 constantine/curves/bn254_snarks_sqrt_fp2.nim create mode 100644 constantine/curves/zoo_square_roots_fp2.nim create mode 100644 constantine/tower_field_extensions/square_root_fp2.nim create mode 100644 constantine/tower_field_extensions/tower_instantiation.nim create mode 100644 sage/derive_square_root.sage diff --git a/constantine/arithmetic/finite_fields_square_root.nim b/constantine/arithmetic/finite_fields_square_root.nim index 7a23cf3..bc2dcfe 100644 --- a/constantine/arithmetic/finite_fields_square_root.nim +++ b/constantine/arithmetic/finite_fields_square_root.nim @@ -89,9 +89,9 @@ func sqrt_invsqrt_if_square_p3mod4[C](sqrt, invsqrt: var Fp[C], a: Fp[C]): Secre ## ## This assumes that the prime field modulus ``p``: p ≡ 3 (mod 4) sqrt_invsqrt_p3mod4(sqrt, invsqrt, a) - var euler {.noInit.}: Fp[C] - euler.prod(sqrt, invsqrt) - result = not(euler.isMinusOne()) + var test {.noInit.}: Fp[C] + test.square(sqrt) + result = test == a func sqrt_if_square_p3mod4[C](a: var Fp[C]): SecretBool {.inline.} = ## If ``a`` is a square, compute the square root of ``a`` diff --git a/constantine/curves/bls12_377_square_root.nim b/constantine/curves/bls12_377_sqrt.nim similarity index 100% rename from constantine/curves/bls12_377_square_root.nim rename to constantine/curves/bls12_377_sqrt.nim diff --git a/constantine/curves/bls12_381_sqrt_fp2.nim b/constantine/curves/bls12_381_sqrt_fp2.nim new file mode 100644 index 0000000..ceb52ea --- /dev/null +++ b/constantine/curves/bls12_381_sqrt_fp2.nim @@ -0,0 +1,26 @@ +# 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, + ../io/io_towers + +# Square Root Fp2 constants +# ----------------------------------------------------------------- +const BLS12_381_sqrt_fp2_QNR* = Fp2[BLS12_381].fromHex( + "0x0", + "0x1" +) +const BLS12_381_sqrt_fp2_sqrt_QNR* = Fp2[BLS12_381].fromHex( + "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09", + "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2" +) +const BLS12_381_sqrt_fp2_minus_sqrt_QNR* = Fp2[BLS12_381].fromHex( + "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2", + "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2" +) diff --git a/constantine/curves/bn254_nogami_sqrt_fp2.nim b/constantine/curves/bn254_nogami_sqrt_fp2.nim new file mode 100644 index 0000000..37e995d --- /dev/null +++ b/constantine/curves/bn254_nogami_sqrt_fp2.nim @@ -0,0 +1,26 @@ +# 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, + ../io/io_towers + +# Square Root Fp2 constants +# ----------------------------------------------------------------- +const BN254_Nogami_sqrt_fp2_QNR* = Fp2[BN254_Nogami].fromHex( + "0x0", + "0x1" +) +const BN254_Nogami_sqrt_fp2_sqrt_QNR* = Fp2[BN254_Nogami].fromHex( + "0x1439ab09c60b248f398c5d77b755f92b9edc5f19d2873545be471151a747e4e", + "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5" +) +const BN254_Nogami_sqrt_fp2_minus_sqrt_QNR* = Fp2[BN254_Nogami].fromHex( + "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5", + "0x23dfc9d1a39f4db8c69b87a8848aa075a7333a0e62d78cbf4b1b8eeae58b81c5" +) diff --git a/constantine/curves/bn254_snarks_sqrt_fp2.nim b/constantine/curves/bn254_snarks_sqrt_fp2.nim new file mode 100644 index 0000000..70e0c21 --- /dev/null +++ b/constantine/curves/bn254_snarks_sqrt_fp2.nim @@ -0,0 +1,26 @@ +# 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, + ../io/io_towers + +# Square Root Fp2 constants +# ----------------------------------------------------------------- +const BN254_Snarks_sqrt_fp2_QNR* = Fp2[BN254_Snarks].fromHex( + "0x0", + "0x1" +) +const BN254_Snarks_sqrt_fp2_sqrt_QNR* = Fp2[BN254_Snarks].fromHex( + "0x4636956ffd65e421c784f990c3a7533717e614fc6e7f616577d10f6464b5204", + "0x4636956ffd65e421c784f990c3a7533717e614fc6e7f616577d10f6464b5204" +) +const BN254_Snarks_sqrt_fp2_minus_sqrt_QNR* = Fp2[BN254_Snarks].fromHex( + "0x2c00e51be15b41e79bd7f61d7546e32a26030941a189d476e4a37b209231ab43", + "0x4636956ffd65e421c784f990c3a7533717e614fc6e7f616577d10f6464b5204" +) diff --git a/constantine/curves/zoo_square_roots.nim b/constantine/curves/zoo_square_roots.nim index 1b52b53..be937b3 100644 --- a/constantine/curves/zoo_square_roots.nim +++ b/constantine/curves/zoo_square_roots.nim @@ -9,7 +9,7 @@ import std/macros, ../config/curves, - ./bls12_377_square_root + ./bls12_377_sqrt {.experimental: "dynamicBindSym".} diff --git a/constantine/curves/zoo_square_roots_fp2.nim b/constantine/curves/zoo_square_roots_fp2.nim new file mode 100644 index 0000000..2835655 --- /dev/null +++ b/constantine/curves/zoo_square_roots_fp2.nim @@ -0,0 +1,19 @@ +# Constantine +# Copyright (c) 2018-2019 Status Research & Development GmbH +# Copyright (c) 2020-Present Mamy André-Ratsimbazafy +# Licensed and distributed under either of +# * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT). +# * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0). +# at your option. This file may not be copied, modified, or distributed except according to those terms. + +import + std/macros, + ../config/curves, + ./bls12_381_sqrt_fp2, + ./bn254_nogami_sqrt_fp2, + ./bn254_snarks_sqrt_fp2 + +{.experimental: "dynamicBindSym".} +macro sqrt_fp2*(C: static Curve, value: untyped): untyped = + ## Get square_root constants + return bindSym($C & "_sqrt_fp2_" & $value) diff --git a/constantine/io/io_towers.nim b/constantine/io/io_towers.nim index 4dae94f..fcf2034 100644 --- a/constantine/io/io_towers.nim +++ b/constantine/io/io_towers.nim @@ -13,7 +13,9 @@ import # Internal ./io_bigints, ./io_fields, ../arithmetic/finite_fields, - ../towers + ../tower_field_extensions/tower_instantiation + +export tower_instantiation # No exceptions allowed {.push raises: [].} diff --git a/constantine/tower_field_extensions/exponentiations.nim b/constantine/tower_field_extensions/exponentiations.nim index f09e6ab..a5e18b1 100644 --- a/constantine/tower_field_extensions/exponentiations.nim +++ b/constantine/tower_field_extensions/exponentiations.nim @@ -24,11 +24,6 @@ import # Square root should be implemented in constant-time for hash-to-curve: # https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-05#section-4 -# -# Further non-constant-time optimization may be used -# - Square Root Computation over Even Extension Fields -# Gora Adj, Francisco Rodríguez-Henríquez, 2012 -# https://eprint.iacr.org/2012/685 # No exceptions allowed {.push raises: [].} @@ -194,83 +189,3 @@ func powUnsafeExponent*[F; bits: static int]( var expBE {.noInit.}: array[(bits + 7) div 8, byte] expBE.exportRawUint(exponent, bigEndian) a.powUnsafeExponent(expBE, window) - -# Square root -# ----------------------------------------------------------- -# -# Warning ⚠️: -# p the characteristic, i.e. the prime modulus of the base field -# in extension field we require q = p^m be of special form -# i.e. q ≡ 3 (mod 4) or q ≡ 9 (mod 16) -# -# In Fp2 in particular p² ≡ 1 (mod 4) always hold -# and p² ≡ 5 (mod 8) is not possible -# if Fp2 = Fp[v]/(v² − β) with β a quadratic non-residue in Fp - -func isSquare*(a: QuadraticExt): SecretBool = - ## Returns true if ``a`` is a square (quadratic residue) in 𝔽p2 - ## - ## Assumes that the prime modulus ``p`` is public. - # Implementation: - # - # (a0, a1) = a in F(p^2) - # is_square(a) = is_square(|a|) over F(p) - # where |a| = a0^2 + a1^2 - # - # This can be done recursively in an extension tower - # - # https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-08#appendix-G.5 - # https://eprint.iacr.org/2012/685 - mixin fromComplexExtension - - var tv1{.noInit.}, tv2{.noInit.}: typeof(a.c0) - - tv1.square(a.c0) # a0² - tv2.square(a.c1) # - β a1² with β = 𝑖² in a complex extension field - when a.fromComplexExtension(): - tv1 += tv2 # a0 - (-1) a1² - else: - tv2 *= NonResidue - tv1 -= tv2 - - result = tv1.isSquare() - -func sqrt_if_square*(a: var QuadraticExt): SecretBool = - ## If ``a`` is a square, compute the square root of ``a`` - ## if not, ``a`` is unmodified. - ## - ## The square root, if it exist is multivalued, - ## i.e. both x² == (-x)² - ## This procedure returns a deterministic result - # - # Implementation via the complex method (which confusingly does not require a complex field) - # We make it constant-time via conditional copies - mixin fromComplexExtension - - var t1{.noInit.}, t2{.noInit.}, t3{.noInit.}: typeof(a.c0) - - t1.square(a.c0) # a0² - t2.square(a.c1) # - β a1² with β = 𝑖² in a complex extension field - when a.fromComplexExtension(): - t1 += t2 # a0 - (-1) a1² - else: - t2 *= NonResidue - t1 -= t2 - - result = t1.sqrt_if_square() - - t2.sum(a.c0, t1) - t2.div2() - - t3.diff(a.c0, t1) - t3.div2() - - let quadResidTest = t2.isSquare() - t2.ccopy(t3, not quadResidTest) - - sqrt_invsqrt(sqrt = t1, invsqrt = t3, t2) - a.c0.ccopy(t1, result) - - t3.div2() - t3 *= a.c1 - a.c1.ccopy(t3, result) diff --git a/constantine/tower_field_extensions/square_root_fp2.nim b/constantine/tower_field_extensions/square_root_fp2.nim new file mode 100644 index 0000000..7180cbd --- /dev/null +++ b/constantine/tower_field_extensions/square_root_fp2.nim @@ -0,0 +1,205 @@ +# 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 + ./tower_instantiation, + ../arithmetic, + ../primitives, + ../config/[common, curves], + ../curves/zoo_square_roots_fp2 + +# Square root +# ----------------------------------------------------------- + +func isSquare*(a: Fp2): SecretBool = + ## Returns true if ``a`` is a square (quadratic residue) in 𝔽p2 + ## + ## Assumes that the prime modulus ``p`` is public. + # Implementation: + # + # (a0, a1) = a in F(p^2) + # is_square(a) = is_square(|a|) over F(p) + # where |a| = a0^2 + a1^2 + # + # This can be done recursively in an extension tower + # + # https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-08#appendix-G.5 + # https://eprint.iacr.org/2012/685 + var tv1{.noInit.}, tv2{.noInit.}: typeof(a.c0) + + tv1.square(a.c0) # a0² + tv2.square(a.c1) # - β a1² with β = 𝑖² in a complex extension field + when a.fromComplexExtension(): + tv1 += tv2 # a0 - (-1) a1² + else: + tv2 *= NonResidue + tv1 -= tv2 + + result = tv1.isSquare() + +func sqrt_rotate_extension( + out_sqrt: var QuadraticExt, + candidate_sqrt: QuadraticExt, + a: QuadraticExt + ): SecretBool = + ## From a field element `a` and a candidate Fp2 square root + ## Search the actual square root by rotating candidate solution + ## in the extension field by 90° + ## + ## if there is one, update out_sqrt with it and return true + ## return false otherwise, out_sqrt is undefined in this case + ## + ## This avoids expensive trial "isSquare" checks (450+ field multiplications) + ## This requires the sqrt of sqrt of the quadratic non-residue + ## to be in Fp2 + var coeff{.noInit.}, cand2{.noInit.}, t{.noInit.}: QuadraticExt + const Curve = typeof(a.c0).C + + # We name µ² the quadratic non-residue + # if p ≡ 3 (mod 4), we have µ = 𝑖 = √-1 and µ² = -1 + # However for BLS12-377 we have µ = √-5 + + # sqrt(cand)² = (a0 + µ a1)² = (a0²-a1²) + (2 a0a1) µ + cand2.square(candidate_sqrt) + + block: # Test 1: (a0²-a1²) + (2 a0a1) µ == cand ? candidate is correct + t.diff(cand2, a) + result = t.isZero() + coeff.setOne() + + block: # Test 2: -((a0²-a1²) + (2 a0a1) µ) == candidate ? candidate must be rotated by 90° + t.sum(cand2, a) + let isSol = t.isZero() + coeff.ccopy(Curve.sqrt_fp2(QNR), isSol) + result = result or isSol + + block: # Test 3: µ((a0²-a1²) + (2 a0a1) µ) == candidate ? candidate must be rotated by 135° + t.c0.diff(cand2.c0, a.c1) + t.c1.sum( cand2.c1, a.c0) + let isSol = t.isZero() + coeff.ccopy(Curve.sqrt_fp2(sqrt_QNR), isSol) + result = result or isSol + + block: # Test 4: -µ((a0²-a1²) + (2 a0a1) µ) == candidate ? candidate must be rotated by 45° + t.c0.sum( cand2.c0, a.c1) + t.c1.diff(cand2.c1, a.c0) + let isSol = t.isZero() + coeff.ccopy(Curve.sqrt_fp2(minus_sqrt_QNR), isSol) + result = result or isSol + + # Rotate the candidate + out_sqrt.prod(candidate_sqrt, coeff) + # result is set + +func sqrt_if_square_opt(a: var Fp2): SecretBool = + ## If ``a`` is a square, compute the square root of ``a`` + ## if not, ``a`` is unmodified. + ## + ## The square root, if it exist is multivalued, + ## i.e. both x² == (-x)² + ## This procedure returns a deterministic result + ## + ## This is an optimized version which is + ## requires the sqrt of sqrt of the quadratic non-residue + ## to be in Fp2 + # + # Implementation via the complex method + # Gora Adj, Francisco Rodríguez-Henríquez, 2012, https://eprint.iacr.org/2012/685 + # Made constant-time and optimized to fuse sqrt and inverse sqrt + # and avoid unfused isSquare tests. + # See discussion and optimization with Andy Polyakov + # https://github.com/supranational/blst/issues/2#issuecomment-686656784 + var t1{.noInit.}, t2{.noInit.}: typeof(a.c0) + var cand{.noInit.}: typeof(a) + + t1.square(a.c0) # a0² + t2.square(a.c1) # - β a1² with β = 𝑖² in a complex extension field + when a.fromComplexExtension(): + t1 += t2 # a0 - (-1) a1² + else: + t2 *= NonResidue + t1 -= t2 + + # t1 being an actual sqrt will be tested in sqrt_rotate_extension + t1.sqrt() # sqrt(a0² - β a1²) + + t2.diff(a.c0, t1) + t1 += a.c0 + t1.ccopy(t2, t1.isZero()) + t1.div2() # (a0 ± sqrt(a0² - β a1²))/2 + + # t1 being an actual sqrt will be tested in sqrt_rotate_extension + # 1/sqrt((a0 ± sqrt(a0² - β b²))/2) + sqrt_invsqrt(sqrt = cand.c1, invsqrt = cand.c0, t1) # we only want invsqrt = cand.c0 + + cand.c1 = a.c1 + cand.c1.div2() + cand.c1 *= cand.c0 # a1/(2*sqrt((a0 ± sqrt(a0² - β a1²))/2)) + cand.c0 *= t1 # sqrt((a0 ± sqrt(a0² - β a1²))/2) + + # Now rotate to check if an actual sqrt exists. + return sqrt_rotate_extension(a, cand, a) + +func sqrt_if_square_generic(a: var Fp2): SecretBool = + ## If ``a`` is a square, compute the square root of ``a`` + ## if not, ``a`` is unmodified. + ## + ## The square root, if it exist is multivalued, + ## i.e. both x² == (-x)² + ## This procedure returns a deterministic result + ## + ## This is a generic version + # Implementation via the complex method + # Gora Adj, Francisco Rodríguez-Henríquez, 2012, + # https://eprint.iacr.org/2012/685 + # Made constant-time and optimized to fuse sqrt and inverse sqrt + var t1{.noInit.}, t2{.noInit.}, t3{.noInit.}: typeof(a.c0) + + t1.square(a.c0) # a0² + t2.square(a.c1) # - β a1² with β = 𝑖² in a complex extension field + when a.fromComplexExtension(): + t1 += t2 # a0 - (-1) a1² + else: + t2 *= NonResidue + t1 -= t2 + + result = t1.sqrt_if_square() + + t2.sum(a.c0, t1) + t2.div2() + + t3.diff(a.c0, t1) + t3.div2() + + let quadResidTest = t2.isSquare() + t2.ccopy(t3, not quadResidTest) + + sqrt_invsqrt(sqrt = t1, invsqrt = t3, t2) + a.c0.ccopy(t1, result) + + t3.div2() + t3 *= a.c1 + a.c1.ccopy(t3, result) + +func sqrt_if_square*(a: var Fp2): SecretBool = + ## If ``a`` is a square, compute the square root of ``a`` + ## if not, ``a`` is unmodified. + ## + ## The square root, if it exist is multivalued, + ## i.e. both x² == (-x)² + when Fp2.C == BLS12_377: + # For BLS12_377, + # the solution µ to x² - µ = 0 being a quadratic non-residue + # is also a quadratic non-residue in Fp2, which means + # we can't use the optimized version which saves an `isSquare` + # which is about 33% of processing time + # as isSquare, sqrt and invsqrt + # all requires over 450 Fp multiplications. + result = a.sqrt_if_square_generic() + else: + result = a.sqrt_if_square_opt() diff --git a/constantine/tower_field_extensions/tower_instantiation.nim b/constantine/tower_field_extensions/tower_instantiation.nim new file mode 100644 index 0000000..39c688a --- /dev/null +++ b/constantine/tower_field_extensions/tower_instantiation.nim @@ -0,0 +1,259 @@ +# 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. + +# Instantiate the actual tower extensions +# that were described in a "conceptualized" way +# ---------------------------------------------------------------- + +import + ../arithmetic, + ../config/curves, + ../io/io_fields, + ./tower_common, + ./quadratic_extensions, + ./cubic_extensions, + ./exponentiations + +export tower_common, quadratic_extensions, cubic_extensions, exponentiations + +# 𝔽p2 +# ---------------------------------------------------------------- + +type + Fp2*[C: static Curve] = object + c0*, c1*: Fp[C] + + β = NonResidue + # Quadratic or Cubic non-residue + + SexticNonResidue* = object + +template fromComplexExtension*[F](elem: F): static bool = + ## Returns true if the input is a complex extension + ## i.e. the irreducible polynomial chosen is + ## x² - µ with µ = -1 + ## and so 𝔽p2 = 𝔽p[x] / x² - µ = 𝔽p[𝑖] + when F is Fp2 and F.C.get_QNR_Fp() == -1: + true + else: + false + +func `*=`*(a: var Fp, _: typedesc[β]) {.inline.} = + ## Multiply an element of 𝔽p by the quadratic non-residue + ## chosen to construct 𝔽p2 + static: doAssert Fp.C.get_QNR_Fp() != -1, "𝔽p2 should be specialized for complex extension" + a *= Fp.C.get_QNR_Fp() + +func `*`*(_: typedesc[β], a: Fp): Fp {.inline, noInit.} = + ## Multiply an element of 𝔽p by the quadratic non-residue + ## chosen to construct 𝔽p2 + result = a + result *= β + +# TODO: rework the quad/cube/sextic non residue declaration + +func `*=`*(a: var Fp, _: typedesc[SexticNonResidue]) {.inline.} = + ## Multiply an element of 𝔽p by the sextic non-residue + ## chosen to construct 𝔽p6 + a *= Fp.C.get_QNR_Fp() + +func `*`*(_: typedesc[SexticNonResidue], a: Fp): Fp {.inline, noInit.} = + ## Multiply an element of 𝔽p by the sextic non-residue + ## chosen to construct 𝔽p6 + result = a + result *= SexticNonResidue + +func `*=`*(a: var Fp2, _: typedesc[SexticNonResidue]) {.inline.} = + ## Multiply an element of 𝔽p2 by the sextic non-residue + ## chosen to construct the sextic twist + # Yet another const tuple unpacking bug + const u = Fp2.C.get_SNR_Fp2()[0] # Sextic non-residue to construct 𝔽p12 + const v = Fp2.C.get_SNR_Fp2()[1] + const Beta {.used.} = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 + # ξ = u + v x + # and x² = β + # + # (c0 + c1 x) (u + v x) => u c0 + (u c0 + u c1)x + v c1 x² + # => u c0 + β v c1 + (v c0 + u c1) x + when a.fromComplexExtension() and u == 1 and v == 1: + let t = a.c0 + a.c0 -= a.c1 + a.c1 += t + else: + let a0 = a.c0 + let a1 = a.c1 + when a.fromComplexExtension(): + a.c0.diff(u * a0, v * a1) + else: + a.c0.sum(u * a0, (Beta * v) * a1) + a.c1.sum(v * a0, u * a1) + +func `/=`*(a: var Fp2, _: typedesc[SexticNonResidue]) {.inline.} = + ## Multiply an element of 𝔽p by the quadratic non-residue + ## chosen to construct sextic twist + # Yet another const tuple unpacking bug + const u = Fp2.C.get_SNR_Fp2()[0] # Sextic non-residue to construct 𝔽p12 + const v = Fp2.C.get_SNR_Fp2()[1] + const Beta = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 + # ξ = u + v x + # and x² = β + # + # (c0 + c1 x) / (u + v x) => (c0 + c1 x)(u - v x) / ((u + vx)(u-vx)) + # => u c0 - v c1 x² + (u c1 - v c0) x / (u² - x²v²) + # => 1/(u² - βv²) * (uc0 - β v c1, u c1 - v c0) + # With β = 𝑖 = √-1 + # 1/(u² + v²) * (u c0 + v c1, u c1 - v c0) + # + # With β = 𝑖 = √-1 and ξ = 1 + 𝑖 + # 1/2 * (c0 + c1, c1 - c0) + + when a.fromComplexExtension() and u == 1 and v == 1: + let t = a.c0 + a.c0 += a.c1 + a.c1 -= t + a.div2() + else: + let a0 = a.c0 + let a1 = a.c1 + const u2v2 = u*u - Beta*v*v # (u² - βv²) + # TODO can be precomputed (or precompute b/µ the twist coefficient) + # and use faster non-constant-time inversion in the VM + var u2v2inv {.noInit.}: a.c0.typeof + u2v2inv.fromUint(u2v2) + u2v2inv.inv() + + a.c0.diff(u * a0, (Beta * v) * a1) + a.c1.diff(u * a1, v * a0) + a.c0 *= u2v2inv + a.c1 *= u2v2inv + +# 𝔽p6 +# ---------------------------------------------------------------- + +type + Fp4*[C: static Curve] = object + c0*, c1*: Fp2[C] + + Fp6*[C: static Curve] = object + c0*, c1*, c2*: Fp2[C] + + ξ* = NonResidue + # We call the non-residue ξ on 𝔽p4/𝔽p6 to avoid confusion + # between non-residue + # of different tower level + +func `*`*(_: typedesc[ξ], a: Fp2): Fp2 {.inline, noInit.} = + ## Multiply an element of 𝔽p2 by the quadratic and cubic non-residue + ## chosen to construct 𝔽p4/𝔽p6 + # Yet another const tuple unpacking bug + const u = Fp2.C.get_CNR_Fp2()[0] # Quadratic & Cubic non-residue to construct 𝔽p4/𝔽p6 + const v = Fp2.C.get_CNR_Fp2()[1] + const Beta {.used.} = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 + # ξ = u + v x + # and x² = β + # + # (c0 + c1 x) (u + v x) => u c0 + (u c0 + u c1)x + v c1 x² + # => u c0 + β v c1 + (v c0 + u c1) x + + # TODO: check code generated when ξ = 1 + 𝑖 + # The mul by constant are inline but + # since we don't have __restrict tag + # and we use arrays (which decay into pointer) + # the compiler might not elide the temporary + when a.fromComplexExtension(): + result.c0.diff(u * a.c0, v * a.c1) + else: + result.c0.sum(u * a.c0, (Beta * v) * a.c1) + result.c1.sum(v * a.c0, u * a.c1 ) + +func `*=`*(a: var Fp2, _: typedesc[ξ]) {.inline.} = + ## Multiply an element of 𝔽p2 by the quadratic non-residue + ## chosen to construct 𝔽p6 + # Yet another const tuple unpacking bug + const u = Fp2.C.get_CNR_Fp2()[0] # Cubic non-residue to construct 𝔽p6 + const v = Fp2.C.get_CNR_Fp2()[1] + const Beta {.used.} = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 + # ξ = u + v x + # and x² = β + # + # (c0 + c1 x) (u + v x) => u c0 + (u c0 + u c1)x + v c1 x² + # => u c0 + β v c1 + (v c0 + u c1) x + when a.fromComplexExtension() and u == 1 and v == 1: + let t = a.c0 + a.c0 -= a.c1 + a.c1 += t + else: # TODO: optim for inline + a = ξ * a + +# 𝔽p12 +# ---------------------------------------------------------------- + +type + Fp12*[C: static Curve] = object + c0*, c1*, c2*: Fp4[C] + + γ = NonResidue + # We call the non-residue γ (Gamma) on 𝔽p6 to avoid confusion between non-residue + # of different tower level + +func `*`*(_: typedesc[γ], a: Fp4): Fp4 {.noInit, inline.} = + ## Multiply an element of 𝔽p6 by the cubic non-residue + ## chosen to construct 𝔽p12 + ## For all curves γ = v with v the factor for 𝔽p6 coordinate + ## and v³ = ξ + ## (c0 + c1 v + c2 v²) v => ξ c2 + c0 v + c1 v² + result.c0 = ξ * a.c1 + result.c1 = a.c0 + +func `*=`*(a: var Fp4, _: typedesc[γ]) {.inline.} = + a = γ * a + +# Sparse functions +# ---------------------------------------------------------------- + +func `*=`*(a: var Fp2, b: Fp) = + ## Multiply an element of Fp2 by an element of Fp + a.c0 *= b + a.c1 *= b + +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₁) + r.c0.prod(a.c1, b) + r.c0 *= NonResidue + r.c1.prod(a.c0, b) + +func mul_sparse_by_0y0*[C: static Curve](r: var Fp6[C], a: Fp6[C], b: Fp2[C]) = + ## Sparse multiplication of an Fp6 element + ## with coordinates (a₀, a₁, a₂) by (0, b₁, 0) + # TODO: make generic and move to tower_field_extensions + + # v0 = a0 b0 = 0 + # v1 = a1 b1 + # v2 = a2 b2 = 0 + # + # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 + # = ξ (a1 b1 + a2 b1 - v1) + # = ξ a2 b1 + # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 + # = a0 b1 + a1 b1 - v1 + # = a0 b1 + # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 + # = v1 + # = a1 b1 + + r.c0.prod(a.c2, b) + r.c0 *= ξ + r.c1.prod(a.c0, b) + r.c2.prod(a.c1, b) diff --git a/constantine/towers.nim b/constantine/towers.nim index 78e1355..993c1c0 100644 --- a/constantine/towers.nim +++ b/constantine/towers.nim @@ -7,251 +7,9 @@ # at your option. This file may not be copied, modified, or distributed except according to those terms. import - ./arithmetic, - ./config/curves, - ./io/io_fields, - ./tower_field_extensions/[ - tower_common, - quadratic_extensions, - cubic_extensions, - exponentiations + tower_field_extensions/[ + tower_instantiation, + square_root_fp2 ] -export tower_common, quadratic_extensions, cubic_extensions, exponentiations - -# 𝔽p2 -# ---------------------------------------------------------------- - -type - Fp2*[C: static Curve] = object - c0*, c1*: Fp[C] - - β = NonResidue - # Quadratic or Cubic non-residue - - SexticNonResidue* = object - -template fromComplexExtension*[F](elem: F): static bool = - ## Returns true if the input is a complex extension - ## i.e. the irreducible polynomial chosen is - ## x² - µ with µ = -1 - ## and so 𝔽p2 = 𝔽p[x] / x² - µ = 𝔽p[𝑖] - when F is Fp2 and F.C.get_QNR_Fp() == -1: - true - else: - false - -func `*=`*(a: var Fp, _: typedesc[β]) {.inline.} = - ## Multiply an element of 𝔽p by the quadratic non-residue - ## chosen to construct 𝔽p2 - static: doAssert Fp.C.get_QNR_Fp() != -1, "𝔽p2 should be specialized for complex extension" - a *= Fp.C.get_QNR_Fp() - -func `*`*(_: typedesc[β], a: Fp): Fp {.inline, noInit.} = - ## Multiply an element of 𝔽p by the quadratic non-residue - ## chosen to construct 𝔽p2 - result = a - result *= β - -# TODO: rework the quad/cube/sextic non residue declaration - -func `*=`*(a: var Fp, _: typedesc[SexticNonResidue]) {.inline.} = - ## Multiply an element of 𝔽p by the sextic non-residue - ## chosen to construct 𝔽p6 - a *= Fp.C.get_QNR_Fp() - -func `*`*(_: typedesc[SexticNonResidue], a: Fp): Fp {.inline, noInit.} = - ## Multiply an element of 𝔽p by the sextic non-residue - ## chosen to construct 𝔽p6 - result = a - result *= SexticNonResidue - -func `*=`*(a: var Fp2, _: typedesc[SexticNonResidue]) {.inline.} = - ## Multiply an element of 𝔽p2 by the sextic non-residue - ## chosen to construct the sextic twist - # Yet another const tuple unpacking bug - const u = Fp2.C.get_SNR_Fp2()[0] # Sextic non-residue to construct 𝔽p12 - const v = Fp2.C.get_SNR_Fp2()[1] - const Beta {.used.} = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 - # ξ = u + v x - # and x² = β - # - # (c0 + c1 x) (u + v x) => u c0 + (u c0 + u c1)x + v c1 x² - # => u c0 + β v c1 + (v c0 + u c1) x - when a.fromComplexExtension() and u == 1 and v == 1: - let t = a.c0 - a.c0 -= a.c1 - a.c1 += t - else: - let a0 = a.c0 - let a1 = a.c1 - when a.fromComplexExtension(): - a.c0.diff(u * a0, v * a1) - else: - a.c0.sum(u * a0, (Beta * v) * a1) - a.c1.sum(v * a0, u * a1) - -func `/=`*(a: var Fp2, _: typedesc[SexticNonResidue]) {.inline.} = - ## Multiply an element of 𝔽p by the quadratic non-residue - ## chosen to construct sextic twist - # Yet another const tuple unpacking bug - const u = Fp2.C.get_SNR_Fp2()[0] # Sextic non-residue to construct 𝔽p12 - const v = Fp2.C.get_SNR_Fp2()[1] - const Beta = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 - # ξ = u + v x - # and x² = β - # - # (c0 + c1 x) / (u + v x) => (c0 + c1 x)(u - v x) / ((u + vx)(u-vx)) - # => u c0 - v c1 x² + (u c1 - v c0) x / (u² - x²v²) - # => 1/(u² - βv²) * (uc0 - β v c1, u c1 - v c0) - # With β = 𝑖 = √-1 - # 1/(u² + v²) * (u c0 + v c1, u c1 - v c0) - # - # With β = 𝑖 = √-1 and ξ = 1 + 𝑖 - # 1/2 * (c0 + c1, c1 - c0) - - when a.fromComplexExtension() and u == 1 and v == 1: - let t = a.c0 - a.c0 += a.c1 - a.c1 -= t - a.div2() - else: - let a0 = a.c0 - let a1 = a.c1 - const u2v2 = u*u - Beta*v*v # (u² - βv²) - # TODO can be precomputed (or precompute b/µ the twist coefficient) - # and use faster non-constant-time inversion in the VM - var u2v2inv {.noInit.}: a.c0.typeof - u2v2inv.fromUint(u2v2) - u2v2inv.inv() - - a.c0.diff(u * a0, (Beta * v) * a1) - a.c1.diff(u * a1, v * a0) - a.c0 *= u2v2inv - a.c1 *= u2v2inv - -# 𝔽p6 -# ---------------------------------------------------------------- - -type - Fp4*[C: static Curve] = object - c0*, c1*: Fp2[C] - - Fp6*[C: static Curve] = object - c0*, c1*, c2*: Fp2[C] - - ξ* = NonResidue - # We call the non-residue ξ on 𝔽p4/𝔽p6 to avoid confusion - # between non-residue - # of different tower level - -func `*`*(_: typedesc[ξ], a: Fp2): Fp2 {.inline, noInit.} = - ## Multiply an element of 𝔽p2 by the quadratic and cubic non-residue - ## chosen to construct 𝔽p4/𝔽p6 - # Yet another const tuple unpacking bug - const u = Fp2.C.get_CNR_Fp2()[0] # Quadratic & Cubic non-residue to construct 𝔽p4/𝔽p6 - const v = Fp2.C.get_CNR_Fp2()[1] - const Beta {.used.} = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 - # ξ = u + v x - # and x² = β - # - # (c0 + c1 x) (u + v x) => u c0 + (u c0 + u c1)x + v c1 x² - # => u c0 + β v c1 + (v c0 + u c1) x - - # TODO: check code generated when ξ = 1 + 𝑖 - # The mul by constant are inline but - # since we don't have __restrict tag - # and we use arrays (which decay into pointer) - # the compiler might not elide the temporary - when a.fromComplexExtension(): - result.c0.diff(u * a.c0, v * a.c1) - else: - result.c0.sum(u * a.c0, (Beta * v) * a.c1) - result.c1.sum(v * a.c0, u * a.c1 ) - -func `*=`*(a: var Fp2, _: typedesc[ξ]) {.inline.} = - ## Multiply an element of 𝔽p2 by the quadratic non-residue - ## chosen to construct 𝔽p6 - # Yet another const tuple unpacking bug - const u = Fp2.C.get_CNR_Fp2()[0] # Cubic non-residue to construct 𝔽p6 - const v = Fp2.C.get_CNR_Fp2()[1] - const Beta {.used.} = Fp2.C.get_QNR_Fp() # Quadratic non-residue to construct 𝔽p2 - # ξ = u + v x - # and x² = β - # - # (c0 + c1 x) (u + v x) => u c0 + (u c0 + u c1)x + v c1 x² - # => u c0 + β v c1 + (v c0 + u c1) x - when a.fromComplexExtension() and u == 1 and v == 1: - let t = a.c0 - a.c0 -= a.c1 - a.c1 += t - else: # TODO: optim for inline - a = ξ * a - -# 𝔽p12 -# ---------------------------------------------------------------- - -type - Fp12*[C: static Curve] = object - c0*, c1*, c2*: Fp4[C] - - γ = NonResidue - # We call the non-residue γ (Gamma) on 𝔽p6 to avoid confusion between non-residue - # of different tower level - -func `*`*(_: typedesc[γ], a: Fp4): Fp4 {.noInit, inline.} = - ## Multiply an element of 𝔽p6 by the cubic non-residue - ## chosen to construct 𝔽p12 - ## For all curves γ = v with v the factor for 𝔽p6 coordinate - ## and v³ = ξ - ## (c0 + c1 v + c2 v²) v => ξ c2 + c0 v + c1 v² - result.c0 = ξ * a.c1 - result.c1 = a.c0 - -func `*=`*(a: var Fp4, _: typedesc[γ]) {.inline.} = - a = γ * a - -# Sparse functions -# ---------------------------------------------------------------- - -func `*=`*(a: var Fp2, b: Fp) = - ## Multiply an element of Fp2 by an element of Fp - a.c0 *= b - a.c1 *= b - -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₁) - r.c0.prod(a.c1, b) - r.c0 *= NonResidue - r.c1.prod(a.c0, b) - -func mul_sparse_by_0y0*[C: static Curve](r: var Fp6[C], a: Fp6[C], b: Fp2[C]) = - ## Sparse multiplication of an Fp6 element - ## with coordinates (a₀, a₁, a₂) by (0, b₁, 0) - # TODO: make generic and move to tower_field_extensions - - # v0 = a0 b0 = 0 - # v1 = a1 b1 - # v2 = a2 b2 = 0 - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b1 + a2 b1 - v1) - # = ξ a2 b1 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = a0 b1 + a1 b1 - v1 - # = a0 b1 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = v1 - # = a1 b1 - - r.c0.prod(a.c2, b) - r.c0 *= ξ - r.c1.prod(a.c0, b) - r.c2.prod(a.c1, b) +export tower_instantiation, square_root_fp2 diff --git a/sage/derive_square_root.sage b/sage/derive_square_root.sage new file mode 100644 index 0000000..0589df0 --- /dev/null +++ b/sage/derive_square_root.sage @@ -0,0 +1,151 @@ +#!/usr/bin/sage +# vim: syntax=python +# vim: set ts=2 sw=2 et: + +# 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. + +# ############################################################ +# +# Frobenius constants +# +# ############################################################ + +# Imports +# --------------------------------------------------------- + +import os +import inspect, textwrap + +# Working directory +# --------------------------------------------------------- + +os.chdir(os.path.dirname(__file__)) + +# Sage imports +# --------------------------------------------------------- +# Accelerate arithmetic by accepting probabilistic proofs +from sage.structure.proof.all import arithmetic +arithmetic(False) + +load('curves.sage') + +# Utilities +# --------------------------------------------------------- + +def fp2_to_hex(a): + v = vector(a) + return '0x' + Integer(v[0]).hex() + ' + β * ' + '0x' + Integer(v[1]).hex() + +def field_to_nim(value, field, curve, prefix = "", comment_above = "", comment_right = ""): + result = '# ' + comment_above + '\n' if comment_above else '' + comment_right = ' # ' + comment_right if comment_right else '' + + if field == 'Fp2': + v = vector(value) + + result += inspect.cleandoc(f""" + {prefix}Fp2[{curve}].fromHex( {comment_right} + "0x{Integer(v[0]).hex()}", + "0x{Integer(v[1]).hex()}" + )""") + elif field == 'Fp': + result += inspect.cleandoc(f""" + {prefix}Fp[{curve}].fromHex( {comment_right} + "0x{Integer(value).hex()}") + """) + else: + raise NotImplementedError() + + return result + +# Code generators +# --------------------------------------------------------- + +def genSqrtFp2Constants(curve_name, curve_config): + embdeg = curve_config[curve_name]['tower']['embedding_degree'] + twdeg = curve_config[curve_name]['tower']['twist_degree'] + g2field = f'Fp{embdeg//twdeg}' if (embdeg//twdeg) > 1 else 'Fp' + + p = curve_config[curve_name]['field']['modulus'] + Fp = GF(p) + K. = PolynomialRing(Fp) + if g2field == 'Fp2': + QNR_Fp = curve_config[curve_name]['tower']['QNR_Fp'] + Fp2. = Fp.extension(u^2 - QNR_Fp) + else: + SNR_Fp = curve_config[curve_name]['tower']['SNR_Fp'] + Fp2. = Fp.extension(u^2 - SNR_Fp) + + sqrt_QNR = Fp2([0, 1]) + sqrt_sqrt_QNR = sqrt_QNR.sqrt() + sqrt_minus_sqrt_QNR = (-sqrt_QNR).sqrt() + + print('\n----> Square root on Fp2 constants <----\n') + buf = inspect.cleandoc(f""" + # Square Root Fp2 constants + # ----------------------------------------------------------------- + """) + buf += '\n' + + buf += f'const {curve_name}_sqrt_QNR* = ' + buf += field_to_nim(sqrt_QNR, 'Fp2', curve_name) + buf += '\n' + + buf += f'const {curve_name}_sqrt_sqrt_QNR* = ' + buf += field_to_nim(sqrt_sqrt_QNR, 'Fp2', curve_name) + buf += '\n' + + buf += f'const {curve_name}_sqrt_minus_sqrt_QNR* = ' + buf += field_to_nim(sqrt_minus_sqrt_QNR, 'Fp2', curve_name) + buf += '\n' + + return buf + +# CLI +# --------------------------------------------------------- + +if __name__ == "__main__": + # Usage + # BLS12-381 + # sage sage/derive_sqrt.sage BLS12_381 + + from argparse import ArgumentParser + + parser = ArgumentParser() + parser.add_argument("curve",nargs="+") + args = parser.parse_args() + + curve = args.curve[0] + + if curve not in Curves: + raise ValueError( + curve + + ' is not one of the available curves: ' + + str(Curves.keys()) + ) + else: + sqrt = genSqrtFp2Constants(curve, Curves) + + with open(f'{curve.lower()}_square_root.nim', 'w') as f: + f.write(copyright()) + f.write('\n\n') + + embdeg = Curves[curve]['tower']['embedding_degree'] + twdeg = Curves[curve]['tower']['twist_degree'] + + f.write(inspect.cleandoc(""" + import + ../config/curves, + ../io/io_towers + """)) + + f.write('\n\n') + f.write(sqrt) + + print(f'Successfully created {curve}_sqrt_fp2.nim')