diff --git a/ecc/bw6-633/g1.go b/ecc/bw6-633/g1.go
index 282352cc91..b911740a77 100644
--- a/ecc/bw6-633/g1.go
+++ b/ecc/bw6-633/g1.go
@@ -995,51 +995,6 @@ func (p *g1Proj) FromAffine(Q *G1Affine) *g1Proj {
 	return p
 }
 
-// BatchProjectiveToAffineG1 converts points in Projective coordinates to Affine coordinates
-// performing a single field inversion (Montgomery batch inversion trick).
-func BatchProjectiveToAffineG1(points []g1Proj) []G1Affine {
-	result := make([]G1Affine, len(points))
-	zeroes := make([]bool, len(points))
-	accumulator := fp.One()
-
-	// batch invert all points[].Z coordinates with Montgomery batch inversion trick
-	// (stores points[].Z^-1 in result[i].X to avoid allocating a slice of fr.Elements)
-	for i := 0; i < len(points); i++ {
-		if points[i].z.IsZero() {
-			zeroes[i] = true
-			continue
-		}
-		result[i].X = accumulator
-		accumulator.Mul(&accumulator, &points[i].z)
-	}
-
-	var accInverse fp.Element
-	accInverse.Inverse(&accumulator)
-
-	for i := len(points) - 1; i >= 0; i-- {
-		if zeroes[i] {
-			// do nothing, (X=0, Y=0) is infinity point in affine
-			continue
-		}
-		result[i].X.Mul(&result[i].X, &accInverse)
-		accInverse.Mul(&accInverse, &points[i].z)
-	}
-
-	// batch convert to affine.
-	parallel.Execute(len(points), func(start, end int) {
-		for i := start; i < end; i++ {
-			if zeroes[i] {
-				// do nothing, (X=0, Y=0) is infinity point in affine
-				continue
-			}
-			a := result[i].X
-			result[i].X.Mul(&points[i].x, &a)
-			result[i].Y.Mul(&points[i].y, &a)
-		}
-	})
-	return result
-}
-
 // BatchJacobianToAffineG1 converts points in Jacobian coordinates to Affine coordinates
 // performing a single field inversion (Montgomery batch inversion trick).
 func BatchJacobianToAffineG1(points []G1Jac) []G1Affine {
diff --git a/ecc/bw6-756/g1.go b/ecc/bw6-756/g1.go
index 49437c7d8f..3858f8ce95 100644
--- a/ecc/bw6-756/g1.go
+++ b/ecc/bw6-756/g1.go
@@ -40,11 +40,6 @@ type g1JacExtended struct {
 	X, Y, ZZ, ZZZ fp.Element
 }
 
-// g1Proj point in projective coordinates
-type g1Proj struct {
-	x, y, z fp.Element
-}
-
 // -------------------------------------------------------------------------------------------------
 // Affine
 
@@ -965,81 +960,6 @@ func (p *g1JacExtended) doubleMixed(q *G1Affine) *g1JacExtended {
 	return p
 }
 
-// -------------------------------------------------------------------------------------------------
-// Homogenous projective
-
-// Set sets p to the provided point
-func (p *g1Proj) Set(a *g1Proj) *g1Proj {
-	p.x, p.y, p.z = a.x, a.y, a.z
-	return p
-}
-
-// Neg computes -G
-func (p *g1Proj) Neg(a *g1Proj) *g1Proj {
-	*p = *a
-	p.y.Neg(&a.y)
-	return p
-}
-
-// FromAffine sets p = Q, p in homogenous projective, Q in affine
-func (p *g1Proj) FromAffine(Q *G1Affine) *g1Proj {
-	if Q.X.IsZero() && Q.Y.IsZero() {
-		p.z.SetZero()
-		p.x.SetOne()
-		p.y.SetOne()
-		return p
-	}
-	p.z.SetOne()
-	p.x.Set(&Q.X)
-	p.y.Set(&Q.Y)
-	return p
-}
-
-// BatchProjectiveToAffineG1 converts points in Projective coordinates to Affine coordinates
-// performing a single field inversion (Montgomery batch inversion trick).
-func BatchProjectiveToAffineG1(points []g1Proj) []G1Affine {
-	result := make([]G1Affine, len(points))
-	zeroes := make([]bool, len(points))
-	accumulator := fp.One()
-
-	// batch invert all points[].Z coordinates with Montgomery batch inversion trick
-	// (stores points[].Z^-1 in result[i].X to avoid allocating a slice of fr.Elements)
-	for i := 0; i < len(points); i++ {
-		if points[i].z.IsZero() {
-			zeroes[i] = true
-			continue
-		}
-		result[i].X = accumulator
-		accumulator.Mul(&accumulator, &points[i].z)
-	}
-
-	var accInverse fp.Element
-	accInverse.Inverse(&accumulator)
-
-	for i := len(points) - 1; i >= 0; i-- {
-		if zeroes[i] {
-			// do nothing, (X=0, Y=0) is infinity point in affine
-			continue
-		}
-		result[i].X.Mul(&result[i].X, &accInverse)
-		accInverse.Mul(&accInverse, &points[i].z)
-	}
-
-	// batch convert to affine.
-	parallel.Execute(len(points), func(start, end int) {
-		for i := start; i < end; i++ {
-			if zeroes[i] {
-				// do nothing, (X=0, Y=0) is infinity point in affine
-				continue
-			}
-			a := result[i].X
-			result[i].X.Mul(&points[i].x, &a)
-			result[i].Y.Mul(&points[i].y, &a)
-		}
-	})
-	return result
-}
-
 // BatchJacobianToAffineG1 converts points in Jacobian coordinates to Affine coordinates
 // performing a single field inversion (Montgomery batch inversion trick).
 func BatchJacobianToAffineG1(points []G1Jac) []G1Affine {
diff --git a/ecc/bw6-756/g2.go b/ecc/bw6-756/g2.go
index 012542286c..20b0c193c6 100644
--- a/ecc/bw6-756/g2.go
+++ b/ecc/bw6-756/g2.go
@@ -40,6 +40,11 @@ type g2JacExtended struct {
 	X, Y, ZZ, ZZZ fp.Element
 }
 
+// g2Proj point in projective coordinates
+type g2Proj struct {
+	x, y, z fp.Element
+}
+
 // -------------------------------------------------------------------------------------------------
 // Affine
 
@@ -868,6 +873,36 @@ func (p *g2JacExtended) doubleMixed(q *G2Affine) *g2JacExtended {
 	return p
 }
 
+// -------------------------------------------------------------------------------------------------
+// Homogenous projective
+
+// Set sets p to the provided point
+func (p *g2Proj) Set(a *g2Proj) *g2Proj {
+	p.x, p.y, p.z = a.x, a.y, a.z
+	return p
+}
+
+// Neg computes -G
+func (p *g2Proj) Neg(a *g2Proj) *g2Proj {
+	*p = *a
+	p.y.Neg(&a.y)
+	return p
+}
+
+// FromAffine sets p = Q, p in homogenous projective, Q in affine
+func (p *g2Proj) FromAffine(Q *G2Affine) *g2Proj {
+	if Q.X.IsZero() && Q.Y.IsZero() {
+		p.z.SetZero()
+		p.x.SetOne()
+		p.y.SetOne()
+		return p
+	}
+	p.z.SetOne()
+	p.x.Set(&Q.X)
+	p.y.Set(&Q.Y)
+	return p
+}
+
 // BatchScalarMultiplicationG2 multiplies the same base by all scalars
 // and return resulting points in affine coordinates
 // uses a simple windowed-NAF like exponentiation algorithm
diff --git a/ecc/bw6-756/internal/fptower/e3.go b/ecc/bw6-756/internal/fptower/e3.go
index b2fa2859d0..15fc2d4bbe 100644
--- a/ecc/bw6-756/internal/fptower/e3.go
+++ b/ecc/bw6-756/internal/fptower/e3.go
@@ -140,6 +140,34 @@ func (z *E3) MulByElement(x *E3, y *fp.Element) *E3 {
 	return z
 }
 
+// MulBy12 multiplication by sparse element (0,b1,b2)
+func (x *E3) MulBy12(b1, b2 *fp.Element) *E3 {
+	var t1, t2, c0, tmp, c1, c2 fp.Element
+	t1.Mul(&x.A1, b1)
+	t2.Mul(&x.A2, b2)
+	c0.Add(&x.A1, &x.A2)
+	tmp.Add(b1, b2)
+	c0.Mul(&c0, &tmp)
+	c0.Sub(&c0, &t1)
+	c0.Sub(&c0, &t2)
+	c0.MulByNonResidue(&c0)
+	c1.Add(&x.A0, &x.A1)
+	c1.Mul(&c1, b1)
+	c1.Sub(&c1, &t1)
+	tmp.MulByNonResidue(&t2)
+	c1.Add(&c1, &tmp)
+	tmp.Add(&x.A0, &x.A2)
+	c2.Mul(b2, &tmp)
+	c2.Sub(&c2, &t2)
+	c2.Add(&c2, &t1)
+
+	x.A0 = c0
+	x.A1 = c1
+	x.A2 = c2
+
+	return x
+}
+
 // MulBy01 multiplication by sparse element (c0,c1,0)
 func (z *E3) MulBy01(c0, c1 *fp.Element) *E3 {
 
diff --git a/ecc/bw6-756/internal/fptower/e6_pairing.go b/ecc/bw6-756/internal/fptower/e6_pairing.go
index f4af22821d..fe7eee34c5 100644
--- a/ecc/bw6-756/internal/fptower/e6_pairing.go
+++ b/ecc/bw6-756/internal/fptower/e6_pairing.go
@@ -126,65 +126,68 @@ func (z *E6) ExptMinus1Div3(x *E6) *E6 {
 	return z
 }
 
-// MulBy034 multiplication by sparse element (c0,0,0,c3,c4,0)
-func (z *E6) MulBy034(c0, c3, c4 *fp.Element) *E6 {
+// MulBy014 multiplication by sparse element (c0,c1,0,0,c4,0)
+func (z *E6) MulBy014(c0, c1, c4 *fp.Element) *E6 {
 
-	var a, b, d E3
+	var a, b E3
+	var d fp.Element
 
-	a.MulByElement(&z.B0, c0)
+	a.Set(&z.B0)
+	a.MulBy01(c0, c1)
 
 	b.Set(&z.B1)
-	b.MulBy01(c3, c4)
+	b.MulBy1(c4)
+	d.Add(c1, c4)
 
-	c0.Add(c0, c3)
-	d.Add(&z.B0, &z.B1)
-	d.MulBy01(c0, c4)
-
-	z.B1.Add(&a, &b).Neg(&z.B1).Add(&z.B1, &d)
-	z.B0.MulByNonResidue(&b).Add(&z.B0, &a)
+	z.B1.Add(&z.B1, &z.B0)
+	z.B1.MulBy01(c0, &d)
+	z.B1.Sub(&z.B1, &a)
+	z.B1.Sub(&z.B1, &b)
+	z.B0.MulByNonResidue(&b)
+	z.B0.Add(&z.B0, &a)
 
 	return z
 }
 
-// Mul034By034 multiplication of sparse element (c0,0,0,c3,c4,0) by sparse element (d0,0,0,d3,d4,0)
-func Mul034By034(d0, d3, d4, c0, c3, c4 *fp.Element) [5]fp.Element {
-	var z00, tmp, x0, x3, x4, x04, x03, x34 fp.Element
+// Mul014By014 multiplication of sparse element (c0,c1,0,0,c4,0) by sparse element (d0,d1,0,0,d4,0)
+func Mul014By014(d0, d1, d4, c0, c1, c4 *fp.Element) [5]fp.Element {
+	var z00, tmp, x0, x1, x4, x04, x01, x14 fp.Element
 	x0.Mul(c0, d0)
-	x3.Mul(c3, d3)
+	x1.Mul(c1, d1)
 	x4.Mul(c4, d4)
 	tmp.Add(c0, c4)
 	x04.Add(d0, d4).
 		Mul(&x04, &tmp).
 		Sub(&x04, &x0).
 		Sub(&x04, &x4)
-	tmp.Add(c0, c3)
-	x03.Add(d0, d3).
-		Mul(&x03, &tmp).
-		Sub(&x03, &x0).
-		Sub(&x03, &x3)
-	tmp.Add(c3, c4)
-	x34.Add(d3, d4).
-		Mul(&x34, &tmp).
-		Sub(&x34, &x3).
-		Sub(&x34, &x4)
+	tmp.Add(c0, c1)
+	x01.Add(d0, d1).
+		Mul(&x01, &tmp).
+		Sub(&x01, &x0).
+		Sub(&x01, &x1)
+	tmp.Add(c1, c4)
+	x14.Add(d1, d4).
+		Mul(&x14, &tmp).
+		Sub(&x14, &x1).
+		Sub(&x14, &x4)
 
 	z00.MulByNonResidue(&x4).
 		Add(&z00, &x0)
 
-	return [5]fp.Element{z00, x3, x34, x03, x04}
+	return [5]fp.Element{z00, x01, x1, x04, x14}
 }
 
-// MulBy01234 multiplies z by an E12 sparse element of the form (x0, x1, x2, x3, x4, 0)
-func (z *E6) MulBy01234(x *[5]fp.Element) *E6 {
+// MulBy01245 multiplies z by an E12 sparse element of the form (x0, x1, x2, 0, x4, x5)
+func (z *E6) MulBy01245(x *[5]fp.Element) *E6 {
 	var c1, a, b, c, z0, z1 E3
 	c0 := &E3{A0: x[0], A1: x[1], A2: x[2]}
-	c1.A0 = x[3]
-	c1.A1 = x[4]
+	c1.A1 = x[3]
+	c1.A2 = x[4]
 	a.Add(&z.B0, &z.B1)
 	b.Add(c0, &c1)
 	a.Mul(&a, &b)
 	b.Mul(&z.B0, c0)
-	c.Set(&z.B1).MulBy01(&x[3], &x[4])
+	c.Set(&z.B1).MulBy12(&x[3], &x[4])
 	z1.Sub(&a, &b)
 	z1.Sub(&z1, &c)
 	z0.MulByNonResidue(&c)
diff --git a/ecc/bw6-756/pairing.go b/ecc/bw6-756/pairing.go
index bd77856426..d4a65e0f0d 100644
--- a/ecc/bw6-756/pairing.go
+++ b/ecc/bw6-756/pairing.go
@@ -124,13 +124,12 @@ func FinalExponentiation(z *GT, _z ...*GT) GT {
 	return result
 }
 
-// MillerLoop Optimal Tate alternative (or twisted ate or Eta revisited)
-// computes the multi-Miller loop
+// MillerLoop computes the multi-Miller loop
 // ∏ᵢ MillerLoop(Pᵢ, Qᵢ) =
-// ∏ᵢ { fᵢ_{x₀+1+λ(x₀³-x₀²-x₀),Pᵢ}(Qᵢ) }
+// ∏ᵢ { fᵢ_{x₀+1+λ(x₀³-x₀²-x₀),Qᵢ}(Pᵢ) }
 //
 // Alg.2 in https://eprint.iacr.org/2021/1359.pdf
-// Eq. (6) in https://hackmd.io/@gnark/BW6-761-changes
+// Eq. (6') in https://hackmd.io/@gnark/BW6-761-changes
 func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	// check input size match
 	n := len(P)
@@ -139,29 +138,34 @@ func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	}
 
 	// filter infinity points
-	p0 := make([]G1Affine, 0, n)
-	q := make([]G2Affine, 0, n)
+	p := make([]G1Affine, 0, n)
+	q0 := make([]G2Affine, 0, n)
 
 	for k := 0; k < n; k++ {
 		if P[k].IsInfinity() || Q[k].IsInfinity() {
 			continue
 		}
-		p0 = append(p0, P[k])
-		q = append(q, Q[k])
+		p = append(p, P[k])
+		q0 = append(q0, Q[k])
 	}
 
-	n = len(q)
+	n = len(p)
 
 	// precomputations
-	pProj1 := make([]g1Proj, n)
-	p1 := make([]G1Affine, n)
+	qProj1 := make([]g2Proj, n)
+	q1 := make([]G2Affine, n)
+	q1Neg := make([]G2Affine, n)
+	q0Neg := make([]G2Affine, n)
 	for k := 0; k < n; k++ {
-		p1[k].Y.Neg(&p0[k].Y)
-		p1[k].X.Mul(&p0[k].X, &thirdRootOneG2)
-		pProj1[k].FromAffine(&p1[k])
+		q1[k].Y.Neg(&q0[k].Y)
+		q0Neg[k].X.Set(&q0[k].X)
+		q0Neg[k].Y.Set(&q1[k].Y)
+		q1[k].X.Mul(&q0[k].X, &thirdRootOneG1)
+		qProj1[k].FromAffine(&q1[k])
+		q1Neg[k].Neg(&q1[k])
 	}
 
-	// f_{a0+λ*a1,P}(Q)
+	// f_{a0+λ*a1,Q}(P)
 	var result GT
 	result.SetOne()
 	var l, l0 lineEvaluation
@@ -170,113 +174,83 @@ func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	var j int8
 
 	if n >= 1 {
-		// i = len(loopCounter0) - 2, separately to avoid an E12 Square
+		// i = 189, separately to avoid an E12 Square
 		// (Square(res) = 1² = 1)
 		// j = 0
-		// k = 0, separately to avoid MulBy034 (res × ℓ)
+		// k = 0, separately to avoid MulBy014 (res × ℓ)
 		// (assign line to res)
-
-		// pProj1[0] ← 2pProj1[0] and l0 the tangent ℓ passing 2pProj1[0]
-		pProj1[0].doubleStep(&l0)
+		// qProj1[0] ← 2qProj1[0] and l0 the tangent ℓ passing 2qProj1[0]
+		qProj1[0].doubleStep(&l0)
 		// line evaluation at Q[0] (assign)
-		result.B1.A0.Mul(&l0.r1, &q[0].X)
-		result.B0.A0.Mul(&l0.r0, &q[0].Y)
-		result.B1.A1.Set(&l0.r2)
+		result.B0.A0.Set(&l0.r0)
+		result.B0.A1.Mul(&l0.r1, &p[0].X)
+		result.B1.A1.Mul(&l0.r2, &p[0].Y)
 	}
 
 	// k = 1
 	if n >= 2 {
-		// pProj1[1] ← 2pProj1[1] and l0 the tangent ℓ passing 2pProj1[1]
-		pProj1[1].doubleStep(&l0)
-		// line evaluation at Q[0]
-		l0.r1.Mul(&l0.r1, &q[1].X)
-		l0.r0.Mul(&l0.r0, &q[1].Y)
-		// ℓ × res
-		prodLines = fptower.Mul034By034(&l0.r0, &l0.r1, &l0.r2, &result.B0.A0, &result.B1.A0, &result.B1.A1)
+		// qProj1[1] ← 2qProj1[1] and l0 the tangent ℓ passing 2qProj1[1]
+		qProj1[1].doubleStep(&l0)
+		// line evaluation at Q[1]
+		l0.r1.Mul(&l0.r1, &p[1].X)
+		l0.r2.Mul(&l0.r2, &p[1].Y)
+		prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &result.B0.A0, &result.B0.A1, &result.B1.A1)
 		result.B0.A0 = prodLines[0]
 		result.B0.A1 = prodLines[1]
 		result.B0.A2 = prodLines[2]
-		result.B1.A0 = prodLines[3]
-		result.B1.A1 = prodLines[4]
+		result.B1.A1 = prodLines[3]
+		result.B1.A2 = prodLines[4]
 	}
 
 	// k >= 2
 	for k := 2; k < n; k++ {
-		// pProj1[1] ← 2pProj1[1] and l0 the tangent ℓ passing 2pProj1[1]
-		pProj1[k].doubleStep(&l0)
+		// qProj1[k] ← 2qProj1[k] and l0 the tangent ℓ passing 2qProj1[k]
+		qProj1[k].doubleStep(&l0)
 		// line evaluation at Q[k]
-		l0.r1.Mul(&l0.r1, &q[k].X)
-		l0.r0.Mul(&l0.r0, &q[k].Y)
+		l0.r1.Mul(&l0.r1, &p[k].X)
+		l0.r2.Mul(&l0.r2, &p[k].Y)
 		// ℓ × res
-		result.MulBy034(&l0.r0, &l0.r1, &l0.r2)
+		result.MulBy014(&l0.r0, &l0.r1, &l0.r2)
 	}
 
-	var tmp G1Affine
-	for i := len(loopCounter0) - 3; i >= 1; i-- {
-		// (∏ᵢfᵢ)²
-		// mutualize the square among n Miller loops
+	for i := 188; i >= 1; i-- {
 		result.Square(&result)
 
 		j = loopCounter1[i]*3 + loopCounter0[i]
 
 		for k := 0; k < n; k++ {
-			// pProj1[1] ← 2pProj1[1] and l0 the tangent ℓ passing 2pProj1[1]
-			pProj1[k].doubleStep(&l0)
-			// line evaluation at Q[k]
-			l0.r1.Mul(&l0.r1, &q[k].X)
-			l0.r0.Mul(&l0.r0, &q[k].Y)
+			qProj1[k].doubleStep(&l0)
+			l0.r1.Mul(&l0.r1, &p[k].X)
+			l0.r2.Mul(&l0.r2, &p[k].Y)
 
 			switch j {
 			// cases -4, -2, 2, 4 do not occur, given the static loopCounters
 			case -3:
-				tmp.Neg(&p1[k])
-				// pProj1[k] ← pProj1[k]-p1[k] and
-				// l the line ℓ passing pProj1[k] and -p1[k]
-				pProj1[k].addMixedStep(&l, &tmp)
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q1Neg[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			case -1:
-				tmp.Neg(&p0[k])
-				// pProj1[k] ← pProj1[k]-p0[k] and
-				// l the line ℓ passing pProj1[k] and -p0[k]
-				pProj1[k].addMixedStep(&l, &tmp)
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q0Neg[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			case 0:
-				// ℓ × res
-				result.MulBy034(&l0.r0, &l0.r1, &l0.r2)
+				result.MulBy014(&l0.r0, &l0.r1, &l0.r2)
 			case 1:
-				// pProj1[k] ← pProj1[k]+p0[k] and
-				// l the line ℓ passing pProj1[k] and p0[k]
-				pProj1[k].addMixedStep(&l, &p0[k])
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q0[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			case 3:
-				// pProj1[k] ← pProj1[k]+p1[k] and
-				// l the line ℓ passing pProj1[k] and p1[k]
-				pProj1[k].addMixedStep(&l, &p1[k])
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q1[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			default:
 				return GT{}, errors.New("invalid loopCounter")
 			}
@@ -287,30 +261,23 @@ func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	// j = -3
 	result.Square(&result)
 	for k := 0; k < n; k++ {
-		// pProj1[k] ← 2pProj1[k] and l0 the tangent ℓ passing 2pProj1[k]
-		pProj1[k].doubleStep(&l0)
-		// line evaluation at Q[k]
-		l0.r1.Mul(&l0.r1, &q[k].X)
-		l0.r0.Mul(&l0.r0, &q[k].Y)
-
-		tmp.Neg(&p1[k])
-		// l the line passing pProj1[k] and -p1
-		pProj1[k].lineCompute(&l, &tmp)
-		// line evaluation at Q[k]
-		l.r1.Mul(&l.r1, &q[k].X)
-		l.r0.Mul(&l.r0, &q[k].Y)
-		// ℓ × ℓ
-		prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-		// (ℓ × ℓ) × res
-		result.MulBy01234(&prodLines)
+		qProj1[k].doubleStep(&l0)
+		l0.r1.Mul(&l0.r1, &p[k].X)
+		l0.r2.Mul(&l0.r2, &p[k].Y)
+		qProj1[k].lineCompute(&l, &q1Neg[k])
+		l.r1.Mul(&l.r1, &p[k].X)
+		l.r2.Mul(&l.r2, &p[k].Y)
+		prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+		result.MulBy01245(&prodLines)
 	}
 
 	return result, nil
+
 }
 
 // doubleStep doubles a point in Homogenous projective coordinates, and evaluates the line in Miller loop
 // https://eprint.iacr.org/2013/722.pdf (Section 4.3)
-func (p *g1Proj) doubleStep(evaluations *lineEvaluation) {
+func (p *g2Proj) doubleStep(evaluations *lineEvaluation) {
 
 	// get some Element from our pool
 	var t1, A, B, C, D, E, EE, F, G, H, I, J, K fp.Element
@@ -320,8 +287,7 @@ func (p *g1Proj) doubleStep(evaluations *lineEvaluation) {
 	C.Square(&p.z)
 	D.Double(&C).
 		Add(&D, &C)
-	// E.Mul(&D, &bCurveCoeff)
-	E.Set(&D)
+	E.Mul(&D, &bTwistCurveCoeff)
 	F.Double(&E).
 		Add(&F, &E)
 	G.Add(&B, &F)
@@ -344,15 +310,15 @@ func (p *g1Proj) doubleStep(evaluations *lineEvaluation) {
 	p.z.Mul(&B, &H)
 
 	// Line evaluation
-	evaluations.r0.Neg(&H)
+	evaluations.r0.Set(&I)
 	evaluations.r1.Double(&J).
 		Add(&evaluations.r1, &J)
-	evaluations.r2.Set(&I)
+	evaluations.r2.Neg(&H)
 }
 
 // addMixedStep point addition in Mixed Homogenous projective and Affine coordinates
 // https://eprint.iacr.org/2013/722.pdf (Section 4.3)
-func (p *g1Proj) addMixedStep(evaluations *lineEvaluation, a *G1Affine) {
+func (p *g2Proj) addMixedStep(evaluations *lineEvaluation, a *G2Affine) {
 
 	// get some Element from our pool
 	var Y2Z1, X2Z1, O, L, C, D, E, F, G, H, t0, t1, t2, J fp.Element
@@ -382,14 +348,14 @@ func (p *g1Proj) addMixedStep(evaluations *lineEvaluation, a *G1Affine) {
 		Sub(&J, &t2)
 
 	// Line evaluation
-	evaluations.r0.Set(&L)
+	evaluations.r0.Set(&J)
 	evaluations.r1.Neg(&O)
-	evaluations.r2.Set(&J)
+	evaluations.r2.Set(&L)
 }
 
 // lineCompute computes the line through p in Homogenous projective coordinates
 // and a in affine coordinates. It does not compute the resulting point p+a.
-func (p *g1Proj) lineCompute(evaluations *lineEvaluation, a *G1Affine) {
+func (p *g2Proj) lineCompute(evaluations *lineEvaluation, a *G2Affine) {
 
 	// get some Element from our pool
 	var Y2Z1, X2Z1, O, L, t2, J fp.Element
@@ -402,7 +368,7 @@ func (p *g1Proj) lineCompute(evaluations *lineEvaluation, a *G1Affine) {
 		Sub(&J, &t2)
 
 	// Line evaluation
-	evaluations.r0.Set(&L)
+	evaluations.r0.Set(&J)
 	evaluations.r1.Neg(&O)
-	evaluations.r2.Set(&J)
+	evaluations.r2.Set(&L)
 }
diff --git a/ecc/bw6-761/g1.go b/ecc/bw6-761/g1.go
index 70a51e6989..aaded0f483 100644
--- a/ecc/bw6-761/g1.go
+++ b/ecc/bw6-761/g1.go
@@ -40,11 +40,6 @@ type g1JacExtended struct {
 	X, Y, ZZ, ZZZ fp.Element
 }
 
-// g1Proj point in projective coordinates
-type g1Proj struct {
-	x, y, z fp.Element
-}
-
 // -------------------------------------------------------------------------------------------------
 // Affine
 
@@ -976,81 +971,6 @@ func (p *g1JacExtended) doubleMixed(q *G1Affine) *g1JacExtended {
 	return p
 }
 
-// -------------------------------------------------------------------------------------------------
-// Homogenous projective
-
-// Set sets p to the provided point
-func (p *g1Proj) Set(a *g1Proj) *g1Proj {
-	p.x, p.y, p.z = a.x, a.y, a.z
-	return p
-}
-
-// Neg computes -G
-func (p *g1Proj) Neg(a *g1Proj) *g1Proj {
-	*p = *a
-	p.y.Neg(&a.y)
-	return p
-}
-
-// FromAffine sets p = Q, p in homogenous projective, Q in affine
-func (p *g1Proj) FromAffine(Q *G1Affine) *g1Proj {
-	if Q.X.IsZero() && Q.Y.IsZero() {
-		p.z.SetZero()
-		p.x.SetOne()
-		p.y.SetOne()
-		return p
-	}
-	p.z.SetOne()
-	p.x.Set(&Q.X)
-	p.y.Set(&Q.Y)
-	return p
-}
-
-// BatchProjectiveToAffineG1 converts points in Projective coordinates to Affine coordinates
-// performing a single field inversion (Montgomery batch inversion trick).
-func BatchProjectiveToAffineG1(points []g1Proj) []G1Affine {
-	result := make([]G1Affine, len(points))
-	zeroes := make([]bool, len(points))
-	accumulator := fp.One()
-
-	// batch invert all points[].Z coordinates with Montgomery batch inversion trick
-	// (stores points[].Z^-1 in result[i].X to avoid allocating a slice of fr.Elements)
-	for i := 0; i < len(points); i++ {
-		if points[i].z.IsZero() {
-			zeroes[i] = true
-			continue
-		}
-		result[i].X = accumulator
-		accumulator.Mul(&accumulator, &points[i].z)
-	}
-
-	var accInverse fp.Element
-	accInverse.Inverse(&accumulator)
-
-	for i := len(points) - 1; i >= 0; i-- {
-		if zeroes[i] {
-			// do nothing, (X=0, Y=0) is infinity point in affine
-			continue
-		}
-		result[i].X.Mul(&result[i].X, &accInverse)
-		accInverse.Mul(&accInverse, &points[i].z)
-	}
-
-	// batch convert to affine.
-	parallel.Execute(len(points), func(start, end int) {
-		for i := start; i < end; i++ {
-			if zeroes[i] {
-				// do nothing, (X=0, Y=0) is infinity point in affine
-				continue
-			}
-			a := result[i].X
-			result[i].X.Mul(&points[i].x, &a)
-			result[i].Y.Mul(&points[i].y, &a)
-		}
-	})
-	return result
-}
-
 // BatchJacobianToAffineG1 converts points in Jacobian coordinates to Affine coordinates
 // performing a single field inversion (Montgomery batch inversion trick).
 func BatchJacobianToAffineG1(points []G1Jac) []G1Affine {
diff --git a/ecc/bw6-761/g2.go b/ecc/bw6-761/g2.go
index 590621d9f8..e8bd0d4a31 100644
--- a/ecc/bw6-761/g2.go
+++ b/ecc/bw6-761/g2.go
@@ -40,6 +40,11 @@ type g2JacExtended struct {
 	X, Y, ZZ, ZZZ fp.Element
 }
 
+// g2Proj point in projective coordinates
+type g2Proj struct {
+	x, y, z fp.Element
+}
+
 // -------------------------------------------------------------------------------------------------
 // Affine
 
@@ -882,6 +887,36 @@ func (p *g2JacExtended) doubleMixed(q *G2Affine) *g2JacExtended {
 	return p
 }
 
+// -------------------------------------------------------------------------------------------------
+// Homogenous projective
+
+// Set sets p to the provided point
+func (p *g2Proj) Set(a *g2Proj) *g2Proj {
+	p.x, p.y, p.z = a.x, a.y, a.z
+	return p
+}
+
+// Neg computes -G
+func (p *g2Proj) Neg(a *g2Proj) *g2Proj {
+	*p = *a
+	p.y.Neg(&a.y)
+	return p
+}
+
+// FromAffine sets p = Q, p in homogenous projective, Q in affine
+func (p *g2Proj) FromAffine(Q *G2Affine) *g2Proj {
+	if Q.X.IsZero() && Q.Y.IsZero() {
+		p.z.SetZero()
+		p.x.SetOne()
+		p.y.SetOne()
+		return p
+	}
+	p.z.SetOne()
+	p.x.Set(&Q.X)
+	p.y.Set(&Q.Y)
+	return p
+}
+
 // BatchScalarMultiplicationG2 multiplies the same base by all scalars
 // and return resulting points in affine coordinates
 // uses a simple windowed-NAF like exponentiation algorithm
diff --git a/ecc/bw6-761/internal/fptower/e3.go b/ecc/bw6-761/internal/fptower/e3.go
index 7434ebd59c..29e7fd7ff7 100644
--- a/ecc/bw6-761/internal/fptower/e3.go
+++ b/ecc/bw6-761/internal/fptower/e3.go
@@ -140,6 +140,34 @@ func (z *E3) MulByElement(x *E3, y *fp.Element) *E3 {
 	return z
 }
 
+// MulBy12 multiplication by sparse element (0,b1,b2)
+func (x *E3) MulBy12(b1, b2 *fp.Element) *E3 {
+	var t1, t2, c0, tmp, c1, c2 fp.Element
+	t1.Mul(&x.A1, b1)
+	t2.Mul(&x.A2, b2)
+	c0.Add(&x.A1, &x.A2)
+	tmp.Add(b1, b2)
+	c0.Mul(&c0, &tmp)
+	c0.Sub(&c0, &t1)
+	c0.Sub(&c0, &t2)
+	c0.MulByNonResidue(&c0)
+	c1.Add(&x.A0, &x.A1)
+	c1.Mul(&c1, b1)
+	c1.Sub(&c1, &t1)
+	tmp.MulByNonResidue(&t2)
+	c1.Add(&c1, &tmp)
+	tmp.Add(&x.A0, &x.A2)
+	c2.Mul(b2, &tmp)
+	c2.Sub(&c2, &t2)
+	c2.Add(&c2, &t1)
+
+	x.A0 = c0
+	x.A1 = c1
+	x.A2 = c2
+
+	return x
+}
+
 // MulBy01 multiplication by sparse element (c0,c1,0)
 func (z *E3) MulBy01(c0, c1 *fp.Element) *E3 {
 
diff --git a/ecc/bw6-761/internal/fptower/e6_pairing.go b/ecc/bw6-761/internal/fptower/e6_pairing.go
index 21178257e0..6f9de4819d 100644
--- a/ecc/bw6-761/internal/fptower/e6_pairing.go
+++ b/ecc/bw6-761/internal/fptower/e6_pairing.go
@@ -149,65 +149,68 @@ func (z *E6) Expc2(x *E6) *E6 {
 	return z
 }
 
-// MulBy034 multiplication by sparse element (c0,0,0,c3,c4,0)
-func (z *E6) MulBy034(c0, c3, c4 *fp.Element) *E6 {
+// MulBy014 multiplication by sparse element (c0,c1,0,0,c4,0)
+func (z *E6) MulBy014(c0, c1, c4 *fp.Element) *E6 {
 
-	var a, b, d E3
+	var a, b E3
+	var d fp.Element
 
-	a.MulByElement(&z.B0, c0)
+	a.Set(&z.B0)
+	a.MulBy01(c0, c1)
 
 	b.Set(&z.B1)
-	b.MulBy01(c3, c4)
+	b.MulBy1(c4)
+	d.Add(c1, c4)
 
-	c0.Add(c0, c3)
-	d.Add(&z.B0, &z.B1)
-	d.MulBy01(c0, c4)
-
-	z.B1.Add(&a, &b).Neg(&z.B1).Add(&z.B1, &d)
-	z.B0.MulByNonResidue(&b).Add(&z.B0, &a)
+	z.B1.Add(&z.B1, &z.B0)
+	z.B1.MulBy01(c0, &d)
+	z.B1.Sub(&z.B1, &a)
+	z.B1.Sub(&z.B1, &b)
+	z.B0.MulByNonResidue(&b)
+	z.B0.Add(&z.B0, &a)
 
 	return z
 }
 
-// Mul034By034 multiplication of sparse element (c0,0,0,c3,c4,0) by sparse element (d0,0,0,d3,d4,0)
-func Mul034By034(d0, d3, d4, c0, c3, c4 *fp.Element) [5]fp.Element {
-	var z00, tmp, x0, x3, x4, x04, x03, x34 fp.Element
+// Mul014By014 multiplication of sparse element (c0,c1,0,0,c4,0) by sparse element (d0,d1,0,0,d4,0)
+func Mul014By014(d0, d1, d4, c0, c1, c4 *fp.Element) [5]fp.Element {
+	var z00, tmp, x0, x1, x4, x04, x01, x14 fp.Element
 	x0.Mul(c0, d0)
-	x3.Mul(c3, d3)
+	x1.Mul(c1, d1)
 	x4.Mul(c4, d4)
 	tmp.Add(c0, c4)
 	x04.Add(d0, d4).
 		Mul(&x04, &tmp).
 		Sub(&x04, &x0).
 		Sub(&x04, &x4)
-	tmp.Add(c0, c3)
-	x03.Add(d0, d3).
-		Mul(&x03, &tmp).
-		Sub(&x03, &x0).
-		Sub(&x03, &x3)
-	tmp.Add(c3, c4)
-	x34.Add(d3, d4).
-		Mul(&x34, &tmp).
-		Sub(&x34, &x3).
-		Sub(&x34, &x4)
+	tmp.Add(c0, c1)
+	x01.Add(d0, d1).
+		Mul(&x01, &tmp).
+		Sub(&x01, &x0).
+		Sub(&x01, &x1)
+	tmp.Add(c1, c4)
+	x14.Add(d1, d4).
+		Mul(&x14, &tmp).
+		Sub(&x14, &x1).
+		Sub(&x14, &x4)
 
 	z00.MulByNonResidue(&x4).
 		Add(&z00, &x0)
 
-	return [5]fp.Element{z00, x3, x34, x03, x04}
+	return [5]fp.Element{z00, x01, x1, x04, x14}
 }
 
-// MulBy01234 multiplies z by an E12 sparse element of the form (x0, x1, x2, x3, x4, 0)
-func (z *E6) MulBy01234(x *[5]fp.Element) *E6 {
+// MulBy01245 multiplies z by an E12 sparse element of the form (x0, x1, x2, 0, x4, x5)
+func (z *E6) MulBy01245(x *[5]fp.Element) *E6 {
 	var c1, a, b, c, z0, z1 E3
 	c0 := &E3{A0: x[0], A1: x[1], A2: x[2]}
-	c1.A0 = x[3]
-	c1.A1 = x[4]
+	c1.A1 = x[3]
+	c1.A2 = x[4]
 	a.Add(&z.B0, &z.B1)
 	b.Add(c0, &c1)
 	a.Mul(&a, &b)
 	b.Mul(&z.B0, c0)
-	c.Set(&z.B1).MulBy01(&x[3], &x[4])
+	c.Set(&z.B1).MulBy12(&x[3], &x[4])
 	z1.Sub(&a, &b)
 	z1.Sub(&z1, &c)
 	z0.MulByNonResidue(&c)
diff --git a/ecc/bw6-761/pairing.go b/ecc/bw6-761/pairing.go
index 12ca16ba03..d1160bbe76 100644
--- a/ecc/bw6-761/pairing.go
+++ b/ecc/bw6-761/pairing.go
@@ -121,13 +121,12 @@ func FinalExponentiation(z *GT, _z ...*GT) GT {
 	return result
 }
 
-// MillerLoop Optimal Tate alternative (or twisted ate or Eta revisited)
-// computes the multi-Miller loop
+// MillerLoop computes the multi-Miller loop
 // ∏ᵢ MillerLoop(Pᵢ, Qᵢ) =
-// ∏ᵢ { fᵢ_{x₀+1+λ(x₀³-x₀²-x₀),Pᵢ}(Qᵢ) }
+// ∏ᵢ { fᵢ_{x₀+1+λ(x₀³-x₀²-x₀),Qᵢ}(Pᵢ) }
 //
 // Alg.2 in https://eprint.iacr.org/2021/1359.pdf
-// Eq. (6) in https://hackmd.io/@gnark/BW6-761-changes
+// Eq. (6') in https://hackmd.io/@gnark/BW6-761-changes
 func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	// check input size match
 	n := len(P)
@@ -136,29 +135,34 @@ func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	}
 
 	// filter infinity points
-	p0 := make([]G1Affine, 0, n)
-	q := make([]G2Affine, 0, n)
+	p := make([]G1Affine, 0, n)
+	q0 := make([]G2Affine, 0, n)
 
 	for k := 0; k < n; k++ {
 		if P[k].IsInfinity() || Q[k].IsInfinity() {
 			continue
 		}
-		p0 = append(p0, P[k])
-		q = append(q, Q[k])
+		p = append(p, P[k])
+		q0 = append(q0, Q[k])
 	}
 
-	n = len(q)
+	n = len(p)
 
 	// precomputations
-	pProj1 := make([]g1Proj, n)
-	p1 := make([]G1Affine, n)
+	qProj1 := make([]g2Proj, n)
+	q1 := make([]G2Affine, n)
+	q1Neg := make([]G2Affine, n)
+	q0Neg := make([]G2Affine, n)
 	for k := 0; k < n; k++ {
-		p1[k].Y.Neg(&p0[k].Y)
-		p1[k].X.Mul(&p0[k].X, &thirdRootOneG2)
-		pProj1[k].FromAffine(&p1[k])
+		q1[k].Y.Neg(&q0[k].Y)
+		q0Neg[k].X.Set(&q0[k].X)
+		q0Neg[k].Y.Set(&q1[k].Y)
+		q1[k].X.Mul(&q0[k].X, &thirdRootOneG1)
+		qProj1[k].FromAffine(&q1[k])
+		q1Neg[k].Neg(&q1[k])
 	}
 
-	// f_{a0+λ*a1,P}(Q)
+	// f_{a0+λ*a1,Q}(P)
 	var result GT
 	result.SetOne()
 	var l, l0 lineEvaluation
@@ -167,113 +171,83 @@ func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	var j int8
 
 	if n >= 1 {
-		// i = len(loopCounter0) - 2, separately to avoid an E12 Square
+		// i = 188, separately to avoid an E12 Square
 		// (Square(res) = 1² = 1)
 		// j = 0
-		// k = 0, separately to avoid MulBy034 (res × ℓ)
+		// k = 0, separately to avoid MulBy014 (res × ℓ)
 		// (assign line to res)
-
-		// pProj1[0] ← 2pProj1[0] and l0 the tangent ℓ passing 2pProj1[0]
-		pProj1[0].doubleStep(&l0)
+		// qProj1[0] ← 2qProj1[0] and l0 the tangent ℓ passing 2qProj1[0]
+		qProj1[0].doubleStep(&l0)
 		// line evaluation at Q[0] (assign)
-		result.B1.A0.Mul(&l0.r1, &q[0].X)
-		result.B0.A0.Mul(&l0.r0, &q[0].Y)
-		result.B1.A1.Set(&l0.r2)
+		result.B0.A0.Set(&l0.r0)
+		result.B0.A1.Mul(&l0.r1, &p[0].X)
+		result.B1.A1.Mul(&l0.r2, &p[0].Y)
 	}
 
 	// k = 1
 	if n >= 2 {
-		// pProj1[1] ← 2pProj1[1] and l0 the tangent ℓ passing 2pProj1[1]
-		pProj1[1].doubleStep(&l0)
-		// line evaluation at Q[0]
-		l0.r1.Mul(&l0.r1, &q[1].X)
-		l0.r0.Mul(&l0.r0, &q[1].Y)
-		// ℓ × res
-		prodLines = fptower.Mul034By034(&l0.r0, &l0.r1, &l0.r2, &result.B0.A0, &result.B1.A0, &result.B1.A1)
+		// qProj1[1] ← 2qProj1[1] and l0 the tangent ℓ passing 2qProj1[1]
+		qProj1[1].doubleStep(&l0)
+		// line evaluation at Q[1]
+		l0.r1.Mul(&l0.r1, &p[1].X)
+		l0.r2.Mul(&l0.r2, &p[1].Y)
+		prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &result.B0.A0, &result.B0.A1, &result.B1.A1)
 		result.B0.A0 = prodLines[0]
 		result.B0.A1 = prodLines[1]
 		result.B0.A2 = prodLines[2]
-		result.B1.A0 = prodLines[3]
-		result.B1.A1 = prodLines[4]
+		result.B1.A1 = prodLines[3]
+		result.B1.A2 = prodLines[4]
 	}
 
 	// k >= 2
 	for k := 2; k < n; k++ {
-		// pProj1[1] ← 2pProj1[1] and l0 the tangent ℓ passing 2pProj1[1]
-		pProj1[k].doubleStep(&l0)
+		// qProj1[k] ← 2qProj1[k] and l0 the tangent ℓ passing 2qProj1[k]
+		qProj1[k].doubleStep(&l0)
 		// line evaluation at Q[k]
-		l0.r1.Mul(&l0.r1, &q[k].X)
-		l0.r0.Mul(&l0.r0, &q[k].Y)
+		l0.r1.Mul(&l0.r1, &p[k].X)
+		l0.r2.Mul(&l0.r2, &p[k].Y)
 		// ℓ × res
-		result.MulBy034(&l0.r0, &l0.r1, &l0.r2)
+		result.MulBy014(&l0.r0, &l0.r1, &l0.r2)
 	}
 
-	var tmp G1Affine
-	for i := len(loopCounter0) - 3; i >= 1; i-- {
-		// (∏ᵢfᵢ)²
-		// mutualize the square among n Miller loops
+	for i := 187; i >= 1; i-- {
 		result.Square(&result)
 
 		j = loopCounter1[i]*3 + loopCounter0[i]
 
 		for k := 0; k < n; k++ {
-			// pProj1[1] ← 2pProj1[1] and l0 the tangent ℓ passing 2pProj1[1]
-			pProj1[k].doubleStep(&l0)
-			// line evaluation at Q[k]
-			l0.r1.Mul(&l0.r1, &q[k].X)
-			l0.r0.Mul(&l0.r0, &q[k].Y)
+			qProj1[k].doubleStep(&l0)
+			l0.r1.Mul(&l0.r1, &p[k].X)
+			l0.r2.Mul(&l0.r2, &p[k].Y)
 
 			switch j {
 			// cases -4, -2, 2, 4 do not occur, given the static loopCounters
 			case -3:
-				tmp.Neg(&p1[k])
-				// pProj1[k] ← pProj1[k]-p1[k] and
-				// l the line ℓ passing pProj1[k] and -p1[k]
-				pProj1[k].addMixedStep(&l, &tmp)
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q1Neg[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			case -1:
-				tmp.Neg(&p0[k])
-				// pProj1[k] ← pProj1[k]-p0[k] and
-				// l the line ℓ passing pProj1[k] and -p0[k]
-				pProj1[k].addMixedStep(&l, &tmp)
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q0Neg[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			case 0:
-				// ℓ × res
-				result.MulBy034(&l0.r0, &l0.r1, &l0.r2)
+				result.MulBy014(&l0.r0, &l0.r1, &l0.r2)
 			case 1:
-				// pProj1[k] ← pProj1[k]+p0[k] and
-				// l the line ℓ passing pProj1[k] and p0[k]
-				pProj1[k].addMixedStep(&l, &p0[k])
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q0[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			case 3:
-				// pProj1[k] ← pProj1[k]+p1[k] and
-				// l the line ℓ passing pProj1[k] and p1[k]
-				pProj1[k].addMixedStep(&l, &p1[k])
-				// line evaluation at Q[k]
-				l.r1.Mul(&l.r1, &q[k].X)
-				l.r0.Mul(&l.r0, &q[k].Y)
-				// ℓ × ℓ
-				prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-				// (ℓ × ℓ) × res
-				result.MulBy01234(&prodLines)
+				qProj1[k].addMixedStep(&l, &q1[k])
+				l.r1.Mul(&l.r1, &p[k].X)
+				l.r2.Mul(&l.r2, &p[k].Y)
+				prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+				result.MulBy01245(&prodLines)
 			default:
 				return GT{}, errors.New("invalid loopCounter")
 			}
@@ -284,30 +258,23 @@ func MillerLoop(P []G1Affine, Q []G2Affine) (GT, error) {
 	// j = -3
 	result.Square(&result)
 	for k := 0; k < n; k++ {
-		// pProj1[k] ← 2pProj1[k] and l0 the tangent ℓ passing 2pProj1[k]
-		pProj1[k].doubleStep(&l0)
-		// line evaluation at Q[k]
-		l0.r1.Mul(&l0.r1, &q[k].X)
-		l0.r0.Mul(&l0.r0, &q[k].Y)
-
-		tmp.Neg(&p1[k])
-		// l the line passing pProj1[k] and -p1
-		pProj1[k].lineCompute(&l, &tmp)
-		// line evaluation at Q[k]
-		l.r1.Mul(&l.r1, &q[k].X)
-		l.r0.Mul(&l.r0, &q[k].Y)
-		// ℓ × ℓ
-		prodLines = fptower.Mul034By034(&l.r0, &l.r1, &l.r2, &l0.r0, &l0.r1, &l0.r2)
-		// (ℓ × ℓ) × res
-		result.MulBy034(&l0.r0, &l0.r1, &l0.r2)
+		qProj1[k].doubleStep(&l0)
+		l0.r1.Mul(&l0.r1, &p[k].X)
+		l0.r2.Mul(&l0.r2, &p[k].Y)
+		qProj1[k].lineCompute(&l, &q1Neg[k])
+		l.r1.Mul(&l.r1, &p[k].X)
+		l.r2.Mul(&l.r2, &p[k].Y)
+		prodLines = fptower.Mul014By014(&l0.r0, &l0.r1, &l0.r2, &l.r0, &l.r1, &l.r2)
+		result.MulBy01245(&prodLines)
 	}
 
 	return result, nil
+
 }
 
 // doubleStep doubles a point in Homogenous projective coordinates, and evaluates the line in Miller loop
 // https://eprint.iacr.org/2013/722.pdf (Section 4.3)
-func (p *g1Proj) doubleStep(evaluations *lineEvaluation) {
+func (p *g2Proj) doubleStep(evaluations *lineEvaluation) {
 
 	// get some Element from our pool
 	var t1, A, B, C, D, E, EE, F, G, H, I, J, K fp.Element
@@ -317,10 +284,7 @@ func (p *g1Proj) doubleStep(evaluations *lineEvaluation) {
 	C.Square(&p.z)
 	D.Double(&C).
 		Add(&D, &C)
-
-		// E.Mul(&D, &bCurveCoeff)
-	E.Neg(&D)
-
+	E.Mul(&D, &bTwistCurveCoeff)
 	F.Double(&E).
 		Add(&F, &E)
 	G.Add(&B, &F)
@@ -343,15 +307,15 @@ func (p *g1Proj) doubleStep(evaluations *lineEvaluation) {
 	p.z.Mul(&B, &H)
 
 	// Line evaluation
-	evaluations.r0.Neg(&H)
+	evaluations.r0.Set(&I)
 	evaluations.r1.Double(&J).
 		Add(&evaluations.r1, &J)
-	evaluations.r2.Set(&I)
+	evaluations.r2.Neg(&H)
 }
 
 // addMixedStep point addition in Mixed Homogenous projective and Affine coordinates
 // https://eprint.iacr.org/2013/722.pdf (Section 4.3)
-func (p *g1Proj) addMixedStep(evaluations *lineEvaluation, a *G1Affine) {
+func (p *g2Proj) addMixedStep(evaluations *lineEvaluation, a *G2Affine) {
 
 	// get some Element from our pool
 	var Y2Z1, X2Z1, O, L, C, D, E, F, G, H, t0, t1, t2, J fp.Element
@@ -381,14 +345,14 @@ func (p *g1Proj) addMixedStep(evaluations *lineEvaluation, a *G1Affine) {
 		Sub(&J, &t2)
 
 	// Line evaluation
-	evaluations.r0.Set(&L)
+	evaluations.r0.Set(&J)
 	evaluations.r1.Neg(&O)
-	evaluations.r2.Set(&J)
+	evaluations.r2.Set(&L)
 }
 
 // lineCompute computes the line through p in Homogenous projective coordinates
 // and a in affine coordinates. It does not compute the resulting point p+a.
-func (p *g1Proj) lineCompute(evaluations *lineEvaluation, a *G1Affine) {
+func (p *g2Proj) lineCompute(evaluations *lineEvaluation, a *G2Affine) {
 
 	// get some Element from our pool
 	var Y2Z1, X2Z1, O, L, t2, J fp.Element
@@ -401,7 +365,7 @@ func (p *g1Proj) lineCompute(evaluations *lineEvaluation, a *G1Affine) {
 		Sub(&J, &t2)
 
 	// Line evaluation
-	evaluations.r0.Set(&L)
+	evaluations.r0.Set(&J)
 	evaluations.r1.Neg(&O)
-	evaluations.r2.Set(&J)
+	evaluations.r2.Set(&L)
 }
diff --git a/internal/generator/config/bw6-756.go b/internal/generator/config/bw6-756.go
index 4657c6cc4b..21ab364ccf 100644
--- a/internal/generator/config/bw6-756.go
+++ b/internal/generator/config/bw6-756.go
@@ -13,7 +13,6 @@ var BW6_756 = Curve{
 		GLV:              true,
 		CofactorCleaning: true,
 		CRange:           []int{4, 5, 8, 16},
-		Projective:       true,
 	},
 	G2: Point{
 		CoordType:        "fp.Element",
@@ -22,6 +21,7 @@ var BW6_756 = Curve{
 		GLV:              true,
 		CofactorCleaning: true,
 		CRange:           []int{4, 5, 8, 16},
+		Projective:       true,
 	},
 	// 2-isogeny
 	HashE1: &HashSuiteSswu{
diff --git a/internal/generator/config/bw6-761.go b/internal/generator/config/bw6-761.go
index 33f34ae699..aa2ec6c363 100644
--- a/internal/generator/config/bw6-761.go
+++ b/internal/generator/config/bw6-761.go
@@ -13,7 +13,6 @@ var BW6_761 = Curve{
 		GLV:              true,
 		CofactorCleaning: true,
 		CRange:           []int{4, 5, 8, 16},
-		Projective:       true,
 	},
 	G2: Point{
 		CoordType:        "fp.Element",
@@ -22,6 +21,7 @@ var BW6_761 = Curve{
 		GLV:              true,
 		CofactorCleaning: true,
 		CRange:           []int{4, 5, 8, 16},
+		Projective:       true,
 	},
 	// 2-isogeny
 	HashE1: &HashSuiteSswu{
diff --git a/internal/generator/ecc/template/point.go.tmpl b/internal/generator/ecc/template/point.go.tmpl
index dfabf475fd..411f14d489 100644
--- a/internal/generator/ecc/template/point.go.tmpl
+++ b/internal/generator/ecc/template/point.go.tmpl
@@ -1521,52 +1521,6 @@ func (p *{{ $TProjective }}) FromAffine(Q *{{ $TAffine }}) *{{ $TProjective }} {
 	return p
 }
 
-{{- if eq .PointName "g1"}}
-// BatchProjectiveToAffine{{ toUpper .PointName }} converts points in Projective coordinates to Affine coordinates
-// performing a single field inversion (Montgomery batch inversion trick).
-func BatchProjectiveToAffine{{ toUpper .PointName }}(points []{{ $TProjective }}) []{{ $TAffine }} {
-	result := make([]{{ $TAffine }}, len(points))
-	zeroes := make([]bool, len(points))
-	accumulator := fp.One()
-
-	// batch invert all points[].Z coordinates with Montgomery batch inversion trick
-	// (stores points[].Z^-1 in result[i].X to avoid allocating a slice of fr.Elements)
-	for i := 0; i < len(points); i++ {
-		if points[i].z.IsZero() {
-			zeroes[i] = true
-			continue
-		}
-		result[i].X = accumulator
-		accumulator.Mul(&accumulator, &points[i].z)
-	}
-
-	var accInverse fp.Element
-	accInverse.Inverse(&accumulator)
-
-	for i := len(points) - 1; i >= 0; i-- {
-		if zeroes[i] {
-			// do nothing, (X=0, Y=0) is infinity point in affine
-			continue
-		}
-		result[i].X.Mul(&result[i].X, &accInverse)
-		accInverse.Mul(&accInverse, &points[i].z)
-	}
-
-	// batch convert to affine.
-	parallel.Execute(len(points), func(start, end int) {
-		for i := start; i < end; i++ {
-			if zeroes[i] {
-                // do nothing, (X=0, Y=0) is infinity point in affine
-				continue
-			}
-			a := result[i].X
-			result[i].X.Mul(&points[i].x, &a)
-			result[i].Y.Mul(&points[i].y, &a)
-		}
-	})
-    return result
-}
-{{end }}
 {{end }}