diff --git a/internal/stats/latest.stats b/internal/stats/latest.stats
index c156102ed..3278097e3 100644
Binary files a/internal/stats/latest.stats and b/internal/stats/latest.stats differ
diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go
index 3763c46d5..90fea92a7 100644
--- a/std/algebra/emulated/sw_bls12381/pairing.go
+++ b/std/algebra/emulated/sw_bls12381/pairing.go
@@ -442,7 +442,8 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 	return res, nil
 }
 
-// doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates, and evaluates the line in Miller loop
+// doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates.
+// Then evaluates the lines going through p1 and p2 (line1) and p1 and p1+p2 (line2).
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *lineEvaluation) {
 
@@ -452,47 +453,47 @@ func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *l
 	// compute λ1 = (y2-y1)/(x2-x1)
 	n := pr.Ext2.Sub(&p1.Y, &p2.Y)
 	d := pr.Ext2.Sub(&p1.X, &p2.X)
-	l1 := pr.Ext2.DivUnchecked(n, d)
+	λ1 := pr.Ext2.DivUnchecked(n, d)
 
 	// compute x3 =λ1²-x1-x2
-	x3 := pr.Ext2.Square(l1)
+	x3 := pr.Ext2.Square(λ1)
 	x3 = pr.Ext2.Sub(x3, pr.Ext2.Add(&p1.X, &p2.X))
 
 	// omit y3 computation
 
 	// compute line1
-	line1.R0 = *l1
-	line1.R1 = *pr.Ext2.Mul(l1, &p1.X)
+	line1.R0 = *λ1
+	line1.R1 = *pr.Ext2.Mul(λ1, &p1.X)
 	line1.R1 = *pr.Ext2.Sub(&line1.R1, &p1.Y)
 
 	// compute λ2 = -λ1-2y1/(x3-x1)
 	n = pr.Ext2.Double(&p1.Y)
 	d = pr.Ext2.Sub(x3, &p1.X)
-	l2 := pr.Ext2.DivUnchecked(n, d)
-	l2 = pr.Ext2.Add(l2, l1)
-	l2 = pr.Ext2.Neg(l2)
+	λ2 := pr.Ext2.DivUnchecked(n, d)
+	λ2 = pr.Ext2.Add(λ2, λ1)
+	λ2 = pr.Ext2.Neg(λ2)
 
 	// compute x4 = λ2²-x1-x3
-	x4 := pr.Ext2.Square(l2)
+	x4 := pr.Ext2.Square(λ2)
 	x4 = pr.Ext2.Sub(x4, pr.Ext2.Add(&p1.X, x3))
 
 	// compute y4 = λ2(x1 - x4)-y1
 	y4 := pr.Ext2.Sub(&p1.X, x4)
-	y4 = pr.Ext2.Mul(l2, y4)
+	y4 = pr.Ext2.Mul(λ2, y4)
 	y4 = pr.Ext2.Sub(y4, &p1.Y)
 
 	p.X = *x4
 	p.Y = *y4
 
 	// compute line2
-	line2.R0 = *l2
-	line2.R1 = *pr.Ext2.Mul(l2, &p1.X)
+	line2.R0 = *λ2
+	line2.R1 = *pr.Ext2.Mul(λ2, &p1.X)
 	line2.R1 = *pr.Ext2.Sub(&line2.R1, &p1.Y)
 
 	return &p, &line1, &line2
 }
 
-// doubleStep doubles a point in affine coordinates, and evaluates the line in Miller loop
+// doubleStep doubles p1 in affine coordinates, and evaluates the tangent line to p1.
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
@@ -575,7 +576,7 @@ func (pr Pairing) tripleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation, *lineEvaluat
 	return &res, &line1, &line2
 }
 
-// tangentCompute computes the line that goes through p1 and p2 but does not compute p1+p2
+// tangentCompute computes the tangent line to p1, but does not compute [2]p1.
 func (pr Pairing) tangentCompute(p1 *g2AffP) *lineEvaluation {
 
 	// λ = 3x²/2y
diff --git a/std/algebra/emulated/sw_bn254/pairing.go b/std/algebra/emulated/sw_bn254/pairing.go
index a2e998ad3..a2b6c941c 100644
--- a/std/algebra/emulated/sw_bn254/pairing.go
+++ b/std/algebra/emulated/sw_bn254/pairing.go
@@ -401,7 +401,7 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 		// The point (x,0) is of order 2. But this function does not check
 		// subgroup membership.
 		yInv[k] = pr.curveF.Inverse(&P[k].Y)
-		xNegOverY[k] = pr.curveF.MulMod(&P[k].X, yInv[k])
+		xNegOverY[k] = pr.curveF.Mul(&P[k].X, yInv[k])
 		xNegOverY[k] = pr.curveF.Neg(xNegOverY[k])
 	}
 
