Scalar mul tests (#28)

* Add sage script for BN254

* Implement (failing) scalar multiplication tests

* Add a first test against sagemath

* Finish the tests against SAGE for BN254

* Add significant test coverage of scalar multiplication with reference checks for BN254_Snarks and BLS12_381

benchmarks/bench_ec_swei_proj_g1.nim

@@ -27,7 +27,7 @@ import
 
 
 const Iters = 1_000_000
-const InvIters = 1000
+const MulIters = 1000
 const AvailableCurves = [
   # P224,
   # BN254_Nogami,
@@ -51,10 +51,19 @@ proc main() =
     separator()
     doublingBench(ECP_SWei_Proj[Fp[curve]], Iters)
     separator()
+    scalarMulBench(ECP_SWei_Proj[Fp[curve]], scratchSpaceSize = 1 shl 2, MulIters)
+    separator()
+    scalarMulBench(ECP_SWei_Proj[Fp[curve]], scratchSpaceSize = 1 shl 3, MulIters)
+    separator()
+    scalarMulBench(ECP_SWei_Proj[Fp[curve]], scratchSpaceSize = 1 shl 4, MulIters)
+    separator()
+    # scalarMulUnsafeDoubleAddBench(ECP_SWei_Proj[Fp[curve]], MulIters)
+    # separator()
+  separator()
 
 main()
 
-echo "Notes:"
+echo "\nNotes:"
 echo "  - GCC is significantly slower than Clang on multiprecision arithmetic."
 echo "  - The simplest operations might be optimized away by the compiler."
 echo "  - Fast Squaring and Fast Multiplication are possible if there are spare bits in the prime representation (i.e. the prime uses 254 bits out of 256 bits)"

benchmarks/bench_elliptic_template.nim

@@ -15,6 +15,8 @@
 import
   # Internals
   ../constantine/config/curves,
+  ../constantine/arithmetic,
+  ../constantine/io/io_bigints,
   # Helpers
   ../helpers/[prng_unsafe, static_for],
   ./platforms,
@@ -69,12 +71,12 @@ when SupportsGetTicks:
 echo "\n=================================================================================================================\n"
 
 proc separator*() =
-  echo "-".repeat(132)
+  echo "-".repeat(157)
 
 proc report(op, elliptic: string, start, stop: MonoTime, startClk, stopClk: int64, iters: int) =
   let ns = inNanoseconds((stop-start) div iters)
   let throughput = 1e9 / float64(ns)
-  echo &"{op:<15} {elliptic:<40} {throughput:>15.3f} ops/s     {ns:>9} ns/op     {(stopClk - startClk) div iters:>9} CPU cycles (approx)"
+  echo &"{op:<40} {elliptic:<40} {throughput:>15.3f} ops/s     {ns:>9} ns/op     {(stopClk - startClk) div iters:>9} CPU cycles (approx)"
 
 macro fixEllipticDisplay(T: typedesc): untyped =
   # At compile-time, enums are integers and their display is buggy
@@ -108,3 +110,35 @@ proc doublingBench*(T: typedesc, iters: int) =
   let P = rng.random_unsafe(T)
   bench("EC Double G1", T, iters):
     r.double(P)
+
+proc scalarMulBench*(T: typedesc, scratchSpaceSize: static int, iters: int) =
+  const bits = T.F.C.getCurveOrderBitwidth()
+
+  var r {.noInit.}: T
+  let P = rng.random_unsafe(T)
+
+  let exponent = rng.random_unsafe(BigInt[bits])
+  var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+  exponentCanonical.exportRawUint(exponent, bigEndian)
+
+  var scratchSpace{.noInit.}: array[scratchSpaceSize, T]
+
+  bench("EC ScalarMul G1 (scratchsize = " & $scratchSpaceSize & ')', T, iters):
+    r = P
+    r.scalarMul(exponentCanonical, scratchSpace)
+
+# import ../tests/support/ec_reference_scalar_mult
+#
+# proc scalarMulUnsafeDoubleAddBench*(T: typedesc, iters: int) =
+#   const bits = T.F.C.getCurveOrderBitwidth()
+#
+#   var r {.noInit.}: T
+#   let P = rng.random_unsafe(T)
+#
+#   let exponent = rng.random_unsafe(BigInt[bits])
+#   var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+#   exponentCanonical.exportRawUint(exponent, bigEndian)
+#
+#   bench("EC ScalarMul G1 (unsafe DoubleAdd)", T, iters):
+#     r = P
+#     r.unsafe_ECmul_double_add(exponentCanonical)

constantine.nimble

@@ -67,6 +67,8 @@ task test, "Run all tests":
 
   # Elliptic curve arithmetic
   test "", "tests/test_ec_weierstrass_projective_g1.nim"
+  test "", "tests/test_ec_bn254.nim"
+  test "", "tests/test_ec_bls12_381.nim"
 
   if sizeof(int) == 8: # 32-bit tests on 64-bit arch
     # Primitives
@@ -95,6 +97,8 @@ task test, "Run all tests":
 
     # Elliptic curve arithmetic
     test "-d:Constantine32", "tests/test_ec_weierstrass_projective_g1.nim"
+    test "-d:Constantine32", "tests/test_ec_bn254.nim"
+    test "-d:Constantine32", "tests/test_ec_bls12_381.nim"
 
   # Benchmarks compile and run
   # ignore Windows 32-bit for the moment
@@ -131,6 +135,8 @@ task test_no_gmp, "Run tests that don't require GMP":
 
   # Elliptic curve arithmetic
   test "", "tests/test_ec_weierstrass_projective_g1.nim"
+  test "", "tests/test_ec_bn254.nim"
+  test "", "tests/test_ec_bls12_381.nim"
 
   if sizeof(int) == 8: # 32-bit tests
     # Primitives
@@ -155,6 +161,8 @@ task test_no_gmp, "Run tests that don't require GMP":
 
     # Elliptic curve arithmetic
     test "-d:Constantine32", "tests/test_ec_weierstrass_projective_g1.nim"
+    test "-d:Constantine32", "tests/test_ec_bn254.nim"
+    test "-d:Constantine32", "tests/test_ec_bls12_381.nim"
 
   # Benchmarks compile and run
   # ignore Windows 32-bit for the moment

constantine/config/curves.nim

@@ -27,7 +27,7 @@ macro Mod*(C: static Curve): untyped =
   ## Get the Modulus associated to a curve
   result = bindSym($C & "_Modulus")
 
-func getCurveBitSize*(C: static Curve): static int =
+func getCurveBitwidth*(C: static Curve): static int =
   ## Returns the number of bits taken by the curve modulus
   result = static(CurveBitWidth[C])
 
@@ -43,8 +43,22 @@ func family*(C: static Curve): CurveFamily =
 #
 # ############################################################
 
+macro getCurveOrder*(C: static Curve): untyped =
+  ## Get the curve order `r`
+  ## i.e. the number of points on the elliptic curve
+  result = bindSym($C & "_Order")
+
+macro getCurveOrderBitwidth*(C: static Curve): untyped =
+  ## Get the curve order `r`
+  ## i.e. the number of points on the elliptic curve
+  result = nnkDotExpr.newTree(
+    getAST(getCurveOrder(C)),
+    ident"bits"
+  )
+
 macro getEquationForm*(C: static Curve): untyped =
   ## Returns the equation form
+  ## (ShortWeierstrass, Montgomery, Twisted Edwards, Weierstrass, ...)
   result = bindSym($C & "_equation_form")
 
 macro getCoefA*(C: static Curve): untyped =

constantine/config/curves_declaration.nim

@@ -94,6 +94,9 @@ declareCurves:
 
     # G1 Equation: Y^2 = X^3 + 3
     # G2 Equation: Y^2 = X^3 + 3/(9+𝑖)
+    order: "0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001"
+    orderBitwidth: 254
+    cofactor: 1
     eq_form: ShortWeierstrass
     coef_a: 0
     coef_b: 3
@@ -141,6 +144,9 @@ declareCurves:
 
     # G1 Equation: y² = x³ + 4
     # G2 Equation: y² = x³ + 4 (1+i)
+    order: "0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001"
+    orderBitwidth: 255
+    cofactor: "0x396c8c005555e1568c00aaab0000aaab"
     eq_form: ShortWeierstrass
     coef_a: 0
     coef_b: 4

constantine/config/curves_parser.nim

@@ -102,6 +102,8 @@ type
     eq_form: CurveEquationForm
     coef_A: CurveCoef
     coef_B: CurveCoef
+    order: NimNode # nnkStrLit (hex)
+    orderBitwidth: NimNode # nnkIntLit
 
     sexticTwist: SexticTwist
     sexticNonResidue_fp2: NimNode # nnkPar(nnkIntLit, nnkIntLit)
@@ -191,6 +193,12 @@ proc parseCurveDecls(defs: var seq[CurveParams], curves: NimNode) =
           params.coef_B = CurveCoef(kind: Small, coef: sectionVal.intVal.int)
         else:
           params.coef_B = CurveCoef(kind: Large, coefHex: sectionVal.strVal)
+      elif sectionId.eqIdent"order":
+        params.order = sectionVal
+      elif sectionId.eqIdent"orderBitwidth":
+        params.orderBitwidth = sectionVal
+      elif sectionId.eqIdent"cofactor":
+        discard "TODO"
       elif sectionId.eqIdent"nonresidue_quad_fp":
         params.nonresidue_quad_fp = sectionVal
       elif sectionId.eqIdent"nonresidue_cube_fp2":
@@ -269,6 +277,16 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
       exported($curve & "_equation_form"),
       newLit curveDef.eq_form
     )
+    if not curveDef.order.isNil:
+      curveDef.orderBitwidth.expectKind(nnkIntLit)
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_Order"),
+        newCall(
+          bindSym"fromHex",
+          nnkBracketExpr.newTree(bindSym"BigInt", curveDef.orderBitwidth),
+          curveDef.order
+        )
+      )
     if curveDef.coef_A.kind != NoCoef and curveDef.coef_B.kind != NoCoef:
       curveEllipticStmts.add newConstStmt(
         exported($curve & "_coef_A"),

constantine/elliptic/ec_weierstrass_projective.nim

@@ -30,7 +30,7 @@ type ECP_SWei_Proj*[F] = object
   ## corresponding to (x, y) with X = xZ and Y = yZ
   ##
   ## Note that projective coordinates are not unique
-  x, y, z: F
+  x*, y*, z*: F
 
 func `==`*[F](P, Q: ECP_SWei_Proj[F]): SecretBool =
   ## Constant-time equality check
@@ -62,6 +62,14 @@ func setInf*(P: var ECP_SWei_Proj) =
   P.y.setOne()
   P.z.setZero()
 
+func ccopy*(P: var ECP_SWei_Proj, Q: ECP_SWei_Proj, ctl: SecretBool) =
+  ## Constant-time conditional copy
+  ## If ctl is true: Q is copied into P
+  ## if ctl is false: Q is not copied and P is unmodified
+  ## Time and memory accesses are the same whether a copy occurs or not
+  for fP, fQ in fields(P, Q):
+    ccopy(fP, fQ, ctl)
+
 func trySetFromCoordsXandZ*[F](P: var ECP_SWei_Proj[F], x, z: F): SecretBool =
   ## Try to create a point the elliptic curve
   ## Y²Z = X³ + aXZ² + bZ³ (projective coordinates)
@@ -319,7 +327,7 @@ func scalarMulPrologue(
      ): uint =
   ## Setup the scratchspace
   ## Returns the fixed-window size for scalar mul with window optimization
-  result = result.scratchspace.len.getWindowLen()
+  result = scratchspace.len.getWindowLen()
   # Precompute window content, special case for window = 1
   # (i.e scratchspace has only space for 2 temporaries)
   # The content scratchspace[2+k] is set at [k]P
@@ -341,11 +349,11 @@ func scalarMulDoubling(
        window: uint,
        acc, acc_len: var uint,
        e: var int
-     ) =
+     ): tuple[k, bits: uint] {.inline.} =
   ## Doubling steps of doubling and add for scalar multiplication
   ## Get the next k bits in range [1, window)
   ## and double k times
