perf: twisted Edwards lattice-based double-base MSM

internal/stats/latest_stats.csv

@@ -87,6 +87,18 @@ math/emulated/secp256k1_64,bn254,plonk,4025,3923
 math/emulated/secp256k1_64,bls12_377,plonk,4025,3923
 math/emulated/secp256k1_64,bls12_381,plonk,4025,3923
 math/emulated/secp256k1_64,bw6_761,plonk,4025,3923
+msm_G1_bn254_2,bn254,groth16,208925,312617
+msm_G1_bn254_2,bn254,plonk,688811,658743
+msm_P256_2,bn254,groth16,185846,288056
+msm_P256_2,bn254,plonk,635297,608874
+msm_babyjubjub_2,bn254,groth16,5269,5683
+msm_babyjubjub_2,bn254,plonk,12389,11848
+msm_bandersnatch_2,bls12_381,groth16,5016,5791
+msm_bandersnatch_2,bls12_381,plonk,12450,11904
+msm_jubjub_2,bls12_381,groth16,5276,5754
+msm_jubjub_2,bls12_381,plonk,12332,11855
+msm_secp256k1_2,bn254,groth16,208997,312737
+msm_secp256k1_2,bn254,plonk,689104,659028
 pairing_bls12377,bw6_761,groth16,11876,11876
 pairing_bls12377,bw6_761,plonk,41914,40431
 pairing_bls12381,bn254,groth16,756837,1242260

internal/stats/snippet.go

@@ -6,13 +6,15 @@ import (
 
 	"github.com/consensys/gnark"
 	"github.com/consensys/gnark-crypto/ecc"
+	twistededwardsCryptoID "github.com/consensys/gnark-crypto/ecc/twistededwards"
 	"github.com/consensys/gnark/frontend"
 	"github.com/consensys/gnark/std/algebra/algopts"
 	"github.com/consensys/gnark/std/algebra/emulated/sw_bls12381"
 	"github.com/consensys/gnark/std/algebra/emulated/sw_bn254"
 	"github.com/consensys/gnark/std/algebra/emulated/sw_bw6761"
 	"github.com/consensys/gnark/std/algebra/emulated/sw_emulated"
 	"github.com/consensys/gnark/std/algebra/native/sw_bls12377"
+	"github.com/consensys/gnark/std/algebra/native/twistededwards"
 	"github.com/consensys/gnark/std/hash/mimc"
 	"github.com/consensys/gnark/std/math/bits"
 	"github.com/consensys/gnark/std/math/emulated"
@@ -412,6 +414,144 @@ func initSnippets() {
 
 	}, ecc.BN254)
 