@@ -534,59 +534,63 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 	return res, nil
 }
 
-// doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates, and evaluates the line in Miller loop
+// doubleAndAddStep doubles p1 and adds or subs p2 to the result in affine coordinates, based on the isSub boolean.
+// Then evaluates the lines going through p1 and p2 or -p2 (line1) and p1 and p1+p2 or p1-p2 (line2).
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
-func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *lineEvaluation) {
+func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP, isSub bool) (*g2AffP, *lineEvaluation, *lineEvaluation) {
 
 	var line1, line2 lineEvaluation
 	var p g2AffP
 
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n := pr.Ext2.Sub(&p1.Y, &p2.Y)
+	// compute λ1 = (y1-y2)/(x1-x2) or λ1 = (y1+y2)/(x1-x2) if isSub is true
+	var n *fields_bn254.E2
+	if isSub {
+		n = pr.Ext2.Add(&p1.Y, &p2.Y)
+	} else {
+		n = pr.Ext2.Sub(&p1.Y, &p2.Y)
+	}
 	d := pr.Ext2.Sub(&p1.X, &p2.X)
-	l1 := pr.Ext2.DivUnchecked(n, d)
+	λ1 := pr.Ext2.DivUnchecked(n, d)
 
 	// 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.Square(λ1)
+	x3 = pr.Ext2.Sub(x3, pr.Ext2.Add(&p1.X, &p2.X))
 
 	// omit y3 computation
 
 	// compute line1
-	line1.R0 = *l1
-	line1.R1 = *pr.Ext2.Mul(l1, &p1.X)
+	line1.R0 = *λ1
+	line1.R1 = *pr.Ext2.Mul(λ1, &p1.X)
 	line1.R1 = *pr.Ext2.Sub(&line1.R1, &p1.Y)
 
 	// compute λ2 = -λ1-2y1/(x3-x1)
-	n = pr.Ext2.Double(&p1.Y)
+	n = pr.Ext2.MulByConstElement(&p1.Y, big.NewInt(2))
 	d = pr.Ext2.Sub(x3, &p1.X)
-	l2 := pr.Ext2.DivUnchecked(n, d)
-	l2 = pr.Ext2.Add(l2, l1)
-	l2 = pr.Ext2.Neg(l2)
+	λ2 := pr.Ext2.DivUnchecked(n, d)
+	λ2 = pr.Ext2.Add(λ2, λ1)
+	λ2 = pr.Ext2.Neg(λ2)
 
 	// 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.Square(λ2)
+	x4 = pr.Ext2.Sub(x4, pr.Ext2.Add(&p1.X, x3))
 
 	// compute y4 = λ2(x1 - x4)-y1
 	y4 := pr.Ext2.Sub(&p1.X, x4)
-	y4 = pr.Ext2.Mul(l2, y4)
+	y4 = pr.Ext2.Mul(λ2, y4)
 	y4 = pr.Ext2.Sub(y4, &p1.Y)
 
 	p.X = *x4
 	p.Y = *y4
 
 	// compute line2
-	line2.R0 = *l2
-	line2.R1 = *pr.Ext2.Mul(l2, &p1.X)
+	line2.R0 = *λ2
+	line2.R1 = *pr.Ext2.Mul(λ2, &p1.X)
 	line2.R1 = *pr.Ext2.Sub(&line2.R1, &p1.Y)
 
 	return &p, &line1, &line2
 }
 
-// doubleStep doubles a point in affine coordinates, and evaluates the line in Miller loop
+// doubleStep doubles p1 in affine coordinates, and evaluates the tangent line to p1.
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
@@ -595,15 +599,13 @@ func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
 	// λ = 3x²/2y
 	n := pr.Ext2.Square(&p1.X)
