diff --git a/go.mod b/go.mod
index 7ed857bba6..28c4e209fb 100644
--- a/go.mod
+++ b/go.mod
@@ -6,7 +6,7 @@ require (
 	github.com/bits-and-blooms/bitset v1.8.0
 	github.com/blang/semver/v4 v4.0.0
 	github.com/consensys/bavard v0.1.13
-	github.com/consensys/gnark-crypto v0.12.2-0.20231013160410-1f65e75b6dfb
+	github.com/consensys/gnark-crypto v0.12.2-0.20231017134050-6652a8b98254
 	github.com/fxamacker/cbor/v2 v2.5.0
 	github.com/google/go-cmp v0.5.9
 	github.com/google/pprof v0.0.0-20230817174616-7a8ec2ada47b
diff --git a/go.sum b/go.sum
index 5cd759d52b..cfb98a25a5 100644
--- a/go.sum
+++ b/go.sum
@@ -4,10 +4,8 @@ github.com/blang/semver/v4 v4.0.0 h1:1PFHFE6yCCTv8C1TeyNNarDzntLi7wMI5i/pzqYIsAM
 github.com/blang/semver/v4 v4.0.0/go.mod h1:IbckMUScFkM3pff0VJDNKRiT6TG/YpiHIM2yvyW5YoQ=
 github.com/consensys/bavard v0.1.13 h1:oLhMLOFGTLdlda/kma4VOJazblc7IM5y5QPd2A/YjhQ=
 github.com/consensys/bavard v0.1.13/go.mod h1:9ItSMtA/dXMAiL7BG6bqW2m3NdSEObYWoH223nGHukI=
-github.com/consensys/gnark-crypto v0.12.2-0.20231012161402-206544105834 h1:o+Q1/PSZfNkoUAl/Gf5N/u+H6LkOzR3gF0mS8nRSWYM=
-github.com/consensys/gnark-crypto v0.12.2-0.20231012161402-206544105834/go.mod h1:v2Gy7L/4ZRosZ7Ivs+9SfUDr0f5UlG+EM5t7MPHiLuY=
-github.com/consensys/gnark-crypto v0.12.2-0.20231013160410-1f65e75b6dfb h1:f0BMgIjhZy4lSRHCXFbQst85f5agZAjtDMixQqBWNpc=
-github.com/consensys/gnark-crypto v0.12.2-0.20231013160410-1f65e75b6dfb/go.mod h1:v2Gy7L/4ZRosZ7Ivs+9SfUDr0f5UlG+EM5t7MPHiLuY=
+github.com/consensys/gnark-crypto v0.12.2-0.20231017134050-6652a8b98254 h1:21iGClQpXhJ0PA6lRBpRbpmJFtfdV7R9QMCsQJ3A9mA=
+github.com/consensys/gnark-crypto v0.12.2-0.20231017134050-6652a8b98254/go.mod h1:v2Gy7L/4ZRosZ7Ivs+9SfUDr0f5UlG+EM5t7MPHiLuY=
 github.com/coreos/go-systemd/v22 v22.5.0/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.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c=
diff --git a/std/algebra/emulated/fields_bw6761/doc.go b/std/algebra/emulated/fields_bw6761/doc.go
new file mode 100644
index 0000000000..29858f7424
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/doc.go
@@ -0,0 +1,6 @@
+// Package fields_bw6761 implements the fields arithmetic of the Fp6 tower
+// used to compute the pairing over the BW6-761 curve.
+//
+//	𝔽p³[u] = 𝔽p/u³+4
+//	𝔽p⁶[v] = 𝔽p²/v²-u
+package fields_bw6761
diff --git a/std/algebra/emulated/fields_bw6761/e3.go b/std/algebra/emulated/fields_bw6761/e3.go
new file mode 100644
index 0000000000..53361f6ab4
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/e3.go
@@ -0,0 +1,401 @@
+package fields_bw6761
+
+import (
+	"math/big"
+
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+type curveF = emulated.Field[emulated.BW6761Fp]
+type baseEl = emulated.Element[emulated.BW6761Fp]
+
+type E3 struct {
+	A0, A1, A2 baseEl
+}
+
+type Ext3 struct {
+	api frontend.API
+	fp  *curveF
+}
+
+func NewExt3(api frontend.API) *Ext3 {
+	fp, err := emulated.NewField[emulated.BW6761Fp](api)
+	if err != nil {
+		panic(err)
+	}
+	return &Ext3{
+		api: api,
+		fp:  fp,
+	}
+}
+
+func (e Ext3) Reduce(x *E3) *E3 {
+	var z E3
+	z.A0 = *e.fp.Reduce(&x.A0)
+	z.A1 = *e.fp.Reduce(&x.A1)
+	z.A2 = *e.fp.Reduce(&x.A2)
+	return &z
+}
+
+func (e Ext3) Zero() *E3 {
+	zero := e.fp.Zero()
+	return &E3{
+		A0: *zero,
+		A1: *zero,
+		A2: *zero,
+	}
+}
+
+func (e Ext3) One() *E3 {
+	one := e.fp.One()
+	zero := e.fp.Zero()
+	return &E3{
+		A0: *one,
+		A1: *zero,
+		A2: *zero,
+	}
+}
+
+func (e Ext3) Neg(x *E3) *E3 {
+	a0 := e.fp.Neg(&x.A0)
+	a1 := e.fp.Neg(&x.A1)
+	a2 := e.fp.Neg(&x.A2)
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+func (e Ext3) Add(x, y *E3) *E3 {
+	a0 := e.fp.Add(&x.A0, &y.A0)
+	a1 := e.fp.Add(&x.A1, &y.A1)
+	a2 := e.fp.Add(&x.A2, &y.A2)
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+func (e Ext3) Sub(x, y *E3) *E3 {
+	a0 := e.fp.Sub(&x.A0, &y.A0)
+	a1 := e.fp.Sub(&x.A1, &y.A1)
+	a2 := e.fp.Sub(&x.A2, &y.A2)
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+func (e Ext3) Double(x *E3) *E3 {
+	two := big.NewInt(2)
+	a0 := e.fp.MulConst(&x.A0, two)
+	a1 := e.fp.MulConst(&x.A1, two)
+	a2 := e.fp.MulConst(&x.A2, two)
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+func mulFpByNonResidue(fp *curveF, x *baseEl) *baseEl {
+
+	z := fp.Neg(x)
+	z = fp.MulConst(z, big.NewInt(4))
+	return z
+}
+
+func (e Ext3) Conjugate(x *E3) *E3 {
+	a1 := e.fp.Neg(&x.A1)
+	return &E3{
+		A0: x.A0,
+		A1: *a1,
+		A2: x.A2,
+	}
+}
+
+func (e Ext3) MulByElement(x *E3, y *baseEl) *E3 {
+	a0 := e.fp.Mul(&x.A0, y)
+	a1 := e.fp.Mul(&x.A1, y)
+	a2 := e.fp.Mul(&x.A2, y)
+	z := &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+	return z
+}
+
+func (e Ext3) MulByConstElement(x *E3, y *big.Int) *E3 {
+	a0 := e.fp.MulConst(&x.A0, y)
+	a1 := e.fp.MulConst(&x.A1, y)
+	a2 := e.fp.MulConst(&x.A2, y)
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+// MulBy01 multiplication by sparse element (c0,c1,0)
+func (e Ext3) MulBy01(z *E3, c0, c1 *baseEl) *E3 {
+
+	a := e.fp.Mul(&z.A0, c0)
+	b := e.fp.Mul(&z.A1, c1)
+
+	tmp := e.fp.Add(&z.A1, &z.A2)
+	t0 := e.fp.Mul(c1, tmp)
+	t0 = e.fp.Sub(t0, b)
+	t0 = mulFpByNonResidue(e.fp, t0)
+	t0 = e.fp.Add(t0, a)
+
+	tmp = e.fp.Add(&z.A0, &z.A2)
+	t2 := e.fp.Mul(c0, tmp)
+	t2 = e.fp.Sub(t2, a)
+	t2 = e.fp.Add(t2, b)
+
+	t1 := e.fp.Add(c0, c1)
+	tmp = e.fp.Add(&z.A0, &z.A1)
+	t1 = e.fp.Mul(t1, tmp)
+	t1 = e.fp.Sub(t1, a)
+	t1 = e.fp.Sub(t1, b)
+
+	return &E3{
+		A0: *t0,
+		A1: *t1,
+		A2: *t2,
+	}
+}
+
+// MulBy1 multiplication of E6 by sparse element (0, c1, 0)
+func (e Ext3) MulBy1(z *E3, c1 *baseEl) *E3 {
+
+	b := e.fp.Mul(&z.A1, c1)
+
+	tmp := e.fp.Add(&z.A1, &z.A2)
+	t0 := e.fp.Mul(c1, tmp)
+	t0 = e.fp.Sub(t0, b)
+	t0 = mulFpByNonResidue(e.fp, t0)
+
+	tmp = e.fp.Add(&z.A0, &z.A1)
+	t1 := e.fp.Mul(c1, tmp)
+	t1 = e.fp.Sub(t1, b)
+
+	return &E3{
+		A0: *t0,
+		A1: *t1,
+		A2: *b,
+	}
+}
+
+// MulBy12 multiplication by sparse element (0,b1,b2)
+func (e Ext3) MulBy12(x *E3, b1, b2 *baseEl) *E3 {
+	t1 := e.fp.Mul(&x.A1, b1)
+	t2 := e.fp.Mul(&x.A2, b2)
+	c0 := e.fp.Add(&x.A1, &x.A2)
+	tmp := e.fp.Add(b1, b2)
+	c0 = e.fp.Mul(c0, tmp)
+	c0 = e.fp.Sub(c0, t1)
+	c0 = e.fp.Sub(c0, t2)
+	c0 = mulFpByNonResidue(e.fp, c0)
+	c1 := e.fp.Add(&x.A0, &x.A1)
+	c1 = e.fp.Mul(c1, b1)
+	c1 = e.fp.Sub(c1, t1)
+	tmp = mulFpByNonResidue(e.fp, t2)
+	c1 = e.fp.Add(c1, tmp)
+	tmp = e.fp.Add(&x.A0, &x.A2)
+	c2 := e.fp.Mul(b2, tmp)
+	c2 = e.fp.Sub(c2, t2)
+	c2 = e.fp.Add(c2, t1)
+	return &E3{
+		A0: *c0,
+		A1: *c1,
+		A2: *c2,
+	}
+}
+
+// Mul01By01 multiplies two E3 sparse element of the form:
+//
+//	E3{
+//		A0: c0,
+//		A1: c1,
+//		A2: 0,
+//	}
+//
+// and
+//
+//	E3{
+//		A0: d0,
+//		A1: d1,
+//		A2: 0,
+//	}
+func (e Ext3) Mul01By01(c0, c1, d0, d1 *baseEl) *E3 {
+	a := e.fp.Mul(d0, c0)
+	b := e.fp.Mul(d1, c1)
+	t0 := e.fp.Mul(c1, d1)
+	t0 = e.fp.Sub(t0, b)
+	t0 = mulFpByNonResidue(e.fp, t0)
+	t0 = e.fp.Add(t0, a)
+	t2 := e.fp.Mul(c0, d0)
+	t2 = e.fp.Sub(t2, a)
+	t2 = e.fp.Add(t2, b)
+	t1 := e.fp.Add(c0, c1)
+	tmp := e.fp.Add(d0, d1)
+	t1 = e.fp.Mul(t1, tmp)
+	t1 = e.fp.Sub(t1, a)
+	t1 = e.fp.Sub(t1, b)
+	return &E3{
+		A0: *t0,
+		A1: *t1,
+		A2: *t2,
+	}
+}
+
+func (e Ext3) Mul(x, y *E3) *E3 {
+	// Algorithm 13 from https://eprint.iacr.org/2010/354.pdf
+	t0 := e.fp.Mul(&x.A0, &y.A0)
+	t1 := e.fp.Mul(&x.A1, &y.A1)
+	t2 := e.fp.Mul(&x.A2, &y.A2)
+
+	c0 := e.fp.Add(&x.A1, &x.A2)
+	tmp := e.fp.Add(&y.A1, &y.A2)
+	c0 = e.fp.Mul(c0, tmp)
+	c0 = e.fp.Sub(c0, t1)
+	c0 = e.fp.Sub(c0, t2)
+	c0 = mulFpByNonResidue(e.fp, c0)
+
+	tmp = e.fp.Add(&x.A0, &x.A2)
+	c2 := e.fp.Add(&y.A0, &y.A2)
+	c2 = e.fp.Mul(c2, tmp)
+	c2 = e.fp.Sub(c2, t0)
+	c2 = e.fp.Sub(c2, t2)
+
+	c1 := e.fp.Add(&x.A0, &x.A1)
+	tmp = e.fp.Add(&y.A0, &y.A1)
+	c1 = e.fp.Mul(c1, tmp)
+	c1 = e.fp.Sub(c1, t0)
+	c1 = e.fp.Sub(c1, t1)
+	t2 = mulFpByNonResidue(e.fp, t2)
+
+	a0 := e.fp.Add(c0, t0)
+	a1 := e.fp.Add(c1, t2)
+	a2 := e.fp.Add(c2, t1)
+
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+func (e Ext3) Square(x *E3) *E3 {
+
+	// Algorithm 16 from https://eprint.iacr.org/2010/354.pdf
+
+	c6 := e.fp.MulConst(&x.A1, big.NewInt(2))
+	c4 := e.fp.Mul(&x.A0, c6) // x.A0 * xA1 * 2
+	c5 := e.fp.Mul(&x.A2, &x.A2)
+	c1 := mulFpByNonResidue(e.fp, c5)
+	c1 = e.fp.Add(c1, c4)
+	c2 := e.fp.Sub(c4, c5)
+
+	c3 := e.fp.Mul(&x.A0, &x.A0)
+	c4 = e.fp.Sub(&x.A0, &x.A1)
+	c4 = e.fp.Add(c4, &x.A2)
+	c5 = e.fp.Mul(c6, &x.A2) // x.A1 * xA2 * 2
+	c4 = e.fp.Mul(c4, c4)
+	c0 := mulFpByNonResidue(e.fp, c5)
+	c4 = e.fp.Add(c4, c5)
+	c4 = e.fp.Sub(c4, c3)
+
+	a0 := e.fp.Add(c0, c3)
+	a1 := c1
+	a2 := e.fp.Add(c2, c4)
+
+	return &E3{
+		A0: *a0,
+		A1: *a1,
+		A2: *a2,
+	}
+}
+
+func (e Ext3) Inverse(x *E3) *E3 {
+	res, err := e.fp.NewHint(inverseE3Hint, 3, &x.A0, &x.A1, &x.A2)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	inv := E3{
+		A0: *res[0],
+		A1: *res[1],
+		A2: *res[2],
+	}
+	one := e.One()
+
+	// 1 == inv * x
+	_one := e.Mul(&inv, x)
+	e.AssertIsEqual(one, _one)
+
+	return &inv
+
+}
+
+func (e Ext3) DivUnchecked(x, y *E3) *E3 {
+	res, err := e.fp.NewHint(divE3Hint, 6, &x.A0, &x.A1, &x.A2, &y.A0, &y.A1, &y.A2)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	div := E3{
+		A0: *res[0],
+		A1: *res[1],
+		A2: *res[2],
+	}
+
+	// x = div * y
+	_x := e.Mul(&div, y)
+	e.AssertIsEqual(x, _x)
+
+	return &div
+
+}
+
+// MulByNonResidue mul x by (0,1,0)
+func (e Ext3) MulByNonResidue(x *E3) *E3 {
+	z := &E3{
+		A0: x.A2,
+		A1: x.A0,
+		A2: x.A1,
+	}
+	z.A0 = *mulFpByNonResidue(e.fp, &z.A0)
+	return z
+}
+
+func (e Ext3) AssertIsEqual(a, b *E3) {
+	e.fp.AssertIsEqual(&a.A0, &b.A0)
+	e.fp.AssertIsEqual(&a.A1, &b.A1)
+	e.fp.AssertIsEqual(&a.A2, &b.A2)
+}
+
+func (e Ext3) Copy(x *E3) *E3 {
+	return &E3{
+		A0: x.A0,
+		A1: x.A1,
+		A2: x.A2,
+	}
+}
+
+func FromE3(a *bw6761.E3) E3 {
+	return E3{
+		A0: emulated.ValueOf[emulated.BW6761Fp](a.A0),
+		A1: emulated.ValueOf[emulated.BW6761Fp](a.A1),
+		A2: emulated.ValueOf[emulated.BW6761Fp](a.A2),
+	}
+}
diff --git a/std/algebra/emulated/fields_bw6761/e3_test.go b/std/algebra/emulated/fields_bw6761/e3_test.go
new file mode 100644
index 0000000000..707be87b67
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/e3_test.go
@@ -0,0 +1,363 @@
+package fields_bw6761
+
+import (
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark-crypto/ecc/bw6-761/fp"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/math/emulated"
+	"github.com/consensys/gnark/test"
+)
+
+type e3Add struct {
+	A, B, C E3
+}
+
+func (circuit *e3Add) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Add(&circuit.A, &circuit.B)
+	e.AssertIsEqual(expected, &circuit.C)
+	return nil
+}
+
+func TestAddFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E3
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Add(&a, &b)
+
+	witness := e3Add{
+		A: FromE3(&a),
+		B: FromE3(&b),
+		C: FromE3(&c),
+	}
+
+	err := test.IsSolved(&e3Add{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Sub struct {
+	A, B, C E3
+}
+
+func (circuit *e3Sub) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Sub(&circuit.A, &circuit.B)
+	e.AssertIsEqual(expected, &circuit.C)
+	return nil
+}
+
+func TestSubFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E3
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Sub(&a, &b)
+
+	witness := e3Sub{
+		A: FromE3(&a),
+		B: FromE3(&b),
+		C: FromE3(&c),
+	}
+
+	err := test.IsSolved(&e3Sub{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Neg struct {
+	A, B E3
+}
+
+func (circuit *e3Neg) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Neg(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestNegFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	b.Neg(&a)
+
+	witness := e3Neg{
+		A: FromE3(&a),
+		B: FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3Neg{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Double struct {
+	A, B E3
+}
+
+func (circuit *e3Double) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Double(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestDoubleFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	b.Double(&a)
+
+	witness := e3Double{
+		A: FromE3(&a),
+		B: FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3Double{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Mul struct {
+	A, B, C E3
+}
+
+func (circuit *e3Mul) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Mul(&circuit.A, &circuit.B)
+	e.AssertIsEqual(expected, &circuit.C)
+	return nil
+}
+
+func TestMulFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E3
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Mul(&a, &b)
+
+	witness := e3Mul{
+		A: FromE3(&a),
+		B: FromE3(&b),
+		C: FromE3(&c),
+	}
+
+	err := test.IsSolved(&e3Mul{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3MulByNonResidue struct {
+	A, B E3
+}
+
+func (circuit *e3MulByNonResidue) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.MulByNonResidue(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestMulByNonResidueFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	b.Set(&a)
+	b.MulByNonResidue(&a)
+
+	witness := e3MulByNonResidue{
+		A: FromE3(&a),
+		B: FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3MulByNonResidue{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3MulByElement struct {
+	A E3
+	Y baseEl
+	B E3
+}
+
+func (circuit *e3MulByElement) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.MulByElement(&circuit.A, &circuit.Y)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestMulByElementFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	var y fp.Element
+	y.SetRandom()
+	b.Set(&a)
+	b.MulByElement(&a, &y)
+
+	witness := e3MulByElement{
+		A: FromE3(&a),
+		Y: emulated.ValueOf[emulated.BW6761Fp](y),
+		B: FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3MulByElement{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3MulBy01 struct {
+	A      E3
+	C0, C1 baseEl
+	B      E3
+}
+
+func (circuit *e3MulBy01) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.MulBy01(&circuit.A, &circuit.C0, &circuit.C1)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestMulBy01Fp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	var c0, c1 fp.Element
+	c0.SetRandom()
+	c1.SetRandom()
+	b.Set(&a)
+	b.MulBy01(&c0, &c1)
+
+	witness := e3MulBy01{
+		A:  FromE3(&a),
+		C0: emulated.ValueOf[emulated.BW6761Fp](c0),
+		C1: emulated.ValueOf[emulated.BW6761Fp](c1),
+		B:  FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3MulBy01{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Square struct {
+	A, B E3
+}
+
+func (circuit *e3Square) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Square(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestSquareFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	b.Square(&a)
+
+	witness := e3Square{
+		A: FromE3(&a),
+		B: FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3Square{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Inverse struct {
+	A, B E3
+}
+
+func (circuit *e3Inverse) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Inverse(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestInverseFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	b.Inverse(&a)
+
+	witness := e3Inverse{
+		A: FromE3(&a),
+		B: FromE3(&b),
+	}
+
+	// add=50605 equals=769 fromBinary=0 mul=50315 sub=558 toBinary=0
+	err := test.IsSolved(&e3Inverse{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e3Div struct {
+	A, B, C E3
+}
+
+func (circuit *e3Div) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.DivUnchecked(&circuit.A, &circuit.B)
+	e.AssertIsEqual(expected, &circuit.C)
+	return nil
+}
+
+func TestDivFp3(t *testing.T) {
+
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E3
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Inverse(&b)
+	c.Mul(&a, &c)
+
+	witness := e3Div{
+		A: FromE3(&a),
+		B: FromE3(&b),
+		C: FromE3(&c),
+	}
+
+	err := test.IsSolved(&e3Div{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+}
+
+type e3Conjugate struct {
+	A, B E3
+}
+
+func (circuit *e3Conjugate) Define(api frontend.API) error {
+	e := NewExt3(api)
+	expected := e.Conjugate(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestConjugateFp3(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E3
+	_, _ = a.SetRandom()
+	b.Conjugate(&a)
+
+	witness := e3Conjugate{
+		A: FromE3(&a),
+		B: FromE3(&b),
+	}
+
+	err := test.IsSolved(&e3Conjugate{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
diff --git a/std/algebra/emulated/fields_bw6761/e6.go b/std/algebra/emulated/fields_bw6761/e6.go
new file mode 100644
index 0000000000..43340bdcd8
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/e6.go
@@ -0,0 +1,399 @@
+package fields_bw6761
+
+import (
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+type E6 struct {
+	B0, B1 E3
+}
+
+type Ext6 struct {
+	*Ext3
+}
+
+func (e Ext6) Reduce(x *E6) *E6 {
+	var z E6
+	z.B0 = *e.Ext3.Reduce(&x.B0)
+	z.B1 = *e.Ext3.Reduce(&x.B1)
+	return &z
+}
+
+func NewExt6(api frontend.API) *Ext6 {
+	return &Ext6{Ext3: NewExt3(api)}
+}
+
+func (e Ext6) Zero() *E6 {
+	b0 := e.Ext3.Zero()
+	b1 := e.Ext3.Zero()
+	return &E6{
+		B0: *b0,
+		B1: *b1,
+	}
+}
+
+func (e Ext6) One() *E6 {
+	return &E6{
+		B0: *e.Ext3.One(),
+		B1: *e.Ext3.Zero(),
+	}
+}
+
+func (e Ext6) Add(x, y *E6) *E6 {
+	return &E6{
+		B0: *e.Ext3.Add(&x.B0, &y.B0),
+		B1: *e.Ext3.Add(&x.B1, &y.B1),
+	}
+}
+
+func (e Ext6) Sub(x, y *E6) *E6 {
+	return &E6{
+		B0: *e.Ext3.Sub(&x.B0, &y.B0),
+		B1: *e.Ext3.Sub(&x.B1, &y.B1),
+	}
+}
+
+func (e Ext6) Double(x *E6) *E6 {
+	return &E6{
+		B0: *e.Ext3.Double(&x.B0),
+		B1: *e.Ext3.Double(&x.B1),
+	}
+}
+
+func (e Ext6) Mul(x, y *E6) *E6 {
+	x = e.Reduce(x)
+	y = e.Reduce(y)
+
+	a := e.Ext3.Add(&x.B0, &x.B1)
+	b := e.Ext3.Add(&y.B0, &y.B1)
+	a = e.Ext3.Mul(a, b)
+	b = e.Ext3.Mul(&x.B0, &y.B0)
+	c := e.Ext3.Mul(&x.B1, &y.B1)
+	b1 := e.Ext3.Sub(a, b)
+	b1 = e.Ext3.Sub(b1, c)
+	b0 := e.Ext3.MulByNonResidue(c)
+	b0 = e.Ext3.Add(b0, b)
+
+	return &E6{
+		B0: *b0,
+		B1: *b1,
+	}
+}
+
+func (e Ext6) Square(x *E6) *E6 {
+
+	x = e.Reduce(x)
+	//Algorithm 22 from https://eprint.iacr.org/2010/354.pdf
+	c0 := e.Ext3.Sub(&x.B0, &x.B1)
+	c3 := e.Ext3.MulByNonResidue(&x.B1)
+	c3 = e.Ext3.Neg(c3)
+	c3 = e.Ext3.Add(&x.B0, c3)
+	c2 := e.Ext3.Mul(&x.B0, &x.B1)
+	c0 = e.Ext3.Mul(c0, c3)
+	c0 = e.Ext3.Add(c0, c2)
+	b1 := e.Ext3.Double(c2)
+	c2 = e.Ext3.MulByNonResidue(c2)
+	b0 := e.Ext3.Add(c0, c2)
+
+	return &E6{
+		B0: *b0,
+		B1: *b1,
+	}
+}
+
+// Karabina's compressed cyclotomic square
+// https://eprint.iacr.org/2010/542.pdf
+// Th. 3.2 with minor modifications to fit our tower
+func (e Ext6) CyclotomicSquareCompressed(x *E6) *E6 {
+	x = e.Reduce(x)
+	z := e.Copy(x)
+
+	var t [7]*baseEl
+
+	// t0 = g1²
+	t[0] = e.fp.Mul(&x.B0.A1, &x.B0.A1)
+	// t1 = g5²
+	t[1] = e.fp.Mul(&x.B1.A2, &x.B1.A2)
+	// t5 = g1 + g5
+	t[5] = e.fp.Add(&x.B0.A1, &x.B1.A2)
+	// t2 = (g1 + g5)²
+	t[2] = e.fp.Mul(t[5], t[5])
+
+	// t3 = g1² + g5²
+	t[3] = e.fp.Add(t[0], t[1])
+	// t5 = 2 * g1 * g5
+	t[5] = e.fp.Sub(t[2], t[3])
+
+	// t6 = g3 + g2
+	t[6] = e.fp.Add(&x.B1.A0, &x.B0.A2)
+	// t3 = (g3 + g2)²
+	t[3] = e.fp.Mul(t[6], t[6])
+	// t2 = g3²
+	t[2] = e.fp.Mul(&x.B1.A0, &x.B1.A0)
+
+	// t6 = 2 * nr * g1 * g5
+	t[6] = mulFpByNonResidue(e.fp, t[5])
+	// t5 = 4 * nr * g1 * g5 + 2 * g3
+	t[5] = e.fp.Add(t[6], &x.B1.A0)
+	t[5] = e.fp.Add(t[5], t[5])
+	// z3 = 6 * nr * g1 * g5 + 2 * g3
+	z.B1.A0 = *e.fp.Add(t[5], t[6])
+
+	// t4 = nr * g5²
+	t[4] = mulFpByNonResidue(e.fp, t[1])
+	// t5 = nr * g5² + g1²
+	t[5] = e.fp.Add(t[0], t[4])
+	// t6 = nr * g5² + g1² - g2
+	t[6] = e.fp.Sub(t[5], &x.B0.A2)
+
+	// t1 = g2²
+	t[1] = e.fp.Mul(&x.B0.A2, &x.B0.A2)
+
+	// t6 = 2 * nr * g5² + 2 * g1² - 2*g2
+	t[6] = e.fp.Add(t[6], t[6])
+	// z2 = 3 * nr * g5² + 3 * g1² - 2*g2
+	z.B0.A2 = *e.fp.Add(t[6], t[5])
+
+	// t4 = nr * g2²
+	t[4] = mulFpByNonResidue(e.fp, t[1])
+	// t5 = g3² + nr * g2²
+	t[5] = e.fp.Add(t[2], t[4])
+	// t6 = g3² + nr * g2² - g1
+	t[6] = e.fp.Sub(t[5], &x.B0.A1)
+	// t6 = 2 * g3² + 2 * nr * g2² - 2 * g1
+	t[6] = e.fp.Add(t[6], t[6])
+	// z1 = 3 * g3² + 3 * nr * g2² - 2 * g1
+	z.B0.A1 = *e.fp.Add(t[6], t[5])
+
+	// t0 = g2² + g3²
+	t[0] = e.fp.Add(t[2], t[1])
+	// t5 = 2 * g3 * g2
+	t[5] = e.fp.Sub(t[3], t[0])
+	// t6 = 2 * g3 * g2 + g5
+	t[6] = e.fp.Add(t[5], &x.B1.A2)
+	// t6 = 4 * g3 * g2 + 2 * g5
+	t[6] = e.fp.Add(t[6], t[6])
+	// z5 = 6 * g3 * g2 + 2 * g5
+	z.B1.A2 = *e.fp.Add(t[5], t[6])
+
+	return z
+}
+
+// DecompressKarabina Karabina's cyclotomic square result
+// if g3 != 0
+//
+//	g4 = (E * g5^2 + 3 * g1^2 - 2 * g2)/4g3
+//
+// if g3 == 0
+//
+//	g4 = 2g1g5/g2
+//
+// if g3=g2=0 then g4=g5=g1=0 and g0=1 (x=1)
+// Theorem 3.1 is well-defined for all x in Gϕₙ\{1}
+//
+// TODO: handle the edge cases
+func (e Ext6) DecompressKarabina(x *E6) *E6 {
+
+	x = e.Reduce(x)
+
+	var z E6
+
+	var t [3]*baseEl
+	one := e.fp.One()
+
+	// t0 = g1^2
+	t[0] = e.fp.Mul(&x.B0.A1, &x.B0.A1)
+	// t1 = 3 * g1^2 - 2 * g2
+	t[1] = e.fp.Sub(t[0], &x.B0.A2)
+	t[1] = e.fp.Add(t[1], t[1])
+	t[1] = e.fp.Add(t[1], t[0])
+	// t0 = E * g5^2 + t1
+	t[2] = e.fp.Mul(&x.B1.A2, &x.B1.A2)
+	t[0] = mulFpByNonResidue(e.fp, t[2])
+	t[0] = e.fp.Add(t[0], t[1])
+	// t1 = 1/(4 * g3)
+	t[1] = e.fp.Add(&x.B1.A0, &x.B1.A0)
+	t[1] = e.fp.Add(t[1], t[1])
+
+	// z4 = g4
+	z.B1.A1 = *e.fp.Div(t[0], t[1])
+	a1 := z.B1.A1
+
+	// t1 = g2 * g1
+	t[1] = e.fp.Mul(&x.B0.A2, &x.B0.A1)
+	// t2 = 2 * g4² - 3 * g2 * g1
+	t[2] = e.fp.Mul(&a1, &a1)
+	t[2] = e.fp.Sub(t[2], t[1])
+	t[2] = e.fp.Add(t[2], t[2])
+	t[2] = e.fp.Sub(t[2], t[1])
+	// t1 = g3 * g5 (g3 can be 0)
+	t[1] = e.fp.Mul(&x.B1.A0, &x.B1.A2)
+	// c₀ = E * (2 * g4² + g3 * g5 - 3 * g2 * g1) + 1
+	t[2] = e.fp.Add(t[2], t[1])
+
+	z.B0.A0 = *mulFpByNonResidue(e.fp, t[2])
+	z.B0.A0 = *e.fp.Add(&z.B0.A0, one)
+
+	z.B0.A1 = x.B0.A1
+	z.B0.A2 = x.B0.A2
+	z.B1.A0 = x.B1.A0
+	z.B1.A2 = x.B1.A2
+
+	return &z
+}
+
+// Granger-Scott's cyclotomic square
+// https://eprint.iacr.org/2009/565.pdf, 3.2
+func (e Ext6) CyclotomicSquare(x *E6) *E6 {
+	// x=(x0,x1,x2,x3,x4,x5,x6,x7) in E3⁶
+	// cyclosquare(x)=(3*x4²*u + 3*x0² - 2*x0,
+	//					3*x2²*u + 3*x3² - 2*x1,
+	//					3*x5²*u + 3*x1² - 2*x2,
+	//					6*x1*x5*u + 2*x3,
+	//					6*x0*x4 + 2*x4,
+	//					6*x2*x3 + 2*x5)
+
+	x = e.Reduce(x)
+
+	var t [9]*baseEl
+
+	t[0] = e.fp.Mul(&x.B1.A1, &x.B1.A1)
+	t[1] = e.fp.Mul(&x.B0.A0, &x.B0.A0)
+	t[6] = e.fp.Add(&x.B1.A1, &x.B0.A0)
+	t[6] = e.fp.Mul(t[6], t[6])
+	t[6] = e.fp.Sub(t[6], t[0])
+	t[6] = e.fp.Sub(t[6], t[1]) // 2*x4*x0
+	t[2] = e.fp.Mul(&x.B0.A2, &x.B0.A2)
+	t[3] = e.fp.Mul(&x.B1.A0, &x.B1.A0)
+	t[7] = e.fp.Add(&x.B0.A2, &x.B1.A0)
+	t[7] = e.fp.Mul(t[7], t[7])
+	t[7] = e.fp.Sub(t[7], t[2])
+	t[7] = e.fp.Sub(t[7], t[3]) // 2*x2*x3
+	t[4] = e.fp.Mul(&x.B1.A2, &x.B1.A2)
+	t[5] = e.fp.Mul(&x.B0.A1, &x.B0.A1)
+	t[8] = e.fp.Add(&x.B1.A2, &x.B0.A1)
+	t[8] = e.fp.Mul(t[8], t[8])
+	t[8] = e.fp.Sub(t[8], t[4])
+	t[8] = e.fp.Sub(t[8], t[5])
+	t[8] = mulFpByNonResidue(e.fp, t[8]) // 2*x5*x1*u
+
+	t[0] = mulFpByNonResidue(e.fp, t[0])
+	t[0] = e.fp.Add(t[0], t[1]) // x4²*u + x0²
+	t[2] = mulFpByNonResidue(e.fp, t[2])
+	t[2] = e.fp.Add(t[2], t[3]) // x2²*u + x3²
+	t[4] = mulFpByNonResidue(e.fp, t[4])
+	t[4] = e.fp.Add(t[4], t[5]) // x5²*u + x1²
+
+	var z E6
+	z.B0.A0 = *e.fp.Sub(t[0], &x.B0.A0)
+	z.B0.A0 = *e.fp.Add(&z.B0.A0, &z.B0.A0)
+	z.B0.A0 = *e.fp.Add(&z.B0.A0, t[0])
+	z.B0.A1 = *e.fp.Sub(t[2], &x.B0.A1)
+	z.B0.A1 = *e.fp.Add(&z.B0.A1, &z.B0.A1)
+	z.B0.A1 = *e.fp.Add(&z.B0.A1, t[2])
+	z.B0.A2 = *e.fp.Sub(t[4], &x.B0.A2)
+	z.B0.A2 = *e.fp.Add(&z.B0.A2, &z.B0.A2)
+	z.B0.A2 = *e.fp.Add(&z.B0.A2, t[4])
+
+	z.B1.A0 = *e.fp.Add(t[8], &x.B1.A0)
+	z.B1.A0 = *e.fp.Add(&z.B1.A0, &z.B1.A0)
+	z.B1.A0 = *e.fp.Add(&z.B1.A0, t[8])
+	z.B1.A1 = *e.fp.Add(t[6], &x.B1.A1)
+	z.B1.A1 = *e.fp.Add(&z.B1.A1, &z.B1.A1)
+	z.B1.A1 = *e.fp.Add(&z.B1.A1, t[6])
+	z.B1.A2 = *e.fp.Add(t[7], &x.B1.A2)
+	z.B1.A2 = *e.fp.Add(&z.B1.A2, &z.B1.A2)
+	z.B1.A2 = *e.fp.Add(&z.B1.A2, t[7])
+
+	return &z
+}
+
+func (e Ext6) Inverse(x *E6) *E6 {
+	res, err := e.fp.NewHint(inverseE6Hint, 6, &x.B0.A0, &x.B0.A1, &x.B0.A2, &x.B1.A0, &x.B1.A1, &x.B1.A2)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	inv := E6{
+		B0: E3{A0: *res[0], A1: *res[1], A2: *res[2]},
+		B1: E3{A0: *res[3], A1: *res[4], A2: *res[5]},
+	}
+	one := e.One()
+
+	// 1 == inv * x
+	_one := e.Mul(&inv, x)
+	e.AssertIsEqual(one, _one)
+
+	return &inv
+
+}
+
+func (e Ext6) DivUnchecked(x, y *E6) *E6 {
+	res, err := e.fp.NewHint(divE6Hint, 12, &x.B0.A0, &x.B0.A1, &x.B0.A2, &x.B1.A0, &x.B1.A1, &x.B1.A2, &y.B0.A0, &y.B0.A1, &y.B0.A2, &y.B1.A0, &y.B1.A1, &y.B1.A2)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	div := E6{
+		B0: E3{A0: *res[0], A1: *res[1], A2: *res[2]},
+		B1: E3{A0: *res[3], A1: *res[4], A2: *res[5]},
+	}
+
+	// x = div * y
+	_x := e.Mul(&div, y)
+	e.AssertIsEqual(x, _x)
+
+	return &div
+
+}
+
+func (e Ext6) Conjugate(x *E6) *E6 {
+	return &E6{
+		B0: x.B0,
+		B1: *e.Ext3.Neg(&x.B1),
+	}
+}
+
+func (e Ext6) AssertIsEqual(a, b *E6) {
+	e.Ext3.AssertIsEqual(&a.B0, &b.B0)
+	e.Ext3.AssertIsEqual(&a.B1, &b.B1)
+}
+
+func (e Ext6) Copy(x *E6) *E6 {
+	b0 := e.Ext3.Copy(&x.B0)
+	b1 := e.Ext3.Copy(&x.B1)
+	return &E6{
+		B0: *b0,
+		B1: *b1,
+	}
+}
+
+func FromE6(a *bw6761.E6) E6 {
+	return E6{
+		B0: FromE3(&a.B0),
+		B1: FromE3(&a.B1),
+	}
+}
+
+// Frobenius set z in E6 to Frobenius(x), return z
+func (e Ext6) Frobenius(x *E6) *E6 {
+	_frobA := emulated.ValueOf[emulated.BW6761Fp]("4922464560225523242118178942575080391082002530232324381063048548642823052024664478336818169867474395270858391911405337707247735739826664939444490469542109391530482826728203582549674992333383150446779312029624171857054392282775648")
+	_frobB := emulated.ValueOf[emulated.BW6761Fp]("1968985824090209297278610739700577151397666382303825728450741611566800370218827257750865013421937292370006175842381275743914023380727582819905021229583192207421122272650305267822868639090213645505120388400344940985710520836292650")
+	_frobC := emulated.ValueOf[emulated.BW6761Fp]("4922464560225523242118178942575080391082002530232324381063048548642823052024664478336818169867474395270858391911405337707247735739826664939444490469542109391530482826728203582549674992333383150446779312029624171857054392282775649")
+	_frobAC := emulated.ValueOf[emulated.BW6761Fp]("-1")
+	_frobBC := emulated.ValueOf[emulated.BW6761Fp]("1968985824090209297278610739700577151397666382303825728450741611566800370218827257750865013421937292370006175842381275743914023380727582819905021229583192207421122272650305267822868639090213645505120388400344940985710520836292651")
+	var z E6
+	z.B0.A0 = x.B0.A0
+	z.B0.A1 = *e.fp.Mul(&x.B0.A1, &_frobA)
+	z.B0.A2 = *e.fp.Mul(&x.B0.A2, &_frobB)
+
+	z.B1.A0 = *e.fp.Mul(&x.B1.A0, &_frobC)
+	z.B1.A1 = *e.fp.Mul(&x.B1.A1, &_frobAC)
+	z.B1.A2 = *e.fp.Mul(&x.B1.A2, &_frobBC)
+
+	return &z
+}
diff --git a/std/algebra/emulated/fields_bw6761/e6_pairing.go b/std/algebra/emulated/fields_bw6761/e6_pairing.go
new file mode 100644
index 0000000000..1e55d911a1
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/e6_pairing.go
@@ -0,0 +1,282 @@
+package fields_bw6761
+
+func (e Ext6) nSquareCompressed(z *E6, n int) *E6 {
+	for i := 0; i < n; i++ {
+		z = e.CyclotomicSquareCompressed(z)
+	}
+	return z
+}
+
+// ExpX0Minus1 set z to z^{x₀-1} in E6 and return z
+// x₀-1 = 91893752504881257682351033800651177983
+func (e Ext6) ExpX0Minus1(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.Copy(z)
+	result = e.nSquareCompressed(result, 5)
+	result = e.DecompressKarabina(result)
+	result = e.Mul(result, z)
+	z33 := e.Copy(result)
+	result = e.nSquareCompressed(result, 7)
+	result = e.DecompressKarabina(result)
+	result = e.Mul(result, z33)
+	result = e.nSquareCompressed(result, 4)
+	result = e.DecompressKarabina(result)
+	result = e.Mul(result, z)
+	result = e.CyclotomicSquare(result)
+	result = e.Mul(result, z)
+	result = e.nSquareCompressed(result, 46)
+	result = e.DecompressKarabina(result)
+
+	return result
+}
+
+// ExpX0Minus1Square set z to z^{(x₀-1)²} in E6 and return z
+// (x₀-1)² = 91893752504881257682351033800651177984
+func (e Ext6) ExpX0Minus1Square(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.Copy(z)
+	result = e.CyclotomicSquare(result)
+	t0 := e.Mul(z, result)
+	t1 := e.CyclotomicSquare(t0)
+	t0 = e.Mul(t0, t1)
+	result = e.Mul(result, t0)
+	t1 = e.Mul(t1, result)
+	t0 = e.Mul(t0, t1)
+	t2 := e.CyclotomicSquare(t0)
+	t2 = e.Mul(t1, t2)
+	t0 = e.Mul(t0, t2)
+	t2 = e.nSquareCompressed(t2, 7)
+	t2 = e.DecompressKarabina(t2)
+	t1 = e.Mul(t1, t2)
+	t1 = e.nSquareCompressed(t1, 11)
+	t1 = e.DecompressKarabina(t1)
+	t1 = e.Mul(t0, t1)
+	t1 = e.nSquareCompressed(t1, 9)
+	t1 = e.DecompressKarabina(t1)
+	t0 = e.Mul(t0, t1)
+	t0 = e.CyclotomicSquare(t0)
+	result = e.Mul(result, t0)
+	result = e.nSquareCompressed(result, 92)
+	result = e.DecompressKarabina(result)
+
+	return result
+
+}
+
+// ExpX0Plus1 set z to z^(x₀+1) in E6 and return z
+// x₀+1 = 91893752504881257682351033800651177985
+func (e Ext6) ExpX0Plus1(z *E6) *E6 {
+	result := e.ExpX0Minus1(z)
+	t := e.CyclotomicSquare(z)
+	result = e.Mul(result, t)
+	return result
+}
+
+// ExpX0Minus1Div3 set z to z^(x₀-1)/3 in E6 and return z
+// (x₀-1)/3 = 3195374304363544576
+func (e Ext6) ExptMinus1Div3(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.Copy(z)
+	result = e.CyclotomicSquare(result)
+	result = e.Mul(result, z)
+	t0 := e.Mul(result, z)
+	t0 = e.CyclotomicSquare(t0)
+	result = e.Mul(result, t0)
+	t0 = result
+	t0 = e.nSquareCompressed(t0, 7)
+	t0 = e.DecompressKarabina(t0)
+	result = e.Mul(result, t0)
+	result = e.nSquareCompressed(result, 5)
+	result = e.DecompressKarabina(result)
+	result = e.Mul(result, z)
+	result = e.nSquareCompressed(result, 46)
+	result = e.DecompressKarabina(result)
+
+	return result
+}
+
+// ExpC1 set z to z^C1 in E6 and return z
+// ht, hy = 13, 9
+// C1 = (ht+hy)/2 = 11
+func (e Ext6) ExpC1(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.CyclotomicSquare(z)
+	result = e.Mul(result, z)
+	t0 := e.Mul(z, result)
+	t0 = e.CyclotomicSquare(t0)
+	result = e.Mul(result, t0)
+
+	return result
+}
+
+// ExpC2 set z to z^C2 in E6 and return z
+// ht, hy = 13, 9
+// C2 = (ht**2+3*hy**2)/4 = 103
+func (e Ext6) ExpC2(z *E6) *E6 {
+	z = e.Reduce(z)
+
+	result := e.CyclotomicSquare(z)
+	result = e.Mul(result, z)
+	t0 := result
+	t0 = e.nSquareCompressed(t0, 4)
+	t0 = e.DecompressKarabina(t0)
+	result = e.Mul(result, t0)
+	result = e.CyclotomicSquare(result)
+	result = e.Mul(result, z)
+
+	return result
+}
+
+// Square034 squares an E6 sparse element of the form
+//
+//	E6{
+//		B0: E3{A0: 1,  A1: 0,  A2: 0},
+//		B1: E3{A0: c3, A1: c4, A2: 0},
+//	}
+func (e *Ext6) Square034(x *E6) *E6 {
+	x = e.Reduce(x)
+	c0 := &E3{
+		A0: *e.fp.Sub(&x.B0.A0, &x.B1.A0),
+		A1: *e.fp.Neg(&x.B1.A1),
+		A2: *e.fp.Zero(),
+	}
+
+	c3 := &E3{
+		A0: x.B0.A0,
+		A1: *e.fp.Neg(&x.B1.A0),
+		A2: *e.fp.Neg(&x.B1.A1),
+	}
+
+	c2 := &E3{
+		A0: x.B1.A0,
+		A1: x.B1.A1,
+		A2: *e.fp.Zero(),
+	}
+	c3 = e.MulBy01(c3, &c0.A0, &c0.A1)
+	c3 = e.Ext3.Add(c3, c2)
+
+	var z E6
+	z.B1.A0 = *e.fp.Add(&c2.A0, &c2.A0)
+	z.B1.A1 = *e.fp.Add(&c2.A1, &c2.A1)
+
+	z.B0.A0 = c3.A0
+	z.B0.A1 = *e.fp.Add(&c3.A1, &c2.A0)
+	z.B0.A2 = *e.fp.Add(&c3.A2, &c2.A1)
+
+	return &z
+}
+
+// MulBy034 multiplies z by an E6 sparse element of the form
+//
+//	E6{
+//		B0: E3{A0: 1,  A1: 0,  A2: 0},
+//		B1: E3{A0: c3, A1: c4, A2: 0},
+//	}
+func (e *Ext6) MulBy034(z *E6, c3, c4 *baseEl) *E6 {
+
+	z = e.Reduce(z)
+	a := z.B0
+	b := e.MulBy01(&z.B1, c3, c4)
+	c3 = e.fp.Add(e.fp.One(), c3)
+	d := e.Ext3.Add(&z.B0, &z.B1)
+	d = e.MulBy01(d, c3, c4)
+
+	zC1 := e.Ext3.Add(&a, b)
+	zC1 = e.Ext3.Neg(zC1)
+	zC1 = e.Ext3.Add(zC1, d)
+	zC0 := e.Ext3.MulByNonResidue(b)
+	zC0 = e.Ext3.Add(zC0, &a)
+
+	return &E6{
+		B0: *zC0,
+		B1: *zC1,
+	}
+}
+
+//	multiplies two E6 sparse element of the form:
+//
+//	E6{
+//		C0: E6{B0: 1, B1: 0, B2: 0},
+//		C1: E6{B0: c3, B1: c4, B2: 0},
+//	}
+//
+// and
+//
+//	E6{
+//		C0: E6{B0: 1, B1: 0, B2: 0},
+//		C1: E6{B0: d3, B1: d4, B2: 0},
+//	}
+func (e *Ext6) Mul034By034(d3, d4, c3, c4 *baseEl) [5]*baseEl {
+	x3 := e.fp.Mul(c3, d3)
+	x4 := e.fp.Mul(c4, d4)
+	x04 := e.fp.Add(c4, d4)
+	x03 := e.fp.Add(c3, d3)
+	tmp := e.fp.Add(c3, c4)
+	x34 := e.fp.Add(d3, d4)
+	x34 = e.fp.Mul(x34, tmp)
+	x34 = e.fp.Sub(x34, x3)
+	x34 = e.fp.Sub(x34, x4)
+
+	zC0B0 := mulFpByNonResidue(e.fp, x4)
+	zC0B0 = e.fp.Add(zC0B0, e.fp.One())
+	zC0B1 := x3
+	zC0B2 := x34
+	zC1B0 := x03
+	zC1B1 := x04
+
+	return [5]*baseEl{zC0B0, zC0B1, zC0B2, zC1B0, zC1B1}
+}
+
+// MulBy01234 multiplies z by an E6 sparse element of the form
+//
+//	E6{
+//		C0: E3{A0: c0, A1: c1, A2: c2},
+//		C1: E3{A0: c3, A1: c4, A2: 0},
+//	}
+func (e *Ext6) MulBy01234(z *E6, x [5]*baseEl) *E6 {
+	c0 := &E3{A0: *x[0], A1: *x[1], A2: *x[2]}
+	c1 := &E3{A0: *x[3], A1: *x[4], A2: *e.fp.Zero()}
+	a := e.Ext3.Add(&z.B0, &z.B1)
+	b := e.Ext3.Add(c0, c1)
+	a = e.Ext3.Mul(a, b)
+	b = e.Ext3.Mul(&z.B0, c0)
+	c := e.Ext3.MulBy01(&z.B1, x[3], x[4])
+	z1 := e.Ext3.Sub(a, b)
+	z1 = e.Ext3.Sub(z1, c)
+	z0 := e.Ext3.MulByNonResidue(c)
+	z0 = e.Ext3.Add(z0, b)
+	return &E6{
+		B0: *z0,
+		B1: *z1,
+	}
+}
+
+//	multiplies two E6 sparse element of the form:
+//
+//	E6{
+//		C0: E2{A0: x0, A1: x1, A2: x2},
+//		C1: E2{A0: x3, A1: x4, A2: 0},
+//	}
+//
+// and
+//
+//	E6{
+//		C0: E3{A0: 1,  A1: 0,  A2: 0},
+//		C1: E3{A0: z3, A1: z4, A2: 0},
+//	}
+func (e *Ext6) Mul01234By034(x [5]*baseEl, z3, z4 *baseEl) *E6 {
+	c0 := &E3{A0: *x[0], A1: *x[1], A2: *x[2]}
+	c1 := &E3{A0: *x[3], A1: *x[4], A2: *e.fp.Zero()}
+	a := e.Ext3.Add(e.Ext3.One(), &E3{A0: *z3, A1: *z4, A2: *e.fp.Zero()})
+	b := e.Ext3.Add(c0, c1)
+	a = e.Ext3.Mul(a, b)
+	c := e.Ext3.Mul01By01(z3, z4, x[3], x[4])
+	z1 := e.Ext3.Sub(a, c0)
+	z1 = e.Ext3.Sub(z1, c)
+	z0 := e.Ext3.MulByNonResidue(c)
+	z0 = e.Ext3.Add(z0, c0)
+	return &E6{
+		B0: *z0,
+		B1: *z1,
+	}
+}
diff --git a/std/algebra/emulated/fields_bw6761/e6_test.go b/std/algebra/emulated/fields_bw6761/e6_test.go
new file mode 100644
index 0000000000..77033ba24c
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/e6_test.go
@@ -0,0 +1,404 @@
+package fields_bw6761
+
+import (
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark-crypto/ecc/bw6-761/fp"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/math/emulated"
+	"github.com/consensys/gnark/test"
+)
+
+type e6Add struct {
+	A, B, C E6
+}
+
+func (circuit *e6Add) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Add(&circuit.A, &circuit.B)
+	e.AssertIsEqual(&expected, &circuit.C)
+	return nil
+}
+
+func TestAddFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E6
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Add(&a, &b)
+
+	witness := e6Add{
+		A: FromE6(&a),
+		B: FromE6(&b),
+		C: FromE6(&c),
+	}
+
+	err := test.IsSolved(&e6Add{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Sub struct {
+	A, B, C E6
+}
+
+func (circuit *e6Sub) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Sub(&circuit.A, &circuit.B)
+	e.AssertIsEqual(&expected, &circuit.C)
+	return nil
+}
+
+func TestSubFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E6
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Sub(&a, &b)
+
+	witness := e6Sub{
+		A: FromE6(&a),
+		B: FromE6(&b),
+		C: FromE6(&c),
+	}
+
+	err := test.IsSolved(&e6Sub{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Double struct {
+	A, B E6
+}
+
+func (circuit *e6Double) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Double(&circuit.A)
+	e.AssertIsEqual(&expected, &circuit.B)
+	return nil
+}
+
+func TestDoubleFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Double(&a)
+
+	witness := e6Double{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6Double{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Mul struct {
+	A, B, C E6
+}
+
+func (circuit *e6Mul) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Mul(&circuit.A, &circuit.B)
+	e.AssertIsEqual(&expected, &circuit.C)
+	return nil
+}
+
+func TestMulFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E6
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Mul(&a, &b)
+
+	witness := e6Mul{
+		A: FromE6(&a),
+		B: FromE6(&b),
+		C: FromE6(&c),
+	}
+
+	err := test.IsSolved(&e6Mul{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Square struct {
+	A, B E6
+}
+
+func (circuit *e6Square) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Square(&circuit.A)
+	e.AssertIsEqual(&expected, &circuit.B)
+	return nil
+}
+
+func TestSquareFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Square(&a)
+
+	witness := e6Square{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6Square{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Inverse struct {
+	A, B E6
+}
+
+func (circuit *e6Inverse) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Inverse(&circuit.A)
+	e.AssertIsEqual(&expected, &circuit.B)
+	return nil
+}
+
+func TestInverseFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Inverse(&a)
+
+	witness := e6Inverse{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6Inverse{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Div struct {
+	A, B, C E6
+}
+
+func (circuit *e6Div) Define(api frontend.API) error {
+	e := NewExt6(api)
+	expected := e.DivUnchecked(&circuit.A, &circuit.B)
+	e.AssertIsEqual(expected, &circuit.C)
+	return nil
+}
+
+func TestDivFp6(t *testing.T) {
+
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b, c bw6761.E6
+	_, _ = a.SetRandom()
+	_, _ = b.SetRandom()
+	c.Inverse(&b)
+	c.Mul(&a, &c)
+
+	witness := e6Div{
+		A: FromE6(&a),
+		B: FromE6(&b),
+		C: FromE6(&c),
+	}
+
+	err := test.IsSolved(&e6Div{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+}
+
+type e6Conjugate struct {
+	A, B E6
+}
+
+func (circuit *e6Conjugate) Define(api frontend.API) error {
+	var expected E6
+	e := NewExt6(api)
+	expected = *e.Conjugate(&circuit.A)
+	e.AssertIsEqual(&expected, &circuit.B)
+	return nil
+}
+
+func TestConjugateFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Conjugate(&a)
+
+	witness := e6Conjugate{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6Conjugate{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6CyclotomicSquareCompressed struct {
+	A, B E6
+}
+
+func (circuit *e6CyclotomicSquareCompressed) Define(api frontend.API) error {
+	e := NewExt6(api)
+	expected := e.CyclotomicSquareCompressed(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestCyclotomicSquareCompressedFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Set(&a)
+	b.CyclotomicSquareCompressed(&a)
+
+	witness := e6CyclotomicSquareCompressed{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6CyclotomicSquareCompressed{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6DecompressKarabina struct {
+	A, B E6
+}
+
+func (circuit *e6DecompressKarabina) Define(api frontend.API) error {
+	e := NewExt6(api)
+	expected := e.DecompressKarabina(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestDecompressKarabinaFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Set(&a)
+	a.DecompressKarabina(&a)
+
+	witness := e6DecompressKarabina{
+		A: FromE6(&b),
+		B: FromE6(&a),
+	}
+
+	err := test.IsSolved(&e6DecompressKarabina{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6CyclotomicSquare struct {
+	A, B E6
+}
+
+func (circuit *e6CyclotomicSquare) Define(api frontend.API) error {
+	e := NewExt6(api)
+	expected := e.CyclotomicSquare(&circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestCyclotomicSquareFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+	b.Set(&a)
+	b.CyclotomicSquare(&a)
+
+	witness := e6CyclotomicSquare{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6CyclotomicSquare{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6Expt struct {
+	A, B E6
+}
+
+func (circuit *e6Expt) Define(api frontend.API) error {
+	e := NewExt6(api)
+	expected := e.ExpX0Minus1(&circuit.A)
+	expected = e.Mul(expected, &circuit.A)
+	e.AssertIsEqual(expected, &circuit.B)
+	return nil
+}
+
+func TestExptFp6(t *testing.T) {
+	assert := test.NewAssert(t)
+	// witness values
+	var a, b bw6761.E6
+	_, _ = a.SetRandom()
+
+	// put a in the cyclotomic subgroup
+	var tmp bw6761.E6
+	tmp.Conjugate(&a)
+	a.Inverse(&a)
+	tmp.Mul(&tmp, &a)
+	a.Frobenius(&tmp).Mul(&a, &tmp)
+
+	b.Expt(&a)
+
+	witness := e6Expt{
+		A: FromE6(&a),
+		B: FromE6(&b),
+	}
+
+	err := test.IsSolved(&e6Expt{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type e6MulBy034 struct {
+	A    E6 `gnark:",public"`
+	W    E6
+	B, C baseEl
+}
+
+func (circuit *e6MulBy034) Define(api frontend.API) error {
+	e := NewExt6(api)
+	res := e.MulBy034(&circuit.A, &circuit.B, &circuit.C)
+	e.AssertIsEqual(res, &circuit.W)
+	return nil
+}
+
+func TestFp12MulBy034(t *testing.T) {
+
+	assert := test.NewAssert(t)
+	// witness values
+	var a, w bw6761.E6
+	_, _ = a.SetRandom()
+	var one, b, c fp.Element
+	one.SetOne()
+	_, _ = b.SetRandom()
+	_, _ = c.SetRandom()
+	w.Set(&a)
+	w.MulBy034(&one, &b, &c)
+
+	witness := e6MulBy034{
+		A: FromE6(&a),
+		B: emulated.ValueOf[emulated.BW6761Fp](&b),
+		C: emulated.ValueOf[emulated.BW6761Fp](&c),
+		W: FromE6(&w),
+	}
+
+	err := test.IsSolved(&e6MulBy034{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+}
diff --git a/std/algebra/emulated/fields_bw6761/hints.go b/std/algebra/emulated/fields_bw6761/hints.go
new file mode 100644
index 0000000000..5c100d085c
--- /dev/null
+++ b/std/algebra/emulated/fields_bw6761/hints.go
@@ -0,0 +1,124 @@
+package fields_bw6761
+
+import (
+	"math/big"
+
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"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{
+		// E3
+		divE3Hint,
+		inverseE3Hint,
+		// E6
+		divE6Hint,
+		inverseE6Hint,
+	}
+}
+
+// E3
+func inverseE3Hint(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 bw6761.E3
+
+			a.A0.SetBigInt(inputs[0])
+			a.A1.SetBigInt(inputs[1])
+			a.A2.SetBigInt(inputs[2])
+
+			c.Inverse(&a)
+
+			c.A0.BigInt(outputs[0])
+			c.A1.BigInt(outputs[1])
+			c.A2.BigInt(outputs[2])
+
+			return nil
+		})
+}
+
+func divE3Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var a, b, c bw6761.E3
+
+			a.A0.SetBigInt(inputs[0])
+			a.A1.SetBigInt(inputs[1])
+			a.A2.SetBigInt(inputs[2])
+			b.A0.SetBigInt(inputs[3])
+			b.A1.SetBigInt(inputs[4])
+			b.A2.SetBigInt(inputs[5])
+
+			c.Inverse(&b).Mul(&c, &a)
+
+			c.A0.BigInt(outputs[0])
+			c.A1.BigInt(outputs[1])
+			c.A2.BigInt(outputs[2])
+
+			return nil
+		})
+}
+
+// E6
+func inverseE6Hint(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 bw6761.E6
+
+			a.B0.A0.SetBigInt(inputs[0])
+			a.B0.A1.SetBigInt(inputs[1])
+			a.B0.A2.SetBigInt(inputs[2])
+			a.B1.A0.SetBigInt(inputs[3])
+			a.B1.A1.SetBigInt(inputs[4])
+			a.B1.A2.SetBigInt(inputs[5])
+
+			c.Inverse(&a)
+
+			c.B0.A0.BigInt(outputs[0])
+			c.B0.A1.BigInt(outputs[1])
+			c.B0.A2.BigInt(outputs[2])
+			c.B1.A0.BigInt(outputs[3])
+			c.B1.A1.BigInt(outputs[4])
+			c.B1.A2.BigInt(outputs[5])
+
+			return nil
+		})
+}
+
+func divE6Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var a, b, c bw6761.E6
+
+			a.B0.A0.SetBigInt(inputs[0])
+			a.B0.A1.SetBigInt(inputs[1])
+			a.B0.A2.SetBigInt(inputs[2])
+			a.B1.A0.SetBigInt(inputs[3])
+			a.B1.A1.SetBigInt(inputs[4])
+			a.B1.A2.SetBigInt(inputs[5])
+			b.B0.A0.SetBigInt(inputs[6])
+			b.B0.A1.SetBigInt(inputs[7])
+			b.B0.A2.SetBigInt(inputs[8])
+			b.B1.A0.SetBigInt(inputs[9])
+			b.B1.A1.SetBigInt(inputs[10])
+			b.B1.A2.SetBigInt(inputs[11])
+
+			c.Inverse(&b).Mul(&c, &a)
+
+			c.B0.A0.BigInt(outputs[0])
+			c.B0.A1.BigInt(outputs[1])
+			c.B0.A2.BigInt(outputs[2])
+			c.B1.A0.BigInt(outputs[3])
+			c.B1.A1.BigInt(outputs[4])
+			c.B1.A2.BigInt(outputs[5])
+
+			return nil
+		})
+}
diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go
index 9f4996d457..7969a7dcbe 100644
--- a/std/algebra/emulated/sw_bls12381/pairing.go
+++ b/std/algebra/emulated/sw_bls12381/pairing.go
@@ -346,7 +346,7 @@ func (pr Pairing) MillerLoop(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
 	// i = 62, separately to avoid an E12 Square
 	// (Square(res) = 1² = 1)
 
-	// k = 0, separately to avoid MulBy034 (res × ℓ)
+	// k = 0, separately to avoid MulBy014 (res × ℓ)
 
 	// Qacc[k] ← 3Qacc[k],
 	// l1 the tangent ℓ to 2Q[k]
diff --git a/std/algebra/emulated/sw_bw6761/doc.go b/std/algebra/emulated/sw_bw6761/doc.go
new file mode 100644
index 0000000000..e608cc7d1c
--- /dev/null
+++ b/std/algebra/emulated/sw_bw6761/doc.go
@@ -0,0 +1,7 @@
+// Package sw_bw6761 implements G1 and G2 arithmetics and pairing computation over BW6-761 curve.
+//
+// The implementation follows [Housni22]: "Pairings in Rank-1 Constraint Systems" and [BW6-761-hackmd].
+//
+// [Housni22]: https://eprint.iacr.org/2022/1162
+// [BW6-761-hackmd]: https://hackmd.io/@gnark/BW6-761-changes
+package sw_bw6761
diff --git a/std/algebra/emulated/sw_bw6761/doc_test.go b/std/algebra/emulated/sw_bw6761/doc_test.go
new file mode 100644
index 0000000000..e524186f3a
--- /dev/null
+++ b/std/algebra/emulated/sw_bw6761/doc_test.go
@@ -0,0 +1,103 @@
+package sw_bw6761_test
+
+import (
+	"crypto/rand"
+	"fmt"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/backend/groth16"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/frontend/cs/r1cs"
+	"github.com/consensys/gnark/std/algebra/emulated/sw_bw6761"
+)
+
+type PairCircuit struct {
+	InG1 sw_bw6761.G1Affine
+	InG2 sw_bw6761.G2Affine
+	Res  sw_bw6761.GTEl
+}
+
+func (c *PairCircuit) Define(api frontend.API) error {
+	pairing, err := sw_bw6761.NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	// Pair method does not check that the points are in the proper groups.
+	// Compute the pairing
+	res, err := pairing.Pair([]*sw_bw6761.G1Affine{&c.InG1}, []*sw_bw6761.G2Affine{&c.InG2})
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	pairing.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func ExamplePairing() {
+	p, q, err := randomG1G2Affines()
+	if err != nil {
+		panic(err)
+	}
+	res, err := bw6761.Pair([]bw6761.G1Affine{p}, []bw6761.G2Affine{q})
+	if err != nil {
+		panic(err)
+	}
+	circuit := PairCircuit{}
+	witness := PairCircuit{
+		InG1: sw_bw6761.NewG1Affine(p),
+		InG2: sw_bw6761.NewG2Affine(q),
+		Res:  sw_bw6761.NewGTEl(res),
+	}
+	ccs, err := frontend.Compile(ecc.BN254.ScalarField(), r1cs.NewBuilder, &circuit)
+	if err != nil {
+		panic(err)
+	} else {
+		fmt.Println("compiled")
+	}
+	pk, vk, err := groth16.Setup(ccs)
+	if err != nil {
+		panic(err)
+	} else {
+		fmt.Println("setup done")
+	}
+	secretWitness, err := frontend.NewWitness(&witness, ecc.BN254.ScalarField())
+	if err != nil {
+		panic(err)
+	} else {
+		fmt.Println("secret witness")
+	}
+	publicWitness, err := secretWitness.Public()
+	if err != nil {
+		panic(err)
+	} else {
+		fmt.Println("public witness")
+	}
+	proof, err := groth16.Prove(ccs, pk, secretWitness)
+	if err != nil {
+		panic(err)
+	} else {
+		fmt.Println("proof")
+	}
+	err = groth16.Verify(proof, vk, publicWitness)
+	if err != nil {
+		panic(err)
+	} else {
+		fmt.Println("verify")
+	}
+}
+
+func randomG1G2Affines() (p bw6761.G1Affine, q bw6761.G2Affine, err error) {
+	_, _, G1AffGen, G2AffGen := bw6761.Generators()
+	mod := bw6761.ID.ScalarField()
+	s1, err := rand.Int(rand.Reader, mod)
+	if err != nil {
+		return p, q, err
+	}
+	s2, err := rand.Int(rand.Reader, mod)
+	if err != nil {
+		return p, q, err
+	}
+	p.ScalarMultiplication(&G1AffGen, s1)
+	q.ScalarMultiplication(&G2AffGen, s2)
+	return
+}
diff --git a/std/algebra/emulated/sw_bw6761/g1.go b/std/algebra/emulated/sw_bw6761/g1.go
new file mode 100644
index 0000000000..0fe5ff2f40
--- /dev/null
+++ b/std/algebra/emulated/sw_bw6761/g1.go
@@ -0,0 +1,16 @@
+package sw_bw6761
+
+import (
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/std/algebra/emulated/sw_emulated"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+type G1Affine = sw_emulated.AffinePoint[emulated.BW6761Fp]
+
+func NewG1Affine(v bw6761.G1Affine) G1Affine {
+	return G1Affine{
+		X: emulated.ValueOf[emulated.BW6761Fp](v.X),
+		Y: emulated.ValueOf[emulated.BW6761Fp](v.Y),
+	}
+}
diff --git a/std/algebra/emulated/sw_bw6761/g2.go b/std/algebra/emulated/sw_bw6761/g2.go
new file mode 100644
index 0000000000..8de701405a
--- /dev/null
+++ b/std/algebra/emulated/sw_bw6761/g2.go
@@ -0,0 +1,16 @@
+package sw_bw6761
+
+import (
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/std/algebra/emulated/sw_emulated"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+type G2Affine = sw_emulated.AffinePoint[emulated.BW6761Fp]
+
+func NewG2Affine(v bw6761.G2Affine) G2Affine {
+	return G2Affine{
+		X: emulated.ValueOf[emulated.BW6761Fp](v.X),
+		Y: emulated.ValueOf[emulated.BW6761Fp](v.Y),
+	}
+}
diff --git a/std/algebra/emulated/sw_bw6761/pairing.go b/std/algebra/emulated/sw_bw6761/pairing.go
new file mode 100644
index 0000000000..80d719f7dc
--- /dev/null
+++ b/std/algebra/emulated/sw_bw6761/pairing.go
@@ -0,0 +1,514 @@
+package sw_bw6761
+
+import (
+	"errors"
+	"fmt"
+	"math/big"
+
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/emulated/fields_bw6761"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+type Pairing struct {
+	api frontend.API
+	*fields_bw6761.Ext6
+	curveF *emulated.Field[emulated.BW6761Fp]
+}
+
+type GTEl = fields_bw6761.E6
+
+func NewGTEl(v bw6761.GT) GTEl {
+	return GTEl{
+		B0: fields_bw6761.E3{
+			A0: emulated.ValueOf[emulated.BW6761Fp](v.B0.A0),
+			A1: emulated.ValueOf[emulated.BW6761Fp](v.B0.A1),
+			A2: emulated.ValueOf[emulated.BW6761Fp](v.B0.A2),
+		},
+		B1: fields_bw6761.E3{
+			A0: emulated.ValueOf[emulated.BW6761Fp](v.B1.A0),
+			A1: emulated.ValueOf[emulated.BW6761Fp](v.B1.A1),
+			A2: emulated.ValueOf[emulated.BW6761Fp](v.B1.A2),
+		},
+	}
+}
+
+func NewPairing(api frontend.API) (*Pairing, error) {
+	ba, err := emulated.NewField[emulated.BW6761Fp](api)
+	if err != nil {
+		return nil, fmt.Errorf("new base api: %w", err)
+	}
+	return &Pairing{
+		api:    api,
+		Ext6:   fields_bw6761.NewExt6(api),
+		curveF: ba,
+	}, nil
+}
+
+// FinalExponentiation computes the exponentiation zᵈ where
+//
+// d = (p⁶-1)/r = (p⁶-1)/Φ₆(p) ⋅ Φ₆(p)/r = (p³-1)(p+1)(p²-p+1)/r
+//
+// we use instead d = s⋅(p³-1)(p+1)(p²-p+1)/r
+// where s is the cofactor (x₀+1)
+func (pr Pairing) FinalExponentiation(z *GTEl) *GTEl {
+
+	z = pr.Reduce(z)
+	result := pr.Copy(z)
+
+	// 1. Easy part
+	// (p³-1)(p+1)
+	buf := pr.Conjugate(result)
+	buf = pr.DivUnchecked(buf, result)
+	result = pr.Frobenius(buf)
+	result = pr.Mul(result, buf)
+
+	// 2. Hard part (up to permutation)
+	// (x₀+1)(p²-p+1)/r
+	// Algorithm 4.4 from https://yelhousni.github.io/phd.pdf
+	a := pr.ExpX0Minus1Square(result)
+	a = pr.Mul(a, pr.Frobenius(result))
+	b := pr.ExpX0Plus1(a)
+	b = pr.Mul(b, pr.Conjugate(result))
+	t := pr.CyclotomicSquare(a)
+	a = pr.Mul(a, t)
+	c := pr.ExptMinus1Div3(b)
+	d := pr.ExpX0Minus1(c)
+	e := pr.ExpX0Minus1Square(d)
+	e = pr.Mul(e, d)
+	d = pr.Conjugate(d)
+	f := pr.Mul(d, b)
+	g := pr.ExpX0Plus1(e)
+	g = pr.Mul(g, f)
+	h := pr.Mul(g, c)
+	i := pr.Mul(g, d)
+	i = pr.ExpX0Plus1(i)
+	i = pr.Mul(i, pr.Conjugate(f))
+	j := pr.ExpC1(h)
+	j = pr.Mul(j, e)
+	k := pr.CyclotomicSquare(j)
+	k = pr.Mul(k, j)
+	k = pr.Mul(k, b)
+	t = pr.ExpC2(i)
+	k = pr.Mul(k, t)
+	result = pr.Mul(a, k)
+
+	return result
+}
+
+// lineEvaluation represents a sparse Fp6 Elmt (result of the line evaluation)
+// line: 1 + R0(x/y) + R1(1/y) = 0 instead of R0'*y + R1'*x + R2' = 0 This
+// makes the multiplication by lines (MulBy034)
+type lineEvaluation struct {
+	R0, R1 emulated.Element[emulated.BW6761Fp]
+}
+
+// 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 subgroup. See IsInSubGroup.
+func (pr Pairing) Pair(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
+	f, err := pr.MillerLoop(P, Q)
+	if err != nil {
+		return nil, err
+	}
+	return pr.FinalExponentiation(f), nil
+}
+
+// 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.
+func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
+	f, err := pr.Pair(P, Q)
+	if err != nil {
+		return err
+
+	}
+	one := pr.One()
+	pr.AssertIsEqual(f, one)
+
+	return nil
+}
+
+func (pr Pairing) AssertIsEqual(x, y *GTEl) {
+	pr.Ext6.AssertIsEqual(x, y)
+}
+
+// seed x₀=9586122913090633729
+//
+// x₀+1 in binary (64 bits) padded with 0s
+var loopCounter1 = [190]int8{
+	0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0,
+	0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+}
+
+// x₀³-x₀²-x₀ in 2-NAF
+var loopCounter2 = [190]int8{
+	-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, -1,
+	0, 1, 0, -1, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1,
+	0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0, 0, -1, 0, 0,
+	1, 0, 0, 0, -1, 0, 0, -1, 0, 1, 0, -1, 0, 0, 0, 1, 0, 0, 1, 0, -1, 0, 1, 0,
+	1, 0, 0, 0, 1, 0, -1, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 1,
+}
+
+// thirdRootOne² + thirdRootOne + 1 = 0 in BW6761Fp
+var thirdRootOne = emulated.ValueOf[emulated.BW6761Fp]("4922464560225523242118178942575080391082002530232324381063048548642823052024664478336818169867474395270858391911405337707247735739826664939444490469542109391530482826728203582549674992333383150446779312029624171857054392282775648")
+
+// MillerLoop computes the optimal Tate multi-Miller loop
+// (or twisted ate or Eta revisited)
+//
+// ∏ᵢ { fᵢ_{x₀+1+λ(x₀³-x₀²-x₀),Pᵢ}(Qᵢ) }
+//
+// Alg.2 in https://eprint.iacr.org/2021/1359.pdf
+// Eq. (6) in https://hackmd.io/@gnark/BW6-761-changes
+func (pr Pairing) MillerLoop(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
+
+	// check input size match
+	n := len(P)
+	if n == 0 || n != len(Q) {
+		return nil, errors.New("invalid inputs sizes")
+	}
+
+	// precomputations
+	negP := make([]*G1Affine, n)
+	imP := make([]*G1Affine, n)
+	imPneg := make([]*G1Affine, n)
+	accP := make([]*G1Affine, n)
+	yInv := make([]*emulated.Element[emulated.BW6761Fp], n)
+	xNegOverY := make([]*emulated.Element[emulated.BW6761Fp], n)
+
+	for k := 0; k < n; k++ {
+		// P and Q are supposed to be on G1 and G2 respectively of prime order r.
+		// The point (x,0) is of order 2. But this function does not check
+		// subgroup membership.
+		yInv[k] = pr.curveF.Inverse(&Q[k].Y)
+		xNegOverY[k] = pr.curveF.MulMod(&Q[k].X, yInv[k])
+		xNegOverY[k] = pr.curveF.Neg(xNegOverY[k])
+		// negP = -P = (x, -y)
+		negP[k] = &G1Affine{X: P[k].X, Y: *pr.curveF.Neg(&P[k].Y)}
+		// imP = (w*x, -y)
+		imP[k] = &G1Affine{X: *pr.curveF.MulMod(&P[k].X, &thirdRootOne), Y: negP[k].Y}
+		// imPneg = (w*x, y)
+		imPneg[k] = &G1Affine{X: imP[k].X, Y: P[k].Y}
+		// point accumulator initialized to imP
+		accP[k] = imP[k]
+	}
+
+	// f_{x₀+1+λ(x₀³-x₀²-x₀),P}(Q)
+	result := pr.Ext6.One()
+	var prodLines [5]*emulated.Element[emulated.BW6761Fp]
+	var l0, l1 *lineEvaluation
+
+	// i = 188, separately to avoid an E6 Square
+	// (Square(res) = 1² = 1)
+	// k = 0, separately to avoid MulBy034 (res × ℓ)
+	// (assign line to res)
+	accP[0], l0 = pr.doubleStep(imP[0])
+	result.B1 = fields_bw6761.E3{
+		A0: *pr.curveF.MulMod(&l0.R0, xNegOverY[0]),
+		A1: *pr.curveF.MulMod(&l0.R1, yInv[0]),
+		A2: result.B1.A2,
+	}
+
+	if n >= 2 {
+		// k = 1, separately to avoid MulBy034 (res × ℓ)
+		// (res is also a line at this point, so we use Mul034By034 ℓ × ℓ)
+		accP[1], l0 = pr.doubleStep(accP[1])
+		l0 = &lineEvaluation{
+			R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[1]),
+			R1: *pr.curveF.MulMod(&l0.R1, yInv[1]),
+		}
+		prodLines = pr.Mul034By034(&l0.R0, &l0.R1, &result.B1.A0, &result.B1.A1)
+		result = &fields_bw6761.E6{
+			B0: fields_bw6761.E3{
+				A0: *prodLines[0],
+				A1: *prodLines[1],
+				A2: *prodLines[2],
+			},
+			B1: fields_bw6761.E3{
+				A0: *prodLines[3],
+				A1: *prodLines[4],
+				A2: result.B1.A2,
+			},
+		}
+	}
+
+	if n >= 3 {
+		// k = 2, separately to avoid MulBy034 (res × ℓ)
+		// (res has a zero E2 element, so we use Mul01234By034)
+		accP[2], l0 = pr.doubleStep(accP[2])
+		l0 = &lineEvaluation{
+			R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[2]),
+			R1: *pr.curveF.MulMod(&l0.R1, yInv[2]),
+		}
+		result = pr.Mul01234By034(prodLines, &l0.R0, &l0.R1)
+
+		// k >= 3
+		for k := 3; k < n; k++ {
+			accP[k], l0 = pr.doubleStep(accP[k])
+			l0 = &lineEvaluation{
+				R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+				R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+			}
+			result = pr.MulBy034(result, &l0.R0, &l0.R1)
+		}
+	}
+
+	// i = 187
+	if n == 1 {
+		result = pr.Square034(result)
+	} else {
+		result = pr.Square(result)
+	}
+	for k := 0; k < n; k++ {
+		accP[k], l0 = pr.doubleStep(accP[k])
+		l0 = &lineEvaluation{
+			R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+			R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+		}
+		result = pr.MulBy034(result, &l0.R0, &l0.R1)
+	}
+
+	for i := 186; i >= 1; i-- {
+		// mutualize the square among n Miller loops
+		// (∏ᵢfᵢ)²
+		result = pr.Square(result)
+
+		j := loopCounter2[i]*3 + loopCounter1[i]
+
+		for k := 0; k < n; k++ {
+			switch j {
+			// cases -4, -2, 2 and 4 are omitted as they do not occur given the
+			// static loop counters.
+			case -3:
+				accP[k], l0, l1 = pr.doubleAndAddStep(accP[k], imPneg[k])
+				l0 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l0.R0, &l0.R1)
+				l1 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l1.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l1.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l1.R0, &l1.R1)
+			case -1:
+				accP[k], l0, l1 = pr.doubleAndAddStep(accP[k], negP[k])
+				l0 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l0.R0, &l0.R1)
+				l1 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l1.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l1.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l1.R0, &l1.R1)
+			case 0:
+				accP[k], l0 = pr.doubleStep(accP[k])
+				l0 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l0.R0, &l0.R1)
+			case 1:
+				accP[k], l0, l1 = pr.doubleAndAddStep(accP[k], P[k])
+				l0 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l0.R0, &l0.R1)
+				l1 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l1.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l1.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l1.R0, &l1.R1)
+			case 3:
+				accP[k], l0, l1 = pr.doubleAndAddStep(accP[k], imP[k])
+				l0 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l0.R0, &l0.R1)
+				l1 = &lineEvaluation{
+					R0: *pr.curveF.MulMod(&l1.R0, xNegOverY[k]),
+					R1: *pr.curveF.MulMod(&l1.R1, yInv[k]),
+				}
+				result = pr.MulBy034(result, &l1.R0, &l1.R1)
+			default:
+				return nil, errors.New("invalid loopCounter")
+			}
+		}
+	}
+
+	// i = 0, j = -3
+	// The resulting accumulator point is the infinity point because
+	// [(x₀+1) + λ(x₀³-x₀²-x₀)]P = [3(x₀-1)² ⋅ r]P = ∞
+	// since we're using affine coordinates, the addStep in the last iteration
+	// (j=-3) will fail as the slope of a vertical line in indefinite. But in
+	// projective coordinates, vertinal lines meet at (0:1:0) so the result
+	// should be unchanged if we ommit the addStep in this case. Moreover we
+	// just compute before the tangent line and not the full doubleStep as we
+	// only care about the Miller loop result in Fp6 and not the point itself.
+	result = pr.Square(result)
+	for k := 0; k < n; k++ {
+		l0 = pr.tangentCompute(accP[k])
+		l0 = &lineEvaluation{
+			R0: *pr.curveF.MulMod(&l0.R0, xNegOverY[k]),
+			R1: *pr.curveF.MulMod(&l0.R1, yInv[k]),
+		}
+		result = pr.MulBy034(result, &l0.R0, &l0.R1)
+	}
+
+	return result, nil
+
+}
+
+// addStep adds two points in affine coordinates, and evaluates the line in Miller loop
+// https://eprint.iacr.org/2022/1162 (Section 6.1)
+func (pr Pairing) addStep(p1, p2 *G1Affine) (*G1Affine, *lineEvaluation) {
+
+	// compute λ = (y2-y1)/(x2-x1)
+	p2ypy := pr.curveF.Sub(&p2.Y, &p1.Y)
+	p2xpx := pr.curveF.Sub(&p2.X, &p1.X)
+	λ := pr.curveF.Div(p2ypy, p2xpx)
+
+	// xr = λ²-x1-x2
+	λλ := pr.curveF.Mul(λ, λ)
+	p2xpx = pr.curveF.Add(&p1.X, &p2.X)
+	xr := pr.curveF.Sub(λλ, p2xpx)
+
+	// yr = λ(x1-xr) - y1
+	pxrx := pr.curveF.Sub(&p1.X, xr)
+	λpxrx := pr.curveF.Mul(λ, pxrx)
+	yr := pr.curveF.Sub(λpxrx, &p1.Y)
+
+	var res G2Affine
+	res.X = *xr
+	res.Y = *yr
+
+	var line lineEvaluation
+	line.R0 = *λ
+	line.R1 = *pr.curveF.Mul(λ, &p1.X)
+	line.R1 = *pr.curveF.Sub(&line.R1, &p1.Y)
+
+	return &res, &line
+
+}
+
+// doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates, and evaluates the line in Miller loop
+// https://eprint.iacr.org/2022/1162 (Section 6.1)
+func (pr Pairing) doubleAndAddStep(p1, p2 *G1Affine) (*G1Affine, *lineEvaluation, *lineEvaluation) {
+
+	var line1, line2 lineEvaluation
+	var p G1Affine
+
+	// compute λ1 = (y2-y1)/(x2-x1)
+	n := pr.curveF.Sub(&p1.Y, &p2.Y)
+	d := pr.curveF.Sub(&p1.X, &p2.X)
+	l1 := pr.curveF.Div(n, d)
+
+	// compute x3 =λ1²-x1-x2
+	x3 := pr.curveF.Mul(l1, l1)
+	x3 = pr.curveF.Sub(x3, &p1.X)
+	x3 = pr.curveF.Sub(x3, &p2.X)
+
+	// omit y3 computation
+
+	// compute line1
+	line1.R0 = *l1
+	line1.R1 = *pr.curveF.Mul(l1, &p1.X)
+	line1.R1 = *pr.curveF.Sub(&line1.R1, &p1.Y)
+
+	// compute λ2 = -λ1-2y1/(x3-x1)
+	n = pr.curveF.Add(&p1.Y, &p1.Y)
+	d = pr.curveF.Sub(x3, &p1.X)
+	l2 := pr.curveF.Div(n, d)
+	l2 = pr.curveF.Add(l2, l1)
+	l2 = pr.curveF.Neg(l2)
+
+	// compute x4 = λ2²-x1-x3
+	x4 := pr.curveF.Mul(l2, l2)
+	x4 = pr.curveF.Sub(x4, &p1.X)
+	x4 = pr.curveF.Sub(x4, x3)
+
+	// compute y4 = λ2(x1 - x4)-y1
+	y4 := pr.curveF.Sub(&p1.X, x4)
+	y4 = pr.curveF.Mul(l2, y4)
+	y4 = pr.curveF.Sub(y4, &p1.Y)
+
+	p.X = *x4
+	p.Y = *y4
+
+	// compute line2
+	line2.R0 = *l2
+	line2.R1 = *pr.curveF.Mul(l2, &p1.X)
+	line2.R1 = *pr.curveF.Sub(&line2.R1, &p1.Y)
+
+	return &p, &line1, &line2
+}
+
+// doubleStep doubles a point in affine coordinates, and evaluates the line in Miller loop
+// https://eprint.iacr.org/2022/1162 (Section 6.1)
+func (pr Pairing) doubleStep(p1 *G1Affine) (*G1Affine, *lineEvaluation) {
+
+	var p G1Affine
+	var line lineEvaluation
+
+	// λ = 3x²/2y
+	n := pr.curveF.Mul(&p1.X, &p1.X)
+	three := big.NewInt(3)
+	n = pr.curveF.MulConst(n, three)
+	d := pr.curveF.Add(&p1.Y, &p1.Y)
+	λ := pr.curveF.Div(n, d)
+
+	// xr = λ²-2x
+	xr := pr.curveF.Mul(λ, λ)
+	xr = pr.curveF.Sub(xr, &p1.X)
+	xr = pr.curveF.Sub(xr, &p1.X)
+
+	// yr = λ(x-xr)-y
+	yr := pr.curveF.Sub(&p1.X, xr)
+	yr = pr.curveF.Mul(λ, yr)
+	yr = pr.curveF.Sub(yr, &p1.Y)
+
+	p.X = *xr
+	p.Y = *yr
+
+	line.R0 = *λ
+	line.R1 = *pr.curveF.Mul(λ, &p1.X)
+	line.R1 = *pr.curveF.Sub(&line.R1, &p1.Y)
+
+	return &p, &line
+
+}
+
+// tangentCompute computes the line that goes through p1 and p2 but does not compute p1+p2
+func (pr Pairing) tangentCompute(p1 *G1Affine) *lineEvaluation {
+
+	// λ = 3x²/2y
+	n := pr.curveF.Mul(&p1.X, &p1.X)
+	three := big.NewInt(3)
+	n = pr.curveF.MulConst(n, three)
+	d := pr.curveF.Add(&p1.Y, &p1.Y)
+	λ := pr.curveF.Div(n, d)
+
+	var line lineEvaluation
+	line.R0 = *λ
+	line.R1 = *pr.curveF.Mul(λ, &p1.X)
+	line.R1 = *pr.curveF.Sub(&line.R1, &p1.Y)
+
+	return &line
+
+}
diff --git a/std/algebra/emulated/sw_bw6761/pairing_test.go b/std/algebra/emulated/sw_bw6761/pairing_test.go
new file mode 100644
index 0000000000..5e9ffbf867
--- /dev/null
+++ b/std/algebra/emulated/sw_bw6761/pairing_test.go
@@ -0,0 +1,245 @@
+package sw_bw6761
+
+import (
+	"bytes"
+	"crypto/rand"
+	"fmt"
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark/constraint"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/frontend/cs/r1cs"
+	"github.com/consensys/gnark/frontend/cs/scs"
+	"github.com/consensys/gnark/test"
+)
+
+func randomG1G2Affines() (bw6761.G1Affine, bw6761.G2Affine) {
+	_, _, G1AffGen, G2AffGen := bw6761.Generators()
+	mod := bw6761.ID.ScalarField()
+	s1, err := rand.Int(rand.Reader, mod)
+	if err != nil {
+		panic(err)
+	}
+	s2, err := rand.Int(rand.Reader, mod)
+	if err != nil {
+		panic(err)
+	}
+	var p bw6761.G1Affine
+	p.ScalarMultiplication(&G1AffGen, s1)
+	var q bw6761.G2Affine
+	q.ScalarMultiplication(&G2AffGen, s2)
+	return p, q
+}
+
+type FinalExponentiationCircuit struct {
+	InGt GTEl
+	Res  GTEl
+}
+
+func (c *FinalExponentiationCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	res := pairing.FinalExponentiation(&c.InGt)
+	pairing.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestFinalExponentiationTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	var gt bw6761.GT
+	gt.SetRandom()
+	res := bw6761.FinalExponentiation(&gt)
+	witness := FinalExponentiationCircuit{
+		InGt: NewGTEl(gt),
+		Res:  NewGTEl(res),
+	}
+	err := test.IsSolved(&FinalExponentiationCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type PairCircuit struct {
+	InG1 G1Affine
+	InG2 G2Affine
+	Res  GTEl
+}
+
+func (c *PairCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	res, err := pairing.Pair([]*G1Affine{&c.InG1}, []*G2Affine{&c.InG2})
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	pairing.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestPairTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p, q := randomG1G2Affines()
+	res, err := bw6761.Pair([]bw6761.G1Affine{p}, []bw6761.G2Affine{q})
+	assert.NoError(err)
+	witness := PairCircuit{
+		InG1: NewG1Affine(p),
+		InG2: NewG2Affine(q),
+		Res:  NewGTEl(res),
+	}
+	err = test.IsSolved(&PairCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type MultiPairCircuit struct {
+	InG1 G1Affine
+	InG2 G2Affine
+	Res  GTEl
+	n    int
+}
+
+func (c *MultiPairCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	P, Q := []*G1Affine{}, []*G2Affine{}
+	for i := 0; i < c.n; i++ {
+		P = append(P, &c.InG1)
+		Q = append(Q, &c.InG2)
+	}
+	res, err := pairing.Pair(P, Q)
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	pairing.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiPairTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p1, q1 := randomG1G2Affines()
+	p := make([]bw6761.G1Affine, 4)
+	q := make([]bw6761.G2Affine, 4)
+	for i := 0; i < 4; i++ {
+		p[i] = p1
+		q[i] = q1
+	}
+
+	for i := 2; i < 4; i++ {
+		res, err := bw6761.Pair(p[:i], q[:i])
+		assert.NoError(err)
+		witness := MultiPairCircuit{
+			InG1: NewG1Affine(p1),
+			InG2: NewG2Affine(q1),
+			Res:  NewGTEl(res),
+		}
+		err = test.IsSolved(&MultiPairCircuit{n: i}, &witness, ecc.BN254.ScalarField())
+		assert.NoError(err)
+	}
+}
+
+/*
+type PairingCheckCircuit struct {
+	In1G1 G1Affine
+	In2G1 G1Affine
+	In1G2 G2Affine
+	In2G2 G2Affine
+}
+
+func (c *PairingCheckCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	err = pairing.PairingCheck([]*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)
+	}
+	return nil
+}
+
+func TestPairingCheckTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p1, q1 := randomG1G2Affines()
+	_, q2 := randomG1G2Affines()
+	var p2 bw6761.G1Affine
+	p2.Neg(&p1)
+	witness := PairingCheckCircuit{
+		In1G1: NewG1Affine(p1),
+		In1G2: NewG2Affine(q1),
+		In2G1: NewG1Affine(p2),
+		In2G2: NewG2Affine(q2),
+	}
+	err := test.IsSolved(&PairingCheckCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+*/
+
+// bench
+func BenchmarkPairing(b *testing.B) {
+
+	p, q := randomG1G2Affines()
+	res, err := bw6761.Pair([]bw6761.G1Affine{p}, []bw6761.G2Affine{q})
+	if err != nil {
+		b.Fatal(err)
+	}
+	witness := PairCircuit{
+		InG1: NewG1Affine(p),
+		InG2: NewG2Affine(q),
+		Res:  NewGTEl(res),
+	}
+	w, err := frontend.NewWitness(&witness, ecc.BN254.ScalarField())
+	if err != nil {
+		b.Fatal(err)
+	}
+	var ccs constraint.ConstraintSystem
+	b.Run("compile scs", func(b *testing.B) {
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			if ccs, err = frontend.Compile(ecc.BN254.ScalarField(), scs.NewBuilder, &PairCircuit{}); err != nil {
+				b.Fatal(err)
+			}
+		}
+	})
+	var buf bytes.Buffer
+	_, err = ccs.WriteTo(&buf)
+	if err != nil {
+		b.Fatal(err)
+	}
+	b.Logf("scs size: %d (bytes), nb constraints %d, nbInstructions: %d", buf.Len(), ccs.GetNbConstraints(), ccs.GetNbInstructions())
+	b.Run("solve scs", func(b *testing.B) {
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			if _, err := ccs.Solve(w); err != nil {
+				b.Fatal(err)
+			}
+		}
+	})
+	b.Run("compile r1cs", func(b *testing.B) {
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			if ccs, err = frontend.Compile(ecc.BN254.ScalarField(), r1cs.NewBuilder, &PairCircuit{}); err != nil {
+				b.Fatal(err)
+			}
+		}
+	})
+	buf.Reset()
+	_, err = ccs.WriteTo(&buf)
+	if err != nil {
+		b.Fatal(err)
+	}
+	b.Logf("r1cs size: %d (bytes), nb constraints %d, nbInstructions: %d", buf.Len(), ccs.GetNbConstraints(), ccs.GetNbInstructions())
+
+	b.Run("solve r1cs", func(b *testing.B) {
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			if _, err := ccs.Solve(w); err != nil {
+				b.Fatal(err)
+			}
+		}
+	})
+}
diff --git a/std/math/emulated/emparams/emparams.go b/std/math/emulated/emparams/emparams.go
index 9b888d6a4d..01088ee435 100644
--- a/std/math/emulated/emparams/emparams.go
+++ b/std/math/emulated/emparams/emparams.go
@@ -31,6 +31,12 @@ func (sixLimbPrimeField) NbLimbs() uint     { return 6 }
 func (sixLimbPrimeField) BitsPerLimb() uint { return 64 }
 func (sixLimbPrimeField) IsPrime() bool     { return true }
 
+type twelveLimbPrimeField struct{}
+
+func (twelveLimbPrimeField) NbLimbs() uint     { return 12 }
+func (twelveLimbPrimeField) BitsPerLimb() uint { return 64 }
+func (twelveLimbPrimeField) IsPrime() bool     { return true }
+
 // Goldilocks provides type parametrization for field emulation:
 //   - limbs: 1
 //   - limb width: 64 bits
@@ -199,3 +205,31 @@ func (P384Fp) Modulus() *big.Int { return elliptic.P384().Params().P }
 type P384Fr struct{ sixLimbPrimeField }
 
 func (P384Fr) Modulus() *big.Int { return elliptic.P384().Params().N }
+
+// BW6761Fp provides type parametrization for field emulation:
+//   - limbs: 12
+//   - limb width: 64 bits
+//
+// The prime modulus for type parametrisation is:
+//
+//	0x122e824fb83ce0ad187c94004faff3eb926186a81d14688528275ef8087be41707ba638e584e91903cebaff25b423048689c8ed12f9fd9071dcd3dc73ebff2e98a116c25667a8f8160cf8aeeaf0a437e6913e6870000082f49d00000000008b (base 16)
+//	6891450384315732539396789682275657542479668912536150109513790160209623422243491736087683183289411687640864567753786613451161759120554247759349511699125301598951605099378508850372543631423596795951899700429969112842764913119068299 (base 10)
+//
+// This is the base field of the BW6-761 curve.
+type BW6761Fp struct{ twelveLimbPrimeField }
+
+func (fp BW6761Fp) Modulus() *big.Int { return ecc.BW6_761.BaseField() }
+
+// BW6761Fr provides type parametrization for field emulation:
+//   - limbs: 6
+//   - limb width: 64 bits
+//
+// The prime modulus for type parametrisation is:
+//
+//	0x1ae3a4617c510eac63b05c06ca1493b1a22d9f300f5138f1ef3622fba094800170b5d44300000008508c00000000001 (base 16)
+//	258664426012969094010652733694893533536393512754914660539884262666720468348340822774968888139573360124440321458177 (base 10)
+//
+// This is the scalar field of the BW6-761 curve.
+type BW6761Fr struct{ sixLimbPrimeField }
+
+func (fp BW6761Fr) Modulus() *big.Int { return ecc.BW6_761.ScalarField() }
diff --git a/std/math/emulated/params.go b/std/math/emulated/params.go
index fc3d8abfc1..892d141d38 100644
--- a/std/math/emulated/params.go
+++ b/std/math/emulated/params.go
@@ -37,4 +37,6 @@ type (
 	P256Fr      = emparams.P256Fr
 	P384Fp      = emparams.P384Fp
 	P384Fr      = emparams.P384Fr
+	BW6761Fp    = emparams.BW6761Fp
+	BW6761Fr    = emparams.BW6761Fr
 )