Feat: emulated subgroup check

go.mod

@@ -6,7 +6,7 @@ require (
 	github.com/bits-and-blooms/bitset v1.5.0
 	github.com/blang/semver/v4 v4.0.0
 	github.com/consensys/bavard v0.1.13
-	github.com/consensys/gnark-crypto v0.9.2-0.20230329155745-a57dcc3b53de
+	github.com/consensys/gnark-crypto v0.10.1-0.20230411152222-ebbf6920ad59
 	github.com/fxamacker/cbor/v2 v2.4.0
 	github.com/google/go-cmp v0.5.9
 	github.com/google/pprof v0.0.0-20230309165930-d61513b1440d

go.sum

@@ -6,6 +6,8 @@ github.com/consensys/bavard v0.1.13 h1:oLhMLOFGTLdlda/kma4VOJazblc7IM5y5QPd2A/Yj
 github.com/consensys/bavard v0.1.13/go.mod h1:9ItSMtA/dXMAiL7BG6bqW2m3NdSEObYWoH223nGHukI=
 github.com/consensys/gnark-crypto v0.9.2-0.20230329155745-a57dcc3b53de h1:W5lRxU8Rk8CDLHMTeyNst0VESbcU5RZ3U1TS9MNGgCQ=
 github.com/consensys/gnark-crypto v0.9.2-0.20230329155745-a57dcc3b53de/go.mod h1:Iq/P3HHl0ElSjsg2E1gsMwhAyxnxoKK5nVyZKd+/KhU=
+github.com/consensys/gnark-crypto v0.10.1-0.20230411152222-ebbf6920ad59 h1:KYAxBDzdiFQWO+lptSBb+7/UUHKjObS87Bx7pBGX1DU=
+github.com/consensys/gnark-crypto v0.10.1-0.20230411152222-ebbf6920ad59/go.mod h1:Iq/P3HHl0ElSjsg2E1gsMwhAyxnxoKK5nVyZKd+/KhU=
 github.com/coreos/go-systemd/v22 v22.3.3-0.20220203105225-a9a7ef127534/go.mod h1:Y58oyj3AT4RCenI/lSvhwexgC+NSVTIJ3seZv2GcEnc=
 github.com/creack/pty v1.1.9/go.mod h1:oKZEueFk5CKHvIhNR5MUki03XCEU+Q6VDXinZuGJ33E=
 github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=

std/algebra/emulated/sw_bls12381/hints.go

@@ -0,0 +1,68 @@
+package sw_bls12381
+
+import (
+	"math/big"
+
+	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
+	"github.com/consensys/gnark/constraint/solver"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+func init() {
+	solver.RegisterHint(GetHints()...)
+}
+
+// GetHints returns all hint functions used in the package.
+func GetHints() []solver.Hint {
+	return []solver.Hint{
+		subgroupG1Hint,
+		subgroupG2Hint,
+	}
+}
+
+func subgroupG1Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var a, c bls12381.G1Affine
+
+			a.X.SetBigInt(inputs[0])
+			a.Y.SetBigInt(inputs[1])
+
+			// c = -[x²]ϕ(p)
+			x0, _ := new(big.Int).SetString("15132376222941642752", 10) // negative
+			var jac bls12381.G1Jac
+			jac.FromAffine(&a)
+			jac.Phi(&jac).ScalarMultiplication(&jac, x0).
+				ScalarMultiplication(&jac, x0).
+				Neg(&jac)
+			c.FromJacobian(&jac)
+
+			c.X.BigInt(outputs[0])
+			c.Y.BigInt(outputs[1])
+
+			return nil
+		})
+}
+
+func subgroupG2Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var a, c bls12381.G2Affine
+
+			a.X.A0.SetBigInt(inputs[0])
+			a.X.A1.SetBigInt(inputs[1])
+			a.Y.A0.SetBigInt(inputs[2])
+			a.Y.A1.SetBigInt(inputs[3])
+
+			// c = [x₀]a
+			x0, _ := new(big.Int).SetString("15132376222941642752", 10) // negative
+			c.ScalarMultiplication(&a, x0).Neg(&c)
+
+			c.X.A0.BigInt(outputs[0])
+			c.X.A1.BigInt(outputs[1])
+			c.Y.A0.BigInt(outputs[2])
+			c.Y.A1.BigInt(outputs[3])
+
+			return nil
+		})
+}

std/algebra/emulated/sw_bls12381/pairing.go

@@ -227,6 +227,100 @@ func (pr Pairing) AssertIsEqual(x, y *GTEl) {
 	pr.Ext12.AssertIsEqual(x, y)
 }
 
