diff --git a/internal/stats/latest.stats b/internal/stats/latest.stats
index d961ee475d..f005c61d24 100644
Binary files a/internal/stats/latest.stats and b/internal/stats/latest.stats differ
diff --git a/std/algebra/emulated/fields_bls12381/e12.go b/std/algebra/emulated/fields_bls12381/e12.go
index d452c4f0cb..84728728c2 100644
--- a/std/algebra/emulated/fields_bls12381/e12.go
+++ b/std/algebra/emulated/fields_bls12381/e12.go
@@ -49,8 +49,8 @@ func (e Ext12) Mul(x, y *E12) *E12 {
 	a = e.Ext6.Mul(a, b)
 	b = e.Ext6.Mul(&x.C0, &y.C0)
 	c := e.Ext6.Mul(&x.C1, &y.C1)
-	z1 := e.Ext6.Sub(a, b)
-	z1 = e.Ext6.Sub(z1, c)
+	d := e.Ext6.Add(c, b)
+	z1 := e.Ext6.Sub(a, d)
 	z0 := e.Ext6.MulByNonResidue(c)
 	z0 = e.Ext6.Add(z0, b)
 	return &E12{
@@ -101,8 +101,7 @@ func (e Ext12) IsZero(z *E12) frontend.Variable {
 func (e Ext12) Square(x *E12) *E12 {
 	c0 := e.Ext6.Sub(&x.C0, &x.C1)
 	c3 := e.Ext6.MulByNonResidue(&x.C1)
-	c3 = e.Ext6.Neg(c3)
-	c3 = e.Ext6.Add(&x.C0, c3)
+	c3 = e.Ext6.Sub(&x.C0, c3)
 	c2 := e.Ext6.Mul(&x.C0, &x.C1)
 	c0 = e.Ext6.Mul(c0, c3)
 	c0 = e.Ext6.Add(c0, c2)
diff --git a/std/algebra/emulated/fields_bls12381/e12_pairing.go b/std/algebra/emulated/fields_bls12381/e12_pairing.go
index 69b0e9e37b..4cd3dc5613 100644
--- a/std/algebra/emulated/fields_bls12381/e12_pairing.go
+++ b/std/algebra/emulated/fields_bls12381/e12_pairing.go
@@ -93,8 +93,8 @@ func (e *Ext12) MulBy014(z *E12, c0, c1 *E2) *E12 {
 
 	zC1 := e.Ext6.Add(&z.C1, &z.C0)
 	zC1 = e.Ext6.MulBy01(zC1, c0, d)
-	zC1 = e.Ext6.Sub(zC1, a)
-	zC1 = e.Ext6.Sub(zC1, &b)
+	tmp := e.Ext6.Add(&b, a)
+	zC1 = e.Ext6.Sub(zC1, tmp)
 	zC0 := e.Ext6.MulByNonResidue(&b)
 	zC0 = e.Ext6.Add(zC0, a)
 
@@ -124,8 +124,8 @@ func (e Ext12) Mul014By014(d0, d1, c0, c1 *E2) [5]*E2 {
 	tmp := e.Ext2.Add(c0, c1)
 	x01 := e.Ext2.Add(d0, d1)
 	x01 = e.Ext2.Mul(x01, tmp)
-	x01 = e.Ext2.Sub(x01, x0)
-	x01 = e.Ext2.Sub(x01, x1)
+	tmp = e.Ext2.Add(x1, x0)
+	x01 = e.Ext2.Sub(x01, tmp)
 	x14 := e.Ext2.Add(c1, d1)
 
 	zC0B0 := e.Ext2.NonResidue()
@@ -148,8 +148,8 @@ func (e *Ext12) MulBy01245(z *E12, x [5]*E2) *E12 {
 	a = e.Ext6.Mul(a, b)
 	b = e.Ext6.Mul(&z.C0, c0)
 	c := e.Ext6.MulBy12(&z.C1, x[3], x[4])
-	z1 := e.Ext6.Sub(a, b)
-	z1 = e.Ext6.Sub(z1, c)
+	d := e.Ext6.Add(c, b)
+	z1 := e.Ext6.Sub(a, d)
 	z0 := e.Ext6.MulByNonResidue(c)
 	z0 = e.Ext6.Add(z0, b)
 	return &E12{
diff --git a/std/algebra/emulated/fields_bls12381/e2.go b/std/algebra/emulated/fields_bls12381/e2.go
index e36178b215..30b5ebda49 100644
--- a/std/algebra/emulated/fields_bls12381/e2.go
+++ b/std/algebra/emulated/fields_bls12381/e2.go
@@ -59,8 +59,8 @@ func NewExt2(api frontend.API) *Ext2 {
 }
 
 func (e Ext2) MulByElement(x *E2, y *baseEl) *E2 {
-	z0 := e.fp.MulMod(&x.A0, y)
-	z1 := e.fp.MulMod(&x.A1, y)
+	z0 := e.fp.Mul(&x.A0, y)
+	z1 := e.fp.Mul(&x.A1, y)
 	return &E2{
 		A0: *z0,
 		A1: *z1,
@@ -110,9 +110,9 @@ func (e Ext2) MulByNonResidue1Power1(x *E2) *E2 {
 // MulByNonResidue1Power2 returns x*(1+u)^(2*(p^1-1)/6)
 func (e Ext2) MulByNonResidue1Power2(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939436")
-	a := e.fp.MulMod(&x.A1, &element)
+	a := e.fp.Mul(&x.A1, &element)
 	a = e.fp.Neg(a)
-	b := e.fp.MulMod(&x.A0, &element)
+	b := e.fp.Mul(&x.A0, &element)
 	return &E2{
 		A0: *a,
 		A1: *b,
@@ -128,8 +128,8 @@ func (e Ext2) MulByNonResidue1Power3(x *E2) *E2 {
 func (e Ext2) MulByNonResidue1Power4(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939437")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -142,8 +142,8 @@ func (e Ext2) MulByNonResidue1Power5(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power1(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("793479390729215512621379701633421447060886740281060493010456487427281649075476305620758731620351")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -151,8 +151,8 @@ func (e Ext2) MulByNonResidue2Power1(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power2(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("793479390729215512621379701633421447060886740281060493010456487427281649075476305620758731620350")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -160,8 +160,8 @@ func (e Ext2) MulByNonResidue2Power2(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power3(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559786")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -169,8 +169,8 @@ func (e Ext2) MulByNonResidue2Power3(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power4(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939436")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -178,23 +178,26 @@ func (e Ext2) MulByNonResidue2Power4(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power5(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BLS12381Fp]("4002409555221667392624310435006688643935503118305586438271171395842971157480381377015405980053539358417135540939437")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
 func (e Ext2) Mul(x, y *E2) *E2 {
-	a := e.fp.Add(&x.A0, &x.A1)
-	b := e.fp.Add(&y.A0, &y.A1)
-	a = e.fp.MulMod(a, b)
-	b = e.fp.MulMod(&x.A0, &y.A0)
-	c := e.fp.MulMod(&x.A1, &y.A1)
-	z1 := e.fp.Sub(a, b)
-	z1 = e.fp.Sub(z1, c)
-	z0 := e.fp.Sub(b, c)
+
+	v0 := e.fp.Mul(&x.A0, &y.A0)
+	v1 := e.fp.Mul(&x.A1, &y.A1)
+
+	b0 := e.fp.Sub(v0, v1)
+	b1 := e.fp.Add(&x.A0, &x.A1)
+	tmp := e.fp.Add(&y.A0, &y.A1)
+	b1 = e.fp.Mul(b1, tmp)
+	tmp = e.fp.Add(v0, v1)
+	b1 = e.fp.Sub(b1, tmp)
+
 	return &E2{
-		A0: *z0,
-		A1: *z1,
+		A0: *b0,
+		A1: *b1,
 	}
 }
 
@@ -242,6 +245,7 @@ func (e Ext2) Zero() *E2 {
 		A1: *z1,
 	}
 }
+
 func (e Ext2) IsZero(z *E2) frontend.Variable {
 	a0 := e.fp.IsZero(&z.A0)
 	a1 := e.fp.IsZero(&z.A1)
@@ -260,8 +264,8 @@ func (e Ext2) NonResidue() *E2 {
 func (e Ext2) Square(x *E2) *E2 {
 	a := e.fp.Add(&x.A0, &x.A1)
 	b := e.fp.Sub(&x.A0, &x.A1)
-	a = e.fp.MulMod(a, b)
-	b = e.fp.MulMod(&x.A0, &x.A1)
+	a = e.fp.Mul(a, b)
+	b = e.fp.Mul(&x.A0, &x.A1)
 	b = e.fp.MulConst(b, big.NewInt(2))
 	return &E2{
 		A0: *a,
diff --git a/std/algebra/emulated/fields_bls12381/e6.go b/std/algebra/emulated/fields_bls12381/e6.go
index 406c317c0a..815af2262a 100644
--- a/std/algebra/emulated/fields_bls12381/e6.go
+++ b/std/algebra/emulated/fields_bls12381/e6.go
@@ -86,22 +86,22 @@ func (e Ext6) Mul(x, y *E6) *E6 {
 	c0 := e.Ext2.Add(&x.B1, &x.B2)
 	tmp := e.Ext2.Add(&y.B1, &y.B2)
 	c0 = e.Ext2.Mul(c0, tmp)
-	c0 = e.Ext2.Sub(c0, t1)
-	c0 = e.Ext2.Sub(c0, t2)
+	tmp = e.Ext2.Add(t2, t1)
+	c0 = e.Ext2.Sub(c0, tmp)
 	c0 = e.Ext2.MulByNonResidue(c0)
 	c0 = e.Ext2.Add(c0, t0)
 	c1 := e.Ext2.Add(&x.B0, &x.B1)
 	tmp = e.Ext2.Add(&y.B0, &y.B1)
 	c1 = e.Ext2.Mul(c1, tmp)
-	c1 = e.Ext2.Sub(c1, t0)
-	c1 = e.Ext2.Sub(c1, t1)
+	tmp = e.Ext2.Add(t0, t1)
+	c1 = e.Ext2.Sub(c1, tmp)
 	tmp = e.Ext2.MulByNonResidue(t2)
 	c1 = e.Ext2.Add(c1, tmp)
 	tmp = e.Ext2.Add(&x.B0, &x.B2)
 	c2 := e.Ext2.Add(&y.B0, &y.B2)
 	c2 = e.Ext2.Mul(c2, tmp)
-	c2 = e.Ext2.Sub(c2, t0)
-	c2 = e.Ext2.Sub(c2, t2)
+	tmp = e.Ext2.Add(t0, t2)
+	c2 = e.Ext2.Sub(c2, tmp)
 	c2 = e.Ext2.Add(c2, t1)
 	return &E6{
 		B0: *c0,
@@ -166,8 +166,8 @@ func (e Ext6) MulBy12(x *E6, b1, b2 *E2) *E6 {
 	c0 := e.Ext2.Add(&x.B1, &x.B2)
 	tmp := e.Ext2.Add(b1, b2)
 	c0 = e.Ext2.Mul(c0, tmp)
-	c0 = e.Ext2.Sub(c0, t1)
-	c0 = e.Ext2.Sub(c0, t2)
+	tmp = e.Ext2.Add(t1, t2)
+	c0 = e.Ext2.Sub(c0, tmp)
 	c0 = e.Ext2.MulByNonResidue(c0)
 	c1 := e.Ext2.Add(&x.B0, &x.B1)
 	c1 = e.Ext2.Mul(c1, b1)
@@ -207,7 +207,13 @@ func (e Ext6) MulBy0(z *E6, c0 *E2) *E6 {
 	}
 }
 
-// MulBy01 multiplication by sparse element (c0,c1,0)
+// MulBy01 multiplies z by an E6 sparse element of the form
+//
+//	E6{
+//		B0: c0,
+//		B1: c1,
+//		B2: 0,
+//	}
 func (e Ext6) MulBy01(z *E6, c0, c1 *E2) *E6 {
 	a := e.Ext2.Mul(&z.B0, c0)
 	b := e.Ext2.Mul(&z.B1, c1)
@@ -224,8 +230,8 @@ func (e Ext6) MulBy01(z *E6, c0, c1 *E2) *E6 {
 	t1 := e.Ext2.Add(c0, c1)
 	tmp = e.Ext2.Add(&z.B0, &z.B1)
 	t1 = e.Ext2.Mul(t1, tmp)
-	t1 = e.Ext2.Sub(t1, a)
-	t1 = e.Ext2.Sub(t1, b)
+	tmp = e.Ext2.Add(a, b)
+	t1 = e.Ext2.Sub(t1, tmp)
 	return &E6{
 		B0: *t0,
 		B1: *t1,
diff --git a/std/algebra/emulated/fields_bn254/e12.go b/std/algebra/emulated/fields_bn254/e12.go
index abe8bd9a9e..fb50d82676 100644
--- a/std/algebra/emulated/fields_bn254/e12.go
+++ b/std/algebra/emulated/fields_bn254/e12.go
@@ -49,8 +49,8 @@ func (e Ext12) Mul(x, y *E12) *E12 {
 	a = e.Ext6.Mul(a, b)
 	b = e.Ext6.Mul(&x.C0, &y.C0)
 	c := e.Ext6.Mul(&x.C1, &y.C1)
-	z1 := e.Ext6.Sub(a, b)
-	z1 = e.Ext6.Sub(z1, c)
+	d := e.Ext6.Add(c, b)
+	z1 := e.Ext6.Sub(a, d)
 	z0 := e.Ext6.MulByNonResidue(c)
 	z0 = e.Ext6.Add(z0, b)
 	return &E12{
@@ -101,8 +101,7 @@ func (e Ext12) IsZero(z *E12) frontend.Variable {
 func (e Ext12) Square(x *E12) *E12 {
 	c0 := e.Ext6.Sub(&x.C0, &x.C1)
 	c3 := e.Ext6.MulByNonResidue(&x.C1)
-	c3 = e.Ext6.Neg(c3)
-	c3 = e.Ext6.Add(&x.C0, c3)
+	c3 = e.Ext6.Sub(&x.C0, c3)
 	c2 := e.Ext6.Mul(&x.C0, &x.C1)
 	c0 = e.Ext6.Mul(c0, c3)
 	c0 = e.Ext6.Add(c0, c2)
diff --git a/std/algebra/emulated/fields_bn254/e12_pairing.go b/std/algebra/emulated/fields_bn254/e12_pairing.go
index 257f7d640f..9f896630ba 100644
--- a/std/algebra/emulated/fields_bn254/e12_pairing.go
+++ b/std/algebra/emulated/fields_bn254/e12_pairing.go
@@ -154,8 +154,8 @@ func (e *Ext12) Mul034By034(d3, d4, c3, c4 *E2) [5]*E2 {
 	tmp := e.Ext2.Add(c3, c4)
 	x34 := e.Ext2.Add(d3, d4)
 	x34 = e.Ext2.Mul(x34, tmp)
-	x34 = e.Ext2.Sub(x34, x3)
-	x34 = e.Ext2.Sub(x34, x4)
+	tmp = e.Ext2.Add(x4, x3)
+	x34 = e.Ext2.Sub(x34, tmp)
 
 	zC0B0 := e.Ext2.MulByNonResidue(x4)
 	zC0B0 = e.Ext2.Add(zC0B0, e.Ext2.One())
@@ -181,8 +181,8 @@ func (e *Ext12) MulBy01234(z *E12, x [5]*E2) *E12 {
 	a = e.Ext6.Mul(a, b)
 	b = e.Ext6.Mul(&z.C0, c0)
 	c := e.Ext6.MulBy01(&z.C1, x[3], x[4])
-	z1 := e.Ext6.Sub(a, b)
-	z1 = e.Ext6.Sub(z1, c)
+	d := e.Ext6.Add(c, b)
+	z1 := e.Ext6.Sub(a, d)
 	z0 := e.Ext6.MulByNonResidue(c)
 	z0 = e.Ext6.Add(z0, b)
 	return &E12{
@@ -211,8 +211,8 @@ func (e *Ext12) Mul01234By034(x [5]*E2, z3, z4 *E2) *E12 {
 	b := e.Ext6.Add(c0, c1)
 	a = e.Ext6.Mul(a, b)
 	c := e.Ext6.Mul01By01(z3, z4, x[3], x[4])
-	z1 := e.Ext6.Sub(a, c0)
-	z1 = e.Ext6.Sub(z1, c)
+	b = e.Ext6.Add(c0, c)
+	z1 := e.Ext6.Sub(a, b)
 	z0 := e.Ext6.MulByNonResidue(c)
 	z0 = e.Ext6.Add(z0, c0)
 	return &E12{
diff --git a/std/algebra/emulated/fields_bn254/e2.go b/std/algebra/emulated/fields_bn254/e2.go
index 161712f197..7533d7bf83 100644
--- a/std/algebra/emulated/fields_bn254/e2.go
+++ b/std/algebra/emulated/fields_bn254/e2.go
@@ -60,8 +60,8 @@ func NewExt2(api frontend.API) *Ext2 {
 }
 
 func (e Ext2) MulByElement(x *E2, y *baseEl) *E2 {
-	z0 := e.fp.MulMod(&x.A0, y)
-	z1 := e.fp.MulMod(&x.A1, y)
+	z0 := e.fp.Mul(&x.A0, y)
+	z1 := e.fp.Mul(&x.A1, y)
 	return &E2{
 		A0: *z0,
 		A1: *z1,
@@ -134,8 +134,8 @@ func (e Ext2) MulByNonResidue1Power5(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power1(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BN254Fp]("21888242871839275220042445260109153167277707414472061641714758635765020556617")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -143,8 +143,8 @@ func (e Ext2) MulByNonResidue2Power1(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power2(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BN254Fp]("21888242871839275220042445260109153167277707414472061641714758635765020556616")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -152,8 +152,8 @@ func (e Ext2) MulByNonResidue2Power2(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power3(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BN254Fp]("21888242871839275222246405745257275088696311157297823662689037894645226208582")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -161,8 +161,8 @@ func (e Ext2) MulByNonResidue2Power3(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power4(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BN254Fp]("2203960485148121921418603742825762020974279258880205651966")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -170,8 +170,8 @@ func (e Ext2) MulByNonResidue2Power4(x *E2) *E2 {
 func (e Ext2) MulByNonResidue2Power5(x *E2) *E2 {
 	element := emulated.ValueOf[emulated.BN254Fp]("2203960485148121921418603742825762020974279258880205651967")
 	return &E2{
-		A0: *e.fp.MulMod(&x.A0, &element),
-		A1: *e.fp.MulMod(&x.A1, &element),
+		A0: *e.fp.Mul(&x.A0, &element),
+		A1: *e.fp.Mul(&x.A1, &element),
 	}
 }
 
@@ -201,17 +201,20 @@ func (e Ext2) MulByNonResidue3Power5(x *E2) *E2 {
 }
 
 func (e Ext2) Mul(x, y *E2) *E2 {
-	a := e.fp.Add(&x.A0, &x.A1)
-	b := e.fp.Add(&y.A0, &y.A1)
-	a = e.fp.MulMod(a, b)
-	b = e.fp.MulMod(&x.A0, &y.A0)
-	c := e.fp.MulMod(&x.A1, &y.A1)
-	z1 := e.fp.Sub(a, b)
-	z1 = e.fp.Sub(z1, c)
-	z0 := e.fp.Sub(b, c)
+
+	v0 := e.fp.Mul(&x.A0, &y.A0)
+	v1 := e.fp.Mul(&x.A1, &y.A1)
+
+	b0 := e.fp.Sub(v0, v1)
+	b1 := e.fp.Add(&x.A0, &x.A1)
+	tmp := e.fp.Add(&y.A0, &y.A1)
+	b1 = e.fp.Mul(b1, tmp)
+	tmp = e.fp.Add(v0, v1)
+	b1 = e.fp.Sub(b1, tmp)
+
 	return &E2{
-		A0: *z0,
-		A1: *z1,
+		A0: *b0,
+		A1: *b1,
 	}
 }
 
@@ -269,8 +272,8 @@ func (e Ext2) IsZero(z *E2) frontend.Variable {
 func (e Ext2) Square(x *E2) *E2 {
 	a := e.fp.Add(&x.A0, &x.A1)
 	b := e.fp.Sub(&x.A0, &x.A1)
-	a = e.fp.MulMod(a, b)
-	b = e.fp.MulMod(&x.A0, &x.A1)
+	a = e.fp.Mul(a, b)
+	b = e.fp.Mul(&x.A0, &x.A1)
 	b = e.fp.MulConst(b, big.NewInt(2))
 	return &E2{
 		A0: *a,
diff --git a/std/algebra/emulated/fields_bn254/e6.go b/std/algebra/emulated/fields_bn254/e6.go
index 31ba7e0515..f7bf6158e2 100644
--- a/std/algebra/emulated/fields_bn254/e6.go
+++ b/std/algebra/emulated/fields_bn254/e6.go
@@ -86,22 +86,22 @@ func (e Ext6) Mul(x, y *E6) *E6 {
 	c0 := e.Ext2.Add(&x.B1, &x.B2)
 	tmp := e.Ext2.Add(&y.B1, &y.B2)
 	c0 = e.Ext2.Mul(c0, tmp)
-	c0 = e.Ext2.Sub(c0, t1)
-	c0 = e.Ext2.Sub(c0, t2)
+	tmp = e.Ext2.Add(t2, t1)
+	c0 = e.Ext2.Sub(c0, tmp)
 	c0 = e.Ext2.MulByNonResidue(c0)
 	c0 = e.Ext2.Add(c0, t0)
 	c1 := e.Ext2.Add(&x.B0, &x.B1)
 	tmp = e.Ext2.Add(&y.B0, &y.B1)
 	c1 = e.Ext2.Mul(c1, tmp)
-	c1 = e.Ext2.Sub(c1, t0)
-	c1 = e.Ext2.Sub(c1, t1)
+	tmp = e.Ext2.Add(t0, t1)
+	c1 = e.Ext2.Sub(c1, tmp)
 	tmp = e.Ext2.MulByNonResidue(t2)
 	c1 = e.Ext2.Add(c1, tmp)
 	tmp = e.Ext2.Add(&x.B0, &x.B2)
 	c2 := e.Ext2.Add(&y.B0, &y.B2)
 	c2 = e.Ext2.Mul(c2, tmp)
-	c2 = e.Ext2.Sub(c2, t0)
-	c2 = e.Ext2.Sub(c2, t2)
+	tmp = e.Ext2.Add(t0, t2)
+	c2 = e.Ext2.Sub(c2, tmp)
 	c2 = e.Ext2.Add(c2, t1)
 	return &E6{
 		B0: *c0,
@@ -204,8 +204,8 @@ func (e Ext6) MulBy01(z *E6, c0, c1 *E2) *E6 {
 	t1 := e.Ext2.Add(c0, c1)
 	tmp = e.Ext2.Add(&z.B0, &z.B1)
 	t1 = e.Ext2.Mul(t1, tmp)
-	t1 = e.Ext2.Sub(t1, a)
-	t1 = e.Ext2.Sub(t1, b)
+	tmp = e.Ext2.Add(a, b)
+	t1 = e.Ext2.Sub(t1, tmp)
 	return &E6{
 		B0: *t0,
 		B1: *t1,
@@ -234,8 +234,8 @@ func (e Ext6) Mul01By01(c0, c1, d0, d1 *E2) *E6 {
 	t1 := e.Ext2.Add(c0, c1)
 	tmp := e.Ext2.Add(d0, d1)
 	t1 = e.Ext2.Mul(t1, tmp)
-	t1 = e.Ext2.Sub(t1, a)
-	t1 = e.Ext2.Sub(t1, b)
+	tmp = e.Ext2.Add(a, b)
+	t1 = e.Ext2.Sub(t1, tmp)
 	return &E6{
 		B0: *a,
 		B1: *t1,
diff --git a/std/algebra/emulated/fields_bw6761/e3.go b/std/algebra/emulated/fields_bw6761/e3.go
index a2dc7530e9..b4ed19bddf 100644
--- a/std/algebra/emulated/fields_bw6761/e3.go
+++ b/std/algebra/emulated/fields_bw6761/e3.go
@@ -163,8 +163,8 @@ func (e Ext3) MulBy01(z *E3, c0, c1 *baseEl) *E3 {
 	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)
+	tmp = e.fp.Add(b, a)
+	t1 = e.fp.Sub(t1, tmp)
 
 	return &E3{
 		A0: *t0,
@@ -239,8 +239,8 @@ func (e Ext3) Mul01By01(c0, c1, d0, d1 *baseEl) *E3 {
 	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)
+	tmp = e.fp.Add(b, a)
+	t1 = e.fp.Sub(t1, tmp)
 	return &E3{
 		A0: *a,
 		A1: *t1,
@@ -249,6 +249,10 @@ func (e Ext3) Mul01By01(c0, c1, d0, d1 *baseEl) *E3 {
 }
 
 func (e Ext3) Mul(x, y *E3) *E3 {
+	return e.MulKaratsuba(x, y)
+}
+
+func (e Ext3) MulKaratsuba(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)
@@ -257,36 +261,118 @@ func (e Ext3) Mul(x, y *E3) *E3 {
 	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(t2, c0)
-	c0 = e.fp.MulConst(c0, big.NewInt(4))
-
-	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)
+	tmp = e.fp.Add(t2, t1)
+	c0 = e.fp.Sub(c0, tmp)
+	c0 = mulFpByNonResidue(e.fp, c0)
+	c0 = e.fp.Add(c0, t0)
 
 	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)
+	tmp = e.fp.MulConst(t2, big.NewInt(4))
+	tmp = e.fp.Add(tmp, t1)
+	tmp = e.fp.Add(tmp, t0)
+	c1 = e.fp.Sub(c1, tmp)
 
-	a0 := e.fp.Add(c0, t0)
-	a1 := e.fp.Add(c1, t2)
-	a2 := e.fp.Add(c2, t1)
+	c2 := e.fp.Add(&y.A0, &y.A2)
+	tmp = e.fp.Add(&x.A0, &x.A2)
+	c2 = e.fp.Mul(c2, tmp)
+	c2 = e.fp.Add(c2, t1)
+	tmp = e.fp.Add(t2, t0)
+	c2 = e.fp.Sub(c2, tmp)
 
 	return &E3{
-		A0: *a0,
-		A1: *a1,
-		A2: *a2,
+		A0: *c0,
+		A1: *c1,
+		A2: *c2,
 	}
 }
 
-func (e Ext3) Square(x *E3) *E3 {
+func (e Ext3) MulToomCook3(x, y *E3) *E3 {
+	// Toom-Cook-3x:
+	// We start by computing five interpolation points – these are evaluations of
+	// the product x(u)y(u) with u ∈ {0, ±1, 2, ∞}:
+	//
+	// v0 = x(0)y(0) = x.A0 * y.A0
+	// v1 = x(1)y(1) = (x.A0 + x.A1 + x.A2)(y.A0 + y.A1 + y.A2)
+	// v2 = x(−1)y(−1) = (x.A0 − x.A1 + x.A2)(y.A0 − y.A1 + y.A2)
+	// v3 = x(2)y(2) = (x.A0 + 2x.A1 + 4x.A2)(y.A0 + 2y.A1 + 4y.A2)
+	// v4 = x(∞)y(∞) = x.A2 * y.A2
+	//
+	// Then the interpolation is performed as:
+	//
+	// a0 = v0 + β((1/2)v0 − (1/2)v1 − (1/6)v2 + (1/6)v3 − 2v4)
+	// a1 = −(1/2)v0 + v1 − (1/3)v2 − (1/6)v3 + 2v4 + βv4
+	// a2 = −v0 + (1/2)v1 + (1/2)v2 − v4
+	//
+	// where is β=-4 the cubic non-residue (mulFpByNonResidue).
+	//
+	// In-circuit, we compute 6*x*y as
+	// a0 = 6v0 - β(3(v1 - v0 + 4v4) + v2 - v3)
+	// a1 = -(3v0 + 2v2 + v3) + 6(v1 + 2v4 + βv4)
+	// a2 = 3(v1 + v2 - 2(v0 + v4))
+	//
+	// and then divide a0, a1 and a2 by 6 using a hint.
+	//
+	// This costs 5M + 22A.
+
+	two := big.NewInt(2)
+	three := big.NewInt(3)
+	four := big.NewInt(4)
+	six := big.NewInt(6)
+
+	v0 := e.fp.Mul(&x.A0, &y.A0)
+	t1 := e.fp.Add(&x.A0, &x.A2)
+	t2 := e.fp.Add(&y.A0, &y.A2)
+	t3 := e.fp.Add(t2, &y.A1)
+	v1 := e.fp.Add(t1, &x.A1)
+	v1 = e.fp.Mul(v1, t3)
+	t3 = e.fp.Sub(t2, &y.A1)
+	v2 := e.fp.Sub(t1, &x.A1)
+	v2 = e.fp.Mul(v2, t3)
+	t1 = e.fp.MulConst(&x.A1, two)
+	t2 = e.fp.MulConst(&x.A2, four)
+	v3 := e.fp.Add(t1, t2)
+	v3 = e.fp.Add(v3, &x.A0)
+	t1 = e.fp.MulConst(&y.A1, two)
+	t2 = e.fp.MulConst(&y.A2, four)
+	t3 = e.fp.Add(t1, t2)
+	t3 = e.fp.Add(t3, &y.A0)
+	v3 = e.fp.Mul(v3, t3)
+	v4 := e.fp.Mul(&x.A2, &y.A2)
+
+	a0 := e.fp.Sub(v1, v0)
+	t1 = e.fp.MulConst(v4, four)
+	a0 = e.fp.Add(a0, t1)
+	a0 = e.fp.MulConst(a0, three)
+	a0 = e.fp.Sub(a0, v3)
+	a0 = e.fp.Add(a0, v2)
+	a0 = e.fp.MulConst(a0, four)
+	t1 = e.fp.MulConst(v0, six)
+	a0 = e.fp.Add(a0, t1)
+
+	t1 = e.fp.MulConst(v0, three)
+	t2 = e.fp.MulConst(v2, two)
+	t1 = e.fp.Add(t1, t2)
+	t1 = e.fp.Add(t1, v3)
+	a1 := e.fp.MulConst(v4, two)
+	a1 = e.fp.Sub(v1, a1)
+	a1 = e.fp.MulConst(a1, six)
+	a1 = e.fp.Sub(a1, t1)
+
+	t1 = e.fp.Add(v0, v4)
+	t1 = e.fp.MulConst(t1, two)
+	a2 := e.fp.Add(v1, v2)
+	a2 = e.fp.Sub(a2, t1)
+	a2 = e.fp.MulConst(a2, three)
+
+	return e.divE3By6(
+		&E3{A0: *a0, A1: *a1, A2: *a2},
+	)
+}
 
+func (e Ext3) Square(x *E3) *E3 {
+	// Chung-Hasan (SQR2)
 	// Algorithm 16 from https://eprint.iacr.org/2010/354.pdf
 
 	c6 := e.fp.MulConst(&x.A1, big.NewInt(2))
@@ -359,6 +445,26 @@ func (e Ext3) DivUnchecked(x, y *E3) *E3 {
 
 }
 
+func (e Ext3) divE3By6(x *E3) *E3 {
+	res, err := e.fp.NewHint(divE3By6Hint, 3, &x.A0, &x.A1, &x.A2)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	y := E3{
+		A0: *res[0],
+		A1: *res[1],
+		A2: *res[2],
+	}
+
+	// x == 6 * y
+	_x := e.MulByConstElement(&y, big.NewInt(6))
+	e.AssertIsEqual(x, _x)
+
+	return &y
+}
+
 // MulByNonResidue mul x by (0,1,0)
 func (e Ext3) MulByNonResidue(x *E3) *E3 {
 	z := &E3{
diff --git a/std/algebra/emulated/fields_bw6761/e6.go b/std/algebra/emulated/fields_bw6761/e6.go
index be781adcb4..2a78781e09 100644
--- a/std/algebra/emulated/fields_bw6761/e6.go
+++ b/std/algebra/emulated/fields_bw6761/e6.go
@@ -68,15 +68,16 @@ 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)
+	v0 := e.Ext3.Mul(&x.B0, &y.B0)
+	v1 := e.Ext3.Mul(&x.B1, &y.B1)
+
+	b0 := e.Ext3.MulByNonResidue(v1)
+	b0 = e.Ext3.Add(b0, v0)
+	b1 := e.Ext3.Add(&x.B0, &x.B1)
+	tmp := e.Ext3.Add(&y.B0, &y.B1)
+	b1 = e.Ext3.Mul(b1, tmp)
+	tmp = e.Ext3.Add(v0, v1)
+	b1 = e.Ext3.Sub(b1, tmp)
 
 	return &E6{
 		B0: *b0,
@@ -85,13 +86,18 @@ func (e Ext6) Mul(x, y *E6) *E6 {
 }
 
 func (e Ext6) Square(x *E6) *E6 {
-
-	x = e.Reduce(x)
 	//Algorithm 22 from https://eprint.iacr.org/2010/354.pdf
+	x = e.Reduce(x)
 	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)
+	c3 := &E3{
+		A0: x.B1.A2,
+		A1: x.B1.A0,
+		A2: x.B1.A1,
+	}
+	c3.A0 = *e.fp.MulConst(&c3.A0, big.NewInt(4))
+	c3.A0 = *e.fp.Add(&x.B0.A0, &c3.A0)
+	c3.A1 = *e.fp.Sub(&x.B0.A1, &c3.A1)
+	c3.A2 = *e.fp.Sub(&x.B0.A2, &c3.A2)
 	c2 := e.Ext3.Mul(&x.B0, &x.B1)
 	c0 = e.Ext3.Mul(c0, c3)
 	c0 = e.Ext3.Add(c0, c2)
@@ -370,14 +376,14 @@ func (e Ext6) CyclotomicSquare(x *E6) *E6 {
 	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
+	tmp := e.fp.Add(t[0], t[1])
+	t[6] = e.fp.Sub(t[6], tmp) // 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
+	tmp = e.fp.Add(t[2], t[3])
+	t[7] = e.fp.Sub(t[7], tmp) // 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)
diff --git a/std/algebra/emulated/fields_bw6761/e6_pairing.go b/std/algebra/emulated/fields_bw6761/e6_pairing.go
index 20e8a9c444..8e5f677b61 100644
--- a/std/algebra/emulated/fields_bw6761/e6_pairing.go
+++ b/std/algebra/emulated/fields_bw6761/e6_pairing.go
@@ -188,8 +188,8 @@ func (e Ext6) Mul014By014(d0, d1, c0, c1 *baseEl) [5]*baseEl {
 	tmp := e.fp.Add(c0, c1)
 	x01 := e.fp.Add(d0, d1)
 	x01 = e.fp.Mul(x01, tmp)
-	x01 = e.fp.Sub(x01, x0)
-	x01 = e.fp.Sub(x01, x1)
+	tmp = e.fp.Add(x0, x1)
+	x01 = e.fp.Sub(x01, tmp)
 	x14 := e.fp.Add(c1, d1)
 
 	four := emulated.ValueOf[emulated.BW6761Fp](big.NewInt(4))
diff --git a/std/algebra/emulated/fields_bw6761/hints.go b/std/algebra/emulated/fields_bw6761/hints.go
index 5c100d085c..879ea22206 100644
--- a/std/algebra/emulated/fields_bw6761/hints.go
+++ b/std/algebra/emulated/fields_bw6761/hints.go
@@ -4,6 +4,7 @@ import (
 	"math/big"
 
 	bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761"
+	"github.com/consensys/gnark-crypto/ecc/bw6-761/fp"
 	"github.com/consensys/gnark/constraint/solver"
 	"github.com/consensys/gnark/std/math/emulated"
 )
@@ -18,6 +19,7 @@ func GetHints() []solver.Hint {
 		// E3
 		divE3Hint,
 		inverseE3Hint,
+		divE3By6Hint,
 		// E6
 		divE6Hint,
 		inverseE6Hint,
@@ -44,6 +46,28 @@ func inverseE3Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) e
 		})
 }
 
+func divE3By6Hint(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])
+
+			var sixInv fp.Element
+			sixInv.SetString("6")
+			sixInv.Inverse(&sixInv)
+			c.MulByElement(&a, &sixInv)
+
+			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 {
diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go
index e180add6dc..bfe908cbd7 100644
--- a/std/algebra/emulated/sw_bls12381/pairing.go
+++ b/std/algebra/emulated/sw_bls12381/pairing.go
@@ -447,8 +447,7 @@ func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *l
 
 	// compute x3 =λ1²-x1-x2
 	x3 := pr.Ext2.Square(l1)
-	x3 = pr.Ext2.Sub(x3, &p1.X)
-	x3 = pr.Ext2.Sub(x3, &p2.X)
+	x3 = pr.Ext2.Sub(x3, pr.Ext2.Add(&p1.X, &p2.X))
 
 	// omit y3 computation
 
@@ -466,8 +465,7 @@ func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *l
 
 	// compute x4 = λ2²-x1-x3
 	x4 := pr.Ext2.Square(l2)
-	x4 = pr.Ext2.Sub(x4, &p1.X)
-	x4 = pr.Ext2.Sub(x4, x3)
+	x4 = pr.Ext2.Sub(x4, pr.Ext2.Add(&p1.X, x3))
 
 	// compute y4 = λ2(x1 - x4)-y1
 	y4 := pr.Ext2.Sub(&p1.X, x4)
@@ -501,8 +499,7 @@ func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
 	// xr = λ²-2x
 	xr := pr.Ext2.Square(λ)
-	xr = pr.Ext2.Sub(xr, &p1.X)
-	xr = pr.Ext2.Sub(xr, &p1.X)
+	xr = pr.Ext2.Sub(xr, pr.Ext2.MulByConstElement(&p1.X, big.NewInt(2)))
 
 	// yr = λ(x-xr)-y
 	yr := pr.Ext2.Sub(&p1.X, xr)
@@ -540,8 +537,7 @@ func (pr Pairing) tripleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation, *lineEvaluat
 
 	// x2 = λ1²-2x
 	x2 := pr.Ext2.Square(λ1)
-	x2 = pr.Ext2.Sub(x2, &p1.X)
-	x2 = pr.Ext2.Sub(x2, &p1.X)
+	x2 = pr.Ext2.Sub(x2, pr.Ext2.MulByConstElement(&p1.X, big.NewInt(2)))
 
 	// ommit yr computation, and
 	// compute λ2 = 2y/(x2 − x) − λ1.