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.}=
+func random_unsafe[T](rng: var RngState, a: var T, C: static Curve) =
-  ## Recursively initialize a BigInt or Field element
+  ## 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 =

sage/testgen_bls12_381.sage

@@ -0,0 +1,69 @@
+# 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.
+
+# ############################################################
+#
+#                    BN254 test generator
+#
+# ############################################################
+
+# Parameters
+u = -(2^63 + 2^62 + 2^60 + 2^57 + 2^48 + 2^16)
+p = (u - 1)^2 * (u^4 - u^2 + 1)//3 + u
+r = u^4 - u^2 + 1
+cofactor = Integer('0x396c8c005555e1568c00aaab0000aaab')
+
+# Finite fields
+F       = GF(p)
+K2.<u>  = PolynomialRing(F)
+# F2.<beta>  = F.extension(u^2+1)
+# K6.<v>  = PolynomialRing(F2)
+# F6.<eta>  = F2.extension(v^3-beta)
+# K12.<w> = PolynomialRing(F6)
+# K12.<gamma> = F6.extension(w^2-eta)
+
+# Curves
+b = 4
+G1 = EllipticCurve(F, [0, b])
+# G2 = EllipticCurve(F2, [0, b*beta])
+
+# Test generator
+set_random_seed(1337)
+
+for i in range(10):
+    print('---------------------------------------')
+    P = G1.random_point()
+    (Px, Py, Pz) = P
+    print('Px: ' + Integer(Px).hex())
+    print('Py: ' + Integer(Py).hex())
+    print('Pz: ' + Integer(Pz).hex())
+    exponent = randrange(r) # Pick an integer below curve order
+    print('scalar: ' + Integer(exponent).hex())
+
+    Q = exponent * P
+    (Qx, Qy, Qz) = Q
+    print('Qx: ' + Integer(Qx).hex())
+    print('Qy: ' + Integer(Qy).hex())
+    print('Qz: ' + Integer(Qz).hex())
+
+
+# CurveOrder sanity check
+#
+# P = G1.random_point()
+# (Px, Py, Pz) = P
+# print('Px: ' + Integer(Px).hex())
+# print('Py: ' + Integer(Py).hex())
+# print('Pz: ' + Integer(Pz).hex())
+#
+# print('order: ' + Integer(r).hex())
+#
+# Q = (r * cofactor) * P
+# (Qx, Qy, Qz) = Q
+# print('Qx: ' + Integer(Qx).hex())
+# print('Qy: ' + Integer(Qy).hex())
+# print('Qz: ' + Integer(Qz).hex())

sage/testgen_bn254_snarks.sage

@@ -0,0 +1,68 @@
+# 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.
+
+# ############################################################
+#
+#                    BN254 test generator
+#
+# ############################################################
+
+# Parameters
+u = Integer('0x44E992B44A6909F1')
+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
+cofactor = 1
+
+# Finite fields
+F       = GF(p)
+K2.<u>  = PolynomialRing(F)
+# F2.<beta>  = F.extension(u^2+9)
+# K6.<v>  = PolynomialRing(F2)
+# F6.<eta>  = F2.extension(v^3-beta)
+# K12.<w> = PolynomialRing(F6)
+# K12.<gamma> = F6.extension(w^2-eta)
+
+# Curves
+b = 3
+G1 = EllipticCurve(F, [0, b])
+# G2 = EllipticCurve(F2, [0, b/beta])
+
+# Test generator
+set_random_seed(1337)
+
+for i in range(10):
+    print('---------------------------------------')
+    P = G1.random_point()
+    (Px, Py, Pz) = P
+    print('Px: ' + Integer(Px).hex())
+    print('Py: ' + Integer(Py).hex())
+    print('Pz: ' + Integer(Pz).hex())
+    exponent = randrange(r) # Pick an integer below curve order
+    print('scalar: ' + Integer(exponent).hex())
+
+    Q = exponent * P
+    (Qx, Qy, Qz) = Q
+    print('Qx: ' + Integer(Qx).hex())
+    print('Qy: ' + Integer(Qy).hex())
+    print('Qz: ' + Integer(Qz).hex())
+
+# CurveOrder sanity check
+#
+# P = G1.random_point()
+# (Px, Py, Pz) = P
+# print('Px: ' + Integer(Px).hex())
+# print('Py: ' + Integer(Py).hex())
+# print('Pz: ' + Integer(Pz).hex())
+#
+# print('order: ' + Integer(r).hex())
+#
+# Q = (r * cofactor) * P
+# (Qx, Qy, Qz) = Q
+# print('Qx: ' + Integer(Qx).hex())
+# print('Qy: ' + Integer(Qy).hex())
+# print('Qz: ' + Integer(Qz).hex())

