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 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..d5cebdf438 --- /dev/null +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2.go @@ -0,0 +1,575 @@ +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 + +// 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). 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 { + if _, ok := g2.api.Compiler().(frontend.PlonkAPI); ok { + return g2CombPlonkWindow + } + return g2CombWindow +} + +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). + // For the lower windows d(j) = 2j − 2^w + 1; for the top window w is + // replaced by tw. + windows [][][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 + +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 + } + 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") + } + r := fr_bls.Modulus() + rBits := r.BitLen() + nw := (rBits + w - 1) / w + 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-- { + 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++ + } + 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 + } + 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++ { + if t > 0 { + for k := 0; k < w; k++ { + 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) + 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 + dc, err := toDoubleBig(&e) + if err != nil { + return nil, fmt.Errorf("double comb table point: %w", err) + } + dbl[j] = dc + } + windows[t] = tab + doubles[t] = dbl + } + var negP bls12381.G2Jac + var negPAff bls12381.G2Affine + negPAff.Neg(&P) + negP.FromAffine(&negPAff) + topEven := 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) + } + + 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]} + } + for t := nbInc + 1; t <= nw-1; t++ { + res = g2.g2CombTailAdd(res, &pts[t], &dbls[t]) + } + return res +} + +// g2CombTryConst returns the comb tables when Q is a compile-time constant +// finite G2 subgroup point and the comb is supported in the current +// compilation context; otherwise ok is false. +func (g2 *G2) g2CombTryConst(Q *G2Affine) (*g2CombData, bool) { + if smallfields.IsSmallField(g2.api.Compiler().Field()) { + return nil, false + } + x0, ok := g2.fp.ConstantValue(&Q.P.X.A0) + if !ok { + return nil, false + } + x1, ok := g2.fp.ConstantValue(&Q.P.X.A1) + if !ok { + return nil, false + } + y0, ok := g2.fp.ConstantValue(&Q.P.Y.A0) + if !ok { + return nil, false + } + y1, ok := g2.fp.ConstantValue(&Q.P.Y.A1) + if !ok { + return nil, false + } + d, err := g2CombDataFor(x0, x1, y0, y1, g2.g2CombWindowFor()) + if err != nil { + return nil, false + } + return d, true +} + +// g2CombRecodeHint computes the comb recoding of the scalar (see the G1 +// version in sw_emulated). +func g2CombRecodeHint(_ *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)) + } + 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 + }) +} + +// 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/fixedbase_g2_test.go b/std/algebra/emulated/sw_bls12381/fixedbase_g2_test.go new file mode 100644 index 0000000000..71ed496534 --- /dev/null +++ b/std/algebra/emulated/sw_bls12381/fixedbase_g2_test.go @@ -0,0 +1,77 @@ +package sw_bls12381 + +import ( + "math/big" + "testing" + + "github.com/consensys/gnark-crypto/ecc" + 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/frontend/cs/r1cs" + "github.com/consensys/gnark/frontend/cs/scs" + "github.com/consensys/gnark/std/math/emulated" + "github.com/consensys/gnark/test" +) + +type g2FixedMulCircuit struct { + S Scalar + Q G2Affine +} + +func (c *g2FixedMulCircuit) Define(api frontend.API) error { + g2, err := NewG2(api) + if err != nil { + return err + } + _, _, _, gen := bls12381.Generators() + P := NewG2Affine(gen) + res := g2.ScalarMul(&P, &c.S) + g2.AssertIsEqual(res, &c.Q) + return nil +} + +func TestG2FixedMulCount(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + circuit := g2FixedMulCircuit{} + ccs, err := frontend.Compile(ecc.BN254.ScalarField(), r1cs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("emulated G2 ScalarMul (const gen) r1cs constraints =", ccs.GetNbConstraints()) + ccs, err = frontend.Compile(ecc.BN254.ScalarField(), scs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("emulated G2 ScalarMul (const gen) scs constraints =", ccs.GetNbConstraints()) +} + +func TestG2FixedMulCorrectness(t *testing.T) { + assert := test.NewAssert(t) + _, _, _, gen := bls12381.Generators() + r := fr_bls.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 i := 0; i < 2; i++ { + var rnd fr_bls.Element + _, _ = rnd.SetRandom() + scalars = append(scalars, rnd.BigInt(new(big.Int))) + } + for _, s := range scalars { + var S bls12381.G2Affine + S.ScalarMultiplication(&gen, s) + circuit := g2FixedMulCircuit{} + witness := g2FixedMulCircuit{ + S: emulated.ValueOf[ScalarField](s), + Q: NewG2Affine(S), + } + err := test.IsSolved(&circuit, &witness, ecc.BN254.ScalarField()) + assert.NoError(err, "s=%s", s.String()) + } +} diff --git a/std/algebra/emulated/sw_bls12381/g2.go b/std/algebra/emulated/sw_bls12381/g2.go index a81fdc5ab2..0819469f0b 100644 --- a/std/algebra/emulated/sw_bls12381/g2.go +++ b/std/algebra/emulated/sw_bls12381/g2.go @@ -666,6 +666,14 @@ func (g2 *G2) scalarMulGeneric(p *G2Affine, s *Scalar, opts ...algopts.AlgebraOp // // [EEMP25]: https://eprint.iacr.org/2025/933 func (g2 *G2) ScalarMul(Q *G2Affine, s *Scalar, opts ...algopts.AlgebraOption) *G2Affine { + // when Q is a compile-time constant subgroup point (for example a fixed + // verification-key or SRS point), use the fixed-base comb method with + // precomputed tables (see fixedbase_g2.go). It uses complete arithmetic + // and handles the zero scalar; [algopts.WithIncompleteArithmetic] is a + // no-op on this path. + if d, ok := g2.g2CombTryConst(Q); ok { + return g2.scalarMulComb(d, s) + } return g2.scalarMulGLVAndFakeGLV(Q, s, opts...) } 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_bls12381/hints.go b/std/algebra/emulated/sw_bls12381/hints.go index 051f7d0d73..9e2abb6d9d 100644 --- a/std/algebra/emulated/sw_bls12381/hints.go +++ b/std/algebra/emulated/sw_bls12381/hints.go @@ -24,6 +24,8 @@ func GetHints() []solver.Hint { pairingCheckHint, 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 new file mode 100644 index 0000000000..a66d46129d --- /dev/null +++ b/std/algebra/emulated/sw_bn254/fixedbase_g2.go @@ -0,0 +1,575 @@ +package sw_bn254 + +import ( + "errors" + "fmt" + "math/big" + "sync" + + bn254 "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/internal/smallfields" + "github.com/consensys/gnark/std/algebra/emulated/fields_bn254" + 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 + +// 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). 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 { + if _, ok := g2.api.Compiler().(frontend.PlonkAPI); ok { + return g2CombPlonkWindow + } + return g2CombWindow +} + +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). + // For the lower windows d(j) = 2j − 2^w + 1; for the top window w is + // replaced by tw. + windows [][][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 + +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 + } + var P bn254.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") + } + r := fr_bn.Modulus() + rBits := r.BitLen() + nw := (rBits + w - 1) / w + 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-- { + 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++ + } + 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 + } + 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++ { + if t > 0 { + for k := 0; k < w; k++ { + 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) + 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 + dc, err := toDoubleBig(&e) + if err != nil { + return nil, fmt.Errorf("double comb table point: %w", err) + } + dbl[j] = dc + } + windows[t] = tab + doubles[t] = dbl + } + var negP bn254.G2Jac + var negPAff bn254.G2Affine + negPAff.Neg(&P) + negP.FromAffine(&negPAff) + topEven := 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) + } + + 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]} + } + for t := nbInc + 1; t <= nw-1; t++ { + res = g2.g2CombTailAdd(res, &pts[t], &dbls[t]) + } + return res +} + +// g2CombTryConst returns the comb tables when Q is a compile-time constant +// finite G2 subgroup point and the comb is supported in the current +// compilation context; otherwise ok is false. +func (g2 *G2) g2CombTryConst(Q *G2Affine) (*g2CombData, bool) { + if smallfields.IsSmallField(g2.api.Compiler().Field()) { + return nil, false + } + x0, ok := g2.fp.ConstantValue(&Q.P.X.A0) + if !ok { + return nil, false + } + x1, ok := g2.fp.ConstantValue(&Q.P.X.A1) + if !ok { + return nil, false + } + y0, ok := g2.fp.ConstantValue(&Q.P.Y.A0) + if !ok { + return nil, false + } + y1, ok := g2.fp.ConstantValue(&Q.P.Y.A1) + if !ok { + return nil, false + } + d, err := g2CombDataFor(x0, x1, y0, y1, g2.g2CombWindowFor()) + if err != nil { + return nil, false + } + return d, true +} + +// g2CombRecodeHint computes the comb recoding of the scalar (see the G1 +// version in sw_emulated). +func g2CombRecodeHint(_ *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)) + } + 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 + }) +} + +// 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/fixedbase_g2_test.go b/std/algebra/emulated/sw_bn254/fixedbase_g2_test.go new file mode 100644 index 0000000000..86a7335e18 --- /dev/null +++ b/std/algebra/emulated/sw_bn254/fixedbase_g2_test.go @@ -0,0 +1,77 @@ +package sw_bn254 + +import ( + "math/big" + "testing" + + "github.com/consensys/gnark-crypto/ecc" + bn254c "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/frontend/cs/r1cs" + "github.com/consensys/gnark/frontend/cs/scs" + "github.com/consensys/gnark/std/math/emulated" + "github.com/consensys/gnark/test" +) + +type g2FixedMulCircuit struct { + S Scalar + Q G2Affine +} + +func (c *g2FixedMulCircuit) Define(api frontend.API) error { + g2, err := NewG2(api) + if err != nil { + return err + } + _, _, _, gen := bn254c.Generators() + P := NewG2Affine(gen) + res := g2.ScalarMul(&P, &c.S) + g2.AssertIsEqual(res, &c.Q) + return nil +} + +func TestG2FixedMulCount(t *testing.T) { + if testing.Short() { + t.Skip() + } + assert := test.NewAssert(t) + circuit := g2FixedMulCircuit{} + ccs, err := frontend.Compile(ecc.BN254.ScalarField(), r1cs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("bn254 emulated G2 ScalarMul (const gen) r1cs constraints =", ccs.GetNbConstraints()) + ccs, err = frontend.Compile(ecc.BN254.ScalarField(), scs.NewBuilder, &circuit) + assert.NoError(err) + t.Log("bn254 emulated G2 ScalarMul (const gen) scs constraints =", ccs.GetNbConstraints()) +} + +func TestG2FixedMulCorrectness(t *testing.T) { + assert := test.NewAssert(t) + _, _, _, gen := bn254c.Generators() + 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 i := 0; i < 2; i++ { + var rnd fr_bn.Element + _, _ = rnd.SetRandom() + scalars = append(scalars, rnd.BigInt(new(big.Int))) + } + for _, s := range scalars { + var S bn254c.G2Affine + S.ScalarMultiplication(&gen, s) + circuit := g2FixedMulCircuit{} + witness := g2FixedMulCircuit{ + S: emulated.ValueOf[ScalarField](s), + Q: NewG2Affine(S), + } + err := test.IsSolved(&circuit, &witness, ecc.BN254.ScalarField()) + assert.NoError(err, "s=%s", s.String()) + } +} diff --git a/std/algebra/emulated/sw_bn254/g2.go b/std/algebra/emulated/sw_bn254/g2.go index 208d704098..6353945e01 100644 --- a/std/algebra/emulated/sw_bn254/g2.go +++ b/std/algebra/emulated/sw_bn254/g2.go @@ -451,6 +451,14 @@ func (g2 *G2) Select(b frontend.Variable, p, q *G2Affine) *G2Affine { // // [EEMP25]: https://eprint.iacr.org/2025/933 func (g2 *G2) ScalarMul(Q *G2Affine, s *Scalar, opts ...algopts.AlgebraOption) *G2Affine { + // when Q is a compile-time constant subgroup point (for example a fixed + // verification-key or SRS point), use the fixed-base comb method with + // precomputed tables (see fixedbase_g2.go). It uses complete arithmetic + // and handles the zero scalar; [algopts.WithIncompleteArithmetic] is a + // no-op on this path. + if d, ok := g2.g2CombTryConst(Q); ok { + return g2.scalarMulComb(d, s) + } return g2.scalarMulGLVAndFakeGLV(Q, s, opts...) } 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 diff --git a/std/algebra/emulated/sw_bn254/hints.go b/std/algebra/emulated/sw_bn254/hints.go index 329f2b908c..25f557a07a 100644 --- a/std/algebra/emulated/sw_bn254/hints.go +++ b/std/algebra/emulated/sw_bn254/hints.go @@ -22,6 +22,8 @@ func GetHints() []solver.Hint { pairingCheckHint, millerLoopAndCheckFinalExpHint, scalarMulG2Hint, + g2CombRecodeHint, + g2CombChainLambdaHint, rationalReconstructExtG2, } } diff --git a/std/algebra/emulated/sw_emulated/fixedbase.go b/std/algebra/emulated/sw_emulated/fixedbase.go new file mode 100644 index 0000000000..871dee7600 --- /dev/null +++ b/std/algebra/emulated/sw_emulated/fixedbase.go @@ -0,0 +1,598 @@ +package sw_emulated + +import ( + "errors" + "fmt" + "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" +) + +// 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= 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, 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 { + return nil, fmt.Errorf("window width %d too large for %d-bit scalar field", w, rBits) + } + n := rBits + if n > scalarCap { + return nil, fmt.Errorf("recoded scalar needs %d bits, scalar field emulation has capacity %d", n, scalarCap) + } + 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-- { + 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++ + } + nbUnified = max(nbUnified, 1) + + G := &combAffine{x: new(big.Int).Set(gx), y: new(big.Int).Set(gy)} + 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 + } + } + } + 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++ { + 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..4cb347839e --- /dev/null +++ b/std/algebra/emulated/sw_emulated/fixedbase_test.go @@ -0,0 +1,474 @@ +package sw_emulated + +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" + 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, 5, 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, 5, 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()) +} + +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()) + } +} + +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") +} + +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()) +} + +// 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/hints.go b/std/algebra/emulated/sw_emulated/hints.go index e436f992cc..9526608b82 100644 --- a/std/algebra/emulated/sw_emulated/hints.go +++ b/std/algebra/emulated/sw_emulated/hints.go @@ -34,6 +34,15 @@ func GetHints() []solver.Hint { scalarMulHint, rationalReconstruct, rationalReconstructExt, + combRecodeHint, + combChainLambdaHint, + ratioHint, + tangentHint, + tangentHintA, + unifiedSlopeHint, + bjSlopeHint, + implicitTangentHint, + implicitChordHint, } } 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..9a47bb0a10 --- /dev/null +++ b/std/algebra/emulated/sw_emulated/joint_comb_count_test.go @@ -0,0 +1,121 @@ +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()) + } +} + +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()) +} + +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 9ed6f6b662..b1dfbcee7e 100644 --- a/std/algebra/emulated/sw_emulated/point.go +++ b/std/algebra/emulated/sw_emulated/point.go @@ -76,6 +76,10 @@ 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 base point and window + // width. + combCache map[string]*combData } // Generator returns the base point of the curve. The method does not copy and @@ -191,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}) @@ -273,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 λ @@ -320,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(λ, λ) @@ -372,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}) @@ -432,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) } @@ -448,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) } @@ -492,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) } @@ -510,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) } @@ -541,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}) @@ -554,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 @@ -643,7 +653,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, c.combWindow()); err == nil { + return c.scalarMulComb(d, s) + } + } + } if c.eigenvalue != nil && c.thirdRootOne != nil { return c.scalarMulGLVAndFakeGLV(p, s, opts...) @@ -1172,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( @@ -1198,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 @@ -1233,20 +1274,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(c.combWindow()); err == nil { + return c.scalarMulBaseComb(s, c.combWindow()) + } if c.eigenvalue != nil && c.thirdRootOne != nil { return c.scalarMulGLVAndFakeGLV(c.Generator(), s, opts...) @@ -1283,6 +1326,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(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 + // 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, c.combWindow()) + sm2 := c.ScalarMul(p, s2, opts...) + return c.AddUnified(sm1, sm2) + } return c.jointScalarMul(c.Generator(), p, s1, s2, opts...) } @@ -1319,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, c.combWindow()); 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 { @@ -1848,6 +1947,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( @@ -1877,10 +1983,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 new file mode 100644 index 0000000000..dbeafd3893 --- /dev/null +++ b/std/algebra/emulated/sw_emulated/slopes.go @@ -0,0 +1,335 @@ +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 + }) +} + +// 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 + }) +} diff --git a/std/algebra/native/sw_bls12377/fixedbase.go b/std/algebra/native/sw_bls12377/fixedbase.go new file mode 100644 index 0000000000..f49372d8fb --- /dev/null +++ b/std/algebra/native/sw_bls12377/fixedbase.go @@ -0,0 +1,282 @@ +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. 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 +// 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<= 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 := 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/g2.go b/std/algebra/native/sw_bls12377/g2.go index ab7901fc9b..51a50050fa 100644 --- a/std/algebra/native/sw_bls12377/g2.go +++ b/std/algebra/native/sw_bls12377/g2.go @@ -203,9 +203,25 @@ func (p *g2AffP) scalarMulBySeed(api frontend.API, Q *g2AffP) *g2AffP { func (p *g2AffP) ScalarMul(api frontend.API, Q g2AffP, s interface{}, opts ...algopts.AlgebraOption) *g2AffP { if n, ok := api.Compiler().ConstantValue(s); ok { return p.constScalarMul(api, Q, n, opts...) - } else { - return p.varScalarMul(api, Q, s, opts...) } + // when Q is a compile-time constant point, use the fixed-base comb + // method with precomputed tables (complete arithmetic, handles the zero + // scalar). See fixedbase_g2.go. + if x0, ok := api.Compiler().ConstantValue(Q.X.A0); ok { + if x1, ok := api.Compiler().ConstantValue(Q.X.A1); ok { + if y0, ok := api.Compiler().ConstantValue(Q.Y.A0); ok { + if y1, ok := api.Compiler().ConstantValue(Q.Y.A1); ok { + if d, err := g2CombDataFor(x0, x1, y0, y1); err == nil && api.Compiler().FieldBitLen() > 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() 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, diff --git a/std/math/emulated/field.go b/std/math/emulated/field.go index e23a29ca1f..c22f4d8f9b 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() { @@ -298,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, 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) }