perf: replace BN254 final exp by a class equivalence check

go.mod

@@ -9,7 +9,7 @@ require (
 	github.com/blang/semver/v4 v4.0.0
 	github.com/consensys/bavard v0.1.13
 	github.com/consensys/compress v0.2.5
-	github.com/consensys/gnark-crypto v0.12.2-0.20240504013751-564b6f724c3b
+	github.com/consensys/gnark-crypto v0.12.2-0.20240703135258-5d8b5fab1afb
 	github.com/fxamacker/cbor/v2 v2.5.0
 	github.com/google/go-cmp v0.5.9
 	github.com/google/pprof v0.0.0-20230817174616-7a8ec2ada47b

go.sum

@@ -6,8 +6,8 @@ github.com/consensys/bavard v0.1.13 h1:oLhMLOFGTLdlda/kma4VOJazblc7IM5y5QPd2A/Yj
 github.com/consensys/bavard v0.1.13/go.mod h1:9ItSMtA/dXMAiL7BG6bqW2m3NdSEObYWoH223nGHukI=
 github.com/consensys/compress v0.2.5 h1:gJr1hKzbOD36JFsF1AN8lfXz1yevnJi1YolffY19Ntk=
 github.com/consensys/compress v0.2.5/go.mod h1:pyM+ZXiNUh7/0+AUjUf9RKUM6vSH7T/fsn5LLS0j1Tk=
-github.com/consensys/gnark-crypto v0.12.2-0.20240504013751-564b6f724c3b h1:tu0NaVr64o6vXzy9rYSK/LCZXmS+u/k9eP1F8OtRUWQ=
-github.com/consensys/gnark-crypto v0.12.2-0.20240504013751-564b6f724c3b/go.mod h1:wKqwsieaKPThcFkHe0d0zMsbHEUWFmZcG7KBCse210o=
+github.com/consensys/gnark-crypto v0.12.2-0.20240703135258-5d8b5fab1afb h1:LMfC1GSeYv1TKp3zNIuwYNEqb9RCVrMkCV4Y9k5ZJ6o=
+github.com/consensys/gnark-crypto v0.12.2-0.20240703135258-5d8b5fab1afb/go.mod h1:wKqwsieaKPThcFkHe0d0zMsbHEUWFmZcG7KBCse210o=
 github.com/coreos/go-systemd/v22 v22.5.0/go.mod h1:Y58oyj3AT4RCenI/lSvhwexgC+NSVTIJ3seZv2GcEnc=
 github.com/davecgh/go-spew v1.1.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c=
 github.com/davecgh/go-spew v1.1.1/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=

std/algebra/emulated/fields_bn254/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_bn254/e12_pairing.go

@@ -4,6 +4,76 @@ 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
+}
+
+// Exponentiation by U=6u+2 where t is the seed u=4965661367192848881
+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
+	//	_1100    = 2*_110
+	//	_1111    = _11 + _1100
+	//	_1100000 = _1100 << 3
+	//	_1100111 = _111 + _1100000
+	//	i22      = ((_1100111 << 2 + 1) << 5 + _1111) << 3
+	//	i38      = ((1 + i22) << 4 + _111) << 9 + _111
+	//	i50      = 2*((i38 << 4 + _11) << 5 + _1111)
+	//	i61      = ((1 + i50) << 5 + _111) << 3 + _11
+	//	i75      = ((i61 << 6 + _111) << 4 + _111) << 2
+	//	return     ((1 + i75) << 2 + 1) << 3
+	//
+	// Operations: 64 squares 18 multiplies
+	//
+	// Generated by github.com/mmcloughlin/addchain v0.4.0.
+
+	z := e.Square(x)
+	t0 := e.Mul(x, z)
+	t1 := e.Square(t0)
+	z = e.Mul(x, t1)
+	t2 := e.Square(t1)
+	t1 = e.Mul(t0, t2)
+	t2 = e.nSquareGS(t2, 3)
+	t2 = e.Mul(z, t2)
+	t2 = e.nSquareGS(t2, 2)
+	t2 = e.Mul(x, t2)
+	t2 = e.nSquareGS(t2, 5)
+	t2 = e.Mul(t1, t2)
+	t2 = e.nSquareGS(t2, 3)
+	t2 = e.Mul(x, t2)
+	t2 = e.nSquareGS(t2, 4)
+	t2 = e.Mul(z, t2)
+	t2 = e.nSquareGS(t2, 9)
+	t2 = e.Mul(z, t2)
+	t2 = e.nSquareGS(t2, 4)
+	t2 = e.Mul(t0, t2)
+	t2 = e.nSquareGS(t2, 5)
+	t1 = e.Mul(t1, t2)
+	t1 = e.Square(t1)
+	t1 = e.Mul(x, t1)
+	t1 = e.nSquareGS(t1, 5)
+	t1 = e.Mul(z, t1)
+	t1 = e.nSquareGS(t1, 3)
+	t0 = e.Mul(t0, t1)
+	t0 = e.nSquareGS(t0, 6)
+	t0 = e.Mul(z, t0)
+	t0 = e.nSquareGS(t0, 4)
+	z = e.Mul(z, t0)
+	z = e.nSquareGS(z, 2)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 2)
+	z = e.Mul(x, z)
+	z = e.nSquareGS(z, 3)
+
+	return z
+}
+
 func (e Ext12) nSquareTorus(z *E6, n int) *E6 {
 	for i := 0; i < n; i++ {
 		z = e.SquareTorus(z)
@@ -14,7 +84,7 @@ func (e Ext12) nSquareTorus(z *E6, n int) *E6 {
 // Exponentiation by the seed t=4965661367192848881
 // The computations are performed on E6 compressed form using Torus-based arithmetic.
 func (e Ext12) ExptTorus(x *E6) *E6 {
-	// Expt computation is derived from the addition chain:
+	// ExptTorus computation is derived from the addition chain:
 	//
 	//	_10     = 2*1
 	//	_100    = 2*_10
@@ -351,3 +421,127 @@ func (e Ext12) FrobeniusCubeTorus(y *E6) *E6 {
 
 	return res
 }
+
+// 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, 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
+		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]},
+		},
+	}
+	cubicNonResiduePower := 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: E6{
+			B0: E2{A0: *res[18], A1: *res[19]},
+			B1: E2{A0: *res[20], A1: *res[21]},
+			B2: E2{A0: *res[22], A1: *res[23]},
+		},
+	}
+
+	// Check that  x * cubicNonResiduePower == residueWitness^λ
+	// where λ = 6u + 2 + q^3 - q^2 + q, with u the BN254 seed
+	// and residueWitness, cubicNonResiduePower from the hint.
+	t2 := e.Mul(&cubicNonResiduePower, x)
+
+	t1 := e.FrobeniusCube(&residueWitness)
+	t0 := e.FrobeniusSquare(&residueWitness)
+	t1 = e.DivUnchecked(t1, t0)
+	t0 = e.Frobenius(&residueWitness)
+	t1 = e.Mul(t1, t0)
+
+	// exponentiation by U=6u+2
+	t0 = e.ExpByU(&residueWitness)
+
+	t0 = e.Mul(t0, t1)
+
+	e.AssertIsEqual(t0, t2)
+
+	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},
+	}
+}
+
+func (e Ext12) FrobeniusCube(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.MulByNonResidue3Power2(t1)
+	t2 = e.Ext2.MulByNonResidue3Power4(t2)
+	t3 = e.Ext2.MulByNonResidue3Power1(t3)
+	t4 = e.Ext2.MulByNonResidue3Power3(t4)
+	t5 = e.Ext2.MulByNonResidue3Power5(t5)
+	return &E12{
+		C0: E6{
+			B0: *t0,
+			B1: *t1,
+			B2: *t2,
+		},
+		C1: E6{
+			B0: *t3,
+			B1: *t4,
+			B2: *t5,
+		},
+	}
+}