diff --git a/internal/stats/latest_stats.csv b/internal/stats/latest_stats.csv
index 94706f3c86..18197f8970 100644
--- a/internal/stats/latest_stats.csv
+++ b/internal/stats/latest_stats.csv
@@ -195,56 +195,56 @@ pairing_bn254,bls24_315,plonk,0,0
 pairing_bn254,bls24_317,plonk,0,0
 pairing_bn254,bw6_761,plonk,0,0
 pairing_bn254,bw6_633,plonk,0,0
-pairing_bw6761,bn254,groth16,1794795,3003881
+pairing_bw6761,bn254,groth16,1782130,2981326
 pairing_bw6761,bls12_377,groth16,0,0
 pairing_bw6761,bls12_381,groth16,0,0
 pairing_bw6761,bls24_315,groth16,0,0
 pairing_bw6761,bls24_317,groth16,0,0
 pairing_bw6761,bw6_761,groth16,0,0
 pairing_bw6761,bw6_633,groth16,0,0
-pairing_bw6761,bn254,plonk,6779434,6155114
+pairing_bw6761,bn254,plonk,6696994,6075840
 pairing_bw6761,bls12_377,plonk,0,0
 pairing_bw6761,bls12_381,plonk,0,0
 pairing_bw6761,bls24_315,plonk,0,0
 pairing_bw6761,bls24_317,plonk,0,0
 pairing_bw6761,bw6_761,plonk,0,0
 pairing_bw6761,bw6_633,plonk,0,0
-scalar_mul_G1_bn254,bn254,groth16,69013,108022
+scalar_mul_G1_bn254,bn254,groth16,59287,91432
 scalar_mul_G1_bn254,bls12_377,groth16,0,0
 scalar_mul_G1_bn254,bls12_381,groth16,0,0
 scalar_mul_G1_bn254,bls24_315,groth16,0,0
 scalar_mul_G1_bn254,bls24_317,groth16,0,0
 scalar_mul_G1_bn254,bw6_761,groth16,0,0
 scalar_mul_G1_bn254,bw6_633,groth16,0,0
-scalar_mul_G1_bn254,bn254,plonk,260289,244439
+scalar_mul_G1_bn254,bn254,plonk,220730,207236
 scalar_mul_G1_bn254,bls12_377,plonk,0,0
 scalar_mul_G1_bn254,bls12_381,plonk,0,0
 scalar_mul_G1_bn254,bls24_315,plonk,0,0
 scalar_mul_G1_bn254,bls24_317,plonk,0,0
 scalar_mul_G1_bn254,bw6_761,plonk,0,0
 scalar_mul_G1_bn254,bw6_633,plonk,0,0
-scalar_mul_P256,bn254,groth16,93170,148354
+scalar_mul_P256,bn254,groth16,78854,124732
 scalar_mul_P256,bls12_377,groth16,0,0
 scalar_mul_P256,bls12_381,groth16,0,0
 scalar_mul_P256,bls24_315,groth16,0,0
 scalar_mul_P256,bls24_317,groth16,0,0
 scalar_mul_P256,bw6_761,groth16,0,0
 scalar_mul_P256,bw6_633,groth16,0,0
-scalar_mul_P256,bn254,plonk,355345,331788
+scalar_mul_P256,bn254,plonk,294014,274427
 scalar_mul_P256,bls12_377,plonk,0,0
 scalar_mul_P256,bls12_381,plonk,0,0
 scalar_mul_P256,bls24_315,plonk,0,0
 scalar_mul_P256,bls24_317,plonk,0,0
 scalar_mul_P256,bw6_761,plonk,0,0
 scalar_mul_P256,bw6_633,plonk,0,0
-scalar_mul_secp256k1,bn254,groth16,69860,109339
+scalar_mul_secp256k1,bn254,groth16,60025,92562
 scalar_mul_secp256k1,bls12_377,groth16,0,0
 scalar_mul_secp256k1,bls12_381,groth16,0,0
 scalar_mul_secp256k1,bls24_315,groth16,0,0
 scalar_mul_secp256k1,bls24_317,groth16,0,0
 scalar_mul_secp256k1,bw6_761,groth16,0,0
 scalar_mul_secp256k1,bw6_633,groth16,0,0
-scalar_mul_secp256k1,bn254,plonk,263180,247131
+scalar_mul_secp256k1,bn254,plonk,223490,209823
 scalar_mul_secp256k1,bls12_377,plonk,0,0
 scalar_mul_secp256k1,bls12_381,plonk,0,0
 scalar_mul_secp256k1,bls24_315,plonk,0,0
diff --git a/std/algebra/emulated/fields_bw6761/e6.go b/std/algebra/emulated/fields_bw6761/e6.go
index f43cdd28e3..125c902d26 100644
--- a/std/algebra/emulated/fields_bw6761/e6.go
+++ b/std/algebra/emulated/fields_bw6761/e6.go
@@ -136,24 +136,6 @@ func (e Ext6) Double(x *E6) *E6 {
 	}
 }
 
