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,5267,5629
+msm_babyjubjub_2,bn254,plonk,12318,11795
+msm_bandersnatch_2,bls12_381,groth16,5532,6301
+msm_bandersnatch_2,bls12_381,plonk,13470,12918
+msm_jubjub_2,bls12_381,groth16,5273,5699
+msm_jubjub_2,bls12_381,plonk,12258,11799
+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,45565,45565
 pairing_bls12381,bn254,groth16,756837,1242260
@@ -95,18 +107,18 @@ pairing_bn254,bn254,groth16,506052,823961
 pairing_bn254,bn254,plonk,1646819,1573151
 pairing_bw6761,bn254,groth16,1589471,2646707
 pairing_bw6761,bn254,plonk,5318762,5097941
-scalar_mul_G1_bn254,bn254,groth16,115934,175413
-scalar_mul_G1_bn254,bn254,plonk,381171,365027
-scalar_mul_G1_bn254_incomplete,bn254,groth16,55441,87984
-scalar_mul_G1_bn254_incomplete,bn254,plonk,200004,192882
+scalar_mul_G1_bn254,bn254,groth16,108168,163915
+scalar_mul_G1_bn254,bn254,plonk,355353,340385
+scalar_mul_G1_bn254_incomplete,bn254,groth16,51579,81902
+scalar_mul_G1_bn254_incomplete,bn254,plonk,185916,179316
 scalar_mul_P256,bn254,groth16,96724,151768
 scalar_mul_P256,bn254,plonk,328895,315729
 scalar_mul_P256_incomplete,bn254,groth16,75542,121798
 scalar_mul_P256_incomplete,bn254,plonk,263160,253523
-scalar_mul_secp256k1,bn254,groth16,117264,177389
-scalar_mul_secp256k1,bn254,plonk,385623,369279
-scalar_mul_secp256k1_incomplete,bn254,groth16,56125,89066
-scalar_mul_secp256k1_incomplete,bn254,plonk,202518,195302
+scalar_mul_secp256k1,bn254,groth16,108204,163975
+scalar_mul_secp256k1,bn254,plonk,355502,340530
+scalar_mul_secp256k1_incomplete,bn254,groth16,51619,81970
+scalar_mul_secp256k1_incomplete,bn254,plonk,186082,179475
 selector/binaryMux_4,bn254,groth16,5,3
 selector/binaryMux_4,bls12_377,groth16,5,3
 selector/binaryMux_4,bls12_381,groth16,5,3

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)
+
 	registerSnippet("selector/mux_3", func(api frontend.API, newVariable func() frontend.Variable) {
 		selector.Mux(api, newVariable(), newVariable(), newVariable(), newVariable())
 	})

std/algebra/emulated/sw_emulated/hints.go

@@ -6,7 +6,7 @@ import (
 	"fmt"
 	"math/big"
 
-	"github.com/consensys/gnark-crypto/algebra/eisenstein"
+	"github.com/consensys/gnark-crypto/algebra/lattice"
 	"github.com/consensys/gnark-crypto/ecc"
 	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
 	bls12381_fp "github.com/consensys/gnark-crypto/ecc/bls12-381/fp"
@@ -32,8 +32,8 @@ func GetHints() []solver.Hint {
 	return []solver.Hint{
 		decomposeScalarG1,
 		scalarMulHint,
-		halfGCD,
-		halfGCDEisenstein,
+		rationalReconstruct,
+		rationalReconstructExt,
 	}
 }
 
@@ -160,7 +160,15 @@ func scalarMulHint(field *big.Int, inputs []*big.Int, outputs []*big.Int) error
 	})
 }
 
