From 1df8f25445a0e38cf9a01333da3c3788ffd27266 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Thu, 23 Jul 2026 20:38:06 -0400 Subject: [PATCH 01/15] perf: comb fixed-base scalar mul --- std/algebra/emulated/sw_emulated/fixedbase.go | 486 ++++++++++++++++++ .../emulated/sw_emulated/fixedbase_test.go | 187 +++++++ std/algebra/emulated/sw_emulated/hints.go | 2 + std/algebra/emulated/sw_emulated/point.go | 3 + std/math/emulated/field.go | 15 + std/math/emulated/field_mul.go | 168 ++++-- 6 files changed, 816 insertions(+), 45 deletions(-) create mode 100644 std/algebra/emulated/sw_emulated/fixedbase.go create mode 100644 std/algebra/emulated/sw_emulated/fixedbase_test.go diff --git a/std/algebra/emulated/sw_emulated/fixedbase.go b/std/algebra/emulated/sw_emulated/fixedbase.go new file mode 100644 index 0000000000..8eaff714bf --- /dev/null +++ b/std/algebra/emulated/sw_emulated/fixedbase.go @@ -0,0 +1,486 @@ +package sw_emulated + +import ( + "errors" + "fmt" + "math/big" + + "github.com/consensys/gnark/frontend" + limbs "github.com/consensys/gnark/std/internal/limbcomposition" + "github.com/consensys/gnark/std/math/emulated" +) + +// This file implements a fixed-base scalar multiplication using a signed-digit +// comb method with compile-time constant window tables. +// +// The scalar s is recoded into an odd integer k' = s + 1 − b0 (b0 the parity +// of s) represented in signed binary digits d_i ∈ {−1, +1}: +// +// k' = Σ_{i 14 { + return nil, fmt.Errorf("unsupported window width %d", w) + } + rBits := r.BitLen() + nw := (rBits + w - 1) / w + if nw < 2 { + return nil, fmt.Errorf("window width %d too large for %d-bit scalar field", w, rBits) + } + n := w * nw + if n > scalarCap { + return nil, fmt.Errorf("recoded scalar needs %d bits, scalar field emulation has capacity %d", n, scalarCap) + } + // adding window t (t ≥ 1) with incomplete formulas is safe iff + // 2^{w·(t+1)} ≤ r; the final addition is always complete as it may cancel + // to the point at infinity or double. + nbUnified := 0 + for t := nw - 1; t >= 1; t-- { + if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + break + } + nbUnified++ + } + nbUnified = max(nbUnified, 1) + + G := &combAffine{x: new(big.Int).Set(gx), y: new(big.Int).Set(gy)} + half := 1 << (w - 1) + windows := make([][][2]*big.Int, nw) + Bt := G + var err error + for t := 0; t < nw; t++ { + if t > 0 { + for k := 0; k < w; k++ { + if Bt, err = combDouble(Bt, a, prime); err != nil { + return nil, err + } + } + } + // odd multiples odd[m] = [(2m+1)·2^{w·t}]G + D, err := combDouble(Bt, a, prime) + if err != nil { + return nil, err + } + odd := make([]*combAffine, half) + odd[0] = Bt + for m := 1; m < half; m++ { + if odd[m], err = combAdd(odd[m-1], D, prime); err != nil { + return nil, err + } + } + tab := make([][2]*big.Int, 1< 0 { + pt = odd[(d-1)/2] + } else { + pt = combNeg(odd[(-d-1)/2], prime) + } + tab[j] = [2]*big.Int{pt.x, pt.y} + } + windows[t] = tab + } + // parity-folded top window: topEven[j] = windows[nw−1][j] + (−G) + negG := combNeg(G, prime) + topEven := make([][2]*big.Int, 1<= 1 { + yAcc := c.baseApi.Eval([][]*emulated.Element[B]{{lamPrev, c.baseApi.Sub(xTPrev, xAcc)}, {yTPrev}}, []int{1, -1}) + acc = &AffinePoint[B]{X: *xAcc, Y: *yAcc} + } else { + acc = &AffinePoint[B]{X: *xT[0], Y: *yT[0]} + } + for t := nbInc + 1; t <= nw-1; t++ { + acc = c.AddUnified(acc, &AffinePoint[B]{X: *xT[t], Y: *yT[t]}) + } + return acc +} + +// combRecodeHint computes the comb recoding of the scalar: given the scalar s +// (nonnative), it returns as native outputs the parity bit b0 = s mod 2 +// followed by the n bits of c = (k' + 2^n − 1)/2 where k' = s + 1 − b0 and n +// is the number of signed digits (inferred from the output count). +func combRecodeHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHintWithNativeOutput(inputs, outputs, func(r *big.Int, in, out []*big.Int) error { + if len(in) != 1 { + return errors.New("expecting one input") + } + if len(out) < 2 { + return errors.New("expecting at least two outputs") + } + n := len(out) - 1 + s := new(big.Int).Mod(in[0], r) + b0 := s.Bit(0) + kp := new(big.Int).Set(s) + if b0 == 0 { + kp.Add(kp, big.NewInt(1)) + } + // c = (k' + 2^n − 1)/2 + cv := new(big.Int).Lsh(big.NewInt(1), uint(n)) + cv.Sub(cv, big.NewInt(1)).Add(cv, kp).Rsh(cv, 1) + out[0].SetUint64(uint64(b0)) + for i := 0; i < n; i++ { + out[1+i].SetUint64(uint64(cv.Bit(i))) + } + return nil + }) +} + +// combChainLambdaHint computes the chord slope of the next comb chain +// addition. Inputs (nonnative, base field): λprev, x (the accumulator +// x-coordinate), xTprev, yTprev (the previously added table point), xT, yT +// (the table point being added). The accumulator y-coordinate is recomputed +// in its implicit form y = λprev·(xTprev − x) − yTprev and the output is +// λ = (yT − y) / (xT − x). +func combChainLambdaHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 6 { + return errors.New("expecting six inputs") + } + if len(out) != 1 { + return errors.New("expecting one output") + } + lamPrev, x, xTPrev, yTPrev, xT, yT := in[0], in[1], in[2], in[3], in[4], in[5] + y := new(big.Int).Sub(xTPrev, x) + y.Mul(y, lamPrev).Sub(y, yTPrev).Mod(y, p) + den := new(big.Int).Sub(xT, x) + den.Mod(den, p) + if den.Sign() == 0 { + return errors.New("comb chain: x-coordinate collision") + } + den.ModInverse(den, p) + out[0].Sub(yT, y).Mul(out[0], den).Mod(out[0], p) + return nil + }) +} diff --git a/std/algebra/emulated/sw_emulated/fixedbase_test.go b/std/algebra/emulated/sw_emulated/fixedbase_test.go new file mode 100644 index 0000000000..a1dd0b1418 --- /dev/null +++ b/std/algebra/emulated/sw_emulated/fixedbase_test.go @@ -0,0 +1,187 @@ +package sw_emulated + +import ( + "math/big" + "testing" + + "github.com/consensys/gnark-crypto/ecc/bn254" + fr_bn "github.com/consensys/gnark-crypto/ecc/bn254/fr" + "github.com/consensys/gnark-crypto/ecc/secp256k1" + fr_secp "github.com/consensys/gnark-crypto/ecc/secp256k1/fr" + "github.com/consensys/gnark/frontend" + "github.com/consensys/gnark/frontend/cs/r1cs" + "github.com/consensys/gnark/frontend/cs/scs" + "github.com/consensys/gnark/std/math/emulated" + "github.com/consensys/gnark/test" +) + +type ScalarMulBaseCombTest[T, S emulated.FieldParams] struct { + Q AffinePoint[T] + S emulated.Element[S] + w int +} + +func (c *ScalarMulBaseCombTest[T, S]) Define(api frontend.API) error { + cr, err := New[T, S](api, GetCurveParams[T]()) + if err != nil { + return err + } + res := cr.scalarMulBaseComb(&c.S, c.w) + cr.AssertIsEqual(res, &c.Q) + return nil +} + +// edge scalars: 0, 1, 2, r−2, r−1 and a 2-power, plus random ones. The point +// [0]G is represented as (0,0) following the package convention. +func combTestScalars(r *big.Int, nbRandom int, randFn func() *big.Int) []*big.Int { + scalars := []*big.Int{ + big.NewInt(0), + big.NewInt(1), + big.NewInt(2), + big.NewInt(3), + new(big.Int).Sub(r, big.NewInt(1)), + new(big.Int).Sub(r, big.NewInt(2)), + new(big.Int).Lsh(big.NewInt(1), 128), + } + for i := 0; i < nbRandom; i++ { + scalars = append(scalars, randFn()) + } + return scalars +} + +func TestScalarMulBaseCombSecp256k1(t *testing.T) { + assert := test.NewAssert(t) + _, g := secp256k1.Generators() + r := fr_secp.Modulus() + randFn := func() *big.Int { + var rnd fr_secp.Element + _, _ = rnd.SetRandom() + return rnd.BigInt(new(big.Int)) + } + for _, w := range []int{4, 8} { + for _, s := range combTestScalars(r, 3, randFn) { + var S secp256k1.G1Affine + S.ScalarMultiplication(&g, s) + circuit := ScalarMulBaseCombTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{w: w} + witness := ScalarMulBaseCombTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{ + S: emulated.ValueOf[emulated.Secp256k1Fr](s), + Q: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](S.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](S.Y), + }, + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err, "w=%d s=%s", w, s.String()) + } + } +} + +func TestScalarMulBaseCombBN254(t *testing.T) { + assert := test.NewAssert(t) + _, _, g, _ := bn254.Generators() + r := fr_bn.Modulus() + randFn := func() *big.Int { + var rnd fr_bn.Element + _, _ = rnd.SetRandom() + return rnd.BigInt(new(big.Int)) + } + for _, w := range []int{4, 8} { + for _, s := range combTestScalars(r, 3, randFn) { + var S bn254.G1Affine + S.ScalarMultiplication(&g, s) + circuit := ScalarMulBaseCombTest[emulated.BN254Fp, emulated.BN254Fr]{w: w} + witness := ScalarMulBaseCombTest[emulated.BN254Fp, emulated.BN254Fr]{ + S: emulated.ValueOf[emulated.BN254Fr](s), + Q: AffinePoint[emulated.BN254Fp]{ + X: emulated.ValueOf[emulated.BN254Fp](S.X), + Y: emulated.ValueOf[emulated.BN254Fp](S.Y), + }, + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err, "w=%d s=%s", w, s.String()) + } + } +} + +func TestScalarMulBaseCombP256(t *testing.T) { + assert := test.NewAssert(t) + // scalar field order of P-256 + r, _ := new(big.Int).SetString("ffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551", 16) + randFn := func() *big.Int { + rnd, _ := new(big.Int).SetString("3d6a4c9e1b5f2a7d8e9c0b1a2f3e4d5c6b7a8901234567890abcdef012345678", 16) + return rnd.Mod(rnd, r) + } + for _, w := range []int{8} { + for _, s := range combTestScalars(r, 1, randFn) { + sr := new(big.Int).Mod(s, r) + // compute the reference with the generic big.Int arithmetic used + // for table computation (P-256 has no gnark-crypto counterpart + // with the same API) + var Sx, Sy *big.Int + { + params := GetP256Params() + var fpp emulated.P256Fp + prime := fpp.Modulus() + if sr.Sign() == 0 { + Sx, Sy = big.NewInt(0), big.NewInt(0) + } else { + acc := &combAffine{x: params.Gx, y: params.Gy} + var err error + for i := sr.BitLen() - 2; i >= 0; i-- { + if acc, err = combDouble(acc, params.A, prime); err != nil { + t.Fatal(err) + } + if sr.Bit(i) == 1 { + if acc, err = combAdd(acc, &combAffine{x: params.Gx, y: params.Gy}, prime); err != nil { + t.Fatal(err) + } + } + } + Sx, Sy = acc.x, acc.y + } + } + circuit := ScalarMulBaseCombTest[emulated.P256Fp, emulated.P256Fr]{w: w} + witness := ScalarMulBaseCombTest[emulated.P256Fp, emulated.P256Fr]{ + S: emulated.ValueOf[emulated.P256Fr](sr), + Q: AffinePoint[emulated.P256Fp]{ + X: emulated.ValueOf[emulated.P256Fp](Sx), + Y: emulated.ValueOf[emulated.P256Fp](Sy), + }, + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err, "w=%d s=%s", w, sr.String()) + } + } +} + +// TestScalarMulBaseCombConstraints reports the constraint counts of the comb +// fixed-base scalar multiplication against the current ScalarMulBase. +func TestScalarMulBaseCombConstraints(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + for _, w := range []int{4, 6, 8, 10} { + circuit := ScalarMulBaseCombTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{w: w} + ccs, err := frontend.Compile(testCurve.ScalarField(), r1cs.NewBuilder, &circuit) + if err != nil { + t.Log("w =", w, "compile error:", err) + continue + } + assert.NoError(err) + t.Log("comb r1cs", "w =", w, "constraints =", ccs.GetNbConstraints()) + } + baseline := ScalarMulBaseTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{} + ccs, err := frontend.Compile(testCurve.ScalarField(), r1cs.NewBuilder, &baseline) + assert.NoError(err) + t.Log("baseline ScalarMulBase r1cs constraints =", ccs.GetNbConstraints()) + + // PLONKish counts + circuit := ScalarMulBaseCombTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{w: 8} + scsCcs, err := frontend.Compile(testCurve.ScalarField(), scs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("comb scs w=8 constraints =", scsCcs.GetNbConstraints()) + scsBase, err := frontend.Compile(testCurve.ScalarField(), scs.NewBuilder, &baseline) + assert.NoError(err) + t.Log("baseline ScalarMulBase scs constraints =", scsBase.GetNbConstraints()) +} diff --git a/std/algebra/emulated/sw_emulated/hints.go b/std/algebra/emulated/sw_emulated/hints.go index e436f992cc..f7b8c44db2 100644 --- a/std/algebra/emulated/sw_emulated/hints.go +++ b/std/algebra/emulated/sw_emulated/hints.go @@ -34,6 +34,8 @@ func GetHints() []solver.Hint { scalarMulHint, rationalReconstruct, rationalReconstructExt, + combRecodeHint, + combChainLambdaHint, } } diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 9ed6f6b662..4a9ca6539d 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -76,6 +76,9 @@ type Curve[Base, Scalars emulated.FieldParams] struct { addA bool eigenvalue *emulated.Element[Scalars] thirdRootOne *emulated.Element[Base] + + // combCache caches the fixed-base comb tables per window width. + combCache map[int]*combData } // Generator returns the base point of the curve. The method does not copy and diff --git a/std/math/emulated/field.go b/std/math/emulated/field.go index e23a29ca1f..65b0ef619e 100644 --- a/std/math/emulated/field.go +++ b/std/math/emulated/field.go @@ -173,6 +173,21 @@ func (f *Field[T]) NewElement(v any) *Element[T] { return c } +// UnsafeFromLimbs constructs an element from the given limbs without +// enforcing any range checks on them. The caller MUST guarantee that every +// limb value is strictly less than 2^BitsPerLimb, otherwise the result is +// unsound. It is useful in gadgets where the limbs are bounded by +// construction, for example when they are convex combinations of constant +// limbs selected by a one-hot vector. +func (f *Field[T]) UnsafeFromLimbs(limbVals []frontend.Variable) *Element[T] { + if len(limbVals) != int(f.fParams.NbLimbs()) { + panic("limb count mismatch") + } + cp := make([]frontend.Variable, len(limbVals)) + copy(cp, limbVals) + return f.newInternalElement(cp, 0) +} + // Zero returns zero as a constant. func (f *Field[T]) Zero() *Element[T] { f.zeroConstOnce.Do(func() { diff --git a/std/math/emulated/field_mul.go b/std/math/emulated/field_mul.go index a434090194..9c69896aef 100644 --- a/std/math/emulated/field_mul.go +++ b/std/math/emulated/field_mul.go @@ -920,6 +920,98 @@ type multivariate[T FieldParams] struct { // // The method returns the result of the evaluation. func (f *Field[T]) Eval(at [][]*Element[T], coefs []int) *Element[T] { + // it is the obvious case - when we don't have any inputs then we need to + // evaluate the zero polynomial which is always zero. + if len(at) == 0 && len(coefs) == 0 { + return f.Zero() + } + mv, allElems := f.polyMvPrep(at, coefs) + + // we call the hint to compute the result. The hint returns the reduced + // result, the quotient and the carries. + k, r, c, kNeg, err := f.callPolyMvHint(mv, allElems, false) + if err != nil { + panic(err) + } + + // finally, we store the deferred check which is performed later. The + // `mvCheck` implements the deferredChecker interface, so that we use the + // generic deferred check method. + mvc := mvCheck[T]{ + f: f, + mv: mv, + vals: allElems, + r: r, + k: k, + c: c, + kNeg: kNeg, + } + + f.deferredChecks = append(f.deferredChecks, &mvc) + + // Record operation for profiling + nbLimbs := 0 + for i := range allElems { + nbLimbs += len(allElems[i].Limbs) + } + nbLimbs += len(r.Limbs) + len(k.Limbs) + len(c.Limbs) + profile.RecordOperation("emulated.Eval", nbLimbs) + return r +} + +// AssertEvalIsZero asserts that the multivariate polynomial given by the terms +// at and coefficients coefs evaluates to zero modulo the field modulus. The +// interface is as in [Field.Eval]: the elements of the inner slices of at are +// multiplied together and summed with the corresponding coefficient. +// +// It is functionally equivalent to asserting that the result of [Field.Eval] +// is zero, but it is cheaper: the reduced remainder is never materialized as a +// witness (saving its allocation and range checks) and no separate equality +// check is needed. +// +// NB! This is experimental API. It does not check that computing the term +// wouldn't overflow the field. +func (f *Field[T]) AssertEvalIsZero(at [][]*Element[T], coefs []int) { + // zero polynomial is always zero + if len(at) == 0 && len(coefs) == 0 { + return + } + mv, allElems := f.polyMvPrep(at, coefs) + + // we call the hint to compute the quotient and the carries. As the + // remainder is asserted to be zero, the hint does not return it and we use + // the zero-limb constant zero element in the deferred check instead. + k, _, c, kNeg, err := f.callPolyMvHint(mv, allElems, true) + if err != nil { + panic(err) + } + + mvc := mvCheck[T]{ + f: f, + mv: mv, + vals: allElems, + r: f.Zero(), // constant zero on zero limbs + k: k, + c: c, + kNeg: kNeg, + } + + f.deferredChecks = append(f.deferredChecks, &mvc) + + // Record operation for profiling + nbLimbs := 0 + for i := range allElems { + nbLimbs += len(allElems[i].Limbs) + } + nbLimbs += len(k.Limbs) + len(c.Limbs) + profile.RecordOperation("emulated.AssertEvalIsZero", nbLimbs) +} + +// polyMvPrep prepares the multivariate polynomial evaluation of the terms at +// with coefficients coefs: it deduplicates the elements appearing in the +// terms, converts the terms into exponent form and ensures that all elements +// have their limb widths enforced. +func (f *Field[T]) polyMvPrep(at [][]*Element[T], coefs []int) (*multivariate[T], []*Element[T]) { if len(at) != len(coefs) { panic("terms and coefficients mismatch") } @@ -928,11 +1020,6 @@ func (f *Field[T]) Eval(at [][]*Element[T], coefs []int) *Element[T] { panic("coefficient math.MinInt overflows on negation") } } - // it is the obvious case - when we don't have any inputs then we need to - // evaluate the zero polynomial which is always zero. - if len(at) == 0 { - return f.Zero() - } // initialize the multivariate struct from the inputs. The current method // takes as input references to the elements. However, the hint function // works with solved values. So it would be better to work with the exact @@ -979,44 +1066,19 @@ func (f *Field[T]) Eval(at [][]*Element[T], coefs []int) *Element[T] { Terms: terms, Coefficients: coefs, } - - // we call the hint to compute the result. The hint returns the reduced - // result, the quotient and the carries. - k, r, c, kNeg, err := f.callPolyMvHint(mv, allElems) - if err != nil { - panic(err) - } - - // finally, we store the deferred check which is performed later. The - // `mvCheck` implements the deferredChecker interface, so that we use the - // generic deferred check method. - mvc := mvCheck[T]{ - f: f, - mv: mv, - vals: allElems, - r: r, - k: k, - c: c, - kNeg: kNeg, - } - - f.deferredChecks = append(f.deferredChecks, &mvc) - - // Record operation for profiling - nbLimbs := 0 - for i := range allElems { - nbLimbs += len(allElems[i].Limbs) - } - nbLimbs += len(r.Limbs) + len(k.Limbs) + len(c.Limbs) - profile.RecordOperation("emulated.Eval", nbLimbs) - return r + return mv, allElems } // callPolyMvHint computes the multivariate evaluation given by mv at at. It // returns the remainder (reduced result), the quotient and the carries. The // computation is performed inside a hint, so it is the callers responsibility to // perform the deferred multiplication check. -func (f *Field[T]) callPolyMvHint(mv *multivariate[T], at []*Element[T]) (quo, rem, carries *Element[T], kNeg frontend.Variable, err error) { +// +// When assertZero is set, the evaluation is asserted to be zero modulo the +// field modulus: the hint does not output the remainder limbs (the returned +// rem is the zero-limb constant zero) and it errors at solving time if the +// evaluation is not divisible by the modulus. +func (f *Field[T]) callPolyMvHint(mv *multivariate[T], at []*Element[T], assertZero bool) (quo, rem, carries *Element[T], kNeg frontend.Variable, err error) { // first compute the length of the result so that we know how many bits we need for the quotient. nbLimbs, nbBits := f.fParams.NbLimbs(), f.fParams.BitsPerLimb() modBits := uint(f.fParams.Modulus().BitLen()) @@ -1026,14 +1088,17 @@ func (f *Field[T]) callPolyMvHint(mv *multivariate[T], at []*Element[T]) (quo, r nbQuoLimbs = (quoSize - modBits + nbBits) / nbBits } nbRemLimbs := nbLimbs + if assertZero { + nbRemLimbs = 0 + } nbCarryLimbs := nbMultiplicationResLimbs(int(nbQuoLimbs), int(nbLimbs)) - 1 - nbHintInputs := 6 + len(mv.Coefficients) + len(at)*len(mv.Terms) + len(mv.Coefficients) + len(f.Modulus().Limbs) + nbHintInputs := 7 + len(mv.Coefficients) + len(at)*len(mv.Terms) + len(mv.Coefficients) + len(f.Modulus().Limbs) for i := range at { nbHintInputs += len(at[i].Limbs) + 1 } hintInputs := make([]frontend.Variable, 0, nbHintInputs) - hintInputs = append(hintInputs, nbBits, nbLimbs, len(mv.Terms), len(at), nbQuoLimbs, nbCarryLimbs) + hintInputs = append(hintInputs, nbBits, nbLimbs, len(mv.Terms), len(at), nbQuoLimbs, nbRemLimbs, nbCarryLimbs) // store per-coefficient signs: 0 = positive, 1 = negative for _, c := range mv.Coefficients { if c < 0 { @@ -1071,7 +1136,11 @@ func (f *Field[T]) callPolyMvHint(mv *multivariate[T], at []*Element[T]) (quo, r return } quo = f.packLimbs(ret[:nbQuoLimbs], false) - rem = f.packLimbs(ret[nbQuoLimbs:nbQuoLimbs+nbRemLimbs], true) + if assertZero { + rem = f.Zero() + } else { + rem = f.packLimbs(ret[nbQuoLimbs:nbQuoLimbs+nbRemLimbs], true) + } carries = f.newInternalElement(ret[nbQuoLimbs+nbRemLimbs:nbQuoLimbs+nbRemLimbs+uint(nbCarryLimbs)], 0) kNeg = ret[nbQuoLimbs+nbRemLimbs+uint(nbCarryLimbs)] f.api.AssertIsBoolean(kNeg) @@ -1248,7 +1317,7 @@ func (f *Field[T]) polyMvEvalQuoSize(mv *multivariate[T], at []*Element[T]) (quo // called directly, but rather through [Field.callPolyMvHint] method which // handles the input packing and output unpacking. func polyMvHint(mod *big.Int, inputs, outputs []*big.Int) error { - if len(inputs) < 7 { + if len(inputs) < 8 { return errors.New("not enough inputs") } var ( @@ -1257,9 +1326,12 @@ func polyMvHint(mod *big.Int, inputs, outputs []*big.Int) error { nbTerms = int(inputs[2].Int64()) nbVars = int(inputs[3].Int64()) nbQuoLimbs = int(inputs[4].Int64()) - nbRemLimbs = nbLimbs - nbCarryLimbs = int(inputs[5].Int64()) + nbRemLimbs = int(inputs[5].Int64()) + nbCarryLimbs = int(inputs[6].Int64()) ) + // nbRemLimbs == 0 indicates that the caller asserts the evaluation to be + // zero modulo the modulus: no remainder limbs are output and a non-zero + // remainder is a solving error. if len(outputs) != nbQuoLimbs+nbRemLimbs+nbCarryLimbs+1 { return errors.New("output length mismatch") } @@ -1272,7 +1344,7 @@ func polyMvHint(mod *big.Int, inputs, outputs []*big.Int) error { outPtr += nbCarryLimbs kNegOut := outputs[outPtr] // read per-coefficient signs: 0 = positive, 1 = negative - ptr := 6 + ptr := 7 signs := make([]int, nbTerms) for i := range signs { signs[i] = int(inputs[ptr].Int64()) @@ -1367,7 +1439,13 @@ func polyMvHint(mod *big.Int, inputs, outputs []*big.Int) error { if err := limbs.Decompose(quo, uint(nbBits), quoLimbs); err != nil { return fmt.Errorf("decompose quo: %w", err) } - if err := limbs.Decompose(rem, uint(nbBits), remLimbs); err != nil { + if nbRemLimbs == 0 { + // the caller asserts the evaluation to be zero modulo the modulus. A + // non-zero remainder means the assertion cannot be satisfied. + if rem.Sign() != 0 { + return errors.New("asserted zero evaluation has non-zero remainder") + } + } else if err := limbs.Decompose(rem, uint(nbBits), remLimbs); err != nil { return fmt.Errorf("decompose rem: %w", err) } From 26b465c63443df88c6ca07c31da3d2f9cfe8d955 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Thu, 23 Jul 2026 21:20:41 -0400 Subject: [PATCH 02/15] perf: use comb in JointScalarMulBase --- std/algebra/emulated/sw_emulated/fixedbase.go | 6 ++ .../emulated/sw_emulated/fixedbase_test.go | 53 ++++++++++++++++++ .../sw_emulated/joint_comb_count_test.go | 55 +++++++++++++++++++ std/algebra/emulated/sw_emulated/point.go | 41 +++++++++++--- 4 files changed, 146 insertions(+), 9 deletions(-) create mode 100644 std/algebra/emulated/sw_emulated/joint_comb_count_test.go diff --git a/std/algebra/emulated/sw_emulated/fixedbase.go b/std/algebra/emulated/sw_emulated/fixedbase.go index 8eaff714bf..741aad9719 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase.go +++ b/std/algebra/emulated/sw_emulated/fixedbase.go @@ -46,6 +46,12 @@ import ( // checks) plus a y materialization (one check). The y-coordinate is // materialized once, before the final complete addition. +// combDefaultWindow is the default window width of the fixed-base comb. With +// 64-bit limb emulation it is supported by all built-in curves (the recoded +// scalar fits the scalar-field limb capacity) and is close to the +// constraint-count optimum in R1CS. +const combDefaultWindow = 8 + // combData holds the compile-time data of the comb: the constant window // tables and the derived parameters. type combData struct { diff --git a/std/algebra/emulated/sw_emulated/fixedbase_test.go b/std/algebra/emulated/sw_emulated/fixedbase_test.go index a1dd0b1418..e7d7fbb4a5 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase_test.go +++ b/std/algebra/emulated/sw_emulated/fixedbase_test.go @@ -185,3 +185,56 @@ func TestScalarMulBaseCombConstraints(t *testing.T) { assert.NoError(err) t.Log("baseline ScalarMulBase scs constraints =", scsBase.GetNbConstraints()) } + +type jointScalarMulBaseCompleteTest[T, S emulated.FieldParams] struct { + P AffinePoint[T] + S1, S2 emulated.Element[S] + Q AffinePoint[T] +} + +func (c *jointScalarMulBaseCompleteTest[T, S]) Define(api frontend.API) error { + cr, err := New[T, S](api, GetCurveParams[T]()) + if err != nil { + return err + } + res := cr.JointScalarMulBase(&c.P, &c.S2, &c.S1) + cr.AssertIsEqual(res, &c.Q) + return nil +} + +// TestJointScalarMulBaseComplete exercises the comb-based complete path of +// JointScalarMulBase, including the zero fixed-base scalar. +func TestJointScalarMulBaseComplete(t *testing.T) { + assert := test.NewAssert(t) + _, g := secp256k1.Generators() + var p secp256k1.G1Affine + p.Double(&g) + r := fr_secp.Modulus() + randFn := func() *big.Int { + var rnd fr_secp.Element + _, _ = rnd.SetRandom() + return rnd.BigInt(new(big.Int)) + } + s2 := randFn() + for _, s1 := range combTestScalars(r, 2, randFn) { + var sm1, sm2, S secp256k1.G1Affine + sm1.ScalarMultiplication(&g, s1) + sm2.ScalarMultiplication(&p, s2) + S.Add(&sm1, &sm2) + circuit := jointScalarMulBaseCompleteTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{} + witness := jointScalarMulBaseCompleteTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{ + S1: emulated.ValueOf[emulated.Secp256k1Fr](s1), + S2: emulated.ValueOf[emulated.Secp256k1Fr](s2), + P: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](p.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](p.Y), + }, + Q: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](S.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](S.Y), + }, + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err, "s1=%s", s1.String()) + } +} diff --git a/std/algebra/emulated/sw_emulated/joint_comb_count_test.go b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go new file mode 100644 index 0000000000..6fda4793de --- /dev/null +++ b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go @@ -0,0 +1,55 @@ +package sw_emulated + +import ( + "testing" + + "github.com/consensys/gnark/frontend" + "github.com/consensys/gnark/frontend/cs/r1cs" + "github.com/consensys/gnark/std/algebra/algopts" + "github.com/consensys/gnark/std/math/emulated" + "github.com/consensys/gnark/test" +) + +type jointCombCount[T, S emulated.FieldParams] struct { + P AffinePoint[T] + S1, S2 emulated.Element[S] + Q AffinePoint[T] + mode int // 0: baseline complete, 1: baseline incomplete, 2: comb-composed +} + +func (c *jointCombCount[T, S]) Define(api frontend.API) error { + cr, err := New[T, S](api, GetCurveParams[T]()) + if err != nil { + return err + } + var res *AffinePoint[T] + switch c.mode { + case 0: + res = cr.JointScalarMulBase(&c.P, &c.S2, &c.S1) + case 1: + res = cr.JointScalarMulBase(&c.P, &c.S2, &c.S1, algopts.WithIncompleteArithmetic()) + case 2: + sm1 := cr.scalarMulBaseComb(&c.S1, 8) + sm2 := cr.ScalarMul(&c.P, &c.S2) + res = cr.AddUnified(sm1, sm2) + case 3: + sm1 := cr.scalarMulBaseComb(&c.S1, 8) + sm2 := cr.ScalarMul(&c.P, &c.S2, algopts.WithIncompleteArithmetic()) + res = cr.AddUnified(sm1, sm2) + } + cr.AssertIsEqual(res, &c.Q) + return nil +} + +func TestJointCombCount(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + for mode := 0; mode <= 3; mode++ { + circuit := jointCombCount[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{mode: mode} + ccs, err := frontend.Compile(testCurve.ScalarField(), r1cs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("mode", mode, "r1cs constraints =", ccs.GetNbConstraints()) + } +} diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 4a9ca6539d..39791ab02c 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -1236,20 +1236,22 @@ func (c *Curve[B, S]) jointScalarMulGLVUnsafe(Q, R *AffinePoint[B], s, t *emulat // ScalarMulBase computes [s]g and returns it where g is the fixed curve generator. It doesn't modify p nor s. // -// By default, uses complete arithmetic. +// It uses the fixed-base comb method with compile-time constant window tables +// (see [Curve.scalarMulBaseComb]), which uses complete arithmetic and +// correctly handles the zero scalar. The [algopts.WithIncompleteArithmetic] +// option is a no-op for this method as the comb method is both complete and +// cheaper than the incomplete variable-base fallbacks. // -// ⚠️ When [algopts.WithIncompleteArithmetic] is set, the exact exceptional set -// depends on the scalar-multiplication algorithm selected for the current -// curve: +// For custom curve parameters where the comb tables cannot be constructed, it +// falls back to the variable-base scalar multiplication with the generator: // - curves without an efficient endomorphism inherit the documented // exceptional set of [Curve.scalarMulFakeGLV] -// and currently include P-256, P-384 and STARK curve // - curves with an efficient endomorphism inherit the documented exceptional -// set of [Curve.scalarMulGLVAndFakeGLV] and currently include BN254, -// BLS12-381, BW6-761 and secp256k1. -// -// ScalarMul calls scalarMulBaseGeneric or scalarMulGLVAndFakeGLV depending on whether an efficient endomorphism is available. +// set of [Curve.scalarMulGLVAndFakeGLV]. func (c *Curve[B, S]) ScalarMulBase(s *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] { + if _, err := c.combData(combDefaultWindow); err == nil { + return c.scalarMulBaseComb(s, combDefaultWindow) + } if c.eigenvalue != nil && c.thirdRootOne != nil { return c.scalarMulGLVAndFakeGLV(c.Generator(), s, opts...) @@ -1286,6 +1288,27 @@ func (c *Curve[B, S]) ScalarMulBase(s *emulated.Element[S], opts ...algopts.Alge // The [EVM] specifies these checks, which are performed on the zkEVM // arithmetization side before calling the circuit that uses this method. func (c *Curve[B, S]) JointScalarMulBase(p *AffinePoint[B], s2, s1 *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] { + cfg, err := algopts.NewConfig(opts...) + if err != nil { + panic(fmt.Sprintf("parse opts: %v", err)) + } + if _, cerr := c.combData(combDefaultWindow); cerr == nil && !cfg.IncompleteArithmetic { + // In complete mode, compute the fixed-base part with the comb method + // (complete, handles s1 = 0) and the variable-base part separately, + // and merge with the complete addition. This is cheaper than two + // variable-base scalar multiplications. + // + // We do NOT take this path in incomplete mode: composing the comb + // with the incomplete variable-base [Curve.ScalarMul] would be + // cheaper still, but the incomplete scalar multiplication has a + // non-negligible exceptional set when p has a small known relation + // to the generator (e.g. p = [2]g fails for a noticeable fraction of + // scalars), whereas the joint Shamir-based algorithm below handles + // those points. + sm1 := c.scalarMulBaseComb(s1, combDefaultWindow) + sm2 := c.ScalarMul(p, s2, opts...) + return c.AddUnified(sm1, sm2) + } return c.jointScalarMul(c.Generator(), p, s1, s2, opts...) } From b474f11bc3c8661cc2a984517fdb9131df60c066 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 10:18:07 -0400 Subject: [PATCH 03/15] perf: generalize the comb beyond the generator --- std/algebra/emulated/sw_emulated/fixedbase.go | 81 ++++++++++++- .../emulated/sw_emulated/fixedbase_test.go | 113 ++++++++++++++++++ .../sw_emulated/joint_comb_count_test.go | 31 +++++ std/algebra/emulated/sw_emulated/point.go | 23 +++- std/math/emulated/field.go | 12 ++ 5 files changed, 252 insertions(+), 8 deletions(-) diff --git a/std/algebra/emulated/sw_emulated/fixedbase.go b/std/algebra/emulated/sw_emulated/fixedbase.go index 741aad9719..09fd3111bf 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase.go +++ b/std/algebra/emulated/sw_emulated/fixedbase.go @@ -120,14 +120,65 @@ func combDouble(p *combAffine, a, prime *big.Int) (*combAffine, error) { return &combAffine{x: xr, y: yr}, nil } +// combCheckPoint checks that the base point (gx, gy) is a finite point on the +// curve y² = x³ + ax + b of prime order r: it verifies the curve equation and +// that [r−1](gx, gy) = −(gx, gy). Both are required for the comb soundness +// argument: the window tables must contain no point at infinity and the +// partial-sum collision analysis works modulo the order of the base point. +func combCheckPoint(gx, gy, a, b, prime, r *big.Int) error { + if gx.Sign() == 0 && gy.Sign() == 0 { + return errors.New("base point is the point at infinity") + } + lhs := new(big.Int).Mul(gy, gy) + lhs.Mod(lhs, prime) + rhs := new(big.Int).Mul(gx, gx) + rhs.Mul(rhs, gx) + if a != nil && a.Sign() != 0 { + rhs.Add(rhs, new(big.Int).Mul(a, gx)) + } + if b != nil { + rhs.Add(rhs, b) + } + rhs.Mod(rhs, prime) + if lhs.Cmp(rhs) != 0 { + return errors.New("base point is not on the curve") + } + // check [r−1]P = −P by double-and-add. The intermediate partial sums are + // [m]P with 0 < m < r−1, so if ord(P) = r the chain never encounters the + // point at infinity nor an x-collision; conversely any such failure means + // ord(P) ≠ r and we reject. + P := &combAffine{x: gx, y: gy} + e := new(big.Int).Sub(r, big.NewInt(1)) + acc := P + var err error + for i := e.BitLen() - 2; i >= 0; i-- { + if acc, err = combDouble(acc, a, prime); err != nil { + return fmt.Errorf("base point order check: %w", err) + } + if e.Bit(i) == 1 { + if acc, err = combAdd(acc, P, prime); err != nil { + return fmt.Errorf("base point order check: %w", err) + } + } + } + negP := combNeg(P, prime) + if acc.x.Cmp(negP.x) != 0 || acc.y.Cmp(negP.y) != 0 { + return errors.New("base point does not have prime order r") + } + return nil +} + // computeCombData computes the comb tables for the curve y² = x³ + ax + b // over the prime field of modulus prime, with base point (gx, gy) of prime // order r, window width w and a recoded-scalar capacity of scalarCap bits // (the recomposition capacity of the scalar field emulation). -func computeCombData(gx, gy, a, prime, r *big.Int, w int, scalarCap int) (*combData, error) { +func computeCombData(gx, gy, a, b, prime, r *big.Int, w int, scalarCap int) (*combData, error) { if w < 2 || w > 14 { return nil, fmt.Errorf("unsupported window width %d", w) } + if err := combCheckPoint(gx, gy, a, b, prime, r); err != nil { + return nil, err + } rBits := r.BitLen() nw := (rBits + w - 1) / w if nw < 2 { @@ -208,21 +259,32 @@ func computeCombData(gx, gy, a, prime, r *big.Int, w int, scalarCap int) (*combD }, nil } -// combData returns the (cached) comb tables for the given window width. +// combData returns the (cached) comb tables for the generator and the given +// window width. func (c *Curve[B, S]) combData(w int) (*combData, error) { - if d, ok := c.combCache[w]; ok { + return c.combDataFor(c.params.Gx, c.params.Gy, w) +} + +// combDataFor returns the (cached) comb tables for the given constant base +// point and window width. It returns an error when the tables cannot be +// constructed: unsupported window width, recoded scalar exceeding the scalar +// field emulation capacity, or a base point which is not a finite curve point +// of prime order r. +func (c *Curve[B, S]) combDataFor(gx, gy *big.Int, w int) (*combData, error) { + key := fmt.Sprintf("%d|%s|%s", w, gx.Text(16), gy.Text(16)) + if d, ok := c.combCache[key]; ok { return d, nil } var fp B var fr S - d, err := computeCombData(c.params.Gx, c.params.Gy, c.params.A, fp.Modulus(), fr.Modulus(), w, int(fr.NbLimbs()*fr.BitsPerLimb())) + d, err := computeCombData(gx, gy, c.params.A, c.params.B, fp.Modulus(), fr.Modulus(), w, int(fr.NbLimbs()*fr.BitsPerLimb())) if err != nil { return nil, err } if c.combCache == nil { - c.combCache = make(map[int]*combData) + c.combCache = make(map[string]*combData) } - c.combCache[w] = d + c.combCache[key] = d return d, nil } @@ -346,6 +408,13 @@ func (c *Curve[B, S]) scalarMulBaseComb(s *emulated.Element[S], w int) *AffinePo if err != nil { panic(fmt.Sprintf("comb data: %v", err)) } + return c.scalarMulComb(d, s) +} + +// scalarMulComb computes [s]P where P is the compile-time constant base point +// of the comb tables d. It returns (0,0) when s ≡ 0 (mod r). +func (c *Curve[B, S]) scalarMulComb(d *combData, s *emulated.Element[S]) *AffinePoint[B] { + w := d.w n, nw := d.n, d.nw // scalar recode: b0 = parity of s, c = (k' + 2^n − 1)/2 with diff --git a/std/algebra/emulated/sw_emulated/fixedbase_test.go b/std/algebra/emulated/sw_emulated/fixedbase_test.go index e7d7fbb4a5..6e95baf5b9 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase_test.go +++ b/std/algebra/emulated/sw_emulated/fixedbase_test.go @@ -4,6 +4,8 @@ import ( "math/big" "testing" + bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381" + fr_bls381 "github.com/consensys/gnark-crypto/ecc/bls12-381/fr" "github.com/consensys/gnark-crypto/ecc/bn254" fr_bn "github.com/consensys/gnark-crypto/ecc/bn254/fr" "github.com/consensys/gnark-crypto/ecc/secp256k1" @@ -238,3 +240,114 @@ func TestJointScalarMulBaseComplete(t *testing.T) { assert.NoError(err, "s1=%s", s1.String()) } } + +type scalarMulConstPointTest[T, S emulated.FieldParams] struct { + S emulated.Element[S] + Q AffinePoint[T] + px *big.Int + py *big.Int +} + +func (c *scalarMulConstPointTest[T, S]) Define(api frontend.API) error { + cr, err := New[T, S](api, GetCurveParams[T]()) + if err != nil { + return err + } + P := AffinePoint[T]{ + X: emulated.ValueOf[T](c.px), + Y: emulated.ValueOf[T](c.py), + } + res := cr.ScalarMul(&P, &c.S) + cr.AssertIsEqual(res, &c.Q) + return nil +} + +// TestScalarMulConstPoint exercises the automatic comb dispatch in ScalarMul +// for compile-time constant points. +func TestScalarMulConstPoint(t *testing.T) { + assert := test.NewAssert(t) + _, g := secp256k1.Generators() + // constant point P = [12345]G + var P secp256k1.G1Affine + P.ScalarMultiplication(&g, big.NewInt(12345)) + px, py := P.X.BigInt(new(big.Int)), P.Y.BigInt(new(big.Int)) + r := fr_secp.Modulus() + randFn := func() *big.Int { + var rnd fr_secp.Element + _, _ = rnd.SetRandom() + return rnd.BigInt(new(big.Int)) + } + for _, s := range combTestScalars(r, 2, randFn) { + var S secp256k1.G1Affine + S.ScalarMultiplication(&P, s) + circuit := scalarMulConstPointTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{px: px, py: py} + witness := scalarMulConstPointTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{ + px: px, py: py, + S: emulated.ValueOf[emulated.Secp256k1Fr](s), + Q: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](S.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](S.Y), + }, + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err, "s=%s", s.String()) + } +} + +// TestScalarMulConstPointBLS12381 checks the comb dispatch on a cofactor +// curve: a subgroup point uses the comb, and a curve point outside the +// r-torsion is rejected by the order check and falls back to the generic +// variable-base path (compilation must succeed). +func TestScalarMulConstPointBLS12381(t *testing.T) { + assert := test.NewAssert(t) + _, _, g, _ := bls12381.Generators() + var P bls12381.G1Affine + P.ScalarMultiplication(&g, big.NewInt(987654321)) + px, py := P.X.BigInt(new(big.Int)), P.Y.BigInt(new(big.Int)) + var rnd fr_bls381.Element + _, _ = rnd.SetRandom() + s := rnd.BigInt(new(big.Int)) + var S bls12381.G1Affine + S.ScalarMultiplication(&P, s) + circuit := scalarMulConstPointTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{px: px, py: py} + witness := scalarMulConstPointTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{ + px: px, py: py, + S: emulated.ValueOf[emulated.BLS12381Fr](s), + Q: AffinePoint[emulated.BLS12381Fp]{ + X: emulated.ValueOf[emulated.BLS12381Fp](S.X), + Y: emulated.ValueOf[emulated.BLS12381Fp](S.Y), + }, + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err) + + // non-r-torsion curve point: search a valid x with y² = x³ + 4 a QR and + // check it is rejected by the comb order check (cofactor > 1 makes a + // random curve point land outside the subgroup w.h.p.). + var fpp emulated.BLS12381Fp + prime := fpp.Modulus() + exp := new(big.Int).Add(prime, big.NewInt(1)) + exp.Rsh(exp, 2) // (p+1)/4, p ≡ 3 mod 4 + found := false + for x := int64(1); x < 50 && !found; x++ { + xx := big.NewInt(x) + rhs := new(big.Int).Exp(xx, big.NewInt(3), prime) + rhs.Add(rhs, big.NewInt(4)).Mod(rhs, prime) + y := new(big.Int).Exp(rhs, exp, prime) + check := new(big.Int).Mul(y, y) + check.Mod(check, prime) + if check.Cmp(rhs) != 0 { + continue + } + // on curve; must not be in the r-torsion for this test to be + // meaningful + var frr emulated.BLS12381Fr + if err := combCheckPoint(xx, y, big.NewInt(0), big.NewInt(4), prime, frr.Modulus()); err == nil { + continue + } + // combCheckPoint rejecting the point is exactly what makes + // combDataFor fall back to the generic path for it. + found = true + } + assert.True(found, "expected to find a non-subgroup curve point") +} diff --git a/std/algebra/emulated/sw_emulated/joint_comb_count_test.go b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go index 6fda4793de..ca2492b6d1 100644 --- a/std/algebra/emulated/sw_emulated/joint_comb_count_test.go +++ b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go @@ -53,3 +53,34 @@ func TestJointCombCount(t *testing.T) { t.Log("mode", mode, "r1cs constraints =", ccs.GetNbConstraints()) } } + +type constPointCountCircuit struct { + S emulated.Element[emulated.Secp256k1Fr] + Q AffinePoint[emulated.Secp256k1Fp] +} + +func (c *constPointCountCircuit) Define(api frontend.API) error { + cr, err := New[emulated.Secp256k1Fp, emulated.Secp256k1Fr](api, GetCurveParams[emulated.Secp256k1Fp]()) + if err != nil { + return err + } + // constant point (not the generator) + P := AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp]("89565891926547004231252920425935692360644145829622209833684329913297188986597"), + Y: emulated.ValueOf[emulated.Secp256k1Fp]("12158399299693830322967808612713398636155367887041628176798871954788371653930"), + } + res := cr.ScalarMul(&P, &c.S) + cr.AssertIsEqual(res, &c.Q) + return nil +} + +func TestConstPointScalarMulCount(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + circuit := constPointCountCircuit{} + ccs, err := frontend.Compile(testCurve.ScalarField(), r1cs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("constant-point ScalarMul r1cs constraints =", ccs.GetNbConstraints()) +} diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 39791ab02c..5046f47056 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -77,8 +77,9 @@ type Curve[Base, Scalars emulated.FieldParams] struct { eigenvalue *emulated.Element[Scalars] thirdRootOne *emulated.Element[Base] - // combCache caches the fixed-base comb tables per window width. - combCache map[int]*combData + // combCache caches the fixed-base comb tables per base point and window + // width. + combCache map[string]*combData } // Generator returns the base point of the curve. The method does not copy and @@ -646,7 +647,25 @@ func (c *Curve[B, S]) muxY8Signed(signBit frontend.Variable, selector frontend.V // N.B. For scalarMulGLVAndFakeGLV, the result is undefined when the input point is // not on the prime order subgroup. For scalarMulFakeGLV the result is well // defined for any point on the curve +// +// When p is a compile-time constant point of prime order r (for example a +// point from a fixed verification key or SRS), the method automatically uses +// the fixed-base comb method with precomputed tables (see +// [Curve.scalarMulComb]), which uses complete arithmetic and is significantly +// cheaper. In this case the [algopts.WithIncompleteArithmetic] option is a +// no-op. func (c *Curve[B, S]) ScalarMul(p *AffinePoint[B], s *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] { + if px, ok := c.baseApi.ConstantValue(&p.X); ok { + if py, ok := c.baseApi.ConstantValue(&p.Y); ok { + // constant point: try the fixed-base comb. combDataFor verifies + // at compile time that (px, py) is a finite curve point of prime + // order r; otherwise we fall back to the variable-base methods + // below which have no such requirement. + if d, err := c.combDataFor(px, py, combDefaultWindow); err == nil { + return c.scalarMulComb(d, s) + } + } + } if c.eigenvalue != nil && c.thirdRootOne != nil { return c.scalarMulGLVAndFakeGLV(p, s, opts...) diff --git a/std/math/emulated/field.go b/std/math/emulated/field.go index 65b0ef619e..c22f4d8f9b 100644 --- a/std/math/emulated/field.go +++ b/std/math/emulated/field.go @@ -313,6 +313,18 @@ func (f *Field[T]) enforceWidthConditional(a *Element[T]) (didConstrain bool) { return } +// ConstantValue returns the constant value of the element modulo the field +// modulus and a boolean indicating if the element is in fact a compile-time +// constant. It allows gadgets to specialize for constant inputs (for example +// scalar multiplication by a constant point can use precomputed tables). +func (f *Field[T]) ConstantValue(v *Element[T]) (*big.Int, bool) { + c, ok := f.constantValue(v) + if !ok { + return nil, false + } + return c.Mod(c, f.fParams.Modulus()), true +} + func (f *Field[T]) constantValue(v *Element[T]) (*big.Int, bool) { // this case happens when we have called [ValueOf] inside a circuit as // [Element.Initialize] has not been called (Limbs are nil). In this case, From 3ef6b3e39e56881bae9bca5e2e394111d043a59b Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 10:43:11 -0400 Subject: [PATCH 04/15] perf: optimize variable-base scalar mul --- std/algebra/emulated/sw_emulated/fixedbase.go | 7 + std/algebra/emulated/sw_emulated/hints.go | 5 + .../sw_emulated/joint_comb_count_test.go | 35 ++++ std/algebra/emulated/sw_emulated/point.go | 74 ++++--- std/algebra/emulated/sw_emulated/slopes.go | 198 ++++++++++++++++++ 5 files changed, 285 insertions(+), 34 deletions(-) create mode 100644 std/algebra/emulated/sw_emulated/slopes.go diff --git a/std/algebra/emulated/sw_emulated/fixedbase.go b/std/algebra/emulated/sw_emulated/fixedbase.go index 09fd3111bf..6874b99018 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase.go +++ b/std/algebra/emulated/sw_emulated/fixedbase.go @@ -6,6 +6,7 @@ import ( "math/big" "github.com/consensys/gnark/frontend" + "github.com/consensys/gnark/internal/smallfields" limbs "github.com/consensys/gnark/std/internal/limbcomposition" "github.com/consensys/gnark/std/math/emulated" ) @@ -271,6 +272,12 @@ func (c *Curve[B, S]) combData(w int) (*combData, error) { // field emulation capacity, or a base point which is not a finite curve point // of prime order r. func (c *Curve[B, S]) combDataFor(gx, gy *big.Int, w int) (*combData, error) { + if smallfields.IsSmallField(c.api.Compiler().Field()) { + // in small-field mode the emulated field instance uses a different + // limb layout than the static field parameters which the comb + // selector and recode assume; fall back to the generic methods. + return nil, errors.New("fixed-base comb unsupported over small fields") + } key := fmt.Sprintf("%d|%s|%s", w, gx.Text(16), gy.Text(16)) if d, ok := c.combCache[key]; ok { return d, nil diff --git a/std/algebra/emulated/sw_emulated/hints.go b/std/algebra/emulated/sw_emulated/hints.go index f7b8c44db2..b15568e87c 100644 --- a/std/algebra/emulated/sw_emulated/hints.go +++ b/std/algebra/emulated/sw_emulated/hints.go @@ -36,6 +36,11 @@ func GetHints() []solver.Hint { rationalReconstructExt, combRecodeHint, combChainLambdaHint, + ratioHint, + tangentHint, + tangentHintA, + unifiedSlopeHint, + bjSlopeHint, } } diff --git a/std/algebra/emulated/sw_emulated/joint_comb_count_test.go b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go index ca2492b6d1..9a47bb0a10 100644 --- a/std/algebra/emulated/sw_emulated/joint_comb_count_test.go +++ b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go @@ -84,3 +84,38 @@ func TestConstPointScalarMulCount(t *testing.T) { assert.NoError(err) t.Log("constant-point ScalarMul r1cs constraints =", ccs.GetNbConstraints()) } + +type varMulCountCircuit struct { + P AffinePoint[emulated.Secp256k1Fp] + S emulated.Element[emulated.Secp256k1Fr] + Q AffinePoint[emulated.Secp256k1Fp] + incomplete bool +} + +func (c *varMulCountCircuit) Define(api frontend.API) error { + cr, err := New[emulated.Secp256k1Fp, emulated.Secp256k1Fr](api, GetCurveParams[emulated.Secp256k1Fp]()) + if err != nil { + return err + } + var res *AffinePoint[emulated.Secp256k1Fp] + if c.incomplete { + res = cr.ScalarMul(&c.P, &c.S, algopts.WithIncompleteArithmetic()) + } else { + res = cr.ScalarMul(&c.P, &c.S) + } + cr.AssertIsEqual(res, &c.Q) + return nil +} + +func TestVarMulCount(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + for _, inc := range []bool{false, true} { + circuit := varMulCountCircuit{incomplete: inc} + ccs, err := frontend.Compile(testCurve.ScalarField(), r1cs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("var-base ScalarMul incomplete =", inc, "r1cs constraints =", ccs.GetNbConstraints()) + } +} diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 5046f47056..7f34929445 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -195,7 +195,7 @@ func (c *Curve[B, S]) add(p, q *AffinePoint[B]) *AffinePoint[B] { // compute λ = (q.y-p.y)/(q.x-p.x) qypy := c.baseApi.Sub(&q.Y, &p.Y) qxpx := c.baseApi.Sub(&q.X, &p.X) - λ := c.baseApi.Div(qypy, qxpx) + λ := c.assertedRatio(qypy, qxpx) // xr = λ²-p.x-q.x xr := c.baseApi.Eval([][]*emulated.Element[B]{{λ, λ}, {c.baseApi.Add(&p.X, &q.X)}}, []int{1, -1}) @@ -277,15 +277,28 @@ func (c *Curve[B, S]) AddUnified(p, q *AffinePoint[B]) *AffinePoint[B] { // so we can safely divide by a dummy 1 and force λ to 0. numChord := c.baseApi.Sub(&q.Y, &p.Y) denChord := xDiff - xx := c.baseApi.MulMod(&p.X, &p.X) - numTangent := c.baseApi.MulConst(xx, big.NewInt(3)) denTangent := c.baseApi.MulConst(&p.Y, big.NewInt(2)) - num := c.baseApi.Select(xEqual, numTangent, numChord) den := c.baseApi.Select(xEqual, denTangent, denChord) denIsZero := c.baseApi.IsZero(den) denSafe := c.baseApi.Select(denIsZero, c.baseApi.One(), den) - λ := c.baseApi.Div(num, denSafe) + // witness the unified slope and certify it with a single deferred + // zero-assertion, blending the chord and tangent numerators with the + // xEqual indicator instead of materializing 3x² and selecting: + // λ·denSafe − xEqual·3x² − (1−xEqual)·(y2−y1) ≡ 0 + // λ remains pinned in all cases; when denIsZero the assertion is + // satisfied by the hint value and λ is discarded by the select below. + lams, err := c.baseApi.NewHint(unifiedSlopeHint, 1, &p.X, &p.Y, &q.X, &q.Y) + if err != nil { + panic(fmt.Sprintf("unified slope hint: %v", err)) + } + λ := lams[0] + zx := c.baseApi.FromBits(xEqual) + nzx := c.baseApi.Sub(c.baseApi.One(), zx) + c.baseApi.AssertEvalIsZero( + [][]*emulated.Element[B]{{λ, denSafe}, {zx, &p.X, &p.X}, {nzx, numChord}}, + []int{1, -3, -1}, + ) λ = c.baseApi.Select(denIsZero, c.baseApi.Zero(), λ) // compute the result point from λ @@ -324,17 +337,23 @@ func (c *Curve[B, S]) AddUnified(p, q *AffinePoint[B]) *AffinePoint[B] { // the j=0 branch, to avoid turning O + Q into O when Q.Y = 0. // --------------------------------------------------------------- - // λ = ((p.x+q.x)² - p.x*q.x + a)/(p.y + q.y) - pxqx := c.baseApi.MulMod(&p.X, &q.X) + // λ = ((p.x+q.x)² - p.x*q.x + a)/(p.y + q.y), certified by a single + // deferred zero-assertion without materializing the numerator: + // λ·denum − (p.x+q.x)² + p.x·q.x − a ≡ 0 pxplusqx := c.baseApi.Add(&p.X, &q.X) - num := c.baseApi.MulMod(pxplusqx, pxplusqx) - num = c.baseApi.Sub(num, pxqx) - num = c.baseApi.Add(num, &c.a) denum := c.baseApi.Add(&p.Y, &q.Y) // if p.y + q.y = 0, assign dummy 1 to denum and continue isYSumZero := c.baseApi.IsZero(denum) denum = c.baseApi.Select(isYSumZero, c.baseApi.One(), denum) - λ := c.baseApi.Div(num, denum) + lams, err := c.baseApi.NewHint(bjSlopeHint, 1, &p.X, &p.Y, &q.X, &q.Y, &c.a) + if err != nil { + panic(fmt.Sprintf("bj slope hint: %v", err)) + } + λ := lams[0] + c.baseApi.AssertEvalIsZero( + [][]*emulated.Element[B]{{λ, denum}, {pxplusqx, pxplusqx}, {&p.X, &q.X}, {&c.a}}, + []int{1, -1, 1, -1}, + ) // x = λ^2 - p.x - q.x xr := c.baseApi.MulMod(λ, λ) @@ -376,23 +395,10 @@ func (c *Curve[B, S]) double(p *AffinePoint[B]) *AffinePoint[B] { } func (c *Curve[B, S]) doubleGeneric(p *AffinePoint[B], unified bool) *AffinePoint[B] { - // compute λ = (3p.x²+a)/2*p.y, here we assume a=0 (j invariant 0 curve) - xx3a := c.baseApi.MulMod(&p.X, &p.X) - xx3a = c.baseApi.MulConst(xx3a, big.NewInt(3)) - if c.addA { - xx3a = c.baseApi.Add(xx3a, &c.a) - } - y2 := c.baseApi.MulConst(&p.Y, big.NewInt(2)) - var isDoubleYZero frontend.Variable = 0 - if unified { - // if 2*p.y = 0, assign dummy 1 to y2 and continue - isDoubleYZero = c.baseApi.IsZero(y2) - y2 = c.baseApi.Select(isDoubleYZero, c.baseApi.One(), y2) - } - λ := c.baseApi.Div(xx3a, y2) - if unified { - λ = c.baseApi.Select(isDoubleYZero, c.baseApi.Zero(), λ) - } + // compute λ = (3p.x²+a)/2*p.y via a hinted witness and a single deferred + // zero-assertion, without materializing x². In unified mode y ≡ 0 forces + // λ = 0 and keeps the assertion satisfiable (see tangentSlope). + λ := c.tangentSlope(p, unified) // xr = λ²-2p.x xr := c.baseApi.Eval([][]*emulated.Element[B]{{λ, λ}, {&p.X}}, []int{1, -2}) @@ -436,7 +442,7 @@ func (c *Curve[B, S]) tripleGeneric(p *AffinePoint[B], unified bool) *AffinePoin isDoubleYZero = c.baseApi.IsZero(y2) y2 = c.baseApi.Select(isDoubleYZero, c.baseApi.One(), y2) } - λ1 := c.baseApi.Div(xx, y2) + λ1 := c.assertedRatio(xx, y2) if unified { λ1 = c.baseApi.Select(isDoubleYZero, c.baseApi.Zero(), λ1) } @@ -452,7 +458,7 @@ func (c *Curve[B, S]) tripleGeneric(p *AffinePoint[B], unified bool) *AffinePoin isSecondSlopeDenominatorZero = c.baseApi.IsZero(x1x2) x1x2 = c.baseApi.Select(isSecondSlopeDenominatorZero, c.baseApi.One(), x1x2) } - λ2 := c.baseApi.Div(y2, x1x2) + λ2 := c.assertedRatio(y2, x1x2) if unified { λ2 = c.baseApi.Select(isSecondSlopeDenominatorZero, c.baseApi.Zero(), λ2) } @@ -496,7 +502,7 @@ func (c *Curve[B, S]) doubleAndAddGeneric(p, q *AffinePoint[B], unified bool) *A isChordDenominatorZero = c.baseApi.IsZero(xqxp) xqxp = c.baseApi.Select(isChordDenominatorZero, c.baseApi.One(), xqxp) } - λ1 := c.baseApi.Div(yqyp, xqxp) + λ1 := c.assertedRatio(yqyp, xqxp) if unified { λ1 = c.baseApi.Select(isChordDenominatorZero, c.baseApi.Zero(), λ1) } @@ -514,7 +520,7 @@ func (c *Curve[B, S]) doubleAndAddGeneric(p, q *AffinePoint[B], unified bool) *A isSecondSlopeDenominatorZero = c.baseApi.IsZero(x2xp) x2xp = c.baseApi.Select(isSecondSlopeDenominatorZero, c.baseApi.One(), x2xp) } - λ2 := c.baseApi.Div(ypyp, x2xp) + λ2 := c.assertedRatio(ypyp, x2xp) if unified { λ2 = c.baseApi.Select(isSecondSlopeDenominatorZero, c.baseApi.Zero(), λ2) } @@ -545,7 +551,7 @@ func (c *Curve[B, S]) doubleAndAddSelect(b frontend.Variable, p, q *AffinePoint[ // compute λ1 = (q.y-p.y)/(q.x-p.x) yqyp := c.baseApi.Sub(&q.Y, &p.Y) xqxp := c.baseApi.Sub(&q.X, &p.X) - λ1 := c.baseApi.Div(yqyp, xqxp) + λ1 := c.assertedRatio(yqyp, xqxp) // compute x2 = λ1²-p.x-q.x x2 := c.baseApi.Eval([][]*emulated.Element[B]{{λ1, λ1}, {&p.X}, {&q.X}}, []int{1, -1, -1}) @@ -558,7 +564,7 @@ func (c *Curve[B, S]) doubleAndAddSelect(b frontend.Variable, p, q *AffinePoint[ // compute -λ2 = λ1+2*t.y/(x2-t.x) ypyp := c.baseApi.MulConst(&t.Y, big.NewInt(2)) x2xp := c.baseApi.Sub(x2, &t.X) - λ2 := c.baseApi.Div(ypyp, x2xp) + λ2 := c.assertedRatio(ypyp, x2xp) λ2 = c.baseApi.Add(λ1, λ2) // compute x3 = (-λ2)²-t.x-x2 diff --git a/std/algebra/emulated/sw_emulated/slopes.go b/std/algebra/emulated/sw_emulated/slopes.go new file mode 100644 index 0000000000..2a4a0421a7 --- /dev/null +++ b/std/algebra/emulated/sw_emulated/slopes.go @@ -0,0 +1,198 @@ +package sw_emulated + +import ( + "errors" + "fmt" + "math/big" + + "github.com/consensys/gnark/std/math/emulated" +) + +// This file implements the slope computations used by the point addition and +// doubling formulas with hinted witnesses certified by single deferred +// zero-assertions, instead of the [emulated.Field.Div] pattern which costs +// two deferred checks (a multiplication check and an equality check). +// +// The constraint content is unchanged: Div asserts λ·den = num (mod p), and +// so does the zero-assertion here. In particular the exceptional behavior is +// identical: when den ≡ 0 and num ≢ 0 the assertion is unsatisfiable, and +// when den ≡ num ≡ 0 the slope is unconstrained (callers exclude or handle +// both cases exactly as they did with Div). + +// assertedRatio returns λ = num/den certified by the single deferred +// assertion λ·den − num ≡ 0 (mod p). +func (c *Curve[B, S]) assertedRatio(num, den *emulated.Element[B]) *emulated.Element[B] { + lams, err := c.baseApi.NewHint(ratioHint, 1, num, den) + if err != nil { + panic(fmt.Sprintf("ratio hint: %v", err)) + } + lam := lams[0] + c.baseApi.AssertEvalIsZero( + [][]*emulated.Element[B]{{lam, den}, {num}}, + []int{1, -1}, + ) + return lam +} + +// tangentSlope returns the tangent slope λ = (3x² + a)/(2y) at p, certified +// by a single deferred zero-assertion without materializing x². +// +// When unified is set, y ≡ 0 forces λ = 0 while keeping the assertion +// satisfiable (matching the dummy-denominator Select pattern of the unified +// formulas): the certified relation becomes λ·(2y + z) − (1−z)·(3x² + a) ≡ 0 +// with z the y ≡ 0 indicator bit. Otherwise y ≡ 0 makes the circuit +// unsatisfiable, as with the previous Div-based tangent. +func (c *Curve[B, S]) tangentSlope(p *AffinePoint[B], unified bool) *emulated.Element[B] { + var lams []*emulated.Element[B] + var err error + if c.addA { + lams, err = c.baseApi.NewHint(tangentHintA, 1, &p.X, &p.Y, &c.a) + } else { + lams, err = c.baseApi.NewHint(tangentHint, 1, &p.X, &p.Y) + } + if err != nil { + panic(fmt.Sprintf("tangent hint: %v", err)) + } + lam := lams[0] + if !unified { + // λ·2y − 3x² − a ≡ 0 + terms := [][]*emulated.Element[B]{{lam, &p.Y}, {&p.X, &p.X}} + coefs := []int{2, -3} + if c.addA { + terms = append(terms, []*emulated.Element[B]{&c.a}) + coefs = append(coefs, -1) + } + c.baseApi.AssertEvalIsZero(terms, coefs) + return lam + } + // unified: z = 1 iff y ≡ 0. The multiplier 2y + z never vanishes: it is + // 2y ≠ 0 when y ≢ 0 and 1 when y ≡ 0, so λ is always pinned; in the + // latter case the right-hand side is zeroed by 1−z and λ = 0. + isYZero := c.baseApi.IsZero(&p.Y) + zEl := c.baseApi.FromBits(isYZero) + nzEl := c.baseApi.Sub(c.baseApi.One(), zEl) + terms := [][]*emulated.Element[B]{{lam, &p.Y}, {lam, zEl}, {nzEl, &p.X, &p.X}} + coefs := []int{2, 1, -3} + if c.addA { + terms = append(terms, []*emulated.Element[B]{nzEl, &c.a}) + coefs = append(coefs, -1) + } + c.baseApi.AssertEvalIsZero(terms, coefs) + return lam +} + +// ratioHint computes num/den modulo the emulated modulus, or 0 when +// den ≡ 0 (the zero-assertion then decides satisfiability). +func ratioHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 2 || len(out) != 1 { + return errors.New("expecting two inputs and one output") + } + den := new(big.Int).Mod(in[1], p) + if den.Sign() == 0 { + out[0].SetInt64(0) + return nil + } + den.ModInverse(den, p) + out[0].Mod(in[0], p) + out[0].Mul(out[0], den).Mod(out[0], p) + return nil + }) +} + +// tangentHint computes 3x²/(2y) modulo the emulated modulus (curves with +// a = 0), or 0 when y ≡ 0. +func tangentHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 2 || len(out) != 1 { + return errors.New("expecting two inputs and one output") + } + return tangentSlopeVal(p, in[0], in[1], nil, out[0]) + }) +} + +// tangentHintA computes (3x² + a)/(2y) modulo the emulated modulus, or 0 +// when y ≡ 0. +func tangentHintA(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 3 || len(out) != 1 { + return errors.New("expecting three inputs and one output") + } + return tangentSlopeVal(p, in[0], in[1], in[2], out[0]) + }) +} + +func tangentSlopeVal(p, x, y, a, out *big.Int) error { + den := new(big.Int).Lsh(y, 1) + den.Mod(den, p) + if den.Sign() == 0 { + out.SetInt64(0) + return nil + } + den.ModInverse(den, p) + num := new(big.Int).Mul(x, x) + num.Mod(num, p) + num.Mul(num, big.NewInt(3)) + if a != nil { + num.Add(num, a) + } + num.Mod(num, p) + out.Mul(num, den).Mod(out, p) + return nil +} + +// unifiedSlopeHint computes the slope of the j-invariant-0 unified addition: +// the tangent 3x1²/(2y1) when x1 ≡ x2, the chord (y2−y1)/(x2−x1) otherwise, +// and 0 when the selected denominator vanishes. Inputs: x1, y1, x2, y2. +func unifiedSlopeHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 4 || len(out) != 1 { + return errors.New("expecting four inputs and one output") + } + x1 := new(big.Int).Mod(in[0], p) + y1 := new(big.Int).Mod(in[1], p) + x2 := new(big.Int).Mod(in[2], p) + y2 := new(big.Int).Mod(in[3], p) + dx := new(big.Int).Sub(x2, x1) + dx.Mod(dx, p) + if dx.Sign() == 0 { + return tangentSlopeVal(p, x1, y1, nil, out[0]) + } + dx.ModInverse(dx, p) + out[0].Sub(y2, y1) + out[0].Mul(out[0], dx).Mod(out[0], p) + return nil + }) +} + +// bjSlopeHint computes the Brier-Joye unified slope +// ((x1+x2)² − x1·x2 + a)/(y1 + y2) used on j ≠ 0 curves. When y1 + y2 ≡ 0 it +// returns the numerator itself, matching the dummy-1-denominator semantics of +// the caller (the result is then discarded by a select). Inputs: x1, y1, x2, +// y2, a. +func bjSlopeHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 5 || len(out) != 1 { + return errors.New("expecting five inputs and one output") + } + x1 := new(big.Int).Mod(in[0], p) + y1 := new(big.Int).Mod(in[1], p) + x2 := new(big.Int).Mod(in[2], p) + y2 := new(big.Int).Mod(in[3], p) + num := new(big.Int).Add(x1, x2) + num.Mul(num, num) + tmp := new(big.Int).Mul(x1, x2) + num.Sub(num, tmp) + num.Add(num, in[4]) + num.Mod(num, p) + den := new(big.Int).Add(y1, y2) + den.Mod(den, p) + if den.Sign() == 0 { + out[0].Set(num) + return nil + } + den.ModInverse(den, p) + out[0].Mul(num, den).Mod(out[0], p) + return nil + }) +} From 1deb0d17fac8554de8057501dd5d2e560314981a Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 11:02:45 -0400 Subject: [PATCH 05/15] perf: apply optims to native --- std/algebra/native/sw_bls12377/fixedbase.go | 278 ++++++++++++++++++ .../sw_bls12377/fixedbase_count_test.go | 69 +++++ std/algebra/native/sw_bls12377/g1.go | 15 + std/algebra/native/sw_bls12377/hints.go | 1 + 4 files changed, 363 insertions(+) create mode 100644 std/algebra/native/sw_bls12377/fixedbase.go create mode 100644 std/algebra/native/sw_bls12377/fixedbase_count_test.go diff --git a/std/algebra/native/sw_bls12377/fixedbase.go b/std/algebra/native/sw_bls12377/fixedbase.go new file mode 100644 index 0000000000..d48b5b32fb --- /dev/null +++ b/std/algebra/native/sw_bls12377/fixedbase.go @@ -0,0 +1,278 @@ +package sw_bls12377 + +import ( + "errors" + "fmt" + "math/big" + "sync" + + bls12377 "github.com/consensys/gnark-crypto/ecc/bls12-377" + fr_bls "github.com/consensys/gnark-crypto/ecc/bls12-377/fr" + "github.com/consensys/gnark/frontend" +) + +// This file implements a fixed-base scalar multiplication for compile-time +// constant G1 points using a signed-digit comb method, following the same +// construction as the emulated sw_emulated comb (see the package comment +// there for the algorithm and soundness argument): +// +// - the scalar is recoded into the odd k' = s + 1 − b0 represented by n +// signed binary digits, witnessed as the bits of c = (k' + 2^n − 1)/2 and +// pinned by the exact native identity 2c + b0 = s + 2^n (no wrap-around: +// 2^{n+1} is far below the native modulus); +// - windows of w digits select from compile-time constant tables +// [d(j)·2^{w·t}]P. With constant tables the selection is a free affine +// combination of the one-hot flags, so only the flag products cost +// constraints; +// - the parity correction −(1−b0) is folded into the top window table; +// - all partial sums are odd non-zero multiples of P, so the chain uses the +// incomplete AddAssign (3 constraints) except for the final complete +// AddUnified addition(s). +// +// In the native setting a small window width is optimal: the incomplete +// addition costs only 3 constraints, so wide windows are dominated by the +// one-hot flag products. With w = 2 the flags cost 2 constraints per window. +const g1CombWindow = 2 + +type g1CombData struct { + w int + nw int + n int + nbUnified int + // windows[t][j] = [d(j)·2^{w·t} mod r]P with d(j) = 2j − 2^w + 1, + // coordinates as big.Int + windows [][][2]*big.Int + // topEven[j] = [d(j)·2^{w·(nw−1)} − 1 mod r]P + topEven [][2]*big.Int +} + +// g1CombCache caches tables per constant base point. Keyed by the hex +// coordinates; access is concurrent-safe as circuits may compile in parallel. +var g1CombCache sync.Map + +// g1CombDataFor returns the comb tables for the constant point (gx, gy). It +// errors when the point is not a finite G1 (prime-order subgroup) point, in +// which case callers fall back to the generic scalar multiplication. +func g1CombDataFor(gx, gy *big.Int) (*g1CombData, error) { + key := gx.Text(16) + "|" + gy.Text(16) + if v, ok := g1CombCache.Load(key); ok { + return v.(*g1CombData), nil + } + var P bls12377.G1Affine + P.X.SetBigInt(gx) + P.Y.SetBigInt(gy) + if P.IsInfinity() { + return nil, errors.New("base point is the point at infinity") + } + if !P.IsOnCurve() { + return nil, errors.New("base point is not on the curve") + } + if !P.IsInSubGroup() { + return nil, errors.New("base point is not in the prime-order subgroup") + } + w := g1CombWindow + r := fr_bls.Modulus() + rBits := r.BitLen() + nw := (rBits + w - 1) / w + n := w * nw + // adding window t (t ≥ 1) with incomplete formulas is safe iff + // 2^{w·(t+1)} ≤ r; the final addition is always complete. + nbUnified := 0 + for t := nw - 1; t >= 1; t-- { + if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + break + } + nbUnified++ + } + nbUnified = max(nbUnified, 1) + + toBig := func(a *bls12377.G1Affine) [2]*big.Int { + if a.IsInfinity() { + // cannot happen for r > 2^w (window digits are odd non-zero and + // small), but guard anyway + return [2]*big.Int{nil, nil} + } + return [2]*big.Int{a.X.BigInt(new(big.Int)), a.Y.BigInt(new(big.Int))} + } + + half := 1 << (w - 1) + windows := make([][][2]*big.Int, nw) + var Bt bls12377.G1Jac + Bt.FromAffine(&P) + for t := 0; t < nw; t++ { + if t > 0 { + for k := 0; k < w; k++ { + Bt.DoubleAssign() + } + } + var D bls12377.G1Jac + D.Set(&Bt).DoubleAssign() + odd := make([]bls12377.G1Affine, half) + acc := new(bls12377.G1Jac).Set(&Bt) + odd[0].FromJacobian(acc) + for m := 1; m < half; m++ { + acc.AddAssign(&D) + odd[m].FromJacobian(acc) + } + tab := make([][2]*big.Int, 1< 0 { + e = odd[(d-1)/2] + } else { + e.Neg(&odd[(-d-1)/2]) + } + c := toBig(&e) + if c[0] == nil { + return nil, errors.New("unexpected point at infinity in comb table") + } + tab[j] = c + } + windows[t] = tab + } + // parity-folded top window: topEven[j] = windows[nw−1][j] + (−P) + var negP bls12377.G1Jac + var negPAff bls12377.G1Affine + negPAff.Neg(&P) + negP.FromAffine(&negPAff) + topEven := make([][2]*big.Int, 1< d.n+2 { + res := g1CombScalarMul(api, d, s) + p.X = res.X + p.Y = res.Y + return p + } + } + } return p.scalarMulGLV(api, Q, s, opts...) } diff --git a/std/algebra/native/sw_bls12377/hints.go b/std/algebra/native/sw_bls12377/hints.go index 7e07a2961d..b269adb4c0 100644 --- a/std/algebra/native/sw_bls12377/hints.go +++ b/std/algebra/native/sw_bls12377/hints.go @@ -13,6 +13,7 @@ import ( func GetHints() []solver.Hint { return []solver.Hint{ decomposeScalarG1Simple, + g1CombRecodeHint, scalarMulGLVG1Hint, rationalReconstructExt, pairingCheckHint, From 012999654663c400ae601c7173eede648ca3531f Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 11:19:26 -0400 Subject: [PATCH 06/15] perf: port optims to G2 --- .../emulated/sw_bls12381/fixedbase_g2.go | 367 ++++++++++++++++++ .../emulated/sw_bls12381/fixedbase_g2_test.go | 73 ++++ std/algebra/emulated/sw_bls12381/g2.go | 8 + std/algebra/emulated/sw_bls12381/hints.go | 1 + std/algebra/emulated/sw_bn254/fixedbase_g2.go | 367 ++++++++++++++++++ .../emulated/sw_bn254/fixedbase_g2_test.go | 73 ++++ std/algebra/emulated/sw_bn254/g2.go | 8 + std/algebra/emulated/sw_bn254/hints.go | 1 + .../native/sw_bls12377/fixedbase_g2.go | 214 ++++++++++ .../native/sw_bls12377/fixedbase_g2_test.go | 71 ++++ std/algebra/native/sw_bls12377/g2.go | 35 +- 11 files changed, 1216 insertions(+), 2 deletions(-) create mode 100644 std/algebra/emulated/sw_bls12381/fixedbase_g2.go create mode 100644 std/algebra/emulated/sw_bls12381/fixedbase_g2_test.go create mode 100644 std/algebra/emulated/sw_bn254/fixedbase_g2.go create mode 100644 std/algebra/emulated/sw_bn254/fixedbase_g2_test.go create mode 100644 std/algebra/native/sw_bls12377/fixedbase_g2.go create mode 100644 std/algebra/native/sw_bls12377/fixedbase_g2_test.go diff --git a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go new file mode 100644 index 0000000000..e091a0adf2 --- /dev/null +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go @@ -0,0 +1,367 @@ +package sw_bls12381 + +import ( + "errors" + "fmt" + "math/big" + "sync" + + bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381" + fr_bls "github.com/consensys/gnark-crypto/ecc/bls12-381/fr" + "github.com/consensys/gnark/frontend" + "github.com/consensys/gnark/internal/smallfields" + "github.com/consensys/gnark/std/algebra/emulated/fields_bls12381" + limbs "github.com/consensys/gnark/std/internal/limbcomposition" + "github.com/consensys/gnark/std/math/emulated" +) + +// Fixed-base signed-digit comb for compile-time constant G2 points, following +// the same construction and soundness argument as the G1 comb in +// sw_emulated/fixedbase.go: odd recode k' = s + 1 − b0 pinned by a single +// mod-r assertion, constant window tables [d(j)·2^{w·t}]P selected by one-hot +// flags shared across all coordinate limbs, provably collision-free +// incomplete chord additions over E2, parity fold in the top window, and a +// complete AddUnified tail. +const g2CombWindow = 8 + +type g2CombData struct { + w int + nw int + n int + nbUnified int + // windows[t][j] = [d(j)·2^{w·t} mod r]P as (X.A0, X.A1, Y.A0, Y.A1) + windows [][][4]*big.Int + topEven [][4]*big.Int +} + +var g2CombCache sync.Map + +func g2CombDataFor(x0, x1, y0, y1 *big.Int) (*g2CombData, error) { + key := x0.Text(16) + "|" + x1.Text(16) + "|" + y0.Text(16) + "|" + y1.Text(16) + if v, ok := g2CombCache.Load(key); ok { + return v.(*g2CombData), nil + } + var P bls12381.G2Affine + P.X.A0.SetBigInt(x0) + P.X.A1.SetBigInt(x1) + P.Y.A0.SetBigInt(y0) + P.Y.A1.SetBigInt(y1) + if P.IsInfinity() { + return nil, errors.New("base point is the point at infinity") + } + if !P.IsOnCurve() { + return nil, errors.New("base point is not on the twist") + } + if !P.IsInSubGroup() { + return nil, errors.New("base point is not in the prime-order subgroup") + } + w := g2CombWindow + r := fr_bls.Modulus() + rBits := r.BitLen() + nw := (rBits + w - 1) / w + n := w * nw + var frParams ScalarField + if n > int(frParams.NbLimbs()*frParams.BitsPerLimb()) { + return nil, errors.New("recoded scalar exceeds scalar field emulation capacity") + } + nbUnified := 0 + for t := nw - 1; t >= 1; t-- { + if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + break + } + nbUnified++ + } + nbUnified = max(nbUnified, 1) + + toBig := func(a *bls12381.G2Affine) ([4]*big.Int, error) { + if a.IsInfinity() { + return [4]*big.Int{}, errors.New("unexpected point at infinity in comb table") + } + return [4]*big.Int{ + a.X.A0.BigInt(new(big.Int)), + a.X.A1.BigInt(new(big.Int)), + a.Y.A0.BigInt(new(big.Int)), + a.Y.A1.BigInt(new(big.Int)), + }, nil + } + + half := 1 << (w - 1) + windows := make([][][4]*big.Int, nw) + var Bt bls12381.G2Jac + Bt.FromAffine(&P) + for t := 0; t < nw; t++ { + if t > 0 { + for k := 0; k < w; k++ { + Bt.DoubleAssign() + } + } + var D bls12381.G2Jac + D.Set(&Bt).DoubleAssign() + odd := make([]bls12381.G2Affine, half) + acc := new(bls12381.G2Jac).Set(&Bt) + odd[0].FromJacobian(acc) + for m := 1; m < half; m++ { + acc.AddAssign(&D) + odd[m].FromJacobian(acc) + } + tab := make([][4]*big.Int, 1< 0 { + e = odd[(d-1)/2] + } else { + e.Neg(&odd[(-d-1)/2]) + } + c, err := toBig(&e) + if err != nil { + return nil, err + } + tab[j] = c + } + windows[t] = tab + } + var negP bls12381.G2Jac + var negPAff bls12381.G2Affine + negPAff.Neg(&P) + negP.FromAffine(&negPAff) + topEven := make([][4]*big.Int, 1< int(frParams.NbLimbs()*frParams.BitsPerLimb()) { + return nil, errors.New("recoded scalar exceeds scalar field emulation capacity") + } + nbUnified := 0 + for t := nw - 1; t >= 1; t-- { + if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + break + } + nbUnified++ + } + nbUnified = max(nbUnified, 1) + + toBig := func(a *bn254.G2Affine) ([4]*big.Int, error) { + if a.IsInfinity() { + return [4]*big.Int{}, errors.New("unexpected point at infinity in comb table") + } + return [4]*big.Int{ + a.X.A0.BigInt(new(big.Int)), + a.X.A1.BigInt(new(big.Int)), + a.Y.A0.BigInt(new(big.Int)), + a.Y.A1.BigInt(new(big.Int)), + }, nil + } + + half := 1 << (w - 1) + windows := make([][][4]*big.Int, nw) + var Bt bn254.G2Jac + Bt.FromAffine(&P) + for t := 0; t < nw; t++ { + if t > 0 { + for k := 0; k < w; k++ { + Bt.DoubleAssign() + } + } + var D bn254.G2Jac + D.Set(&Bt).DoubleAssign() + odd := make([]bn254.G2Affine, half) + acc := new(bn254.G2Jac).Set(&Bt) + odd[0].FromJacobian(acc) + for m := 1; m < half; m++ { + acc.AddAssign(&D) + odd[m].FromJacobian(acc) + } + tab := make([][4]*big.Int, 1< 0 { + e = odd[(d-1)/2] + } else { + e.Neg(&odd[(-d-1)/2]) + } + c, err := toBig(&e) + if err != nil { + return nil, err + } + tab[j] = c + } + windows[t] = tab + } + var negP bn254.G2Jac + var negPAff bn254.G2Affine + negPAff.Neg(&P) + negP.FromAffine(&negPAff) + topEven := make([][4]*big.Int, 1<= 1; t-- { + if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + break + } + nbUnified++ + } + nbUnified = max(nbUnified, 1) + + toBig := func(a *bls12377.G2Affine) ([4]*big.Int, error) { + if a.IsInfinity() { + return [4]*big.Int{}, errors.New("unexpected point at infinity in comb table") + } + return [4]*big.Int{ + a.X.A0.BigInt(new(big.Int)), + a.X.A1.BigInt(new(big.Int)), + a.Y.A0.BigInt(new(big.Int)), + a.Y.A1.BigInt(new(big.Int)), + }, nil + } + + half := 1 << (w - 1) + windows := make([][][4]*big.Int, nw) + var Bt bls12377.G2Jac + Bt.FromAffine(&P) + for t := 0; t < nw; t++ { + if t > 0 { + for k := 0; k < w; k++ { + Bt.DoubleAssign() + } + } + var D bls12377.G2Jac + D.Set(&Bt).DoubleAssign() + odd := make([]bls12377.G2Affine, half) + acc := new(bls12377.G2Jac).Set(&Bt) + odd[0].FromJacobian(acc) + for m := 1; m < half; m++ { + acc.AddAssign(&D) + odd[m].FromJacobian(acc) + } + tab := make([][4]*big.Int, 1< 0 { + e = odd[(d-1)/2] + } else { + e.Neg(&odd[(-d-1)/2]) + } + c, err := toBig(&e) + if err != nil { + return nil, err + } + tab[j] = c + } + windows[t] = tab + } + var negP bls12377.G2Jac + var negPAff bls12377.G2Affine + negPAff.Neg(&P) + negP.FromAffine(&negPAff) + topEven := make([][4]*big.Int, 1< d.n+2 { + res := g2CombScalarMul(api, d, s) + p.X = res.X + p.Y = res.Y + return p + } + } + } + } + } + return p.varScalarMul(api, Q, s, opts...) } // varScalarMul sets P = [s]Q and returns P. It doesn't modify Q nor s. @@ -503,6 +519,21 @@ func (p *g2AffP) DoubleAndAdd(api frontend.API, p1, p2 *g2AffP) *g2AffP { // does not support complete arithmetic and will produce incorrect results for // s=0. func (p *g2AffP) ScalarMulBase(api frontend.API, s frontend.Variable) *g2AffP { + // use the fixed-base comb on the generator when supported (complete + // arithmetic, handles the zero scalar); fall back to the legacy + // precomputed-table double-and-add otherwise. + _, _, _, g2gen := bls12377.Generators() + if d, err := g2CombDataFor( + g2gen.X.A0.BigInt(new(big.Int)), + g2gen.X.A1.BigInt(new(big.Int)), + g2gen.Y.A0.BigInt(new(big.Int)), + g2gen.Y.A1.BigInt(new(big.Int)), + ); err == nil && api.Compiler().FieldBitLen() > d.n+2 { + res := g2CombScalarMul(api, d, s) + p.X = res.X + p.Y = res.Y + return p + } points := getTwistPoints() From 43679124de9a3ffad6b767fde4afd8208c4ef4be Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 11:34:39 -0400 Subject: [PATCH 07/15] perf: carry y acc unmaterialized --- std/algebra/emulated/sw_emulated/hints.go | 2 + std/algebra/emulated/sw_emulated/point.go | 27 +++- std/algebra/emulated/sw_emulated/slopes.go | 137 +++++++++++++++++++++ 3 files changed, 165 insertions(+), 1 deletion(-) diff --git a/std/algebra/emulated/sw_emulated/hints.go b/std/algebra/emulated/sw_emulated/hints.go index b15568e87c..9526608b82 100644 --- a/std/algebra/emulated/sw_emulated/hints.go +++ b/std/algebra/emulated/sw_emulated/hints.go @@ -41,6 +41,8 @@ func GetHints() []solver.Hint { tangentHintA, unifiedSlopeHint, bjSlopeHint, + implicitTangentHint, + implicitChordHint, } } diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 7f34929445..972b2d849d 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -1200,6 +1200,12 @@ func (c *Curve[B, S]) jointScalarMulGLVUnsafe(Q, R *AffinePoint[B], s, t *emulat // hence have the same X coordinates. var Bi *AffinePoint[B] + // run the loop with the accumulator y-coordinate in implicit form on + // j = 0 curves (see implicitDoubleAndAdd). + var iAcc *implicitAcc[B] + if !c.addA { + iAcc = c.implicitFromAffine(Acc) + } for i := nbits - 1; i > 0; i-- { // selectorY takes values in [0,15] selectorY := c.api.Add( @@ -1226,7 +1232,14 @@ func (c *Curve[B, S]) jointScalarMulGLVUnsafe(Q, R *AffinePoint[B], s, t *emulat ), } // Acc = [2]Acc + Bi - Acc = c.doubleAndAdd(Acc, Bi) + if iAcc != nil { + c.implicitDoubleAndAdd(iAcc, Bi) + } else { + Acc = c.doubleAndAdd(Acc, Bi) + } + } + if iAcc != nil { + Acc = c.implicitToAffine(iAcc) } // i = 0 @@ -1899,6 +1912,13 @@ func (c *Curve[B, S]) scalarMulGLVAndFakeGLV(P *AffinePoint[B], s *emulated.Elem // hence have the same X coordinates. var Bi *AffinePoint[B] + // in incomplete mode on j = 0 curves, run the loop with the accumulator + // y-coordinate in implicit form (see implicitDoubleAndAdd): one deferred + // check saved per iteration, y materialized once after the loop. + var iAcc *implicitAcc[B] + if cfg.IncompleteArithmetic && !c.addA { + iAcc = c.implicitFromAffine(Acc) + } for i := nbits - 1; i > 0; i-- { // selectorY takes values in [0,15] selectorY := c.api.Add( @@ -1928,10 +1948,15 @@ func (c *Curve[B, S]) scalarMulGLVAndFakeGLV(P *AffinePoint[B], s *emulated.Elem if !cfg.IncompleteArithmetic { Acc = c.doubleGeneric(Acc, true) Acc = addFn(Acc, Bi) + } else if iAcc != nil { + c.implicitDoubleAndAdd(iAcc, Bi) } else { Acc = c.doubleAndAdd(Acc, Bi) } } + if iAcc != nil { + Acc = c.implicitToAffine(iAcc) + } // i = 0 // subtract the P, Q, Φ(P), Φ(Q) if the first bits are 0 diff --git a/std/algebra/emulated/sw_emulated/slopes.go b/std/algebra/emulated/sw_emulated/slopes.go index 2a4a0421a7..dbeafd3893 100644 --- a/std/algebra/emulated/sw_emulated/slopes.go +++ b/std/algebra/emulated/sw_emulated/slopes.go @@ -196,3 +196,140 @@ func bjSlopeHint(_ *big.Int, inputs, outputs []*big.Int) error { return nil }) } + +// implicitAcc carries the accumulator of a double-and-add chain with its +// y-coordinate in implicit form +// +// y = lam·(xT − x) − yT +// +// where (xT, yT) is the last added point and lam the slope of that addition. +// The expression is degree-1 in materialized values, so consecutive chain +// steps can consume it inside their deferred zero-assertions instead of +// materializing y at every iteration (one Eval saved per iteration). The +// encoding resets at every addition — the implicit form always references +// only the most recent slope — so expressions do not grow with chain length. +type implicitAcc[B emulated.FieldParams] struct { + x *emulated.Element[B] + lam *emulated.Element[B] + xT *emulated.Element[B] + yT *emulated.Element[B] +} + +// implicitFromAffine encodes a materialized point: with xT = x the lam term +// vanishes identically for any lam, and yT = −y recovers y. +func (c *Curve[B, S]) implicitFromAffine(p *AffinePoint[B]) *implicitAcc[B] { + return &implicitAcc[B]{x: &p.X, lam: &p.X, xT: &p.X, yT: c.baseApi.Neg(&p.Y)} +} + +// implicitToAffine materializes the accumulator y-coordinate. +func (c *Curve[B, S]) implicitToAffine(acc *implicitAcc[B]) *AffinePoint[B] { + y := c.baseApi.Eval( + [][]*emulated.Element[B]{{acc.lam, c.baseApi.Sub(acc.xT, acc.x)}, {acc.yT}}, + []int{1, -1}, + ) + return &AffinePoint[B]{X: *acc.x, Y: *y} +} + +// implicitDoubleAndAdd sets acc = 2·acc + q using incomplete formulas with +// the accumulator y kept implicit throughout: the doubling tangent and the +// addition chord are each certified by a single deferred zero-assertion which +// consumes the implicit y, and only the two x-coordinates are materialized. +// +// ⚠️ Incomplete: requires j-invariant 0 (a = 0), acc not 2-torsion and +// q ≠ ±[2]acc — the same exceptional envelope as the ELM-based doubleAndAdd +// it replaces in the incomplete scalar-multiplication loops, where the +// accumulator anchoring excludes these cases for honest witnesses. +func (c *Curve[B, S]) implicitDoubleAndAdd(acc *implicitAcc[B], q *AffinePoint[B]) { + // tangent: λd·2y = 3x² with y = lam·(xT − x) − yT: + // 2·λd·lam·(xT − x) − 2·λd·yT − 3·x² ≡ 0 + lamDs, err := c.baseApi.NewHint(implicitTangentHint, 1, acc.lam, acc.x, acc.xT, acc.yT) + if err != nil { + panic(fmt.Sprintf("implicit tangent hint: %v", err)) + } + lamD := lamDs[0] + dxPrev := c.baseApi.Sub(acc.xT, acc.x) + c.baseApi.AssertEvalIsZero( + [][]*emulated.Element[B]{{lamD, acc.lam, dxPrev}, {lamD, acc.yT}, {acc.x, acc.x}}, + []int{2, -2, -3}, + ) + // xd = λd² − 2x + xd := c.baseApi.Eval([][]*emulated.Element[B]{{lamD, lamD}, {acc.x}}, []int{1, -2}) + + // chord through q: λa·(xq − xd) = yq − yd with + // yd = λd·(x − xd) − y = λd·(x − xd) − lam·(xT − x) + yT: + // λa·(xq − xd) + λd·(x − xd) − lam·(xT − x) + yT − yq ≡ 0 + lamAs, err := c.baseApi.NewHint(implicitChordHint, 1, acc.lam, acc.x, acc.xT, acc.yT, lamD, xd, &q.X, &q.Y) + if err != nil { + panic(fmt.Sprintf("implicit chord hint: %v", err)) + } + lamA := lamAs[0] + dxq := c.baseApi.Sub(&q.X, xd) + dxd := c.baseApi.Sub(acc.x, xd) + c.baseApi.AssertEvalIsZero( + [][]*emulated.Element[B]{{lamA, dxq}, {lamD, dxd}, {acc.lam, dxPrev}, {acc.yT}, {&q.Y}}, + []int{1, 1, -1, 1, -1}, + ) + // xa = λa² − xd − xq + xa := c.baseApi.Eval([][]*emulated.Element[B]{{lamA, lamA}, {xd}, {&q.X}}, []int{1, -1, -1}) + + // the new implicit y references only this addition: y = λa·(xq − xa) − yq + acc.x = xa + acc.lam = lamA + acc.xT = &q.X + acc.yT = &q.Y +} + +// implicitYVal recomputes the implicit accumulator y = lam·(xT − x) − yT. +func implicitYVal(p, lam, x, xT, yT *big.Int) *big.Int { + y := new(big.Int).Sub(xT, x) + y.Mul(y, lam).Sub(y, yT).Mod(y, p) + return y +} + +// implicitTangentHint computes the doubling slope 3x²/(2y) (a = 0 curves) +// with y in implicit form. Inputs: lam, x, xT, yT. Returns 0 when y ≡ 0. +func implicitTangentHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 4 || len(out) != 1 { + return errors.New("expecting four inputs and one output") + } + lam := new(big.Int).Mod(in[0], p) + x := new(big.Int).Mod(in[1], p) + xT := new(big.Int).Mod(in[2], p) + yT := new(big.Int).Mod(in[3], p) + y := implicitYVal(p, lam, x, xT, yT) + return tangentSlopeVal(p, x, y, nil, out[0]) + }) +} + +// implicitChordHint computes the second slope of the implicit double-and-add: +// λa = (yq − yd)/(xq − xd) with yd = λd·(x − xd) − y and y in implicit form. +// Inputs: lam, x, xT, yT, λd, xd, xq, yq. Returns 0 when xq ≡ xd. +func implicitChordHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 8 || len(out) != 1 { + return errors.New("expecting eight inputs and one output") + } + lam := new(big.Int).Mod(in[0], p) + x := new(big.Int).Mod(in[1], p) + xT := new(big.Int).Mod(in[2], p) + yT := new(big.Int).Mod(in[3], p) + lamD := new(big.Int).Mod(in[4], p) + xd := new(big.Int).Mod(in[5], p) + xq := new(big.Int).Mod(in[6], p) + yq := new(big.Int).Mod(in[7], p) + y := implicitYVal(p, lam, x, xT, yT) + yd := new(big.Int).Sub(x, xd) + yd.Mul(yd, lamD).Sub(yd, y).Mod(yd, p) + den := new(big.Int).Sub(xq, xd) + den.Mod(den, p) + if den.Sign() == 0 { + out[0].SetInt64(0) + return nil + } + den.ModInverse(den, p) + out[0].Sub(yq, yd) + out[0].Mul(out[0], den).Mod(out[0], p) + return nil + }) +} From 3affdf314f05d4ee099db632a665e85ebe7b564c Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 11:52:05 -0400 Subject: [PATCH 08/15] perf: E2 slope fusion --- .../emulated/sw_bls12381/fixedbase_g2.go | 136 ++++++++++++++++-- std/algebra/emulated/sw_bls12381/hints.go | 1 + std/algebra/emulated/sw_bn254/fixedbase_g2.go | 136 ++++++++++++++++-- std/algebra/emulated/sw_bn254/hints.go | 1 + 4 files changed, 252 insertions(+), 22 deletions(-) diff --git a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go index e091a0adf2..2103730b05 100644 --- a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go @@ -257,13 +257,53 @@ func (g2 *G2) g2CombSelect(table [][4]*big.Int, bs []frontend.Variable) g2AffP { } } -// g2CombAdd is the incomplete chord addition over E2 (complete on the comb -// chain by the collision-freeness argument). -func (g2 *G2) g2CombAdd(p, q *g2AffP) *g2AffP { - lam := g2.Ext2.DivUnchecked(g2.Ext2.Sub(&q.Y, &p.Y), g2.Ext2.Sub(&q.X, &p.X)) - xr := g2.Ext2.Sub(g2.Ext2.Square(lam), g2.Ext2.Add(&p.X, &q.X)) - yr := g2.Ext2.Sub(g2.Ext2.Mul(lam, g2.Ext2.Sub(&p.X, xr)), &p.Y) - return &g2AffP{X: *xr, Y: *yr} +// g2CombChainStep performs one incomplete chord addition of the comb chain +// with the accumulator y-coordinate kept implicit (y = λprev·(xTprev − x) − +// yTprev over E2, as in the G1 comb): the slope λ ∈ Fp² is witnessed by a +// hint and certified by two deferred zero-assertions (one per E2 component, +// with the u² = −1 cross terms), and only the two x-components are +// materialized. This costs 4 deferred checks per addition instead of ~9 for +// the DivUnchecked/Square/Mul formulation. +func (g2 *G2) g2CombChainStep(lamPrev, xAcc, xTPrev, yTPrev *fields_bls12381.E2, q *g2AffP) (lam, xNew *fields_bls12381.E2) { + lams, err := g2.fp.NewHint(g2CombChainLambdaHint, 2, + &lamPrev.A0, &lamPrev.A1, &xAcc.A0, &xAcc.A1, + &xTPrev.A0, &xTPrev.A1, &yTPrev.A0, &yTPrev.A1, + &q.X.A0, &q.X.A1, &q.Y.A0, &q.Y.A1) + if err != nil { + panic(fmt.Sprintf("comb chain hint: %v", err)) + } + lam = &fields_bls12381.E2{A0: *lams[0], A1: *lams[1]} + dxC := g2.Ext2.Sub(&q.X, xAcc) + dxP := g2.Ext2.Sub(xTPrev, xAcc) + // slope identity λ·(xT − x) + λprev·(xTprev − x) − yT − yTprev = 0 (E2), + // component 0 (u² = −1): + g2.fp.AssertEvalIsZero( + [][]*emulated.Element[BaseField]{ + {&lam.A0, &dxC.A0}, {&lam.A1, &dxC.A1}, + {&lamPrev.A0, &dxP.A0}, {&lamPrev.A1, &dxP.A1}, + {&q.Y.A0}, {&yTPrev.A0}, + }, + []int{1, -1, 1, -1, -1, -1}, + ) + // component 1: + g2.fp.AssertEvalIsZero( + [][]*emulated.Element[BaseField]{ + {&lam.A0, &dxC.A1}, {&lam.A1, &dxC.A0}, + {&lamPrev.A0, &dxP.A1}, {&lamPrev.A1, &dxP.A0}, + {&q.Y.A1}, {&yTPrev.A1}, + }, + []int{1, 1, 1, 1, -1, -1}, + ) + // x' = λ² − x − xT componentwise: (λ0² − λ1², 2λ0λ1) − ... + x0 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lam.A0, &lam.A0}, {&lam.A1, &lam.A1}, {&xAcc.A0}, {&q.X.A0}}, + []int{1, -1, -1, -1}, + ) + x1 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lam.A0, &lam.A1}, {&xAcc.A1}, {&q.X.A1}}, + []int{2, -1, -1}, + ) + return lam, &fields_bls12381.E2{A0: *x0, A1: *x1} } // scalarMulComb computes [s]P for the constant base point of the tables d. @@ -298,11 +338,34 @@ func (g2 *G2) scalarMulComb(d *g2CombData, s *Scalar) *G2Affine { pts[nw-1] = g2.g2CombSelect(stacked, topBits) nbInc := nw - 1 - d.nbUnified - acc := &pts[0] - for t := 1; t <= nbInc; t++ { - acc = g2.g2CombAdd(acc, &pts[t]) + var res *G2Affine + if nbInc >= 1 { + // implicit-y chain: at the first step the accumulator is T_0 whose + // y is directly available; encode it as λprev·(xTprev − x) − yTprev + // with xTprev = x = xT0 and yTprev = −yT0, so the λprev term + // vanishes identically for any λprev (we pass xT0 as a dummy). + xAcc := &pts[0].X + lamPrev := &pts[0].X + xTPrev := &pts[0].X + yTPrev := g2.Ext2.Neg(&pts[0].Y) + for t := 1; t <= nbInc; t++ { + lamPrev, xAcc = g2.g2CombChainStep(lamPrev, xAcc, xTPrev, yTPrev, &pts[t]) + xTPrev, yTPrev = &pts[t].X, &pts[t].Y + } + // materialize y = λprev·(xTprev − x) − yTprev once + dy := g2.Ext2.Sub(xTPrev, xAcc) + y0 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lamPrev.A0, &dy.A0}, {&lamPrev.A1, &dy.A1}, {&yTPrev.A0}}, + []int{1, -1, -1}, + ) + y1 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lamPrev.A0, &dy.A1}, {&lamPrev.A1, &dy.A0}, {&yTPrev.A1}}, + []int{1, 1, -1}, + ) + res = &G2Affine{P: g2AffP{X: *xAcc, Y: fields_bls12381.E2{A0: *y0, A1: *y1}}} + } else { + res = &G2Affine{P: pts[0]} } - res := &G2Affine{P: *acc} for t := nbInc + 1; t <= nw-1; t++ { res = g2.AddUnified(res, &G2Affine{P: pts[t]}) } @@ -365,3 +428,54 @@ func g2CombRecodeHint(_ *big.Int, inputs, outputs []*big.Int) error { return nil }) } + +// g2CombChainLambdaHint computes the chord slope of the next comb chain +// addition over Fp² (u² = −1). Inputs (component pairs): λprev, x, xTprev, +// yTprev, xT, yT. The accumulator y is recomputed in its implicit form +// y = λprev·(xTprev − x) − yTprev and the output is λ = (yT − y)/(xT − x). +func g2CombChainLambdaHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 12 || len(out) != 2 { + return errors.New("expecting twelve inputs and two outputs") + } + mod := func(v *big.Int) *big.Int { return new(big.Int).Mod(v, p) } + lamP := [2]*big.Int{mod(in[0]), mod(in[1])} + x := [2]*big.Int{mod(in[2]), mod(in[3])} + xTp := [2]*big.Int{mod(in[4]), mod(in[5])} + yTp := [2]*big.Int{mod(in[6]), mod(in[7])} + xT := [2]*big.Int{mod(in[8]), mod(in[9])} + yT := [2]*big.Int{mod(in[10]), mod(in[11])} + e2Sub := func(a, b [2]*big.Int) [2]*big.Int { + return [2]*big.Int{ + new(big.Int).Mod(new(big.Int).Sub(a[0], b[0]), p), + new(big.Int).Mod(new(big.Int).Sub(a[1], b[1]), p), + } + } + e2Mul := func(a, b [2]*big.Int) [2]*big.Int { + c0 := new(big.Int).Mul(a[0], b[0]) + c0.Sub(c0, new(big.Int).Mul(a[1], b[1])).Mod(c0, p) + c1 := new(big.Int).Mul(a[0], b[1]) + c1.Add(c1, new(big.Int).Mul(a[1], b[0])).Mod(c1, p) + return [2]*big.Int{c0, c1} + } + // y = λprev·(xTprev − x) − yTprev + y := e2Mul(lamP, e2Sub(xTp, x)) + y = e2Sub(y, yTp) + den := e2Sub(xT, x) + // den⁻¹ = (den0 − den1·u)/(den0² + den1²) + nrm := new(big.Int).Mul(den[0], den[0]) + nrm.Add(nrm, new(big.Int).Mul(den[1], den[1])).Mod(nrm, p) + if nrm.Sign() == 0 { + return errors.New("comb chain: x-coordinate collision") + } + nrm.ModInverse(nrm, p) + inv := [2]*big.Int{ + new(big.Int).Mod(new(big.Int).Mul(den[0], nrm), p), + new(big.Int).Mod(new(big.Int).Mul(new(big.Int).Neg(den[1]), nrm), p), + } + lam := e2Mul(e2Sub(yT, y), inv) + out[0].Set(lam[0]) + out[1].Set(lam[1]) + return nil + }) +} diff --git a/std/algebra/emulated/sw_bls12381/hints.go b/std/algebra/emulated/sw_bls12381/hints.go index eb036198ba..9e2abb6d9d 100644 --- a/std/algebra/emulated/sw_bls12381/hints.go +++ b/std/algebra/emulated/sw_bls12381/hints.go @@ -25,6 +25,7 @@ func GetHints() []solver.Hint { millerLoopAndCheckFinalExpHint, scalarMulG2Hint, g2CombRecodeHint, + g2CombChainLambdaHint, rationalReconstructExtG2, g1SqrtRatioHint, g2SqrtRatioHint, diff --git a/std/algebra/emulated/sw_bn254/fixedbase_g2.go b/std/algebra/emulated/sw_bn254/fixedbase_g2.go index 3853de5df2..034f031ebe 100644 --- a/std/algebra/emulated/sw_bn254/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bn254/fixedbase_g2.go @@ -257,13 +257,53 @@ func (g2 *G2) g2CombSelect(table [][4]*big.Int, bs []frontend.Variable) g2AffP { } } -// g2CombAdd is the incomplete chord addition over E2 (complete on the comb -// chain by the collision-freeness argument). -func (g2 *G2) g2CombAdd(p, q *g2AffP) *g2AffP { - lam := g2.Ext2.DivUnchecked(g2.Ext2.Sub(&q.Y, &p.Y), g2.Ext2.Sub(&q.X, &p.X)) - xr := g2.Ext2.Sub(g2.Ext2.Square(lam), g2.Ext2.Add(&p.X, &q.X)) - yr := g2.Ext2.Sub(g2.Ext2.Mul(lam, g2.Ext2.Sub(&p.X, xr)), &p.Y) - return &g2AffP{X: *xr, Y: *yr} +// g2CombChainStep performs one incomplete chord addition of the comb chain +// with the accumulator y-coordinate kept implicit (y = λprev·(xTprev − x) − +// yTprev over E2, as in the G1 comb): the slope λ ∈ Fp² is witnessed by a +// hint and certified by two deferred zero-assertions (one per E2 component, +// with the u² = −1 cross terms), and only the two x-components are +// materialized. This costs 4 deferred checks per addition instead of ~9 for +// the DivUnchecked/Square/Mul formulation. +func (g2 *G2) g2CombChainStep(lamPrev, xAcc, xTPrev, yTPrev *fields_bn254.E2, q *g2AffP) (lam, xNew *fields_bn254.E2) { + lams, err := g2.fp.NewHint(g2CombChainLambdaHint, 2, + &lamPrev.A0, &lamPrev.A1, &xAcc.A0, &xAcc.A1, + &xTPrev.A0, &xTPrev.A1, &yTPrev.A0, &yTPrev.A1, + &q.X.A0, &q.X.A1, &q.Y.A0, &q.Y.A1) + if err != nil { + panic(fmt.Sprintf("comb chain hint: %v", err)) + } + lam = &fields_bn254.E2{A0: *lams[0], A1: *lams[1]} + dxC := g2.Ext2.Sub(&q.X, xAcc) + dxP := g2.Ext2.Sub(xTPrev, xAcc) + // slope identity λ·(xT − x) + λprev·(xTprev − x) − yT − yTprev = 0 (E2), + // component 0 (u² = −1): + g2.fp.AssertEvalIsZero( + [][]*emulated.Element[BaseField]{ + {&lam.A0, &dxC.A0}, {&lam.A1, &dxC.A1}, + {&lamPrev.A0, &dxP.A0}, {&lamPrev.A1, &dxP.A1}, + {&q.Y.A0}, {&yTPrev.A0}, + }, + []int{1, -1, 1, -1, -1, -1}, + ) + // component 1: + g2.fp.AssertEvalIsZero( + [][]*emulated.Element[BaseField]{ + {&lam.A0, &dxC.A1}, {&lam.A1, &dxC.A0}, + {&lamPrev.A0, &dxP.A1}, {&lamPrev.A1, &dxP.A0}, + {&q.Y.A1}, {&yTPrev.A1}, + }, + []int{1, 1, 1, 1, -1, -1}, + ) + // x' = λ² − x − xT componentwise: (λ0² − λ1², 2λ0λ1) − ... + x0 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lam.A0, &lam.A0}, {&lam.A1, &lam.A1}, {&xAcc.A0}, {&q.X.A0}}, + []int{1, -1, -1, -1}, + ) + x1 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lam.A0, &lam.A1}, {&xAcc.A1}, {&q.X.A1}}, + []int{2, -1, -1}, + ) + return lam, &fields_bn254.E2{A0: *x0, A1: *x1} } // scalarMulComb computes [s]P for the constant base point of the tables d. @@ -298,11 +338,34 @@ func (g2 *G2) scalarMulComb(d *g2CombData, s *Scalar) *G2Affine { pts[nw-1] = g2.g2CombSelect(stacked, topBits) nbInc := nw - 1 - d.nbUnified - acc := &pts[0] - for t := 1; t <= nbInc; t++ { - acc = g2.g2CombAdd(acc, &pts[t]) + var res *G2Affine + if nbInc >= 1 { + // implicit-y chain: at the first step the accumulator is T_0 whose + // y is directly available; encode it as λprev·(xTprev − x) − yTprev + // with xTprev = x = xT0 and yTprev = −yT0, so the λprev term + // vanishes identically for any λprev (we pass xT0 as a dummy). + xAcc := &pts[0].X + lamPrev := &pts[0].X + xTPrev := &pts[0].X + yTPrev := g2.Ext2.Neg(&pts[0].Y) + for t := 1; t <= nbInc; t++ { + lamPrev, xAcc = g2.g2CombChainStep(lamPrev, xAcc, xTPrev, yTPrev, &pts[t]) + xTPrev, yTPrev = &pts[t].X, &pts[t].Y + } + // materialize y = λprev·(xTprev − x) − yTprev once + dy := g2.Ext2.Sub(xTPrev, xAcc) + y0 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lamPrev.A0, &dy.A0}, {&lamPrev.A1, &dy.A1}, {&yTPrev.A0}}, + []int{1, -1, -1}, + ) + y1 := g2.fp.Eval( + [][]*emulated.Element[BaseField]{{&lamPrev.A0, &dy.A1}, {&lamPrev.A1, &dy.A0}, {&yTPrev.A1}}, + []int{1, 1, -1}, + ) + res = &G2Affine{P: g2AffP{X: *xAcc, Y: fields_bn254.E2{A0: *y0, A1: *y1}}} + } else { + res = &G2Affine{P: pts[0]} } - res := &G2Affine{P: *acc} for t := nbInc + 1; t <= nw-1; t++ { res = g2.AddUnified(res, &G2Affine{P: pts[t]}) } @@ -365,3 +428,54 @@ func g2CombRecodeHint(_ *big.Int, inputs, outputs []*big.Int) error { return nil }) } + +// g2CombChainLambdaHint computes the chord slope of the next comb chain +// addition over Fp² (u² = −1). Inputs (component pairs): λprev, x, xTprev, +// yTprev, xT, yT. The accumulator y is recomputed in its implicit form +// y = λprev·(xTprev − x) − yTprev and the output is λ = (yT − y)/(xT − x). +func g2CombChainLambdaHint(_ *big.Int, inputs, outputs []*big.Int) error { + return emulated.UnwrapHint(inputs, outputs, func(p *big.Int, in, out []*big.Int) error { + if len(in) != 12 || len(out) != 2 { + return errors.New("expecting twelve inputs and two outputs") + } + mod := func(v *big.Int) *big.Int { return new(big.Int).Mod(v, p) } + lamP := [2]*big.Int{mod(in[0]), mod(in[1])} + x := [2]*big.Int{mod(in[2]), mod(in[3])} + xTp := [2]*big.Int{mod(in[4]), mod(in[5])} + yTp := [2]*big.Int{mod(in[6]), mod(in[7])} + xT := [2]*big.Int{mod(in[8]), mod(in[9])} + yT := [2]*big.Int{mod(in[10]), mod(in[11])} + e2Sub := func(a, b [2]*big.Int) [2]*big.Int { + return [2]*big.Int{ + new(big.Int).Mod(new(big.Int).Sub(a[0], b[0]), p), + new(big.Int).Mod(new(big.Int).Sub(a[1], b[1]), p), + } + } + e2Mul := func(a, b [2]*big.Int) [2]*big.Int { + c0 := new(big.Int).Mul(a[0], b[0]) + c0.Sub(c0, new(big.Int).Mul(a[1], b[1])).Mod(c0, p) + c1 := new(big.Int).Mul(a[0], b[1]) + c1.Add(c1, new(big.Int).Mul(a[1], b[0])).Mod(c1, p) + return [2]*big.Int{c0, c1} + } + // y = λprev·(xTprev − x) − yTprev + y := e2Mul(lamP, e2Sub(xTp, x)) + y = e2Sub(y, yTp) + den := e2Sub(xT, x) + // den⁻¹ = (den0 − den1·u)/(den0² + den1²) + nrm := new(big.Int).Mul(den[0], den[0]) + nrm.Add(nrm, new(big.Int).Mul(den[1], den[1])).Mod(nrm, p) + if nrm.Sign() == 0 { + return errors.New("comb chain: x-coordinate collision") + } + nrm.ModInverse(nrm, p) + inv := [2]*big.Int{ + new(big.Int).Mod(new(big.Int).Mul(den[0], nrm), p), + new(big.Int).Mod(new(big.Int).Mul(new(big.Int).Neg(den[1]), nrm), p), + } + lam := e2Mul(e2Sub(yT, y), inv) + out[0].Set(lam[0]) + out[1].Set(lam[1]) + return nil + }) +} diff --git a/std/algebra/emulated/sw_bn254/hints.go b/std/algebra/emulated/sw_bn254/hints.go index b906cff12d..25f557a07a 100644 --- a/std/algebra/emulated/sw_bn254/hints.go +++ b/std/algebra/emulated/sw_bn254/hints.go @@ -23,6 +23,7 @@ func GetHints() []solver.Hint { millerLoopAndCheckFinalExpHint, scalarMulG2Hint, g2CombRecodeHint, + g2CombChainLambdaHint, rationalReconstructExtG2, } } From e8d89f7be6b9f2acce2620ce3bee28380fe91c8d Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 12:00:42 -0400 Subject: [PATCH 09/15] perf: non-fold MSM to use scalarMulComb --- .../emulated/sw_emulated/fixedbase_test.go | 88 +++++++++++++++++++ std/algebra/emulated/sw_emulated/point.go | 51 +++++++++-- 2 files changed, 131 insertions(+), 8 deletions(-) diff --git a/std/algebra/emulated/sw_emulated/fixedbase_test.go b/std/algebra/emulated/sw_emulated/fixedbase_test.go index 6e95baf5b9..dab3973dcf 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase_test.go +++ b/std/algebra/emulated/sw_emulated/fixedbase_test.go @@ -351,3 +351,91 @@ func TestScalarMulConstPointBLS12381(t *testing.T) { } assert.True(found, "expected to find a non-subgroup curve point") } + +type msmMixedTest[T, S emulated.FieldParams] struct { + P AffinePoint[T] // variable point + S [4]emulated.Element[S] + Q AffinePoint[T] + useOld bool +} + +func (c *msmMixedTest[T, S]) Define(api frontend.API) error { + cr, err := New[T, S](api, GetCurveParams[T]()) + if err != nil { + return err + } + // two constant points (G and 2G), two variable (P twice to keep the + // witness small) + g := cr.Generator() + g2 := AffinePoint[T]{ + X: *cr.baseApi.NewElement(cr.params.Gm[0][0]), + Y: *cr.baseApi.NewElement(cr.params.Gm[0][1]), + } + pts := []*AffinePoint[T]{g, &g2, &c.P, &c.P} + scs := []*emulated.Element[S]{&c.S[0], &c.S[1], &c.S[2], &c.S[3]} + res, err := cr.MultiScalarMul(pts, scs) + if err != nil { + return err + } + cr.AssertIsEqual(res, &c.Q) + return nil +} + +// TestMSMConstRouting checks correctness of the constant-term routing in +// MultiScalarMul against gnark-crypto. +func TestMSMConstRouting(t *testing.T) { + assert := test.NewAssert(t) + _, g := secp256k1.Generators() + var g2, P secp256k1.G1Affine + // params.Gm[0] is [3]G + g2.ScalarMultiplication(&g, big.NewInt(3)) + var rp fr_secp.Element + _, _ = rp.SetRandom() + P.ScalarMultiplication(&g, rp.BigInt(new(big.Int))) + var S [4]*big.Int + var expected secp256k1.G1Jac + pts := []secp256k1.G1Affine{g, g2, P, P} + for i := range S { + var rs fr_secp.Element + _, _ = rs.SetRandom() + S[i] = rs.BigInt(new(big.Int)) + var t secp256k1.G1Jac + var ta secp256k1.G1Affine + ta.ScalarMultiplication(&pts[i], S[i]) + t.FromAffine(&ta) + if i == 0 { + expected = t + } else { + expected.AddAssign(&t) + } + } + var E secp256k1.G1Affine + E.FromJacobian(&expected) + circuit := msmMixedTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{} + witness := msmMixedTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{ + P: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](P.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](P.Y), + }, + Q: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](E.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](E.Y), + }, + } + for i := range S { + witness.S[i] = emulated.ValueOf[emulated.Secp256k1Fr](S[i]) + } + err := test.IsSolved(&circuit, &witness, testCurve.ScalarField()) + assert.NoError(err) +} + +func TestMSMConstRoutingCount(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + circuit := msmMixedTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{} + ccs, err := frontend.Compile(testCurve.ScalarField(), r1cs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("MSM 4 terms (2 const + 2 var) r1cs =", ccs.GetNbConstraints()) +} diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 972b2d849d..488b973cc3 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -1383,16 +1383,51 @@ func (c *Curve[B, S]) MultiScalarMul(p []*AffinePoint[B], s []*emulated.Element[ if len(p) != len(s) { return nil, fmt.Errorf("mismatching points and scalars slice lengths") } - n := len(p) + // route compile-time constant points of prime order through the + // fixed-base comb (complete arithmetic, much cheaper) and fold only + // the remaining variable points through the joint scalar + // multiplications. var res *AffinePoint[B] - if n%2 == 1 { - res = c.ScalarMul(p[n-1], s[n-1], opts...) - } else { - res = c.jointScalarMul(p[n-2], p[n-1], s[n-2], s[n-1], opts...) + varP := make([]*AffinePoint[B], 0, len(p)) + varS := make([]*emulated.Element[S], 0, len(s)) + for i := range p { + var d *combData + if px, ok := c.baseApi.ConstantValue(&p[i].X); ok { + if py, ok := c.baseApi.ConstantValue(&p[i].Y); ok { + if dd, derr := c.combDataFor(px, py, combDefaultWindow); derr == nil { + d = dd + } + } + } + if d == nil { + varP = append(varP, p[i]) + varS = append(varS, s[i]) + continue + } + q := c.scalarMulComb(d, s[i]) + if res == nil { + res = q + } else { + res = addFn(res, q) + } } - for i := 1; i < n-1; i += 2 { - q := c.jointScalarMul(p[i-1], p[i], s[i-1], s[i], opts...) - res = addFn(res, q) + n := len(varP) + if n > 0 { + var vres *AffinePoint[B] + if n%2 == 1 { + vres = c.ScalarMul(varP[n-1], varS[n-1], opts...) + } else { + vres = c.jointScalarMul(varP[n-2], varP[n-1], varS[n-2], varS[n-1], opts...) + } + for i := 1; i < n-1; i += 2 { + q := c.jointScalarMul(varP[i-1], varP[i], varS[i-1], varS[i], opts...) + vres = addFn(vres, q) + } + if res == nil { + res = vres + } else { + res = addFn(res, vres) + } } return res, nil } else { From 7e56e19ccaf3b64bc4fd8b4baf364b94a464ff98 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 12:13:16 -0400 Subject: [PATCH 10/15] perf: tune windows in comb scs --- .../emulated/sw_bls12381/fixedbase_g2.go | 20 ++++++++--- std/algebra/emulated/sw_bn254/fixedbase_g2.go | 20 ++++++++--- std/algebra/emulated/sw_emulated/fixedbase.go | 16 +++++++++ .../emulated/sw_emulated/fixedbase_test.go | 33 +++++++++++++++++++ std/algebra/emulated/sw_emulated/point.go | 12 +++---- 5 files changed, 87 insertions(+), 14 deletions(-) diff --git a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go index 2103730b05..aa1d28b784 100644 --- a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go @@ -24,6 +24,19 @@ import ( // complete AddUnified tail. const g2CombWindow = 8 +// g2CombPlonkWindow is the window width on PLONKish backends, where the +// one-hot selection's constant linear combinations expand into addition +// gates (see the G1 comb in sw_emulated). +const g2CombPlonkWindow = 4 + +// g2CombWindowFor returns the comb window width for the current backend. +func (g2 *G2) g2CombWindowFor() int { + if _, ok := g2.api.Compiler().(frontend.PlonkAPI); ok { + return g2CombPlonkWindow + } + return g2CombWindow +} + type g2CombData struct { w int nw int @@ -36,8 +49,8 @@ type g2CombData struct { var g2CombCache sync.Map -func g2CombDataFor(x0, x1, y0, y1 *big.Int) (*g2CombData, error) { - key := x0.Text(16) + "|" + x1.Text(16) + "|" + y0.Text(16) + "|" + y1.Text(16) +func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { + key := fmt.Sprintf("%d|%s|%s|%s|%s", w, x0.Text(16), x1.Text(16), y0.Text(16), y1.Text(16)) if v, ok := g2CombCache.Load(key); ok { return v.(*g2CombData), nil } @@ -55,7 +68,6 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int) (*g2CombData, error) { if !P.IsInSubGroup() { return nil, errors.New("base point is not in the prime-order subgroup") } - w := g2CombWindow r := fr_bls.Modulus() rBits := r.BitLen() nw := (rBits + w - 1) / w @@ -395,7 +407,7 @@ func (g2 *G2) g2CombTryConst(Q *G2Affine) (*g2CombData, bool) { if !ok { return nil, false } - d, err := g2CombDataFor(x0, x1, y0, y1) + d, err := g2CombDataFor(x0, x1, y0, y1, g2.g2CombWindowFor()) if err != nil { return nil, false } diff --git a/std/algebra/emulated/sw_bn254/fixedbase_g2.go b/std/algebra/emulated/sw_bn254/fixedbase_g2.go index 034f031ebe..74f4f09f86 100644 --- a/std/algebra/emulated/sw_bn254/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bn254/fixedbase_g2.go @@ -24,6 +24,19 @@ import ( // complete AddUnified tail. const g2CombWindow = 8 +// g2CombPlonkWindow is the window width on PLONKish backends, where the +// one-hot selection's constant linear combinations expand into addition +// gates (see the G1 comb in sw_emulated). +const g2CombPlonkWindow = 4 + +// g2CombWindowFor returns the comb window width for the current backend. +func (g2 *G2) g2CombWindowFor() int { + if _, ok := g2.api.Compiler().(frontend.PlonkAPI); ok { + return g2CombPlonkWindow + } + return g2CombWindow +} + type g2CombData struct { w int nw int @@ -36,8 +49,8 @@ type g2CombData struct { var g2CombCache sync.Map -func g2CombDataFor(x0, x1, y0, y1 *big.Int) (*g2CombData, error) { - key := x0.Text(16) + "|" + x1.Text(16) + "|" + y0.Text(16) + "|" + y1.Text(16) +func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { + key := fmt.Sprintf("%d|%s|%s|%s|%s", w, x0.Text(16), x1.Text(16), y0.Text(16), y1.Text(16)) if v, ok := g2CombCache.Load(key); ok { return v.(*g2CombData), nil } @@ -55,7 +68,6 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int) (*g2CombData, error) { if !P.IsInSubGroup() { return nil, errors.New("base point is not in the prime-order subgroup") } - w := g2CombWindow r := fr_bn.Modulus() rBits := r.BitLen() nw := (rBits + w - 1) / w @@ -395,7 +407,7 @@ func (g2 *G2) g2CombTryConst(Q *G2Affine) (*g2CombData, bool) { if !ok { return nil, false } - d, err := g2CombDataFor(x0, x1, y0, y1) + d, err := g2CombDataFor(x0, x1, y0, y1, g2.g2CombWindowFor()) if err != nil { return nil, false } diff --git a/std/algebra/emulated/sw_emulated/fixedbase.go b/std/algebra/emulated/sw_emulated/fixedbase.go index 6874b99018..906dc952e1 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase.go +++ b/std/algebra/emulated/sw_emulated/fixedbase.go @@ -53,6 +53,22 @@ import ( // constraint-count optimum in R1CS. const combDefaultWindow = 8 +// combPlonkWindow is the window width used on PLONKish backends. The one-hot +// selection's wide constant linear combinations are free in R1CS but expand +// into one addition gate per term in PLONK, making the selector cost scale +// with 2^w·nbLimbs per window; a smaller window rebalances selector versus +// chain-addition cost (measured on secp256k1: w=4 is ~35% cheaper than w=8 in +// scs, while w=8 remains ~30% cheaper than w=4 in R1CS). +const combPlonkWindow = 4 + +// combWindow returns the comb window width for the current backend. +func (c *Curve[B, S]) combWindow() int { + if _, ok := c.api.Compiler().(frontend.PlonkAPI); ok { + return combPlonkWindow + } + return combDefaultWindow +} + // combData holds the compile-time data of the comb: the constant window // tables and the derived parameters. type combData struct { diff --git a/std/algebra/emulated/sw_emulated/fixedbase_test.go b/std/algebra/emulated/sw_emulated/fixedbase_test.go index dab3973dcf..65749a317c 100644 --- a/std/algebra/emulated/sw_emulated/fixedbase_test.go +++ b/std/algebra/emulated/sw_emulated/fixedbase_test.go @@ -439,3 +439,36 @@ func TestMSMConstRoutingCount(t *testing.T) { assert.NoError(err) t.Log("MSM 4 terms (2 const + 2 var) r1cs =", ccs.GetNbConstraints()) } + +// TestScalarMulBaseCombPlonkSelector validates the comb on an actual PLONKish +// (scs) compilation, end-to-end through witness solving. +func TestScalarMulBaseCombPlonkSelector(t *testing.T) { + assert := test.NewAssert(t) + _, g := secp256k1.Generators() + r := fr_secp.Modulus() + randFn := func() *big.Int { + var rnd fr_secp.Element + _, _ = rnd.SetRandom() + return rnd.BigInt(new(big.Int)) + } + // ScalarMulBaseTest goes through the public dispatch, which picks the + // PLONK window on scs + circuit := ScalarMulBaseTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{} + ccs, err := frontend.Compile(testCurve.ScalarField(), scs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("ScalarMulBase scs constraints =", ccs.GetNbConstraints()) + for _, s := range []*big.Int{big.NewInt(0), big.NewInt(1), new(big.Int).Sub(r, big.NewInt(1)), randFn()} { + var S secp256k1.G1Affine + S.ScalarMultiplication(&g, s) + witness := ScalarMulBaseTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{ + S: emulated.ValueOf[emulated.Secp256k1Fr](s), + Q: AffinePoint[emulated.Secp256k1Fp]{ + X: emulated.ValueOf[emulated.Secp256k1Fp](S.X), + Y: emulated.ValueOf[emulated.Secp256k1Fp](S.Y), + }, + } + w, err := frontend.NewWitness(&witness, testCurve.ScalarField()) + assert.NoError(err) + assert.NoError(ccs.IsSolved(w), "s=%s", s.String()) + } +} diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go index 488b973cc3..b1dfbcee7e 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -667,7 +667,7 @@ func (c *Curve[B, S]) ScalarMul(p *AffinePoint[B], s *emulated.Element[S], opts // at compile time that (px, py) is a finite curve point of prime // order r; otherwise we fall back to the variable-base methods // below which have no such requirement. - if d, err := c.combDataFor(px, py, combDefaultWindow); err == nil { + if d, err := c.combDataFor(px, py, c.combWindow()); err == nil { return c.scalarMulComb(d, s) } } @@ -1287,8 +1287,8 @@ func (c *Curve[B, S]) jointScalarMulGLVUnsafe(Q, R *AffinePoint[B], s, t *emulat // - curves with an efficient endomorphism inherit the documented exceptional // set of [Curve.scalarMulGLVAndFakeGLV]. func (c *Curve[B, S]) ScalarMulBase(s *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] { - if _, err := c.combData(combDefaultWindow); err == nil { - return c.scalarMulBaseComb(s, combDefaultWindow) + if _, err := c.combData(c.combWindow()); err == nil { + return c.scalarMulBaseComb(s, c.combWindow()) } if c.eigenvalue != nil && c.thirdRootOne != nil { return c.scalarMulGLVAndFakeGLV(c.Generator(), s, opts...) @@ -1330,7 +1330,7 @@ func (c *Curve[B, S]) JointScalarMulBase(p *AffinePoint[B], s2, s1 *emulated.Ele if err != nil { panic(fmt.Sprintf("parse opts: %v", err)) } - if _, cerr := c.combData(combDefaultWindow); cerr == nil && !cfg.IncompleteArithmetic { + if _, cerr := c.combData(c.combWindow()); cerr == nil && !cfg.IncompleteArithmetic { // In complete mode, compute the fixed-base part with the comb method // (complete, handles s1 = 0) and the variable-base part separately, // and merge with the complete addition. This is cheaper than two @@ -1343,7 +1343,7 @@ func (c *Curve[B, S]) JointScalarMulBase(p *AffinePoint[B], s2, s1 *emulated.Ele // to the generator (e.g. p = [2]g fails for a noticeable fraction of // scalars), whereas the joint Shamir-based algorithm below handles // those points. - sm1 := c.scalarMulBaseComb(s1, combDefaultWindow) + sm1 := c.scalarMulBaseComb(s1, c.combWindow()) sm2 := c.ScalarMul(p, s2, opts...) return c.AddUnified(sm1, sm2) } @@ -1394,7 +1394,7 @@ func (c *Curve[B, S]) MultiScalarMul(p []*AffinePoint[B], s []*emulated.Element[ var d *combData if px, ok := c.baseApi.ConstantValue(&p[i].X); ok { if py, ok := c.baseApi.ConstantValue(&p[i].Y); ok { - if dd, derr := c.combDataFor(px, py, combDefaultWindow); derr == nil { + if dd, derr := c.combDataFor(px, py, c.combWindow()); derr == nil { d = dd } } From 23adea923e537875fbb52b88f7f6ab7cb2da7f0c Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 24 Jul 2026 13:34:12 -0400 Subject: [PATCH 11/15] test: up stats --- internal/stats/latest_stats.csv | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/internal/stats/latest_stats.csv b/internal/stats/latest_stats.csv index 9c94a4fc85..4fdd5a3eec 100644 --- a/internal/stats/latest_stats.csv +++ b/internal/stats/latest_stats.csv @@ -111,24 +111,24 @@ pairing_bn254,bn254,groth16,505959,823836 pairing_bn254,bn254,plonk,1646544,1572896 pairing_bw6761,bn254,groth16,1589471,2646707 pairing_bw6761,bn254,plonk,5318762,5097941 -scalar_mul_G1_bn254,bn254,groth16,107499,163170 -scalar_mul_G1_bn254,bn254,plonk,353793,338853 -scalar_mul_G1_bn254_incomplete,bn254,groth16,50892,81121 -scalar_mul_G1_bn254_incomplete,bn254,plonk,184266,177694 +scalar_mul_G1_bn254,bn254,groth16,91587,137256 +scalar_mul_G1_bn254,bn254,plonk,296080,282901 +scalar_mul_G1_bn254_incomplete,bn254,groth16,40447,63364 +scalar_mul_G1_bn254_incomplete,bn254,plonk,143517,138048 scalar_mul_G2_bls12381,bn254,groth16,302753,456022 scalar_mul_G2_bls12381,bn254,plonk,988292,945775 scalar_mul_G2_bn254,bn254,groth16,216495,325402 scalar_mul_G2_bn254,bn254,plonk,719650,688306 scalar_mul_G2_bw6761,bn254,groth16,387972,616049 scalar_mul_G2_bw6761,bn254,plonk,1241833,1189810 -scalar_mul_P256,bn254,groth16,96434,151466 -scalar_mul_P256,bn254,plonk,328264,315107 -scalar_mul_P256_incomplete,bn254,groth16,75252,121496 -scalar_mul_P256_incomplete,bn254,plonk,262529,252901 -scalar_mul_secp256k1,bn254,groth16,107536,163231 -scalar_mul_secp256k1,bn254,plonk,353942,338998 -scalar_mul_secp256k1_incomplete,bn254,groth16,50932,81189 -scalar_mul_secp256k1_incomplete,bn254,plonk,184432,177853 +scalar_mul_P256,bn254,groth16,78648,122092 +scalar_mul_P256,bn254,plonk,261626,250726 +scalar_mul_P256_incomplete,bn254,groth16,60232,95719 +scalar_mul_P256_incomplete,bn254,plonk,209383,201119 +scalar_mul_secp256k1,bn254,groth16,91623,137316 +scalar_mul_secp256k1,bn254,plonk,296229,283046 +scalar_mul_secp256k1_incomplete,bn254,groth16,40487,63432 +scalar_mul_secp256k1_incomplete,bn254,plonk,143683,138207 selector/binaryMux_4,bn254,groth16,5,3 selector/binaryMux_4,bls12_377,groth16,5,3 selector/binaryMux_4,bls12_381,groth16,5,3 From 655a640322d7e4d18463144f4257aa582e22d6ed Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Thu, 30 Jul 2026 11:25:00 -0400 Subject: [PATCH 12/15] perf: partial top window 5 in fixed-base scalarmul --- .../emulated/sw_bls12381/fixedbase_g2.go | 38 +++++++---- .../emulated/sw_bls12381/fixedbase_g2_test.go | 4 ++ std/algebra/emulated/sw_bn254/fixedbase_g2.go | 38 +++++++---- .../emulated/sw_bn254/fixedbase_g2_test.go | 4 ++ std/algebra/emulated/sw_emulated/fixedbase.go | 66 +++++++++++-------- .../emulated/sw_emulated/fixedbase_test.go | 4 +- 6 files changed, 102 insertions(+), 52 deletions(-) diff --git a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go index aa1d28b784..e023c4a4e0 100644 --- a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go @@ -26,8 +26,9 @@ const g2CombWindow = 8 // g2CombPlonkWindow is the window width on PLONKish backends, where the // one-hot selection's constant linear combinations expand into addition -// gates (see the G1 comb in sw_emulated). -const g2CombPlonkWindow = 4 +// gates (see the G1 comb in sw_emulated). Width 5 is the current SCS optimum +// measured on BN254 and BLS12-381 G2. +const g2CombPlonkWindow = 5 // g2CombWindowFor returns the comb window width for the current backend. func (g2 *G2) g2CombWindowFor() int { @@ -41,8 +42,11 @@ type g2CombData struct { w int nw int n int + tw int nbUnified int - // windows[t][j] = [d(j)·2^{w·t} mod r]P as (X.A0, X.A1, Y.A0, Y.A1) + // windows[t][j] = [d(j)·2^{w·t} mod r]P as (X.A0, X.A1, Y.A0, Y.A1). + // For the lower windows d(j) = 2j − 2^w + 1; for the top window w is + // replaced by tw. windows [][][4]*big.Int topEven [][4]*big.Int } @@ -71,14 +75,19 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { r := fr_bls.Modulus() rBits := r.BitLen() nw := (rBits + w - 1) / w - n := w * nw + n := rBits var frParams ScalarField if n > int(frParams.NbLimbs()*frParams.BitsPerLimb()) { return nil, errors.New("recoded scalar exceeds scalar field emulation capacity") } + tw := n - w*(nw-1) nbUnified := 0 for t := nw - 1; t >= 1; t-- { - if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + endBit := w * (t + 1) + if t == nw-1 { + endBit = n + } + if new(big.Int).Lsh(big.NewInt(1), uint(endBit)).Cmp(r) <= 0 { break } nbUnified++ @@ -97,7 +106,6 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { }, nil } - half := 1 << (w - 1) windows := make([][][4]*big.Int, nw) var Bt bls12381.G2Jac Bt.FromAffine(&P) @@ -107,8 +115,13 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { Bt.DoubleAssign() } } + tw := w + if t == nw-1 { + tw = n - w*(nw-1) + } var D bls12381.G2Jac D.Set(&Bt).DoubleAssign() + half := 1 << (tw - 1) odd := make([]bls12381.G2Affine, half) acc := new(bls12381.G2Jac).Set(&Bt) odd[0].FromJacobian(acc) @@ -116,9 +129,9 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { acc.AddAssign(&D) odd[m].FromJacobian(acc) } - tab := make([][4]*big.Int, 1< 0 { e = odd[(d-1)/2] @@ -137,7 +150,7 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { var negPAff bls12381.G2Affine negPAff.Neg(&P) negP.FromAffine(&negPAff) - topEven := make([][4]*big.Int, 1<= 1; t-- { - if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + endBit := w * (t + 1) + if t == nw-1 { + endBit = n + } + if new(big.Int).Lsh(big.NewInt(1), uint(endBit)).Cmp(r) <= 0 { break } nbUnified++ @@ -97,7 +106,6 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { }, nil } - half := 1 << (w - 1) windows := make([][][4]*big.Int, nw) var Bt bn254.G2Jac Bt.FromAffine(&P) @@ -107,8 +115,13 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { Bt.DoubleAssign() } } + tw := w + if t == nw-1 { + tw = n - w*(nw-1) + } var D bn254.G2Jac D.Set(&Bt).DoubleAssign() + half := 1 << (tw - 1) odd := make([]bn254.G2Affine, half) acc := new(bn254.G2Jac).Set(&Bt) odd[0].FromJacobian(acc) @@ -116,9 +129,9 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { acc.AddAssign(&D) odd[m].FromJacobian(acc) } - tab := make([][4]*big.Int, 1< 0 { e = odd[(d-1)/2] @@ -137,7 +150,7 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { var negPAff bn254.G2Affine negPAff.Neg(&P) negP.FromAffine(&negPAff) - topEven := make([][4]*big.Int, 1< scalarCap { return nil, fmt.Errorf("recoded scalar needs %d bits, scalar field emulation has capacity %d", n, scalarCap) } - // adding window t (t ≥ 1) with incomplete formulas is safe iff - // 2^{w·(t+1)} ≤ r; the final addition is always complete as it may cancel - // to the point at infinity or double. + tw := n - w*(nw-1) + // Adding window t (t ≥ 1) with incomplete formulas is safe iff the absolute + // value of every running sum is below r. For a partial top window this + // bound is 2^n on the top window and 2^{w·(t+1)} otherwise. The final + // addition is always complete as it may cancel to the point at infinity or + // double. nbUnified := 0 for t := nw - 1; t >= 1; t-- { - if new(big.Int).Lsh(big.NewInt(1), uint(w*(t+1))).Cmp(r) <= 0 { + endBit := w * (t + 1) + if t == nw-1 { + endBit = n + } + if new(big.Int).Lsh(big.NewInt(1), uint(endBit)).Cmp(r) <= 0 { break } nbUnified++ @@ -218,7 +227,6 @@ func computeCombData(gx, gy, a, b, prime, r *big.Int, w int, scalarCap int) (*co nbUnified = max(nbUnified, 1) G := &combAffine{x: new(big.Int).Set(gx), y: new(big.Int).Set(gy)} - half := 1 << (w - 1) windows := make([][][2]*big.Int, nw) Bt := G var err error @@ -230,11 +238,16 @@ func computeCombData(gx, gy, a, b, prime, r *big.Int, w int, scalarCap int) (*co } } } + tw := w + if t == nw-1 { + tw = n - w*(nw-1) + } // odd multiples odd[m] = [(2m+1)·2^{w·t}]G D, err := combDouble(Bt, a, prime) if err != nil { return nil, err } + half := 1 << (tw - 1) odd := make([]*combAffine, half) odd[0] = Bt for m := 1; m < half; m++ { @@ -242,9 +255,9 @@ func computeCombData(gx, gy, a, b, prime, r *big.Int, w int, scalarCap int) (*co return nil, err } } - tab := make([][2]*big.Int, 1< 0 { pt = odd[(d-1)/2] @@ -257,7 +270,7 @@ func computeCombData(gx, gy, a, b, prime, r *big.Int, w int, scalarCap int) (*co } // parity-folded top window: topEven[j] = windows[nw−1][j] + (−G) negG := combNeg(G, prime) - topEven := make([][2]*big.Int, 1< Date: Fri, 31 Jul 2026 11:26:08 -0400 Subject: [PATCH 13/15] perf: specilize the last add in fixed-base scalarmul --- .../emulated/sw_bls12381/fixedbase_g2.go | 88 ++++++++++++++++--- std/algebra/emulated/sw_bls12381/g2_test.go | 59 +++++++++++++ std/algebra/emulated/sw_bn254/fixedbase_g2.go | 88 ++++++++++++++++--- std/algebra/emulated/sw_bn254/g2_test.go | 61 +++++++++++++ 4 files changed, 276 insertions(+), 20 deletions(-) diff --git a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go index e023c4a4e0..d5cebdf438 100644 --- a/std/algebra/emulated/sw_bls12381/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go @@ -48,7 +48,11 @@ type g2CombData struct { // For the lower windows d(j) = 2j − 2^w + 1; for the top window w is // replaced by tw. windows [][][4]*big.Int - topEven [][4]*big.Int + // doubles[t][j] = 2*windows[t][j], used by the screened complete tail + // when a trailing addition degenerates into a doubling. + doubles [][][4]*big.Int + topEven [][4]*big.Int + topEvenDoubles [][4]*big.Int } var g2CombCache sync.Map @@ -105,8 +109,17 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { a.Y.A1.BigInt(new(big.Int)), }, nil } + toDoubleBig := func(a *bls12381.G2Affine) ([4]*big.Int, error) { + var j bls12381.G2Jac + j.FromAffine(a) + j.DoubleAssign() + var d bls12381.G2Affine + d.FromJacobian(&j) + return toBig(&d) + } windows := make([][][4]*big.Int, nw) + doubles := make([][][4]*big.Int, nw) var Bt bls12381.G2Jac Bt.FromAffine(&P) for t := 0; t < nw; t++ { @@ -130,6 +143,7 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { odd[m].FromJacobian(acc) } tab := make([][4]*big.Int, 1< nbInc { + dbls[t] = g2.g2CombSelect(d.doubles[t], cbits[t*w:(t+1)*w]) + } } stacked := make([][4]*big.Int, 0, 2< nbInc { + dbls[nw-1] = g2.g2CombSelect(stackedDoubles, topBits) + } - nbInc := nw - 1 - d.nbUnified var res *G2Affine if nbInc >= 1 { // implicit-y chain: at the first step the accumulator is T_0 whose @@ -393,7 +461,7 @@ func (g2 *G2) scalarMulComb(d *g2CombData, s *Scalar) *G2Affine { res = &G2Affine{P: pts[0]} } for t := nbInc + 1; t <= nw-1; t++ { - res = g2.AddUnified(res, &G2Affine{P: pts[t]}) + res = g2.g2CombTailAdd(res, &pts[t], &dbls[t]) } return res } diff --git a/std/algebra/emulated/sw_bls12381/g2_test.go b/std/algebra/emulated/sw_bls12381/g2_test.go index 25ec197e9e..9b6412319c 100644 --- a/std/algebra/emulated/sw_bls12381/g2_test.go +++ b/std/algebra/emulated/sw_bls12381/g2_test.go @@ -101,6 +101,65 @@ func TestScalarMulG2EdgeCases(t *testing.T) { } } +type scalarMulConstG2Circuit struct { + S Scalar + Res G2Affine + + px0, px1 *big.Int + py0, py1 *big.Int +} + +func (c *scalarMulConstG2Circuit) Define(api frontend.API) error { + g2, err := NewG2(api) + if err != nil { + return fmt.Errorf("new G2 struct: %w", err) + } + P := G2Affine{P: g2AffP{ + X: fields_bls12381.E2{ + A0: emulated.ValueOf[BaseField](c.px0), + A1: emulated.ValueOf[BaseField](c.px1), + }, + Y: fields_bls12381.E2{ + A0: emulated.ValueOf[BaseField](c.py0), + A1: emulated.ValueOf[BaseField](c.py1), + }, + }} + res := g2.ScalarMul(&P, &c.S) + g2.AssertIsEqual(res, &c.Res) + return nil +} + +func TestScalarMulConstG2Comb(t *testing.T) { + assert := test.NewAssert(t) + _, _, _, gen := bls12381.Generators() + var P bls12381.G2Affine + P.ScalarMultiplication(&gen, big.NewInt(12345)) + px0, px1 := P.X.A0.BigInt(new(big.Int)), P.X.A1.BigInt(new(big.Int)) + py0, py1 := P.Y.A0.BigInt(new(big.Int)), P.Y.A1.BigInt(new(big.Int)) + r := fr_bls12381.Modulus() + scalars := []*big.Int{ + big.NewInt(0), + big.NewInt(1), + big.NewInt(2), + big.NewInt(3), + new(big.Int).Sub(r, big.NewInt(1)), + new(big.Int).Sub(r, big.NewInt(2)), + new(big.Int).Lsh(big.NewInt(1), 128), + } + for _, s := range scalars { + var S bls12381.G2Affine + S.ScalarMultiplication(&P, s) + circuit := scalarMulConstG2Circuit{px0: px0, px1: px1, py0: py0, py1: py1} + witness := scalarMulConstG2Circuit{ + px0: px0, px1: px1, py0: py0, py1: py1, + S: emulated.ValueOf[ScalarField](s), + Res: NewG2Affine(S), + } + err := test.IsSolved(&circuit, &witness, ecc.BN254.ScalarField()) + assert.NoError(err, "s=%s", s.String()) + } +} + type addG2Circuit struct { In1, In2 G2Affine Res G2Affine diff --git a/std/algebra/emulated/sw_bn254/fixedbase_g2.go b/std/algebra/emulated/sw_bn254/fixedbase_g2.go index 0a6bd00e8e..a66d46129d 100644 --- a/std/algebra/emulated/sw_bn254/fixedbase_g2.go +++ b/std/algebra/emulated/sw_bn254/fixedbase_g2.go @@ -48,7 +48,11 @@ type g2CombData struct { // For the lower windows d(j) = 2j − 2^w + 1; for the top window w is // replaced by tw. windows [][][4]*big.Int - topEven [][4]*big.Int + // doubles[t][j] = 2*windows[t][j], used by the screened complete tail + // when a trailing addition degenerates into a doubling. + doubles [][][4]*big.Int + topEven [][4]*big.Int + topEvenDoubles [][4]*big.Int } var g2CombCache sync.Map @@ -105,8 +109,17 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { a.Y.A1.BigInt(new(big.Int)), }, nil } + toDoubleBig := func(a *bn254.G2Affine) ([4]*big.Int, error) { + var j bn254.G2Jac + j.FromAffine(a) + j.DoubleAssign() + var d bn254.G2Affine + d.FromJacobian(&j) + return toBig(&d) + } windows := make([][][4]*big.Int, nw) + doubles := make([][][4]*big.Int, nw) var Bt bn254.G2Jac Bt.FromAffine(&P) for t := 0; t < nw; t++ { @@ -130,6 +143,7 @@ func g2CombDataFor(x0, x1, y0, y1 *big.Int, w int) (*g2CombData, error) { odd[m].FromJacobian(acc) } tab := make([][4]*big.Int, 1< nbInc { + dbls[t] = g2.g2CombSelect(d.doubles[t], cbits[t*w:(t+1)*w]) + } } stacked := make([][4]*big.Int, 0, 2< nbInc { + dbls[nw-1] = g2.g2CombSelect(stackedDoubles, topBits) + } - nbInc := nw - 1 - d.nbUnified var res *G2Affine if nbInc >= 1 { // implicit-y chain: at the first step the accumulator is T_0 whose @@ -393,7 +461,7 @@ func (g2 *G2) scalarMulComb(d *g2CombData, s *Scalar) *G2Affine { res = &G2Affine{P: pts[0]} } for t := nbInc + 1; t <= nw-1; t++ { - res = g2.AddUnified(res, &G2Affine{P: pts[t]}) + res = g2.g2CombTailAdd(res, &pts[t], &dbls[t]) } return res } diff --git a/std/algebra/emulated/sw_bn254/g2_test.go b/std/algebra/emulated/sw_bn254/g2_test.go index 8eb94ad1b2..3a8411d92e 100644 --- a/std/algebra/emulated/sw_bn254/g2_test.go +++ b/std/algebra/emulated/sw_bn254/g2_test.go @@ -6,8 +6,10 @@ import ( "github.com/consensys/gnark-crypto/ecc" "github.com/consensys/gnark-crypto/ecc/bn254" + fr_bn "github.com/consensys/gnark-crypto/ecc/bn254/fr" "github.com/consensys/gnark/frontend" "github.com/consensys/gnark/std/algebra/algopts" + "github.com/consensys/gnark/std/algebra/emulated/fields_bn254" "github.com/consensys/gnark/std/math/emulated" "github.com/consensys/gnark/test" ) @@ -163,6 +165,65 @@ func TestScalarMulG2EdgeCases(t *testing.T) { } } +type scalarMulConstG2Circuit struct { + S Scalar + Res G2Affine + + px0, px1 *big.Int + py0, py1 *big.Int +} + +func (c *scalarMulConstG2Circuit) Define(api frontend.API) error { + g2, err := NewG2(api) + if err != nil { + panic(err) + } + P := G2Affine{P: g2AffP{ + X: fields_bn254.E2{ + A0: emulated.ValueOf[BaseField](c.px0), + A1: emulated.ValueOf[BaseField](c.px1), + }, + Y: fields_bn254.E2{ + A0: emulated.ValueOf[BaseField](c.py0), + A1: emulated.ValueOf[BaseField](c.py1), + }, + }} + res := g2.ScalarMul(&P, &c.S) + g2.AssertIsEqual(res, &c.Res) + return nil +} + +func TestScalarMulConstG2Comb(t *testing.T) { + assert := test.NewAssert(t) + _, _, _, gen := bn254.Generators() + var P bn254.G2Affine + P.ScalarMultiplication(&gen, big.NewInt(12345)) + px0, px1 := P.X.A0.BigInt(new(big.Int)), P.X.A1.BigInt(new(big.Int)) + py0, py1 := P.Y.A0.BigInt(new(big.Int)), P.Y.A1.BigInt(new(big.Int)) + r := fr_bn.Modulus() + scalars := []*big.Int{ + big.NewInt(0), + big.NewInt(1), + big.NewInt(2), + big.NewInt(3), + new(big.Int).Sub(r, big.NewInt(1)), + new(big.Int).Sub(r, big.NewInt(2)), + new(big.Int).Lsh(big.NewInt(1), 128), + } + for _, s := range scalars { + var S bn254.G2Affine + S.ScalarMultiplication(&P, s) + circuit := scalarMulConstG2Circuit{px0: px0, px1: px1, py0: py0, py1: py1} + witness := scalarMulConstG2Circuit{ + px0: px0, px1: px1, py0: py0, py1: py1, + S: emulated.ValueOf[ScalarField](s), + Res: NewG2Affine(S), + } + err := test.IsSolved(&circuit, &witness, ecc.BN254.ScalarField()) + assert.NoError(err, "s=%s", s.String()) + } +} + type scalarMulG2BySeedCircuit struct { In1 G2Affine Res G2Affine From 5c719af92f516fc6b7a4c473b49ded895749beb0 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 31 Jul 2026 11:51:26 -0400 Subject: [PATCH 14/15] fix: constrain the top recode bit in native fixed-base scalar mul --- std/algebra/native/sw_bls12377/fixedbase.go | 16 ++++--- .../sw_bls12377/fixedbase_count_test.go | 45 +++++++++++++++++++ .../native/sw_bls12377/fixedbase_g2.go | 4 +- .../native/sw_bls12377/fixedbase_g2_test.go | 21 +++++++++ 4 files changed, 79 insertions(+), 7 deletions(-) diff --git a/std/algebra/native/sw_bls12377/fixedbase.go b/std/algebra/native/sw_bls12377/fixedbase.go index d48b5b32fb..f49372d8fb 100644 --- a/std/algebra/native/sw_bls12377/fixedbase.go +++ b/std/algebra/native/sw_bls12377/fixedbase.go @@ -18,8 +18,9 @@ import ( // // - the scalar is recoded into the odd k' = s + 1 − b0 represented by n // signed binary digits, witnessed as the bits of c = (k' + 2^n − 1)/2 and -// pinned by the exact native identity 2c + b0 = s + 2^n (no wrap-around: -// 2^{n+1} is far below the native modulus); +// pinned by the exact native identity 2c + b0 = s + 2^n. The top bit of c +// is constrained to 1, proving 2c + b0 ≥ 2^n and excluding native-field +// wrap-around; // - windows of w digits select from compile-time constant tables // [d(j)·2^{w·t}]P. With constant tables the selection is a free affine // combination of the one-hot flags, so only the flag products cost @@ -198,7 +199,8 @@ func g1CombSelect(api frontend.API, table [][2]*big.Int, flags []frontend.Variab } // g1CombScalarMul computes [s]P for the constant base point of the tables d. -// It returns (0,0) when s ≡ 0 (mod r). The scalar must be reduced (s < r). +// It returns (0,0) when s ≡ 0 (mod r). The recoding constraints force s < 2^n; +// callers that need canonical scalar encodings must separately enforce s < r. func g1CombScalarMul(api frontend.API, d *g1CombData, s frontend.Variable) *G1Affine { w, n, nw := d.w, d.n, d.nw rets, err := api.Compiler().NewHint(g1CombRecodeHint, 1+n, s) @@ -208,9 +210,10 @@ func g1CombScalarMul(api frontend.API, d *g1CombData, s frontend.Variable) *G1Af b0 := rets[0] cbits := rets[1:] api.AssertIsBoolean(b0) - // exact native identity 2·c + b0 = s + 2^n: all quantities are below - // 2^{n+1} which is far below the native modulus, so the identity holds - // over the integers and pins k' = 2c − (2^n − 1) = s + 1 − b0. + // The top bit of c proves 2c+b0 ≥ 2^n. Since 2c+b0 < 2^{n+1} + // and 2^{n+1} is far below the native modulus, the equality below cannot + // be satisfied by a wrapped s+2^n value. It therefore holds over the + // integers and pins k' = 2c − (2^n − 1) = s + 1 − b0. cSum := frontend.Variable(0) coef := big.NewInt(2) for i := range cbits { @@ -218,6 +221,7 @@ func g1CombScalarMul(api frontend.API, d *g1CombData, s frontend.Variable) *G1Af cSum = api.Add(cSum, api.Mul(cbits[i], new(big.Int).Set(coef))) coef.Lsh(coef, 1) } + api.AssertIsEqual(cbits[n-1], 1) twoN := new(big.Int).Lsh(big.NewInt(1), uint(n)) api.AssertIsEqual(api.Add(cSum, b0), api.Add(s, twoN)) diff --git a/std/algebra/native/sw_bls12377/fixedbase_count_test.go b/std/algebra/native/sw_bls12377/fixedbase_count_test.go index 35902603e1..9fda68f59b 100644 --- a/std/algebra/native/sw_bls12377/fixedbase_count_test.go +++ b/std/algebra/native/sw_bls12377/fixedbase_count_test.go @@ -7,6 +7,7 @@ import ( "github.com/consensys/gnark-crypto/ecc" bls12377 "github.com/consensys/gnark-crypto/ecc/bls12-377" fr_bls "github.com/consensys/gnark-crypto/ecc/bls12-377/fr" + "github.com/consensys/gnark/constraint/solver" "github.com/consensys/gnark/frontend" "github.com/consensys/gnark/frontend/cs/r1cs" "github.com/consensys/gnark/test" @@ -24,6 +25,33 @@ func (c *nativeBaseMulCount) Define(api frontend.API) error { return nil } +func nativeCombZeroRecodeHint(_ *big.Int, _ []*big.Int, outputs []*big.Int) error { + for _, out := range outputs { + out.SetUint64(0) + } + return nil +} + +func nativeCombWrappedScalar(n int) *big.Int { + twoN := new(big.Int).Lsh(big.NewInt(1), uint(n)) + return new(big.Int).Sub(ecc.BW6_761.ScalarField(), twoN) +} + +func nativeCombNegativeTwoNScalar(n int) *big.Int { + twoN := new(big.Int).Lsh(big.NewInt(1), uint(n)) + twoN.Neg(twoN) + return twoN.Mod(twoN, fr_bls.Modulus()) +} + +func nativeCombSolveWithZeroRecode(circuit, witness frontend.Circuit) error { + return test.IsSolved( + circuit, + witness, + ecc.BW6_761.ScalarField(), + test.WithReplacementHint(solver.GetHintID(g1CombRecodeHint), nativeCombZeroRecodeHint), + ) +} + func TestNativeBaseMulCount(t *testing.T) { if testing.Short() { t.Skip() @@ -67,3 +95,20 @@ func TestNativeCombScalarMulBase(t *testing.T) { assert.NoError(err, "s=%s", s.String()) } } + +func TestNativeCombRejectsWrappedScalarRecode(t *testing.T) { + assert := test.NewAssert(t) + _, _, g, _ := bls12377.Generators() + d, err := g1CombDataFor(g.X.BigInt(new(big.Int)), g.Y.BigInt(new(big.Int))) + assert.NoError(err) + + var wrong bls12377.G1Affine + wrong.ScalarMultiplication(&g, nativeCombNegativeTwoNScalar(d.n)) + witness := nativeBaseMulCount{ + S: nativeCombWrappedScalar(d.n), + Q: G1Affine{X: wrong.X.BigInt(new(big.Int)), Y: wrong.Y.BigInt(new(big.Int))}, + } + + err = nativeCombSolveWithZeroRecode(&nativeBaseMulCount{}, &witness) + assert.Error(err, "wrapped scalar accepted with malicious all-zero comb recode") +} diff --git a/std/algebra/native/sw_bls12377/fixedbase_g2.go b/std/algebra/native/sw_bls12377/fixedbase_g2.go index 68e71a5fb4..a2835dfa55 100644 --- a/std/algebra/native/sw_bls12377/fixedbase_g2.go +++ b/std/algebra/native/sw_bls12377/fixedbase_g2.go @@ -168,7 +168,8 @@ func g2CombSelect(api frontend.API, table [][4]*big.Int, flags []frontend.Variab } // g2CombScalarMul computes [s]P for the constant base point of the tables d. -// It returns (0,0) when s ≡ 0 (mod r). The scalar must be reduced (s < r). +// It returns (0,0) when s ≡ 0 (mod r). The recoding constraints force s < 2^n; +// callers that need canonical scalar encodings must separately enforce s < r. func g2CombScalarMul(api frontend.API, d *g2CombData, s frontend.Variable) *g2AffP { w, n, nw := d.w, d.n, d.nw rets, err := api.Compiler().NewHint(g1CombRecodeHint, 1+n, s) @@ -185,6 +186,7 @@ func g2CombScalarMul(api frontend.API, d *g2CombData, s frontend.Variable) *g2Af cSum = api.Add(cSum, api.Mul(cbits[i], new(big.Int).Set(coef))) coef.Lsh(coef, 1) } + api.AssertIsEqual(cbits[n-1], 1) twoN := new(big.Int).Lsh(big.NewInt(1), uint(n)) api.AssertIsEqual(api.Add(cSum, b0), api.Add(s, twoN)) diff --git a/std/algebra/native/sw_bls12377/fixedbase_g2_test.go b/std/algebra/native/sw_bls12377/fixedbase_g2_test.go index 1f4824ba86..b8e5db4ced 100644 --- a/std/algebra/native/sw_bls12377/fixedbase_g2_test.go +++ b/std/algebra/native/sw_bls12377/fixedbase_g2_test.go @@ -69,3 +69,24 @@ func TestNativeG2CombScalarMulBase(t *testing.T) { assert.NoError(err, "s=%s", s.String()) } } + +func TestNativeG2CombRejectsWrappedScalarRecode(t *testing.T) { + assert := test.NewAssert(t) + _, _, _, g2 := bls12377.Generators() + d, err := g2CombDataFor( + g2.X.A0.BigInt(new(big.Int)), + g2.X.A1.BigInt(new(big.Int)), + g2.Y.A0.BigInt(new(big.Int)), + g2.Y.A1.BigInt(new(big.Int)), + ) + assert.NoError(err) + + var wrong bls12377.G2Affine + wrong.ScalarMultiplication(&g2, nativeCombNegativeTwoNScalar(d.n)) + var witness nativeG2BaseMulCount + witness.S = nativeCombWrappedScalar(d.n) + witness.Q.Assign(&wrong) + + err = nativeCombSolveWithZeroRecode(&nativeG2BaseMulCount{}, &witness) + assert.Error(err, "wrapped scalar accepted with malicious all-zero comb recode") +} From 02acf8bd5f5b2fd9421bb1ab0317e16a27c97fcb Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 31 Jul 2026 11:54:27 -0400 Subject: [PATCH 15/15] refactor: remove dead code --- std/algebra/native/sw_bls12377/fixedbase_g2_test.go | 3 --- 1 file changed, 3 deletions(-) diff --git a/std/algebra/native/sw_bls12377/fixedbase_g2_test.go b/std/algebra/native/sw_bls12377/fixedbase_g2_test.go index b8e5db4ced..1e4904d7f9 100644 --- a/std/algebra/native/sw_bls12377/fixedbase_g2_test.go +++ b/std/algebra/native/sw_bls12377/fixedbase_g2_test.go @@ -58,9 +58,6 @@ func TestNativeG2CombScalarMulBase(t *testing.T) { for _, s := range scalars { var S bls12377.G2Affine S.ScalarMultiplication(&g2, s) - var expected g2AffP - expected.Assign(&S) - _ = expected circuit := nativeG2BaseMulCount{} var w nativeG2BaseMulCount w.S = s