perf: PairingCheck for BN254, BLS12-381, BLS12-377 and BW6-761

std/algebra/emulated/sw_bls12381/hints.go

@@ -14,7 +14,135 @@ func init() {
 
 // GetHints returns all hint functions used in the package.
 func GetHints() []solver.Hint {
-	return []solver.Hint{finalExpHint}
+	return []solver.Hint{
+		finalExpHint,
+		pairingCheckHint,
+	}
+}
+
+func pairingCheckHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	// This is inspired from https://eprint.iacr.org/2024/640.pdf
+	// and based on a personal communication with the author Andrija Novakovic.
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var root, rootPthInverse, root27thInverse, residueWitness, scalingFactor bls12381.E12
+			var order3rd, order3rdPower, exponent, exponentInv, finalExpFactor, polyFactor big.Int
+			var P bls12381.G1Affine
+			var Q bls12381.G2Affine
+			n := len(inputs)
+			p := make([]bls12381.G1Affine, 0, n/6)
+			q := make([]bls12381.G2Affine, 0, n/6)
+			for k := 0; k < n/6+1; k += 2 {
+				P.X.SetBigInt(inputs[k])
+				P.Y.SetBigInt(inputs[k+1])
+				p = append(p, P)
+			}
+			for k := n / 3; k < n/2+3; k += 4 {
+				Q.X.A0.SetBigInt(inputs[k])
+				Q.X.A1.SetBigInt(inputs[k+1])
+				Q.Y.A0.SetBigInt(inputs[k+2])
+				Q.Y.A1.SetBigInt(inputs[k+3])
+				q = append(q, Q)
+			}
+
+			lines := make([][2][len(bls12381.LoopCounter) - 1]bls12381.LineEvaluationAff, 0, len(q))
+			for _, qi := range q {
+				lines = append(lines, bls12381.PrecomputeLines(qi))
+			}
+			millerLoop, err := bls12381.MillerLoopFixedQ(p, lines)
+			if err != nil {
+				return err
+			}
+			millerLoop.Conjugate(&millerLoop)
+
+			// polyFactor = (1-x)/3
+			polyFactor.SetString("5044125407647214251", 10)
+			// finalExpFactor = ((q^12 - 1) / r) / (27 * polyFactor)
+			finalExpFactor.SetString("2366356426548243601069753987687709088104621721678962410379583120840019275952471579477684846670499039076873213559162845121989217658133790336552276567078487633052653005423051750848782286407340332979263075575489766963251914185767058009683318020965829271737924625612375201545022326908440428522712877494557944965298566001441468676802477524234094954960009227631543471415676620753242466901942121887152806837594306028649150255258504417829961387165043999299071444887652375514277477719817175923289019181393803729926249507024121957184340179467502106891835144220611408665090353102353194448552304429530104218473070114105759487413726485729058069746063140422361472585604626055492939586602274983146215294625774144156395553405525711143696689756441298365274341189385646499074862712688473936093315628166094221735056483459332831845007196600723053356837526749543765815988577005929923802636375670820616189737737304893769679803809426304143627363860243558537831172903494450556755190448279875942974830469855835666815454271389438587399739607656399812689280234103023464545891697941661992848552456326290792224091557256350095392859243101357349751064730561345062266850238821755009430903520645523345000326783803935359711318798844368754833295302563158150573540616830138810935344206231367357992991289265295323280", 10)
+
+			// 1. get pth-root inverse
+			exponent.Mul(&finalExpFactor, big.NewInt(27))
+			root.Exp(millerLoop, &exponent)
+			if root.IsOne() {
+				rootPthInverse.SetOne()
+			} else {
+				exponentInv.ModInverse(&exponent, &polyFactor)
+				exponent.Neg(&exponentInv).Mod(&exponent, &polyFactor)
+				rootPthInverse.Exp(root, &exponent)
+			}
+
+			// 2.1. get order of 3rd primitive root
+			var three big.Int
+			three.SetUint64(3)
+			exponent.Mul(&polyFactor, &finalExpFactor)
+			root.Exp(millerLoop, &exponent)
+			if root.IsOne() {
+				order3rdPower.SetUint64(0)
+			}
+			root.Exp(root, &three)
+			if root.IsOne() {
+				order3rdPower.SetUint64(1)
+			}
+			root.Exp(root, &three)
+			if root.IsOne() {
+				order3rdPower.SetUint64(2)
+			}
+			root.Exp(root, &three)
+			if root.IsOne() {
+				order3rdPower.SetUint64(3)
+			}
+
+			// 2.2. get 27th root inverse
+			if order3rdPower.Uint64() == 0 {
+				root27thInverse.SetOne()
+			} else {
+				order3rd.Exp(&three, &order3rdPower, nil)
+				exponent.Mul(&polyFactor, &finalExpFactor)
+				root.Exp(millerLoop, &exponent)
+				exponentInv.ModInverse(&exponent, &order3rd)
+				exponent.Neg(&exponentInv).Mod(&exponent, &order3rd)
+				root27thInverse.Exp(root, &exponent)
+			}
+
+			// 2.3. shift the Miller loop result so that millerLoop * scalingFactor
+			// is of order finalExpFactor
+			scalingFactor.Mul(&rootPthInverse, &root27thInverse)
+			millerLoop.Mul(&millerLoop, &scalingFactor)
+
+			// 3. get the witness residue
+			//
+			// lambda = q - u, the optimal exponent
+			var lambda big.Int
+			lambda.SetString("4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129030796414117214202539", 10)
+			exponent.ModInverse(&lambda, &finalExpFactor)
+			residueWitness.Exp(millerLoop, &exponent)
+
+			// return the witness residue
+			residueWitness.C0.B0.A0.BigInt(outputs[0])
+			residueWitness.C0.B0.A1.BigInt(outputs[1])
+			residueWitness.C0.B1.A0.BigInt(outputs[2])
+			residueWitness.C0.B1.A1.BigInt(outputs[3])
+			residueWitness.C0.B2.A0.BigInt(outputs[4])
+			residueWitness.C0.B2.A1.BigInt(outputs[5])
+			residueWitness.C1.B0.A0.BigInt(outputs[6])
+			residueWitness.C1.B0.A1.BigInt(outputs[7])
+			residueWitness.C1.B1.A0.BigInt(outputs[8])
+			residueWitness.C1.B1.A1.BigInt(outputs[9])
+			residueWitness.C1.B2.A0.BigInt(outputs[10])
+			residueWitness.C1.B2.A1.BigInt(outputs[11])
+
+			// return the scaling factor
+			scalingFactor.C0.B0.A0.BigInt(outputs[12])
+			scalingFactor.C0.B0.A1.BigInt(outputs[13])
+			scalingFactor.C0.B1.A0.BigInt(outputs[14])
+			scalingFactor.C0.B1.A1.BigInt(outputs[15])
+			scalingFactor.C0.B2.A0.BigInt(outputs[16])
+			scalingFactor.C0.B2.A1.BigInt(outputs[17])
+
+			return nil
+
+		})
+
 }
 
 func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {

std/algebra/emulated/sw_bls12381/pairing.go

@@ -100,14 +100,122 @@ func (pr Pairing) Pair(P []*G1Affine, Q []*G2Affine) (*GTEl, error) {
 //
 // This function doesn't check that the inputs are in the correct subgroups.
 func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
-	f, err := pr.MillerLoop(P, Q)
+	// check input size match
+	nP := len(P)
+	nQ := len(Q)
+	if nP == 0 || nP != nQ {
+		return nil
+	}
+	// hint the non-residue witness
+	inputs := make([]*baseEl, 0, 2*nP+4*nQ)
+	for _, p := range P {
+		inputs = append(inputs, &p.X, &p.Y)
+	}
+	for _, q := range Q {
+		inputs = append(inputs, &q.P.X.A0, &q.P.X.A1, &q.P.Y.A0, &q.P.Y.A1)
+	}
+	hint, err := pr.curveF.NewHint(pairingCheckHint, 18, inputs...)
 	if err != nil {
-		return err
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
 
+	residueWitness := pr.FromTower([12]*baseEl{hint[0], hint[1], hint[2], hint[3], hint[4], hint[5], hint[6], hint[7], hint[8], hint[9], hint[10], hint[11]})
+	// constrain cubicNonResiduePower to be in Fp6
+	// that is: a100=a101=a110=a111=a120=a121=0
+	// or
+	//     A0  =  a000 - a001
+	//     A1  =  0
+	//     A2  =  a010 - a011
+	//     A3  =  0
+	//     A4  =  a020 - a021
+	//     A5  =  0
+	//     A6  =  a001
+	//     A7  =  0
+	//     A8  =  a011
+	//     A9  =  0
+	//     A10 =  a021
+	//     A11 =  0
+	scalingFactor := GTEl{
+		A0:  *pr.curveF.Sub(hint[12], hint[13]),
+		A1:  *pr.curveF.Zero(),
+		A2:  *pr.curveF.Sub(hint[14], hint[15]),
+		A3:  *pr.curveF.Zero(),
+		A4:  *pr.curveF.Sub(hint[16], hint[17]),
+		A5:  *pr.curveF.Zero(),
+		A6:  *hint[13],
+		A7:  *pr.curveF.Zero(),
+		A8:  *hint[15],
+		A9:  *pr.curveF.Zero(),
+		A10: *hint[17],
+		A11: *pr.curveF.Zero(),
 	}
 
-	pr.AssertFinalExponentiationIsOne(f)
+	lines := make([]lineEvaluations, nQ)
+	for i := range Q {
+		if Q[i].Lines == nil {
+			Qlines := pr.computeLines(&Q[i].P)
+			Q[i].Lines = &Qlines
+		}
+		lines[i] = *Q[i].Lines
+	}
 
+	// precomputations
+	yInv := make([]*baseEl, nP)
+	xNegOverY := make([]*baseEl, nP)
+
+	for k := 0; k < nP; k++ {
+		// P are supposed to be on G1 respectively of prime order r.
+		// 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.Mul(&P[k].X, yInv[k])
+		xNegOverY[k] = pr.curveF.Neg(xNegOverY[k])
+	}
+
+	// init Miller loop accumulator to residueWitnessInv to share the squarings
+	// of residueWitnessInv^{x₀}
+	residueWitnessInv := pr.Ext12.Inverse(residueWitness)
+	res := residueWitnessInv
+
+	// Compute ∏ᵢ { fᵢ_{x₀,Q}(P) }
+	for i := 62; i >= 0; i-- {
+		// mutualize the square among n Miller loops
+		// (∏ᵢfᵢ)²
+		res = pr.Ext12.Square(res)
+
+		if loopCounter[i] == 0 {
+			for k := 0; k < nP; k++ {
+				res = pr.MulBy02368(res,
+					pr.MulByElement(&lines[k][0][i].R1, yInv[k]),
+					pr.MulByElement(&lines[k][0][i].R0, xNegOverY[k]),
+				)
+			}
+		} else {
+			// multiply by residueWitnessInv when bit=1
+			res = pr.Ext12.Mul(res, residueWitnessInv)
+			for k := 0; k < nP; k++ {
+				res = pr.MulBy02368(res,
+					pr.MulByElement(&lines[k][0][i].R1, yInv[k]),
+					pr.MulByElement(&lines[k][0][i].R0, xNegOverY[k]),
+				)
+				res = pr.MulBy02368(res,
+					pr.MulByElement(&lines[k][1][i].R1, yInv[k]),
+					pr.MulByElement(&lines[k][1][i].R0, xNegOverY[k]),
+				)
+			}
+		}
+	}
+
+	// Check that  res * scalingFactor == residueWitness^(q)
+	// where u=-0xd201000000010000 is the BLS12-381 seed,
+	// and residueWitness, scalingFactor from the hint.
+	// Note that res is already MillerLoop(P,Q) * residueWitnessInv^{-x₀} since
+	// we initialized the Miller loop accumulator with residueWitnessInv.
+	t0 := pr.Frobenius(residueWitness)
+	t1 := pr.Ext12.Mul(res, &scalingFactor)
+
+	pr.AssertIsEqual(t0, t1)
 	return nil
 }
 
@@ -312,7 +420,7 @@ func (pr Pairing) FinalExponentiation(e *GTEl) *GTEl {
 func (pr Pairing) AssertFinalExponentiationIsOne(x *GTEl) {
 	tower := pr.ToTower(x)
 
-	res, err := pr.curveF.NewHint(finalExpHint, 24, tower[0], tower[1], tower[2], tower[3], tower[4], tower[5], tower[6], tower[7], tower[8], tower[9], tower[10], tower[11])
+	res, err := pr.curveF.NewHint(finalExpHint, 18, tower[0], tower[1], tower[2], tower[3], tower[4], tower[5], tower[6], tower[7], tower[8], tower[9], tower[10], tower[11])
 	if err != nil {
 		// err is non-nil only for invalid number of inputs
 		panic(err)

std/algebra/emulated/sw_bls12381/pairing_test.go

@@ -33,6 +33,47 @@ func randomG1G2Affines() (bls12381.G1Affine, bls12381.G2Affine) {
 	return p, q
 }
 
+type MillerLoopCircuit struct {
+	In1G1, In2G1 G1Affine
+	In1G2, In2G2 G2Affine
+	Res          GTEl
+}
+
+func (c *MillerLoopCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	if err != nil {
+		return fmt.Errorf("new pairing: %w", err)
+	}
+	res, err := pairing.MillerLoop([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2})
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	pairing.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMillerLoopTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p1, q1 := randomG1G2Affines()
+	p2, q2 := randomG1G2Affines()
+	lines1 := bls12381.PrecomputeLines(q1)
+	lines2 := bls12381.PrecomputeLines(q2)
+	res, err := bls12381.MillerLoopFixedQ(
+		[]bls12381.G1Affine{p1, p2},
+		[][2][len(bls12381.LoopCounter) - 1]bls12381.LineEvaluationAff{lines1, lines2},
+	)
+	assert.NoError(err)
+	witness := MillerLoopCircuit{
+		In1G1: NewG1Affine(p1),
+		In1G2: NewG2Affine(q1),
+		In2G1: NewG1Affine(p2),
+		In2G2: NewG2Affine(q2),
+		Res:   NewGTEl(res),
+	}
+	err = test.IsSolved(&MillerLoopCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
 type FinalExponentiationCircuit struct {
 	InGt GTEl
 	Res  GTEl