wrong result with scalarMul in G2 curve

[bkomuves]

For a small number of possible exponents, scalarMul on G2 (tested with the BN254 curve) gives the wrong result.

Standalone Nim code to reproduce:

import constantine/math/arithmetic    
import constantine/math/io/io_fields  
import constantine/math/io/io_bigints
import constantine/math/config/curves  
import constantine/math/extension_fields/towers 
import constantine/math/elliptic/ec_shortweierstrass_affine    
import constantine/math/elliptic/ec_shortweierstrass_projective 
import constantine/math/elliptic/ec_scalar_mul                

#-------------------------------------------------------------------------------

type B  = BigInt[254]
type F  = Fp[BN254Snarks]
type F2 = QuadraticExt[F]
type G  = ECP_ShortW_Prj[F2, G2]
  
#-------------------------------------------------------------------------------

# size of the scalar field
let r : B = fromHex( B , "0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001" )

let expo : B = fromHex( B, "0x7b17fcc286b01af79176aa7da3a8615020eacda89a90e4ff5d0a085483f0448" )

let expoA : B = fromHex( B, "0x1234567890123456789001234567890" )
var expoB : B = expo
expoB -= expoA

let zeroF : F = fromHex( F, "0x00" )
let oneF  : F = fromHex( F, "0x01" )

#-------------------------------------------------------------------------------

# standard generator of G2

let gen2_xi : F = fromHex( F, "0x1adcd0ed10df9cb87040f46655e3808f98aa68a570acf5b0bde23fab1f149701" )
let gen2_xu : F = fromHex( F, "0x09e847e9f05a6082c3cd2a1d0a3a82e6fbfbe620f7f31269fa15d21c1c13b23b" )
let gen2_yi : F = fromHex( F, "0x056c01168a5319461f7ca7aa19d4fcfd1c7cdf52dbfc4cbee6f915250b7f6fc8" )
let gen2_yu : F = fromHex( F, "0x0efe500a2d02dd77f5f401329f30895df553b878fc3c0dadaaa86456a623235c" )

let gen2_x : F2 = F2( coords: [gen2_xi, gen2_xu] )
let gen2_y : F2 = F2( coords: [gen2_yi, gen2_yu] )
let gen2_z : F2 = F2( coords: [oneF   , zeroF  ] )

let gen2 : G = G( x: gen2_x, y: gen2_y, z: gen2_z )

#-------------------------------------------------------------------------------

proc printF( x: F ) =
  echo(" = " & x.toDecimal)
  
proc printF2( z: F2) =
  echo("   1 ~> " & z.coords[0].toDecimal )
  echo("   u ~> " & z.coords[1].toDecimal )


proc printG( pt: G ) =
  var aff : ECP_ShortW_Aff[F2, G2];
  aff.affine(pt)
  echo(" affine x coord: ");  printF2( aff.x )
  echo(" affine y coord: ");  printF2( aff.y )

#-------------------------------------------------------------------------------

var p : G
var q : G

echo("")
echo("sanity check: g2^r should be infinity")
p = gen2
p.scalarMul(r)
printG(p)

echo("")
echo("LHS = g2^expo")
p = gen2
p.scalarMul(expo)
printG(p)
let lhs : G = p

echo("")
echo("RHS = g2^expoA * g2^expoB, where expo = expoA + expoB")
p = gen2
q = gen2
p.scalarMul(expoA)
q.scalarMul(expoB)
p += q
printG(p)
let rhs : G = p

echo("")
echo("reference from SageMath")
echo("  sage x coord:")
echo("    1 -> 17216390949661727229956939928583223226083668728437958793715435751523027888005 ")
echo("    u -> 3082945034329785101034278215941854680789766318859358488904629243958221738137 ")
echo("  sage y coord:") 
echo("    1 -> 20108673238932196920264801702661201943173809015346082727725783869161803474440 ")
echo("    u -> 10405477402946058176045590740070709500904395284580129777629727895349459816649 ")

echo("")
echo("LHS - RHS = ")
p =  lhs
p -= rhs
printG(p)

#-------------------------------------------------------------------------------

#[

SageMath code

