perf: optimize pre-computation in fixed-argument pairings by batching double steps

ecc/bls12-377/pairing.go

@@ -367,18 +367,327 @@ func PairingCheckFixedQ(P []G1Affine, lines [][2][len(LoopCounter) - 1]LineEvalu
 	return f.Equal(&one), nil
 }
 
+// manyDoubleSteps performs k consecutive doublings on p and returns the line evaluations.
+// It uses a recurrence to compute 2^k*P with a single batch inversion.
+func (p *G2Affine) manyDoubleSteps(k int, evaluations []LineEvaluationAff) {
+	if k == 0 {
+		return
+	}
+
+	// Step 1: Compute A[i], B[i], C[i] using the recurrence
+	A := make([]fptower.E2, k+1)
+	B := make([]fptower.E2, k+1)
+	C := make([]fptower.E2, k+1)
+
+	var tmp fptower.E2
+	A[0].Set(&p.X)
+	C[0].Neg(&p.Y)
+	tmp.Square(&p.X)
+	B[0].Double(&tmp).Add(&B[0], &tmp) // B[0] = 3x²
+
+	for i := 1; i <= k; i++ {
+		var Csq, ACs, eightACs fptower.E2
+		Csq.Square(&C[i-1])
+		ACs.Mul(&A[i-1], &Csq)
+		eightACs.Double(&ACs).Double(&eightACs).Double(&eightACs)
+		A[i].Square(&B[i-1]).Sub(&A[i], &eightACs)
+
+		tmp.Square(&A[i])
+		B[i].Double(&tmp).Add(&B[i], &tmp)
+
+		var C4, fourACs, diff fptower.E2
+		C4.Square(&Csq)
+		fourACs.Double(&ACs).Double(&fourACs)
+		diff.Sub(&A[i], &fourACs)
+		C[i].Double(&C4).Double(&C[i]).Double(&C[i]) // 8*C[i-1]⁴
+		tmp.Mul(&B[i-1], &diff)
+		C[i].Add(&C[i], &tmp) // C[i] = 8*C[i-1]⁴ + B[i-1]*(A[i] - 4*A[i-1]*C[i-1]²)
+	}
+
+	// Step 2: Compute D[i] = -2*C[i] = 2*y[i] for i = 0..k-1
+	D := make([]fptower.E2, k)
+	for i := 0; i < k; i++ {
+		D[i].Double(&C[i]).Neg(&D[i])
+	}
+
+	// Step 3: Compute T[i] = D[0]*D[1]*...*D[i] for i = 0..k-1
+	T := make([]fptower.E2, k)
+	T[0].Set(&D[0])
+	for i := 1; i < k; i++ {
+		T[i].Mul(&T[i-1], &D[i])
+	}
+
+	// Step 4: Batch invert T
+	invT := fptower.BatchInvertE2(T)
+
+	// Step 5: Compute line evaluations
+	// For i = 0: x[0] = A[0], y[0] = -C[0]
+	// For i > 0: x[i] = A[i] / T[i-1]², y[i] = -C[i] / T[i-1]³
+
+	// Step 0: special case since scaling is 1
+	evaluations[0].R0.Mul(&B[0], &invT[0])
+	evaluations[0].R1.Mul(&B[0], &A[0]).Mul(&evaluations[0].R1, &invT[0]).Add(&evaluations[0].R1, &C[0])
+
+	// Steps 1 to k-1
+	var invT2, invT3 fptower.E2
+	for i := 1; i < k; i++ {
+		// R0 = B[i] / T[i]
+		evaluations[i].R0.Mul(&B[i], &invT[i])
+
+		// R1 = B[i]*A[i]/(T[i]*T[i-1]²) + C[i]/T[i-1]³
+		invT2.Square(&invT[i-1])
+		invT3.Mul(&invT2, &invT[i-1])
+
+		var term1, term2 fptower.E2
+		term1.Mul(&B[i], &A[i]).Mul(&term1, &invT[i]).Mul(&term1, &invT2)
+		term2.Mul(&C[i], &invT3)
+		evaluations[i].R1.Add(&term1, &term2)
+	}
+
+	// Step 6: Final point coordinates
+	// x[k] = A[k] / T[k-1]²
+	// y[k] = -C[k] / T[k-1]³
+	invT2.Square(&invT[k-1])
+	invT3.Mul(&invT2, &invT[k-1])
+	p.X.Mul(&A[k], &invT2)
+	p.Y.Mul(&C[k], &invT3).Neg(&p.Y)
+}
+
+// manyDoublesAndAdd performs k doublings followed by an add with point a,
+// using only a single batch inversion for all operations combined.
+// This fuses manyDoubleSteps(k) + doubleAndAddStep into one inversion.
+func (p *G2Affine) manyDoublesAndAdd(k int, doubEvals []LineEvaluationAff, addEval1, addEval2 *LineEvaluationAff, a *G2Affine) {
+	if k == 0 {
+		p.doubleAndAddStep(addEval1, addEval2, a)
+		return
+	}
+
+	// Step 1: Compute A[i], B[i], C[i] using the recurrence (same as manyDoubleSteps)
+	A := make([]fptower.E2, k+1)
+	B := make([]fptower.E2, k+1)
+	C := make([]fptower.E2, k+1)
+
+	var tmp fptower.E2
+	A[0].Set(&p.X)
+	C[0].Neg(&p.Y) // C[i] = -y[i] (scaled)
+	tmp.Square(&p.X)
+	B[0].Double(&tmp).Add(&B[0], &tmp) // B[0] = 3x²
+
+	for i := 1; i <= k; i++ {
+		var Csq, ACs, eightACs fptower.E2
+		Csq.Square(&C[i-1])
+		ACs.Mul(&A[i-1], &Csq)
+		eightACs.Double(&ACs).Double(&eightACs).Double(&eightACs)
+		A[i].Square(&B[i-1]).Sub(&A[i], &eightACs)
+
+		tmp.Square(&A[i])
+		B[i].Double(&tmp).Add(&B[i], &tmp)
+
+		var C4, fourACs, diff fptower.E2
+		C4.Square(&Csq)
+		fourACs.Double(&ACs).Double(&fourACs)
+		diff.Sub(&A[i], &fourACs)
+		C[i].Double(&C4).Double(&C[i]).Double(&C[i])
+		tmp.Mul(&B[i-1], &diff)
+		C[i].Add(&C[i], &tmp)
+	}
+
+	// Step 2: Compute D[i] = -2*C[i] for i = 0..k-1 (denominators for doubling lines)
+	D := make([]fptower.E2, k)
+	for i := 0; i < k; i++ {
+		D[i].Double(&C[i]).Neg(&D[i])
+	}
+
+	// Step 3: Compute partial products T[i] = D[0]*D[1]*...*D[i]
+	T := make([]fptower.E2, k+1)
+	T[0].Set(&D[0])
+	for i := 1; i < k; i++ {
+		T[i].Mul(&T[i-1], &D[i])
+	}
+
+	// Step 4: Compute ELM formula numerators using scaled coordinates
+	// At position k, we have:
+	//   x_P = A[k] / S² where S = T[k-1] = D[0]*...*D[k-1]
+	//   y_P = -C[k] / S³
+	// For ELM with point a = (a.X, a.Y):
+	//   A_num = x_P - a.X = (A[k] - a.X * S²) / S²
+	//   B_num = y_P - a.Y = (-C[k] - a.Y * S³) / S³
+	var S, S2, S3 fptower.E2
+	S.Set(&T[k-1])
+	S2.Square(&S)
+	S3.Mul(&S2, &S)
+
+	var A_num, B_num fptower.E2
+	A_num.Mul(&a.X, &S2)
+	A_num.Sub(&A[k], &A_num) // A_num = A[k] - a.X * S²
+
+	B_num.Mul(&a.Y, &S3)
+	B_num.Add(&C[k], &B_num)
+	B_num.Neg(&B_num) // B_num = -C[k] - a.Y * S³
+
+	// U_num = B_num² - (2*A[k] + a.X*S²) * A_num²
+	var A_num2, B_num2, twoAk, axS2, coeff, U_num fptower.E2
+	A_num2.Square(&A_num)
+	B_num2.Square(&B_num)
+	twoAk.Double(&A[k])
+	axS2.Mul(&a.X, &S2)
+	coeff.Add(&twoAk, &axS2)
+	U_num.Mul(&coeff, &A_num2)
+	U_num.Sub(&B_num2, &U_num)
+
+	// Step 5: Extend T with the ELM factor
+	// T[k] = S * A_num * U_num for batch inversion
+	T[k].Mul(&S, &A_num)
+	T[k].Mul(&T[k], &U_num)
+
+	// Step 6: Batch invert all of T[0..k]
+	invT := fptower.BatchInvertE2(T)
+
+	// Step 7: Compute doubling line evaluations (same as manyDoubleSteps)
+	// Step 0: special case since scaling is 1
+	doubEvals[0].R0.Mul(&B[0], &invT[0])
+	doubEvals[0].R1.Mul(&B[0], &A[0]).Mul(&doubEvals[0].R1, &invT[0]).Add(&doubEvals[0].R1, &C[0])
+
+	// Steps 1 to k-1
+	var invT2, invT3 fptower.E2
+	for i := 1; i < k; i++ {
+		doubEvals[i].R0.Mul(&B[i], &invT[i])
+
+		invT2.Square(&invT[i-1])
+		invT3.Mul(&invT2, &invT[i-1])
+
+		var term1, term2 fptower.E2
+		term1.Mul(&B[i], &A[i]).Mul(&term1, &invT[i]).Mul(&term1, &invT2)
+		term2.Mul(&C[i], &invT3)
+		doubEvals[i].R1.Add(&term1, &term2)
+	}
+
+	// Step 8: Compute point P coordinates at position k
+	// x_P = A[k] / S², y_P = -C[k] / S³
+	var x_P, y_P fptower.E2
+	invT2.Square(&invT[k-1])
+	invT3.Mul(&invT2, &invT[k-1])
+	x_P.Mul(&A[k], &invT2)
+	y_P.Mul(&C[k], &invT3).Neg(&y_P)
+
+	// Step 9: Compute ELM slopes using the batch inversion result
+	// invT[k] = 1 / (S * A_num * U_num)
+	// invSA = U_num * invT[k] = 1 / (S * A_num)
+	// invSU = A_num * invT[k] = 1 / (S * U_num)
+	var invSA, invSU fptower.E2
+	invSA.Mul(&U_num, &invT[k])
+	invSU.Mul(&A_num, &invT[k])
+
+	// λ1 = B_num / (S * A_num) = B_num * invSA
+	// But we need λ1 = (y_P - a.Y) / (x_P - a.X)
+	// B_num/S³ = y_P - a.Y, A_num/S² = x_P - a.X
+	// So λ1 = (B_num/S³) / (A_num/S²) = B_num / (S * A_num) = B_num * invSA
+	var l1 fptower.E2
+	l1.Mul(&B_num, &invSA)
+
+	// x3 = λ1² - x_P - a.X
+	var x3 fptower.E2
+	x3.Square(&l1)
+	x3.Sub(&x3, &x_P)
+	x3.Sub(&x3, &a.X)
+
+	// Line 1 evaluation: R0 = λ1, R1 = λ1*x_P - y_P
+	addEval1.R0.Set(&l1)
+	addEval1.R1.Mul(&l1, &x_P)
+	addEval1.R1.Sub(&addEval1.R1, &y_P)
+
+	// λ2 = -λ1 - 2*y_P*A²/U
+	// In scaled coordinates:
+	//   y_P = -C[k]/S³, A² = A_num²/S⁴, U = U_num/S⁶
+	//   2*y_P*A²/U = 2*(-C[k]/S³)*(A_num²/S⁴)*(S⁶/U_num)
+	//             = -2*C[k]*A_num²/(S*U_num)
+	//             = -2*C[k]*A_num²*invSU
+	var l2 fptower.E2
+	l2.Double(&C[k])
+	l2.Neg(&l2)          // -2*C[k]
+	l2.Mul(&l2, &A_num2) // -2*C[k]*A_num²
+	l2.Mul(&l2, &invSU)  // -2*C[k]*A_num²/(S*U_num) = 2*y_P*A²/U
+	l2.Add(&l2, &l1)
+	l2.Neg(&l2)
+
+	// x4 = λ2² - x_P - x3
+	var x4, y4 fptower.E2
+	x4.Square(&l2)
+	x4.Sub(&x4, &x_P)
+	x4.Sub(&x4, &x3)
+
+	// y4 = λ2*(x_P - x4) - y_P
+	y4.Sub(&x_P, &x4)
+	y4.Mul(&l2, &y4)
+	y4.Sub(&y4, &y_P)
+
+	// Line 2 evaluation: R0 = λ2, R1 = λ2*x_P - y_P
+	addEval2.R0.Set(&l2)
+	addEval2.R1.Mul(&l2, &x_P)
+	addEval2.R1.Sub(&addEval2.R1, &y_P)
+
+	// Update point
+	p.X.Set(&x4)
+	p.Y.Set(&y4)
+}
+
 // PrecomputeLines precomputes the lines for the fixed-argument Miller loop
 func PrecomputeLines(Q G2Affine) (PrecomputedLines [2][len(LoopCounter) - 1]LineEvaluationAff) {
 	var accQ G2Affine
 	accQ.Set(&Q)
 
-	for i := len(LoopCounter) - 2; i >= 0; i-- {
-		if LoopCounter[i] == 0 {
-			accQ.doubleStep(&PrecomputedLines[0][i])
-		} else {
-			accQ.doubleAndAddStep(&PrecomputedLines[0][i], &PrecomputedLines[1][i], &Q)
+	// LoopCounter non-zero values: 0(1), 46(1), 47(1), 51(1), 56(1), 58(1), 63(1)
+	// Optimization: Use manyDoublesAndAdd to fuse consecutive doublings with the following add
+	// This reduces total inversions from 11 to 7.
+
+	// i=62→58: 4 consecutive zeros followed by add at i=58
+	// Uses manyDoublesAndAdd(4) = 1 inversion instead of manyDoubleSteps(4) + doubleAndAddStep = 2 inversions
+	{
+		var doubEvals [4]LineEvaluationAff
+		accQ.manyDoublesAndAdd(4, doubEvals[:], &PrecomputedLines[0][58], &PrecomputedLines[1][58], &Q)
+		PrecomputedLines[0][62] = doubEvals[0]
+		PrecomputedLines[0][61] = doubEvals[1]
+		PrecomputedLines[0][60] = doubEvals[2]
+		PrecomputedLines[0][59] = doubEvals[3]
+	}
+
+	// i=57: LoopCounter[57]=0 (isolated between two 1s)
+	accQ.doubleStep(&PrecomputedLines[0][57])
+
+	// i=56: LoopCounter[56]=1
+	accQ.doubleAndAddStep(&PrecomputedLines[0][56], &PrecomputedLines[1][56], &Q)
+
+	// i=55→51: 4 consecutive zeros followed by add at i=51
+	{
+		var doubEvals [4]LineEvaluationAff
+		accQ.manyDoublesAndAdd(4, doubEvals[:], &PrecomputedLines[0][51], &PrecomputedLines[1][51], &Q)
+		PrecomputedLines[0][55] = doubEvals[0]
+		PrecomputedLines[0][54] = doubEvals[1]
+		PrecomputedLines[0][53] = doubEvals[2]
+		PrecomputedLines[0][52] = doubEvals[3]
+	}
+
+	// i=50→47: 3 consecutive zeros followed by add at i=47
+	{
+		var doubEvals [3]LineEvaluationAff
+		accQ.manyDoublesAndAdd(3, doubEvals[:], &PrecomputedLines[0][47], &PrecomputedLines[1][47], &Q)
+		PrecomputedLines[0][50] = doubEvals[0]
+		PrecomputedLines[0][49] = doubEvals[1]
+		PrecomputedLines[0][48] = doubEvals[2]
+	}
+
+	// i=46: LoopCounter[46]=1 (consecutive 1 after i=47)
+	accQ.doubleAndAddStep(&PrecomputedLines[0][46], &PrecomputedLines[1][46], &Q)
+
+	// i=45→0: 45 consecutive zeros followed by add at i=0
+	{
+		var doubEvals [45]LineEvaluationAff
+		accQ.manyDoublesAndAdd(45, doubEvals[:], &PrecomputedLines[0][0], &PrecomputedLines[1][0], &Q)
+		for j := 0; j < 45; j++ {
+			PrecomputedLines[0][45-j] = doubEvals[j]
 		}
 	}
+
 	return PrecomputedLines
 }