feat(secp256k1): add endomorphism acceleration (#444)

benchmarks/bench_ec_g1.nim

@@ -34,7 +34,7 @@ const AvailableCurves = [
   BN254_Snarks,
   # Edwards25519,
   # P256,
-  # Secp256k1,
+  Secp256k1,
   Pallas,
   Vesta,
   BLS12_377,

benchmarks/bench_ec_g1_scalar_mul.nim

@@ -31,7 +31,7 @@ const AvailableCurves = [
   BN254_Snarks,
   # Edwards25519,
   # P256,
-  # Secp256k1,
+  Secp256k1,
   Pallas,
   Vesta,
   BLS12_377,
@@ -67,7 +67,7 @@ proc main() =
     scalarMulVartimeWNAFBench(EC_ShortW_Jac[Fp[curve], G1], bits, window = 4, MulIters)
     scalarMulVartimeWNAFBench(EC_ShortW_Jac[Fp[curve], G1], bits, window = 5, MulIters)
     separator()
-    when bits >= EndomorphismThreshold: # All endomorphisms constants are below this threshold
+    when curve.hasEndomorphismAcceleration():
       scalarMulVartimeEndoWNAFBench(EC_ShortW_Prj[Fp[curve], G1], bits, window = 2, MulIters)
       scalarMulVartimeEndoWNAFBench(EC_ShortW_Prj[Fp[curve], G1], bits, window = 3, MulIters)
       scalarMulVartimeEndoWNAFBench(EC_ShortW_Prj[Fp[curve], G1], bits, window = 4, MulIters)

benchmarks/bench_ec_g2_scalar_mul.nim

@@ -66,7 +66,7 @@ proc main() =
     scalarMulVartimeWNAFBench(EC_ShortW_Jac[Fp2[curve], G2], bits, window = 4, MulIters)
     scalarMulVartimeWNAFBench(EC_ShortW_Jac[Fp2[curve], G2], bits, window = 5, MulIters)
     separator()
-    when bits >= EndomorphismThreshold: # All endomorphisms constants are below this threshold
+    when curve.hasEndomorphismAcceleration():
       scalarMulVartimeEndoWNAFBench(EC_ShortW_Prj[Fp2[curve], G2], bits, window = 2, MulIters)
       scalarMulVartimeEndoWNAFBench(EC_ShortW_Prj[Fp2[curve], G2], bits, window = 3, MulIters)
       scalarMulVartimeEndoWNAFBench(EC_ShortW_Prj[Fp2[curve], G2], bits, window = 4, MulIters)

benchmarks/bench_elliptic_parallel_template.nim

@@ -56,14 +56,14 @@ proc multiAddParallelBench*(EC: typedesc, numInputs: int, iters: int) =
 type BenchMsmContext*[EC] = object
   tp: Threadpool
   numInputs: int
-  coefs: seq[getBigInt(EC.getName(), kScalarField)]
+  coefs: seq[BigInt[64]] # seq[getBigInt(EC.getName(), kScalarField)]
   points: seq[affine(EC)]
 
 proc createBenchMsmContext*(EC: typedesc, inputSizes: openArray[int]): BenchMsmContext[EC] =
   result.tp = Threadpool.new()
   let maxNumInputs = inputSizes.max()
 
-  const bits = EC.getScalarField().bits()
+  const bits = 64 # EC.getScalarField().bits()
   type ECaff = affine(EC)
 
   result.numInputs = maxNumInputs
@@ -104,7 +104,7 @@ proc createBenchMsmContext*(EC: typedesc, inputSizes: openArray[int]): BenchMsmC
   stdout.write &"in {float64(inNanoSeconds(stop-start)) / 1e6:6.3f} ms\n"
 
 proc msmParallelBench*[EC](ctx: var BenchMsmContext[EC], numInputs: int, iters: int) =
-  const bits = EC.getScalarField().bits()
+  const bits = 64 # EC.getScalarField().bits()
   type ECaff = affine(EC)
 
   template coefs: untyped = ctx.coefs.toOpenArray(0, numInputs-1)

constantine.nimble

@@ -529,6 +529,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/math_elliptic_curves/t_ec_sage_bls12_381.nim", false),
   ("tests/math_elliptic_curves/t_ec_sage_pallas.nim", false),
   ("tests/math_elliptic_curves/t_ec_sage_vesta.nim", false),
+  ("tests/math_elliptic_curves/t_ec_sage_secp256k1.nim", false),
 
   # Edge cases highlighted by past bugs
   # ----------------------------------------------------------

constantine/math/arithmetic/finite_fields.nim

@@ -472,6 +472,15 @@ func `*=`*(a: var FF, b: static int) =
     a.double(t)  # 6
     a.double()   # 12
     a += t       # 15
+  elif b == 21:
+    var t {.noInit.}: typeof(a)
+    t.double(a)
+    t.double()   # 4
+    t += a       # 5
+    t.double()   # 10
+    t.double()   # 20
+    a += t       # 21
+
   else:
     {.error: "Multiplication by this small int not implemented".}
 

constantine/math/elliptic/ec_multi_scalar_mul_parallel.nim

@@ -520,7 +520,13 @@ template withEndo[coefsBits: static int, EC, ECaff](
            coefs: ptr UncheckedArray[BigInt[coefsBits]],
            points: ptr UncheckedArray[ECaff],
            N: int, c: static int, useParallelBuckets: static bool) =
-  when coefsBits <= EC.getScalarField().bits() and hasEndomorphismAcceleration(EC.getName()):
+  when hasEndomorphismAcceleration(EC.getName()) and
+        EndomorphismThreshold <= coefsBits and
+        coefsBits <= EC.getScalarField().bits() and
+        # computeEndomorphism assumes they can be applied to affine repr
+        # but this is not the case for Bandersnatch/wagon
+        # instead Twisted Edwards MSM should be overloaded for Projective/ProjectiveExtended
+        EC.getName() notin {Bandersnatch, Banderwagon}:
     let (endoCoefs, endoPoints, endoN) = applyEndomorphism_parallel(tp, coefs, points, N)
     # Given that bits and N changed, we are able to use a bigger `c`
     # but it has no significant impact on performance

constantine/named/constants/secp256k1_endomorphisms.nim

@@ -0,0 +1,31 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  constantine/named/algebras,
+  constantine/math/io/[io_bigints, io_fields]
+
+# Secp256k1 G1
+# ------------------------------------------------------------
+
+const Secp256k1_cubicRootOfUnity_mod_p* =
+  Fp[Secp256k1].fromHex"0x851695d49a83f8ef919bb86153cbcb16630fb68aed0a766a3ec693d68e6afa40"
+
+const Secp256k1_Lattice_G1* = (
+  # (BigInt, isNeg)
+  ((BigInt[128].fromHex"0xe4437ed6010e88286f547fa90abfe4c3", false),
+   (BigInt[126].fromHex"0x3086d221a7d46bcde86c90e49284eb15", true)),
+  ((BigInt[126].fromHex"0x3086d221a7d46bcde86c90e49284eb15", false),
+   (BigInt[129].fromHex"0x114ca50f7a8e2f3f657c1108d9d44cfd8", false))
+)
+
+const Secp256k1_Babai_G1* = (
+  # (BigInt, isNeg)
+  (BigInt[129].fromHex"0x114ca50f7a8e2f3f657c1108d9d44cfd9", false),
+  (BigInt[126].fromHex"0x3086d221a7d46bcde86c90e49284eb15", false)
+)

constantine/named/zoo_endomorphisms.nim

@@ -27,6 +27,7 @@ import
   ./constants/bw6_761_endomorphisms,
   ./constants/pallas_endomorphisms,
   ./constants/vesta_endomorphisms,
+  ./constants/secp256k1_endomorphisms,
   ./constants/bandersnatch_endomorphisms,
   ./constants/banderwagon_endomorphisms
 
@@ -112,6 +113,7 @@ func hasEndomorphismAcceleration*(Name: static Algebra): bool {.compileTime.} =
     Banderwagon,
     BN254_Nogami,
     BN254_Snarks,
+    Secp256k1,
     BLS12_377,
     BLS12_381,
     BW6_761,

sage/curves.sage

@@ -174,5 +174,16 @@ Curves = {
       'a': 0,
       'b': 5
     }
-  }
+  },
+  'Secp256k1': {
+    'field': {
+      'modulus':  Integer('0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f'),
+      'order': Integer('0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141'),
+    },
+    'curve': {
+      'form': 'short_weierstrass',
+      'a': 0,
+      'b': 7
+    }
+  },
 }

