perf(bls12-381): eliminate finalexp ~naively

std/algebra/emulated/fields_bls12381/e12.go

@@ -114,6 +114,65 @@ func (e Ext12) Square(x *E12) *E12 {
 	}
 }
 
+// Granger--Scott cyclotomic square
+func (e Ext12) CyclotomicSquare(x *E12) *E12 {
+	t0 := e.Ext2.Square(&x.C1.B1)
+	t1 := e.Ext2.Square(&x.C0.B0)
+	t6 := e.Ext2.Add(&x.C1.B1, &x.C0.B0)
+	t6 = e.Ext2.Square(t6)
+	t6 = e.Ext2.Sub(t6, t0)
+	t6 = e.Ext2.Sub(t6, t1)
+	t2 := e.Ext2.Square(&x.C0.B2)
+	t3 := e.Ext2.Square(&x.C1.B0)
+	t7 := e.Ext2.Add(&x.C0.B2, &x.C1.B0)
+	t7 = e.Ext2.Square(t7)
+	t7 = e.Ext2.Sub(t7, t2)
+	t7 = e.Ext2.Sub(t7, t3)
+	t4 := e.Ext2.Square(&x.C1.B2)
+	t5 := e.Ext2.Square(&x.C0.B1)
+	t8 := e.Ext2.Add(&x.C1.B2, &x.C0.B1)
+	t8 = e.Ext2.Square(t8)
+	t8 = e.Ext2.Sub(t8, t4)
+	t8 = e.Ext2.Sub(t8, t5)
+	t8 = e.Ext2.MulByNonResidue(t8)
+	t0 = e.Ext2.MulByNonResidue(t0)
+	t0 = e.Ext2.Add(t0, t1)
+	t2 = e.Ext2.MulByNonResidue(t2)
+	t2 = e.Ext2.Add(t2, t3)
+	t4 = e.Ext2.MulByNonResidue(t4)
+	t4 = e.Ext2.Add(t4, t5)
+	z00 := e.Ext2.Sub(t0, &x.C0.B0)
+	z00 = e.Ext2.Double(z00)
+	z00 = e.Ext2.Add(z00, t0)
+	z01 := e.Ext2.Sub(t2, &x.C0.B1)
+	z01 = e.Ext2.Double(z01)
+	z01 = e.Ext2.Add(z01, t2)
+	z02 := e.Ext2.Sub(t4, &x.C0.B2)
+	z02 = e.Ext2.Double(z02)
+	z02 = e.Ext2.Add(z02, t4)
+	z10 := e.Ext2.Add(t8, &x.C1.B0)
+	z10 = e.Ext2.Double(z10)
+	z10 = e.Ext2.Add(z10, t8)
+	z11 := e.Ext2.Add(t6, &x.C1.B1)
+	z11 = e.Ext2.Double(z11)
+	z11 = e.Ext2.Add(z11, t6)
+	z12 := e.Ext2.Add(t7, &x.C1.B2)
+	z12 = e.Ext2.Double(z12)
+	z12 = e.Ext2.Add(z12, t7)
+	return &E12{
+		C0: E6{
+			B0: *z00,
+			B1: *z01,
+			B2: *z02,
+		},
+		C1: E6{
+			B0: *z10,
+			B1: *z11,
+			B2: *z12,
+		},
+	}
+}
+
 func (e Ext12) AssertIsEqual(x, y *E12) {
 	e.Ext6.AssertIsEqual(&x.C0, &y.C0)
 	e.Ext6.AssertIsEqual(&x.C1, &y.C1)

std/algebra/emulated/fields_bls12381/e12_pairing.go

@@ -2,6 +2,121 @@ package fields_bls12381
 
 import "github.com/consensys/gnark/std/math/emulated"
 
+func (e Ext12) nSquareGS(z *E12, n int) *E12 {
+	for i := 0; i < n; i++ {
+		z = e.CyclotomicSquare(z)
+	}
+	return z
+}
+
+// Expt sets z to x^t in E12 and return z
+// where t = -u = 15132376222941642752
+func (e Ext12) Expt(x *E12) *E12 {
+	// Expt computation is derived from the addition chain:
+	//
+	//	_10      = 2*1
+	//	_11      = 1 + _10
+	//	_1100    = _11 << 2
+	//	_1101    = 1 + _1100
+	//	_1101000 = _1101 << 3
+	//	_1101001 = 1 + _1101000
+	//	return     ((_1101001 << 9 + 1) << 32 + 1) << 15
+	//
+	// Operations: 62 squares 5 multiplies
+	//
+	// Generated by github.com/mmcloughlin/addchain v0.4.0.
+
+	z := e.CyclotomicSquare(x)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 2)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 3)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 9)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 32)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 15)
+	z = e.CyclotomicSquare(z)
+
+	return z
+}
+
+// ExpByU sets z to x^U in E12 and return z
+// where U = (u-1)^2/3 = 76329603384216526031706109802092473003
+func (e Ext12) ExpByU(x *E12) *E12 {
+	// ExpByU computation is derived from the addition chain:
+	//
+	//	_10       = 2*1
+	//	_11       = 1 + _10
+	//	_110      = 2*_11
+	//	_111      = 1 + _110
+	//	_1000     = 1 + _111
+	//	_100000   = _1000 << 2
+	//	_100011   = _11 + _100000
+	//	_101010   = _111 + _100011
+	//	_1010100  = 2*_101010
+	//	_1010101  = 1 + _1010100
+	//	_1010110  = 1 + _1010101
+	//	_1011001  = _11 + _1010110
+	//	_1100001  = _1000 + _1011001
+	//	_10101011 = _1010101 + _1010110
+	//	_11000010 = 2*_1100001
+	//	_11100101 = _100011 + _11000010
+	//	i49       = ((_11100101 << 7 + _1011001) << 5 + _11) << 18
+	//	i68       = ((_1010101 + i49) << 9 + _10101011) << 7 + _1100001
+	//	i85       = (i68 << 9 + _10 + _10101011) << 5 + _11
+	//	i136      = ((i85 << 17 + _1010101) << 9 + _10101011) << 23
+	//	return      (_1010101 + i136) << 9 + _10101011
+	//
+	// Operations: 124 squares 23 multiplies
+	//
+	// Generated by github.com/mmcloughlin/addchain v0.4.0.
+
+	t2 := e.CyclotomicSquare(x)
+	t1 := e.Mul(x, t2)
+	z := e.CyclotomicSquare(t1)
+	z = e.Mul(x, z)
+	t3 := e.Mul(x, z)
+	t0 := e.CyclotomicSquare(t3)
+	t0 = e.CyclotomicSquare(t0)
+	t5 := e.Mul(t1, t0)
+	z = e.Mul(z, t5)
+	z = e.CyclotomicSquare(z)
+	t0 = e.Mul(x, z)
+	z = e.Mul(x, t0)
+	t4 := e.Mul(t1, z)
+	t3 = e.Mul(t3, t4)
+	z = e.Mul(t0, z)
+	t6 := e.CyclotomicSquare(t3)
+	t5 = e.Mul(t5, t6)
+	t5 = e.nSquareGS(t5, 7)
+	t4 = e.Mul(t4, t5)
+	t4 = e.nSquareGS(t4, 5)
+	t4 = e.Mul(t1, t4)
+	t4 = e.nSquareGS(t4, 18)
+	t4 = e.Mul(t0, t4)
+	t4 = e.nSquareGS(t4, 9)
+	t4 = e.Mul(z, t4)
+	t4 = e.nSquareGS(t4, 7)
+	t3 = e.Mul(t3, t4)
+	t3 = e.nSquareGS(t3, 9)
+	t2 = e.Mul(t2, t3)
+	t2 = e.Mul(z, t2)
+	t2 = e.nSquareGS(t2, 5)
+	t1 = e.Mul(t1, t2)
+	t1 = e.nSquareGS(t1, 17)
+	t1 = e.Mul(t0, t1)
+	t1 = e.nSquareGS(t1, 9)
+	t1 = e.Mul(z, t1)
+	t1 = e.nSquareGS(t1, 23)
+	t0 = e.Mul(t0, t1)
+	t0 = e.nSquareGS(t0, 9)
+	z = e.Mul(z, t0)
+
+	return z
+}
+
 func (e Ext12) nSquareTorus(z *E6, n int) *E6 {
 	for i := 0; i < n; i++ {
 		z = e.SquareTorus(z)
@@ -268,3 +383,84 @@ 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
+// 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 {
+	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)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	residueWitness := E12{
+		C0: E6{
+			B0: E2{A0: *res[0], A1: *res[1]},
+			B1: E2{A0: *res[2], A1: *res[3]},
+			B2: E2{A0: *res[4], A1: *res[5]},
+		},
+		C1: E6{
+			B0: E2{A0: *res[6], A1: *res[7]},
+			B1: E2{A0: *res[8], A1: *res[9]},
+			B2: E2{A0: *res[10], A1: *res[11]},
+		},
+	}
+
+	// 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.
+	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)
+
+	return nil
+}
+
+func (e Ext12) Frobenius(x *E12) *E12 {
+	t0 := e.Ext2.Conjugate(&x.C0.B0)
+	t1 := e.Ext2.Conjugate(&x.C0.B1)
+	t2 := e.Ext2.Conjugate(&x.C0.B2)
+	t3 := e.Ext2.Conjugate(&x.C1.B0)
+	t4 := e.Ext2.Conjugate(&x.C1.B1)
+	t5 := e.Ext2.Conjugate(&x.C1.B2)
+	t1 = e.Ext2.MulByNonResidue1Power2(t1)
+	t2 = e.Ext2.MulByNonResidue1Power4(t2)
+	t3 = e.Ext2.MulByNonResidue1Power1(t3)
+	t4 = e.Ext2.MulByNonResidue1Power3(t4)
+	t5 = e.Ext2.MulByNonResidue1Power5(t5)
+	return &E12{
+		C0: E6{
+			B0: *t0,
+			B1: *t1,
+			B2: *t2,
+		},
+		C1: E6{
+			B0: *t3,
+			B1: *t4,
+			B2: *t5,
+		},
+	}
+}
+
+func (e Ext12) FrobeniusSquare(x *E12) *E12 {
+	z00 := &x.C0.B0
+	z01 := e.Ext2.MulByNonResidue2Power2(&x.C0.B1)
+	z02 := e.Ext2.MulByNonResidue2Power4(&x.C0.B2)
+	z10 := e.Ext2.MulByNonResidue2Power1(&x.C1.B0)
+	z11 := e.Ext2.MulByNonResidue2Power3(&x.C1.B1)
+	z12 := e.Ext2.MulByNonResidue2Power5(&x.C1.B2)
+	return &E12{
+		C0: E6{B0: *z00, B1: *z01, B2: *z02},
+		C1: E6{B0: *z10, B1: *z11, B2: *z12},
+	}
+}

