Sqrt fp2 acceleration (#109)

* Use sqrt.square() == a instead of sqrt * invsqrt = -1 (Euler criterion) for sqrt existence.

* Accelerate sqrt_fp2 by 33%
This commit is contained in:
Mamy Ratsimbazafy 2020-12-13 17:08:32 +01:00 committed by GitHub
parent f0b18ecfe0
commit a022db1c08
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
13 changed files with 723 additions and 336 deletions

View File

@ -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) ## This assumes that the prime field modulus ``p``: p ≡ 3 (mod 4)
sqrt_invsqrt_p3mod4(sqrt, invsqrt, a) sqrt_invsqrt_p3mod4(sqrt, invsqrt, a)
var euler {.noInit.}: Fp[C] var test {.noInit.}: Fp[C]
euler.prod(sqrt, invsqrt) test.square(sqrt)
result = not(euler.isMinusOne()) result = test == a
func sqrt_if_square_p3mod4[C](a: var Fp[C]): SecretBool {.inline.} = func sqrt_if_square_p3mod4[C](a: var Fp[C]): SecretBool {.inline.} =
## If ``a`` is a square, compute the square root of ``a`` ## If ``a`` is a square, compute the square root of ``a``

View File

@ -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"
)

View File

@ -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"
)

View File

@ -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"
)

View File

@ -9,7 +9,7 @@
import import
std/macros, std/macros,
../config/curves, ../config/curves,
./bls12_377_square_root ./bls12_377_sqrt
{.experimental: "dynamicBindSym".} {.experimental: "dynamicBindSym".}

View File

@ -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)

View File

@ -13,7 +13,9 @@ import
# Internal # Internal
./io_bigints, ./io_fields, ./io_bigints, ./io_fields,
../arithmetic/finite_fields, ../arithmetic/finite_fields,
../towers ../tower_field_extensions/tower_instantiation
export tower_instantiation
# No exceptions allowed # No exceptions allowed
{.push raises: [].} {.push raises: [].}

View File

@ -24,11 +24,6 @@ import
# Square root should be implemented in constant-time for hash-to-curve: # 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 # 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 # No exceptions allowed
{.push raises: [].} {.push raises: [].}
@ -194,83 +189,3 @@ func powUnsafeExponent*[F; bits: static int](
var expBE {.noInit.}: array[(bits + 7) div 8, byte] var expBE {.noInit.}: array[(bits + 7) div 8, byte]
expBE.exportRawUint(exponent, bigEndian) expBE.exportRawUint(exponent, bigEndian)
a.powUnsafeExponent(expBE, window) 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)

View File

@ -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()

View File

@ -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)

View File

@ -7,251 +7,9 @@
# at your option. This file may not be copied, modified, or distributed except according to those terms. # at your option. This file may not be copied, modified, or distributed except according to those terms.
import import
./arithmetic, tower_field_extensions/[
./config/curves, tower_instantiation,
./io/io_fields, square_root_fp2
./tower_field_extensions/[
tower_common,
quadratic_extensions,
cubic_extensions,
exponentiations
] ]
export tower_common, quadratic_extensions, cubic_extensions, exponentiations export tower_instantiation, square_root_fp2
# 𝔽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)

View File

@ -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.<u> = PolynomialRing(Fp)
if g2field == 'Fp2':
QNR_Fp = curve_config[curve_name]['tower']['QNR_Fp']
Fp2.<beta> = Fp.extension(u^2 - QNR_Fp)
else:
SNR_Fp = curve_config[curve_name]['tower']['SNR_Fp']
Fp2.<beta> = 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')