diff --git a/std/algebra/emulated/sw_bls12381/doc_test.go b/std/algebra/emulated/sw_bls12381/doc_test.go
index ae87b72997..a1ce0f5ca5 100644
--- a/std/algebra/emulated/sw_bls12381/doc_test.go
+++ b/std/algebra/emulated/sw_bls12381/doc_test.go
@@ -23,6 +23,10 @@ func (c *PairCircuit) Define(api frontend.API) error {
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
+	// Pair method does not check that the points are in the proper groups.
+	pairing.AssertIsOnG1(&c.InG1)
+	pairing.AssertIsOnG2(&c.InG2)
+	// Compute the pairing
 	res, err := pairing.Pair([]*sw_bls12381.G1Affine{&c.InG1}, []*sw_bls12381.G2Affine{&c.InG2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
diff --git a/std/algebra/emulated/sw_bls12381/g1.go b/std/algebra/emulated/sw_bls12381/g1.go
index fd95f4722b..f7908b5ded 100644
--- a/std/algebra/emulated/sw_bls12381/g1.go
+++ b/std/algebra/emulated/sw_bls12381/g1.go
@@ -1,7 +1,10 @@
 package sw_bls12381
 
 import (
+	"fmt"
+
 	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
+	"github.com/consensys/gnark/frontend"
 	"github.com/consensys/gnark/std/algebra/emulated/sw_emulated"
 	"github.com/consensys/gnark/std/math/emulated"
 )
@@ -14,3 +17,29 @@ func NewG1Affine(v bls12381.G1Affine) G1Affine {
 		Y: emulated.ValueOf[emulated.BLS12381Fp](v.Y),
 	}
 }
+
+type G1 struct {
+	curveF *emulated.Field[emulated.BLS12381Fp]
+	w      *emulated.Element[emulated.BLS12381Fp]
+}
+
+func NewG1(api frontend.API) (*G1, error) {
+	ba, err := emulated.NewField[emulated.BLS12381Fp](api)
+	if err != nil {
+		return nil, fmt.Errorf("new base api: %w", err)
+	}
+	w := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939436")
+	return &G1{
+		curveF: ba,
+		w:      &w,
+	}, nil
+}
+
+func (g1 *G1) phi(q *G1Affine) *G1Affine {
+	x := g1.curveF.Mul(&q.X, g1.w)
+
+	return &G1Affine{
+		X: *x,
+		Y: q.Y,
+	}
+}
diff --git a/std/algebra/emulated/sw_bls12381/g2.go b/std/algebra/emulated/sw_bls12381/g2.go
index f96ea4c354..45077a6cba 100644
--- a/std/algebra/emulated/sw_bls12381/g2.go
+++ b/std/algebra/emulated/sw_bls12381/g2.go
@@ -1,15 +1,39 @@
 package sw_bls12381
 
 import (
+	"math/big"
+
 	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
+	"github.com/consensys/gnark/frontend"
 	"github.com/consensys/gnark/std/algebra/emulated/fields_bls12381"
 	"github.com/consensys/gnark/std/math/emulated"
 )
 
+type G2 struct {
+	*fields_bls12381.Ext2
+	u1, w *emulated.Element[emulated.BLS12381Fp]
+	v     *fields_bls12381.E2
+}
+
 type G2Affine struct {
 	X, Y fields_bls12381.E2
 }
 
+func NewG2(api frontend.API) *G2 {
+	w := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939436")
+	u1 := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939437")
+	v := fields_bls12381.E2{
+		A0: emulated.ValueOf[emulated.BLS12381Fp]("2973677408986561043442465346520108879172042883009249989176415018091420807192182638567116318576472649347015917690530"),
+		A1: emulated.ValueOf[emulated.BLS12381Fp]("1028732146235106349975324479215795277384839936929757896155643118032610843298655225875571310552543014690878354869257"),
+	}
+	return &G2{
+		Ext2: fields_bls12381.NewExt2(api),
+		w:    &w,
+		u1:   &u1,
+		v:    &v,
+	}
+}
+
 func NewG2Affine(v bls12381.G2Affine) G2Affine {
 	return G2Affine{
 		X: fields_bls12381.E2{
@@ -22,3 +46,174 @@ func NewG2Affine(v bls12381.G2Affine) G2Affine {
 		},
 	}
 }
+
+func (g2 *G2) psi(q *G2Affine) *G2Affine {
+	x := g2.Ext2.MulByElement(&q.X, g2.u1)
+	y := g2.Ext2.Conjugate(&q.Y)
+	y = g2.Ext2.Mul(y, g2.v)
+
+	return &G2Affine{
+		X: fields_bls12381.E2{A0: x.A1, A1: x.A0},
+		Y: *y,
+	}
+}
+
+func (g2 *G2) scalarMulBySeed(q *G2Affine) *G2Affine {
+
+	z := g2.triple(q)
+	z = g2.double(z)
+	z = g2.doubleAndAdd(z, q)
+	z = g2.doubleN(z, 2)
+	z = g2.doubleAndAdd(z, q)
+	z = g2.doubleN(z, 8)
+	z = g2.doubleAndAdd(z, q)
+	z = g2.doubleN(z, 31)
+	z = g2.doubleAndAdd(z, q)
+	z = g2.doubleN(z, 16)
+
+	return g2.neg(z)
+}
+
+func (g2 G2) add(p, q *G2Affine) *G2Affine {
+	// compute λ = (q.y-p.y)/(q.x-p.x)
+	qypy := g2.Ext2.Sub(&q.Y, &p.Y)
+	qxpx := g2.Ext2.Sub(&q.X, &p.X)
+	λ := g2.Ext2.DivUnchecked(qypy, qxpx)
+
+	// xr = λ²-p.x-q.x
+	λλ := g2.Ext2.Square(λ)
+	qxpx = g2.Ext2.Add(&p.X, &q.X)
+	xr := g2.Ext2.Sub(λλ, qxpx)
+
+	// p.y = λ(p.x-r.x) - p.y
+	pxrx := g2.Ext2.Sub(&p.X, xr)
+	λpxrx := g2.Ext2.Mul(λ, pxrx)
+	yr := g2.Ext2.Sub(λpxrx, &p.Y)
+
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 G2) neg(p *G2Affine) *G2Affine {
+	xr := &p.X
+	yr := g2.Ext2.Neg(&p.Y)
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 G2) sub(p, q *G2Affine) *G2Affine {
+	qNeg := g2.neg(q)
+	return g2.add(p, qNeg)
+}
+
+func (g2 *G2) double(p *G2Affine) *G2Affine {
+	// compute λ = (3p.x²)/2*p.y
+	xx3a := g2.Square(&p.X)
+	xx3a = g2.MulByConstElement(xx3a, big.NewInt(3))
+	y2 := g2.Double(&p.Y)
+	λ := g2.DivUnchecked(xx3a, y2)
+
+	// xr = λ²-2p.x
+	x2 := g2.Double(&p.X)
+	λλ := g2.Square(λ)
+	xr := g2.Sub(λλ, x2)
+
+	// yr = λ(p-xr) - p.y
+	pxrx := g2.Sub(&p.X, xr)
+	λpxrx := g2.Mul(λ, pxrx)
+	yr := g2.Sub(λpxrx, &p.Y)
+
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 *G2) doubleN(p *G2Affine, n int) *G2Affine {
+	pn := p
+	for s := 0; s < n; s++ {
+		pn = g2.double(pn)
+	}
+	return pn
+}
+
+func (g2 G2) triple(p *G2Affine) *G2Affine {
+
+	// compute λ1 = (3p.x²)/2p.y
+	xx := g2.Square(&p.X)
+	xx = g2.MulByConstElement(xx, big.NewInt(3))
+	y2 := g2.Double(&p.Y)
+	λ1 := g2.DivUnchecked(xx, y2)
+
+	// xr = λ1²-2p.x
+	x2 := g2.MulByConstElement(&p.X, big.NewInt(2))
+	λ1λ1 := g2.Square(λ1)
+	x2 = g2.Sub(λ1λ1, x2)
+
+	// ommit y2 computation, and
+	// compute λ2 = 2p.y/(x2 − p.x) − λ1.
+	x1x2 := g2.Sub(&p.X, x2)
+	λ2 := g2.DivUnchecked(y2, x1x2)
+	λ2 = g2.Sub(λ2, λ1)
+
+	// xr = λ²-p.x-x2
+	λ2λ2 := g2.Square(λ2)
+	qxrx := g2.Add(x2, &p.X)
+	xr := g2.Sub(λ2λ2, qxrx)
+
+	// yr = λ(p.x-xr) - p.y
+	pxrx := g2.Sub(&p.X, xr)
+	λ2pxrx := g2.Mul(λ2, pxrx)
+	yr := g2.Sub(λ2pxrx, &p.Y)
+
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 G2) doubleAndAdd(p, q *G2Affine) *G2Affine {
+
+	// compute λ1 = (q.y-p.y)/(q.x-p.x)
+	yqyp := g2.Ext2.Sub(&q.Y, &p.Y)
+	xqxp := g2.Ext2.Sub(&q.X, &p.X)
+	λ1 := g2.Ext2.DivUnchecked(yqyp, xqxp)
+
+	// compute x2 = λ1²-p.x-q.x
+	λ1λ1 := g2.Ext2.Square(λ1)
+	xqxp = g2.Ext2.Add(&p.X, &q.X)
+	x2 := g2.Ext2.Sub(λ1λ1, xqxp)
+
+	// ommit y2 computation
+	// compute λ2 = -λ1-2*p.y/(x2-p.x)
+	ypyp := g2.Ext2.Add(&p.Y, &p.Y)
+	x2xp := g2.Ext2.Sub(x2, &p.X)
+	λ2 := g2.Ext2.DivUnchecked(ypyp, x2xp)
+	λ2 = g2.Ext2.Add(λ1, λ2)
+	λ2 = g2.Ext2.Neg(λ2)
+
+	// compute x3 =λ2²-p.x-x3
+	λ2λ2 := g2.Ext2.Square(λ2)
+	x3 := g2.Ext2.Sub(λ2λ2, &p.X)
+	x3 = g2.Ext2.Sub(x3, x2)
+
+	// compute y3 = λ2*(p.x - x3)-p.y
+	y3 := g2.Ext2.Sub(&p.X, x3)
+	y3 = g2.Ext2.Mul(λ2, y3)
+	y3 = g2.Ext2.Sub(y3, &p.Y)
+
+	return &G2Affine{
+		X: *x3,
+		Y: *y3,
+	}
+}
+
+// AssertIsEqual asserts that p and q are the same point.
+func (g2 *G2) AssertIsEqual(p, q *G2Affine) {
+	g2.Ext2.AssertIsEqual(&p.X, &q.X)
+	g2.Ext2.AssertIsEqual(&p.Y, &q.Y)
+}
diff --git a/std/algebra/emulated/sw_bls12381/g2_test.go b/std/algebra/emulated/sw_bls12381/g2_test.go
new file mode 100644
index 0000000000..76937d9623
--- /dev/null
+++ b/std/algebra/emulated/sw_bls12381/g2_test.go
@@ -0,0 +1,120 @@
+package sw_bls12381
+
+import (
+	"math/big"
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/test"
+)
+
+type addG2Circuit struct {
+	In1, In2 G2Affine
+	Res      G2Affine
+}
+
+func (c *addG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.add(&c.In1, &c.In2)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestAddG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines(assert)
+	_, in2 := randomG1G2Affines(assert)
+	var res bls12381.G2Affine
+	res.Add(&in1, &in2)
+	witness := addG2Circuit{
+		In1: NewG2Affine(in1),
+		In2: NewG2Affine(in2),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&addG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type doubleG2Circuit struct {
+	In1 G2Affine
+	Res G2Affine
+}
+
+func (c *doubleG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.double(&c.In1)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestDoubleG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines(assert)
+	var res bls12381.G2Affine
+	var in1Jac, resJac bls12381.G2Jac
+	in1Jac.FromAffine(&in1)
+	resJac.Double(&in1Jac)
+	res.FromJacobian(&resJac)
+	witness := doubleG2Circuit{
+		In1: NewG2Affine(in1),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&doubleG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type doubleAndAddG2Circuit struct {
+	In1, In2 G2Affine
+	Res      G2Affine
+}
+
+func (c *doubleAndAddG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.doubleAndAdd(&c.In1, &c.In2)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestDoubleAndAddG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines(assert)
+	_, in2 := randomG1G2Affines(assert)
+	var res bls12381.G2Affine
+	res.Double(&in1).
+		Add(&res, &in2)
+	witness := doubleAndAddG2Circuit{
+		In1: NewG2Affine(in1),
+		In2: NewG2Affine(in2),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&doubleAndAddG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type scalarMulG2BySeedCircuit struct {
+	In1 G2Affine
+	Res G2Affine
+}
+
+func (c *scalarMulG2BySeedCircuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.scalarMulBySeed(&c.In1)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestScalarMulG2BySeedTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines(assert)
+	var res bls12381.G2Affine
+	x0, _ := new(big.Int).SetString("15132376222941642752", 10)
+	res.ScalarMultiplication(&in1, x0).Neg(&res)
+	witness := scalarMulG2BySeedCircuit{
+		In1: NewG2Affine(in1),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&scalarMulG2BySeedCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go
index 788f6203b1..58038989c0 100644
--- a/std/algebra/emulated/sw_bls12381/pairing.go
+++ b/std/algebra/emulated/sw_bls12381/pairing.go
@@ -8,6 +8,7 @@ import (
 	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
 	"github.com/consensys/gnark/frontend"
 	"github.com/consensys/gnark/std/algebra/emulated/fields_bls12381"
+	"github.com/consensys/gnark/std/algebra/emulated/sw_emulated"
 	"github.com/consensys/gnark/std/math/emulated"
 )
 
@@ -15,6 +16,10 @@ type Pairing struct {
 	api frontend.API
 	*fields_bls12381.Ext12
 	curveF *emulated.Field[emulated.BLS12381Fp]
+	g2     *G2
+	g1     *G1
+	curve  *sw_emulated.Curve[emulated.BLS12381Fp, emulated.BLS12381Fr]
+	bTwist *fields_bls12381.E2
 }
 
 type GTEl = fields_bls12381.E12
@@ -57,10 +62,26 @@ func NewPairing(api frontend.API) (*Pairing, error) {
 	if err != nil {
 		return nil, fmt.Errorf("new base api: %w", err)
 	}
+	curve, err := sw_emulated.New[emulated.BLS12381Fp, emulated.BLS12381Fr](api, sw_emulated.GetBLS12381Params())
+	if err != nil {
+		return nil, fmt.Errorf("new curve: %w", err)
+	}
+	bTwist := fields_bls12381.E2{
+		A0: emulated.ValueOf[emulated.BLS12381Fp]("4"),
+		A1: emulated.ValueOf[emulated.BLS12381Fp]("4"),
+	}
+	g1, err := NewG1(api)
+	if err != nil {
+		return nil, fmt.Errorf("new G1 struct: %w", err)
+	}
 	return &Pairing{
 		api:    api,
 		Ext12:  fields_bls12381.NewExt12(api),
 		curveF: ba,
+		curve:  curve,
+		g1:     g1,
+		g2:     NewG2(api),
+		bTwist: &bTwist,
 	}, nil
 }
 
@@ -227,6 +248,57 @@ func (pr Pairing) AssertIsEqual(x, y *GTEl) {
 	pr.Ext12.AssertIsEqual(x, y)
 }
 
+func (pr Pairing) AssertIsOnCurve(P *G1Affine) {
+	pr.curve.AssertIsOnCurve(P)
+}
+
+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))
+
+	b := pr.Ext2.Select(selector, pr.Ext2.Zero(), pr.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
+	// TODO: add phi and scalarMulBySeedSquare to g1.go
+	// [x²]ϕ(P)
+	phiP := pr.g1.phi(P)
+	seedSquare := emulated.ValueOf[emulated.BLS12381Fr]("228988810152649578064853576960394133504")
+	// TODO: use addchain to construct a fixed-scalar ScalarMul
+	_P := pr.curve.ScalarMul(phiP, &seedSquare)
+	_P = pr.curve.Neg(_P)
+
+	// [r]Q == 0 <==>  P = -[x²]ϕ(P)
+	pr.curve.AssertIsEqual(_P, 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
+	// [x₀]Q
+	xQ := pr.g2.scalarMulBySeed(Q)
+	// ψ(Q)
+	psiQ := pr.g2.psi(Q)
+
+	// [r]Q == 0 <==>  ψ(Q) == [x₀]Q
+	pr.g2.AssertIsEqual(xQ, psiQ)
+}
+
 // loopCounter = seed in binary
 //
 //	seed=-15132376222941642752
diff --git a/std/algebra/emulated/sw_bls12381/pairing_test.go b/std/algebra/emulated/sw_bls12381/pairing_test.go
index ec1bf2baec..5a104df470 100644
--- a/std/algebra/emulated/sw_bls12381/pairing_test.go
+++ b/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)
+}
diff --git a/std/algebra/emulated/sw_bn254/doc_test.go b/std/algebra/emulated/sw_bn254/doc_test.go
index db095a0210..7d8ef6a6cd 100644
--- a/std/algebra/emulated/sw_bn254/doc_test.go
+++ b/std/algebra/emulated/sw_bn254/doc_test.go
@@ -23,6 +23,10 @@ func (c *PairCircuit) Define(api frontend.API) error {
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
+	// Pair method does not check that the points are in the proper groups.
+	pairing.AssertIsOnG1(&c.InG1)
+	pairing.AssertIsOnG2(&c.InG2)
+	// Compute the pairing
 	res, err := pairing.Pair([]*sw_bn254.G1Affine{&c.InG1}, []*sw_bn254.G2Affine{&c.InG2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
diff --git a/std/algebra/emulated/sw_bn254/g2.go b/std/algebra/emulated/sw_bn254/g2.go
index 9fd94325c5..71026f66f0 100644
--- a/std/algebra/emulated/sw_bn254/g2.go
+++ b/std/algebra/emulated/sw_bn254/g2.go
@@ -1,15 +1,42 @@
 package sw_bn254
 
 import (
+	"math/big"
+
 	"github.com/consensys/gnark-crypto/ecc/bn254"
+	"github.com/consensys/gnark/frontend"
 	"github.com/consensys/gnark/std/algebra/emulated/fields_bn254"
 	"github.com/consensys/gnark/std/math/emulated"
 )
 
+type G2 struct {
+	*fields_bn254.Ext2
+	w    *emulated.Element[emulated.BN254Fp]
+	u, v *fields_bn254.E2
+}
+
 type G2Affine struct {
 	X, Y fields_bn254.E2
 }
 
+func NewG2(api frontend.API) *G2 {
+	w := emulated.ValueOf[emulated.BN254Fp]("21888242871839275220042445260109153167277707414472061641714758635765020556616")
+	u := fields_bn254.E2{
+		A0: emulated.ValueOf[emulated.BN254Fp]("21575463638280843010398324269430826099269044274347216827212613867836435027261"),
+		A1: emulated.ValueOf[emulated.BN254Fp]("10307601595873709700152284273816112264069230130616436755625194854815875713954"),
+	}
+	v := fields_bn254.E2{
+		A0: emulated.ValueOf[emulated.BN254Fp]("2821565182194536844548159561693502659359617185244120367078079554186484126554"),
+		A1: emulated.ValueOf[emulated.BN254Fp]("3505843767911556378687030309984248845540243509899259641013678093033130930403"),
+	}
+	return &G2{
+		Ext2: fields_bn254.NewExt2(api),
+		w:    &w,
+		u:    &u,
+		v:    &v,
+	}
+}
+
 func NewG2Affine(v bn254.G2Affine) G2Affine {
 	return G2Affine{
 		X: fields_bn254.E2{
@@ -22,3 +49,168 @@ func NewG2Affine(v bn254.G2Affine) G2Affine {
 		},
 	}
 }
+
+func (g2 *G2) phi(q *G2Affine) *G2Affine {
+	x := g2.Ext2.MulByElement(&q.X, g2.w)
+
+	return &G2Affine{
+		X: *x,
+		Y: q.Y,
+	}
+}
+
+func (g2 *G2) psi(q *G2Affine) *G2Affine {
+	x := g2.Ext2.Conjugate(&q.X)
+	x = g2.Ext2.Mul(x, g2.u)
+	y := g2.Ext2.Conjugate(&q.Y)
+	y = g2.Ext2.Mul(y, g2.v)
+
+	return &G2Affine{
+		X: *x,
+		Y: *y,
+	}
+}
+
+func (g2 *G2) scalarMulBySeed(q *G2Affine) *G2Affine {
+	z := g2.double(q)
+	t0 := g2.add(q, z)
+	t2 := g2.add(q, t0)
+	t1 := g2.add(z, t2)
+	z = g2.doubleAndAdd(t1, t0)
+	t0 = g2.add(t0, z)
+	t2 = g2.add(t2, t0)
+	t1 = g2.add(t1, t2)
+	t0 = g2.add(t0, t1)
+	t1 = g2.add(t1, t0)
+	t0 = g2.add(t0, t1)
+	t2 = g2.add(t2, t0)
+	t1 = g2.doubleAndAdd(t2, t1)
+	t2 = g2.add(t2, t1)
+	z = g2.add(z, t2)
+	t2 = g2.add(t2, z)
+	z = g2.doubleAndAdd(t2, z)
+	t0 = g2.add(t0, z)
+	t1 = g2.add(t1, t0)
+	t3 := g2.double(t1)
+	t3 = g2.doubleAndAdd(t3, t1)
+	t2 = g2.add(t2, t3)
+	t1 = g2.add(t1, t2)
+	t2 = g2.add(t2, t1)
+	t2 = g2.doubleN(t2, 16)
+	t1 = g2.doubleAndAdd(t2, t1)
+	t1 = g2.doubleN(t1, 13)
+	t0 = g2.doubleAndAdd(t1, t0)
+	t0 = g2.doubleN(t0, 15)
+	z = g2.doubleAndAdd(t0, z)
+
+	return z
+}
+
+func (g2 G2) add(p, q *G2Affine) *G2Affine {
+	// compute λ = (q.y-p.y)/(q.x-p.x)
+	qypy := g2.Ext2.Sub(&q.Y, &p.Y)
+	qxpx := g2.Ext2.Sub(&q.X, &p.X)
+	λ := g2.Ext2.DivUnchecked(qypy, qxpx)
+
+	// xr = λ²-p.x-q.x
+	λλ := g2.Ext2.Square(λ)
+	qxpx = g2.Ext2.Add(&p.X, &q.X)
+	xr := g2.Ext2.Sub(λλ, qxpx)
+
+	// p.y = λ(p.x-r.x) - p.y
+	pxrx := g2.Ext2.Sub(&p.X, xr)
+	λpxrx := g2.Ext2.Mul(λ, pxrx)
+	yr := g2.Ext2.Sub(λpxrx, &p.Y)
+
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 G2) neg(p *G2Affine) *G2Affine {
+	xr := &p.X
+	yr := g2.Ext2.Neg(&p.Y)
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 G2) sub(p, q *G2Affine) *G2Affine {
+	qNeg := g2.neg(q)
+	return g2.add(p, qNeg)
+}
+
+func (g2 *G2) double(p *G2Affine) *G2Affine {
+	// compute λ = (3p.x²)/2*p.y
+	xx3a := g2.Square(&p.X)
+	xx3a = g2.MulByConstElement(xx3a, big.NewInt(3))
+	y2 := g2.Double(&p.Y)
+	λ := g2.DivUnchecked(xx3a, y2)
+
+	// xr = λ²-2p.x
+	x2 := g2.Double(&p.X)
+	λλ := g2.Square(λ)
+	xr := g2.Sub(λλ, x2)
+
+	// yr = λ(p-xr) - p.y
+	pxrx := g2.Sub(&p.X, xr)
+	λpxrx := g2.Mul(λ, pxrx)
+	yr := g2.Sub(λpxrx, &p.Y)
+
+	return &G2Affine{
+		X: *xr,
+		Y: *yr,
+	}
+}
+
+func (g2 *G2) doubleN(p *G2Affine, n int) *G2Affine {
+	pn := p
+	for s := 0; s < n; s++ {
+		pn = g2.double(pn)
+	}
+	return pn
+}
+
+func (g2 G2) doubleAndAdd(p, q *G2Affine) *G2Affine {
+
+	// compute λ1 = (q.y-p.y)/(q.x-p.x)
+	yqyp := g2.Ext2.Sub(&q.Y, &p.Y)
+	xqxp := g2.Ext2.Sub(&q.X, &p.X)
+	λ1 := g2.Ext2.DivUnchecked(yqyp, xqxp)
+
+	// compute x2 = λ1²-p.x-q.x
+	λ1λ1 := g2.Ext2.Square(λ1)
+	xqxp = g2.Ext2.Add(&p.X, &q.X)
+	x2 := g2.Ext2.Sub(λ1λ1, xqxp)
+
+	// ommit y2 computation
+	// compute λ2 = -λ1-2*p.y/(x2-p.x)
+	ypyp := g2.Ext2.Add(&p.Y, &p.Y)
+	x2xp := g2.Ext2.Sub(x2, &p.X)
+	λ2 := g2.Ext2.DivUnchecked(ypyp, x2xp)
+	λ2 = g2.Ext2.Add(λ1, λ2)
+	λ2 = g2.Ext2.Neg(λ2)
+
+	// compute x3 =λ2²-p.x-x3
+	λ2λ2 := g2.Ext2.Square(λ2)
+	x3 := g2.Ext2.Sub(λ2λ2, &p.X)
+	x3 = g2.Ext2.Sub(x3, x2)
+
+	// compute y3 = λ2*(p.x - x3)-p.y
+	y3 := g2.Ext2.Sub(&p.X, x3)
+	y3 = g2.Ext2.Mul(λ2, y3)
+	y3 = g2.Ext2.Sub(y3, &p.Y)
+
+	return &G2Affine{
+		X: *x3,
+		Y: *y3,
+	}
+}
+
+// AssertIsEqual asserts that p and q are the same point.
+func (g2 *G2) AssertIsEqual(p, q *G2Affine) {
+	g2.Ext2.AssertIsEqual(&p.X, &q.X)
+	g2.Ext2.AssertIsEqual(&p.Y, &q.Y)
+}
diff --git a/std/algebra/emulated/sw_bn254/g2_test.go b/std/algebra/emulated/sw_bn254/g2_test.go
new file mode 100644
index 0000000000..3e61f05cc4
--- /dev/null
+++ b/std/algebra/emulated/sw_bn254/g2_test.go
@@ -0,0 +1,144 @@
+package sw_bn254
+
+import (
+	"math/big"
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	"github.com/consensys/gnark-crypto/ecc/bn254"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/test"
+)
+
+type addG2Circuit struct {
+	In1, In2 G2Affine
+	Res      G2Affine
+}
+
+func (c *addG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.add(&c.In1, &c.In2)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestAddG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines()
+	_, in2 := randomG1G2Affines()
+	var res bn254.G2Affine
+	res.Add(&in1, &in2)
+	witness := addG2Circuit{
+		In1: NewG2Affine(in1),
+		In2: NewG2Affine(in2),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&addG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type doubleG2Circuit struct {
+	In1 G2Affine
+	Res G2Affine
+}
+
+func (c *doubleG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.double(&c.In1)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestDoubleG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines()
+	var res bn254.G2Affine
+	var in1Jac, resJac bn254.G2Jac
+	in1Jac.FromAffine(&in1)
+	resJac.Double(&in1Jac)
+	res.FromJacobian(&resJac)
+	witness := doubleG2Circuit{
+		In1: NewG2Affine(in1),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&doubleG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type doubleAndAddG2Circuit struct {
+	In1, In2 G2Affine
+	Res      G2Affine
+}
+
+func (c *doubleAndAddG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.doubleAndAdd(&c.In1, &c.In2)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestDoubleAndAddG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines()
+	_, in2 := randomG1G2Affines()
+	var res bn254.G2Affine
+	res.Double(&in1).
+		Add(&res, &in2)
+	witness := doubleAndAddG2Circuit{
+		In1: NewG2Affine(in1),
+		In2: NewG2Affine(in2),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&doubleAndAddG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type scalarMulG2BySeedCircuit struct {
+	In1 G2Affine
+	Res G2Affine
+}
+
+func (c *scalarMulG2BySeedCircuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res := g2.scalarMulBySeed(&c.In1)
+	g2.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestScalarMulG2BySeedTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines()
+	var res bn254.G2Affine
+	x0, _ := new(big.Int).SetString("4965661367192848881", 10)
+	res.ScalarMultiplication(&in1, x0)
+	witness := scalarMulG2BySeedCircuit{
+		In1: NewG2Affine(in1),
+		Res: NewG2Affine(res),
+	}
+	err := test.IsSolved(&scalarMulG2BySeedCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type endomorphismG2Circuit struct {
+	In1 G2Affine
+}
+
+func (c *endomorphismG2Circuit) Define(api frontend.API) error {
+	g2 := NewG2(api)
+	res1 := g2.phi(&c.In1)
+	res1 = g2.neg(res1)
+	res2 := g2.psi(&c.In1)
+	res2 = g2.psi(res2)
+	g2.AssertIsEqual(res1, res2)
+	return nil
+}
+
+func TestEndomorphismG2TestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, in1 := randomG1G2Affines()
+	witness := endomorphismG2Circuit{
+		In1: NewG2Affine(in1),
+	}
+	err := test.IsSolved(&endomorphismG2Circuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
diff --git a/std/algebra/emulated/sw_bn254/pairing.go b/std/algebra/emulated/sw_bn254/pairing.go
index 88245f1714..5ed4e08f81 100644
--- a/std/algebra/emulated/sw_bn254/pairing.go
+++ b/std/algebra/emulated/sw_bn254/pairing.go
@@ -8,6 +8,7 @@ import (
 	"github.com/consensys/gnark-crypto/ecc/bn254"
 	"github.com/consensys/gnark/frontend"
 	"github.com/consensys/gnark/std/algebra/emulated/fields_bn254"
+	"github.com/consensys/gnark/std/algebra/emulated/sw_emulated"
 	"github.com/consensys/gnark/std/math/emulated"
 )
 
@@ -15,6 +16,9 @@ type Pairing struct {
 	api frontend.API
 	*fields_bn254.Ext12
 	curveF *emulated.Field[emulated.BN254Fp]
+	curve  *sw_emulated.Curve[emulated.BN254Fp, emulated.BN254Fr]
+	g2     *G2
+	bTwist *fields_bn254.E2
 }
 
 type GTEl = fields_bn254.E12
@@ -57,10 +61,21 @@ func NewPairing(api frontend.API) (*Pairing, error) {
 	if err != nil {
 		return nil, fmt.Errorf("new base api: %w", err)
 	}
+	curve, err := sw_emulated.New[emulated.BN254Fp, emulated.BN254Fr](api, sw_emulated.GetBN254Params())
+	if err != nil {
+		return nil, fmt.Errorf("new curve: %w", err)
+	}
+	bTwist := fields_bn254.E2{
+		A0: emulated.ValueOf[emulated.BN254Fp]("19485874751759354771024239261021720505790618469301721065564631296452457478373"),
+		A1: emulated.ValueOf[emulated.BN254Fp]("266929791119991161246907387137283842545076965332900288569378510910307636690"),
+	}
 	return &Pairing{
 		api:    api,
 		Ext12:  fields_bn254.NewExt12(api),
 		curveF: ba,
+		curve:  curve,
+		g2:     NewG2(api),
+		bTwist: &bTwist,
 	}, nil
 }
 
@@ -202,7 +217,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 +230,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 +247,56 @@ func (pr Pairing) AssertIsEqual(x, y *GTEl) {
 	pr.Ext12.AssertIsEqual(x, y)
 }
 
+func (pr Pairing) AssertIsOnCurve(P *G1Affine) {
+	pr.curve.AssertIsOnCurve(P)
+}
+
+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))
+	b := pr.Ext2.Select(selector, pr.Ext2.Zero(), pr.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
+
+	// [x₀]Q
+	xQ := pr.g2.scalarMulBySeed(Q)
+	// ψ([x₀]Q)
+	psixQ := pr.g2.psi(xQ)
+	// ψ²([x₀]Q) = -ϕ([x₀]Q)
+	psi2xQ := pr.g2.phi(xQ)
+	// ψ³([2x₀]Q)
+	psi3xxQ := pr.g2.double(psi2xQ)
+	psi3xxQ = pr.g2.psi(psi3xxQ)
+
+	// _Q = ψ³([2x₀]Q) - ψ²([x₀]Q) - ψ([x₀]Q) - [x₀]Q
+	_Q := pr.g2.sub(psi2xQ, psi3xxQ)
+	_Q = pr.g2.sub(_Q, psixQ)
+	_Q = pr.g2.sub(_Q, xQ)
+
+	// [r]Q == 0 <==>  _Q == Q
+	pr.g2.AssertIsEqual(Q, _Q)
+}
+
 // loopCounter = 6x₀+2 = 29793968203157093288
 //
 // in 2-NAF
@@ -436,15 +501,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)
diff --git a/std/algebra/emulated/sw_bn254/pairing_test.go b/std/algebra/emulated/sw_bn254/pairing_test.go
index 90234b63d2..66a77219fa 100644
--- a/std/algebra/emulated/sw_bn254/pairing_test.go
+++ b/std/algebra/emulated/sw_bn254/pairing_test.go
@@ -74,6 +74,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)
@@ -109,6 +111,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)
@@ -205,6 +211,32 @@ func TestFinalExponentiationSafeCircuit(t *testing.T) {
 	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()
+	witness := GroupMembershipCircuit{
+		InG1: NewG1Affine(p),
+		InG2: NewG2Affine(q),
+	}
+	err := test.IsSolved(&GroupMembershipCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
 func BenchmarkPairing(b *testing.B) {
 
 	p1, q1 := randomG1G2Affines()
diff --git a/std/algebra/emulated/sw_emulated/params.go b/std/algebra/emulated/sw_emulated/params.go
index b9f00ec707..bb513b6e41 100644
--- a/std/algebra/emulated/sw_emulated/params.go
+++ b/std/algebra/emulated/sw_emulated/params.go
@@ -3,6 +3,7 @@ package sw_emulated
 import (
 	"math/big"
 
+	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
 	"github.com/consensys/gnark-crypto/ecc/bn254"
 	"github.com/consensys/gnark-crypto/ecc/secp256k1"
 	"github.com/consensys/gnark/std/math/emulated"
@@ -50,6 +51,20 @@ func GetBN254Params() CurveParams {
 	}
 }
 
+// GetBLS12381Params returns the curve parameters for the curve BLS12-381.
+// When initialising new curve, use the base field [emulated.BLS12381Fp] and scalar
+// field [emulated.BLS12381Fr].
+func GetBLS12381Params() CurveParams {
+	_, _, g1aff, _ := bls12381.Generators()
+	return CurveParams{
+		A:  big.NewInt(0),
+		B:  big.NewInt(4),
+		Gx: g1aff.X.BigInt(new(big.Int)),
+		Gy: g1aff.Y.BigInt(new(big.Int)),
+		Gm: computeBLS12381Table(),
+	}
+}
+
 // GetCurveParams returns suitable curve parameters given the parametric type Base as base field.
 func GetCurveParams[Base emulated.FieldParams]() CurveParams {
 	var t Base
@@ -58,6 +73,8 @@ func GetCurveParams[Base emulated.FieldParams]() CurveParams {
 		return secp256k1Params
 	case "30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47":
 		return bn254Params
+	case "1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab":
+		return bls12381Params
 	default:
 		panic("no stored parameters")
 	}
@@ -66,9 +83,11 @@ func GetCurveParams[Base emulated.FieldParams]() CurveParams {
 var (
 	secp256k1Params CurveParams
 	bn254Params     CurveParams
+	bls12381Params  CurveParams
 )
 
 func init() {
 	secp256k1Params = GetSecp256k1Params()
 	bn254Params = GetBN254Params()
+	bls12381Params = GetBLS12381Params()
 }
diff --git a/std/algebra/emulated/sw_emulated/params_compute.go b/std/algebra/emulated/sw_emulated/params_compute.go
index fbdabea300..58d3e26df8 100644
--- a/std/algebra/emulated/sw_emulated/params_compute.go
+++ b/std/algebra/emulated/sw_emulated/params_compute.go
@@ -3,6 +3,7 @@ package sw_emulated
 import (
 	"math/big"
 
+	bls12381 "github.com/consensys/gnark-crypto/ecc/bls12-381"
 	"github.com/consensys/gnark-crypto/ecc/bn254"
 	"github.com/consensys/gnark-crypto/ecc/secp256k1"
 )
@@ -58,3 +59,29 @@ func computeBN254Table() [][2]*big.Int {
 	}
 	return table
 }
+
+func computeBLS12381Table() [][2]*big.Int {
+	Gjac, _, _, _ := bls12381.Generators()
+	table := make([][2]*big.Int, 256)
+	tmp := new(bls12381.G1Jac).Set(&Gjac)
+	aff := new(bls12381.G1Affine)
+	jac := new(bls12381.G1Jac)
+	for i := 1; i < 256; i++ {
+		tmp = tmp.Double(tmp)
+		switch i {
+		case 1, 2:
+			jac.Set(tmp).AddAssign(&Gjac)
+			aff.FromJacobian(jac)
+			table[i-1] = [2]*big.Int{aff.X.BigInt(new(big.Int)), aff.Y.BigInt(new(big.Int))}
+		case 3:
+			jac.Set(tmp).SubAssign(&Gjac)
+			aff.FromJacobian(jac)
+			table[i-1] = [2]*big.Int{aff.X.BigInt(new(big.Int)), aff.Y.BigInt(new(big.Int))}
+			fallthrough
+		default:
+			aff.FromJacobian(tmp)
+			table[i] = [2]*big.Int{aff.X.BigInt(new(big.Int)), aff.Y.BigInt(new(big.Int))}
+		}
+	}
+	return table
+}
diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go
index c5a56dabd4..7f4dc1dda5 100644
--- a/std/algebra/emulated/sw_emulated/point.go
+++ b/std/algebra/emulated/sw_emulated/point.go
@@ -38,6 +38,7 @@ func New[Base, Scalars emulated.FieldParams](api frontend.API, params CurveParam
 		},
 		gm:   emuGm,
 		a:    emulated.ValueOf[Base](params.A),
+		b:    emulated.ValueOf[Base](params.B),
 		addA: params.A.Cmp(big.NewInt(0)) != 0,
 	}, nil
 }
@@ -60,6 +61,7 @@ type Curve[Base, Scalars emulated.FieldParams] struct {
 	gm []AffinePoint[Base]
 
 	a    emulated.Element[Base]
+	b    emulated.Element[Base]
 	addA bool
 }
 
@@ -125,6 +127,24 @@ func (c *Curve[B, S]) add(p, q *AffinePoint[B]) *AffinePoint[B] {
 	}
 }
 
+// AssertIsOnCurve asserts if p belongs to the curve. It doesn't modify p.
+func (c *Curve[B, S]) AssertIsOnCurve(p *AffinePoint[B]) {
+	// (X,Y) ∈ {Y² == X³ + aX + b} U (0,0)
+
+	// if p=(0,0) we assign b=0 and continue
+	selector := c.api.And(c.baseApi.IsZero(&p.X), c.baseApi.IsZero(&p.Y))
+	b := c.baseApi.Select(selector, c.baseApi.Zero(), &c.b)
+
+	left := c.baseApi.Mul(&p.Y, &p.Y)
+	right := c.baseApi.Mul(&p.X, c.baseApi.Mul(&p.X, &p.X))
+	right = c.baseApi.Add(right, b)
+	if c.addA {
+		ax := c.baseApi.Mul(&c.a, &p.X)
+		right = c.baseApi.Add(right, ax)
+	}
+	c.baseApi.AssertIsEqual(left, right)
+}
+
 // AddUnified adds p and q and returns it. It doesn't modify p nor q.
 //
 // ✅ p can be equal to q, and either or both can be (0,0).
@@ -337,6 +357,7 @@ func (c *Curve[B, S]) Lookup2(b0, b1 frontend.Variable, i0, i1, i2, i3 *AffinePo
 }
 
 // ScalarMul computes s * p and returns it. It doesn't modify p nor s.
+// This function doesn't check that the p is on the curve. See AssertIsOnCurve.
 //
 // ✅ p can can be (0,0) and s can be 0.
 // (0,0) is not on the curve but we conventionally take it as the
diff --git a/std/algebra/emulated/sw_emulated/point_test.go b/std/algebra/emulated/sw_emulated/point_test.go
index 0b0f90fe9e..0543c8bb91 100644
--- a/std/algebra/emulated/sw_emulated/point_test.go
+++ b/std/algebra/emulated/sw_emulated/point_test.go
@@ -5,6 +5,8 @@ import (
 	"testing"
 
 	"github.com/consensys/gnark-crypto/ecc"
+	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"
@@ -393,6 +395,28 @@ func TestScalarMulBase2(t *testing.T) {
 	assert.NoError(err)
 }
 
+func TestScalarMulBase3(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, _, g, _ := bls12381.Generators()
+	var r fr_bn.Element
+	_, _ = r.SetRandom()
+	s := new(big.Int)
+	r.BigInt(s)
+	var S bls12381.G1Affine
+	S.ScalarMultiplication(&g, s)
+
+	circuit := ScalarMulBaseTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{}
+	witness := ScalarMulBaseTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{
+		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)
+}
+
 type ScalarMulTest[T, S emulated.FieldParams] struct {
 	P, Q AffinePoint[T]
 	S    emulated.Element[S]
@@ -518,3 +542,138 @@ func TestScalarMulEdgeCasesEdgeCases(t *testing.T) {
 	err = test.IsSolved(&circuit, &witness2, testCurve.ScalarField())
 	assert.NoError(err)
 }
+
+func TestScalarMul3(t *testing.T) {
+	assert := test.NewAssert(t)
+	var r fr_bls381.Element
+	_, _ = r.SetRandom()
+	s := new(big.Int)
+	r.BigInt(s)
+	var res bls12381.G1Affine
+	_, _, gen, _ := bls12381.Generators()
+	res.ScalarMultiplication(&gen, s)
+
+	circuit := ScalarMulTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{}
+	witness := ScalarMulTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{
+		S: emulated.ValueOf[emulated.BLS12381Fr](s),
+		P: AffinePoint[emulated.BLS12381Fp]{
+			X: emulated.ValueOf[emulated.BLS12381Fp](gen.X),
+			Y: emulated.ValueOf[emulated.BLS12381Fp](gen.Y),
+		},
+		Q: AffinePoint[emulated.BLS12381Fp]{
+			X: emulated.ValueOf[emulated.BLS12381Fp](res.X),
+			Y: emulated.ValueOf[emulated.BLS12381Fp](res.Y),
+		},
+	}
+	err := test.IsSolved(&circuit, &witness, testCurve.ScalarField())
+	assert.NoError(err)
+}
+
+type IsOnCurveTest[T, S emulated.FieldParams] struct {
+	Q AffinePoint[T]
+}
+
+func (c *IsOnCurveTest[T, S]) Define(api frontend.API) error {
+	cr, err := New[T, S](api, GetCurveParams[T]())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsOnCurve(&c.Q)
+	return nil
+}
+
+func TestIsOnCurve(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, g := secp256k1.Generators()
+	var r fr_secp.Element
+	_, _ = r.SetRandom()
+	s := new(big.Int)
+	r.BigInt(s)
+	var Q, infinity secp256k1.G1Affine
+	Q.ScalarMultiplication(&g, s)
+
+	// Q=[s]G is on curve
+	circuit := IsOnCurveTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{}
+	witness1 := IsOnCurveTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{
+		Q: AffinePoint[emulated.Secp256k1Fp]{
+			X: emulated.ValueOf[emulated.Secp256k1Fp](Q.X),
+			Y: emulated.ValueOf[emulated.Secp256k1Fp](Q.Y),
+		},
+	}
+	err := test.IsSolved(&circuit, &witness1, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// (0,0) is on curve
+	witness2 := IsOnCurveTest[emulated.Secp256k1Fp, emulated.Secp256k1Fr]{
+		Q: AffinePoint[emulated.Secp256k1Fp]{
+			X: emulated.ValueOf[emulated.Secp256k1Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.Secp256k1Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness2, testCurve.ScalarField())
+	assert.NoError(err)
+}
+
+func TestIsOnCurve2(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, _, g, _ := bn254.Generators()
+	var r fr_secp.Element
+	_, _ = r.SetRandom()
+	s := new(big.Int)
+	r.BigInt(s)
+	var Q, infinity bn254.G1Affine
+	Q.ScalarMultiplication(&g, s)
+
+	// Q=[s]G is on curve
+	circuit := IsOnCurveTest[emulated.BN254Fp, emulated.BN254Fr]{}
+	witness1 := IsOnCurveTest[emulated.BN254Fp, emulated.BN254Fr]{
+		Q: AffinePoint[emulated.BN254Fp]{
+			X: emulated.ValueOf[emulated.BN254Fp](Q.X),
+			Y: emulated.ValueOf[emulated.BN254Fp](Q.Y),
+		},
+	}
+	err := test.IsSolved(&circuit, &witness1, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// (0,0) is on curve
+	witness2 := IsOnCurveTest[emulated.BN254Fp, emulated.BN254Fr]{
+		Q: AffinePoint[emulated.BN254Fp]{
+			X: emulated.ValueOf[emulated.BN254Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BN254Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness2, testCurve.ScalarField())
+	assert.NoError(err)
+}
+
+func TestIsOnCurve3(t *testing.T) {
+	assert := test.NewAssert(t)
+	_, _, g, _ := bls12381.Generators()
+	var r fr_secp.Element
+	_, _ = r.SetRandom()
+	s := new(big.Int)
+	r.BigInt(s)
+	var Q, infinity bls12381.G1Affine
+	Q.ScalarMultiplication(&g, s)
+
+	// Q=[s]G is on curve
+	circuit := IsOnCurveTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{}
+	witness1 := IsOnCurveTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{
+		Q: AffinePoint[emulated.BLS12381Fp]{
+			X: emulated.ValueOf[emulated.BLS12381Fp](Q.X),
+			Y: emulated.ValueOf[emulated.BLS12381Fp](Q.Y),
+		},
+	}
+	err := test.IsSolved(&circuit, &witness1, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// (0,0) is on curve
+	witness2 := IsOnCurveTest[emulated.BLS12381Fp, emulated.BLS12381Fr]{
+		Q: AffinePoint[emulated.BLS12381Fp]{
+			X: emulated.ValueOf[emulated.BLS12381Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BLS12381Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness2, testCurve.ScalarField())
+	assert.NoError(err)
+}
diff --git a/std/evmprecompiles/06-bnadd.go b/std/evmprecompiles/06-bnadd.go
index 0f819d9833..ff6c397fc9 100644
--- a/std/evmprecompiles/06-bnadd.go
+++ b/std/evmprecompiles/06-bnadd.go
@@ -14,7 +14,8 @@ func ECAdd(api frontend.API, P, Q *sw_emulated.AffinePoint[emulated.BN254Fp]) *s
 	if err != nil {
 		panic(err)
 	}
-	// We use AddUnified because P can be equal to Q
+	// Check that P and Q are on the curve (done in the zkEVM ⚠️ )
+	// We use AddUnified because P can be equal to Q, -Q and either or both can be (0,0)
 	res := curve.AddUnified(P, Q)
 	return res
 }
diff --git a/std/evmprecompiles/07-bnmul.go b/std/evmprecompiles/07-bnmul.go
index 1b6d2d6fa8..eb1c02ccda 100644
--- a/std/evmprecompiles/07-bnmul.go
+++ b/std/evmprecompiles/07-bnmul.go
@@ -14,6 +14,7 @@ func ECMul(api frontend.API, P *sw_emulated.AffinePoint[emulated.BN254Fp], u *em
 	if err != nil {
 		panic(err)
 	}
+	// Check that P is on the curve (done in the zkEVM ⚠️ )
 	res := curve.ScalarMul(P, u)
 	return res
 }
diff --git a/std/evmprecompiles/08-bnpairing.go b/std/evmprecompiles/08-bnpairing.go
index cab5c3ad2b..01dbbab99c 100644
--- a/std/evmprecompiles/08-bnpairing.go
+++ b/std/evmprecompiles/08-bnpairing.go
@@ -13,6 +13,13 @@ func ECPair(api frontend.API, P []*sw_bn254.G1Affine, Q []*sw_bn254.G2Affine) {
 	if err != nil {
 		panic(err)
 	}
+	// 1- Check that Pᵢ are on G1 (done in the zkEVM ⚠️ )
+	// 2- Check that Qᵢ are on G2
+	for i := 0; i < len(Q); i++ {
+		pair.AssertIsOnG2(Q[i])
+	}
+
+	// 3- Check that ∏ᵢ e(Pᵢ, Qᵢ) == 1
 	if err := pair.PairingCheck(P, Q); err != nil {
 		panic(err)
 	}
diff --git a/std/math/emulated/params.go b/std/math/emulated/params.go
index eda8ebddc5..b40b2f2493 100644
--- a/std/math/emulated/params.go
+++ b/std/math/emulated/params.go
@@ -97,3 +97,14 @@ func (fp BLS12381Fp) NbLimbs() uint     { return 6 }
 func (fp BLS12381Fp) BitsPerLimb() uint { return 64 }
 func (fp BLS12381Fp) IsPrime() bool     { return true }
 func (fp BLS12381Fp) Modulus() *big.Int { return ecc.BLS12_381.BaseField() }
+
+// BLS12381Fr provide type parametrization for emulated field on 4 limb of width
+// 64bits for modulus
+// 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001.
+// This is the scalar field of the BLS12-381 curve.
+type BLS12381Fr struct{}
+
+func (fp BLS12381Fr) NbLimbs() uint     { return 4 }
+func (fp BLS12381Fr) BitsPerLimb() uint { return 64 }
+func (fp BLS12381Fr) IsPrime() bool     { return true }
+func (fp BLS12381Fr) Modulus() *big.Int { return ecc.BLS12_381.ScalarField() }