std/algebra/emulated/fields_bls12381/hints.go

@@ -27,6 +27,7 @@ func GetHints() []solver.Hint {
 		// E12
 		divE12Hint,
 		inverseE12Hint,
+		finalExpHint,
 	}
 }
 
@@ -268,3 +269,46 @@ func divE12Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) erro
 			return nil
 		})
 }
+
+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)
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var millerLoop, residueWitness bls12381.E12
+			var rInv big.Int
+
+			millerLoop.C0.B0.A0.SetBigInt(inputs[0])
+			millerLoop.C0.B0.A1.SetBigInt(inputs[1])
+			millerLoop.C0.B1.A0.SetBigInt(inputs[2])
+			millerLoop.C0.B1.A1.SetBigInt(inputs[3])
+			millerLoop.C0.B2.A0.SetBigInt(inputs[4])
+			millerLoop.C0.B2.A1.SetBigInt(inputs[5])
+			millerLoop.C1.B0.A0.SetBigInt(inputs[6])
+			millerLoop.C1.B0.A1.SetBigInt(inputs[7])
+			millerLoop.C1.B1.A0.SetBigInt(inputs[8])
+			millerLoop.C1.B1.A1.SetBigInt(inputs[9])
+			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)
+
+			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 nil
+		})
+}

std/algebra/emulated/sw_bls12381/pairing.go

@@ -245,13 +245,22 @@ 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.Pair(P, Q)
+	f, err := pr.MillerLoop(P, Q)
 	if err != nil {
 		return err
 
 	}
-	one := pr.One()
-	pr.AssertIsEqual(f, one)
+	// We perform the easy part of the final exp to push f to the cyclotomic
+	// subgroup so that FinalExponentiationCheck is carried with optimized
+	// cyclotomic squaring (e.g. Karabina12345).
+	//
+	// f = f^(p⁶-1)(p²+1)
+	buf := pr.Conjugate(f)
+	buf = pr.DivUnchecked(buf, f)
+	f = pr.FrobeniusSquare(buf)
+	f = pr.Mul(f, buf)
+
+	pr.FinalExponentiationCheck(f)
 
 	return nil
 }