perf: optimize class equivalence check for BLS12 final exp

std/algebra/emulated/fields_bls12381/e12_pairing.go

@@ -384,14 +384,15 @@ func (e Ext12) FrobeniusSquareTorus(y *E6) *E6 {
 	return &E6{B0: *t0, B1: *t1, B2: *t2}
 }
 
-// FinalExponentiationCheck checks that a Miller function output x lies in the
+// AssertFinalExponentiationIsOne checks that a Miller function output x lies in the
 // same equivalence class as the reduced pairing. This replaces the final
 // exponentiation step in-circuit.
-// The method follows Section 4 of [On Proving Pairings] paper by A. Novakovic and L. Eagen.
+// The method is inspired from [On Proving Pairings] paper by A. Novakovic and
+// L. Eagen, and is based on a personal communication with A. Novakovic.
 //
 // [On Proving Pairings]: https://eprint.iacr.org/2024/640.pdf
-func (e Ext12) FinalExponentiationCheck(x *E12) *E12 {
-	res, err := e.fp.NewHint(finalExpHint, 12, &x.C0.B0.A0, &x.C0.B0.A1, &x.C0.B1.A0, &x.C0.B1.A1, &x.C0.B2.A0, &x.C0.B2.A1, &x.C1.B0.A0, &x.C1.B0.A1, &x.C1.B1.A0, &x.C1.B1.A1, &x.C1.B2.A0, &x.C1.B2.A1)
+func (e Ext12) AssertFinalExponentiationIsOne(x *E12) {
+	res, err := e.fp.NewHint(finalExpHint, 18, &x.C0.B0.A0, &x.C0.B0.A1, &x.C0.B1.A0, &x.C0.B1.A1, &x.C0.B2.A0, &x.C0.B2.A1, &x.C1.B0.A0, &x.C1.B0.A1, &x.C1.B1.A0, &x.C1.B1.A1, &x.C1.B2.A0, &x.C1.B2.A1)
 	if err != nil {
 		// err is non-nil only for invalid number of inputs
 		panic(err)
@@ -409,21 +410,27 @@ func (e Ext12) FinalExponentiationCheck(x *E12) *E12 {
 			B2: E2{A0: *res[10], A1: *res[11]},
 		},
 	}
+	// constrain cubicNonResiduePower to be in Fp6
+	scalingFactor := E12{
+		C0: E6{
+			B0: E2{A0: *res[12], A1: *res[13]},
+			B1: E2{A0: *res[14], A1: *res[15]},
+			B2: E2{A0: *res[16], A1: *res[17]},
+		},
+		C1: (*e.Ext6.Zero()),
+	}
 
-	// Check that  x == residueWitness^r by checking that:
-	// x^k == residueWitness^(q-u)
-	// where k = (u-1)^2/3, u=-0xd201000000010000 the BLS12-381 seed
-	// and residueWitness from the hint.
+	// Check that  x * scalingFactor == residueWitness^(q-u)
+	// where u=-0xd201000000010000 is the BLS12-381 seed,
+	// and residueWitness, scalingFactor from the hint.
 	t0 := e.Frobenius(&residueWitness)
 	// exponentiation by -u
 	t1 := e.Expt(&residueWitness)
 	t0 = e.Mul(t0, t1)
-	// exponentiation by U=(u-1)^2/3
-	t1 = e.ExpByU(x)
 
-	e.AssertIsEqual(t0, t1)
+	t1 = e.Mul(x, &scalingFactor)
 
-	return nil
+	e.AssertIsEqual(t0, t1)
 }
 
 func (e Ext12) Frobenius(x *E12) *E12 {

std/algebra/emulated/fields_bls12381/hints.go

@@ -271,11 +271,11 @@ func divE12Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) erro
 }
 
 func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
-	// This follows section 4.1 of https://eprint.iacr.org/2024/640.pdf (Th. 1)
+	// 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 millerLoop, residueWitness bls12381.E12
-			var rInv big.Int
+			var millerLoop bls12381.E12
 
 			millerLoop.C0.B0.A0.SetBigInt(inputs[0])
 			millerLoop.C0.B0.A1.SetBigInt(inputs[1])
@@ -290,12 +290,71 @@ func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) er
 			millerLoop.C1.B2.A0.SetBigInt(inputs[10])
 			millerLoop.C1.B2.A1.SetBigInt(inputs[11])
 