-	three := big.NewInt(3)
-	n = pr.Ext2.MulByConstElement(n, three)
-	d := pr.Ext2.Double(&p1.Y)
+	n = pr.Ext2.MulByConstElement(n, big.NewInt(3))
+	d := pr.Ext2.MulByConstElement(&p1.Y, big.NewInt(2))
 	λ := pr.Ext2.DivUnchecked(n, d)
 
 	// 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)
@@ -621,7 +623,7 @@ func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
 }
 
-// addStep adds two points in affine coordinates, and evaluates the line in Miller loop
+// addStep adds p1 and p2 in affine coordinates, and evaluates the line through p1 and p2.
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func (pr Pairing) addStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation) {
 
@@ -631,9 +633,8 @@ func (pr Pairing) addStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation) {
 	λ := pr.Ext2.DivUnchecked(p2ypy, p2xpx)
 
 	// xr = λ²-x1-x2
-	λλ := pr.Ext2.Square(λ)
-	p2xpx = pr.Ext2.Add(&p1.X, &p2.X)
-	xr := pr.Ext2.Sub(λλ, p2xpx)
+	xr := pr.Ext2.Square(λ)
+	xr = pr.Ext2.Sub(xr, pr.Ext2.Add(&p1.X, &p2.X))
 
 	// yr = λ(x1-xr) - y1
 	pxrx := pr.Ext2.Sub(&p1.X, xr)
@@ -653,12 +654,12 @@ func (pr Pairing) addStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation) {
 
 }
 