+func (pr Pairing) AssertIsOnCurve(P *G1Affine) {
+	// Curve: Y² == X³ + aX + b, where a=0 and b=4
+	// (X,Y) ∈ {Y² == X³ + aX + b} U (0,0)
+
+	// if P=(0,0) we assign b=0 otherwise 3, and continue
+	selector := pr.api.And(pr.curveF.IsZero(&P.X), pr.curveF.IsZero(&P.Y))
+	bCurve := emulated.ValueOf[emulated.BLS12381Fp](4)
+	b := pr.curveF.Select(selector, pr.curveF.Zero(), &bCurve)
+
+	left := pr.curveF.Mul(&P.Y, &P.Y)
+	right := pr.curveF.Mul(&P.X, &P.X)
+	right = pr.curveF.Mul(right, &P.X)
+	right = pr.curveF.Add(right, b)
+	pr.curveF.AssertIsEqual(left, right)
+}
+
+func (pr Pairing) AssertIsOnTwist(Q *G2Affine) {
+	// Twist: Y² == X³ + aX + b, where a=0 and b=4(1+u)
+	// (X,Y) ∈ {Y² == X³ + aX + b} U (0,0)
+
+	// if Q=(0,0) we assign b=0 otherwise 3/(9+u), and continue
+	selector := pr.api.And(pr.Ext2.IsZero(&Q.X), pr.Ext2.IsZero(&Q.Y))
+	bTwist := fields_bls12381.E2{
+		A0: emulated.ValueOf[emulated.BLS12381Fp]("4"),
+		A1: emulated.ValueOf[emulated.BLS12381Fp]("4"),
+	}
+	b := pr.Ext2.Select(selector, pr.Ext2.Zero(), &bTwist)
+
+	left := pr.Ext2.Square(&Q.Y)
+	right := pr.Ext2.Square(&Q.X)
+	right = pr.Ext2.Mul(right, &Q.X)
+	right = pr.Ext2.Add(right, b)
+	pr.Ext2.AssertIsEqual(left, right)
+}
+
+func (pr Pairing) AssertIsOnG1(P *G1Affine) {
+	// 1- Check P is on the curve
+	pr.AssertIsOnCurve(P)
+
+	// 2- Check P has the right subgroup order
+	res, err := pr.curveF.NewHint(subgroupG1Hint, 2, &P.X, &P.Y)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	// _P = -[x²]ϕ(p)
+	_P := G1Affine{
+		X: *res[0],
+		Y: *res[1],
+	}
+
+	// [r]Q == 0 <==>  P = _P
+	pr.curveF.AssertIsEqual(&P.X, &_P.X)
+	pr.curveF.AssertIsEqual(&P.Y, &_P.Y)
+}
+
+func (pr Pairing) AssertIsOnG2(Q *G2Affine) {
+	// 1- Check Q is on the curve
+	pr.AssertIsOnTwist(Q)
+
+	// 2- Check Q has the right subgroup order
+	res, err := pr.curveF.NewHint(subgroupG2Hint, 4, &Q.X.A0, &Q.X.A1, &Q.Y.A0, &Q.Y.A1)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	// _Q = [x₀]Q
+	_Q := G2Affine{
+		X: fields_bls12381.E2{A0: *res[0], A1: *res[1]},
+		Y: fields_bls12381.E2{A0: *res[2], A1: *res[3]},
+	}
+
+	// [r]Q == 0 <==>  ψ(Q) = _Q
+	var psiQ G2Affine
+	u1 := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939437")
+	v := fields_bls12381.E2{
+		A0: emulated.ValueOf[emulated.BLS12381Fp]("2973677408986561043442465346520108879172042883009249989176415018091420807192182638567116318576472649347015917690530"),
+		A1: emulated.ValueOf[emulated.BLS12381Fp]("1028732146235106349975324479215795277384839936929757896155643118032610843298655225875571310552543014690878354869257"),
+	}
+
+	// ψ(Q)
+	psiQ.X = *pr.Ext2.Conjugate(&Q.X)
+	tmp := *pr.curveF.Mul(&psiQ.X.A1, &u1)
+	psiQ.X.A1 = *pr.curveF.Mul(&psiQ.X.A0, &u1)
+	psiQ.X.A0 = *pr.curveF.Neg(&tmp)
+	psiQ.Y = *pr.Ext2.Conjugate(&Q.Y)
+	psiQ.Y = *pr.Ext2.Mul(&psiQ.Y, &v)
+
+	pr.Ext2.AssertIsEqual(&psiQ.X, &_Q.X)
+	pr.Ext2.AssertIsEqual(&psiQ.Y, &_Q.Y)
+}
+
 // loopCounter = seed in binary
 //
 //	seed=-15132376222941642752

std/algebra/emulated/sw_bls12381/pairing_test.go