-			// compute r-th root:
-			// Exponentiate to rInv where
-			// rInv = 1/r mod (p^12-1)/r
-			rInv.SetString("169662389312441398885310937191698694666993326870281216192803558492181163400934408837135364582394949149589560242411491538960982200559697133935443307582773537814554128992403254243871087441488619811839498788505657962013599019994544063402394719913759780901881538869078447034832302535303591303383830742161317593225991746471557492001710830538428792119562309446698444646787667517629943447802199824630112988907247336627481159245442124709621313522294197747687500252452962523217400829932174349352696726049683687654879009114460723993703760367089269403767790334911644010940272722630305066645230222732316445557889124653426141642271480304669447694344127599708992364443461893123938202386892312748211835322692697497854107961493711137028209148238339237355911496376520814450515612396561384525661635220451168152178239892009375229296874955612623691164738926395993739297557487207643426168321070539996994036837992284584225139752716615623194417718962478029165908544042568334172107008712033983002554672734519081879196926275059798317879322062358113986901925780890205936071364647548199159506709147492864081514759663116291487638998943660232689862634717010538047493292265992334130695994203833154950619462266484292385471162124464248375625748097868775829652908052615424796255913420292818674303286242639225711610323988077268116737", 10)
-			residueWitness.Exp(millerLoop, &rInv)
-
+			var root, rootPthInverse, root27thInverse, residueWitness, scalingFactor bls12381.E12
+			var order3rd, order3rdPower, exponent, exponentInv, finalExpFactor, polyFactor big.Int
+			// 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])
@@ -309,6 +368,14 @@ func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) er
 			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
 		})
 }

std/algebra/emulated/fields_bn254/e12_pairing.go

@@ -422,13 +422,13 @@ func (e Ext12) FrobeniusCubeTorus(y *E6) *E6 {
 	return res
 }
 
-// FinalExponentiationCheck checks that a Miller function output x lies in the
+// AssertFinalExponentiationIsOne checks that a Miller function output x lies in the
 // same equivalence class as the reduced pairing. This replaces the final
 // exponentiation step in-circuit.
 // The method follows Section 4 of [On Proving Pairings] paper by A. Novakovic and L. Eagen.
 //
 // [On Proving Pairings]: https://eprint.iacr.org/2024/640.pdf