-// lineCompute computes the line that goes through p1 and p2 but does not compute p1+p2
+// lineCompute computes the line through p1 and p2, but does not compute p1+p2.
 func (pr Pairing) lineCompute(p1, p2 *g2AffP) *lineEvaluation {
 
-	// compute λ = (y2-y1)/(x2-x1)
-	qypy := pr.Ext2.Sub(&p2.Y, &p1.Y)
-	qxpx := pr.Ext2.Sub(&p2.X, &p1.X)
+	// compute λ = (y2+y1)/(x2-x1)
+	qypy := pr.Ext2.Add(&p1.Y, &p2.Y)
+	qxpx := pr.Ext2.Sub(&p1.X, &p2.X)
 	λ := pr.Ext2.DivUnchecked(qypy, qxpx)
 
 	var line lineEvaluation
@@ -727,7 +728,7 @@ func (pr Pairing) millerLoopAndFinalExpResult(P *G1Affine, Q *G2Affine, previous
 
 	// precomputations
 	yInv := pr.curveF.Inverse(&P.Y)
-	xNegOverY := pr.curveF.MulMod(&P.X, yInv)
+	xNegOverY := pr.curveF.Mul(&P.X, yInv)
 	xNegOverY = pr.curveF.Neg(xNegOverY)
 
 	// init Miller loop accumulator to residueWitnessInv to share the squarings
diff --git a/std/algebra/emulated/sw_bn254/precomputations.go b/std/algebra/emulated/sw_bn254/precomputations.go
index fe55542f4..84f8ffeab 100644
--- a/std/algebra/emulated/sw_bn254/precomputations.go
+++ b/std/algebra/emulated/sw_bn254/precomputations.go
@@ -33,22 +33,18 @@ func (p *Pairing) computeLines(Q *g2AffP) lineEvaluations {
 
 	var cLines lineEvaluations
 	Qacc := Q
-	QNeg := &g2AffP{
-		X: Q.X,
-		Y: *p.Ext2.Neg(&Q.Y),
-	}
 	n := len(bn254.LoopCounter)
 	Qacc, cLines[0][n-2] = p.doubleStep(Qacc)
-	cLines[1][n-3] = p.lineCompute(Qacc, QNeg)
+	cLines[1][n-3] = p.lineCompute(Qacc, Q)
 	Qacc, cLines[0][n-3] = p.addStep(Qacc, Q)
 	for i := n - 4; i >= 0; i-- {
 		switch loopCounter[i] {
 		case 0:
 			Qacc, cLines[0][i] = p.doubleStep(Qacc)
 		case 1:
-			Qacc, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(Qacc, Q)
+			Qacc, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(Qacc, Q, false)
 		case -1:
-			Qacc, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(Qacc, QNeg)
+			Qacc, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(Qacc, Q, true)
 		default:
 			return lineEvaluations{}
 		}
@@ -64,7 +60,6 @@ func (p *Pairing) computeLines(Q *g2AffP) lineEvaluations {
 	}
 
 	Q2Y := p.Ext2.MulByNonResidue2Power3(&Q.Y)
-	Q2Y = p.Ext2.Neg(Q2Y)
 	Q2 := &g2AffP{
 		X: *p.Ext2.MulByNonResidue2Power2(&Q.X),
 		Y: *Q2Y,
diff --git a/std/algebra/emulated/sw_bw6761/pairing.go b/std/algebra/emulated/sw_bw6761/pairing.go
index 47ad91556..6c9704bcc 100644
--- a/std/algebra/emulated/sw_bw6761/pairing.go
+++ b/std/algebra/emulated/sw_bw6761/pairing.go
@@ -383,15 +383,21 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 
 }
 
-// doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates, and evaluates the line in Miller loop
+// doubleAndAddStep doubles p1 and adds or subs p2 to the result in affine coordinates, based on the isSub boolean.
+// Then evaluates the lines going through p1 and p2 or -p2 (line1) and p1 and p1+p2 or p1-p2 (line2).
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
-func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *lineEvaluation) {
+func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP, isSub bool) (*g2AffP, *lineEvaluation, *lineEvaluation) {
 
 	var line1, line2 lineEvaluation
 	var p g2AffP
 
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n := pr.curveF.Sub(&p1.Y, &p2.Y)
+	// compute λ1 = (y1-y2)/(x1-x2) or λ1 = (y1+y2)/(x1-x2) if isSub is true
+	var n *emulated.Element[BaseField]
+	if isSub {
+		n = pr.curveF.Add(&p1.Y, &p2.Y)
+	} else {
+		n = pr.curveF.Sub(&p1.Y, &p2.Y)
+	}
 	d := pr.curveF.Sub(&p1.X, &p2.X)
 	l1 := pr.curveF.Div(n, d)
 
@@ -432,7 +438,57 @@ func (pr Pairing) doubleAndAddStep(p1, p2 *g2AffP) (*g2AffP, *lineEvaluation, *l
 	return &p, &line1, &line2
 }
 
-// doubleStep doubles a point in affine coordinates, and evaluates the line in Miller loop
+// 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))
+
+	// 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(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, 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
+}
+
+// doubleStep doubles p1 in affine coordinates, and evaluates the tangent line to p1.
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
@@ -465,7 +521,7 @@ func (pr Pairing) doubleStep(p1 *g2AffP) (*g2AffP, *lineEvaluation) {
 
 }
 
-// tangentCompute computes the line that goes through p1 and p2 but does not compute p1+p2
+// tangentCompute computes the tangent line to p1, but does not compute [2]p1.
 func (pr Pairing) tangentCompute(p1 *g2AffP) *lineEvaluation {
 
 	// λ = 3x²/2y
diff --git a/std/algebra/emulated/sw_bw6761/precomputations.go b/std/algebra/emulated/sw_bw6761/precomputations.go
index afb406012..829dd1b2b 100644
--- a/std/algebra/emulated/sw_bw6761/precomputations.go
+++ b/std/algebra/emulated/sw_bw6761/precomputations.go
@@ -35,31 +35,23 @@ func (p *Pairing) computeLines(Q *g2AffP) lineEvaluations {
 		X: *p.curveF.Mul(&Q.X, &thirdRootOne),
 		Y: *p.curveF.Neg(&Q.Y),
 	}
-	negQ := &g2AffP{
-		X: Q.X,
-		Y: imQ.Y,
-	}
 	accQ := &g2AffP{
 		X: imQ.X,
 		Y: imQ.Y,
 	}
-	imQneg := &g2AffP{
-		X: imQ.X,
-		Y: Q.Y,
-	}
 	for i := len(loopCounter2) - 2; i > 0; i-- {
 		switch loopCounter2[i]*3 + loopCounter1[i] {
 		// cases -4, -2, 2, 4 do not occur, given the static LoopCounters
 		case -3:
-			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, imQneg)
+			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, imQ, true)
 		case -1:
-			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, negQ)
+			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, Q, true)
 		case 0:
 			accQ, cLines[0][i] = p.doubleStep(accQ)
 		case 1:
-			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, Q)
+			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, Q, false)
 		case 3:
-			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, imQ)
+			accQ, cLines[0][i], cLines[1][i] = p.doubleAndAddStep(accQ, imQ, false)
 		default:
 			panic("unknown case for loopCounter")
 		}