diff --git a/constantine.nimble b/constantine.nimble index a129b46..a2d74f6 100644 --- a/constantine.nimble +++ b/constantine.nimble @@ -62,6 +62,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[ ("tests/t_ec_wstrass_prj_g2_mul_distri_bls12_381.nim", false), ("tests/t_ec_wstrass_prj_g2_mul_vs_ref_bls12_381.nim", false), # Elliptic curve arithmetic vs Sagemath + ("tests/t_ec_frobenius.nim", false), ("tests/t_ec_sage_bn254.nim", false), ("tests/t_ec_sage_bls12_381.nim", false), # Edge cases highlighted by past bugs diff --git a/constantine/arithmetic/finite_fields.nim b/constantine/arithmetic/finite_fields.nim index e68a936..e4b1983 100644 --- a/constantine/arithmetic/finite_fields.nim +++ b/constantine/arithmetic/finite_fields.nim @@ -32,6 +32,11 @@ import when UseASM_X86_64: import ./limbs_asm_modular_x86 +when nimvm: + from ../config/precompute import montyResidue_precompute +else: + discard + export Fp # No exceptions allowed @@ -43,13 +48,16 @@ export Fp # # ############################################################ -func fromBig*[C: static Curve](T: type Fp[C], src: BigInt): Fp[C] {.noInit, inline.} = - ## Convert a BigInt to its Montgomery form - result.mres.montyResidue(src, C.Mod, C.getR2modP(), C.getNegInvModWord(), C.canUseNoCarryMontyMul()) - func fromBig*[C: static Curve](dst: var Fp[C], src: BigInt) {.inline.}= ## Convert a BigInt to its Montgomery form - dst.mres.montyResidue(src, C.Mod, C.getR2modP(), C.getNegInvModWord(), C.canUseNoCarryMontyMul()) + when nimvm: + dst.mres.montyResidue_precompute(src, C.Mod, C.getR2modP(), C.getNegInvModWord()) + else: + dst.mres.montyResidue(src, C.Mod, C.getR2modP(), C.getNegInvModWord(), C.canUseNoCarryMontyMul()) + +func fromBig*[C: static Curve](T: type Fp[C], src: BigInt): Fp[C] {.noInit, inline.} = + ## Convert a BigInt to its Montgomery form + result.fromBig(src) func toBig*(src: Fp): auto {.noInit, inline.} = ## Convert a finite-field element to a BigInt in natural representation diff --git a/constantine/config/precompute.nim b/constantine/config/precompute.nim index b6c7f6b..055b2d0 100644 --- a/constantine/config/precompute.nim +++ b/constantine/config/precompute.nim @@ -522,7 +522,7 @@ func montyMul_precompute(r: var BigInt, a, b, M: BigInt, m0ninv: BaseType) = t.limbs[N-1] = SecretWord(tmp) addC(carry, tN, tNp1, 0, carry) - discard t.csub(M, (tN != 0) or not(t < M)) + discard t.csub(M, (tN != 0) or not(precompute.`<`(t, M))) r = t func montyResidue_precompute*(r: var BigInt, a, M, r2modM: BigInt, diff --git a/constantine/elliptic/ec_scalar_mul.nim b/constantine/elliptic/ec_scalar_mul.nim index 9b45747..9ba6529 100644 --- a/constantine/elliptic/ec_scalar_mul.nim +++ b/constantine/elliptic/ec_scalar_mul.nim @@ -226,4 +226,4 @@ func scalarMul*( scratchSpace: array[1 shl 4, ECP_SWei_Proj] scalarCanonicalBE: array[(scalar.bits+7)div 8, byte] # canonical big endian representation scalarCanonicalBE.exportRawUint(scalar, bigEndian) # Export is constant-time - P.scalarMulGeneric(scratchSpace) + P.scalarMulGeneric(scalarCanonicalBE, scratchSpace) diff --git a/constantine/io/io_towers.nim b/constantine/io/io_towers.nim index 4b01504..27e0334 100644 --- a/constantine/io/io_towers.nim +++ b/constantine/io/io_towers.nim @@ -53,3 +53,9 @@ func fromHex*(dst: var Fp2, c0, c1: string) {.raises: [ValueError].}= ## β is the quadratic non-residue chosen to construct 𝔽p2 dst.c0.fromHex(c0) dst.c1.fromHex(c1) + +func fromHex*(T: typedesc[Fp2], c0, c1: string): T {.raises: [ValueError].}= + ## Convert 2 coordinates to an element of 𝔽p2 + ## with dst = c0 + β * c1 + ## β is the quadratic non-residue chosen to construct 𝔽p2 + result.fromHex(c0, c1) diff --git a/constantine/isogeny/frobenius.nim b/constantine/isogeny/frobenius.nim new file mode 100644 index 0000000..6ea3164 --- /dev/null +++ b/constantine/isogeny/frobenius.nim @@ -0,0 +1,185 @@ +# 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/[common, curves], + ../io/io_towers, + ../towers, ../arithmetic + +# Frobenius automorphism +# ------------------------------------------------------------ +# +# https://en.wikipedia.org/wiki/Frobenius_endomorphism +# For p prime +# +# a^p (mod p) ≡ a (mod p) +# +# Also +# +# (a + b)^p (mod p) ≡ a^p + b^p (mod p) +# ≡ a + b (mod p) +# +# Because p is prime and all expanded terms (from binomial expansion) +# besides a^p and b^p are divisible by p + +# For 𝔽p2, with `u` the quadratic non-residue (usually the complex 𝑖) +# (a + u b)^p² (mod p²) ≡ a + u^p² b (mod p²) +# +# For 𝔽p2, frobenius acts like the conjugate +# whether u = √-1 = i +# or √-2 or √-5 + +func frobenius*(r: var Fp2, a: Fp2, k: static int = 1) {.inline.} = + ## Computes a^(p^k) + ## The p-power frobenius automorphism on 𝔽p2 + r.c0 = a.c0 + when (k and 1) == 1: + r.c1.neg(a.c1) + else: + r.c1 = a.c1 + +# Frobenius endomorphism +# ------------------------------------------------------------ +# TODO: generate those constants via Sage in a Json file +# and parse at compile-time + +# Constants: +# Assuming embedding degree of 12 and a sextic twist +# with SNR the sextic non-residue +# +# BN254_Snarks is a D-Twist: SNR^((p-1)/6) +const FrobConst_BN254_Snarks_psi1_coef1 = Fp2[BN254_Snarks].fromHex( + "0x1284b71c2865a7dfe8b99fdd76e68b605c521e08292f2176d60b35dadcc9e470", + "0x246996f3b4fae7e6a6327cfe12150b8e747992778eeec7e5ca5cf05f80f362ac" +) +# SNR^((p-1)/3) +const FrobConst_BN254_Snarks_psi1_coef2 = Fp2[BN254_Snarks].fromHex( + "0x2fb347984f7911f74c0bec3cf559b143b78cc310c2c3330c99e39557176f553d", + "0x16c9e55061ebae204ba4cc8bd75a079432ae2a1d0b7c9dce1665d51c640fcba2" +) +# SNR^((p-1)/2) +const FrobConst_BN254_Snarks_psi1_coef3 = Fp2[BN254_Snarks].fromHex( + "0x63cf305489af5dcdc5ec698b6e2f9b9dbaae0eda9c95998dc54014671a0135a", + "0x7c03cbcac41049a0704b5a7ec796f2b21807dc98fa25bd282d37f632623b0e3" +) +# norm(SNR)^((p-1)/3) +const FrobConst_BN254_Snarks_psi2_coef2 = Fp2[BN254_Snarks].fromHex( + "0x30644e72e131a0295e6dd9e7e0acccb0c28f069fbb966e3de4bd44e5607cfd48", + "0x0" +) + +# BLS12_377 is a D-Twist: SNR^((p-1)/6) +const FrobConst_BLS12_377_psi1_coef1 = Fp2[BLS12_377].fromHex( + "0x9a9975399c019633c1e30682567f915c8a45e0f94ebc8ec681bf34a3aa559db57668e558eb0188e938a9d1104f2031", + "0x0" +) +# SNR^((p-1)/3) +const FrobConst_BLS12_377_psi1_coef2 = Fp2[BLS12_377].fromHex( + "0x9b3af05dd14f6ec619aaf7d34594aabc5ed1347970dec00452217cc900000008508c00000000002", + "0x0" +) +# SNR^((p-1)/2) +const FrobConst_BLS12_377_psi1_coef3 = Fp2[BLS12_377].fromHex( + "0x1680a40796537cac0c534db1a79beb1400398f50ad1dec1bce649cf436b0f6299588459bff27d8e6e76d5ecf1391c63", + "0x0" +) +# norm(SNR)^((p-1)/3) +const FrobConst_BLS12_377_psi2_coef2 = Fp2[BLS12_377].fromHex( + "0x9b3af05dd14f6ec619aaf7d34594aabc5ed1347970dec00452217cc900000008508c00000000001", + "0x0" +) + +# BLS12_381 is a M-twist: (1/SNR)^((p-1)/6) +const FrobConst_BLS12_381_psi1_coef1 = Fp2[BLS12_381].fromHex( + "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116", + "0x5b2cfd9013a5fd8df47fa6b48b1e045f39816240c0b8fee8beadf4d8e9c0566c63a3e6e257f87329b18fae980078116" +) +# (1/SNR)^((p-1)/3) +const FrobConst_BLS12_381_psi1_coef2 = Fp2[BLS12_381].fromHex( + "0x0", + "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad" +) +# (1/SNR)^((p-1)/2) +const FrobConst_BLS12_381_psi1_coef3 = Fp2[BLS12_381].fromHex( + "0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2", + "0x6af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09" +) +# norm(SNR)^((p-1)/3) +const FrobConst_BLS12_381_psi2_coef2 = Fp2[BLS12_381].fromHex( + "0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac", + "0x0" +) + +{.experimental: "dynamicBindSym".} + +macro frobConst(C: static Curve, psipow, coefpow: static int): untyped = + return bindSym("FrobConst_" & $C & "_psi" & $psipow & "_coef" & $coefpow) + +template mulCheckSparse[Fp2](a: var Fp2, b: Fp2) = + when b.c0.isZero().bool: + a.mul_sparse_by_0y(b) + elif b.c1.isZero().bool: + a.mul_sparse_by_x0(b) + else: + a *= b + +func frobenius_psi*[PointG2](r: var PointG2, P: PointG2) = + ## "Untwist-Frobenius-Twist" endomorphism + ## r = ψ(P) + for coordR, coordP in fields(r, P): + coordR.frobenius(coordP, 1) + + # With ξ (xi) the sextic non-residue + # c = ξ^((p-1)/6) for D-Twist + # c = (1/ξ)^((p-1)/6) for M-Twist + # + # c1_2 = c² + # c1_3 = c³ + + r.x.mulCheckSparse frobConst(PointG2.F.C, psipow=1, coefpow=2) + r.y.mulCheckSparse frobConst(PointG2.F.C, psipow=1, coefpow=3) + +func frobenius_psi2*[PointG2](r: var PointG2, P: PointG2) = + ## "Untwist-Frobenius-Twist" endomorphism applied twice + ## r = ψ(ψ(P)) + for coordR, coordP in fields(r, P): + coordR.frobenius(coordP, 2) + + # With ξ (xi) the sextic non-residue + # c = ξ for D-Twist + # c = (1/ξ) for M-Twist + # + # frobenius(a) = conj(a) = a^p + # + # c1_2 = (c^((p-1)/6))² = c^((p-1)/3) + # c1_3 = (c^((p-1)/6))³ = c^((p-1)/2) + # + # c2_2 = c1_2 * frobenius(c1_2) = c^((p-1)/3) * c^((p-1)/3)^p + # = c^((p-1)/3) * conj(c)^((p-1)/3) + # = norm(c)^((p-1)/3) + # + # c2_3 = c1_3 * frobenius(c1_3) = c^((p-1)/2) * c^((p-1)/2)^p + # = c^((p-1)/2) * conj(c)^((p-1)/2) + # = norm(c)^((p-1)/2) + # We prove that c2_3 ≡ -1 (mod p²) with the following: + # + # - Whether c = ξ or c = (1/ξ), c is a quadratic non-residue (QNR) in 𝔽p2 + # because: + # - ξ is quadratic non-residue as it is a sextic non-residue + # by construction of the tower extension + # - if a is QNR then 1/a is also a QNR + # - Then c^((p²-1)/2) ≡ -1 (mod p²) from the Legendre symbol in 𝔽p2 + # + # c2_3 = c^((p-1)/2) * c^((p-1)/2)^p = c^((p+1)*(p-1)/2) + # = c^((p²-1)/2) + # c2_3 ≡ -1 (mod p²) + # QED + + r.x.mulCheckSparse frobConst(PointG2.F.C, psipow=2, coefpow=2) + r.y.neg(r.y) diff --git a/constantine/primitives/constant_time.nim b/constantine/primitives/constant_time.nim index 715f8c6..f01887a 100644 --- a/constantine/primitives/constant_time.nim +++ b/constantine/primitives/constant_time.nim @@ -117,7 +117,9 @@ template `-`*[T: Ct](x: T): T = ## Unary minus returns the two-complement representation ## of an unsigned integer # We could use "not(x) + 1" but the codegen is not optimal - block: + when nimvm: + not(x) + T(1) + else: # Use C so that compiler uses the "neg" instructions var neg: T {.emit:[neg, " = -", x, ";"].} neg diff --git a/constantine/tower_field_extensions/quadratic_extensions.nim b/constantine/tower_field_extensions/quadratic_extensions.nim index 48767c7..6a5e65f 100644 --- a/constantine/tower_field_extensions/quadratic_extensions.nim +++ b/constantine/tower_field_extensions/quadratic_extensions.nim @@ -138,6 +138,26 @@ func prod_complex(r: var QuadraticExt, a, b: QuadraticExt) = # - Function calls? # - push/pop stack? +func mul_sparse_complex_by_0y(r: var QuadraticExt, a, sparseB: QuadraticExt) = + ## Multiply `a` by `b` with sparse coordinates (0, y) + ## On a complex quadratic extension field 𝔽p2 = 𝔽p[𝑖] + # + # r0 = a0 b0 - a1 b1 + # r1 = (a0 + a1) (b0 + b1) - a0 b0 - a1 b1 (Karatsuba) + # + # with b0 = 0, hence + # + # r0 = - a1 b1 + # r1 = (a0 + a1) b1 - a1 b1 = a0 b1 + mixin fromComplexExtension + static: doAssert r.fromComplexExtension() + + template b(): untyped = sparseB + + r.c0.prod(a.c1, b.c1) + r.c0.neg(r.c0) + r.c1.prod(a.c0, b.c1) + # Commutative ring implementation for generic quadratic extension fields # ------------------------------------------------------------------- @@ -205,6 +225,24 @@ func prod_generic(r: var QuadraticExt, a, b: QuadraticExt) = # r0 <- a0 b0 + β a1 b1 r.c0 += NonResidue * t +func mul_sparse_generic_by_x0(r: var QuadraticExt, a, sparseB: QuadraticExt) = + ## Multiply `a` by `b` with sparse coordinates (x, 0) + ## On a generic quadratic extension field + # Algorithm (with β the non-residue in the base field) + # + # r0 = a0 b0 + β a1 b1 + # r1 = (a0 + a1) (b0 + b1) - a0 b0 - a1 b1 (Karatsuba) + # + # with b1 = 0, hence + # + # r0 = a0 b0 + # r1 = (a0 + a1) b0 - a0 b0 = a1 b0 + mixin prod + template b(): untyped = sparseB + + r.c0.prod(a.c0, b.c0) + r.c1.prod(a.c1, b.c0) + # Exported symbols # ------------------------------------------------------------------- @@ -263,6 +301,20 @@ func `*=`*(a: var QuadraticExt, b: QuadraticExt) {.inline.} = let t = a a.prod(t, b) -func square*(a: var QuadraticExt){.inline.} = +func square*(a: var QuadraticExt) {.inline.} = let t = a a.square(t) + +func mul_sparse_by_0y*(a: var QuadraticExt, sparseB: QuadraticExt) {.inline.} = + ## Sparse in-place multiplication + mixin fromComplexExtension + when a.fromComplexExtension(): + let t = a + a.mul_sparse_complex_by_0y(t, sparseB) + else: + {.error: "Not implemented".} + +func mul_sparse_by_x0*(a: var QuadraticExt, sparseB: QuadraticExt) {.inline.} = + ## Sparse in-place multiplication + let t = a + a.mul_sparse_generic_by_x0(t, sparseB) diff --git a/sage/curve_family_bls12.sage b/sage/curve_family_bls12.sage index 6f34e46..8a8f61c 100644 --- a/sage/curve_family_bls12.sage +++ b/sage/curve_family_bls12.sage @@ -14,21 +14,23 @@ # ############################################################ # # This module derives a BLS12 curve parameters from -# its base parameter u +# its base parameter x -def compute_curve_characteristic(u_str): - u = sage_eval(u_str) - p = (u - 1)^2 * (u^4 - u^2 + 1)//3 + u - r = u^4 - u^2 + 1 +def compute_curve_characteristic(x_str): + x = sage_eval(x_str) + p = (x - 1)^2 * (x^4 - x^2 + 1)//3 + x + r = x^4 - x^2 + 1 + t = x + 1 print(f'BLS12 family - {p.nbits()} bits') - print(' Prime modulus: 0x' + p.hex()) - print(' Curve order: 0x' + r.hex()) - print(' Parameter u: ' + u_str) - if u < 0: - print(' Parameter u (hex): -0x' + (-u).hex()) + print(' Prime modulus p: 0x' + p.hex()) + print(' Curve order r: 0x' + r.hex()) + print(' trace t: 0x' + t.hex()) + print(' Parameter x: ' + x_str) + if x < 0: + print(' Parameter x (hex): -0x' + (-x).hex()) else: - print(' Parameter u (hex): 0x' + u.hex()) + print(' Parameter x (hex): 0x' + x.hex()) print() @@ -61,11 +63,11 @@ def compute_curve_characteristic(u_str): print(f' GLV-2 decomposition of s into (k1, k2) on G1') print(f' (k1, k2) = (s, 0) - 𝛼1 b1 - 𝛼2 b2') print(f' 𝛼i = 𝛼\u0302i * s / r') - print(f' Lattice b1: ' + str(['0x' + b.hex() for b in [u^2-1, -1]])) - print(f' Lattice b2: ' + str(['0x' + b.hex() for b in [1, u^2]])) + print(f' Lattice b1: ' + str(['0x' + b.hex() for b in [x^2-1, -1]])) + print(f' Lattice b2: ' + str(['0x' + b.hex() for b in [1, x^2]])) # Babai rounding - ahat1 = u^2 + ahat1 = x^2 ahat2 = 1 # We want a1 = ahat1 * s/r with m = 2 (for a 2-dim decomposition) and r the curve order # To handle rounding errors we instead multiply by diff --git a/sage/curve_family_bn.sage b/sage/curve_family_bn.sage index 1b3a006..b353a8d 100644 --- a/sage/curve_family_bn.sage +++ b/sage/curve_family_bn.sage @@ -14,21 +14,23 @@ # ############################################################ # # This module derives a BN curve parameters from -# its base parameter u +# its base parameterx -def compute_curve_characteristic(u_str): - u = sage_eval(u_str) - p = 36*u^4 + 36*u^3 + 24*u^2 + 6*u + 1 - r = 36*u^4 + 36*u^3 + 18*u^2 + 6*u + 1 +def compute_curve_characteristic(x_str): + x = sage_eval(x_str) + p = 36*x^4 + 36*x^3 + 24*x^2 + 6*x + 1 + r = 36*x^4 + 36*x^3 + 18*x^2 + 6*x + 1 + t = 6*x^2 + 1 print(f'BN family - {p.nbits()} bits') print(' Prime modulus p: 0x' + p.hex()) print(' Curve order r: 0x' + r.hex()) - print(' Parameter u: ' + u_str) - if u < 0: - print(' Parameter u (hex): -0x' + (-u).hex()) + print(' trace t: 0x' + t.hex()) + print(' Parameterx: ' + x_str) + if x < 0: + print(' Parameterx (hex): -0x' + (-x).hex()) else: - print(' Parameter u (hex): 0x' + u.hex()) + print(' Parameterx (hex): 0x' +x.hex()) print(f' p mod 3: ' + str(p % 3)) print(f' p mod 4: ' + str(p % 4)) @@ -61,12 +63,12 @@ def compute_curve_characteristic(u_str): print(f' GLV-2 decomposition of s into (k1, k2) on G1') print(f' (k1, k2) = (s, 0) - 𝛼1 b1 - 𝛼2 b2') print(f' 𝛼i = 𝛼\u0302i * s / r') - print(f' Lattice b1: ' + str(['0x' + b.hex() for b in [2*u+1, 6*u^2+4*u+1]])) - print(f' Lattice b2: ' + str(['0x' + b.hex() for b in [6*u^2+2*u, -2*u-1]])) + print(f' Lattice b1: ' + str(['0x' + b.hex() for b in [2*x+1, 6*x^2+4*x+1]])) + print(f' Lattice b2: ' + str(['0x' + b.hex() for b in [6*x^2+2*x, -2*x-1]])) # Babai rounding - ahat1 = 2*u+1 - ahat2 = 6*u^2+4*u+1 + ahat1 = 2*x+1 + ahat2 = 6*x^2+4*x+1 # We want a1 = ahat1 * s/r with m = 2 (for a 2-dim decomposition) and r the curve order # To handle rounding errors we instead multiply by # 𝜈 = (2^WordBitWidth)^w (i.e. the same as the R magic constant for Montgomery arithmetic) diff --git a/sage/frobenius_bls12_377.sage b/sage/frobenius_bls12_377.sage new file mode 100644 index 0000000..473f17a --- /dev/null +++ b/sage/frobenius_bls12_377.sage @@ -0,0 +1,147 @@ +# 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. + +# ############################################################ +# +# BLS12-377 +# Frobenius Endomorphism +# Untwist-Frobenius-Twist isogeny +# +# ############################################################ + +# Parameters +x = 3 * 2^46 * (7 * 13 * 499) + 1 +p = (x - 1)^2 * (x^4 - x^2 + 1)//3 + x +r = x^4 - x^2 + 1 +t = x + 1 +print('p : ' + p.hex()) +print('r : ' + r.hex()) +print('t : ' + t.hex()) + +# Finite fields +Fp = GF(p) +K2. = PolynomialRing(Fp) +Fp2. = Fp.extension(u^2+5) # √-5 quadratic non-residue +# K6. = PolynomialRing(F2) +# Fp6. = Fp2.extension(v^3-Fp2([0, 1]) +# K12. = PolynomialRing(Fp6) +# Fp12. = Fp6.extension(w^2-eta) + +# Curves +b = 1 +SNR = Fp2([0, 1]) # √-5 sextic non-residue +G1 = EllipticCurve(Fp, [0, b]) +G2 = EllipticCurve(Fp2, [0, b/SNR]) + +# Utilities +def fp2_to_hex(a): + v = vector(a) + return Integer(v[0]).hex() + ' + β * ' + Integer(v[1]).hex() + +# Frobenius constants (D type: use SNR, M type use 1/SNR) +FrobConst_psi = SNR^((p-1)/6) +FrobConst_psi_2 = FrobConst_psi * FrobConst_psi +FrobConst_psi_3 = FrobConst_psi_2 * FrobConst_psi +print('FrobConst_psi : ' + fp2_to_hex(FrobConst_psi)) +print('FrobConst_psi_2 : ' + fp2_to_hex(FrobConst_psi_2)) +print('FrobConst_psi_3 : ' + fp2_to_hex(FrobConst_psi_3)) + +print('') +FrobConst_psi2_2 = FrobConst_psi_2 * FrobConst_psi_2**p +FrobConst_psi2_3 = FrobConst_psi_3 * FrobConst_psi_3**p +print('FrobConst_psi2_2 : ' + fp2_to_hex(FrobConst_psi2_2)) +print('FrobConst_psi2_3 : ' + fp2_to_hex(FrobConst_psi2_3)) + +print('') +FrobConst_psi3_2 = FrobConst_psi_2 * FrobConst_psi2_2**p +FrobConst_psi3_3 = FrobConst_psi_3 * FrobConst_psi2_3**p +print('FrobConst_psi3_2 : ' + fp2_to_hex(FrobConst_psi3_2)) +print('FrobConst_psi3_3 : ' + fp2_to_hex(FrobConst_psi3_3)) + +# Recap, with ξ (xi) the sextic non-residue +# psi_2 = (ξ^((p-1)/6))^2 = ξ^((p-1)/3) +# psi_3 = psi_2 * ξ^((p-1)/6) = ξ^((p-1)/3) * ξ^((p-1)/6) = ξ^((p-1)/2) +# +# Reminder, in 𝔽p2, frobenius(a) = a^p = conj(a) +# psi2_2 = psi_2 * psi_2^p = ξ^((p-1)/3) * ξ^((p-1)/3)^p = ξ^((p-1)/3) * frobenius(ξ)^((p-1)/3) +# = norm(ξ)^((p-1)/3) +# psi2_3 = psi_3 * psi_3^p = ξ^((p-1)/2) * ξ^((p-1)/2)^p = ξ^((p-1)/2) * frobenius(ξ)^((p-1)/2) +# = norm(ξ)^((p-1)/2) +# +# In Fp²: +# - quadratic non-residues respect the equation a^((p²-1)/2) ≡ -1 (mod p²) by the Legendre symbol +# - sextic non-residues are also quadratic non-residues so ξ^((p²-1)/2) ≡ -1 (mod p²) +# +# We have norm(ξ)^((p-1)/2) = (ξ*frobenius(ξ))^((p-1)/2) = (ξ*(ξ^p))^((p-1)/2) = ξ^(p+1)^(p-1)/2 +# = ξ^((p²-1)/2) +# And ξ^((p²-1)/2) ≡ -1 (mod p²) +# So psi2_3 ≡ -1 (mod p²) +# +# TODO: explain why psi3_2 = [0, -1] + +# Frobenius Fp2 +A = Fp2([5, 7]) +Aconj = Fp2([5, -7]) +AF = A.frobenius(1) # or pth_power(1) +AF2 = A.frobenius(2) +AF3 = A.frobenius(3) +print('') +print('A : ' + fp2_to_hex(A)) +print('A conjugate: ' + fp2_to_hex(Aconj)) +print('') +print('AF1 : ' + fp2_to_hex(AF)) +print('AF2 : ' + fp2_to_hex(AF2)) +print('AF3 : ' + fp2_to_hex(AF3)) + +def psi(P): + (Px, Py, Pz) = P + return G2([ + FrobConst_psi_2 * Px.frobenius(), + FrobConst_psi_3 * Py.frobenius() + # Pz.frobenius() - Always 1 after extract + ]) + +def psi2(P): + (Px, Py, Pz) = P + return G2([ + FrobConst_psi2_2 * Px.frobenius(2), + FrobConst_psi2_3 * Py.frobenius(2) + # Pz - Always 1 after extract + ]) + +# Test generator +set_random_seed(1337) + +# Vectors +print('\nTest vectors:') +for i in range(4): + P = G2.random_point() + + (Px, Py, Pz) = P + vPx = vector(Px) + vPy = vector(Py) + # Pz = vector(Pz) + print(f'\nTest {i}') + print(' Px: ' + Integer(vPx[0]).hex() + ' + β * ' + Integer(vPx[1]).hex()) + print(' Py: ' + Integer(vPy[0]).hex() + ' + β * ' + Integer(vPy[1]).hex()) + + # Galbraith-Lin-Scott, 2008, Theorem 1 + # Fuentes-Castaneda et al, 2011, Equation (2) + assert psi(psi(P)) - t*psi(P) + p*P == G2([0, 1, 0]) + + # Galbraith-Scott, 2008, Lemma 1 + # k-th cyclotomic polynomial with k = 12 + assert psi2(psi2(P)) - psi2(P) + P == G2([0, 1, 0]) + + assert psi(psi(P)) == psi2(P) + + (Qx, Qy, Qz) = psi(P) + vQx = vector(Qx) + vQy = vector(Qy) + print(' Qx: ' + Integer(vQx[0]).hex() + ' + β * ' + Integer(vQx[1]).hex()) + print(' Qy: ' + Integer(vQy[0]).hex() + ' + β * ' + Integer(vQy[1]).hex()) diff --git a/sage/frobenius_bls12_381.sage b/sage/frobenius_bls12_381.sage new file mode 100644 index 0000000..e3ad1a3 --- /dev/null +++ b/sage/frobenius_bls12_381.sage @@ -0,0 +1,150 @@ +# 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. + +# ############################################################ +# +# BLS12-381 +# Frobenius Endomorphism +# Untwist-Frobenius-Twist isogeny +# +# ############################################################ + +# Parameters +x = -(2^63 + 2^62 + 2^60 + 2^57 + 2^48 + 2^16) +p = (x - 1)^2 * (x^4 - x^2 + 1)//3 + x +r = x^4 - x^2 + 1 +t = x + 1 +print('p : ' + p.hex()) +print('r : ' + r.hex()) +print('t : ' + t.hex()) + +# Finite fields +Fp = GF(p) +K2. = PolynomialRing(Fp) +Fp2. = Fp.extension(u^2+1) +# K6. = PolynomialRing(F2) +# Fp6. = Fp2.extension(v^3-Fp2([1, 1]) +# K12. = PolynomialRing(Fp6) +# Fp12. = Fp6.extension(w^2-eta) + +# Curves +b = 4 +SNR = Fp2([1, 1]) +G1 = EllipticCurve(Fp, [0, b]) +G2 = EllipticCurve(Fp2, [0, b*SNR]) + +# Utilities +def fp2_to_hex(a): + v = vector(a) + return Integer(v[0]).hex() + ' + β * ' + Integer(v[1]).hex() + +# Frobenius constants (D type: use SNR, M type use 1/SNR) +print('1/sextic_non_residue: ' + fp2_to_hex(1/SNR)) +FrobConst_psi = (1/SNR)^((p-1)/6) +FrobConst_psi_2 = FrobConst_psi * FrobConst_psi +FrobConst_psi_3 = FrobConst_psi_2 * FrobConst_psi +print('FrobConst_psi : ' + fp2_to_hex(FrobConst_psi)) +print('FrobConst_psi_2 : ' + fp2_to_hex(FrobConst_psi_2)) +print('FrobConst_psi_3 : ' + fp2_to_hex(FrobConst_psi_3)) + +print('') +FrobConst_psi2_2 = FrobConst_psi_2 * FrobConst_psi_2**p +FrobConst_psi2_3 = FrobConst_psi_3 * FrobConst_psi_3**p +print('FrobConst_psi2_2 : ' + fp2_to_hex(FrobConst_psi2_2)) +print('FrobConst_psi2_3 : ' + fp2_to_hex(FrobConst_psi2_3)) + +print('') +FrobConst_psi3_2 = FrobConst_psi_2 * FrobConst_psi2_2**p +FrobConst_psi3_3 = FrobConst_psi_3 * FrobConst_psi2_3**p +print('FrobConst_psi3_2 : ' + fp2_to_hex(FrobConst_psi3_2)) +print('FrobConst_psi3_3 : ' + fp2_to_hex(FrobConst_psi3_3)) + +# Recap, with ξ (xi) the sextic non-residue +# psi_2 = ((1/ξ)^((p-1)/6))^2 = (1/ξ)^((p-1)/3) +# psi_3 = psi_2 * (1/ξ)^((p-1)/6) = (1/ξ)^((p-1)/3) * (1/ξ)^((p-1)/6) = (1/ξ)^((p-1)/2) +# +# Reminder, in 𝔽p2, frobenius(a) = a^p = conj(a) +# psi2_2 = psi_2 * psi_2^p = (1/ξ)^((p-1)/3) * (1/ξ)^((p-1)/3)^p = (1/ξ)^((p-1)/3) * frobenius((1/ξ))^((p-1)/3) +# = norm(1/ξ)^((p-1)/3) +# psi2_3 = psi_3 * psi_3^p = (1/ξ)^((p-1)/2) * (1/ξ)^((p-1)/2)^p = (1/ξ)^((p-1)/2) * frobenius((1/ξ))^((p-1)/2) +# = norm(1/ξ)^((p-1)/2) +# +# In Fp²: +# - quadratic non-residues respect the equation a^((p²-1)/2) ≡ -1 (mod p²) by the Legendre symbol +# - sextic non-residues are also quadratic non-residues so ξ^((p²-1)/2) ≡ -1 (mod p²) +# - QRT(1/a) = QRT(a) with QRT the quadratic residuosity test +# +# We have norm(ξ)^((p-1)/2) = (ξ*frobenius(ξ))^((p-1)/2) = (ξ*(ξ^p))^((p-1)/2) = ξ^(p+1)^(p-1)/2 +# = ξ^((p²-1)/2) +# And ξ^((p²-1)/2) ≡ -1 (mod p²) +# So psi2_3 ≡ -1 (mod p²) +# +# TODO: explain why psi3_2 = [0, -1] + +# Frobenius Fp2 +A = Fp2([5, 7]) +Aconj = Fp2([5, -7]) +AF = A.frobenius(1) # or pth_power(1) +AF2 = A.frobenius(2) +AF3 = A.frobenius(3) +print('') +print('A : ' + fp2_to_hex(A)) +print('A conjugate: ' + fp2_to_hex(Aconj)) +print('') +print('AF1 : ' + fp2_to_hex(AF)) +print('AF2 : ' + fp2_to_hex(AF2)) +print('AF3 : ' + fp2_to_hex(AF3)) + +def psi(P): + (Px, Py, Pz) = P + return G2([ + FrobConst_psi_2 * Px.frobenius(), + FrobConst_psi_3 * Py.frobenius() + # Pz.frobenius() - Always 1 after extract + ]) + +def psi2(P): + (Px, Py, Pz) = P + return G2([ + FrobConst_psi2_2 * Px.frobenius(2), + FrobConst_psi2_3 * Py.frobenius(2) + # Pz - Always 1 after extract + ]) + +# Test generator +set_random_seed(1337) + +# Vectors +print('\nTest vectors:') +for i in range(4): + P = G2.random_point() + + (Px, Py, Pz) = P + vPx = vector(Px) + vPy = vector(Py) + # Pz = vector(Pz) + print(f'\nTest {i}') + print(' Px: ' + Integer(vPx[0]).hex() + ' + β * ' + Integer(vPx[1]).hex()) + print(' Py: ' + Integer(vPy[0]).hex() + ' + β * ' + Integer(vPy[1]).hex()) + + # Galbraith-Lin-Scott, 2008, Theorem 1 + # Fuentes-Castaneda et al, 2011, Equation (2) + assert psi(psi(P)) - t*psi(P) + p*P == G2([0, 1, 0]) + + # Galbraith-Scott, 2008, Lemma 1 + # k-th cyclotomic polynomial with k = 12 + assert psi2(psi2(P)) - psi2(P) + P == G2([0, 1, 0]) + + assert psi(psi(P)) == psi2(P) + + + (Qx, Qy, Qz) = psi(P) + vQx = vector(Qx) + vQy = vector(Qy) + print(' Qx: ' + Integer(vQx[0]).hex() + ' + β * ' + Integer(vQx[1]).hex()) + print(' Qy: ' + Integer(vQy[0]).hex() + ' + β * ' + Integer(vQy[1]).hex()) diff --git a/sage/frobenius_bn254_snarks.sage b/sage/frobenius_bn254_snarks.sage new file mode 100644 index 0000000..8f4c384 --- /dev/null +++ b/sage/frobenius_bn254_snarks.sage @@ -0,0 +1,147 @@ +# 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. + +# ############################################################ +# +# BLS12-381 +# Frobenius Endomorphism +# Untwist-Frobenius-Twist isogeny +# +# ############################################################ + +# Parameters +x = Integer('0x44E992B44A6909F1') +p = 36*x^4 + 36*x^3 + 24*x^2 + 6*x + 1 +r = 36*x^4 + 36*x^3 + 18*x^2 + 6*x + 1 +t = 6*x^2 + 1 +cofactor = 1 +print('p : ' + p.hex()) +print('r : ' + r.hex()) +print('t : ' + t.hex()) + +# Finite fields +Fp = GF(p) +K2. = PolynomialRing(Fp) +Fp2. = Fp.extension(u^2+1) +# K6. = PolynomialRing(F2) +# Fp6. = Fp2.extension(v^3-Fp2([9, 1])) +# K12. = PolynomialRing(Fp6) +# K12. = F6.extension(w^2-eta) + +# Curves +b = 3 +SNR = Fp2([9, 1]) +G1 = EllipticCurve(Fp, [0, b]) +G2 = EllipticCurve(Fp2, [0, b/SNR]) + +# Utilities +def fp2_to_hex(a): + v = vector(a) + return Integer(v[0]).hex() + ' + β * ' + Integer(v[1]).hex() + +# Frobenius constants (D type: use SNR, M type use 1/SNR) +FrobConst_psi = SNR^((p-1)/6) +FrobConst_psi_2 = FrobConst_psi * FrobConst_psi +FrobConst_psi_3 = FrobConst_psi_2 * FrobConst_psi +print('FrobConst_psi : ' + fp2_to_hex(FrobConst_psi)) +print('FrobConst_psi_2 : ' + fp2_to_hex(FrobConst_psi_2)) +print('FrobConst_psi_3 : ' + fp2_to_hex(FrobConst_psi_3)) + +print('') +FrobConst_psi2_2 = FrobConst_psi_2 * FrobConst_psi_2^p +FrobConst_psi2_3 = FrobConst_psi_3 * FrobConst_psi_3^p +print('FrobConst_psi2_2 : ' + fp2_to_hex(FrobConst_psi2_2)) +print('FrobConst_psi2_3 : ' + fp2_to_hex(FrobConst_psi2_3)) + +print('') +FrobConst_psi3_2 = FrobConst_psi_2 * FrobConst_psi2_2^p +FrobConst_psi3_3 = FrobConst_psi_3 * FrobConst_psi2_3^p +print('FrobConst_psi3_2 : ' + fp2_to_hex(FrobConst_psi3_2)) +print('FrobConst_psi3_3 : ' + fp2_to_hex(FrobConst_psi3_3)) + +# Recap, with ξ (xi) the sextic non-residue +# psi_2 = (ξ^((p-1)/6))^2 = ξ^((p-1)/3) +# psi_3 = psi_2 * ξ^((p-1)/6) = ξ^((p-1)/3) * ξ^((p-1)/6) = ξ^((p-1)/2) +# +# Reminder, in 𝔽p2, frobenius(a) = a^p = conj(a) +# psi2_2 = psi_2 * psi_2^p = ξ^((p-1)/3) * ξ^((p-1)/3)^p = ξ^((p-1)/3) * frobenius(ξ)^((p-1)/3) +# = norm(ξ)^((p-1)/3) +# psi2_3 = psi_3 * psi_3^p = ξ^((p-1)/2) * ξ^((p-1)/2)^p = ξ^((p-1)/2) * frobenius(ξ)^((p-1)/2) +# = norm(ξ)^((p-1)/2) +# +# In Fp²: +# - quadratic non-residues respect the equation a^((p²-1)/2) ≡ -1 (mod p²) by the Legendre symbol +# - sextic non-residues are also quadratic non-residues so ξ^((p²-1)/2) ≡ -1 (mod p²) +# +# We have norm(ξ)^((p-1)/2) = (ξ*frobenius(ξ))^((p-1)/2) = (ξ*(ξ^p))^((p-1)/2) = ξ^(p+1)^(p-1)/2 +# = ξ^((p²-1)/2) +# And ξ^((p²-1)/2) ≡ -1 (mod p²) +# So psi2_3 ≡ -1 (mod p²) + +# Frobenius Fp2 +A = Fp2([5, 7]) +Aconj = Fp2([5, -7]) +AF = A.frobenius(1) # or pth_power(1) +AF2 = A.frobenius(2) +AF3 = A.frobenius(3) +print('') +print('A : ' + fp2_to_hex(A)) +print('A conjugate: ' + fp2_to_hex(Aconj)) +print('') +print('AF1 : ' + fp2_to_hex(AF)) +print('AF2 : ' + fp2_to_hex(AF2)) +print('AF3 : ' + fp2_to_hex(AF3)) + +def psi(P): + (Px, Py, Pz) = P + return G2([ + FrobConst_psi_2 * Px.frobenius(), + FrobConst_psi_3 * Py.frobenius() + # Pz.frobenius() - Always 1 after extract + ]) + +def psi2(P): + (Px, Py, Pz) = P + return G2([ + FrobConst_psi2_2 * Px.frobenius(2), + FrobConst_psi2_3 * Py.frobenius(2) + # Pz - Always 1 after extract + ]) + +# Test generator +set_random_seed(1337) + +# Vectors +print('\nTest vectors:') +for i in range(4): + P = G2.random_point() + + (Px, Py, Pz) = P + vPx = vector(Px) + vPy = vector(Py) + # vPz = vector(Pz) + print(f'\nTest {i}') + print(' Px: ' + Integer(vPx[0]).hex() + ' + β * ' + Integer(vPx[1]).hex()) + print(' Py: ' + Integer(vPy[0]).hex() + ' + β * ' + Integer(vPy[1]).hex()) + # print(' Pz: ' + Integer(vPz[0]).hex() + ' + β * ' + Integer(vPz[1]).hex()) + + # Galbraith-Lin-Scott, 2008, Theorem 1 + # Fuentes-Castaneda et al, 2011, Equation (2) + assert psi(psi(P)) - t*psi(P) + p*P == G2([0, 1, 0]) + + # Galbraith-Scott, 2008, Lemma 1 + # k-th cyclotomic polynomial with k = 12 + assert psi2(psi2(P)) - psi2(P) + P == G2([0, 1, 0]) + + assert psi(psi(P)) == psi2(P) + + (Qx, Qy, Qz) = psi(P) + vQx = vector(Qx) + vQy = vector(Qy) + print(' Qx: ' + Integer(vQx[0]).hex() + ' + β * ' + Integer(vQx[1]).hex()) + print(' Qy: ' + Integer(vQy[0]).hex() + ' + β * ' + Integer(vQy[1]).hex()) diff --git a/tests/t_bigints.nim b/tests/t_bigints.nim index be43128..b338321 100644 --- a/tests/t_bigints.nim +++ b/tests/t_bigints.nim @@ -26,15 +26,30 @@ proc mainArith() = check: x.isZero().bool test "isZero for non-zero": block: - var x = fromHex(BigInt[128], "0x00000000_00000000_00000000_00000001") + let x = fromHex(BigInt[128], "0x00000000_00000000_00000000_00000001") check: not x.isZero().bool block: - var x = fromHex(BigInt[128], "0x00000000_00000001_00000000_00000000") + let x = fromHex(BigInt[128], "0x00000000_00000001_00000000_00000000") check: not x.isZero().bool block: - var x = fromHex(BigInt[128], "0xFFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFF") + let x = fromHex(BigInt[128], "0xFFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFF") check: not x.isZero().bool + test "isZero for zero (compile-time)": + const x = BigInt[128]() + check: static(x.isZero().bool) + test "isZero for non-zero (compile-time)": + block: + const x = fromHex(BigInt[128], "0x00000000_00000000_00000000_00000001") + check: static(not x.isZero().bool) + block: + const x = fromHex(BigInt[128], "0x00000000_00000001_00000000_00000000") + check: static(not x.isZero().bool) + block: + const x = fromHex(BigInt[128], "0xFFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFF") + check: static(not x.isZero().bool) + + suite "Arithmetic operations - Addition" & " [" & $WordBitwidth & "-bit mode]": test "Adding 2 zeros": var a = fromHex(BigInt[128], "0x00000000_00000000_00000000_00000000") @@ -622,9 +637,9 @@ proc mainModularInverse() = check: bool(r == expected) mainArith() -mainMul() -mainMulHigh() -mainModular() -mainNeg() -mainCopySwap() -mainModularInverse() +# mainMul() +# mainMulHigh() +# mainModular() +# mainNeg() +# mainCopySwap() +# mainModularInverse() diff --git a/tests/t_ec_frobenius.nim b/tests/t_ec_frobenius.nim new file mode 100644 index 0000000..a8c5f26 --- /dev/null +++ b/tests/t_ec_frobenius.nim @@ -0,0 +1,324 @@ +# 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 + # Standard library + std/[times, unittest], + # Internals + ../constantine/config/[common, curves], + ../constantine/[arithmetic, towers], + ../constantine/io/[io_bigints, io_ec], + ../constantine/elliptic/[ec_weierstrass_affine, ec_weierstrass_projective, ec_scalar_mul], + ../constantine/isogeny/frobenius, + # Tests + ../helpers/prng_unsafe, + ./t_ec_template + +echo "\n------------------------------------------------------\n" + +# Random seed for reproducibility +var rng: RngState +let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32 +rng.seed(seed) +echo "frobenius xoshiro512** seed: ", seed + +proc test( + id: int, + EC: typedesc[ECP_SWei_Proj], + Px0, Px1, Py0, Py1: string, + Qx0, Qx1, Qy0, Qy1: string + ) = + + test "test " & $EC & " - " & $id: + var P: EC + let pOK = P.fromHex(Px0, Px1, Py0, Py1) + doAssert pOK + + var Q: EC + let qOK = Q.fromHex(Qx0, Qx1, Qy0, Qy1) + + var R{.noInit.}: EC + + R.frobenius_psi(P) + doAssert: bool(R == Q) + +suite "ψ (Psi) - Untwist-Frobenius-Twist Endomorphism on G2 vs SageMath" & " [" & $WordBitwidth & "-bit mode]": + # Generated via + # - sage sage/frobenius_bn254_snarks.sage + # - sage sage/frobenius_bls12_377.sage + # - sage sage/frobenius_bls12_381.sage + test( + id = 0, + EC = ECP_SWei_Proj[Fp2[BN254_Snarks]], + Px0 = "598e4c8c14c24c90834f2debedee4db3d31fed98a5134177704bfec14f46cb5", + Px1 = "c6fffa61daeb7caaf96983e70f164931d958c6820b205cdde19f2fa1eaaa7b1", + Py0 = "2f5fa252a27df56f5ca2e9c3382c17e531d317d50396f3fe952704304946a5a", + Py1 = "25f58d3c1a91b83ad645c6cf43f63b5cc36c08d920eaa9c2d7e5a8f8d9f01c6a", + Qx0 = "14ce20d36d5ad6f84f3f13a4db7dc2e4cc1b568af633f0c1164fbb7804738071", + Qx1 = "1689941c634aca12aecb1e20e04e9d6b22fc012b534e30436bf0cec03ff0ddde", + Qy0 = "2eb451e1a4860bcc16a0f898c4fa5b9638975ae20fa3096ed9f3127ca6f0b6de", + Qy1 = "317e8196e9e77ae0139c6fa6b56f7df1dd88e1f472a6c60e1c1c63065ebc71f" + ) + + test( + id = 1, + EC = ECP_SWei_Proj[Fp2[BN254_Snarks]], + Px0 = "21014830dd88a0e7961e704cea531200866c5df46cb25aa3e2aac8d4fec64c6e", + Px1 = "1db17d8364def10443beab6e4a055c210d3e49c7c3af31e9cfb66d829938dca7", + Py0 = "1394ab8c346ad3eba14fa14789d3bbfc2deed5a7a510da8e9418580515d27bda", + Py1 = "1157a26b2639a5258d62577efad8eadffd5814ada3c9b0ce6cae8d29ccafbec9", + Qx0 = "2099f3e18a043417f28d7eac1ba96856df748201ed667467bb40b14e389e4eb1", + Qx1 = "2870f4936ec8a10b45b407e493a883da8fb2ebd68dbb35f998bfe7a63d512498", + Qy0 = "1dfe671f9556e11d771b0a3f0d3884883416866c889ef59cb26a1b9d74c24bb", + Qy1 = "a348314ec66afb05f015b6f6868e3f483bb03b591d4e53d0c9d9d165b4a9d8" + ) + + test( + id = 2, + EC = ECP_SWei_Proj[Fp2[BN254_Snarks]], + Px0 = "46f2a2be9a3e19c1bb484fc37703ff64c3d7379de22249ccf0881037948beec", + Px1 = "10a5aaae14cb028f4ff4b81d41b712038b9f620a99e208c23504887e56831806", + Py0 = "2e6c3ebe0f3dada0063dc59f85fe2264dc3502bf65206336106a8d39d838a7b2", + Py1 = "1fc7d880dc104e05dfdc7f8b96a6c1f2486c6228d13a0f04d4e10b15f0c77e96", + Qx0 = "284e28ea45121a3ec0fffd4b8d3f0c470eff76c914e6fc7527f6035cc2a4bf12", + Qx1 = "6dd235191328f1a7245b968396a04e3f6a569491bd6bdc651092cbb95ff65c9", + Qy0 = "24823a8704bb36a2f014c93ba78a2152e8f90ebd19e9196a4f61b15b04409fcc", + Qy1 = "303d0b6bf9db34bdf3beb16dcdbff3a0822b6f241d27b06b3ac1f8707941b4e1" + ) + + test( + id = 3, + EC = ECP_SWei_Proj[Fp2[BN254_Snarks]], + Px0 = "1cf3af1d41e89d8df378aa81463a978c021f27f4a48387e74655ce2cf5c1f298", + Px1 = "36553e80e5c7c7360c7a2ae6bf1b8f68eb48804fc7eba7d2f56f09e87bbb0b1", + Py0 = "25f03e551d74b6be3268bf001905dfbe0bcbe43a2d1aac645a3ca8650b52e551", + Py1 = "2c24c71b843695a003dd5657dba745ce44f7708d9c5c4e0fd1f905751724a57a", + Qx0 = "2c2381c01df71c0762db9458cc369d43a7b2a4f28861580d543010959a8790c1", + Qx1 = "15c9e53568fe7c0d260224d5bc179ecd1bc16b09f421ed56609809c5c5c3bf9b", + Qy0 = "1e5d31d9cda2dd3344a585ffa3273fbed22a1fdf33b45025a480f2e4c07c10ec", + Qy1 = "a8e13d82a6d1f503ce2437733b55b17452d5cc2ff7219f684f343b9c4b09a81" + ) + + # -------------------------------------------------------------------------- + + test( + id = 0, + EC = ECP_SWei_Proj[Fp2[BLS12_377]], + Px0 = "112de13b7cd42bccdb005f2d4dc2726f360243103335ef6cf5e217e777554ae7c1deff5ddb5bcbb581fc9f13728a439", + Px1 = "10d1a8963e5c6854d5e610ece9914f9b5619c27652be1e9ec3e87687d63ed5d45b449bf59c2481e18ac6159f75966ac", + Py0 = "8aaf3a8660cf0edd6e97a2cd7837af1c63ec89e18f9bf4c64638662a661636b928a4f8097e6a2e8dfa11e13c51b075", + Py1 = "163eeb32f275bc5e17546382180b0baefeea482d4da1f7d4938670c66167c7912f571ab3e0426266247b102f8351b3c", + Qx0 = "2ffc357b6f63a3a040b9f1113d1806d35897abcc38fc7617354b9ea834f4c66dcd87e459ab6cafdcdfe2ae44f8bc5", + Qx1 = "17f1a16aa1cfead79134b075300cb5999015f4314d82656c04871289a476451221adeb202754a4d21a57fb03f5c39aa", + Qy0 = "1904a3f203c94c832388cda6c10aa38243c44ae61ee31472d9197d9668e37d58cc6f2181004a9520b27bddc7e523e3a", + Qy1 = "98505fd21506437c605d28f809bd6215431154fa8175f0eaa02928130437e44840acb284f09bf3973350dac32a6dae" + ) + + test( + id = 1, + EC = ECP_SWei_Proj[Fp2[BLS12_377]], + Px0 = "2f9318360b53c2d706061f527571e91679e6086a72ce8203ba1a04850f83bb192b29307e9b2d63feb1d23979e3f632", + Px1 = "3cbab0789968a3a35fa5d2e2326baa40c34d11a4af05a4109350944300ce32eef74dc5e47ba46717bd8bf87604696d", + Py0 = "14ea84922f76f2681fec869dce26141392975dcdb4f21d5fa8aec06b37bf71ba6249c219ecbaef4a266196dafb4ad19", + Py1 = "187cac5daa215b608daab087a9c5ba4364424bb4770c4c5e33112efe931c8a87253f90db38948f3094eb71f3ba593e5", + Qx0 = "3fdbf071e73c2b81d6c3bb7d0deb03460a6fbc13d488644023a16fa0e7bc992f9304f62c37cbd67af0d7f5ef00891c", + Qx1 = "14c3aa3f9ee2a1a9dfd629611907842444fa49212be5732f5594f26c6bcc455d51ba78f593d344729e1eb1a32646149", + Qy0 = "150af74166adf3ae9645178b525e3db9d0141c1b39d9259cf17ba066ea71a3b1163a8d1299a4afce6ca4ee5eb2b9ec3", + Qy1 = "11a869dea20407de3dcfa6288d54edc4afe799c84f8308a7fd9164ea8edcc12d5eda436aafc82cd4ada3b73c9481568" + ) + + test( + id = 2, + EC = ECP_SWei_Proj[Fp2[BLS12_377]], + Px0 = "833ca23630be463c388ea6cfcff5b0e3b055065702a84310d2c726aee14d9e140cba05be79b5cb0441816d9e8c8370", + Px1 = "264a9755524baac8d9e53b0a45789e9dafcb6b453e965061fcfa20bb12a27d9b9417d5277ae2a499b1cfe567d75e2d", + Py0 = "5b670b9789825e2b48101b5b6e660cf9117e29c521dad54640cb356b674b3946c98cb43909c3495fb6d6d231891b7e", + Py1 = "8d794bfd3f87b76ef28af168999e89e6b4fe95da0a539e94a0d0215a7abb756c4b479de5d05a950edf720fd0a20d2d", + Qx0 = "14ab7d6cb634f741384133209d3944cd455ba6abfb7568052f79dedce2a39bbb3a7e4f038a23e8c1b186444997b2f3f", + Qx1 = "15c2626a4b919778485be345bc044717616ee7b9b97c111c2082cf0f249fc946d0aebd7d523234afd5ee40549a61978", + Qy0 = "a783f13eeda1f3d58229008ef4b131135cf5363963dbb712be14186d83e281452b1e0a257ab4851feee7efa150247f", + Qy1 = "11ece9f651b6ca92fe75b277fe3d98ef63cf1a8d467e786be66bb8a705c9cdf2208c715f57540f72b38a3f4cedfc066" + ) + + test( + id = 3, + EC = ECP_SWei_Proj[Fp2[BLS12_377]], + Px0 = "14cd89e2e2755ddc086f63fd62e1f9904c3c1497243455c578a963e81b389f04e95ceafc4f47dc777579cdc82eca79b", + Px1 = "ba8801beba0654f20ccb78783efa7a911d182ec0eb99abe10f9a3d26b46fb7f90552e4ff6beb4df4611a9072be648b", + Py0 = "12e23bc97d891f2a047bac9c90e728cb89760c812156f96c95e36c40f1c830cf6ecbb5d407b189070d48a92eb461ea6", + Py1 = "3b5b911592b3b4110f1690afc0c334b15dccb0eaf6d68b4361c19dfad31e55bbdc219e4328026dc31a4ec122235579", + Qx0 = "151806c45fdb293c57f07e905ea7cbc6f3c9d3803103f7ffec8b499365925f80e63a648a02523e504238ba1417aeef3", + Qx1 = "159018bfc69cf11ac1666efdace397117d0016e70fdff1f095ea5a5340b8e1d2e2a99ade477c8cd56184994ff8860b4", + Qy0 = "1756e93f4b477ab09524a52e4a4551f8ffc8216fc29ec04fbdc31a236750535bd91f369171f763f73196cc5d8d1ea52", + Qy1 = "17aad3c2d6c1918fd733c25be872c445a7afea2191a90b451f85edf7e44a8f9ae4012f7388dad8fc5ccb8e0cb0967b2" + ) + + # -------------------------------------------------------------------------- + + test( + id = 0, + EC = ECP_SWei_Proj[Fp2[BLS12_381]], + Px0 = "d6904be428a0310dbd6e15a744a774bcf9800abe27536267a5383f1ddbd7783e1dc20098a8e045e3cca66b83f6d7f0f", + Px1 = "12107f6ef71d0d1e3bcba9e00a0675d3080519dd1b6c086bd660eb2d2bca8f276e283a891b5c0615064d7886af625cf2", + Py0 = "c592a3546d2d61d671070909e97860822db0a389e351c1744bdbb2c472cf52f3ca3e94068b0b6f3b0121923659131f5", + Py1 = "f3c8aef3a00761f30948689a45dfa0d48ccda74981147e3e8f4877e1784c6bec49e180be98a139e2ed9dcd36ea31c67", + Qx0 = "1030dc4b36a818bdafe99b645008e815abe7e75826d62067787cd76e300a1a7a40e7b53d8b4e74fe78e6d9812376d294", + Qx1 = "17fd5bb7be865765282839efa48fecbcab8be35711d2c878e2d9e874f3328dced52d2a55b1f305f231601f8e700ee6bd", + Qy0 = "b3359ac6f39e703ceabd522c9e472c60db04086bbd8d446f11c2ebb7fb578dba1b3e064207d26ca17270fc46a008fac", + Qy1 = "ce5fe60372147ef37f6fb2ce0fbadb1d1e911bacd1d7ffb367183a4ec475fa93552f3016942c5d0838c572d2ae6c32b" + ) + + test( + id = 1, + EC = ECP_SWei_Proj[Fp2[BLS12_381]], + Px0 = "112de130b7cd42bccdb005f2d4dc2726f360243103335ef6cf5e217e777554ae7c1deff5ddb5bcbb581fc9f13728a439", + Px1 = "10d1a89a63e5c6854d5e610ece9914f9b5619c27652be1e9ec3e87687d63ed5d45b449bf59c2481e18ac6159f75966ac", + Py0 = "11261c8fcb0f4f560479547fe6b2a1c1e8b648d87e54c39f299eba8729294e99b415851d134ca31e8bb861c42e6f1022", + Py1 = "1674a925228b822022bf721c9be8946825b9776c2c06158b330831856d5e05c5a454271d3ada3cd882cd385e2732db4d", + Qx0 = "19a679660d4079c1ed073081951e1a6ce2c06efdd16cb782c227badcf412a61692ba9a89fb23c77a00256cc48505e2b2", + Qx1 = "27442be7b2676bfcda0d655156695f391e6d480b5105b5cd7cee64749bb77b15ce8862d423bfca47f429e0746a5964a", + Qy0 = "ae1d1f52e7b1445112906d8c7c6c7044f08ddd56c24d054e0cd3f3cbc0591a5164ea362c939b09bbdf3357f31dcd93e", + Qy1 = "13662a83e7e9b26bbf5df7e1323a09770eb7345435f23207138012be9acfa6b3ee32da7df4830cb1b1cf00ca4bdccce6" + ) + + test( + id = 2, + EC = ECP_SWei_Proj[Fp2[BLS12_381]], + Px0 = "2f93183360b53c2d706061f527571e91679e6086a72ce8203ba1a04850f83bb192b29307e9b2d63feb1d23979e3f632", + Px1 = "3cbab0c789968a3a35fa5d2e2326baa40c34d11a4af05a4109350944300ce32eef74dc5e47ba46717bd8bf87604696d", + Py0 = "2b8d995b0f2114442b7bbdbe5732fbf94430d6d413e1f388031f3abb956e598cb6764275a75832c1670868c458378b6", + Py1 = "63743343c86cd84b8396413c23fd851144d607cef8b39a1eeb4dc8d4420739574e0be8c498cc26d74096201f7c40580", + Qx0 = "1482a7b3dbdd5810c2167cfea50f6fcc076a23d78f9471205f056ce53fc93266a96969f2c627f8da2190760163d69d4e", + Qx1 = "1106e16aa7dd71959f6a1469de3c5978e6f2177e668faf3d8616ee56a9729b671839754290cece62056ffa9bef7aec55", + Qy0 = "6826a449ed06645081c1d9212231c4f63c98398a514d2950712bb0f4f9dfd0f83327f8623b2d9f32a9d9d3f891519cf", + Qy1 = "19424b599275b17e6c8cfdb9c08100a30a64a0b961e91bf2a6ecb5bcd4e698fa4ef823386d488bd368213dc9b23d2ebd" + ) + + test( + id = 3, + EC = ECP_SWei_Proj[Fp2[BLS12_381]], + Px0 = "d7d1c55ddf8bd03b7a15c3ea4f8f69aee37bf282d4aac82b7bd1fd47139250b9c708997a7ff8f603e48f0471c2cfe03", + Px1 = "d145a91934a6ad865d24ab556ae1e6c42decdd05d676b80e53365a6ff7536332859c9682e7200e40515f675415d71a3", + Py0 = "6de67fa12af93813a42612b1e9449c7b1f160c5de004ec26ea61010e48ba38dcf158d2692f347fdc6c6332bbec7106f", + Py1 = "8177897afa93cc9aca82eb9df71652ef7360ce3ef484041450ae14b090246f1bb46293d658540c7652df8ac48842c0c", + Qx0 = "175a00d5aedc82a684a1af83cedb126112af904f9d84c1c2b0175f739a2805f720c9b97454c57e3c09b8231e9fe5ab7b", + Qx1 = "fbee1da493f45154e2a9bc43b15386757c8dc0db390f9a00f1ccc49ebb14838fa391b8d10b8369294885deb2ae033ed", + Qy0 = "169027ec8c608a49325edc87b726cec5f768acf48e80bba04d2810a7259d6b834043e65e1775c4d9181aecc4ff1430ec", + Qy1 = "77ef6850d4a8f181a10196398cd344011a44c50dce00e18578f3526301263492086d44c7c3d1db5b12499b4033116e1" + ) + +suite "ψ - psi(psi(P)) == psi2(P) - (Untwist-Frobenius-Twist Endomorphism)" & " [" & $WordBitwidth & "-bit mode]": + const Iters = 8 + proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) = + for i in 0 ..< Iters: + let P = rng.random_point(EC, randZ, gen) + + var Q1 {.noInit.}: EC + Q1.frobenius_psi(P) + Q1.frobenius_psi(Q1) + + var Q2 {.noInit.}: EC + Q2.frobenius_psi2(P) + + doAssert bool(Q1 == Q2), "\nIters: " & $i & "\n" & + "P: " & P.toHex() & "\n" & + "Q1: " & Q1.toHex() & "\n" & + "Q2: " & Q2.toHex() + + proc testAll(EC: typedesc) = + test "psi(psi(P)) == psi2(P) for " & $EC: + test(EC, randZ = false, gen = Uniform) + test(EC, randZ = true, gen = Uniform) + test(EC, randZ = false, gen = HighHammingWeight) + test(EC, randZ = true, gen = HighHammingWeight) + test(EC, randZ = false, gen = Long01Sequence) + test(EC, randZ = true, gen = Long01Sequence) + + testAll(ECP_SWei_Proj[Fp2[BN254_Snarks]]) + # testAll(ECP_SWei_Proj[Fp2[BLS12_377]]) + testAll(ECP_SWei_Proj[Fp2[BLS12_381]]) + +suite "ψ²(P) - [t]ψ(P) + [p]P = Inf" & " [" & $WordBitwidth & "-bit mode]": + const Iters = 10 + proc trace(C: static Curve): auto = + # Returns (abs(trace), isNegativeSign) + when C == BN254_Snarks: + # x = "0x44E992B44A6909F1" + # t = 6x²+1 + return (BigInt[127].fromHex"0x6f4d8248eeb859fbf83e9682e87cfd47", false) + elif C == BLS12_381: + # x = "-(2^63 + 2^62 + 2^60 + 2^57 + 2^48 + 2^16)" + # t = x+1 + return (BigInt[64].fromHex"0xd20100000000ffff", true) + else: + {.error: "Not implemented".} + + proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) = + let trace = trace(EC.F.C) + + for i in 0 ..< Iters: + let P = rng.random_point(EC, randZ, gen) + + var r {.noInit.}, psi2 {.noInit.}, tpsi {.noInit.}, pP {.noInit.}: EC + + psi2.frobenius_psi2(P) + tpsi.frobenius_psi(P) + tpsi.scalarMul(trace[0]) # Should be valid for GLS scalar mul even if cofactor isn't cleared + if trace[1]: # negative trace + tpsi.neg() + pP = P + pP.scalarMul(EC.F.C.Mod) # Should be valid for GLS scalar mul even if cofactor isn't cleared + + # ψ²(P) - [t]ψ(P) + [p]P = InfinityPoint + r.diff(psi2, tpsi) + r += pP + + doAssert bool(r.isInf()) + + proc testAll(EC: typedesc) = + test "ψ²(P) - [t]ψ(P) + [p]P = Inf for " & $EC: + test(EC, randZ = false, gen = Uniform) + test(EC, randZ = true, gen = Uniform) + test(EC, randZ = false, gen = HighHammingWeight) + test(EC, randZ = true, gen = HighHammingWeight) + test(EC, randZ = false, gen = Long01Sequence) + test(EC, randZ = true, gen = Long01Sequence) + + testAll(ECP_SWei_Proj[Fp2[BN254_Snarks]]) + # testAll(ECP_SWei_Proj[Fp2[BLS12_377]]) + testAll(ECP_SWei_Proj[Fp2[BLS12_381]]) + +suite "ψ⁴(P) - ψ²(P) + P = Inf (k-th cyclotomic polynomial with embedding degree k=12)" & " [" & $WordBitwidth & "-bit mode]": + const Iters = 10 + + proc test(EC: typedesc, randZ: static bool, gen: static RandomGen) = + for i in 0 ..< Iters: + let P = rng.random_point(EC, randZ, gen) + + var r {.noInit.}, psi4 {.noInit.}, psi2 {.noInit.}: EC + + psi2.frobenius_psi2(P) + psi4.frobenius_psi2(psi2) + r.diff(psi4, psi2) + r += P + + doAssert bool(r.isInf()) + + proc testAll(EC: typedesc) = + test "ψ⁴(P) - ψ²(P) + P = Inf for " & $EC: + test(EC, randZ = false, gen = Uniform) + test(EC, randZ = true, gen = Uniform) + test(EC, randZ = false, gen = HighHammingWeight) + test(EC, randZ = true, gen = HighHammingWeight) + test(EC, randZ = false, gen = Long01Sequence) + test(EC, randZ = true, gen = Long01Sequence) + + testAll(ECP_SWei_Proj[Fp2[BN254_Snarks]]) + # testAll(ECP_SWei_Proj[Fp2[BLS12_377]]) + testAll(ECP_SWei_Proj[Fp2[BLS12_381]]) diff --git a/tests/t_ec_sage_bn254.nim b/tests/t_ec_sage_bn254.nim index d2094a1..ef981d0 100644 --- a/tests/t_ec_sage_bn254.nim +++ b/tests/t_ec_sage_bn254.nim @@ -8,7 +8,7 @@ import # Standard library - std/[unittest, times], + std/unittest, # Internals ../constantine/config/[common, curves], ../constantine/arithmetic, diff --git a/tests/t_ec_template.nim b/tests/t_ec_template.nim index e78a66a..bb8fdd6 100644 --- a/tests/t_ec_template.nim +++ b/tests/t_ec_template.nim @@ -26,12 +26,12 @@ import ./support/ec_reference_scalar_mult type - RandomGen = enum + RandomGen* = enum Uniform HighHammingWeight Long01Sequence -func random_point(rng: var RngState, F: typedesc, randZ: static bool, gen: static RandomGen): F {.inline, noInit.} = +func random_point*(rng: var RngState, F: typedesc, randZ: static bool, gen: static RandomGen): F {.inline, noInit.} = when not randZ: when gen == Uniform: result = rng.random_unsafe(F)