# BN128 elliptic curve
p  = 21888242871839275222246405745257275088696311157297823662689037894645226208583
r  = 21888242871839275222246405745257275088548364400416034343698204186575808495617
h  = 1
Fp = GF(p)
Fr = GF(r)
A  = Fp(0)
B  = Fp(3)
E  = EllipticCurve(Fp,[A,B])
gx = Fp(1)
gy = Fp(2)
gen = E(gx,gy)  # subgroup generator
print("scalar field check: ", gen.additive_order() == r )
print("cofactor check:     ", E.cardinality() == r*h )

# extension field
R.<x>   = Fp[]
Fp2.<u> = Fp.extension(x^2+1)

# twisted curve
B_twist = Fp2(19485874751759354771024239261021720505790618469301721065564631296452457478373 + 266929791119991161246907387137283842545076965332900288569378510910307636690*u )
E2 = EllipticCurve(Fp2,[0,B_twist])
size_E2     = E2.cardinality();
cofactor_E2 = size_E2 / r;

gen2_xi = Fp( 0x1adcd0ed10df9cb87040f46655e3808f98aa68a570acf5b0bde23fab1f149701 )
gen2_xu = Fp( 0x09e847e9f05a6082c3cd2a1d0a3a82e6fbfbe620f7f31269fa15d21c1c13b23b )
gen2_yi = Fp( 0x056c01168a5319461f7ca7aa19d4fcfd1c7cdf52dbfc4cbee6f915250b7f6fc8 )
gen2_yu = Fp( 0x0efe500a2d02dd77f5f401329f30895df553b878fc3c0dadaaa86456a623235c )

gen2_x = gen2_xi + u * gen2_xu
gen2_y = gen2_yi + u * gen2_yu

gen2 = E2(gen2_x, gen2_y)

print("g2^r: ", gen2*r )

expo = 0x7b17fcc286b01af79176aa7da3a8615020eacda89a90e4ff5d0a085483f0448

print("g2^expo: ")
print(gen2*expo)

]#

#-------------------------------------------------------------------------------

[bkomuves]

btw scalarMulGeneric gives the correct result, at least for this particular example (thanks @advaita-saha for noticing this!)

[mratsim]

Thanks for raising this,

After discussing with @advaita-saha apparently this is used for Groth16, you can use scalarMul_vartime for your use-case.

This is likely related to this comment:

TODO, min amount of bits for endomorphisms?

constantine/constantine/math/elliptic/ec_scalar_mul.nim

Lines 243 to 251 in c4be82d

	when BigInt.bits <= EC.F.C.getCurveOrderBitwidth() and
	EC.F.C.hasEndomorphismAcceleration():
	# TODO, min amount of bits for endomorphisms?
	when EC.F is Fp:
	P.scalarMulGLV_m2w2(scalar)
	elif EC.F is Fp2:
	P.scalarMulEndo(scalar)
	else: # Curves defined on Fp^m with m > 2
	{.error: "Unreachable".}

[mratsim]

Deeper Analysis

scalarMul on G2 gives the wrong result but:

• both scalarMulGeneric (no endomorphism)
• and scalarMul_vartime (endomorphism, same decomposition into k0, k1, k2, k3 scalars but different recoding instead of GLV-SAC)

gives the correct result.

The base GLV-SAC algorithm with only positive numbers uses length of ceildiv(bitlength(order), m)+1

But to deal with negative mini-scalars that may arise from lattice reduction:

constantine/constantine/math/constants/bn254_snarks_endomorphisms.nim

