diff --git a/ecc/bls12-377/pairing.go b/ecc/bls12-377/pairing.go
index bdef731e1c..a9c27df112 100644
--- a/ecc/bls12-377/pairing.go
+++ b/ecc/bls12-377/pairing.go
@@ -563,43 +563,73 @@ func (p *G2Affine) addStep(evaluations *LineEvaluationAff, a *G2Affine) {
 }
 
 func (p *G2Affine) doubleAndAddStep(evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
-	var n, d, l1, x3, l2, x4, y4 fptower.E2
+	var A, B, A2, B2, X2A2, t, U, AU, invAU, invA, invU, l1, x3, l2, x4, y4 fptower.E2
 
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n.Sub(&p.Y, &a.Y)
-	d.Sub(&p.X, &a.X)
-	l1.Div(&n, &d)
+	// The Eisenträger-Lauter-Montgomery formula for 2P+Q (https://eprint.iacr.org/2003/257)
+	// computes both slopes λ1 and λ2 using a single field inversion via batch inversion.
+	//
+	// Given P = (x1, y1) and Q = (x2, y2), let:
+	//   A = x1 - x2
+	//   B = y1 - y2
+	//   U = B² - (2x1 + x2)·A²
+	//
+	// Then:
+	//   λ1 = B/A                    (slope for P + Q)
+	//   λ2 = -λ1 - 2y1·A²/U         (slope for P + (P+Q))
+	//
+	// We compute 1/A and 1/U using Montgomery's batch inversion:
+	//   1/A = U/(A·U) and 1/U = A/(A·U) with a single inversion of A·U.
+
+	// Compute A = x1 - x2 and B = y1 - y2
+	A.Sub(&p.X, &a.X)
+	B.Sub(&p.Y, &a.Y)
+
+	// Compute A² and B²
+	A2.Square(&A)
+	B2.Square(&B)
 
-	// compute x3 =λ1²-x1-x2
+	// Compute U = B² - (2x1 + x2)·A²
+	t.Double(&p.X).Add(&t, &a.X)
+	X2A2.Mul(&t, &A2)
+	U.Sub(&B2, &X2A2)
+
+	// Batch inversion: compute 1/A and 1/U with a single inversion
+	AU.Mul(&A, &U)
+	invAU.Inverse(&AU)
+	invA.Mul(&U, &invAU)
+	invU.Mul(&A, &invAU)
+
+	// λ1 = B/A = B·(1/A)
+	l1.Mul(&B, &invA)
+
+	// x3 = λ1² - x1 - x2
 	x3.Square(&l1)
 	x3.Sub(&x3, &p.X)
 	x3.Sub(&x3, &a.X)
 
-	// omit y3 computation
-
-	// compute line1
+	// line1 evaluation
 	evaluations1.R0.Set(&l1)
 	evaluations1.R1.Mul(&l1, &p.X)
 	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
 
-	// compute λ2 = -λ1-2y1/(x3-x1)
-	n.Double(&p.Y)
-	d.Sub(&x3, &p.X)
-	l2.Div(&n, &d)
+	// λ2 = -λ1 - 2y1·A²/U = -λ1 - 2y1·A²·(1/U)
+	l2.Double(&p.Y)
+	l2.Mul(&l2, &A2)
+	l2.Mul(&l2, &invU)
 	l2.Add(&l2, &l1)
 	l2.Neg(&l2)
 
-	// compute x4 = λ2²-x1-x3
+	// x4 = λ2² - x1 - x3
 	x4.Square(&l2)
 	x4.Sub(&x4, &p.X)
 	x4.Sub(&x4, &x3)
 
-	// compute y4 = λ2(x1 - x4)-y1
+	// y4 = λ2·(x1 - x4) - y1
 	y4.Sub(&p.X, &x4)
 	y4.Mul(&l2, &y4)
 	y4.Sub(&y4, &p.Y)
 
-	// compute line2
+	// line2 evaluation
 	evaluations2.R0.Set(&l2)
 	evaluations2.R1.Mul(&l2, &p.X)
 	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
diff --git a/ecc/bls12-377/pairing_compatibility_test.go b/ecc/bls12-377/pairing_compatibility_test.go
new file mode 100644
index 0000000000..3fe36d68fa
--- /dev/null
+++ b/ecc/bls12-377/pairing_compatibility_test.go
@@ -0,0 +1,64 @@
+// Copyright 2020-2025 Consensys Software Inc.
+// Licensed under the Apache License, Version 2.0. See the LICENSE file for details.
+
+package bls12377
+
+import (
+	"github.com/consensys/gnark-crypto/ecc/bls12-377/internal/fptower"
+)
+
+// doubleAndAddStepRef is the reference (pre-optimization) implementation
+// of the doubleAndAddStep function. It computes 2P+Q using two field inversions.
+//
+// This version uses the standard chord-tangent method:
+//   - λ1 = (y2-y1)/(x2-x1) for P + Q
+//   - λ2 = -λ1 - 2y1/(x3-x1) for doubling and adding back
+//
+// The optimized version uses the Eisenträger-Lauter-Montgomery formula
+// (https://eprint.iacr.org/2003/257) which computes both slopes with a single
+// field inversion via Montgomery's batch inversion trick.
+func doubleAndAddStepRef(p *G2Affine, evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
+	var n, d, l1, x3, l2, x4, y4 fptower.E2
+
+	// compute λ1 = (y2-y1)/(x2-x1)
+	n.Sub(&p.Y, &a.Y)
+	d.Sub(&p.X, &a.X)
+	l1.Div(&n, &d)
+
+	// compute x3 =λ1²-x1-x2
+	x3.Square(&l1)
+	x3.Sub(&x3, &p.X)
+	x3.Sub(&x3, &a.X)
+
+	// omit y3 computation
+
+	// compute line1
+	evaluations1.R0.Set(&l1)
+	evaluations1.R1.Mul(&l1, &p.X)
+	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
+
+	// compute λ2 = -λ1-2y1/(x3-x1)
+	n.Double(&p.Y)
+	d.Sub(&x3, &p.X)
+	l2.Div(&n, &d)
+	l2.Add(&l2, &l1)
+	l2.Neg(&l2)
+
+	// compute x4 = λ2²-x1-x3
+	x4.Square(&l2)
+	x4.Sub(&x4, &p.X)
+	x4.Sub(&x4, &x3)
+
+	// compute y4 = λ2(x1 - x4)-y1
+	y4.Sub(&p.X, &x4)
+	y4.Mul(&l2, &y4)
+	y4.Sub(&y4, &p.Y)
+
+	// compute line2
+	evaluations2.R0.Set(&l2)
+	evaluations2.R1.Mul(&l2, &p.X)
+	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
+
+	p.X.Set(&x4)
+	p.Y.Set(&y4)
+}
diff --git a/ecc/bls12-377/pairing_test.go b/ecc/bls12-377/pairing_test.go
index 0bf62d72d5..0c81a1d92c 100644
--- a/ecc/bls12-377/pairing_test.go
+++ b/ecc/bls12-377/pairing_test.go
@@ -189,6 +189,67 @@ func TestPairing(t *testing.T) {
 	properties.TestingRun(t, gopter.ConsoleReporter(false))
 }
 
+func TestDoubleAndAddStepEquivalence(t *testing.T) {
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genR1 := GenFr()
+	genR2 := GenFr()
+
+	properties.Property("[BLS12-377] doubleAndAddStep: optimized (Eisenträger-Lauter-Montgomery) matches reference implementation", prop.ForAll(
+		func(s1, s2 fr.Element) bool {
+			var p1, p2, p1Ref G2Affine
+
+			// Get generator
+			_, _, _, g2Gen := Generators()
+
+			// Scale by random values
+			var s1Int, s2Int big.Int
+			s1.BigInt(&s1Int)
+			s2.BigInt(&s2Int)
+			p1.ScalarMultiplication(&g2Gen, &s1Int)
+			p2.ScalarMultiplication(&g2Gen, &s2Int)
+
+			// Fail if points are the same - this indicates a problem with random generation
+			if p1.X.Equal(&p2.X) {
+				return false
+			}
+
+			p1Ref.Set(&p1)
+
+			// Compute using optimized implementation
+			var eval1Opt, eval2Opt LineEvaluationAff
+			p1.doubleAndAddStep(&eval1Opt, &eval2Opt, &p2)
+
+			// Compute using reference implementation
+			var eval1Ref, eval2Ref LineEvaluationAff
+			doubleAndAddStepRef(&p1Ref, &eval1Ref, &eval2Ref, &p2)
+
+			// Compare results: resulting point coordinates
+			pointsEqual := p1.X.Equal(&p1Ref.X) && p1.Y.Equal(&p1Ref.Y)
+
+			// Compare results: line evaluation 1 (from first addition)
+			eval1Equal := eval1Opt.R0.Equal(&eval1Ref.R0) && eval1Opt.R1.Equal(&eval1Ref.R1)
+
+			// Compare results: line evaluation 2 (from doubling step)
+			eval2Equal := eval2Opt.R0.Equal(&eval2Ref.R0) && eval2Opt.R1.Equal(&eval2Ref.R1)
+
+			return pointsEqual && eval1Equal && eval2Equal
+		},
+		genR1,
+		genR2,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
 func TestMillerLoop(t *testing.T) {
 
 	t.Parallel()
@@ -404,6 +465,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/ecc/bls12-381/pairing.go b/ecc/bls12-381/pairing.go
index 3379bd25b1..33bc812b30 100644
--- a/ecc/bls12-381/pairing.go
+++ b/ecc/bls12-381/pairing.go
@@ -533,43 +533,73 @@ func (p *G2Affine) addStep(evaluations *LineEvaluationAff, a *G2Affine) {
 }
 
 func (p *G2Affine) doubleAndAddStep(evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
-	var n, d, l1, x3, l2, x4, y4 fptower.E2
+	var A, B, A2, B2, X2A2, t, U, AU, invAU, invA, invU, l1, x3, l2, x4, y4 fptower.E2
 
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n.Sub(&p.Y, &a.Y)
-	d.Sub(&p.X, &a.X)
-	l1.Div(&n, &d)
+	// The Eisenträger-Lauter-Montgomery formula for 2P+Q (https://eprint.iacr.org/2003/257)
+	// computes both slopes λ1 and λ2 using a single field inversion via batch inversion.
+	//
+	// Given P = (x1, y1) and Q = (x2, y2), let:
+	//   A = x1 - x2
+	//   B = y1 - y2
+	//   U = B² - (2x1 + x2)·A²
+	//
+	// Then:
+	//   λ1 = B/A                    (slope for P + Q)
+	//   λ2 = -λ1 - 2y1·A²/U         (slope for P + (P+Q))
+	//
+	// We compute 1/A and 1/U using Montgomery's batch inversion:
+	//   1/A = U/(A·U) and 1/U = A/(A·U) with a single inversion of A·U.
+
+	// Compute A = x1 - x2 and B = y1 - y2
+	A.Sub(&p.X, &a.X)
+	B.Sub(&p.Y, &a.Y)
+
+	// Compute A² and B²
+	A2.Square(&A)
+	B2.Square(&B)
 
-	// compute x3 =λ1²-x1-x2
+	// Compute U = B² - (2x1 + x2)·A²
+	t.Double(&p.X).Add(&t, &a.X)
+	X2A2.Mul(&t, &A2)
+	U.Sub(&B2, &X2A2)
+
+	// Batch inversion: compute 1/A and 1/U with a single inversion
+	AU.Mul(&A, &U)
+	invAU.Inverse(&AU)
+	invA.Mul(&U, &invAU)
+	invU.Mul(&A, &invAU)
+
+	// λ1 = B/A = B·(1/A)
+	l1.Mul(&B, &invA)
+
+	// x3 = λ1² - x1 - x2
 	x3.Square(&l1)
 	x3.Sub(&x3, &p.X)
 	x3.Sub(&x3, &a.X)
 
-	// omit y3 computation
-
-	// compute line1
+	// line1 evaluation
 	evaluations1.R0.Set(&l1)
 	evaluations1.R1.Mul(&l1, &p.X)
 	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
 
-	// compute λ2 = -λ1-2y1/(x3-x1)
-	n.Double(&p.Y)
-	d.Sub(&x3, &p.X)
-	l2.Div(&n, &d)
+	// λ2 = -λ1 - 2y1·A²/U = -λ1 - 2y1·A²·(1/U)
+	l2.Double(&p.Y)
+	l2.Mul(&l2, &A2)
+	l2.Mul(&l2, &invU)
 	l2.Add(&l2, &l1)
 	l2.Neg(&l2)
 
-	// compute x4 = λ2²-x1-x3
+	// x4 = λ2² - x1 - x3
 	x4.Square(&l2)
 	x4.Sub(&x4, &p.X)
 	x4.Sub(&x4, &x3)
 
-	// compute y4 = λ2(x1 - x4)-y1
+	// y4 = λ2·(x1 - x4) - y1
 	y4.Sub(&p.X, &x4)
 	y4.Mul(&l2, &y4)
 	y4.Sub(&y4, &p.Y)
 
-	// compute line2
+	// line2 evaluation
 	evaluations2.R0.Set(&l2)
 	evaluations2.R1.Mul(&l2, &p.X)
 	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
diff --git a/ecc/bls12-381/pairing_compatibility_test.go b/ecc/bls12-381/pairing_compatibility_test.go
new file mode 100644
index 0000000000..b710419fec
--- /dev/null
+++ b/ecc/bls12-381/pairing_compatibility_test.go
@@ -0,0 +1,64 @@
+// Copyright 2020-2025 Consensys Software Inc.
+// Licensed under the Apache License, Version 2.0. See the LICENSE file for details.
+
+package bls12381
+
+import (
+	"github.com/consensys/gnark-crypto/ecc/bls12-381/internal/fptower"
+)
+
+// doubleAndAddStepRef is the reference (pre-optimization) implementation
+// of the doubleAndAddStep function. It computes 2P+Q using two field inversions.
+//
+// This version uses the standard chord-tangent method:
+//   - λ1 = (y2-y1)/(x2-x1) for P + Q
+//   - λ2 = -λ1 - 2y1/(x3-x1) for doubling and adding back
+//
+// The optimized version uses the Eisenträger-Lauter-Montgomery formula
+// (https://eprint.iacr.org/2003/257) which computes both slopes with a single
+// field inversion via Montgomery's batch inversion trick.
+func doubleAndAddStepRef(p *G2Affine, evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
+	var n, d, l1, x3, l2, x4, y4 fptower.E2
+
+	// compute λ1 = (y2-y1)/(x2-x1)
+	n.Sub(&p.Y, &a.Y)
+	d.Sub(&p.X, &a.X)
+	l1.Div(&n, &d)
+
+	// compute x3 =λ1²-x1-x2
+	x3.Square(&l1)
+	x3.Sub(&x3, &p.X)
+	x3.Sub(&x3, &a.X)
+
+	// omit y3 computation
+
+	// compute line1
+	evaluations1.R0.Set(&l1)
+	evaluations1.R1.Mul(&l1, &p.X)
+	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
+
+	// compute λ2 = -λ1-2y1/(x3-x1)
+	n.Double(&p.Y)
+	d.Sub(&x3, &p.X)
+	l2.Div(&n, &d)
+	l2.Add(&l2, &l1)
+	l2.Neg(&l2)
+
+	// compute x4 = λ2²-x1-x3
+	x4.Square(&l2)
+	x4.Sub(&x4, &p.X)
+	x4.Sub(&x4, &x3)
+
+	// compute y4 = λ2(x1 - x4)-y1
+	y4.Sub(&p.X, &x4)
+	y4.Mul(&l2, &y4)
+	y4.Sub(&y4, &p.Y)
+
+	// compute line2
+	evaluations2.R0.Set(&l2)
+	evaluations2.R1.Mul(&l2, &p.X)
+	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
+
+	p.X.Set(&x4)
+	p.Y.Set(&y4)
+}
diff --git a/ecc/bls12-381/pairing_test.go b/ecc/bls12-381/pairing_test.go
index 065b07d75d..851c5a326d 100644
--- a/ecc/bls12-381/pairing_test.go
+++ b/ecc/bls12-381/pairing_test.go
@@ -189,6 +189,67 @@ func TestPairing(t *testing.T) {
 	properties.TestingRun(t, gopter.ConsoleReporter(false))
 }
 
+func TestDoubleAndAddStepEquivalence(t *testing.T) {
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genR1 := GenFr()
+	genR2 := GenFr()
+
+	properties.Property("[BLS12-381] doubleAndAddStep: optimized (Eisenträger-Lauter-Montgomery) matches reference implementation", prop.ForAll(
+		func(s1, s2 fr.Element) bool {
+			var p1, p2, p1Ref G2Affine
+
+			// Get generator
+			_, _, _, g2Gen := Generators()
+
+			// Scale by random values
+			var s1Int, s2Int big.Int
+			s1.BigInt(&s1Int)
+			s2.BigInt(&s2Int)
+			p1.ScalarMultiplication(&g2Gen, &s1Int)
+			p2.ScalarMultiplication(&g2Gen, &s2Int)
+
+			// Fail if points are the same - this indicates a problem with random generation
+			if p1.X.Equal(&p2.X) {
+				return false
+			}
+
+			p1Ref.Set(&p1)
+
+			// Compute using optimized implementation
+			var eval1Opt, eval2Opt LineEvaluationAff
+			p1.doubleAndAddStep(&eval1Opt, &eval2Opt, &p2)
+
+			// Compute using reference implementation
+			var eval1Ref, eval2Ref LineEvaluationAff
+			doubleAndAddStepRef(&p1Ref, &eval1Ref, &eval2Ref, &p2)
+
+			// Compare results: resulting point coordinates
+			pointsEqual := p1.X.Equal(&p1Ref.X) && p1.Y.Equal(&p1Ref.Y)
+
+			// Compare results: line evaluation 1 (from first addition)
+			eval1Equal := eval1Opt.R0.Equal(&eval1Ref.R0) && eval1Opt.R1.Equal(&eval1Ref.R1)
+
+			// Compare results: line evaluation 2 (from doubling step)
+			eval2Equal := eval2Opt.R0.Equal(&eval2Ref.R0) && eval2Opt.R1.Equal(&eval2Ref.R1)
+
+			return pointsEqual && eval1Equal && eval2Equal
+		},
+		genR1,
+		genR2,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
 func TestMillerLoop(t *testing.T) {
 
 	t.Parallel()
@@ -404,6 +465,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/ecc/bls24-315/pairing_test.go b/ecc/bls24-315/pairing_test.go
index 3070b3958a..c4af521e01 100644
--- a/ecc/bls24-315/pairing_test.go
+++ b/ecc/bls24-315/pairing_test.go
@@ -405,6 +405,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/ecc/bls24-317/pairing_test.go b/ecc/bls24-317/pairing_test.go
index 99d2a92aed..3a36b77b36 100644
--- a/ecc/bls24-317/pairing_test.go
+++ b/ecc/bls24-317/pairing_test.go
@@ -405,6 +405,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/ecc/bn254/pairing.go b/ecc/bn254/pairing.go
index 080db40f88..af4b3451cd 100644
--- a/ecc/bn254/pairing.go
+++ b/ecc/bn254/pairing.go
@@ -666,43 +666,73 @@ func (p *G2Affine) addStep(evaluations *LineEvaluationAff, a *G2Affine) {
 }
 
 func (p *G2Affine) doubleAndAddStep(evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
-	var n, d, l1, x3, l2, x4, y4 fptower.E2
+	var A, B, A2, B2, X2A2, t, U, AU, invAU, invA, invU, l1, x3, l2, x4, y4 fptower.E2
 
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n.Sub(&p.Y, &a.Y)
-	d.Sub(&p.X, &a.X)
-	l1.Div(&n, &d)
+	// The Eisenträger-Lauter-Montgomery formula for 2P+Q (https://eprint.iacr.org/2003/257)
+	// computes both slopes λ1 and λ2 using a single field inversion via batch inversion.
+	//
+	// Given P = (x1, y1) and Q = (x2, y2), let:
+	//   A = x1 - x2
+	//   B = y1 - y2
+	//   U = B² - (2x1 + x2)·A²
+	//
+	// Then:
+	//   λ1 = B/A                    (slope for P + Q)
+	//   λ2 = -λ1 - 2y1·A²/U         (slope for P + (P+Q))
+	//
+	// We compute 1/A and 1/U using Montgomery's batch inversion:
+	//   1/A = U/(A·U) and 1/U = A/(A·U) with a single inversion of A·U.
+
+	// Compute A = x1 - x2 and B = y1 - y2
+	A.Sub(&p.X, &a.X)
+	B.Sub(&p.Y, &a.Y)
+
+	// Compute A² and B²
+	A2.Square(&A)
+	B2.Square(&B)
 
-	// compute x3 =λ1²-x1-x2
+	// Compute U = B² - (2x1 + x2)·A²
+	t.Double(&p.X).Add(&t, &a.X)
+	X2A2.Mul(&t, &A2)
+	U.Sub(&B2, &X2A2)
+
+	// Batch inversion: compute 1/A and 1/U with a single inversion
+	AU.Mul(&A, &U)
+	invAU.Inverse(&AU)
+	invA.Mul(&U, &invAU)
+	invU.Mul(&A, &invAU)
+
+	// λ1 = B/A = B·(1/A)
+	l1.Mul(&B, &invA)
+
+	// x3 = λ1² - x1 - x2
 	x3.Square(&l1)
 	x3.Sub(&x3, &p.X)
 	x3.Sub(&x3, &a.X)
 
-	// omit y3 computation
-
-	// compute line1
+	// line1 evaluation
 	evaluations1.R0.Set(&l1)
 	evaluations1.R1.Mul(&l1, &p.X)
 	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
 
-	// compute λ2 = -λ1-2y1/(x3-x1)
-	n.Double(&p.Y)
-	d.Sub(&x3, &p.X)
-	l2.Div(&n, &d)
+	// λ2 = -λ1 - 2y1·A²/U = -λ1 - 2y1·A²·(1/U)
+	l2.Double(&p.Y)
+	l2.Mul(&l2, &A2)
+	l2.Mul(&l2, &invU)
 	l2.Add(&l2, &l1)
 	l2.Neg(&l2)
 
-	// compute x4 = λ2²-x1-x3
+	// x4 = λ2² - x1 - x3
 	x4.Square(&l2)
 	x4.Sub(&x4, &p.X)
 	x4.Sub(&x4, &x3)
 
-	// compute y4 = λ2(x1 - x4)-y1
+	// y4 = λ2·(x1 - x4) - y1
 	y4.Sub(&p.X, &x4)
 	y4.Mul(&l2, &y4)
 	y4.Sub(&y4, &p.Y)
 
-	// compute line2
+	// line2 evaluation
 	evaluations2.R0.Set(&l2)
 	evaluations2.R1.Mul(&l2, &p.X)
 	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
diff --git a/ecc/bn254/pairing_compatibility_test.go b/ecc/bn254/pairing_compatibility_test.go
new file mode 100644
index 0000000000..37a3d0d348
--- /dev/null
+++ b/ecc/bn254/pairing_compatibility_test.go
@@ -0,0 +1,64 @@
+// Copyright 2020-2025 Consensys Software Inc.
+// Licensed under the Apache License, Version 2.0. See the LICENSE file for details.
+
+package bn254
+
+import (
+	"github.com/consensys/gnark-crypto/ecc/bn254/internal/fptower"
+)
+
+// doubleAndAddStepRef is the reference (pre-optimization) implementation
+// of the doubleAndAddStep function. It computes 2P+Q using two field inversions.
+//
+// This version uses the standard chord-tangent method:
+//   - λ1 = (y2-y1)/(x2-x1) for P + Q
+//   - λ2 = -λ1 - 2y1/(x3-x1) for doubling and adding back
+//
+// The optimized version uses the Eisenträger-Lauter-Montgomery formula
+// (https://eprint.iacr.org/2003/257) which computes both slopes with a single
+// field inversion via Montgomery's batch inversion trick.
+func doubleAndAddStepRef(p *G2Affine, evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
+	var n, d, l1, x3, l2, x4, y4 fptower.E2
+
+	// compute λ1 = (y2-y1)/(x2-x1)
+	n.Sub(&p.Y, &a.Y)
+	d.Sub(&p.X, &a.X)
+	l1.Div(&n, &d)
+
+	// compute x3 =λ1²-x1-x2
+	x3.Square(&l1)
+	x3.Sub(&x3, &p.X)
+	x3.Sub(&x3, &a.X)
+
+	// omit y3 computation
+
+	// compute line1
+	evaluations1.R0.Set(&l1)
+	evaluations1.R1.Mul(&l1, &p.X)
+	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
+
+	// compute λ2 = -λ1-2y1/(x3-x1)
+	n.Double(&p.Y)
+	d.Sub(&x3, &p.X)
+	l2.Div(&n, &d)
+	l2.Add(&l2, &l1)
+	l2.Neg(&l2)
+
+	// compute x4 = λ2²-x1-x3
+	x4.Square(&l2)
+	x4.Sub(&x4, &p.X)
+	x4.Sub(&x4, &x3)
+
+	// compute y4 = λ2(x1 - x4)-y1
+	y4.Sub(&p.X, &x4)
+	y4.Mul(&l2, &y4)
+	y4.Sub(&y4, &p.Y)
+
+	// compute line2
+	evaluations2.R0.Set(&l2)
+	evaluations2.R1.Mul(&l2, &p.X)
+	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
+
+	p.X.Set(&x4)
+	p.Y.Set(&y4)
+}
diff --git a/ecc/bn254/pairing_test.go b/ecc/bn254/pairing_test.go
index 2744f8a387..19f3713224 100644
--- a/ecc/bn254/pairing_test.go
+++ b/ecc/bn254/pairing_test.go
@@ -189,6 +189,67 @@ func TestPairing(t *testing.T) {
 	properties.TestingRun(t, gopter.ConsoleReporter(false))
 }
 
+func TestDoubleAndAddStepEquivalence(t *testing.T) {
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genR1 := GenFr()
+	genR2 := GenFr()
+
+	properties.Property("[BN254] doubleAndAddStep: optimized (Eisenträger-Lauter-Montgomery) matches reference implementation", prop.ForAll(
+		func(s1, s2 fr.Element) bool {
+			var p1, p2, p1Ref G2Affine
+
+			// Get generator
+			_, _, _, g2Gen := Generators()
+
+			// Scale by random values
+			var s1Int, s2Int big.Int
+			s1.BigInt(&s1Int)
+			s2.BigInt(&s2Int)
+			p1.ScalarMultiplication(&g2Gen, &s1Int)
+			p2.ScalarMultiplication(&g2Gen, &s2Int)
+
+			// Fail if points are the same - this indicates a problem with random generation
+			if p1.X.Equal(&p2.X) {
+				return false
+			}
+
+			p1Ref.Set(&p1)
+
+			// Compute using optimized implementation
+			var eval1Opt, eval2Opt LineEvaluationAff
+			p1.doubleAndAddStep(&eval1Opt, &eval2Opt, &p2)
+
+			// Compute using reference implementation
+			var eval1Ref, eval2Ref LineEvaluationAff
+			doubleAndAddStepRef(&p1Ref, &eval1Ref, &eval2Ref, &p2)
+
+			// Compare results: resulting point coordinates
+			pointsEqual := p1.X.Equal(&p1Ref.X) && p1.Y.Equal(&p1Ref.Y)
+
+			// Compare results: line evaluation 1 (from first addition)
+			eval1Equal := eval1Opt.R0.Equal(&eval1Ref.R0) && eval1Opt.R1.Equal(&eval1Ref.R1)
+
+			// Compare results: line evaluation 2 (from doubling step)
+			eval2Equal := eval2Opt.R0.Equal(&eval2Ref.R0) && eval2Opt.R1.Equal(&eval2Ref.R1)
+
+			return pointsEqual && eval1Equal && eval2Equal
+		},
+		genR1,
+		genR2,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
 func TestMillerLoop(t *testing.T) {
 
 	t.Parallel()
@@ -404,6 +465,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/ecc/bw6-633/pairing_test.go b/ecc/bw6-633/pairing_test.go
index 36089c2cef..925fae1fe1 100644
--- a/ecc/bw6-633/pairing_test.go
+++ b/ecc/bw6-633/pairing_test.go
@@ -405,6 +405,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/ecc/bw6-761/pairing.go b/ecc/bw6-761/pairing.go
index 802cc9e374..bb025b8136 100644
--- a/ecc/bw6-761/pairing.go
+++ b/ecc/bw6-761/pairing.go
@@ -581,43 +581,73 @@ func (p *G2Affine) addStep(evaluations *LineEvaluationAff, a *G2Affine) {
 }
 
 func (p *G2Affine) doubleAndAddStep(evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
-	var n, d, l1, x3, l2, x4, y4 fp.Element
+	var A, B, A2, B2, X2A2, t, U, AU, invAU, invA, invU, l1, x3, l2, x4, y4 fp.Element
 
-	// compute λ1 = (y2-y1)/(x2-x1)
-	n.Sub(&p.Y, &a.Y)
-	d.Sub(&p.X, &a.X)
-	l1.Div(&n, &d)
+	// The Eisenträger-Lauter-Montgomery formula for 2P+Q (https://eprint.iacr.org/2003/257)
+	// computes both slopes λ1 and λ2 using a single field inversion via batch inversion.
+	//
+	// Given P = (x1, y1) and Q = (x2, y2), let:
+	//   A = x1 - x2
+	//   B = y1 - y2
+	//   U = B² - (2x1 + x2)·A²
+	//
+	// Then:
+	//   λ1 = B/A                    (slope for P + Q)
+	//   λ2 = -λ1 - 2y1·A²/U         (slope for P + (P+Q))
+	//
+	// We compute 1/A and 1/U using Montgomery's batch inversion:
+	//   1/A = U/(A·U) and 1/U = A/(A·U) with a single inversion of A·U.
+
+	// Compute A = x1 - x2 and B = y1 - y2
+	A.Sub(&p.X, &a.X)
+	B.Sub(&p.Y, &a.Y)
+
+	// Compute A² and B²
+	A2.Square(&A)
+	B2.Square(&B)
 
-	// compute x3 =λ1²-x1-x2
+	// Compute U = B² - (2x1 + x2)·A²
+	t.Double(&p.X).Add(&t, &a.X)
+	X2A2.Mul(&t, &A2)
+	U.Sub(&B2, &X2A2)
+
+	// Batch inversion: compute 1/A and 1/U with a single inversion
+	AU.Mul(&A, &U)
+	invAU.Inverse(&AU)
+	invA.Mul(&U, &invAU)
+	invU.Mul(&A, &invAU)
+
+	// λ1 = B/A = B·(1/A)
+	l1.Mul(&B, &invA)
+
+	// x3 = λ1² - x1 - x2
 	x3.Square(&l1)
 	x3.Sub(&x3, &p.X)
 	x3.Sub(&x3, &a.X)
 
-	// omit y3 computation
-
-	// compute line1
+	// line1 evaluation
 	evaluations1.R0.Set(&l1)
 	evaluations1.R1.Mul(&l1, &p.X)
 	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
 
-	// compute λ2 = -λ1-2y1/(x3-x1)
-	n.Double(&p.Y)
-	d.Sub(&x3, &p.X)
-	l2.Div(&n, &d)
+	// λ2 = -λ1 - 2y1·A²/U = -λ1 - 2y1·A²·(1/U)
+	l2.Double(&p.Y)
+	l2.Mul(&l2, &A2)
+	l2.Mul(&l2, &invU)
 	l2.Add(&l2, &l1)
 	l2.Neg(&l2)
 
-	// compute x4 = λ2²-x1-x3
+	// x4 = λ2² - x1 - x3
 	x4.Square(&l2)
 	x4.Sub(&x4, &p.X)
 	x4.Sub(&x4, &x3)
 
-	// compute y4 = λ2(x1 - x4)-y1
+	// y4 = λ2·(x1 - x4) - y1
 	y4.Sub(&p.X, &x4)
 	y4.Mul(&l2, &y4)
 	y4.Sub(&y4, &p.Y)
 
-	// compute line2
+	// line2 evaluation
 	evaluations2.R0.Set(&l2)
 	evaluations2.R1.Mul(&l2, &p.X)
 	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
diff --git a/ecc/bw6-761/pairing_compatibility_test.go b/ecc/bw6-761/pairing_compatibility_test.go
new file mode 100644
index 0000000000..d9cbffe9c8
--- /dev/null
+++ b/ecc/bw6-761/pairing_compatibility_test.go
@@ -0,0 +1,64 @@
+// Copyright 2020-2025 Consensys Software Inc.
+// Licensed under the Apache License, Version 2.0. See the LICENSE file for details.
+
+package bw6761
+
+import (
+	"github.com/consensys/gnark-crypto/ecc/bw6-761/fp"
+)
+
+// doubleAndAddStepRef is the reference (pre-optimization) implementation
+// of the doubleAndAddStep function. It computes 2P+Q using two field inversions.
+//
+// This version uses the standard chord-tangent method:
+//   - λ1 = (y2-y1)/(x2-x1) for P + Q
+//   - λ2 = -λ1 - 2y1/(x3-x1) for doubling and adding back
+//
+// The optimized version uses the Eisenträger-Lauter-Montgomery formula
+// (https://eprint.iacr.org/2003/257) which computes both slopes with a single
+// field inversion via Montgomery's batch inversion trick.
+func doubleAndAddStepRef(p *G2Affine, evaluations1, evaluations2 *LineEvaluationAff, a *G2Affine) {
+	var n, d, l1, x3, l2, x4, y4 fp.Element
+
+	// compute λ1 = (y2-y1)/(x2-x1)
+	n.Sub(&p.Y, &a.Y)
+	d.Sub(&p.X, &a.X)
+	l1.Div(&n, &d)
+
+	// compute x3 =λ1²-x1-x2
+	x3.Square(&l1)
+	x3.Sub(&x3, &p.X)
+	x3.Sub(&x3, &a.X)
+
+	// omit y3 computation
+
+	// compute line1
+	evaluations1.R0.Set(&l1)
+	evaluations1.R1.Mul(&l1, &p.X)
+	evaluations1.R1.Sub(&evaluations1.R1, &p.Y)
+
+	// compute λ2 = -λ1-2y1/(x3-x1)
+	n.Double(&p.Y)
+	d.Sub(&x3, &p.X)
+	l2.Div(&n, &d)
+	l2.Add(&l2, &l1)
+	l2.Neg(&l2)
+
+	// compute x4 = λ2²-x1-x3
+	x4.Square(&l2)
+	x4.Sub(&x4, &p.X)
+	x4.Sub(&x4, &x3)
+
+	// compute y4 = λ2(x1 - x4)-y1
+	y4.Sub(&p.X, &x4)
+	y4.Mul(&l2, &y4)
+	y4.Sub(&y4, &p.Y)
+
+	// compute line2
+	evaluations2.R0.Set(&l2)
+	evaluations2.R1.Mul(&l2, &p.X)
+	evaluations2.R1.Sub(&evaluations2.R1, &p.Y)
+
+	p.X.Set(&x4)
+	p.Y.Set(&y4)
+}
diff --git a/ecc/bw6-761/pairing_test.go b/ecc/bw6-761/pairing_test.go
index b60093f204..f010265263 100644
--- a/ecc/bw6-761/pairing_test.go
+++ b/ecc/bw6-761/pairing_test.go
@@ -218,6 +218,67 @@ func TestPairing(t *testing.T) {
 	properties.TestingRun(t, gopter.ConsoleReporter(false))
 }
 
+func TestDoubleAndAddStepEquivalence(t *testing.T) {
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genR1 := GenFr()
+	genR2 := GenFr()
+
+	properties.Property("[BW6-761] doubleAndAddStep: optimized (Eisenträger-Lauter-Montgomery) matches reference implementation", prop.ForAll(
+		func(s1, s2 fr.Element) bool {
+			var p1, p2, p1Ref G2Affine
+
+			// Get generator
+			_, _, _, g2Gen := Generators()
+
+			// Scale by random values
+			var s1Int, s2Int big.Int
+			s1.BigInt(&s1Int)
+			s2.BigInt(&s2Int)
+			p1.ScalarMultiplication(&g2Gen, &s1Int)
+			p2.ScalarMultiplication(&g2Gen, &s2Int)
+
+			// Fail if points are the same - this indicates a problem with random generation
+			if p1.X.Equal(&p2.X) {
+				return false
+			}
+
+			p1Ref.Set(&p1)
+
+			// Compute using optimized implementation
+			var eval1Opt, eval2Opt LineEvaluationAff
+			p1.doubleAndAddStep(&eval1Opt, &eval2Opt, &p2)
+
+			// Compute using reference implementation
+			var eval1Ref, eval2Ref LineEvaluationAff
+			doubleAndAddStepRef(&p1Ref, &eval1Ref, &eval2Ref, &p2)
+
+			// Compare results: resulting point coordinates
+			pointsEqual := p1.X.Equal(&p1Ref.X) && p1.Y.Equal(&p1Ref.Y)
+
+			// Compare results: line evaluation 1 (from first addition)
+			eval1Equal := eval1Opt.R0.Equal(&eval1Ref.R0) && eval1Opt.R1.Equal(&eval1Ref.R1)
+
+			// Compare results: line evaluation 2 (from doubling step)
+			eval2Equal := eval2Opt.R0.Equal(&eval2Ref.R0) && eval2Opt.R1.Equal(&eval2Ref.R1)
+
+			return pointsEqual && eval1Equal && eval2Equal
+		},
+		genR1,
+		genR2,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
 func TestMillerLoop(t *testing.T) {
 
 	t.Parallel()
@@ -433,6 +494,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine
diff --git a/internal/generator/pairing/template/tests/pairing.go.tmpl b/internal/generator/pairing/template/tests/pairing.go.tmpl
index 712386fb96..5f60b90887 100644
--- a/internal/generator/pairing/template/tests/pairing.go.tmpl
+++ b/internal/generator/pairing/template/tests/pairing.go.tmpl
@@ -230,6 +230,68 @@ func TestPairing(t *testing.T) {
 	properties.TestingRun(t, gopter.ConsoleReporter(false))
 }
 
+{{if or (eq .Name "bn254") (eq .Name "bls12-381") (eq .Name "bls12-377") (eq .Name "bw6-761")}}
+func TestDoubleAndAddStepEquivalence(t *testing.T) {
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genR1 := GenFr()
+	genR2 := GenFr()
+
+	properties.Property("[{{ toUpper .Name}}] doubleAndAddStep: optimized (Eisenträger-Lauter-Montgomery) matches reference implementation", prop.ForAll(
+		func(s1, s2 fr.Element) bool {
+			var p1, p2, p1Ref G2Affine
+
+			// Get generator
+			_, _, _, g2Gen := Generators()
+
+			// Scale by random values
+			var s1Int, s2Int big.Int
+			s1.BigInt(&s1Int)
+			s2.BigInt(&s2Int)
+			p1.ScalarMultiplication(&g2Gen, &s1Int)
+			p2.ScalarMultiplication(&g2Gen, &s2Int)
+
+			// Fail if points are the same - this indicates a problem with random generation
+			if p1.X.Equal(&p2.X) {
+				return false
+			}
+
+			p1Ref.Set(&p1)
+
+			// Compute using optimized implementation
+			var eval1Opt, eval2Opt LineEvaluationAff
+			p1.doubleAndAddStep(&eval1Opt, &eval2Opt, &p2)
+
+			// Compute using reference implementation
+			var eval1Ref, eval2Ref LineEvaluationAff
+			doubleAndAddStepRef(&p1Ref, &eval1Ref, &eval2Ref, &p2)
+
+			// Compare results: resulting point coordinates
+			pointsEqual := p1.X.Equal(&p1Ref.X) && p1.Y.Equal(&p1Ref.Y)
+
+			// Compare results: line evaluation 1 (from first addition)
+			eval1Equal := eval1Opt.R0.Equal(&eval1Ref.R0) && eval1Opt.R1.Equal(&eval1Ref.R1)
+
+			// Compare results: line evaluation 2 (from doubling step)
+			eval2Equal := eval2Opt.R0.Equal(&eval2Ref.R0) && eval2Opt.R1.Equal(&eval2Ref.R1)
+
+			return pointsEqual && eval1Equal && eval2Equal
+		},
+		genR1,
+		genR2,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+{{- end}}
 
 func TestMillerLoop(t *testing.T) {
 
@@ -464,6 +526,17 @@ func BenchmarkFinalExponentiation(b *testing.B) {
 
 }
 
+func BenchmarkPrecomputeLines(b *testing.B) {
+
+	var g2GenAff G2Affine
+	g2GenAff.FromJacobian(&g2Gen)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		PrecomputeLines(g2GenAff)
+	}
+}
+
 func BenchmarkMultiMiller(b *testing.B) {
 
 	var g1GenAff G1Affine