Add `dualScalarMul` implementation for $R = α A + β B$

[Vindaar]

In ECDSA we encounter operations of the form:

$R = α A + β B$

where the coefficients are scalars in the field $α, β ∈ 𝔽ᵣ$ and $A, B ∈ 𝔾$ elliptic curve points. These appear both in the verification and public key recovery.

Currently they are represented as:

  var
    point1 {.noinit.}: EC_ShortW_Jac[Fp[Name], G1]
    point2 {.noinit.}: EC_ShortW_Jac[Fp[Name], G1]
  point1.scalarMul(α, A)
  point2.scalarMul(β, B)
  var R {.noinit.}: EC_ShortW_Jac[Fp[Name], G1]
  R.sum(point1, point2)

but there is possible room for performance improvements. Refactoring these (and for potentially other future cases) is a good idea. Potential avenues for optimization:

constantine/constantine/math/elliptic/ec_scalar_mul_vartime.nim

Lines 252 to 331 in 6b65b0e

	func scalarMulEndo_wNAF_vartime*[scalBits: static int; EC](
	P: var EC,
	scalar: BigInt[scalBits],
	window: static int) {.tags:[VarTime, Alloca], meter.} =
	## Endomorphism-accelerated windowed vartime scalar multiplication
	##
	## P <- [k] P
	##
	## This uses windowed-NAF (wNAF)
	## This MUST NOT be used with secret data.
	##
	## This is highly VULNERABLE to timing attacks and power analysis attacks
	# Signed digits divides precomputation table size by 2
	# Odd-only divides precomputation table size by another 2
	const precompSize = 1 shl (window - 2)
	static: doAssert window < 8, "Window of size " & $window & " is too large and precomputation would use " & $(precompSize * sizeof(EC)) & " stack space."
	# 1. Compute endomorphisms
	const M = when P.F is Fp: 2
	elif P.F is Fp2: 4
	else: {.error: "Unconfigured".}
	const G = when EC isnot EC_ShortW_Aff|EC_ShortW_Jac|EC_ShortW_Prj: G1
	else: EC.G
	var endos {.noInit.}: array[M-1, EC]
	endos.computeEndomorphisms(P)
	# 2. Decompose scalar into mini-scalars
	const L = EC.getScalarField().bits().ceilDiv_vartime(M) + 1
	var miniScalars {.noInit.}: array[M, BigInt[L]]
	var negatePoints {.noInit.}: array[M, SecretBool]
	miniScalars.decomposeEndo(negatePoints, scalar, EC.getScalarField().bits(), EC.getName(), G)
	# 3. Handle negative mini-scalars
	if negatePoints[0].bool:
	P.neg()
	for m in 1 ..< M:
	if negatePoints[m].bool:
	endos[m-1].neg()
	# 4. EC precomputed table
	var tabEC {.noinit.}: array[M, array[precompSize, EC]]
	for m in 0 ..< M:
	var P2{.noInit.}: EC
	if m == 0:
	tabEC[0][0] = P
	P2.double(P)
	else:
	tabEC[m][0] = endos[m-1]
	P2.double(endos[m-1])
	for i in 1 ..< tabEC[m].len:
	tabEC[m][i].sum_vartime(tabEC[m][i-1], P2)
	var tab {.noinit.}: array[M, array[precompSize, affine(EC)]]
	tab.batchAffine(tabEC)
	# 5. wNAF precomputed tables
	const NafLen = L+1
	var tabNaf {.noinit.}: array[M, array[NafLen, int8]]
	for m in 0 ..< M:
	# tabNaf returns NAF from least-significant to most significant bits
	let miniScalarLen = tabNaf[m].recode_r2l_signed_window_vartime(miniScalars[m], window)
	# We compute from most significant to least significant
	# so we pad with 0
	for i in miniScalarLen ..< NafLen:
	tabNaf[m][i] = 0
	# 6. Compute
	var isInit = false
	for i in 0 ..< NafLen:
	if isInit:
	P.double()
	for m in 0 ..< M:
	if isInit:
	P.accumNAF(tab[m], tabNaf[m], NafLen, i)
	else:
	isInit = P.initNAF(tab[m], tabNaf[m], NafLen, i)

The scalarMulEndo implementation handles such cases internally efficiently. We may be able to extract those internals.

The two locations where this code currently appears:

• constantine/constantine/signatures/ecdsa.nim

  Lines 283 to 292 in 9642ca6

  	# 4. Compute u₁G + u₂Q
  	var
  	point1 {.noinit.}: EC_ShortW_Jac[Fp[Name], G1]
  	point2 {.noinit.}: EC_ShortW_Jac[Fp[Name], G1]
  	# Generator of the curve
  	const G = publicKey.F.Name.getGenerator($publicKey.G)
  	point1.scalarMul(u1, G)
  	point2.scalarMul(u2, publicKey)
  	var R {.noinit.}: EC_ShortW_Jac[Fp[Name], G1]
  	R.sum(point1, point2)

• constantine/constantine/signatures/ecdsa.nim

  Lines 371 to 375 in 664d985

  	var Q {.noinit.}: ECJac # the potential public key
  	var point1 {.noinit.}, point2 {.noinit.}: ECJac
  	point1.scalarMul(u1, G) # `p₁ = u₁ * G`
  	point2.scalarMul(u2, R) # `p₂ = u₂ * R`
  	Q.sum(point1, point2) # `Q = p₁ + p₂`