tests/support/ec_reference_scalar_mult.nim

@@ -0,0 +1,52 @@
+# 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/config/[common, curves],
+  ../../constantine/arithmetic,
+  ../../constantine/elliptic/[ec_weierstrass_affine, ec_weierstrass_projective]
+
+# Support files for testing Elliptic Curve arithmetic
+# ------------------------------------------------------------------------------
+
+iterator unpack(scalarByte: byte): bool =
+  yield bool((scalarByte and 0b10000000) shr 7)
+  yield bool((scalarByte and 0b01000000) shr 6)
+  yield bool((scalarByte and 0b00100000) shr 5)
+  yield bool((scalarByte and 0b00010000) shr 4)
+  yield bool((scalarByte and 0b00001000) shr 3)
+  yield bool((scalarByte and 0b00000100) shr 2)
+  yield bool((scalarByte and 0b00000010) shr 1)
+  yield bool( scalarByte and 0b00000001)
+
+func unsafe_ECmul_double_add*(
+       P: var ECP_SWei_Proj,
+       scalar: openArray[byte],
+     ) =
+  ## **Unsafe** Elliptic Curve Scalar Multiplication
+  ##
+  ##   P <- [k] P
+  ##
+  ## This uses the double-and-add algorithm to verify the constant-time production implementation
+  ## This is UNSAFE to use in production and only intended for testing purposes.
+  ##
+  ## This is highly VULNERABLE to timing attacks and power analysis attacks
+  ##
+  ## `scalar` is in canonical representation in BigEndian (octet string)
+  var t0{.noInit.}, t1{.noInit.}: typeof(P)
+  t0.setInf()
+  t1.setInf()
+  for scalarByte in scalar:
+    for bit in unpack(scalarByte):
+      t1.double(t0)
+      if bit:
+        t0.sum(t1, P)
+      else:
+        t0 = t1
+  P = t0

tests/test_ec_bls12_381.nim