+	// MSM(2, n) snippets for the four curve classes — used to evaluate which
+	// MSM-size variant is best in complete-arithmetic mode (Phase 4 of plan).
+	// Baselines: existing scalar_mul_* divided by 2 gives lower bound for
+	// MSM(2, n) via two ScalarMul + Add.
+	registerSnippet("msm_secp256k1_2", func(api frontend.API, newVariable func() frontend.Variable) {
+		cr, err := sw_emulated.New[emulated.Secp256k1Fp, emulated.Secp256k1Fr](api, sw_emulated.GetCurveParams[emulated.Secp256k1Fp]())
+		if err != nil {
+			panic(err)
+		}
+		fr, _ := emulated.NewField[emulated.Secp256k1Fr](api)
+		newFr := func() *emulated.Element[emulated.Secp256k1Fr] {
+			n, _ := emulated.GetEffectiveFieldParams[emulated.Secp256k1Fr](api.Compiler().Field())
+			limbs := make([]frontend.Variable, n)
+			for i := range limbs {
+				limbs[i] = newVariable()
+			}
+			return fr.NewElement(limbs)
+		}
+		fp, _ := emulated.NewField[emulated.Secp256k1Fp](api)
+		newPoint := func() *sw_emulated.AffinePoint[emulated.Secp256k1Fp] {
+			n, _ := emulated.GetEffectiveFieldParams[emulated.Secp256k1Fp](api.Compiler().Field())
+			x := make([]frontend.Variable, n)
+			y := make([]frontend.Variable, n)
+			for i := range x {
+				x[i] = newVariable()
+				y[i] = newVariable()
+			}
+			return &sw_emulated.AffinePoint[emulated.Secp256k1Fp]{X: *fp.NewElement(x), Y: *fp.NewElement(y)}
+		}
+		_, _ = cr.MultiScalarMul(
+			[]*sw_emulated.AffinePoint[emulated.Secp256k1Fp]{newPoint(), newPoint()},
+			[]*emulated.Element[emulated.Secp256k1Fr]{newFr(), newFr()},
+		)
+	}, ecc.BN254)
+
+	registerSnippet("msm_P256_2", func(api frontend.API, newVariable func() frontend.Variable) {
+		cr, err := sw_emulated.New[emulated.P256Fp, emulated.P256Fr](api, sw_emulated.GetCurveParams[emulated.P256Fp]())
+		if err != nil {
+			panic(err)
+		}
+		fr, _ := emulated.NewField[emulated.P256Fr](api)
+		newFr := func() *emulated.Element[emulated.P256Fr] {
+			n, _ := emulated.GetEffectiveFieldParams[emulated.P256Fr](api.Compiler().Field())
+			limbs := make([]frontend.Variable, n)
+			for i := range limbs {
+				limbs[i] = newVariable()
+			}
+			return fr.NewElement(limbs)
+		}
+		fp, _ := emulated.NewField[emulated.P256Fp](api)
+		newPoint := func() *sw_emulated.AffinePoint[emulated.P256Fp] {
+			n, _ := emulated.GetEffectiveFieldParams[emulated.P256Fp](api.Compiler().Field())
+			x := make([]frontend.Variable, n)
+			y := make([]frontend.Variable, n)
+			for i := range x {
+				x[i] = newVariable()
+				y[i] = newVariable()
+			}
+			return &sw_emulated.AffinePoint[emulated.P256Fp]{X: *fp.NewElement(x), Y: *fp.NewElement(y)}
+		}
+		_, _ = cr.MultiScalarMul(
+			[]*sw_emulated.AffinePoint[emulated.P256Fp]{newPoint(), newPoint()},
+			[]*emulated.Element[emulated.P256Fr]{newFr(), newFr()},
+		)
+	}, ecc.BN254)
+
+	registerSnippet("msm_G1_bn254_2", func(api frontend.API, newVariable func() frontend.Variable) {
+		cr, err := sw_emulated.New[emulated.BN254Fp, emulated.BN254Fr](api, sw_emulated.GetCurveParams[emulated.BN254Fp]())
+		if err != nil {
+			panic(err)
+		}
+		fr, _ := emulated.NewField[emulated.BN254Fr](api)
+		newFr := func() *emulated.Element[emulated.BN254Fr] {
+			n, _ := emulated.GetEffectiveFieldParams[emulated.BN254Fr](api.Compiler().Field())
+			limbs := make([]frontend.Variable, n)
+			for i := range limbs {
+				limbs[i] = newVariable()
+			}
+			return fr.NewElement(limbs)
+		}
+		fp, _ := emulated.NewField[emulated.BN254Fp](api)
+		newPoint := func() *sw_emulated.AffinePoint[emulated.BN254Fp] {
+			n, _ := emulated.GetEffectiveFieldParams[emulated.BN254Fp](api.Compiler().Field())
+			x := make([]frontend.Variable, n)
+			y := make([]frontend.Variable, n)
+			for i := range x {
+				x[i] = newVariable()
+				y[i] = newVariable()
+			}
+			return &sw_emulated.AffinePoint[emulated.BN254Fp]{X: *fp.NewElement(x), Y: *fp.NewElement(y)}
+		}
+		_, _ = cr.MultiScalarMul(
+			[]*sw_emulated.AffinePoint[emulated.BN254Fp]{newPoint(), newPoint()},
+			[]*emulated.Element[emulated.BN254Fr]{newFr(), newFr()},
+		)
+	}, ecc.BN254)
+
+	// Twisted Edwards DoubleBaseScalarMul snippets — exercise the new
+	// MSM(3, 2n/3) (no GLV) and MSM(6, n/3) (GLV) variants from PR #1697.
+	registerSnippet("msm_babyjubjub_2", func(api frontend.API, newVariable func() frontend.Variable) {
+		curve, err := twistededwards.NewEdCurve(api, twistededwardsCryptoID.BN254)
+		if err != nil {
+			panic(err)
+		}
+		var P1, P2 twistededwards.Point
+		P1.X = newVariable()
+		P1.Y = newVariable()
+		P2.X = newVariable()
+		P2.Y = newVariable()
+		_ = curve.DoubleBaseScalarMulNonZero(P1, P2, newVariable(), newVariable())
+	}, ecc.BN254)
+
+	registerSnippet("msm_jubjub_2", func(api frontend.API, newVariable func() frontend.Variable) {
+		curve, err := twistededwards.NewEdCurve(api, twistededwardsCryptoID.BLS12_381)
+		if err != nil {
+			panic(err)
+		}
+		var P1, P2 twistededwards.Point
+		P1.X = newVariable()
+		P1.Y = newVariable()
+		P2.X = newVariable()
+		P2.Y = newVariable()
+		_ = curve.DoubleBaseScalarMulNonZero(P1, P2, newVariable(), newVariable())
+	}, ecc.BLS12_381)
+
+	registerSnippet("msm_bandersnatch_2", func(api frontend.API, newVariable func() frontend.Variable) {
+		curve, err := twistededwards.NewEdCurve(api, twistededwardsCryptoID.BLS12_381_BANDERSNATCH)
+		if err != nil {
+			panic(err)
+		}
+		var P1, P2 twistededwards.Point
+		P1.X = newVariable()
+		P1.Y = newVariable()
+		P2.X = newVariable()
+		P2.Y = newVariable()
+		_ = curve.DoubleBaseScalarMulNonZero(P1, P2, newVariable(), newVariable())
+	}, ecc.BLS12_381)
+
 	// G2 scalar mul snippets — exercise the GLV+FakeGLV path to track its cost.
 	registerSnippet("scalar_mul_G2_bls12381", func(api frontend.API, newVariable func() frontend.Variable) {
 		bls12381fp, _ := emulated.NewField[emulated.BLS12381Fp](api)

std/algebra/native/twistededwards/curve.go

@@ -10,6 +10,7 @@ type curve struct {
 	api    frontend.API
 	id     twistededwards.ID
 	params *CurveParams
+	endo   *EndoParams // non-nil iff the curve has a GLV endomorphism (Bandersnatch)
 }
 
 func (c *curve) Params() *CurveParams {
@@ -44,8 +45,29 @@ func (c *curve) ScalarMul(p1 Point, scalar frontend.Variable) Point {
 	p.scalarMul(c.api, &p1, scalar, c.params)
 	return p
 }
+
+// DoubleBaseScalarMul computes s1*p1 + s2*p2. It is complete for all scalar
+// inputs, including zero, and for identity points.
 func (c *curve) DoubleBaseScalarMul(p1, p2 Point, s1, s2 frontend.Variable) Point {
 	var p Point
 	p.doubleBaseScalarMul(c.api, &p1, &p2, s1, s2, c.params)
 	return p
 }
+
+// DoubleBaseScalarMulNonZero computes s1*p1 + s2*p2 using the most efficient
+// lattice-based MSM variant available for the curve:
+//   - GLV-equipped curves (Bandersnatch): 6-MSM with r^(1/3)-bounded sub-scalars.
+//   - non-GLV curves (Jubjub, BabyJubjub, edBLS12-377, edBW6-761): 3-MSM with
+//     r^(2/3)-bounded sub-scalars and LogUp lookups.
+//
+// The scalars s1, s2 must be nonzero and p1, p2 must not be the TE identity
+// (0, 1). Use DoubleBaseScalarMul for complete edge-case handling.
+func (c *curve) DoubleBaseScalarMulNonZero(p1, p2 Point, s1, s2 frontend.Variable) Point {
+	var p Point
+	if c.endo != nil {
+		p.doubleBaseScalarMul6MSMLogUp(c.api, &p1, &p2, s1, s2, c.params, c.endo)
+	} else {
+		p.doubleBaseScalarMul3MSMLogUp(c.api, &p1, &p2, s1, s2, c.params)
+	}
+	return p
+}

std/algebra/native/twistededwards/curve_test.go

@@ -191,6 +191,22 @@ func (circuit *doubleBaseScalarMulCircuit) Define(api frontend.API) error {
 	return nil
 }
 
+type doubleBaseScalarMulNonZeroCircuit struct {
+	curveID twistededwards.ID
+	P1, P2  Point
+	S1, S2  frontend.Variable
+	Result  Point
+}
+
+func (circuit *doubleBaseScalarMulNonZeroCircuit) Define(api frontend.API) error {
+	curve, err := NewEdCurve(api, circuit.curveID)
+	if err != nil {
+		return err
+	}
+	assertPointEqual(api, curve.DoubleBaseScalarMulNonZero(circuit.P1, circuit.P2, circuit.S1, circuit.S2), circuit.Result)
+	return nil
+}
+
 func TestAdd(t *testing.T) {
 	for _, curveID := range curves {
 		params, err := GetCurveParams(curveID)
@@ -296,6 +312,21 @@ func TestDoubleBaseScalarMul(t *testing.T) {
 	}
 }
 
+func TestDoubleBaseScalarMulNonZero(t *testing.T) {
+	for _, curveID := range curves {
+		params, err := GetCurveParams(curveID)
+		if err != nil {
+			t.Fatalf("%s: get curve params: %v", curveLabel(curveID), err)
+		}
+		data := randomTestData(params, curveID)
+		circuit := &doubleBaseScalarMulNonZeroCircuit{curveID: curveID}
+		witness := &doubleBaseScalarMulNonZeroCircuit{P1: data.P1, P2: data.P2, S1: data.S1, S2: data.S2, Result: data.DoubleScalarMulResult}
+		invalidWitness := *witness
+		invalidWitness.Result = offCurvePoint()
+		checkCircuitForCurve(t, curveID, circuit, witness, &invalidWitness)
+	}
+}
+
 func TestAddEdgeCases(t *testing.T) {
 	for _, curveID := range curves {
 		params, err := GetCurveParams(curveID)
@@ -380,18 +411,24 @@ func TestFixedScalarMulEdgeCases(t *testing.T) {
 	}
 }
 
+// TestDoubleBaseScalarMulEdgeCases covers the complete public method, including
+// zero scalars and identity points. The optimized NonZero variant is tested
+// separately.
 func TestDoubleBaseScalarMulEdgeCases(t *testing.T) {
 	for _, curveID := range curves {
 		params, err := GetCurveParams(curveID)
 		if err != nil {
 			t.Fatalf("%s: get curve params: %v", curveLabel(curveID), err)
 		}
 		data := testDataForScalars(params, curveID, big.NewInt(1), big.NewInt(2))
+		base := Point{X: params.Base[0], Y: params.Base[1]}
 
 		t.Run(curveLabel(curveID), func(t *testing.T) {
-			assertSolvedForCurve(t, curveID, &doubleBaseScalarMulCircuit{curveID: curveID}, &doubleBaseScalarMulCircuit{P1: data.P1, P2: data.P2, S1: 0, S2: 0, Result: identityPoint()})
-			assertSolvedForCurve(t, curveID, &doubleBaseScalarMulCircuit{curveID: curveID}, &doubleBaseScalarMulCircuit{P1: data.P1, P2: data.P2, S1: 1, S2: 0, Result: data.P1})
-			assertSolvedForCurve(t, curveID, &doubleBaseScalarMulCircuit{curveID: curveID}, &doubleBaseScalarMulCircuit{P1: data.P1, P2: data.P2, S1: 0, S2: 1, Result: data.P2})
+			circuit := &doubleBaseScalarMulCircuit{curveID: curveID}
+			assertSolvedForCurve(t, curveID, circuit, &doubleBaseScalarMulCircuit{P1: data.P1, P2: data.P2, S1: 0, S2: 0, Result: identityPoint()})
+			assertSolvedForCurve(t, curveID, circuit, &doubleBaseScalarMulCircuit{P1: data.P1, P2: data.P2, S1: 1, S2: 0, Result: data.P1})
+			assertSolvedForCurve(t, curveID, circuit, &doubleBaseScalarMulCircuit{P1: data.P1, P2: data.P2, S1: 0, S2: 1, Result: data.P2})
+			assertSolvedForCurve(t, curveID, circuit, &doubleBaseScalarMulCircuit{P1: identityPoint(), P2: base, S1: 1, S2: 2, Result: data.P2})
 		})
 	}
 }
@@ -630,8 +667,8 @@ func zeroRationalReconstructHint(_ *big.Int, inputs, outputs []*big.Int) error {
 	if len(inputs) != 2 {
 		return errors.New("expecting two inputs")
 	}
-	if len(outputs) != 4 {
-		return errors.New("expecting four outputs")
+	if len(outputs) != 3 {
+		return errors.New("expecting three outputs")
 	}
 	for i := range outputs {
 		outputs[i].SetUint64(0)
@@ -659,3 +696,96 @@ func TestScalarMulFakeGLVRegressionTrivialDecomposition(t *testing.T) {
 	)
 	assert.Error(err)
 }
+
+func forgedBN254DoubleBaseScalarMulHint(_ *big.Int, inputs, outputs []*big.Int) error {
+	if len(inputs) != 7 {
+		return errors.New("expecting seven inputs")
+	}
+	if len(outputs) != 4 {
+		return errors.New("expecting four outputs")
+	}
+	var p1, p2, q1, q2 tbn254.PointAffine
+	p1.X.SetBigInt(inputs[0])
+	p1.Y.SetBigInt(inputs[1])
+	p2.X.SetBigInt(inputs[3])
+	p2.Y.SetBigInt(inputs[4])
+	q1.ScalarMultiplication(&p1, inputs[2])
+	q2.ScalarMultiplication(&p2, inputs[5])
+
+	var delta tbn254.PointAffine
+	delta.Set(&p1)
+
+	var q1Hint tbn254.PointAffine
+	q1Hint.Add(&q1, &delta)
+
+	q1Hint.X.BigInt(outputs[0])
+	q1Hint.Y.BigInt(outputs[1])
+	q2.X.BigInt(outputs[2])
+	q2.Y.BigInt(outputs[3])
+	return nil
+}
+
+func forgedBN254DoubleBaseResult(params *CurveParams, s1, s2 *big.Int) (Point, error) {
+	var p1, p2 tbn254.PointAffine
+	p1.X.SetBigInt(params.Base[0])
+	p1.Y.SetBigInt(params.Base[1])
+	p2.Set(&p1)
+	p1.ScalarMultiplication(&p1, s1)
+	p2.ScalarMultiplication(&p2, s2)
+
+	p1X, p1Y := new(big.Int), new(big.Int)
+	p2X, p2Y := new(big.Int), new(big.Int)
+	p1.X.BigInt(p1X)
+	p1.Y.BigInt(p1Y)
+	p2.X.BigInt(p2X)
+	p2.Y.BigInt(p2Y)
+
+	inputs := []*big.Int{
+		p1X,
+		p1Y,
+		new(big.Int).Set(s1),
+		p2X,
+		p2Y,
+		new(big.Int).Set(s2),
+		new(big.Int).Set(params.Order),
+	}
+	outputs := []*big.Int{new(big.Int), new(big.Int), new(big.Int), new(big.Int)}
+	if err := forgedBN254DoubleBaseScalarMulHint(nil, inputs, outputs); err != nil {
+		return Point{}, err
+	}
+	var q1, q2, r tbn254.PointAffine
+	q1.X.SetBigInt(outputs[0])
+	q1.Y.SetBigInt(outputs[1])
+	q2.X.SetBigInt(outputs[2])
+	q2.Y.SetBigInt(outputs[3])
+	r.Add(&q1, &q2)
+	rX, rY := new(big.Int), new(big.Int)
+	r.X.BigInt(rX)
+	r.Y.BigInt(rY)
+	return Point{X: rX, Y: rY}, nil
+}
+
+func TestDoubleBaseScalarMulNonZeroRejectsForgedPartialHints(t *testing.T) {
+	assert := require.New(t)
+	params, err := GetCurveParams(twistededwards.BN254)
+	assert.NoError(err)
+
+	data := testDataForScalars(params, twistededwards.BN254, big.NewInt(5), big.NewInt(7))
+	forged, err := forgedBN254DoubleBaseResult(params, data.S1, data.S2)
+	assert.NoError(err)
+
+	witness := doubleBaseScalarMulNonZeroCircuit{
+		P1:     data.P1,
+		P2:     data.P2,
+		S1:     data.S1,
+		S2:     data.S2,
+		Result: forged,
+	}
+	err = test.IsSolved(
+		&doubleBaseScalarMulNonZeroCircuit{curveID: twistededwards.BN254},
+		&witness,
+		ecc.BN254.ScalarField(),
+		test.WithReplacementHint(solver.GetHintID(doubleBaseScalarMulHint), forgedBN254DoubleBaseScalarMulHint),
+	)
+	assert.Error(err)
+}