sage/derive_endomorphisms.sage

@@ -258,16 +258,20 @@ if __name__ == "__main__":
     print('\nPrecomputing G1 - 𝜑 (phi) cubic root endomorphism')
     print('----------------------------------------------------\n')
     cubeRootModP, g1lat, g1babai = genCubicRootEndo(curve, Curves)
-    print('\n\nPrecomputing 𝔾₂ - ψ (Psi) - untwist-Frobenius-twist endomorphism')
-    print('----------------------------------------------------\n')
-    g2lat, g2babai = genPsiEndo(curve, Curves)
+
+    hasG2 = 'tower' in Curves[curve]
+
+    if hasG2:
+        print('\n\nPrecomputing 𝔾₂ - ψ (Psi) - untwist-Frobenius-twist endomorphism')
+        print('----------------------------------------------------\n')
+        g2lat, g2babai = genPsiEndo(curve, Curves)
 
     with open(f'{curve.lower()}_endomorphisms.nim', 'w') as f:
       f.write(copyright())
       f.write('\n\n')
       f.write(inspect.cleandoc(f"""
         import
-          constantine/named/algebra,
+          constantine/named/algebras,
           constantine/math/io/[io_bigints, io_fields]
 
         # {curve} G1
@@ -288,18 +292,21 @@ if __name__ == "__main__":
         f'{curve}_Babai_G1',
         dumpBabai(g1babai)
       ))
-      f.write('\n\n')
-      f.write(inspect.cleandoc(f"""
-        # {curve} 𝔾₂
-        # ------------------------------------------------------------
-      """))
-      f.write('\n\n')
-      f.write(dumpConst(
-        f'{curve}_Lattice_G2',
-        dumpLattice(g2lat)
-      ))
       f.write('\n')
-      f.write(dumpConst(
-        f'{curve}_Babai_G2',
-        dumpBabai(g2babai)
-      ))
+      if hasG2:
+        f.write('\n')
+        f.write(inspect.cleandoc(f"""
+            # {curve} 𝔾₂
+            # ------------------------------------------------------------
+        """))
+        f.write('\n\n')
+        f.write(dumpConst(
+            f'{curve}_Lattice_G2',
+            dumpLattice(g2lat)
+        ))
+        f.write('\n')
+        f.write(dumpConst(
+            f'{curve}_Babai_G2',
+            dumpBabai(g2babai)
+        ))
+        f.write('\n')

sage/testgen_scalar_mul.sage

@@ -263,5 +263,5 @@ if __name__ == "__main__":
     elif group == 'G2':
       out = genScalarMulG2(curve, Curves, count, seed, scalarBits)
 
-    with open(f'tv_{curve}_scalar_mul_{group}_{bits}bits.json', 'w') as f:
+    with open(f'tv_{curve}_scalar_mul_{group}_{bits}bit.json', 'w') as f:
       json.dump(out, f, indent=2)

tests/math_elliptic_curves/t_ec_sage_secp256k1.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
+  # Internals
+  constantine/named/algebras,
+  constantine/math/extension_fields,
+  constantine/math/ec_shortweierstrass,
+  # Test utilities
+  ./t_ec_sage_template
+
+staticFor(bits, [32, 64, 128, Fr[Secp256k1].bits()]):
+  run_scalar_mul_test_vs_sage(
+    EC_ShortW_Prj[Fp[Secp256k1], G1], bits,
+    "t_ec_sage_secp256k1_g1_projective"
+  )
+
+  run_scalar_mul_test_vs_sage(
+    EC_ShortW_Jac[Fp[Secp256k1], G1], bits,
+    "t_ec_sage_secp256k1_g1_jacobian"
+  )

tests/math_elliptic_curves/t_ec_shortw_jac_g1_add_double.nim

@@ -22,6 +22,12 @@ run_EC_addition_tests(
     moduleName = "test_ec_shortweierstrass_jacobian_g1_add_double_" & $BN254_Snarks
   )
 
+run_EC_addition_tests(
+    ec = EC_ShortW_Jac[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_jacobian_g1_add_double_" & $Secp256k1
+  )
+
 run_EC_addition_tests(
     ec = EC_ShortW_Jac[Fp[BLS12_381], G1],
     Iters = Iters,
@@ -58,6 +64,12 @@ run_EC_addition_vartime_tests(
     moduleName = "test_ec_shortweierstrass_jacobian_g1_add_double_vartime_" & $BN254_Snarks
   )
 
+run_EC_addition_vartime_tests(
+    ec = EC_ShortW_Jac[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_jacobian_g1_add_double_vartime_" & $Secp256k1
+  )
+
 run_EC_addition_vartime_tests(
     ec = EC_ShortW_Jac[Fp[BLS12_381], G1],
     Iters = Iters,

tests/math_elliptic_curves/t_ec_shortw_jac_g1_mixed_add.nim

@@ -23,6 +23,12 @@ run_EC_mixed_add_impl(
     moduleName = "test_ec_shortweierstrass_jacobian_mixed_add_" & $BN254_Snarks
   )
 
+run_EC_mixed_add_impl(
+    ec = EC_ShortW_Jac[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_jacobian_mixed_add_" & $Secp256k1
+  )
+
 run_EC_mixed_add_impl(
     ec = EC_ShortW_Jac[Fp[BLS12_381], G1],
     Iters = Iters,

tests/math_elliptic_curves/t_ec_shortw_jac_g1_mul_distri.nim

@@ -23,6 +23,12 @@ run_EC_mul_distributive_tests(
     moduleName = "test_ec_shortweierstrass_jacobian_g1_mul_distributive_" & $BN254_Snarks
   )
 
+run_EC_mul_distributive_tests(
+    ec = EC_ShortW_Jac[Fp[Secp256k1], G1],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_shortweierstrass_jacobian_g1_mul_distributive_" & $Secp256k1
+  )
+
 run_EC_mul_distributive_tests(
     ec = EC_ShortW_Jac[Fp[BLS12_381], G1],
     ItersMul = ItersMul,

tests/math_elliptic_curves/t_ec_shortw_jac_g1_mul_sanity.nim

@@ -68,6 +68,12 @@ suite "Order checks on BN254_Snarks":
       bool not ay.isSquare()
       bool not ay.sqrt_if_square()
 
+run_EC_mul_sanity_tests(
+    ec = EC_ShortW_Jac[Fp[Secp256k1], G1],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_shortweierstrass_jacobian_g1_mul_sanity_" & $Secp256k1
+  )
+
 run_EC_mul_sanity_tests(
     ec = EC_ShortW_Jac[Fp[BLS12_381], G1],
     ItersMul = ItersMul,

tests/math_elliptic_curves/t_ec_shortw_jac_g1_mul_vs_ref.nim

@@ -23,6 +23,12 @@ run_EC_mul_vs_ref_impl(
     moduleName = "test_ec_shortweierstrass_jacobian_g1_mul_vs_ref_" & $BN254_Snarks
   )
 
+run_EC_mul_vs_ref_impl(
+    ec = EC_ShortW_Jac[Fp[Secp256k1], G1],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_shortweierstrass_jacobian_g1_mul_vs_ref_" & $Secp256k1
+  )
+
 run_EC_mul_vs_ref_impl(
     ec = EC_ShortW_Jac[Fp[BLS12_381], G1],
     ItersMul = ItersMul,

tests/math_elliptic_curves/t_ec_shortw_jacext_g1_add_double.nim

@@ -22,6 +22,12 @@ run_EC_addition_tests(
     moduleName = "test_ec_shortweierstrass_jacobian_extended_g1_add_double_" & $BN254_Snarks
   )
 
+run_EC_addition_tests(
+    ec = EC_ShortW_JacExt[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_jacobian_extended_g1_add_double_" & $Secp256k1
+  )
+
 run_EC_addition_tests(
     ec = EC_ShortW_JacExt[Fp[BLS12_381], G1],
     Iters = Iters,

tests/math_elliptic_curves/t_ec_shortw_jacext_g1_mixed_add.nim

@@ -23,6 +23,12 @@ run_EC_mixed_add_impl(
     moduleName = "test_ec_shortweierstrass_jacobian_extendedmixed_add_" & $BN254_Snarks
   )
 
+run_EC_mixed_add_impl(
+    ec = EC_ShortW_JacExt[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_jacobian_extendedmixed_add_" & $Secp256k1
+  )
+
 run_EC_mixed_add_impl(
     ec = EC_ShortW_JacExt[Fp[BLS12_381], G1],
     Iters = Iters,

tests/math_elliptic_curves/t_ec_shortw_prj_g1_add_double.nim

@@ -22,6 +22,12 @@ run_EC_addition_tests(
     moduleName = "test_ec_shortweierstrass_projective_g1_add_double_" & $BN254_Snarks
   )
 
+run_EC_addition_tests(
+    ec = EC_ShortW_Prj[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_projective_g1_add_double_" & $Secp256k1
+  )
+
 run_EC_addition_tests(
     ec = EC_ShortW_Prj[Fp[BLS12_381], G1],
     Iters = Iters,

tests/math_elliptic_curves/t_ec_shortw_prj_g1_mixed_add.nim

@@ -23,6 +23,12 @@ run_EC_mixed_add_impl(
     moduleName = "test_ec_shortweierstrass_projective_mixed_add_" & $BN254_Snarks
   )
 
+run_EC_mixed_add_impl(
+    ec = EC_ShortW_Prj[Fp[Secp256k1], G1],
+    Iters = Iters,
+    moduleName = "test_ec_shortweierstrass_projective_mixed_add_" & $Secp256k1
+  )
+
 run_EC_mixed_add_impl(
     ec = EC_ShortW_Prj[Fp[BLS12_381], G1],
     Iters = Iters,