@@ -0,0 +1,151 @@
+# 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
+  unittest, times,
+  # Internals
+  ../constantine/config/[common, curves],
+  ../constantine/arithmetic,
+  ../constantine/io/[io_bigints, io_ec],
+  ../constantine/elliptic/[ec_weierstrass_projective],
+  # Test utilities
+  ./support/ec_reference_scalar_mult
+
+proc test(
+       id: int,
+       EC: typedesc[ECP_SWei_Proj],
+       Px, Py: string,
+       scalar: string,
+       Qx, Qy: string
+     ) =
+
+  test "test " & $id:
+    var P: EC
+    let pOK = P.fromHex(Px, Py)
+    doAssert pOK
+
+    var Q: EC
+    let qOK = Q.fromHex(Qx, Qy)
+
+    let exponent = EC.F.C.matchingBigInt.fromHex(scalar)
+    var exponentCanonical: array[(exponent.bits+7) div 8, byte]
+    exponentCanonical.exportRawUint(exponent, bigEndian)
+
+    var
+      impl = P
+      reference = P
+      scratchSpace: array[1 shl 4, EC]
+
+    impl.scalarMul(exponentCanonical, scratchSpace)
+    reference.unsafe_ECmul_double_add(exponentCanonical)
+
+    doAssert: bool(Q == reference)
+    doAssert: bool(Q == impl)
+
+suite "BLS12_381 implementation (and unsafe reference impl) vs SageMath":
+  # Generated via sage sage/testgen_bls12_381.sage
+  test(
+    id = 1,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "f21eda282230f72b855d48055e68ab3825da87831fa5147a64fa071bade4c26bddd45e8b602e62df4d907414a6ec1b4",
+    Py = "531b38866cb35c19951f4a1ac62242f11fa714a1b99c6116a630fa75e7f4407fcd1ae9770a821c5899a777d341c915a",
+    scalar = "f7e60a832eb77ac47374bc93251360d6c81c21add62767ff816caf11a20d8db",
+    Qx = "18d7ca3fb93d7300a0484233f3bac9bca00b45595a4b9caf66aa0b2237f6fd51559a24a634f3876451332c5f754438b2",
+    Qy = "edbb203999303fc99ef04368412da4b3555f999c703b425dedff3fdc799317c292751c46275b27990c53d933de2db63"
+  )
+
+  test(
+    id = 2,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "9ca8e33d8a330b04b052af6cf44bf2ed08cc93d83a4eb48cbb0cabfe02ffb2ef910df44862b271354352f15b70e45b5",
+    Py = "102f6d07ef45f51de9a4ecef5ec34eae16833f4761c2ddfbe2b414173c3580721135e5bbb74269ab85ba83cb03020d9b",
+    scalar = "5f10367bdae7aa872d90b5ac209321ce5a15181ce22848d032a8d452055cbfd0",
+    Qx = "a50d49e3d8757f994aae312dedd55205687c432bc9d97efbe69e87bef4256b87af1b665a669d06657cda6ff01ee42df",
+    Qy = "160d50aaa21f9d5b4faada77e4f91d8d4f152a0fcca4d30d271d74b20c1bba8638128f99f52d9603d4a24f8e27219bcd"
+  )
+
+  test(
+    id = 3,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "173c28687d23de83c950131e548485e8e09d4053d32b814d13b618ee4159e8b61bf6320148ddabcedf2b04d3c9787cd4",
+    Py = "277f935b4e0a90155915960c617f395dcadead1c7297cf92916add07308fc3f0493aa6dabf31d1f15953f56ac37d3d9",
+    scalar = "4c321d72220c098fc0fd52306de98f8be9446bf854cf1e4d8dbae62375d18faf",
+    Qx = "16259e878b5921bbe1e5672cccea0f29fedbb93b8ce1bae4d4b602b6dd5708c6d4e5d82ff92868828c46fd333aadf82d",
+    Qy = "16d09713f4fe5705f2e3491aa9a1d5827fb3b280f5a1fdde0b01a2b75f5803d528d5f5603cc0e9da29d6a07b8e14df7c"
+  )
+
+  test(
+    id = 4,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "177d32dfa6e97daf5ea8aee520bc5c33b7bee37cba114fda52dd948030ad81abdffbdda402c930791e1048ad531b2c80",
+    Py = "14e9f915698dadd2220586a8bcc921b1f855273c3c0e41a88569e5a9fd2a4e886eeff9a7a11b02ec286987c5a52d55ce",
+    scalar = "1738857afb76c55f615c2a20b44ca90dcb3267d804ec23fddea431dbee4eb37f",
+    Qx = "a4bfcfc65eb16562752f5c164349ef673477e19fe020de84eddbc2958f6d40bbbba39fc67ee8c8fdf007922fec97f79",
+    Qy = "106ccd382d15773e6097f8ea6f012cbec15184d6f4ea08bac2842ed419f0e555f1a43f7434b2e017f9e02971d07eb59d"
+  )
+
+  test(
+    id = 5,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "1bc38e70e62770489063d3b7a3bebbc6735ab8dc389576ff03272c2883045aa051d74f8407920530b4c719b140cda81",
+    Py = "bd24e4fb09ed4098d61e3d2cb456f03d7818ded79dfba9cfe7956829797b12e10f1766c46c1a2e1cf2957295124c782",
+    scalar = "19c47811813444020c999a2b263940b5054cf45bb8ad8e086ff126bfcd5507e1",
+    Qx = "b310d4688f2c9f8cd4c030b62ed27341f4c71341fe9c56858a949a2d51670eb6ebe1339163bdb833e692b0ee0cf4e92",
+    Qy = "c92300561e1acb1e1ae6a1b75f83b9d2d2cb5f07c3f8ea945990ceb75e7ea12c4aec115227c13a05be92f5caed9268e"
+  )
+
+  test(
+    id = 6,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "48cddafaca93d33caf910a7a6f74dc3489d53da9fa2f940b70b6dcf538cc08da1a369809ab86a8ee49cead0ed6bfef6",
+    Py = "173f8dfb384aea011bed89aaca625085dc2940d0775e5f2647fc7574ce822643d0d7b1b39e9a51b9f5a0dca7486bddd0",
+    scalar = "43ffcda71e45a3e90b7502d92b30a0b06c54c95a91aa21e0438677b1c2714ecb",
+    Qx = "ef1e4967a3eb19318a66d092eada9810bebf301c168cea7c73fad9d98f7d4c2bde1071fd142c3da90830509f22a82b5",
+    Qy = "da537922dcb6bf79e4d09237c1a3c5804e3a83b6f18ccb26991d50d77c81bef76139fa73d39c684c7c1616151b1058b"
+  )
+
+  test(
+    id = 7,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "9d56cb273bdeef945078066192b74d2f3077f00f5bd1a50b338c44f7c640005a614f9c6fc89cb4678140b2a721c69a8",
+    Py = "107b42b9a0c22b9e9cd2191b90fede2ab280532ea26806338a5b28533cf9431bde1a8010677a5078c63482953d4f2451",
+    scalar = "64ad0d6c36dba5368e71f0010aebf860288f54611e5aaf18082bae7a404ebfd8",
+    Qx = "e0c78d1e1ed993fdeb14e4872965bc90014aa39c728c457a720bf3123ebcdcb17ac553a619b9b7073ada436565d4bb4",
+    Qy = "c2d9ba441ed90bae4f1597da90e434f1668fda320e4fa04cddcdce0eacb3bc54185d5f7cde826f5bd0e3d59b2424906"
+  )
+
+  test(
+    id = 8,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "150a83a868fa6a74dbc5658445ea99ec47009572f303ce1d3c76370804c5a8c26d40c8b4b35a6585612d704c5fb090cb",
+    Py = "31e73ed0aedebcf0b58d60c16f2e5ddd2d4eb2a6e34177939efcca0767cde241966b5950c3333c62ccddee51de26fe6",
+    scalar = "b0ac3d0e685583075aa46c03a00859dfbec24ccb36e2cae3806d82275adcc03",
+    Qx = "9c5e69fbd492a64e5811af7cc69e42bc14d8626f6d384d3f479d8e06c20ec5f460a1e3839f33899b4a9e0ada876ac6e",
+    Qy = "16990d7d308897c74b87368f847df3ac0bb6609091c8d39b22d5778a4229f0bb92fea385d27db41e237dcfb0d05bd0e7"
+  )
+
+  test(
+    id = 9,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "69498486a06c18f836a8e9ed507bbb563d6d03545e03e08f628e8fbd2e5d098e58950071d516ffd044d92b3a8b07184",
+    Py = "18a169f06fc94f40cd131bdcd23e48e95b1276c0c8daacf56c3a5e278e89ee1094c94aa516113aa4a2455da149f0f989",
+    scalar = "23941bb3c3659423d6fdafb7cff52e0e02de0ac91e64c537c6203d64905b63d0",
+    Qx = "482e085550f5e514dd98f2d9b119c284ac165514d228c8f7a179f2b442968984873223af2255a499dc931c63543c0ba",
+    Qy = "151ce80ca51dd09243d2b1a7937096d6b7494e89190da5ab7604cd913dc4105c871e48c815fefadee2906b8b401e7e71"
+  )
+
+  test(
+    id = 10,
+    EC = ECP_SWei_Proj[Fp[BLS12_381]],
+    Px = "98cc20aa561769b7ee569304503a94752e236bba52938fed7f3093d5867f65361dc8b48c83bd7db490c26736196e20e",
+    Py = "10a68394358903122bd649bd30b473f4d3b4f0830bfe7da1c48ae87d9429d8fd26f5b4be8d8fd8e4214017044696da29",
+    scalar = "4203156dcf70582ea8cbd0388104f47fd5a18ae336b2fed8458e1e4e74d7baf5",
+    Qx = "18ff1dfd96799b7d0bffaa7480121c3a719047815ae41419f1bd1fdd593288bed8827b3d9e45a3a1e01bf7d603b5ba0",
+    Qy = "49b95ca2c0f75dfb15fc07e5692d23f8eb38cb1cc9c48cd0e93a80adbff135a3945cc7a5d53d2b7510d6ee7cf97308d"
+  )