Lines 39 to 65 in dbd2630

	const BN254_Snarks_Lattice_G2* = (
	# (BigInt, isNeg)
	((BigInt[64].fromHex"0x89d3256894d213e2", false),
	(BigInt[63].fromHex"0x44e992b44a6909f2", false),
	(BigInt[63].fromHex"0x44e992b44a6909f1", true),
	(BigInt[63].fromHex"0x44e992b44a6909f1", false)),
	((BigInt[63].fromHex"0x44e992b44a6909f1", true),
	(BigInt[63].fromHex"0x44e992b44a6909f1", false),
	(BigInt[63].fromHex"0x44e992b44a6909f1", true),
	(BigInt[64].fromHex"0x89d3256894d213e3", true)),
	((BigInt[63].fromHex"0x44e992b44a6909f2", false),
	(BigInt[63].fromHex"0x44e992b44a6909f1", false),
	(BigInt[63].fromHex"0x44e992b44a6909f1", false),
	(BigInt[64].fromHex"0x89d3256894d213e2", true)),
	((BigInt[64].fromHex"0x89d3256894d213e3", false),
	(BigInt[63].fromHex"0x44e992b44a6909f1", true),
	(BigInt[63].fromHex"0x44e992b44a6909f2", true),
	(BigInt[63].fromHex"0x44e992b44a6909f1", true))
	)
	const BN254_Snarks_Babai_G2* = (
	# (BigInt, isNeg)
	(BigInt[192].fromHex"0x9e80318ab0d92b9308e5da66fc7184ae46f4bda995d51bb1", false),
	(BigInt[192].fromHex"0x9e80318ab0d92b9555b4ca7ba3e5577f2dff291532e42728", true),
	(BigInt[192].fromHex"0x9e80318ab0d92b9555b4ca7ba3e55782071c4c43fac4daff", false),
	(BigInt[192].fromHex"0x9e80318ab0d92b9555b4ca7ba3e5577dc170977dcef3cd3f", false)
	)

we need to:

• either posit the scalar and negate the matching curve point
• or use +2 instead of +1
  

Speculation

I assume the change to +2 guarantees an extra "sign" bit.

The numbers used are

order: 0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001
expo : 0x07b17fcc286b01af79176aa7da3a8615020eacda89a90e4ff5d0a085483f0448

The exponent is very small compared to the order. Hence I assume the lattice reduction which is supposed to reduce a 254-bit scalar into 4 64-bit scalars with substractions underflowed.

TODO

In that case, it may be possible that this avoids this very restrictive condition in the scalarMul_vartime to only use endomorphism on exact bitlength match:

constantine/constantine/math/elliptic/ec_scalar_mul_vartime.nim

Lines 353 to 377 in dbd2630

	const L = scalBits.ceilDiv_vartime(M) + 1
	let usedBits = scalar.limbs.getBits_LE_vartime()
	when scalBits == EC.F.C.getCurveOrderBitwidth() and
	EC.F.C.hasEndomorphismAcceleration():
	if usedBits >= L:
	when EC.F is Fp:
	P.scalarMulEndo_minHammingWeight_windowed_vartime(scalar, window = 4)
	elif EC.F is Fp2:
	P.scalarMulEndo_minHammingWeight_windowed_vartime(scalar, window = 3)
	else: # Curves defined on Fp^m with m > 2
	{.error: "Unreachable".}
	return
	if 64 < usedBits:
	# With a window of 5, we precompute 2^3 = 8 points
	P.scalarMul_minHammingWeight_windowed_vartime(scalar, window = 5)
	elif 16 < usedBits:
	# With a window of 3, we precompute 2^1 = 2 points
	P.scalarMul_minHammingWeight_windowed_vartime(scalar, window = 3)
	elif 4 < usedBits:
	P.scalarMul_doubleAdd_vartime(scalar)
	else:
	P.scalarMul_addchain_4bit_vartime(scalar)

Cost

The cost has been detailed in the code:

constantine/constantine/math/elliptic/ec_endomorphism_accel.nim

Lines 50 to 59 in dbd2630

	## A scalar decomposition might lead to negative miniscalar(s).
	## For proper handling it requires either:
	## 1. Negating it and then negating the corresponding curve point P
	## 2. Adding an extra bit to the recoding, which will do the right thing™
	##
	## For implementation solution 1 is faster:
	## - Double + Add is about 5000~8000 cycles on 6 64-bits limbs (BLS12-381)
	## - Conditional negate is about 10 cycles per Fp, on G2 projective we have 3 (coords) * 2 (Fp2) * 10 (cycles) ~= 60 cycles
	## We need to test the mini scalar, which is 65 bits so 2 Fp so about 2 cycles
	## and negate it as well.

An extra double+addition, on G2 there are 64 minimum add+dbl so cost is 1.56% perf

Alternative

We can decompose those scalars with an alternative method that doesn't produce negative mini-scalars, for example with euclidean division.