-  ## Returns the niumber of doubling done and the corresponding bits.
+  ## Returns the number of doubling done and the corresponding bits.
   ##
   ## Updates iteration variables and accumulators
   #

constantine/io/io_ec.nim

@@ -0,0 +1,57 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  ./io_bigints, ./io_fields,
+  ../config/curves,
+  ../arithmetic/[bigints, finite_fields],
+  ../elliptic/[
+    ec_weierstrass_affine,
+    ec_weierstrass_projective
+  ]
+
+# No exceptions allowed
+{.push raises: [].}
+{.push inline.}
+
+# ############################################################
+#
+#   Parsing from canonical inputs to internal representation
+#
+# ############################################################
+
+func toHex*(P: ECP_SWei_Proj): string =
+  ## Stringify an elliptic curve point to Hex
+  ## Note. Leading zeros are not removed.
+  ## Result is prefixed with 0x
+  ##
+  ## Output will be padded with 0s to maintain constant-time.
+  ##
+  ## CT:
+  ##   - no leaks
+  ##
+  ## TODO: only normalize and don't display the Z coordinate
+  ##
+  ## This proc output may change format in the future
+  result = $P.F.C & "(x: "
+  result &= P.x.tohex(bigEndian)
+  result &= ", y: "
+  result &= P.y.tohex(bigEndian)
+  result &= ", z: "
+  result &= P.y.tohex(bigEndian)
+  result &= ')'
+
+func fromHex*(dst: var ECP_SWei_Proj, x, y: string): bool {.raises: [ValueError].}=
+  ## Convert hex strings to a curve point
+  ## Returns `false`
+  ## if there is no point with coordinates (`x`, `y`) on the curve
+  ## In that case, `dst` content is undefined.
+  dst.x.fromHex(x)
+  dst.y.fromHex(y)
+  dst.z.setOne()
+  return bool(isOnCurve(dst.x, dst.y))

