From bdcdb275bd1f3c0b2eca5026dcb17a5979fbca80 Mon Sep 17 00:00:00 2001
From: Youssef El Housni <youssef.housni21@gmail.com>
Date: Mon, 15 Apr 2024 18:31:58 -0400
Subject: [PATCH] perf(pairing/bw6): optimize Miller loop

---
 .../emulated/fields_bw6761/e6_pairing.go      | 27 ++++++++++++++++
 std/algebra/emulated/sw_bw6761/pairing.go     | 32 ++++++++++++++-----
 2 files changed, 51 insertions(+), 8 deletions(-)

diff --git a/std/algebra/emulated/fields_bw6761/e6_pairing.go b/std/algebra/emulated/fields_bw6761/e6_pairing.go
index e4926c064f..20e8a9c444 100644
--- a/std/algebra/emulated/fields_bw6761/e6_pairing.go
+++ b/std/algebra/emulated/fields_bw6761/e6_pairing.go
@@ -198,6 +198,33 @@ func (e Ext6) Mul014By014(d0, d1, c0, c1 *baseEl) [5]*baseEl {
 	return [5]*baseEl{zC0B0, x01, x1, x04, x14}
 }
 
+// MulBy01245 multiplies z by an E6 sparse element of the form
+//
+//	E6{
+//		B0: E3{A0: c0, A1: c1, A2: c2},
+//		B1: E3{A0: 0, A1: c4, A2: c5},
+//	}
+func (e *Ext6) MulBy01245(z *E6, x [5]*baseEl) *E6 {
+	c0 := &E3{A0: *x[0], A1: *x[1], A2: *x[2]}
+	a := e.Ext3.Add(&z.B0, &z.B1)
+	b := &E3{
+		A0: c0.A0,
+		A1: *e.fp.Add(&c0.A1, x[3]),
+		A2: *e.fp.Add(&c0.A2, x[4]),
+	}
+	a = e.Ext3.Mul(a, b)
+	b = e.Ext3.Mul(&z.B0, c0)
+	c := e.Ext3.MulBy12(&z.B1, x[3], x[4])
+	z1 := e.Ext3.Sub(a, b)
+	z1 = e.Ext3.Sub(z1, c)
+	z0 := e.Ext3.MulByNonResidue(c)
+	z0 = e.Ext3.Add(z0, b)
+	return &E6{
+		B0: *z0,
+		B1: *z1,
+	}
+}
+
 // Mul01245By014 multiplies two E6 sparse element of the form
 //
 //	E6{
diff --git a/std/algebra/emulated/sw_bw6761/pairing.go b/std/algebra/emulated/sw_bw6761/pairing.go
index c032b8f6b7..d5b0ac4fa7 100644
--- a/std/algebra/emulated/sw_bw6761/pairing.go
+++ b/std/algebra/emulated/sw_bw6761/pairing.go
@@ -358,19 +358,35 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 		// (∏ᵢfᵢ)²
 		result = pr.Square(result)
 
-		for k := 0; k < n; k++ {
-			result = pr.MulBy014(result,
-				pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]),
-				pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]),
-			)
-		}
-
 		if i > 0 && loopCounter2[i]*3+loopCounter1[i] != 0 {
 			for k := 0; k < n; k++ {
-				result = pr.MulBy014(result,
+				prodLines = pr.Mul014By014(
+					pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]),
+					pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]),
 					pr.curveF.Mul(&lines[k][1][i].R1, yInv[k]),
 					pr.curveF.Mul(&lines[k][1][i].R0, xNegOverY[k]),
 				)
+				result = pr.MulBy01245(result, prodLines)
+			}
+		} else {
+			// if number of lines is odd, mul last line by res
+			// works for n=1 as well
+			if n%2 != 0 {
+				// ℓ × res
+				result = pr.MulBy014(result,
+					pr.curveF.Mul(&lines[n-1][0][i].R1, yInv[n-1]),
+					pr.curveF.Mul(&lines[n-1][0][i].R0, xNegOverY[n-1]),
+				)
+			}
+			// mul lines 2-by-2
+			for k := 1; k < n; k += 2 {
+				prodLines = pr.Mul014By014(
+					pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]),
+					pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]),
+					pr.curveF.Mul(&lines[k-1][0][i].R1, yInv[k-1]),
+					pr.curveF.Mul(&lines[k-1][0][i].R0, xNegOverY[k-1]),
+				)
+				result = pr.MulBy01245(result, prodLines)
 			}
 		}
 	}