-func halfGCD(mod *big.Int, inputs []*big.Int, outputs []*big.Int) error {
+// rationalReconstruct decomposes a scalar s ∈ Fr into (s1, |s2|, signBit) such
+// that s1 ≡ s2·s (mod r), with |s1|, |s2| < γ₂·√r ≈ 1.15·√r (proven LLL/Hermite
+// bound from gnark-crypto/algebra/lattice). Replaces the older heuristic
+// HalfGCD-based decomposition.
+//
+// In-circuit: 1 native sign bit + 2 emulated outputs (s1, |s2|). The caller
+// reconstructs the signed s2 as ±|s2| based on the sign bit and asserts
+// s1 + s·s2 ≡ 0 (mod r).
+func rationalReconstruct(mod *big.Int, inputs []*big.Int, outputs []*big.Int) error {
 	return emulated.UnwrapHintContext(mod, inputs, outputs, func(hc emulated.HintContext) error {
 		moduli := hc.EmulatedModuli()
 		if len(moduli) != 1 {
@@ -177,25 +185,38 @@ func halfGCD(mod *big.Int, inputs []*big.Int, outputs []*big.Int) error {
 		if len(emuOutputs) != 2 {
 			return fmt.Errorf("expecting two outputs, got %d", len(emuOutputs))
 		}
-		glvBasis := new(ecc.Lattice)
-		ecc.PrecomputeLattice(moduli[0], emuInputs[0], glvBasis)
-		emuOutputs[0].Set(&glvBasis.V1[0])
-		emuOutputs[1].Set(&glvBasis.V1[1])
-		// we need the absolute values for the in-circuit computations,
-		// otherwise the negative values will be reduced modulo the SNARK scalar
-		// field and not the emulated field.
-		// 		output0 = |s0| mod r
-		// 		output1 = |s1| mod r
+		// lattice.RationalReconstruct returns (x, z) with x ≡ z·s (mod r),
+		// i.e., x − z·s ≡ 0 (mod r). The circuit expects: s1 + s·_s2 ≡ 0
+		// (mod r), so s1 = x and _s2 = −z.
+		rc := lattice.NewReconstructor(moduli[0])
+		res := rc.RationalReconstruct(emuInputs[0])
+		x, z := new(big.Int).Set(res[0]), new(big.Int).Set(res[1])
+
+		// Normalise so s1 ≥ 0; flipping (x, z) preserves x ≡ z·s mod r.
+		if x.Sign() < 0 {
+			x.Neg(x)
+			z.Neg(z)
+		}
+		emuOutputs[0].Set(x)
+		emuOutputs[1].Abs(z)
+
+		// signBit = 1 iff −z < 0 iff z > 0 (so the in-circuit code negates
+		// |z| to recover s2 = −z).
 		nativeOutputs[0].SetUint64(0)
-		if emuOutputs[1].Sign() == -1 {
-			emuOutputs[1].Neg(emuOutputs[1])
-			nativeOutputs[0].SetUint64(1) // we return the sign of the second subscalar
+		if z.Sign() > 0 {
+			nativeOutputs[0].SetUint64(1)
 		}
 		return nil
 	})
 }
 
-func halfGCDEisenstein(mod *big.Int, inputs []*big.Int, outputs []*big.Int) error {
+// rationalReconstructExt is the 4-D Eisenstein-style decomposition: given a
+// scalar s and GLV eigenvalue λ, finds (u1, u2, v1, v2) such that
+// s·(v1 + λ·v2) + u1 + λ·u2 ≡ 0 (mod r), with |u_i|, |v_i| < γ₄·r^(1/4) ≈
+// 1.25·r^(1/4) (proven LLL bound). Replaces the older Eisenstein HalfGCD.
+//
+// In-circuit: 4 native sign bits + 4 emulated absolute values.
+func rationalReconstructExt(mod *big.Int, inputs []*big.Int, outputs []*big.Int) error {
 	return emulated.UnwrapHintContext(mod, inputs, outputs, func(hc emulated.HintContext) error {
 		moduli := hc.EmulatedModuli()
 		if len(moduli) != 1 {
@@ -213,47 +234,43 @@ func halfGCDEisenstein(mod *big.Int, inputs []*big.Int, outputs []*big.Int) erro
 			return fmt.Errorf("expecting four outputs, got %d", len(emuOutputs))
 		}
 
-		glvBasis := new(ecc.Lattice)
-		ecc.PrecomputeLattice(moduli[0], emuInputs[1], glvBasis)
-		r := eisenstein.ComplexNumber{
-			A0: glvBasis.V1[0],
-			A1: glvBasis.V1[1],
-		}
-		sp := ecc.SplitScalar(emuInputs[0], glvBasis)
-		// in-circuit we check that Q - [s]P = 0 or equivalently Q + [-s]P = 0
-		// so here we return -s instead of s.
-		s := eisenstein.ComplexNumber{
-			A0: sp[0],
-			A1: sp[1],
-		}
-		s.Neg(&s)
+		// Inputs: emuInputs[0] = s, emuInputs[1] = λ.
+		// In-circuit we check Q − [s]P = 0, equivalently [−s]P + Q = 0, so we
+		// negate the scalar before reconstruction (matches the previous
+		// halfGCDEisenstein convention).
+		k := new(big.Int).Neg(emuInputs[0])
+		k.Mod(k, moduli[0])
+
+		rc := lattice.NewReconstructor(moduli[0]).SetLambda(emuInputs[1])
+		res := rc.RationalReconstructExt(k)
+		// res = (x, y, z, t) with k = (x + λ·y)/(z + λ·t) mod r,
+		// i.e., (x + λ·y) − k·(z + λ·t) ≡ 0 (mod r).
+		// Mapping onto our convention u1 + λ·u2 + s·(v1 + λ·v2) ≡ 0 with k = −s:
+		// u1 = x, u2 = y, v1 = z, v2 = t.
+		u1 := new(big.Int).Set(res[0])
+		u2 := new(big.Int).Set(res[1])
+		v1 := new(big.Int).Set(res[2])
+		v2 := new(big.Int).Set(res[3])
+
+		emuOutputs[0].Abs(u1)
+		emuOutputs[1].Abs(u2)
+		emuOutputs[2].Abs(v1)
+		emuOutputs[3].Abs(v2)
 
-		res := eisenstein.HalfGCD(&r, &s)
-		// values
-		emuOutputs[0].Set(&res[0].A0)
-		emuOutputs[1].Set(&res[0].A1)
-		emuOutputs[2].Set(&res[1].A0)
-		emuOutputs[3].Set(&res[1].A1)
-		// signs
 		nativeOutputs[0].SetUint64(0)
 		nativeOutputs[1].SetUint64(0)
 		nativeOutputs[2].SetUint64(0)
 		nativeOutputs[3].SetUint64(0)
-
-		if res[0].A0.Sign() == -1 {
-			emuOutputs[0].Neg(emuOutputs[0])
+		if u1.Sign() < 0 {
 			nativeOutputs[0].SetUint64(1)
 		}
-		if res[0].A1.Sign() == -1 {
-			emuOutputs[1].Neg(emuOutputs[1])
+		if u2.Sign() < 0 {
 			nativeOutputs[1].SetUint64(1)
 		}
-		if res[1].A0.Sign() == -1 {
-			emuOutputs[2].Neg(emuOutputs[2])
+		if v1.Sign() < 0 {
 			nativeOutputs[2].SetUint64(1)
 		}
-		if res[1].A1.Sign() == -1 {
-			emuOutputs[3].Neg(emuOutputs[3])
+		if v2.Sign() < 0 {
 			nativeOutputs[3].SetUint64(1)
 		}
 		return nil