tests/test_ec_bn254.nim

@@ -0,0 +1,151 @@
+# 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
+  unittest, times,
+  # Internals
+  ../constantine/config/[common, curves],
+  ../constantine/arithmetic,
+  ../constantine/io/[io_bigints, io_ec],
+  ../constantine/elliptic/[ec_weierstrass_projective],
+  # Test utilities
+  ./support/ec_reference_scalar_mult
+
+proc test(
+       id: int,
+       EC: typedesc[ECP_SWei_Proj],
+       Px, Py: string,
+       scalar: string,
+       Qx, Qy: string
+     ) =
+
+  test "test " & $id:
+    var P: EC
+    let pOK = P.fromHex(Px, Py)
+    doAssert pOK
+
+    var Q: EC
+    let qOK = Q.fromHex(Qx, Qy)
+
+    let exponent = EC.F.C.matchingBigInt.fromHex(scalar)
+    var exponentCanonical: array[(exponent.bits+7) div 8, byte]
+    exponentCanonical.exportRawUint(exponent, bigEndian)
+
+    var
+      impl = P
+      reference = P
+      scratchSpace: array[1 shl 4, EC]
+
+    impl.scalarMul(exponentCanonical, scratchSpace)
+    reference.unsafe_ECmul_double_add(exponentCanonical)
+
+    doAssert: bool(Q == reference)
+    doAssert: bool(Q == impl)
+
+suite "BN254 implementation (and unsafe reference impl) vs SageMath":
+  # Generated via sage sage/testgen_bn254_snarks.sage
+  test(
+    id = 1,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "22d3af0f3ee310df7fc1a2a204369ac13eb4a48d969a27fcd2861506b2dc0cd7",
+    Py = "1c994169687886ccd28dd587c29c307fb3cab55d796d73a5be0bbf9aab69912e",
+    scalar = "e08a292f940cfb361cc82bc24ca564f51453708c9745a9cf8707b11c84bc448",
+    Qx = "267c05cd49d681c5857124876748365313b9c285e783206f48513ce06d3df931",
+    Qy = "2fa00719ce37465dbe7037f723ed5df08c76b9a27a4dd80d86c0ee5157349b96"
+  )
+
+  test(
+    id = 2,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "2724750abe620fce759b6f18729e40f891a514160d477811a44b222372cc4ea3",
+    Py = "105cdcbe363921790a56bf2696e73642447c60b814827ca4dba86c814912c98a",
+    scalar = "2f5c2960850eabadab1e5595ff0bf841206885653e7f2024248b281a86744790",
+    Qx = "57d2dcbc665fb93fd5119bb982c29700d025423d60a42b5fe17210fd5a868fd",
+    Qy = "2abad564ff78fbc266dfb77bdd110b22271136b33ce5049fb3ca05107787abc"
+  )
+
+  test(
+    id = 3,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "39bc19c41835082f86ca046b71875b051575072e4d6a4aeedac31eee34b07df",
+    Py = "1fdbf42fc20421e1e775fd93ed1888d614f7e39067e7443f21b6a4817481c346",
+    scalar = "29e140c33f706c0111443699b0b8396d8ead339a3d6f3c212b08749cf2a16f6b",
+    Qx = "83895d1c7a2b15a5dfe9371983196591415182978e8ff0e83262e32d768c712",
+    Qy = "2ed8b88e1cd08814ce1d1929d0e4bba6fb5897f915b3525cf12349256da95499"
+  )
+
+  test(
+    id = 4,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "157a3e1ff9dabccced9746e19855a9438098be6d734f07d1c069aa1bd05b8d87",
+    Py = "1c96bf3e48bc1a6635d93d4f1302a0eba39bd907c5d861f2a9d0c714ee60f04d",
+    scalar = "29b05bd55963e262e0fa458c76297fb5be3ec1421fdb1354789f68fdce81dc2c",
+    Qx = "196aeca74447934eeaba0f2263177fcb7eb239985814f8ef2d7bf08677108c9",
+    Qy = "1f5aa4c7df4a9855113c63d8fd55c512c7e919b8ae0352e280bdb1009299c3b2"
+  )
+
+  test(
+    id = 5,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "2f260967d4cd5d15f98c0a0a9d5abaae0c70d3b8d83e1e884586cd6ece395fe7",
+    Py = "2a102c7aebdfaa999d5a99984148ada142f72f5d4158c10368a2e13dded886f6",
+    scalar = "1796de74c1edac90d102e7c33f3fad94304eaff4a67a018cae678774d377f6cd",
+    Qx = "28c73e276807863ecf4ae60b1353790f10f176ca8c55b3db774e33c569ef39d5",
+    Qy = "c386e24828cead255ec7657698559b23a26fc9bd5db70a1fe20b48ecfbd6db9"
+  )
+
+  test(
+    id = 6,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "1b4ccef57f4411360a02b8228e4251896c9492ff93a69ba3720da0cd46a04e83",
+    Py = "1fabcb215bd7c06ead2e6b0167497efc2cdd3dbacf69bcb0244142fd63c1e405",
+    scalar = "116741cd19dac61c5e77877fc6fef40f363b164b501dfbdbc09e17ea51d6beb0",
+    Qx = "192ca2e120b0f5296baf7cc47bfebbbc74748c8847bbdbe485bcb796de2622aa",
+    Qy = "8bc6b1aa4532c727be8fd21a8176d55bc721c727af327f601f7a8dff655b0b9"
+  )
+
+  test(
+    id = 7,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "2807c88d6759280d6bd83a54d349a533d1a66dc32f72cab8114ab707f10e829b",
+    Py = "dbf0d486aeed3d303880f324faa2605aa0219e35661bc88150470c7df1c0b61",
+    scalar = "2a5976268563870739ced3e6efd8cf53887e8e4426803377095708509dd156ca",
+    Qx = "2841f67de361436f64e582a134fe36ab7196334c758a07e732e1cf1ccb35a476",
+    Qy = "21fb9b8311e53832044be5ff024f737aee474bc504c7c158fe760cc999da8612"
+  )
+
+  test(
+    id = 8,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "2754a174a33a55f2a31573767e9bf5381b47dca1cbebc8b68dd4df58b3f1cc2",
+    Py = "f222f59c8893ad87c581dacb3f8b6e7c20e7a13bc5fb6e24262a3436d663b1",
+    scalar = "25d596bf6caf4565fbfd22d81f9cef40c8f89b1e5939f20caa1b28056e0e4f58",
+    Qx = "2b48dd3ace8e403c2905f00cdf13814f0dbecb0c0465e6455fe390cc9730f5a",
+    Qy = "fe65f0cd4ae0d2e459daa4163f32deed1250b5c384eb5aeb933162a41793d25"
+  )
+
+  test(
+    id = 9,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "273bf6c679d8e880034590d16c007bbabc6c65ed870a263b5d1ce7375c18fd7",
+    Py = "2904086cb9e33657999229b082558a74c19b2b619a0499afb2e21d804d8598ee",
+    scalar = "67a499a389129f3902ba6140660c431a56811b53de01d043e924711bd341e53",
+    Qx = "1d827e4569f17f068457ffc52f1c6ed7e2ec89b8b520efae48eff41827f79128",
+    Qy = "be8c488bb9587bcb0faba916277974afe12511e54fbd749e27d3d7efd998713"
+  )
+
+  test(
+    id = 10,
+    EC = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Px = "ec892c09a5f1c68c1bfec7780a1ebd279739383f2698eeefbba745b3e717fd5",
+    Py = "23d273a1b9750fe1d4ebd4b7c25f4a8d7d94f6662c436305cca8ff2cdbd3f736",
+    scalar = "d2f09ceaa2638b7ac3d7d4aa9eff7a12e93dc85db0f9676e5f19fb86d6273e9",
+    Qx = "305d7692b141962a4a92038adfacc0d2691e5589ed097a1c661cc48c84e2b64e",
+    Qy = "bafa230a0f5cc2fa3cf07fa46312cb724fc944b097890fa60f2cf42a1be7963"
+  )