@@ -66,6 +66,8 @@ func (c *PairCircuit) Define(api frontend.API) error {
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
+	pairing.AssertIsOnG1(&c.InG1)
+	pairing.AssertIsOnG2(&c.InG2)
 	res, err := pairing.Pair([]*G1Affine{&c.InG1}, []*G2Affine{&c.InG2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
@@ -101,6 +103,10 @@ func (c *MultiPairCircuit) Define(api frontend.API) error {
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
+	pairing.AssertIsOnG1(&c.In1G1)
+	pairing.AssertIsOnG1(&c.In2G1)
+	pairing.AssertIsOnG2(&c.In1G2)
+	pairing.AssertIsOnG2(&c.In2G2)
 	res, err := pairing.Pair([]*G1Affine{&c.In1G1, &c.In1G1, &c.In2G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2, &c.In1G2, &c.In2G2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
@@ -196,3 +202,29 @@ func TestFinalExponentiationSafeCircuit(t *testing.T) {
 	}, ecc.BN254.ScalarField())
 	assert.NoError(err)
 }
+
+type GroupMembershipCircuit struct {
+	InG1 G1Affine
+	InG2 G2Affine
+}
+
+func (c *GroupMembershipCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	pairing.AssertIsOnG1(&c.InG1)
+	pairing.AssertIsOnG2(&c.InG2)
+	return nil
+}
+
+func TestGroupMembershipSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p, q := randomG1G2Affines(assert)
+	witness := GroupMembershipCircuit{
+		InG1: NewG1Affine(p),
+		InG2: NewG2Affine(q),
+	}
+	err := test.IsSolved(&GroupMembershipCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}

std/algebra/emulated/sw_bn254/hints.go

@@ -0,0 +1,43 @@
+package sw_bn254
+
+import (
+	"math/big"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254"
+	"github.com/consensys/gnark/constraint/solver"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+func init() {
+	solver.RegisterHint(subgroupG2Hint)
+}
+
+func subgroupG2Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var a, c bn254.G2Affine
+
+			a.X.A0.SetBigInt(inputs[0])
+			a.X.A1.SetBigInt(inputs[1])
+			a.Y.A0.SetBigInt(inputs[2])
+			a.Y.A1.SetBigInt(inputs[3])
+
+			// c = ψ³([2x₀]a) - ψ²([x₀]a) - ψ([x₀]a) - [x₀]a
+			x0, _ := new(big.Int).SetString("4965661367192848881", 10)
+			var aJac, t1, t2, t3 bn254.G2Jac
+			aJac.FromAffine(&a)
+			aJac.ScalarMultiplication(&aJac, x0)
+			t1.Psi(&aJac)
+			t2.Psi(&t1)
+			t3.Psi(&t2).Double(&t3).Neg(&t3)
+			aJac.AddAssign(&t1).AddAssign(&t2).AddAssign(&t3).Neg(&aJac)
+			c.FromJacobian(&aJac)
+
+			c.X.A0.BigInt(outputs[0])
+			c.X.A1.BigInt(outputs[1])
+			c.Y.A0.BigInt(outputs[2])
+			c.Y.A1.BigInt(outputs[3])
+
+			return nil
+		})
+}

std/algebra/emulated/sw_bn254/pairing.go

@@ -202,7 +202,7 @@ func (pr Pairing) finalExponentiation(e *GTEl, unsafe bool) *GTEl {
 // Pair calculates the reduced pairing for a set of points
 // ∏ᵢ e(Pᵢ, Qᵢ).
 //
-// This function doesn't check that the inputs are in the correct subgroups.
+// This function doesn't check that the inputs are in the correct subgroups. See AssertIsOnG1 and AssertIsOnG2.
 func (pr Pairing) Pair(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
 	res, err := pr.MillerLoop(P, Q)
 	if err != nil {
@@ -215,7 +215,7 @@ func (pr Pairing) Pair(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
 // PairingCheck calculates the reduced pairing for a set of points and asserts if the result is One
 // ∏ᵢ e(Pᵢ, Qᵢ) =? 1
 //
-// This function doesn't check that the inputs are in the correct subgroups.
+// This function doesn't check that the inputs are in the correct subgroups. See AssertIsOnG1 and AssertIsOnG2.
 func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	f, err := pr.Pair(P, Q)
 	if err != nil {
@@ -232,6 +232,69 @@ func (pr Pairing) AssertIsEqual(x, y *GTEl) {
 	pr.Ext12.AssertIsEqual(x, y)
 }
 
+func (pr Pairing) AssertIsOnCurve(P *G1Affine) {
+	// Curve: Y² == X³ + aX + b, where a=0 and b=3
+	// (X,Y) ∈ {Y² == X³ + aX + b} U (0,0)
+
+	// if P=(0,0) we assign b=0 otherwise 3, and continue
+	selector := pr.api.And(pr.curveF.IsZero(&P.X), pr.curveF.IsZero(&P.Y))
+	bCurve := emulated.ValueOf[emulated.BN254Fp](3)
+	b := pr.curveF.Select(selector, pr.curveF.Zero(), &bCurve)
+
+	left := pr.curveF.Mul(&P.Y, &P.Y)
+	right := pr.curveF.Mul(&P.X, &P.X)
+	right = pr.curveF.Mul(right, &P.X)
+	right = pr.curveF.Add(right, b)
+	pr.curveF.AssertIsEqual(left, right)
+}
+
+func (pr Pairing) AssertIsOnTwist(Q *G2Affine) {
+	// Twist: Y² == X³ + aX + b, where a=0 and b=3/(9+u)
+	// (X,Y) ∈ {Y² == X³ + aX + b} U (0,0)
+
+	// if Q=(0,0) we assign b=0 otherwise 3/(9+u), and continue
+	selector := pr.api.And(pr.Ext2.IsZero(&Q.X), pr.Ext2.IsZero(&Q.Y))
+	bTwist := fields_bn254.E2{
+		A0: emulated.ValueOf[emulated.BN254Fp]("19485874751759354771024239261021720505790618469301721065564631296452457478373"),
+		A1: emulated.ValueOf[emulated.BN254Fp]("266929791119991161246907387137283842545076965332900288569378510910307636690"),
+	}
+	b := pr.Ext2.Select(selector, pr.Ext2.Zero(), &bTwist)
+
+	left := pr.Ext2.Square(&Q.Y)
+	right := pr.Ext2.Square(&Q.X)
+	right = pr.Ext2.Mul(right, &Q.X)
+	right = pr.Ext2.Add(right, b)
+	pr.Ext2.AssertIsEqual(left, right)
+}
+
+func (pr Pairing) AssertIsOnG1(P *G1Affine) {
+	// BN254 has a prime order, so we only
+	// 1- Check P is on the curve
+	pr.AssertIsOnCurve(P)
+}
+
+func (pr Pairing) AssertIsOnG2(Q *G2Affine) {
+	// 1- Check Q is on the curve
+	pr.AssertIsOnTwist(Q)
+
+	// 2- Check Q has the right subgroup order
+	res, err := pr.curveF.NewHint(subgroupG2Hint, 4, &Q.X.A0, &Q.X.A1, &Q.Y.A0, &Q.Y.A1)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	// _Q = ψ³([2x₀]Q) - ψ²([x₀]Q) - ψ([x₀]Q) - [x₀]Q
+	_Q := G2Affine{
+		X: fields_bn254.E2{A0: *res[0], A1: *res[1]},
+		Y: fields_bn254.E2{A0: *res[2], A1: *res[3]},
+	}
+
+	// [r]Q == 0 <==>  Q = _Q
+	pr.Ext2.AssertIsEqual(&Q.X, &_Q.X)
+	pr.Ext2.AssertIsEqual(&Q.Y, &_Q.Y)
+}
+
 // loopCounter = 6x₀+2 = 29793968203157093288
 //
 // in 2-NAF
@@ -436,15 +499,15 @@ func (pr Pairing) MillerLoop(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
 	Q1, Q2 := new(G2Affine), new(G2Affine)
 	for k := 0; k < n; k++ {
 		//Q1 = π(Q)
-		Q1.X = *pr.Ext12.Ext2.Conjugate(&Q[k].X)
-		Q1.X = *pr.Ext12.Ext2.MulByNonResidue1Power2(&Q1.X)
-		Q1.Y = *pr.Ext12.Ext2.Conjugate(&Q[k].Y)
-		Q1.Y = *pr.Ext12.Ext2.MulByNonResidue1Power3(&Q1.Y)
+		Q1.X = *pr.Ext2.Conjugate(&Q[k].X)
+		Q1.X = *pr.Ext2.MulByNonResidue1Power2(&Q1.X)
+		Q1.Y = *pr.Ext2.Conjugate(&Q[k].Y)
+		Q1.Y = *pr.Ext2.MulByNonResidue1Power3(&Q1.Y)
 
 		// Q2 = -π²(Q)
-		Q2.X = *pr.Ext12.Ext2.MulByNonResidue2Power2(&Q[k].X)
-		Q2.Y = *pr.Ext12.Ext2.MulByNonResidue2Power3(&Q[k].Y)
-		Q2.Y = *pr.Ext12.Ext2.Neg(&Q2.Y)
+		Q2.X = *pr.Ext2.MulByNonResidue2Power2(&Q[k].X)
+		Q2.Y = *pr.Ext2.MulByNonResidue2Power3(&Q[k].Y)
+		Q2.Y = *pr.Ext2.Neg(&Q2.Y)
 
 		// Qacc[k] ← Qacc[k]+π(Q) and
 		// l1 the line passing Qacc[k] and π(Q)