-func (e Ext6) MulByElement(x *E6, y *baseEl) *E6 {
-	a0 := e.fp.Mul(&x.A0, y)
-	a1 := e.fp.Mul(&x.A1, y)
-	a2 := e.fp.Mul(&x.A2, y)
-	a3 := e.fp.Mul(&x.A3, y)
-	a4 := e.fp.Mul(&x.A4, y)
-	a5 := e.fp.Mul(&x.A5, y)
-	z := &E6{
-		A0: *a0,
-		A1: *a1,
-		A2: *a2,
-		A3: *a3,
-		A4: *a4,
-		A5: *a5,
-	}
-	return z
-}
-
 func (e Ext6) MulByConstElement(x *E6, y *big.Int) *E6 {
 	a0 := e.fp.MulConst(&x.A0, y)
 	a1 := e.fp.MulConst(&x.A1, y)
diff --git a/std/algebra/emulated/fields_bw6761/e6_pairing.go b/std/algebra/emulated/fields_bw6761/e6_pairing.go
index 5750f7e63f..a4a40e840f 100644
--- a/std/algebra/emulated/fields_bw6761/e6_pairing.go
+++ b/std/algebra/emulated/fields_bw6761/e6_pairing.go
@@ -253,7 +253,7 @@ func (e Ext6) mul023by023Direct(d0, d1, c0, c1 *baseEl) [5]*baseEl {
 	// c3 = d0 + c0
 	z3 := e.fp.Add(d0, c0)
 	// c4 = d1c1
-	z4 := e.fp.Eval([][]*baseEl{{d1, c1}}, []int{1})
+	z4 := e.fp.Mul(d1, c1)
 	// c5 = d1 + c1,
 	z5 := e.fp.Add(d1, c1)
 
diff --git a/std/algebra/emulated/sw_bw6761/g1.go b/std/algebra/emulated/sw_bw6761/g1.go
index c31d3159ec..02278f42be 100644
--- a/std/algebra/emulated/sw_bw6761/g1.go
+++ b/std/algebra/emulated/sw_bw6761/g1.go
@@ -72,14 +72,11 @@ func (g1 *G1) double(p *G1Affine) *G1Affine {
 	λ := g1.curveF.Div(xx3a, y1)
 
 	// xr = λ²-2p.x
-	x1 := g1.curveF.MulConst(&p.X, big.NewInt(2))
-	λλ := g1.curveF.Mul(λ, λ)
-	xr := g1.curveF.Sub(λλ, x1)
+	mone := g1.curveF.NewElement(-1)
+	xr := g1.curveF.Eval([][]*baseEl{{λ, λ}, {mone, &p.X}}, []int{1, 2})
 
-	// yr = λ(p-xr) - p.y
-	pxrx := g1.curveF.Sub(&p.X, xr)
-	λpxrx := g1.curveF.Mul(λ, pxrx)
-	yr := g1.curveF.Sub(λpxrx, &p.Y)
+	// yr = λ(p.x-xr) - p.y
+	yr := g1.curveF.Eval([][]*baseEl{{λ, &p.X}, {mone, λ, xr}, {mone, &p.Y}}, []int{1, 1, 1})
 
 	return &G1Affine{
 		X: *xr,
@@ -102,14 +99,11 @@ func (g1 G1) add(p, q *G1Affine) *G1Affine {
 	λ := g1.curveF.Div(qypy, qxpx)
 
 	// xr = λ²-p.x-q.x
-	λλ := g1.curveF.Mul(λ, λ)
-	qxpx = g1.curveF.Add(&p.X, &q.X)
-	xr := g1.curveF.Sub(λλ, qxpx)
+	mone := g1.curveF.NewElement(-1)
+	xr := g1.curveF.Eval([][]*baseEl{{λ, λ}, {mone, &p.X}, {mone, &q.X}}, []int{1, 1, 1})
 
 	// p.y = λ(p.x-r.x) - p.y
-	pxrx := g1.curveF.Sub(&p.X, xr)
-	λpxrx := g1.curveF.Mul(λ, pxrx)
-	yr := g1.curveF.Sub(λpxrx, &p.Y)
+	yr := g1.curveF.Eval([][]*baseEl{{λ, &p.X}, {mone, λ, xr}, {mone, &p.Y}}, []int{1, 1, 1})
 
 	return &G1Affine{
 		X: *xr,
@@ -138,10 +132,9 @@ func (g1 G1) doubleAndAdd(p, q *G1Affine) *G1Affine {
 	xqxp := g1.curveF.Sub(&q.X, &p.X)
 	λ1 := g1.curveF.Div(yqyp, xqxp)
 
-	// compute x1 = λ1²-p.x-q.x
-	λ1λ1 := g1.curveF.Mul(λ1, λ1)
-	xqxp = g1.curveF.Add(&p.X, &q.X)
-	x2 := g1.curveF.Sub(λ1λ1, xqxp)
+	// compute x2 = λ1²-p.x-q.x
+	mone := g1.curveF.NewElement(-1)
+	x2 := g1.curveF.Eval([][]*baseEl{{λ1, λ1}, {mone, &p.X}, {mone, &q.X}}, []int{1, 1, 1})
 
 	// omit y1 computation
 	// compute λ1 = -λ1-1*p.y/(x1-p.x)
@@ -152,14 +145,10 @@ func (g1 G1) doubleAndAdd(p, q *G1Affine) *G1Affine {
 	λ2 = g1.curveF.Neg(λ2)
 
 	// compute x3 =λ2²-p.x-x3
-	λ2λ2 := g1.curveF.Mul(λ2, λ2)
-	x3 := g1.curveF.Sub(λ2λ2, &p.X)
-	x3 = g1.curveF.Sub(x3, x2)
+	x3 := g1.curveF.Eval([][]*baseEl{{λ2, λ2}, {mone, &p.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// compute y3 = λ2*(p.x - x3)-p.y
-	y3 := g1.curveF.Sub(&p.X, x3)
-	y3 = g1.curveF.Mul(λ2, y3)
-	y3 = g1.curveF.Sub(y3, &p.Y)
+	y3 := g1.curveF.Eval([][]*baseEl{{λ2, &p.X}, {mone, λ2, x3}, {mone, &p.Y}}, []int{1, 1, 1})
 
 	return &G1Affine{
 		X: *x3,
@@ -176,9 +165,8 @@ func (g1 G1) triple(p *G1Affine) *G1Affine {
 	λ1 := g1.curveF.Div(xx, y2)
 
 	// xr = λ²-2p.x
-	x2 := g1.curveF.MulConst(&p.X, big.NewInt(2))
-	λ1λ1 := g1.curveF.Mul(λ1, λ1)
-	x2 = g1.curveF.Sub(λ1λ1, x2)
+	mone := g1.curveF.NewElement(-1)
+	x2 := g1.curveF.Eval([][]*baseEl{{λ1, λ1}, {mone, &p.X}}, []int{1, 2})
 
 	// omit y2 computation, and
 	// compute λ2 = 2p.y/(x2 − p.x) − λ1.
@@ -187,14 +175,10 @@ func (g1 G1) triple(p *G1Affine) *G1Affine {
 	λ2 = g1.curveF.Sub(λ2, λ1)
 
 	// xr = λ²-p.x-x2
-	λ2λ2 := g1.curveF.Mul(λ2, λ2)
-	qxrx := g1.curveF.Add(x2, &p.X)
-	xr := g1.curveF.Sub(λ2λ2, qxrx)
+	xr := g1.curveF.Eval([][]*baseEl{{λ2, λ2}, {mone, &p.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// yr = λ(p.x-xr) - p.y
-	pxrx := g1.curveF.Sub(&p.X, xr)
-	λ2pxrx := g1.curveF.Mul(λ2, pxrx)
-	yr := g1.curveF.Sub(λ2pxrx, &p.Y)
+	yr := g1.curveF.Eval([][]*baseEl{{λ2, &p.X}, {mone, λ2, xr}, {mone, &p.Y}}, []int{1, 1, 1})
 
 	return &G1Affine{
 		X: *xr,
diff --git a/std/algebra/emulated/sw_bw6761/g2.go b/std/algebra/emulated/sw_bw6761/g2.go
index d1b2a4f752..0709eb5db9 100644
--- a/std/algebra/emulated/sw_bw6761/g2.go
+++ b/std/algebra/emulated/sw_bw6761/g2.go
@@ -110,14 +110,11 @@ func (g2 G2) add(p, q *G2Affine) *G2Affine {
 	λ := g2.curveF.Div(qypy, qxpx)
 
 	// xr = λ²-p.x-q.x
-	λλ := g2.curveF.Mul(λ, λ)
-	qxpx = g2.curveF.Add(&p.P.X, &q.P.X)
-	xr := g2.curveF.Sub(λλ, qxpx)
+	mone := g2.curveF.NewElement(-1)
+	xr := g2.curveF.Eval([][]*baseEl{{λ, λ}, {mone, &p.P.X}, {mone, &q.P.X}}, []int{1, 1, 1})
 
 	// p.y = λ(p.x-r.x) - p.y
-	pxrx := g2.curveF.Sub(&p.P.X, xr)
-	λpxrx := g2.curveF.Mul(λ, pxrx)
-	yr := g2.curveF.Sub(λpxrx, &p.P.Y)
+	yr := g2.curveF.Eval([][]*baseEl{{λ, &p.P.X}, {mone, λ, xr}, {mone, &p.P.Y}}, []int{1, 1, 1})
 
 	return &G2Affine{
 		P: g2AffP{
@@ -151,14 +148,11 @@ func (g2 *G2) double(p *G2Affine) *G2Affine {
 	λ := g2.curveF.Div(xx3a, y2)
 
 	// xr = λ²-2p.x
-	x2 := g2.curveF.MulConst(&p.P.X, big.NewInt(2))
-	λλ := g2.curveF.Mul(λ, λ)
-	xr := g2.curveF.Sub(λλ, x2)
+	mone := g2.curveF.NewElement(-1)
+	xr := g2.curveF.Eval([][]*baseEl{{λ, λ}, {mone, &p.P.X}}, []int{1, 2})
 
-	// yr = λ(p-xr) - p.y
-	pxrx := g2.curveF.Sub(&p.P.X, xr)
-	λpxrx := g2.curveF.Mul(λ, pxrx)
-	yr := g2.curveF.Sub(λpxrx, &p.P.Y)
+	// yr = λ(p.x-xr) - p.y
+	yr := g2.curveF.Eval([][]*baseEl{{λ, &p.P.X}, {mone, λ, xr}, {mone, &p.P.Y}}, []int{1, 1, 1})
 
 	return &G2Affine{
 		P: g2AffP{
@@ -184,9 +178,8 @@ func (g2 G2) doubleAndAdd(p, q *G2Affine) *G2Affine {
 	λ1 := g2.curveF.Div(yqyp, xqxp)
 
 	// compute x2 = λ1²-p.x-q.x
-	λ1λ1 := g2.curveF.Mul(λ1, λ1)
-	xqxp = g2.curveF.Add(&p.P.X, &q.P.X)
-	x2 := g2.curveF.Sub(λ1λ1, xqxp)
+	mone := g2.curveF.NewElement(-1)
+	x2 := g2.curveF.Eval([][]*baseEl{{λ1, λ1}, {mone, &p.P.X}, {mone, &q.P.X}}, []int{1, 1, 1})
 
 	// omit y2 computation
 	// compute λ2 = -λ1-2*p.y/(x2-p.x)
@@ -197,14 +190,10 @@ func (g2 G2) doubleAndAdd(p, q *G2Affine) *G2Affine {
 	λ2 = g2.curveF.Neg(λ2)
 
 	// compute x3 =λ2²-p.x-x3
-	λ2λ2 := g2.curveF.Mul(λ2, λ2)
-	x3 := g2.curveF.Sub(λ2λ2, &p.P.X)
-	x3 = g2.curveF.Sub(x3, x2)
+	x3 := g2.curveF.Eval([][]*baseEl{{λ2, λ2}, {mone, &p.P.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// compute y3 = λ2*(p.x - x3)-p.y
-	y3 := g2.curveF.Sub(&p.P.X, x3)
-	y3 = g2.curveF.Mul(λ2, y3)
-	y3 = g2.curveF.Sub(y3, &p.P.Y)
+	y3 := g2.curveF.Eval([][]*baseEl{{λ2, &p.P.X}, {mone, λ2, x3}, {mone, &p.P.Y}}, []int{1, 1, 1})
 
 	return &G2Affine{
 		P: g2AffP{
@@ -223,9 +212,8 @@ func (g2 G2) triple(p *G2Affine) *G2Affine {
 	λ1 := g2.curveF.Div(xx, y2)
 
 	// xr = λ²-2p.x
-	x2 := g2.curveF.MulConst(&p.P.X, big.NewInt(2))
-	λ1λ1 := g2.curveF.Mul(λ1, λ1)
-	x2 = g2.curveF.Sub(λ1λ1, x2)
+	mone := g2.curveF.NewElement(-1)
+	x2 := g2.curveF.Eval([][]*baseEl{{λ1, λ1}, {mone, &p.P.X}}, []int{1, 2})
 
 	// omit y2 computation, and
 	// compute λ2 = 2p.y/(x2 − p.x) − λ1.
@@ -234,14 +222,10 @@ func (g2 G2) triple(p *G2Affine) *G2Affine {
 	λ2 = g2.curveF.Sub(λ2, λ1)
 
 	// xr = λ²-p.x-x2
-	λ2λ2 := g2.curveF.Mul(λ2, λ2)
-	qxrx := g2.curveF.Add(x2, &p.P.X)
-	xr := g2.curveF.Sub(λ2λ2, qxrx)
+	xr := g2.curveF.Eval([][]*baseEl{{λ2, λ2}, {mone, &p.P.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// yr = λ(p.x-xr) - p.y
-	pxrx := g2.curveF.Sub(&p.P.X, xr)
-	λ2pxrx := g2.curveF.Mul(λ2, pxrx)
-	yr := g2.curveF.Sub(λ2pxrx, &p.P.Y)
+	yr := g2.curveF.Eval([][]*baseEl{{λ2, &p.P.X}, {mone, λ2, xr}, {mone, &p.P.Y}}, []int{1, 1, 1})
 
 	return &G2Affine{
 		P: g2AffP{
diff --git a/std/algebra/emulated/sw_bw6761/pairing.go b/std/algebra/emulated/sw_bw6761/pairing.go
index 6c9704bcc4..3b690be6fa 100644
--- a/std/algebra/emulated/sw_bw6761/pairing.go
+++ b/std/algebra/emulated/sw_bw6761/pairing.go
@@ -19,10 +19,10 @@ type Pairing struct {
 	curve  *sw_emulated.Curve[BaseField, ScalarField]
 	g1     *G1
 	g2     *G2
-	g2gen  *G2Affine
 }
 
 type GTEl = fields_bw6761.E6
+type baseEl = emulated.Element[BaseField]
 
 func NewGTEl(v bw6761.GT) GTEl {
 	return GTEl{
@@ -62,15 +62,6 @@ func NewPairing(api frontend.API) (*Pairing, error) {
 	}, nil
 }
 
-func (pr Pairing) generators() *G2Affine {
-	if pr.g2gen == nil {
-		_, _, _, g2gen := bw6761.Generators()
-		cg2gen := NewG2AffineFixed(g2gen)
-		pr.g2gen = &cg2gen
-	}
-	return pr.g2gen
-}
-
 // FinalExponentiation computes the exponentiation zᵈ where
 //
 // d = (p⁶-1)/r = (p⁶-1)/Φ₆(p) ⋅ Φ₆(p)/r = (p³-1)(p+1)(p²-p+1)/r
@@ -288,8 +279,8 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 	}
 
 	// precomputations
-	yInv := make([]*emulated.Element[BaseField], n)
-	xNegOverY := make([]*emulated.Element[BaseField], n)
+	yInv := make([]*baseEl, n)
+	xNegOverY := make([]*baseEl, n)
 
 	for k := 0; k < n; k++ {
 		// P are supposed to be on G1 respectively of prime order r.
@@ -301,7 +292,7 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 	}
 
 	// f_{x₀+1+λ(x₀³-x₀²-x₀),Q}(P), Q is known in advance
-	var prodLines [5]*emulated.Element[BaseField]
+	var prodLines [5]*baseEl
 	result := pr.Ext6.One()
 
 	// i = 188
@@ -390,6 +381,7 @@ func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP, isSub bool) (*g2AffP, *lineEv
 
 	var line1, line2 lineEvaluation
 	var p g2AffP
+	mone := pr.curveF.NewElement(-1)
 
 	// compute λ1 = (y1-y2)/(x1-x2) or λ1 = (y1+y2)/(x1-x2) if isSub is true
 	var n *emulated.Element[BaseField]
@@ -402,88 +394,32 @@ func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP, isSub bool) (*g2AffP, *lineEv
 	l1 := pr.curveF.Div(n, d)
 
 	// compute x3 =λ1²-x1-x2
-	x3 := pr.curveF.Mul(l1, l1)
-	x3 = pr.curveF.Sub(x3, pr.curveF.Add(&p1.X, &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.MulConst(&p1.Y, big.NewInt(2))
-	d = pr.curveF.Sub(&p1.X, x3)
-	l2 := pr.curveF.Div(n, d)
-	l2 = pr.curveF.Sub(l2, l1)
-
-	// compute x4 = λ2²-x1-x3
-	x4 := pr.curveF.Mul(l2, l2)
-	x4 = pr.curveF.Sub(x4, pr.curveF.Add(&p1.X, 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
-}
-
-// doubleAndSubStep doubles p1 and subs 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) doubleAndSubStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *lineEvaluation) {
-
-	var line1, line2 lineEvaluation
-	var p g2AffP
-
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n := pr.curveF.Add(&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, pr.curveF.Add(&p1.X, &p2.X))
+	x3 := pr.curveF.Eval([][]*baseEl{{l1, l1}, {mone, &p1.X}, {mone, &p2.X}}, []int{1, 1, 1})
 
 	// 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)
+	line1.R1 = *pr.curveF.Eval([][]*baseEl{{l1, &p1.X}, {mone, &p1.Y}}, []int{1, 1})
 
-	// compute λ2 = -λ1-2y1/(x3-x1)
-	n = pr.curveF.MulConst(&p1.Y, big.NewInt(2))
-	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 -λ2 = λ1+2y1/(x3-x1)
+	ypyp := pr.curveF.MulConst(&p1.Y, big.NewInt(2))
+	x2xp := pr.curveF.Sub(x3, &p1.X)
+	l2 := pr.curveF.Div(ypyp, x2xp)
+	l2 = pr.curveF.Add(l1, l2)
 
-	// compute x4 = λ2²-x1-x3
-	x4 := pr.curveF.Mul(l2, l2)
-	x4 = pr.curveF.Sub(x4, pr.curveF.Add(&p1.X, x3))
+	// compute x4 = (-λ2)²-x1-x3
+	x4 := pr.curveF.Eval([][]*baseEl{{l2, l2}, {mone, &p1.X}, {mone, x3}}, []int{1, 1, 1})
 
-	// 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)
+	// compute y4 = -λ2(-x1 + x4)-y1
+	y4 := pr.curveF.Eval([][]*baseEl{{l2, pr.curveF.Sub(x4, &p1.X)}, {mone, &p1.Y}}, []int{1, 1})
 
 	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)
+	line2.R0 = *pr.curveF.Neg(l2)
+	line2.R1 = *pr.curveF.Eval([][]*baseEl{{mone, l2, &p1.X}, {mone, &p1.Y}}, []int{1, 1})
 
 	return &p, &line1, &line2
 }
@@ -494,6 +430,7 @@ func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
 	var p g2AffP
 	var line lineEvaluation
+	mone := pr.curveF.NewElement(-1)
 
 	// λ = 3x²/2y
 	n := pr.curveF.Mul(&p1.X, &p1.X)
@@ -502,20 +439,16 @@ func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 	λ := pr.curveF.Div(n, d)
 
 	// xr = λ²-2x
-	xr := pr.curveF.Mul(λ, λ)
-	xr = pr.curveF.Sub(xr, pr.curveF.MulConst(&p1.X, big.NewInt(2)))
+	xr := pr.curveF.Eval([][]*baseEl{{λ, λ}, {mone, &p1.X}}, []int{1, 2})
 
 	// yr = λ(x-xr)-y
-	yr := pr.curveF.Sub(&p1.X, xr)
-	yr = pr.curveF.Mul(λ, yr)
-	yr = pr.curveF.Sub(yr, &p1.Y)
+	yr := pr.curveF.Eval([][]*baseEl{{λ, pr.curveF.Sub(&p1.X, xr)}, {mone, &p1.Y}}, []int{1, 1})
 
 	p.X = *xr
 	p.Y = *yr
 
 	line.R0 = *λ
-	line.R1 = *pr.curveF.Mul(λ, &p1.X)
-	line.R1 = *pr.curveF.Sub(&line.R1, &p1.Y)
+	line.R1 = *pr.curveF.Eval([][]*baseEl{{λ, &p1.X}, {mone, &p1.Y}}, []int{1, 1})
 
 	return &p, &line
 
@@ -532,8 +465,7 @@ func (pr Pairing) tangentCompute(p1 *g2AffP) *lineEvaluation {
 
 	var line lineEvaluation
 	line.R0 = *λ
-	line.R1 = *pr.curveF.Mul(λ, &p1.X)
-	line.R1 = *pr.curveF.Sub(&line.R1, &p1.Y)
+	line.R1 = *pr.curveF.Eval([][]*baseEl{{λ, &p1.X}, {pr.curveF.NewElement(-1), &p1.Y}}, []int{1, 1})
 
 	return &line
 
diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go
index 433300ccda..2cd9e71a6c 100644
--- a/std/algebra/emulated/sw_emulated/point.go
+++ b/std/algebra/emulated/sw_emulated/point.go
@@ -188,20 +188,17 @@ func (c *Curve[B, S]) AssertIsEqual(p, q *AffinePoint[B]) {
 //
 // It uses incomplete formulas in affine coordinates.
 func (c *Curve[B, S]) add(p, q *AffinePoint[B]) *AffinePoint[B] {
+	mone := c.baseApi.NewElement(-1)
 	// compute λ = (q.y-p.y)/(q.x-p.x)
 	qypy := c.baseApi.Sub(&q.Y, &p.Y)
 	qxpx := c.baseApi.Sub(&q.X, &p.X)
 	λ := c.baseApi.Div(qypy, qxpx)
 
 	// xr = λ²-p.x-q.x
-	λλ := c.baseApi.MulMod(λ, λ)
-	qxpx = c.baseApi.Add(&p.X, &q.X)
-	xr := c.baseApi.Sub(λλ, qxpx)
+	xr := c.baseApi.Eval([][]*emulated.Element[B]{{λ, λ}, {mone, c.baseApi.Add(&p.X, &q.X)}}, []int{1, 1})
 
 	// p.y = λ(p.x-r.x) - p.y
-	pxrx := c.baseApi.Sub(&p.X, xr)
-	λpxrx := c.baseApi.MulMod(λ, pxrx)
-	yr := c.baseApi.Sub(λpxrx, &p.Y)
+	yr := c.baseApi.Eval([][]*emulated.Element[B]{{λ, c.baseApi.Sub(&p.X, xr)}, {mone, &p.Y}}, []int{1, 1})
 
 	return &AffinePoint[B]{
 		X: *c.baseApi.Reduce(xr),
@@ -296,6 +293,7 @@ func (c *Curve[B, S]) Add(p, q *AffinePoint[B]) *AffinePoint[B] {
 // It uses affine coordinates.
 func (c *Curve[B, S]) double(p *AffinePoint[B]) *AffinePoint[B] {
 
+	mone := c.baseApi.NewElement(-1)
 	// compute λ = (3p.x²+a)/2*p.y, here we assume a=0 (j invariant 0 curve)
 	xx3a := c.baseApi.MulMod(&p.X, &p.X)
 	xx3a = c.baseApi.MulConst(xx3a, big.NewInt(3))
@@ -306,14 +304,10 @@ func (c *Curve[B, S]) double(p *AffinePoint[B]) *AffinePoint[B] {
 	λ := c.baseApi.Div(xx3a, y2)
 
 	// xr = λ²-2p.x
-	x2 := c.baseApi.MulConst(&p.X, big.NewInt(2))
-	λλ := c.baseApi.MulMod(λ, λ)
-	xr := c.baseApi.Sub(λλ, x2)
+	xr := c.baseApi.Eval([][]*emulated.Element[B]{{λ, λ}, {mone, &p.X}}, []int{1, 2})
 
 	// yr = λ(p-xr) - p.y
-	pxrx := c.baseApi.Sub(&p.X, xr)
-	λpxrx := c.baseApi.MulMod(λ, pxrx)
-	yr := c.baseApi.Sub(λpxrx, &p.Y)
+	yr := c.baseApi.Eval([][]*emulated.Element[B]{{λ, c.baseApi.Sub(&p.X, xr)}, {mone, &p.Y}}, []int{1, 1})
 
 	return &AffinePoint[B]{
 		X: *c.baseApi.Reduce(xr),
@@ -334,6 +328,7 @@ func (c *Curve[B, S]) double(p *AffinePoint[B]) *AffinePoint[B] {
 // [ELM03]: https://arxiv.org/pdf/math/0208038.pdf
 func (c *Curve[B, S]) triple(p *AffinePoint[B]) *AffinePoint[B] {
 
+	mone := c.baseApi.NewElement(-1)
 	// compute λ1 = (3p.x²+a)/2p.y, here we assume a=0 (j invariant 0 curve)
 	xx := c.baseApi.MulMod(&p.X, &p.X)
 	xx = c.baseApi.MulConst(xx, big.NewInt(3))
@@ -344,9 +339,7 @@ func (c *Curve[B, S]) triple(p *AffinePoint[B]) *AffinePoint[B] {
 	λ1 := c.baseApi.Div(xx, y2)
 
 	// xr = λ1²-2p.x
-	x2 := c.baseApi.MulConst(&p.X, big.NewInt(2))
-	λ1λ1 := c.baseApi.MulMod(λ1, λ1)
-	x2 = c.baseApi.Sub(λ1λ1, x2)
+	x2 := c.baseApi.Eval([][]*emulated.Element[B]{{λ1, λ1}, {mone, &p.X}}, []int{1, 2})
 
 	// omit y2 computation, and
 	// compute λ2 = 2p.y/(x2 − p.x) − λ1.
@@ -355,14 +348,10 @@ func (c *Curve[B, S]) triple(p *AffinePoint[B]) *AffinePoint[B] {
 	λ2 = c.baseApi.Sub(λ2, λ1)
 
 	// xr = λ²-p.x-x2
-	λ2λ2 := c.baseApi.MulMod(λ2, λ2)
-	qxrx := c.baseApi.Add(x2, &p.X)
-	xr := c.baseApi.Sub(λ2λ2, qxrx)
+	xr := c.baseApi.Eval([][]*emulated.Element[B]{{λ2, λ2}, {mone, &p.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// yr = λ(p.x-xr) - p.y
-	pxrx := c.baseApi.Sub(&p.X, xr)
-	λ2pxrx := c.baseApi.MulMod(λ2, pxrx)
-	yr := c.baseApi.Sub(λ2pxrx, &p.Y)
+	yr := c.baseApi.Eval([][]*emulated.Element[B]{{λ2, c.baseApi.Sub(&p.X, xr)}, {mone, &p.Y}}, []int{1, 1})
 
 	return &AffinePoint[B]{
 		X: *c.baseApi.Reduce(xr),
@@ -383,15 +372,14 @@ func (c *Curve[B, S]) triple(p *AffinePoint[B]) *AffinePoint[B] {
 // [ELM03]: https://arxiv.org/pdf/math/0208038.pdf
 func (c *Curve[B, S]) doubleAndAdd(p, q *AffinePoint[B]) *AffinePoint[B] {
 
+	mone := c.baseApi.NewElement(-1)
 	// compute λ1 = (q.y-p.y)/(q.x-p.x)
 	yqyp := c.baseApi.Sub(&q.Y, &p.Y)
 	xqxp := c.baseApi.Sub(&q.X, &p.X)
 	λ1 := c.baseApi.Div(yqyp, xqxp)
 
 	// compute x2 = λ1²-p.x-q.x
-	λ1λ1 := c.baseApi.MulMod(λ1, λ1)
-	xqxp = c.baseApi.Add(&p.X, &q.X)
-	x2 := c.baseApi.Sub(λ1λ1, xqxp)
+	x2 := c.baseApi.Eval([][]*emulated.Element[B]{{λ1, λ1}, {mone, c.baseApi.Add(&p.X, &q.X)}}, []int{1, 1})
 
 	// omit y2 computation
 	// compute λ2 = λ1+2*p.y/(x2-p.x)
@@ -401,13 +389,10 @@ func (c *Curve[B, S]) doubleAndAdd(p, q *AffinePoint[B]) *AffinePoint[B] {
 	λ2 = c.baseApi.Add(λ1, λ2)
 
 	// compute x3 =λ2²-p.x-x2
-	λ2λ2 := c.baseApi.MulMod(λ2, λ2)
-	x3 := c.baseApi.Sub(λ2λ2, c.baseApi.Add(&p.X, x2))
+	x3 := c.baseApi.Eval([][]*emulated.Element[B]{{λ2, λ2}, {mone, &p.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// compute y3 = λ2*(-p.x + x3)-p.y
-	y3 := c.baseApi.Sub(x3, &p.X)
-	y3 = c.baseApi.Mul(λ2, y3)
-	y3 = c.baseApi.Sub(y3, &p.Y)
+	y3 := c.baseApi.Eval([][]*emulated.Element[B]{{λ2, c.baseApi.Sub(x3, &p.X)}, {mone, &p.Y}}, []int{1, 1})
 
 	return &AffinePoint[B]{
 		X: *c.baseApi.Reduce(x3),
@@ -425,15 +410,14 @@ func (c *Curve[B, S]) doubleAndAdd(p, q *AffinePoint[B]) *AffinePoint[B] {
 // and then based on a Select adds either p or q.
 func (c *Curve[B, S]) doubleAndAddSelect(b frontend.Variable, p, q *AffinePoint[B]) *AffinePoint[B] {
 
+	mone := c.baseApi.NewElement(-1)
 	// compute λ1 = (q.y-p.y)/(q.x-p.x)
 	yqyp := c.baseApi.Sub(&q.Y, &p.Y)
 	xqxp := c.baseApi.Sub(&q.X, &p.X)
 	λ1 := c.baseApi.Div(yqyp, xqxp)
 
 	// compute x2 = λ1²-p.x-q.x
-	λ1λ1 := c.baseApi.MulMod(λ1, λ1)
-	xqxp = c.baseApi.Add(&p.X, &q.X)
-	x2 := c.baseApi.Sub(λ1λ1, xqxp)
+	x2 := c.baseApi.Eval([][]*emulated.Element[B]{{λ1, λ1}, {mone, &p.X}, {mone, &q.X}}, []int{1, 1, 1})
 
 	// omit y2 computation
 
@@ -447,13 +431,10 @@ func (c *Curve[B, S]) doubleAndAddSelect(b frontend.Variable, p, q *AffinePoint[
 	λ2 = c.baseApi.Add(λ1, λ2)
 
 	// compute x3 =λ2²-t.x-x2
-	λ2λ2 := c.baseApi.MulMod(λ2, λ2)
-	x3 := c.baseApi.Sub(λ2λ2, c.baseApi.Add(&t.X, x2))
+	x3 := c.baseApi.Eval([][]*emulated.Element[B]{{λ2, λ2}, {mone, &t.X}, {mone, x2}}, []int{1, 1, 1})
 
 	// compute y3 = -λ2*(t.x - x3)-t.y
-	y3 := c.baseApi.Sub(x3, &t.X)
-	y3 = c.baseApi.Mul(λ2, y3)
-	y3 = c.baseApi.Sub(y3, &t.Y)
+	y3 := c.baseApi.Eval([][]*emulated.Element[B]{{λ2, x3}, {mone, λ2, &t.X}, {mone, &t.Y}}, []int{1, 1, 1})
 
 	return &AffinePoint[B]{
 		X: *c.baseApi.Reduce(x3),
diff --git a/std/math/emulated/field_mul.go b/std/math/emulated/field_mul.go
index 224ffb4cc9..a1947cb7c3 100644
--- a/std/math/emulated/field_mul.go
+++ b/std/math/emulated/field_mul.go
@@ -838,7 +838,7 @@ func (mc *mvCheck[T]) cleanEvaluations() {
 // As it only depends on the bit-length of the inputs, then we can precompute it
 // regardless of the actual values.
 func (f *Field[T]) polyMvEvalQuoSize(mv *multivariate[T], at []*Element[T]) (quoSize int) {
-	modBits := f.fParams.Modulus().BitLen()
+	var fp T
 	quoSizes := make([]int, len(mv.Terms))
 	for i, term := range mv.Terms {
 		// for every term, the result length is the sum of the lengths of the
@@ -846,11 +846,11 @@ func (f *Field[T]) polyMvEvalQuoSize(mv *multivariate[T], at []*Element[T]) (quo
 		var lengths []int
 		for j, pow := range term {
 			for k := 0; k < pow; k++ {
-				lengths = append(lengths, modBits+int(at[j].overflow))
+				lengths = append(lengths, len(at[j].Limbs)*int(fp.BitsPerLimb())+int(at[j].overflow))
 			}
 		}
 		lengths = append(lengths, bits.Len(uint(mv.Coefficients[i])))
-		quoSizes[i] = sum(lengths...)
+		quoSizes[i] = sum(lengths...) - 1
 	}
 	// and for the full result, it is maximum of the inputs. We also add a bit
 	// for every term for overflow.