-func (e Ext12) FinalExponentiationCheck(x *E12) *E12 {
+func (e Ext12) AssertFinalExponentiationIsOne(x *E12) {
 	res, err := e.fp.NewHint(finalExpHint, 24, &x.C0.B0.A0, &x.C0.B0.A1, &x.C0.B1.A0, &x.C0.B1.A1, &x.C0.B2.A0, &x.C0.B2.A1, &x.C1.B0.A0, &x.C1.B0.A1, &x.C1.B1.A0, &x.C1.B1.A1, &x.C1.B2.A0, &x.C1.B2.A1)
 	if err != nil {
 		// err is non-nil only for invalid number of inputs
@@ -474,8 +474,6 @@ func (e Ext12) FinalExponentiationCheck(x *E12) *E12 {
 	t0 = e.Mul(t0, t1)
 
 	e.AssertIsEqual(t0, t2)
-
-	return nil
 }
 
 func (e Ext12) Frobenius(x *E12) *E12 {

std/algebra/emulated/fields_bw6761/e6_pairing.go

@@ -322,13 +322,13 @@ func (e *Ext6) MulBy02345(z *E6, x [5]*baseEl) *E6 {
 	}
 }
 
-// FinalExponentiationCheck checks that a Miller function output x lies in the
+// AssertFinalExponentiationIsOne checks that a Miller function output x lies in the
 // same equivalence class as the reduced pairing. This replaces the final
 // exponentiation step in-circuit.
 // The method is adapted from Section 4 of [On Proving Pairings] paper by A. Novakovic and L. Eagen.
 //
 // [On Proving Pairings]: https://eprint.iacr.org/2024/640.pdf
-func (e Ext6) FinalExponentiationCheck(x *E6) *E6 {
+func (e Ext6) AssertFinalExponentiationIsOne(x *E6) {
 	res, err := e.fp.NewHint(finalExpHint, 6, &x.A0, &x.A1, &x.A2, &x.A3, &x.A4, &x.A5)
 	if err != nil {
 		// err is non-nil only for invalid number of inputs
@@ -357,8 +357,6 @@ func (e Ext6) FinalExponentiationCheck(x *E6) *E6 {
 	t0 = e.DivUnchecked(t0, t1)
 
 	e.AssertIsEqual(t0, x)
-
-	return nil
 }
 
 // ExpByU2 set z to z^(x₀+1) in E12 and return z

std/algebra/emulated/sw_bls12381/pairing.go

@@ -251,7 +251,7 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 
 	}
 	// We perform the easy part of the final exp to push f to the cyclotomic
-	// subgroup so that FinalExponentiationCheck is carried with optimized
+	// subgroup so that AssertFinalExponentiationIsOne is carried with optimized
 	// cyclotomic squaring (e.g. Karabina12345).
 	//
 	// f = f^(p⁶-1)(p²+1)
@@ -260,7 +260,7 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	f = pr.FrobeniusSquare(buf)
 	f = pr.Mul(f, buf)
 
-	pr.FinalExponentiationCheck(f)
+	pr.AssertFinalExponentiationIsOne(f)
 
 	return nil
 }

std/algebra/emulated/sw_bls12381/pairing_test.go

@@ -177,7 +177,7 @@ func (c *PairingCheckCircuit) Define(api frontend.API) error {
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
-	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In1G1, &c.In2G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2, &c.In1G2, &c.In2G2})
+	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
 	}
@@ -186,10 +186,13 @@ func (c *PairingCheckCircuit) Define(api frontend.API) error {
 
 func TestPairingCheckTestSolve(t *testing.T) {
 	assert := test.NewAssert(t)
+	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
-	_, q2 := randomG1G2Affines()
 	var p2 bls12381.G1Affine
-	p2.Neg(&p1)
+	p2.Double(&p1).Neg(&p2)
+	var q2 bls12381.G2Affine
+	q2.Set(&q1)
+	q1.Double(&q1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),
@@ -228,11 +231,13 @@ func TestGroupMembershipSolve(t *testing.T) {
 
 // bench
 func BenchmarkPairing(b *testing.B) {
-
+	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
-	_, q2 := randomG1G2Affines()
 	var p2 bls12381.G1Affine
-	p2.Neg(&p1)
+	p2.Double(&p1).Neg(&p2)
+	var q2 bls12381.G2Affine
+	q2.Set(&q1)
+	q1.Double(&q1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),

std/algebra/emulated/sw_bn254/pairing.go

@@ -251,7 +251,7 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 
 	}
 	// We perform the easy part of the final exp to push f to the cyclotomic
-	// subgroup so that FinalExponentiationCheck is carried with optimized
+	// subgroup so that AssertFinalExponentiationIsOne is carried with optimized
 	// cyclotomic squaring (e.g. Karabina12345).
 	//
 	// f = f^(p⁶-1)(p²+1)
@@ -260,7 +260,7 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	f = pr.FrobeniusSquare(buf)
 	f = pr.Mul(f, buf)
 
-	pr.FinalExponentiationCheck(f)
+	pr.AssertFinalExponentiationIsOne(f)
 
 	return nil
 }

std/algebra/emulated/sw_bn254/pairing_test.go

@@ -259,7 +259,7 @@ func (c *PairingCheckCircuit) Define(api frontend.API) error {
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
-	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In1G1, &c.In2G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2, &c.In1G2, &c.In2G2})
+	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
 	}
@@ -268,10 +268,13 @@ func (c *PairingCheckCircuit) Define(api frontend.API) error {
 
 func TestPairingCheckTestSolve(t *testing.T) {
 	assert := test.NewAssert(t)
+	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
-	_, q2 := randomG1G2Affines()
 	var p2 bn254.G1Affine
-	p2.Neg(&p1)
+	p2.Double(&p1).Neg(&p2)
+	var q2 bn254.G2Affine
+	q2.Set(&q1)
+	q1.Double(&q1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),
@@ -477,11 +480,13 @@ func TestIsMillerLoopAndFinalExpCircuitTestSolve(t *testing.T) {
 
 // bench
 func BenchmarkPairing(b *testing.B) {
-
+	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
-	_, q2 := randomG1G2Affines()
 	var p2 bn254.G1Affine
-	p2.Neg(&p1)
+	p2.Double(&p1).Neg(&p2)
+	var q2 bn254.G2Affine
+	q2.Set(&q1)
+	q1.Double(&q1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),

std/algebra/emulated/sw_bw6761/pairing.go

@@ -147,7 +147,7 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 
 	}
 	// We perform the easy part of the final exp to push f to the cyclotomic
-	// subgroup so that FinalExponentiationCheck is carried with optimized
+	// subgroup so that AssertFinalExponentiationIsOne is carried with optimized
 	// cyclotomic squaring (e.g. Karabina12345).
 	//
 	// f = f^(p³-1)(p+1)
@@ -156,7 +156,7 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	f = pr.Frobenius(buf)
 	f = pr.Mul(f, buf)
 
-	pr.FinalExponentiationCheck(f)
+	pr.AssertFinalExponentiationIsOne(f)
 
 	return nil
 }