tests/test_ec_weierstrass_projective_g1.nim

@@ -12,11 +12,15 @@ import
   # Internals
   ../constantine/config/[common, curves],
   ../constantine/arithmetic,
+  ../constantine/io/io_bigints,
   ../constantine/elliptic/[ec_weierstrass_affine, ec_weierstrass_projective],
   # Test utilities
-  ../helpers/prng_unsafe
+  ../helpers/prng_unsafe,
+  ./support/ec_reference_scalar_mult
 
-const Iters = 128
+const
+  Iters = 128
+  ItersMul = Iters div 4
 
 var rng: RngState
 let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
@@ -169,3 +173,209 @@ suite "Elliptic curve in Short Weierstrass form y² = x³ + a x + b with project
     test(Fp[BN254_Snarks], randZ = true)
     test(Fp[BLS12_381], randZ = false)
     test(Fp[BLS12_381], randZ = true)
+
+
+  const BN254_Snarks_order_bits = BN254_Snarks.getCurveOrderBitwidth()
+  const BLS12_381_order_bits = BLS12_381.getCurveOrderBitwidth()
+
+  test "EC mul [0]P == Inf":
+    proc test(F: typedesc, bits: static int, randZ: static bool) =
+      for _ in 0 ..< ItersMul:
+        when randZ:
+          let a = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+        else:
+          let a = rng.random_unsafe(ECP_SWei_Proj[F])
+
+        # zeroInit
+        var exponentCanonical: array[(bits+7) div 8, byte]
+
+        var
+          impl = a
+          reference = a
+          scratchSpace{.noInit.}: array[1 shl 4, ECP_SWei_Proj[F]]
+
+        impl.scalarMul(exponentCanonical, scratchSpace)
+        reference.unsafe_ECmul_double_add(exponentCanonical)
+
+        check:
+          bool(impl.isInf())
+          bool(reference.isInf())
+
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = false)
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = true)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = false)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = true)
+
+  test "EC mul [Order]P == Inf":
+    proc test(F: typedesc, bits: static int, randZ: static bool) =
+      for _ in 0 ..< ItersMul:
+        when randZ:
+          let a = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+        else:
+          let a = rng.random_unsafe(ECP_SWei_Proj[F])
+
+        let exponent = F.C.getCurveOrder()
+        var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+        exponentCanonical.exportRawUint(exponent, bigEndian)
+
+        var
+          impl = a
+          reference = a
+          scratchSpace{.noInit.}: array[1 shl 4, ECP_SWei_Proj[F]]
+
+        impl.scalarMul(exponentCanonical, scratchSpace)
+        reference.unsafe_ECmul_double_add(exponentCanonical)
+
+        check:
+          bool(impl.isInf())
+          bool(reference.isInf())
+
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = false)
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = true)
+    # TODO: BLS12 is using a subgroup of order "r" such as r*h = CurveOrder
+    #       with h the curve cofactor
+    #       instead of the full group
+    # test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = false)
+    # test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = true)
+
+  test "EC mul [1]P == P":
+    proc test(F: typedesc, bits: static int, randZ: static bool) =
+      for _ in 0 ..< ItersMul:
+        when randZ:
+          let a = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+        else:
+          let a = rng.random_unsafe(ECP_SWei_Proj[F])
+
+        var exponent{.noInit.}: BigInt[bits]
+        exponent.setOne()
+        var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+        exponentCanonical.exportRawUint(exponent, bigEndian)
+
+        var
+          impl = a
+          reference = a
+          scratchSpace{.noInit.}: array[1 shl 4, ECP_SWei_Proj[F]]
+
+        impl.scalarMul(exponentCanonical, scratchSpace)
+        reference.unsafe_ECmul_double_add(exponentCanonical)
+
+        check:
+          bool(impl == a)
+          bool(reference == a)
+
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = false)
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = true)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = false)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = true)
+
+  test "EC mul [2]P == P.double()":
+    proc test(F: typedesc, bits: static int, randZ: static bool) =
+      for _ in 0 ..< ItersMul:
+        when randZ:
+          let a = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+        else:
+          let a = rng.random_unsafe(ECP_SWei_Proj[F])
+
+        var doubleA{.noInit.}: ECP_SWei_Proj[F]
+        doubleA.double(a)
+
+        let exponent = BigInt[bits].fromUint(2)
+        var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+        exponentCanonical.exportRawUint(exponent, bigEndian)
+
+        var
+          impl = a
+          reference = a
+          scratchSpace{.noInit.}: array[1 shl 4, ECP_SWei_Proj[F]]
+
+        impl.scalarMul(exponentCanonical, scratchSpace)
+        reference.unsafe_ECmul_double_add(exponentCanonical)
+
+        check:
+          bool(impl == doubleA)
+          bool(reference == doubleA)
+
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = false)
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = true)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = false)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = true)
+
+  test "EC mul is distributive over EC add":
+    proc test(F: typedesc, bits: static int, randZ: static bool) =
+      for _ in 0 ..< ItersMul:
+        when randZ:
+          let a = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+          let b = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+        else:
+          let a = rng.random_unsafe(ECP_SWei_Proj[F])
+          let b = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+
+        let exponent = rng.random_unsafe(BigInt[bits])
+        var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+        exponentCanonical.exportRawUint(exponent, bigEndian)
+
+        # [k](a + b) - Factorized
+        var
+          fImpl{.noInit.}: ECP_SWei_Proj[F]
+          fReference{.noInit.}: ECP_SWei_Proj[F]
+          scratchSpace{.noInit.}: array[1 shl 4, ECP_SWei_Proj[F]]
+
+        fImpl.sum(a, b)
+        fReference.sum(a, b)
+
+        fImpl.scalarMul(exponentCanonical, scratchSpace)
+        fReference.unsafe_ECmul_double_add(exponentCanonical)
+
+        # [k]a + [k]b - Distributed
+        var kaImpl = a
+        var kaRef = a
+
+        kaImpl.scalarMul(exponentCanonical, scratchSpace)
+        kaRef.unsafe_ECmul_double_add(exponentCanonical)
+
+        var kbImpl = b
+        var kbRef = b
+
+        kbImpl.scalarMul(exponentCanonical, scratchSpace)
+        kbRef.unsafe_ECmul_double_add(exponentCanonical)
+
+        var kakbImpl{.noInit.}, kakbRef{.noInit.}: ECP_SWei_Proj[F]
+        kakbImpl.sum(kaImpl, kbImpl)
+        kakbRef.sum(kaRef, kbRef)
+
+        check:
+          bool(fImpl == kakbImpl)
+          bool(fReference == kakbRef)
+          bool(fImpl == fReference)
+
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = false)
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = true)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = false)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = true)
+
+  test "EC mul constant-time is equivalent to a simple double-and-add algorithm":
+    proc test(F: typedesc, bits: static int, randZ: static bool) =
+      for _ in 0 ..< ItersMul:
+        when randZ:
+          let a = rng.random_unsafe_with_randZ(ECP_SWei_Proj[F])
+        else:
+          let a = rng.random_unsafe(ECP_SWei_Proj[F])
+
+        let exponent = rng.random_unsafe(BigInt[bits])
+        var exponentCanonical{.noInit.}: array[(bits+7) div 8, byte]
+        exponentCanonical.exportRawUint(exponent, bigEndian)
+
+        var
+          impl = a
+          reference = a
+          scratchSpace{.noInit.}: array[1 shl 4, ECP_SWei_Proj[F]]
+
+        impl.scalarMul(exponentCanonical, scratchSpace)
+        reference.unsafe_ECmul_double_add(exponentCanonical)
+
+        check: bool(impl == reference)
+
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = false)
+    test(Fp[BN254_Snarks], bits = BN254_Snarks_order_bits, randZ = true)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = false)
+    test(Fp[BLS12_381], bits = BLS12_381_order_bits, randZ = true)