helpers/prng_unsafe.nim

@@ -82,8 +82,8 @@ func next(rng: var RngState): uint64 =
 # BigInts and Fields
 # ------------------------------------------------------------
 
-func random_unsafe[T](rng: var RngState, a: var T, C: static Curve) {.noInit.}=
-  ## Recursively initialize a BigInt or Field element
+func random_unsafe[T](rng: var RngState, a: var T, C: static Curve) =
+  ## Recursively initialize a BigInt (part of a field) or Field element
   ## Unsafe: for testing and benchmarking purposes only
   when T is BigInt:
     var reduced, unreduced{.noInit.}: T
@@ -99,6 +99,11 @@ func random_unsafe[T](rng: var RngState, a: var T, C: static Curve) {.noInit.}=
     for field in fields(a):
       rng.random_unsafe(field, C)
 
+func random_unsafe(rng: var RngState, a: var BigInt) =
+  ## Initialize a standalone BigInt
+  for i in 0 ..< a.limbs.len:
+    a.limbs[i] = SecretWord(rng.next())
+
 # Elliptic curves
 # ------------------------------------------------------------
 
@@ -172,7 +177,9 @@ func random_unsafe*(rng: var RngState, T: typedesc): T =
     rng.random_unsafe(result)
   elif T is SomeNumber:
     cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
-  else:
+  elif T is BigInt:
+    rng.random_unsafe(result)
+  else: # Fields
     rng.random_unsafe(result, T.C)
 
 func random_unsafe_with_randZ*(rng: var RngState, T: typedesc[ECP_SWei_Proj]): T =