tests/test_finite_fields_vs_gmp.nim

@@ -56,7 +56,7 @@ proc binary_prologue[C: static Curve, N: static int](
         a, b, p: var mpz_t,
         aTest, bTest: var Fp[C],
         aBuf, bBuf: var array[N, byte]) =
-  const bits = C.getCurveBitSize()
+  const bits = C.getCurveBitwidth()
 
   # Generate random value in the range 0 ..< 2^(bits-1)
   mpz_urandomb(a, gmpRng, uint bits)
@@ -123,7 +123,7 @@ proc addTests(gmpRng: var gmp_randstate_t, a, b, p, r: var mpz_t, C: static Curv
   # echo "Testing: random modular addition on ", $C
 
   const
-    bits = C.getCurveBitSize()
+    bits = C.getCurveBitwidth()
     bufLen = (bits + 7) div 8
   var
     aTest, bTest{.noInit.}: Fp[C]
@@ -146,7 +146,7 @@ proc subTests(gmpRng: var gmp_randstate_t, a, b, p, r: var mpz_t, C: static Curv
   # echo "Testing: random modular substraction on ", $C
 
   const
-    bits = C.getCurveBitSize()
+    bits = C.getCurveBitwidth()
     bufLen = (bits + 7) div 8
   var
     aTest, bTest{.noInit.}: Fp[C]
@@ -169,7 +169,7 @@ proc mulTests(gmpRng: var gmp_randstate_t, a, b, p, r: var mpz_t, C: static Curv
   # echo "Testing: random modular multiplication on ", $C
 
   const
-    bits = C.getCurveBitSize()
+    bits = C.getCurveBitwidth()
     bufLen = (bits + 7) div 8
   var
     aTest, bTest{.noInit.}: Fp[C]
@@ -189,7 +189,7 @@ proc invTests(gmpRng: var gmp_randstate_t, a, b, p, r: var mpz_t, C: static Curv
   # echo "Testing: random modular inversion on ", $C
 
   const
-    bits = C.getCurveBitSize()
+    bits = C.getCurveBitwidth()
     bufLen = (bits + 7) div 8
   var
     aTest, bTest{.noInit.}: Fp[C]