diff --git a/std/algebra/emulated/fields_bw6761/doc.go b/std/algebra/emulated/fields_bw6761/doc.go
index 29858f7424..3779fdb2dc 100644
--- a/std/algebra/emulated/fields_bw6761/doc.go
+++ b/std/algebra/emulated/fields_bw6761/doc.go
@@ -2,5 +2,5 @@
 // used to compute the pairing over the BW6-761 curve.
 //
 //	𝔽p³[u] = 𝔽p/u³+4
-//	𝔽p⁶[v] = 𝔽p²/v²-u
+//	𝔽p⁶[v] = 𝔽p³/v²-u
 package fields_bw6761
diff --git a/std/algebra/emulated/sw_bls12381/hints.go b/std/algebra/emulated/sw_bls12381/hints.go
new file mode 100644
index 0000000000..3280da400a
--- /dev/null
+++ b/std/algebra/emulated/sw_bls12381/hints.go
@@ -0,0 +1,142 @@
+package sw_bls12381
+
+import (
+	"math/big"
+
+	"github.com/consensys/gnark-crypto/ecc/bls12-381"
+	"github.com/consensys/gnark/constraint/solver"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+func init() {
+	solver.RegisterHint(GetHints()...)
+}
+
+// GetHints returns all hint functions used in the package.
+func GetHints() []solver.Hint {
+	return []solver.Hint{
+		doublePairingCheckHint,
+	}
+}
+
+func doublePairingCheckHint(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 P0, P1 bls12381.G1Affine
+			var Q0, Q1 bls12381.G2Affine
+
+			P0.X.SetBigInt(inputs[0])
+			P0.Y.SetBigInt(inputs[1])
+			P1.X.SetBigInt(inputs[2])
+			P1.Y.SetBigInt(inputs[3])
+			Q0.X.A0.SetBigInt(inputs[4])
+			Q0.X.A1.SetBigInt(inputs[5])
+			Q0.Y.A0.SetBigInt(inputs[6])
+			Q0.Y.A1.SetBigInt(inputs[7])
+			Q1.X.A0.SetBigInt(inputs[8])
+			Q1.X.A1.SetBigInt(inputs[9])
+			Q1.Y.A0.SetBigInt(inputs[10])
+			Q1.Y.A1.SetBigInt(inputs[11])
+
+			lines0 := bls12381.PrecomputeLines(Q0)
+			lines1 := bls12381.PrecomputeLines(Q1)
+			millerLoop, err := bls12381.MillerLoopFixedQ(
+				[]bls12381.G1Affine{P0, P1},
+				[][2][len(bls12381.LoopCounter) - 1]bls12381.LineEvaluationAff{lines0, lines1},
+			)
+			millerLoop.Conjugate(&millerLoop)
+			if err != nil {
+				return err
+			}
+
+			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])
+			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
+		})
+}
diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go
index 4bc99671d6..b59d18c8b4 100644
--- a/std/algebra/emulated/sw_bls12381/pairing.go
+++ b/std/algebra/emulated/sw_bls12381/pairing.go
@@ -265,6 +265,121 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck calculates the reduced pairing for a 2 pairs of points and asserts if the result is One
+// e(P0, Q0) * e(P1, Q1) =? 1
+//
+// This function doesn't check that the inputs are in the correct subgroups. See AssertIsOnG1 and AssertIsOnG2.
+func (pr Pairing) DoublePairingCheck(P [2]*G1Affine, Q [2]*G2Affine) error {
+	// hint the non-residue witness
+	hint, err := pr.curveF.NewHint(doublePairingCheckHint, 18, &P[0].X, &P[0].Y, &P[1].X, &P[1].Y, &Q[0].P.X.A0, &Q[0].P.X.A1, &Q[0].P.Y.A0, &Q[0].P.Y.A1, &Q[1].P.X.A0, &Q[1].P.X.A1, &Q[1].P.Y.A0, &Q[1].P.Y.A1)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	residueWitness := fields_bls12381.E12{
+		C0: fields_bls12381.E6{
+			B0: fields_bls12381.E2{A0: *hint[0], A1: *hint[1]},
+			B1: fields_bls12381.E2{A0: *hint[2], A1: *hint[3]},
+			B2: fields_bls12381.E2{A0: *hint[4], A1: *hint[5]},
+		},
+		C1: fields_bls12381.E6{
+			B0: fields_bls12381.E2{A0: *hint[6], A1: *hint[7]},
+			B1: fields_bls12381.E2{A0: *hint[8], A1: *hint[9]},
+			B2: fields_bls12381.E2{A0: *hint[10], A1: *hint[11]},
+		},
+	}
+	// constrain scaling factor to be in Fp6
+	scalingFactor := fields_bls12381.E12{
+		C0: fields_bls12381.E6{
+			B0: fields_bls12381.E2{A0: *hint[12], A1: *hint[13]},
+			B1: fields_bls12381.E2{A0: *hint[14], A1: *hint[15]},
+			B2: fields_bls12381.E2{A0: *hint[16], A1: *hint[17]},
+		},
+		C1: (*pr.Ext6.Zero()),
+	}
+
+	// residueWitnessInv = 1 / residueWitness
+	residueWitnessInv := pr.Inverse(&residueWitness)
+
+	if Q[0].Lines == nil {
+		Q0lines := pr.computeLines(&Q[0].P)
+		Q[0].Lines = &Q0lines
+	}
+	lines0 := *Q[0].Lines
+	if Q[1].Lines == nil {
+		Q1lines := pr.computeLines(&Q[1].P)
+		Q[1].Lines = &Q1lines
+	}
+	lines1 := *Q[1].Lines
+
+	// precomputations
+	y0Inv := pr.curveF.Inverse(&P[0].Y)
+	x0NegOverY0 := pr.curveF.Mul(&P[0].X, y0Inv)
+	x0NegOverY0 = pr.curveF.Neg(x0NegOverY0)
+	y1Inv := pr.curveF.Inverse(&P[1].Y)
+	x1NegOverY1 := pr.curveF.Mul(&P[1].X, y1Inv)
+	x1NegOverY1 = pr.curveF.Neg(x1NegOverY1)
+
+	// init Miller loop accumulator to residueWitnessInv to share the squarings
+	// of residueWitnessInv^{-x₀}
+	res := residueWitnessInv
+
+	// Compute f_{x₀,Q}(P)
+	for i := 62; i >= 0; i-- {
+		res = pr.Square(res)
+
+		if loopCounter[i] == 0 {
+			// ℓ × res
+			res = pr.MulBy014(res,
+				pr.MulByElement(&lines0[0][i].R1, y0Inv),
+				pr.MulByElement(&lines0[0][i].R0, x0NegOverY0),
+			)
+
+			// ℓ × res
+			res = pr.MulBy014(res,
+				pr.MulByElement(&lines1[0][i].R1, y1Inv),
+				pr.MulByElement(&lines1[0][i].R0, x1NegOverY1),
+			)
+		} else {
+			// multiply by residueWitnessInv when bit=1
+			res = pr.Mul(res, residueWitnessInv)
+
+			res = pr.MulBy014(res,
+				pr.MulByElement(&lines0[0][i].R1, y0Inv),
+				pr.MulByElement(&lines0[0][i].R0, x0NegOverY0),
+			)
+			res = pr.MulBy014(res,
+				pr.MulByElement(&lines0[1][i].R1, y0Inv),
+				pr.MulByElement(&lines0[1][i].R0, x0NegOverY0),
+			)
+
+			res = pr.MulBy014(res,
+				pr.MulByElement(&lines1[0][i].R1, y1Inv),
+				pr.MulByElement(&lines1[0][i].R0, x1NegOverY1),
+			)
+			res = pr.MulBy014(res,
+				pr.MulByElement(&lines1[1][i].R1, y1Inv),
+				pr.MulByElement(&lines1[1][i].R0, x1NegOverY1),
+			)
+		}
+	}
+
+	// Check that  res * scalingFactor * residueWitnessInv^λ' == 1
+	// where λ' = q, with u the BLS12-381 seed
+	// and residueWitnessInv, 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(residueWitnessInv)
+	t0 = pr.Mul(t0, res)
+	t0 = pr.Mul(t0, &scalingFactor)
+
+	pr.AssertIsEqual(t0, pr.One())
+
+	return nil
+
+}
+
 func (pr Pairing) AssertIsEqual(x, y *GTEl) {
 	pr.Ext12.AssertIsEqual(x, y)
 }
diff --git a/std/algebra/emulated/sw_bls12381/pairing_test.go b/std/algebra/emulated/sw_bls12381/pairing_test.go
index 9ffdd18d56..2eaf6d08fb 100644
--- a/std/algebra/emulated/sw_bls12381/pairing_test.go
+++ b/std/algebra/emulated/sw_bls12381/pairing_test.go
@@ -165,26 +165,26 @@ func TestMultiPairTestSolve(t *testing.T) {
 	}
 }
 
-type PairingCheckCircuit struct {
+type DoublePairingCheckCircuit struct {
 	In1G1 G1Affine
 	In2G1 G1Affine
 	In1G2 G2Affine
 	In2G2 G2Affine
 }
 
-func (c *PairingCheckCircuit) Define(api frontend.API) error {
+func (c *DoublePairingCheckCircuit) Define(api frontend.API) error {
 	pairing, err := NewPairing(api)
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
-	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2})
+	err = pairing.DoublePairingCheck([2]*G1Affine{&c.In1G1, &c.In2G1}, [2]*G2Affine{&c.In1G2, &c.In2G2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
 	}
 	return nil
 }
 
-func TestPairingCheckTestSolve(t *testing.T) {
+func TestDoublePairingCheckTestSolve(t *testing.T) {
 	assert := test.NewAssert(t)
 	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
@@ -193,6 +193,41 @@ func TestPairingCheckTestSolve(t *testing.T) {
 	var q2 bls12381.G2Affine
 	q2.Set(&q1)
 	q1.Double(&q1)
+	witness := DoublePairingCheckCircuit{
+		In1G1: NewG1Affine(p1),
+		In1G2: NewG2Affine(q1),
+		In2G1: NewG1Affine(p2),
+		In2G2: NewG2Affine(q2),
+	}
+	err := test.IsSolved(&DoublePairingCheckCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type PairingCheckCircuit struct {
+	In1G1 G1Affine
+	In2G1 G1Affine
+	In1G2 G2Affine
+	In2G2 G2Affine
+}
+
+func (c *PairingCheckCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	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})
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	return nil
+}
+
+func TestPairingCheckTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p1, q1 := randomG1G2Affines()
+	_, q2 := randomG1G2Affines()
+	var p2 bls12381.G1Affine
+	p2.Neg(&p1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),
diff --git a/std/algebra/emulated/sw_bn254/hints.go b/std/algebra/emulated/sw_bn254/hints.go
index 3c4d5642fe..8b918eb7ca 100644
--- a/std/algebra/emulated/sw_bn254/hints.go
+++ b/std/algebra/emulated/sw_bn254/hints.go
@@ -17,6 +17,7 @@ func init() {
 func GetHints() []solver.Hint {
 	return []solver.Hint{
 		millerLoopAndCheckFinalExpHint,
+		doublePairingCheckHint,
 	}
 }
 
@@ -135,7 +136,135 @@ func millerLoopAndCheckFinalExpHint(nativeMod *big.Int, nativeInputs, nativeOutp
 
 			// x is now the cube root of residueWitness
 			residueWitness.Set(&x)
+
+			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])
+
+			// we also need to return the cubic non-residue power
+			cubicNonResiduePower.C0.B0.A0.BigInt(outputs[12])
+			cubicNonResiduePower.C0.B0.A1.BigInt(outputs[13])
+			cubicNonResiduePower.C0.B1.A0.BigInt(outputs[14])
+			cubicNonResiduePower.C0.B1.A1.BigInt(outputs[15])
+			cubicNonResiduePower.C0.B2.A0.BigInt(outputs[16])
+			cubicNonResiduePower.C0.B2.A1.BigInt(outputs[17])
+
+			return nil
+		})
+}
+
+func doublePairingCheckHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error {
+	// This follows section 4.3.2 of https://eprint.iacr.org/2024/640.pdf
+	return emulated.UnwrapHint(nativeInputs, nativeOutputs,
+		func(mod *big.Int, inputs, outputs []*big.Int) error {
+			var tmp, x3, cubicNonResiduePower, x, millerLoop, residueWitness, residueWitnessInv, one, root27thOf1 bn254.E12
+			var exp1, exp2, rInv, mInv big.Int
+			var P0, P1 bn254.G1Affine
+			var Q0, Q1 bn254.G2Affine
+
+			P0.X.SetBigInt(inputs[0])
+			P0.Y.SetBigInt(inputs[1])
+			P1.X.SetBigInt(inputs[2])
+			P1.Y.SetBigInt(inputs[3])
+			Q0.X.A0.SetBigInt(inputs[4])
+			Q0.X.A1.SetBigInt(inputs[5])
+			Q0.Y.A0.SetBigInt(inputs[6])
+			Q0.Y.A1.SetBigInt(inputs[7])
+			Q1.X.A0.SetBigInt(inputs[8])
+			Q1.X.A1.SetBigInt(inputs[9])
+			Q1.Y.A0.SetBigInt(inputs[10])
+			Q1.Y.A1.SetBigInt(inputs[11])
+
+			lines0 := bn254.PrecomputeLines(Q0)
+			lines1 := bn254.PrecomputeLines(Q1)
+			millerLoop, err := bn254.MillerLoopFixedQ(
+				[]bn254.G1Affine{P0, P1},
+				[][2][len(bn254.LoopCounter)]bn254.LineEvaluationAff{lines0, lines1},
+			)
+			if err != nil {
+				return err
+			}
+
+			// exp1 = (p^12-1)/3
+			exp1.SetString("4030969696062745741797811005853058291874379204406359442560681893891674450106959530046539719647151210908190211459382793062006703141168852426020468083171325367934590379984666859998399967609544754664110191464072930598755441160008826659219834762354786403012110463250131961575955268597858015384895449311534622125256548620283853223733396368939858981844663598065852816056384933498610930035891058807598891752166582271931875150099691598048016175399382213304673796601585080509443902692818733420199004555566113537482054218823936116647313678747500267068559627206777530424029211671772692598157901876223857571299238046741502089890557442500582300718504160740314926185458079985126192563953772118929726791041828902047546977272656240744693339962973939047279285351052107950250121751682659529260304162131862468322644288196213423232132152125277136333208005221619443705106431645884840489295409272576227859206166894626854018093044908314720", 10)
+			// root27thOf1 = (0, c010, c011, 0, 0, 0, 0, 0, 0, 0, 0, 0)
+			// is a 27-th root of  unity which is necessarily a cubic non-residue
+			// since h/r = (p^12-1)/r = 27·l and 3 does not divide l.
+			// it was computed as w^((p^12-1)/27) = c2 * w^2 + c8 * w^8 where
+			// Fp12 = Fp[w]/w^12-18w^6+82 which is isomorphic to our Fp12 tower
+			// then c010 = (c2 + 9 * c8) % p and c011 = c8
+			root27thOf1.C0.B1.A0.SetString("9483667112135124394372960210728142145589475128897916459350428495526310884707")
+			root27thOf1.C0.B1.A1.SetString("4534159768373982659291990808346042891252278737770656686799127720849666919525")
+
+			if one.Exp(millerLoop, &exp1).IsOne() {
+				// residueWitness = millerLoop is a cubic residue
+				cubicNonResiduePower.SetOne()
+				residueWitness.Set(&millerLoop)
+			} else if one.Exp(*millerLoop.Mul(&millerLoop, &root27thOf1), &exp1).IsOne() {
+				// residueWitness = millerLoop * root27thOf1 is a cubic residue
+				cubicNonResiduePower.Set(&root27thOf1)
+				residueWitness.Set(&millerLoop)
+			} else {
+				// residueWitness = millerLoop * root27thOf1^2 is a cubic residue
+				cubicNonResiduePower.Square(&root27thOf1)
+				residueWitness.Mul(&millerLoop, &root27thOf1)
+			}
+
+			// 1. compute r-th root:
+			// Exponentiate to rInv where
+			// rInv = 1/r mod (p^12-1)/r
+			rInv.SetString("495819184011867778744231927046742333492451180917315223017345540833046880485481720031136878341141903241966521818658471092566752321606779256340158678675679238405722886654128392203338228575623261160538734808887996935946888297414610216445334190959815200956855428635568184508263913274453942864817234480763055154719338281461936129150171789463489422401982681230261920147923652438266934726901346095892093443898852488218812468761027620988447655860644584419583586883569984588067403598284748297179498734419889699245081714359110559679136004228878808158639412436468707589339209058958785568729925402190575720856279605832146553573981587948304340677613460685405477047119496887534881410757668344088436651291444274840864486870663164657544390995506448087189408281061890434467956047582679858345583941396130713046072603335601764495918026585155498301896749919393", 10)
+			residueWitness.Exp(residueWitness, &rInv)
+
+			// 2. compute m-th root:
+			// where m = (6x + 2 + q^3 - q^2 + q)/(3r)
+			// Exponentiate to mInv where
+			// mInv = 1/m mod p^12-1
+			mInv.SetString("17840267520054779749190587238017784600702972825655245554504342129614427201836516118803396948809179149954197175783449826546445899524065131269177708416982407215963288737761615699967145070776364294542559324079147363363059480104341231360692143673915822421222230661528586799190306058519400019024762424366780736540525310403098758015600523609594113357130678138304964034267260758692953579514899054295817541844330584721967571697039986079722203518034173581264955381924826388858518077894154909963532054519350571947910625755075099598588672669612434444513251495355121627496067454526862754597351094345783576387352673894873931328099247263766690688395096280633426669535619271711975898132416216382905928886703963310231865346128293216316379527200971959980873989485521004596686352787540034457467115536116148612884807380187255514888720048664139404687086409399", 10)
+			residueWitness.Exp(residueWitness, &mInv)
+
+			// 3. compute cube root:
+			// since gcd(3, (p^12-1)/r) ≠ 1 we use a modified Toneelli-Shanks algorithm
+			// see Alg.4 of https://eprint.iacr.org/2024/640.pdf
+			// Typo in the paper: p^k-1 = 3^n * s instead of p-1 = 3^r * s
+			// where k=12 and n=3 here and exp2 = (s+1)/3
 			residueWitnessInv.Inverse(&residueWitness)
+			exp2.SetString("149295173928249842288807815031594751550902933496531831205951181255247201855813315927649619246190785589192230054051214557852100116339587126889646966043382421034614458517950624444385183985538694617189266350521219651805757080000326913304438324531658755667115202342597480058368713651772519088329461085612393412046538837788290860138273939590365147475728281409846400594680923462911515927255224400281440435265428973034513894448136725853630228718495637529802733207466114092942366766400693830377740909465411612499335341437923559875826432546203713595131838044695464089778859691547136762894737106526809539677749557286722299625576201574095640767352005953344997266128077036486155280146436004404804695964512181557316554713802082990544197776406442186936269827816744738898152657469728130713344598597476387715653492155415311971560450078713968012341037230430349766855793764662401499603533676762082513303932107208402000670112774382027", 10)
+			x.Exp(residueWitness, &exp2)
+
+			// 3^t is ord(x^3 / residueWitness)
+			x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv)
+			t := 0
+			for !x3.IsOne() {
+				t++
+				tmp.Square(&x3)
+				x3.Mul(&tmp, &x3)
+			}
+
+			for t != 0 {
+				x.Mul(&x, tmp.Exp(root27thOf1, &exp2))
+
+				// 3^t is ord(x^3 / residueWitness)
+				x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv)
+				t = 0
+				for !x3.IsOne() {
+					t++
+					tmp.Square(&x3)
+					x3.Mul(&tmp, &x3)
+				}
+			}
+
+			// x is now the cube root of residueWitness
+			residueWitness.Set(&x)
 
 			residueWitness.C0.B0.A0.BigInt(outputs[0])
 			residueWitness.C0.B0.A1.BigInt(outputs[1])
diff --git a/std/algebra/emulated/sw_bn254/pairing.go b/std/algebra/emulated/sw_bn254/pairing.go
index 61a7e051e6..d0b1c6723c 100644
--- a/std/algebra/emulated/sw_bn254/pairing.go
+++ b/std/algebra/emulated/sw_bn254/pairing.go
@@ -265,6 +265,190 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck calculates the reduced pairing for a 2 pairs of points and asserts if the result is One
+// e(P0, Q0) * e(P1, Q1) =? 1
+//
+// This function doesn't check that the inputs are in the correct subgroups. See AssertIsOnG1 and AssertIsOnG2.
+func (pr Pairing) DoublePairingCheck(P [2]*G1Affine, Q [2]*G2Affine) error {
+	// hint the non-residue witness
+	hint, err := pr.curveF.NewHint(doublePairingCheckHint, 18, &P[0].X, &P[0].Y, &P[1].X, &P[1].Y, &Q[0].P.X.A0, &Q[0].P.X.A1, &Q[0].P.Y.A0, &Q[0].P.Y.A1, &Q[1].P.X.A0, &Q[1].P.X.A1, &Q[1].P.Y.A0, &Q[1].P.Y.A1)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	residueWitness := fields_bn254.E12{
+		C0: fields_bn254.E6{
+			B0: fields_bn254.E2{A0: *hint[0], A1: *hint[1]},
+			B1: fields_bn254.E2{A0: *hint[2], A1: *hint[3]},
+			B2: fields_bn254.E2{A0: *hint[4], A1: *hint[5]},
+		},
+		C1: fields_bn254.E6{
+			B0: fields_bn254.E2{A0: *hint[6], A1: *hint[7]},
+			B1: fields_bn254.E2{A0: *hint[8], A1: *hint[9]},
+			B2: fields_bn254.E2{A0: *hint[10], A1: *hint[11]},
+		},
+	}
+	// constrain cubicNonResiduePower to be in Fp6
+	cubicNonResiduePower := fields_bn254.E12{
+		C0: fields_bn254.E6{
+			B0: fields_bn254.E2{A0: *hint[12], A1: *hint[13]},
+			B1: fields_bn254.E2{A0: *hint[14], A1: *hint[15]},
+			B2: fields_bn254.E2{A0: *hint[16], A1: *hint[17]},
+		},
+		C1: (*pr.Ext6.Zero()),
+	}
+
+	// residueWitnessInv = 1 / residueWitness
+	residueWitnessInv := pr.Inverse(&residueWitness)
+
+	if Q[0].Lines == nil {
+		Q0lines := pr.computeLines(&Q[0].P)
+		Q[0].Lines = &Q0lines
+	}
+	lines0 := *Q[0].Lines
+	if Q[1].Lines == nil {
+		Q1lines := pr.computeLines(&Q[1].P)
+		Q[1].Lines = &Q1lines
+	}
+	lines1 := *Q[1].Lines
+
+	// precomputations
+	y0Inv := pr.curveF.Inverse(&P[0].Y)
+	x0NegOverY0 := pr.curveF.Mul(&P[0].X, y0Inv)
+	x0NegOverY0 = pr.curveF.Neg(x0NegOverY0)
+	y1Inv := pr.curveF.Inverse(&P[1].Y)
+	x1NegOverY1 := pr.curveF.Mul(&P[1].X, y1Inv)
+	x1NegOverY1 = pr.curveF.Neg(x1NegOverY1)
+
+	// init Miller loop accumulator to residueWitnessInv to share the squarings
+	// of residueWitnessInv^{6x₀+2}
+	res := residueWitnessInv
+
+	// Compute f_{6x₀+2,Q}(P)
+	for i := 64; i >= 0; i-- {
+		res = pr.Square(res)
+
+		switch loopCounter[i] {
+		case 0:
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines0[0][i].R0, x0NegOverY0),
+				pr.MulByElement(&lines0[0][i].R1, y0Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines1[0][i].R0, x1NegOverY1),
+				pr.MulByElement(&lines1[0][i].R1, y1Inv),
+			)
+		case 1:
+			// multiply by residueWitnessInv when bit=1
+			res = pr.Mul(res, residueWitnessInv)
+
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines0[0][i].R0, x0NegOverY0),
+				pr.MulByElement(&lines0[0][i].R1, y0Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines0[1][i].R0, x0NegOverY0),
+				pr.MulByElement(&lines0[1][i].R1, y0Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines1[0][i].R0, x1NegOverY1),
+				pr.MulByElement(&lines1[0][i].R1, y1Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines1[1][i].R0, x1NegOverY1),
+				pr.MulByElement(&lines1[1][i].R1, y1Inv),
+			)
+		case -1:
+			// multiply by residueWitness when bit=-1
+			res = pr.Mul(res, &residueWitness)
+
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines0[0][i].R0, x0NegOverY0),
+				pr.MulByElement(&lines0[0][i].R1, y0Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines0[1][i].R0, x0NegOverY0),
+				pr.MulByElement(&lines0[1][i].R1, y0Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines1[0][i].R0, x1NegOverY1),
+				pr.MulByElement(&lines1[0][i].R1, y1Inv),
+			)
+			// ℓ × res
+			res = pr.MulBy034(
+				res,
+				pr.MulByElement(&lines1[1][i].R0, x1NegOverY1),
+				pr.MulByElement(&lines1[1][i].R1, y1Inv),
+			)
+		default:
+			panic(fmt.Sprintf("invalid loop counter value %d", loopCounter[i]))
+		}
+	}
+
+	// Compute  ℓ_{[6x₀+2]Q,π(Q)}(P) · ℓ_{[6x₀+2]Q+π(Q),-π²(Q)}(P)
+	// ℓ × res
+	res = pr.MulBy034(
+		res,
+		pr.MulByElement(&lines0[0][65].R0, x0NegOverY0),
+		pr.MulByElement(&lines0[0][65].R1, y0Inv),
+	)
+	// ℓ × res
+	res = pr.MulBy034(
+		res,
+		pr.MulByElement(&lines0[1][65].R0, x0NegOverY0),
+		pr.MulByElement(&lines0[1][65].R1, y0Inv),
+	)
+	res = pr.MulBy034(
+		res,
+		pr.MulByElement(&lines1[0][65].R0, x1NegOverY1),
+		pr.MulByElement(&lines1[0][65].R1, y1Inv),
+	)
+	// ℓ × res
+	res = pr.MulBy034(
+		res,
+		pr.MulByElement(&lines1[1][65].R0, x1NegOverY1),
+		pr.MulByElement(&lines1[1][65].R1, y1Inv),
+	)
+
+	// Check that  res * cubicNonResiduePower * residueWitnessInv^λ' == 1
+	// where λ' = q^3 - q^2 + q, with u the BN254 seed
+	// and residueWitnessInv, cubicNonResiduePower from the hint.
+	// Note that res is already MillerLoop(P,Q) * residueWitnessInv^{6x₀+2} since
+	// we initialized the Miller loop accumulator with residueWitnessInv.
+	t2 := pr.Mul(&cubicNonResiduePower, res)
+
+	t1 := pr.FrobeniusCube(residueWitnessInv)
+	t0 := pr.FrobeniusSquare(residueWitnessInv)
+	t1 = pr.DivUnchecked(t1, t0)
+	t0 = pr.Frobenius(residueWitnessInv)
+	t1 = pr.Mul(t1, t0)
+
+	t2 = pr.Mul(t2, t1)
+	pr.AssertIsEqual(t2, pr.One())
+
+	return nil
+
+}
+
 func (pr Pairing) IsEqual(x, y *GTEl) frontend.Variable {
 	return pr.Ext12.IsEqual(x, y)
 }
@@ -401,7 +585,7 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl
 		// 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.MulMod(&P[k].X, yInv[k])
+		xNegOverY[k] = pr.curveF.Mul(&P[k].X, yInv[k])
 		xNegOverY[k] = pr.curveF.Neg(xNegOverY[k])
 	}
 
@@ -727,7 +911,7 @@ func (pr Pairing) millerLoopAndFinalExpResult(P *G1Affine, Q *G2Affine, previous
 
 	// precomputations
 	yInv := pr.curveF.Inverse(&P.Y)
-	xNegOverY := pr.curveF.MulMod(&P.X, yInv)
+	xNegOverY := pr.curveF.Mul(&P.X, yInv)
 	xNegOverY = pr.curveF.Neg(xNegOverY)
 
 	// init Miller loop accumulator to residueWitnessInv to share the squarings
diff --git a/std/algebra/emulated/sw_bn254/pairing_test.go b/std/algebra/emulated/sw_bn254/pairing_test.go
index 8b49607441..c6de7af3ca 100644
--- a/std/algebra/emulated/sw_bn254/pairing_test.go
+++ b/std/algebra/emulated/sw_bn254/pairing_test.go
@@ -247,26 +247,26 @@ func TestMultiPairTestSolve(t *testing.T) {
 	}
 }
 
-type PairingCheckCircuit struct {
+type DoublePairingCheckCircuit struct {
 	In1G1 G1Affine
 	In2G1 G1Affine
 	In1G2 G2Affine
 	In2G2 G2Affine
 }
 
-func (c *PairingCheckCircuit) Define(api frontend.API) error {
+func (c *DoublePairingCheckCircuit) Define(api frontend.API) error {
 	pairing, err := NewPairing(api)
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
-	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2})
+	err = pairing.DoublePairingCheck([2]*G1Affine{&c.In1G1, &c.In2G1}, [2]*G2Affine{&c.In1G2, &c.In2G2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
 	}
 	return nil
 }
 
-func TestPairingCheckTestSolve(t *testing.T) {
+func TestDoublePairingCheckTestSolve(t *testing.T) {
 	assert := test.NewAssert(t)
 	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
@@ -275,6 +275,41 @@ func TestPairingCheckTestSolve(t *testing.T) {
 	var q2 bn254.G2Affine
 	q2.Set(&q1)
 	q1.Double(&q1)
+	witness := DoublePairingCheckCircuit{
+		In1G1: NewG1Affine(p1),
+		In1G2: NewG2Affine(q1),
+		In2G1: NewG1Affine(p2),
+		In2G2: NewG2Affine(q2),
+	}
+	err := test.IsSolved(&DoublePairingCheckCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type PairingCheckCircuit struct {
+	In1G1 G1Affine
+	In2G1 G1Affine
+	In1G2 G2Affine
+	In2G2 G2Affine
+}
+
+func (c *PairingCheckCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	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})
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	return nil
+}
+
+func TestPairingCheckTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p1, q1 := randomG1G2Affines()
+	_, q2 := randomG1G2Affines()
+	var p2 bn254.G1Affine
+	p2.Neg(&p1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),
diff --git a/std/algebra/emulated/sw_bw6761/pairing.go b/std/algebra/emulated/sw_bw6761/pairing.go
index 47ad915567..2c481e75b7 100644
--- a/std/algebra/emulated/sw_bw6761/pairing.go
+++ b/std/algebra/emulated/sw_bw6761/pairing.go
@@ -161,6 +161,14 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck calculates the reduced pairing for a 2 pairs of points and asserts if the result is One
+// e(P0, Q0) * e(P1, Q1) =? 1
+//
+// This function doesn't check that the inputs are in the correct subgroups. See AssertIsOnG1 and AssertIsOnG2.
+func (pr Pairing) DoublePairingCheck(P [2]*G1Affine, Q [2]*G2Affine) error {
+	return pr.PairingCheck(P[:], Q[:])
+}
+
 func (pr Pairing) AssertIsEqual(x, y *GTEl) {
 	pr.Ext6.AssertIsEqual(x, y)
 }
diff --git a/std/algebra/emulated/sw_bw6761/pairing_test.go b/std/algebra/emulated/sw_bw6761/pairing_test.go
index 65f087a555..7f2f4683f8 100644
--- a/std/algebra/emulated/sw_bw6761/pairing_test.go
+++ b/std/algebra/emulated/sw_bw6761/pairing_test.go
@@ -163,26 +163,26 @@ func TestMultiPairTestSolve(t *testing.T) {
 	}
 }
 
-type PairingCheckCircuit struct {
+type DoublePairingCheckCircuit struct {
 	In1G1 G1Affine
 	In2G1 G1Affine
 	In1G2 G2Affine
 	In2G2 G2Affine
 }
 
-func (c *PairingCheckCircuit) Define(api frontend.API) error {
+func (c *DoublePairingCheckCircuit) Define(api frontend.API) error {
 	pairing, err := NewPairing(api)
 	if err != nil {
 		return fmt.Errorf("new pairing: %w", err)
 	}
-	err = pairing.PairingCheck([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2})
+	err = pairing.DoublePairingCheck([2]*G1Affine{&c.In1G1, &c.In2G1}, [2]*G2Affine{&c.In1G2, &c.In2G2})
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
 	}
 	return nil
 }
 
-func TestPairingCheckTestSolve(t *testing.T) {
+func TestDoublePairingCheckTestSolve(t *testing.T) {
 	assert := test.NewAssert(t)
 	// e(a,2b) * e(-2a,b) == 1
 	p1, q1 := randomG1G2Affines()
@@ -191,6 +191,41 @@ func TestPairingCheckTestSolve(t *testing.T) {
 	var q2 bw6761.G2Affine
 	q2.Set(&q1)
 	q1.Double(&q1)
+	witness := DoublePairingCheckCircuit{
+		In1G1: NewG1Affine(p1),
+		In1G2: NewG2Affine(q1),
+		In2G1: NewG1Affine(p2),
+		In2G2: NewG2Affine(q2),
+	}
+	err := test.IsSolved(&DoublePairingCheckCircuit{}, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type PairingCheckCircuit struct {
+	In1G1 G1Affine
+	In2G1 G1Affine
+	In1G2 G2Affine
+	In2G2 G2Affine
+}
+
+func (c *PairingCheckCircuit) Define(api frontend.API) error {
+	pairing, err := NewPairing(api)
+	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})
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+	return nil
+}
+
+func TestPairingCheckTestSolve(t *testing.T) {
+	assert := test.NewAssert(t)
+	p1, q1 := randomG1G2Affines()
+	_, q2 := randomG1G2Affines()
+	var p2 bw6761.G1Affine
+	p2.Neg(&p1)
 	witness := PairingCheckCircuit{
 		In1G1: NewG1Affine(p1),
 		In1G2: NewG2Affine(q1),
diff --git a/std/algebra/interfaces.go b/std/algebra/interfaces.go
index 3775d0f922..666c1ecb5d 100644
--- a/std/algebra/interfaces.go
+++ b/std/algebra/interfaces.go
@@ -97,6 +97,10 @@ type Pairing[G1El G1ElementT, G2El G2ElementT, GtEl GtElementT] interface {
 	// when the inputs are of mismatching length. It does not modify the inputs.
 	PairingCheck([]*G1El, []*G2El) error
 
+	// DoublePairingCheck asserts that the pairing result of 2 pairs of points
+	// is 1.  It does not modify the inputs.
+	DoublePairingCheck([2]*G1El, [2]*G2El) error
+
 	// AssertIsEqual asserts the equality of the inputs.
 	AssertIsEqual(*GtEl, *GtEl)
 
diff --git a/std/algebra/native/sw_bls12377/hints.go b/std/algebra/native/sw_bls12377/hints.go
index d59ef955ef..09d2550efd 100644
--- a/std/algebra/native/sw_bls12377/hints.go
+++ b/std/algebra/native/sw_bls12377/hints.go
@@ -5,6 +5,7 @@ import (
 	"math/big"
 
 	"github.com/consensys/gnark-crypto/ecc"
+	bls12377 "github.com/consensys/gnark-crypto/ecc/bls12-377"
 	"github.com/consensys/gnark/constraint/solver"
 )
 
@@ -13,6 +14,7 @@ func GetHints() []solver.Hint {
 		decomposeScalarG1,
 		decomposeScalarG1Simple,
 		decomposeScalarG2,
+		doublePairingCheckHint,
 	}
 }
 
@@ -88,3 +90,88 @@ func decomposeScalarG2(scalarField *big.Int, inputs []*big.Int, outputs []*big.I
 
 	return nil
 }
+
+func doublePairingCheckHint(_ *big.Int, inputs, outputs []*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.
+	var P0, P1 bls12377.G1Affine
+	var Q0, Q1 bls12377.G2Affine
+
+	P0.X.SetBigInt(inputs[0])
+	P0.Y.SetBigInt(inputs[1])
+	P1.X.SetBigInt(inputs[2])
+	P1.Y.SetBigInt(inputs[3])
+	Q0.X.A0.SetBigInt(inputs[4])
+	Q0.X.A1.SetBigInt(inputs[5])
+	Q0.Y.A0.SetBigInt(inputs[6])
+	Q0.Y.A1.SetBigInt(inputs[7])
+	Q1.X.A0.SetBigInt(inputs[8])
+	Q1.X.A1.SetBigInt(inputs[9])
+	Q1.Y.A0.SetBigInt(inputs[10])
+	Q1.Y.A1.SetBigInt(inputs[11])
+
+	lines0 := bls12377.PrecomputeLines(Q0)
+	lines1 := bls12377.PrecomputeLines(Q1)
+	millerLoop, err := bls12377.MillerLoopFixedQ(
+		[]bls12377.G1Affine{P0, P1},
+		[][2][len(bls12377.LoopCounter) - 1]bls12377.LineEvaluationAff{lines0, lines1},
+	)
+	if err != nil {
+		return err
+	}
+
+	var root, rootPthInverse, residueWitness, scalingFactor bls12377.E12
+	var exponent, exponentInv, finalExpFactor, polyFactor big.Int
+	// polyFactor = 12(x-1)
+	polyFactor.SetString("115033474957087604736", 10)
+	// finalExpFactor = ((q^12 - 1) / r) / polyFactor
+	finalExpFactor.SetString("92351561334497520756349650336409370070948672672207914824247073415859727964231807559847070685040742345026775319680739143654748316009031763764029886042408725311062057776702838555815712331129279611544378217895455619058809454575474763035923260395518532422855090028311239234310116353269618927871828693919559964406939845784130633021661399269804065961999062695977580539176029238189119059338698461832966347603096853909366901376879505972606045770762516580639801134008192256366142553202619529638202068488750102055204336502584141399828818871664747496033599618827160583206926869573005874449182200210044444351826855938563862937638034918413235278166699461287943529570559518592586872860190313088429391521694808994276205429071153237122495989095857292965461625387657577981811772819764071512345106346232882471034669258055302790607847924560040527682025558360106509628206144255667203317787586698694011876342903106644003067103035176245790275561392007119121995936066014208972135762663107247939004517852248103325700169848524693333524025685325993207375736519358185783520948988673594976115901587076295116293065682366935313875411927779217584729138600463438806153265891176654957439524358472291492028580820575807385461119025678550977847392818655362610734928283105671242634809807533919011078145", 10)
+
+	// 1. get pth-root inverse
+	exponent.Set(&finalExpFactor)
+	root.Exp(millerLoop, &finalExpFactor)
+	if root.IsOne() {
+		rootPthInverse.SetOne()
+	} else {
+		exponentInv.ModInverse(&exponent, &polyFactor)
+		exponent.Neg(&exponentInv).Mod(&exponent, &polyFactor)
+		rootPthInverse.Exp(root, &exponent)
+	}
+
+	// 3. shift the Miller loop result so that millerLoop * scalingFactor
+	// is of order finalExpFactor
+	scalingFactor.Set(&rootPthInverse)
+	millerLoop.Mul(&millerLoop, &scalingFactor)
+
+	// 4. get the witness residue
+	//
+	// lambda = q - u, the optimal exponent
+	var lambda big.Int
+	lambda.SetString("258664426012969094010652733694893533536393512754914660539884262666720468348340822774968888139563774001527230824448", 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
+}
diff --git a/std/algebra/native/sw_bls12377/pairing.go b/std/algebra/native/sw_bls12377/pairing.go
index 9cd1ce851b..75e2af030a 100644
--- a/std/algebra/native/sw_bls12377/pairing.go
+++ b/std/algebra/native/sw_bls12377/pairing.go
@@ -273,6 +273,130 @@ func PairingCheck(api frontend.API, P []G1Affine, Q []G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck calculates the reduced pairing for 2 pairs of points and asserts if the result is One
+// e(P0, Q0) * e(P1, Q1) =? 1
+//
+// This function doesn't check that the inputs are in the correct subgroups
+func DoublePairingCheck(api frontend.API, P [2]G1Affine, Q [2]G2Affine) error {
+	// hint the non-residue witness
+	hint, err := api.NewHint(doublePairingCheckHint, 18, P[0].X, P[0].Y, P[1].X, P[1].Y, Q[0].P.X.A0, Q[0].P.X.A1, Q[0].P.Y.A0, Q[0].P.Y.A1, Q[1].P.X.A0, Q[1].P.X.A1, Q[1].P.Y.A0, Q[1].P.Y.A1)
+	if err != nil {
+		// err is non-nil only for invalid number of inputs
+		panic(err)
+	}
+
+	var residueWitness, residueWitnessInv, scalingFactor, t0 fields_bls12377.E12
+	residueWitness.C0.B0.A0 = hint[0]
+	residueWitness.C0.B0.A1 = hint[1]
+	residueWitness.C0.B1.A0 = hint[2]
+	residueWitness.C0.B1.A1 = hint[3]
+	residueWitness.C0.B2.A0 = hint[4]
+	residueWitness.C0.B2.A1 = hint[5]
+	residueWitness.C1.B0.A0 = hint[6]
+	residueWitness.C1.B0.A1 = hint[7]
+	residueWitness.C1.B1.A0 = hint[8]
+	residueWitness.C1.B1.A1 = hint[9]
+	residueWitness.C1.B2.A0 = hint[10]
+	residueWitness.C1.B2.A1 = hint[11]
+	// constrain cubicNonResiduePower to be in Fp6
+	scalingFactor.C0.B0.A0 = hint[12]
+	scalingFactor.C0.B0.A1 = hint[13]
+	scalingFactor.C0.B1.A0 = hint[14]
+	scalingFactor.C0.B1.A1 = hint[15]
+	scalingFactor.C0.B2.A0 = hint[16]
+	scalingFactor.C0.B2.A1 = hint[17]
+	scalingFactor.C1.SetZero()
+
+	// residueWitnessInv = 1 / residueWitness
+	residueWitnessInv.Inverse(api, residueWitness)
+
+	if Q[0].Lines == nil {
+		Q0lines := computeLines(api, Q[0].P)
+		Q[0].Lines = Q0lines
+	}
+	lines0 := *Q[0].Lines
+	if Q[1].Lines == nil {
+		Q1lines := computeLines(api, Q[1].P)
+		Q[1].Lines = Q1lines
+	}
+	lines1 := *Q[1].Lines
+
+	// precomputations
+	y0Inv := api.Inverse(P[0].Y)
+	x0NegOverY0 := api.Mul(P[0].X, y0Inv)
+	x0NegOverY0 = api.Neg(x0NegOverY0)
+	y1Inv := api.Inverse(P[1].Y)
+	x1NegOverY1 := api.Mul(P[1].X, y1Inv)
+	x1NegOverY1 = api.Neg(x1NegOverY1)
+
+	// init Miller loop accumulator to residueWitness to share the squarings
+	// of residueWitnessInv^{-x₀}
+	res := residueWitness
+
+	// Compute f_{x₀,Q}(P)
+	var prodLines [5]fields_bls12377.E2
+	var l0, l1 lineEvaluation
+	for i := 62; i >= 0; i-- {
+		// mutualize the square among n Miller loops
+		// (∏ᵢfᵢ)²
+		res.Square(api, res)
+
+		if loopCounter[i] == 0 {
+			// line evaluation at P
+			// ℓ × res
+			res.MulBy034(api,
+				*l0.R0.MulByFp(api, lines0[0][i].R0, x0NegOverY0),
+				*l0.R1.MulByFp(api, lines0[0][i].R1, y0Inv),
+			)
+			// ℓ × res
+			res.MulBy034(api,
+				*l0.R0.MulByFp(api, lines1[0][i].R0, x1NegOverY1),
+				*l0.R1.MulByFp(api, lines1[0][i].R1, y1Inv),
+			)
+		} else {
+			// multiply by residueWitness when bit=1
+			res.Mul(api, res, residueWitness)
+
+			// lines evaluation at P
+			// ℓ × ℓ
+			prodLines = *fields_bls12377.Mul034By034(api,
+				*l0.R0.MulByFp(api, lines0[0][i].R0, x0NegOverY0),
+				*l0.R1.MulByFp(api, lines0[0][i].R1, y0Inv),
+				*l1.R0.MulByFp(api, lines0[1][i].R0, x0NegOverY0),
+				*l1.R1.MulByFp(api, lines0[1][i].R1, y0Inv),
+			)
+			// (ℓ × ℓ) × res
+			res.MulBy01234(api, prodLines)
+			// lines evaluation at P
+			// ℓ × ℓ
+			prodLines = *fields_bls12377.Mul034By034(api,
+				*l0.R0.MulByFp(api, lines1[0][i].R0, x1NegOverY1),
+				*l0.R1.MulByFp(api, lines1[0][i].R1, y1Inv),
+				*l1.R0.MulByFp(api, lines1[1][i].R0, x1NegOverY1),
+				*l1.R1.MulByFp(api, lines1[1][i].R1, y1Inv),
+			)
+			// (ℓ × ℓ) × res
+			res.MulBy01234(api, prodLines)
+		}
+	}
+
+	// Check that  res * scalingFactor * residueWitnessInv^λ' == 1
+	// where λ' = q, with u the BLS12-377 seed
+	// and residueWitnessInv, scalingFactor from the hint.
+	// Note that res is already MillerLoop(P,Q) * residueWitnessInv^{x₀} since
+	// we initialized the Miller loop accumulator with residueWitnessInv.
+	t0.Frobenius(api, residueWitnessInv)
+	t0.Mul(api, t0, res)
+	t0.Mul(api, t0, scalingFactor)
+
+	var one GT
+	one.SetOne()
+	t0.AssertIsEqual(api, one)
+
+	return nil
+
+}
+
 // doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates, and evaluates the line in Miller loop
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func doubleAndAddStep(api frontend.API, p1, p2 *g2AffP) (g2AffP, *lineEvaluation, *lineEvaluation) {
diff --git a/std/algebra/native/sw_bls12377/pairing2.go b/std/algebra/native/sw_bls12377/pairing2.go
index d057ea13d1..b765d731e7 100644
--- a/std/algebra/native/sw_bls12377/pairing2.go
+++ b/std/algebra/native/sw_bls12377/pairing2.go
@@ -330,6 +330,19 @@ func (p *Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck computes the multi-pairing of the 2 input pairs and asserts that
+// the result is an identity element in the target group. It returns an error if
+// there is a mismatch between the lengths of the inputs.
+func (p *Pairing) DoublePairingCheck(P [2]*G1Affine, Q [2]*G2Affine) error {
+	var inP [2]G1Affine
+	inP[0] = *P[0]
+	inP[1] = *P[1]
+	var inQ [2]G2Affine
+	inQ[0] = *Q[0]
+	inQ[1] = *Q[1]
+	return DoublePairingCheck(p.api, inP, inQ)
+}
+
 // AssertIsEqual asserts the equality of the target group elements.
 func (p *Pairing) AssertIsEqual(e1, e2 *GT) {
 	e1.AssertIsEqual(p.api, *e2)
diff --git a/std/algebra/native/sw_bls12377/pairing_test.go b/std/algebra/native/sw_bls12377/pairing_test.go
index 1956764d2a..02e1ccb1e2 100644
--- a/std/algebra/native/sw_bls12377/pairing_test.go
+++ b/std/algebra/native/sw_bls12377/pairing_test.go
@@ -57,6 +57,35 @@ func TestFinalExp(t *testing.T) {
 	assert.CheckCircuit(&finalExp{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
 }
 
+type millerLoopBLS377 struct {
+	P   G1Affine
+	Q   G2Affine
+	Res GT
+}
+
+func (circuit *millerLoopBLS377) Define(api frontend.API) error {
+	pairingRes, _ := MillerLoop(api, []G1Affine{circuit.P}, []G2Affine{circuit.Q})
+	pairingRes.AssertIsEqual(api, circuit.Res)
+
+	return nil
+}
+
+func TestMillerLoopBLS377(t *testing.T) {
+
+	// Miller loop test data
+	P, Q, milRes := millerLoopData()
+
+	// assign values to witness
+	witness := millerLoopBLS377{
+		P:   NewG1Affine(P),
+		Q:   NewG2Affine(Q),
+		Res: NewGTEl(milRes),
+	}
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&millerLoopBLS377{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
+
+}
+
 type pairingBLS377 struct {
 	P   G1Affine
 	Q   G2Affine
@@ -184,12 +213,12 @@ func TestDoublePairingFixedBLS377(t *testing.T) {
 
 }
 
-type pairingCheckBLS377 struct {
+type pairingCheck struct {
 	P1, P2 G1Affine
 	Q1, Q2 G2Affine
 }
 
-func (circuit *pairingCheckBLS377) Define(api frontend.API) error {
+func (circuit *pairingCheck) Define(api frontend.API) error {
 
 	err := PairingCheck(api, []G1Affine{circuit.P1, circuit.P2}, []G2Affine{circuit.Q1, circuit.Q2})
 
@@ -200,19 +229,48 @@ func (circuit *pairingCheckBLS377) Define(api frontend.API) error {
 	return nil
 }
 
-func TestPairingCheckBLS377(t *testing.T) {
+func TestPairingCheck(t *testing.T) {
 
 	// pairing test data
 	P, Q := pairingCheckData()
-	witness := pairingCheckBLS377{
+	witness := pairingCheck{
 		P1: NewG1Affine(P[0]),
 		P2: NewG1Affine(P[1]),
 		Q1: NewG2Affine(Q[0]),
 		Q2: NewG2Affine(Q[1]),
 	}
 	assert := test.NewAssert(t)
-	assert.CheckCircuit(&pairingCheckBLS377{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761), test.NoProverChecks())
+	assert.CheckCircuit(&pairingCheck{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761), test.NoProverChecks())
+}
 
+type doublePairingCheck struct {
+	P1, P2 G1Affine
+	Q1, Q2 G2Affine
+}
+
+func (circuit *doublePairingCheck) Define(api frontend.API) error {
+
+	err := DoublePairingCheck(api, [2]G1Affine{circuit.P1, circuit.P2}, [2]G2Affine{circuit.Q1, circuit.Q2})
+
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+
+	return nil
+}
+
+func TestDoublePairingCheck(t *testing.T) {
+
+	// pairing test data
+	P, Q := pairingCheckData()
+	witness := doublePairingCheck{
+		P1: NewG1Affine(P[0]),
+		P2: NewG1Affine(P[1]),
+		Q1: NewG2Affine(Q[0]),
+		Q2: NewG2Affine(Q[1]),
+	}
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&doublePairingCheck{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761), test.NoProverChecks())
 }
 
 type groupMembership struct {
@@ -243,6 +301,16 @@ func TestGroupMembership(t *testing.T) {
 }
 
 // utils
+func millerLoopData() (P bls12377.G1Affine, Q bls12377.G2Affine, milRes bls12377.GT) {
+	_, _, P, Q = bls12377.Generators()
+	lines := bls12377.PrecomputeLines(Q)
+	milRes, _ = bls12377.MillerLoopFixedQ(
+		[]bls12377.G1Affine{P},
+		[][2][len(bls12377.LoopCounter) - 1]bls12377.LineEvaluationAff{lines},
+	)
+	return
+}
+
 func pairingData() (P bls12377.G1Affine, Q bls12377.G2Affine, milRes, pairingRes bls12377.GT) {
 	_, _, P, Q = bls12377.Generators()
 	milRes, _ = bls12377.MillerLoop([]bls12377.G1Affine{P}, []bls12377.G2Affine{Q})
diff --git a/std/algebra/native/sw_bls24315/pairing.go b/std/algebra/native/sw_bls24315/pairing.go
index 50d741f602..1b8f02f1c0 100644
--- a/std/algebra/native/sw_bls24315/pairing.go
+++ b/std/algebra/native/sw_bls24315/pairing.go
@@ -171,10 +171,10 @@ func FinalExponentiation(api frontend.API, e1 GT) GT {
 	return result
 }
 
-// PairingCheck calculates the reduced pairing for a set of points and returns True if the result is One
-// ∏ᵢ e(Pᵢ, Qᵢ) =? 1
+// Pair calculates the reduced pairing for a set of points
+// ∏ᵢ e(Pᵢ, Qᵢ).
 //
-// This function doesn't check that the inputs are in the correct subgroup. See IsInSubGroup.
+// This function doesn't check that the inputs are in the correct subgroup
 func Pair(api frontend.API, P []G1Affine, Q []G2Affine) (GT, error) {
 	f, err := MillerLoop(api, P, Q)
 	if err != nil {
@@ -199,6 +199,14 @@ func PairingCheck(api frontend.API, P []G1Affine, Q []G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck calculates the reduced pairing for 2 pairs of points and asserts if the result is One
+// e(P0, Q0) * e(P1, Q1) =? 1
+//
+// This function doesn't check that the inputs are in the correct subgroups
+func DoublePairingCheck(api frontend.API, P [2]G1Affine, Q [2]G2Affine) error {
+	return PairingCheck(api, P[:], Q[:])
+}
+
 // doubleAndAddStep doubles p1 and adds p2 to the result in affine coordinates, and evaluates the line in Miller loop
 // https://eprint.iacr.org/2022/1162 (Section 6.1)
 func doubleAndAddStep(api frontend.API, p1, p2 *g2AffP) (g2AffP, *lineEvaluation, *lineEvaluation) {
diff --git a/std/algebra/native/sw_bls24315/pairing2.go b/std/algebra/native/sw_bls24315/pairing2.go
index fc5fac0ba8..b86d06ee60 100644
--- a/std/algebra/native/sw_bls24315/pairing2.go
+++ b/std/algebra/native/sw_bls24315/pairing2.go
@@ -313,6 +313,13 @@ func (p *Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error {
 	return nil
 }
 
+// DoublePairingCheck computes the multi-pairing of the 2 input pairs and asserts that
+// the result is an identity element in the target group. It returns an error if
+// there is a mismatch between the lengths of the inputs.
+func (p *Pairing) DoublePairingCheck(P [2]*G1Affine, Q [2]*G2Affine) error {
+	return p.PairingCheck(P[:], Q[:])
+}
+
 // AssertIsEqual asserts the equality of the target group elements.
 func (p *Pairing) AssertIsEqual(e1, e2 *GT) {
 	e1.AssertIsEqual(p.api, *e2)
diff --git a/std/algebra/native/sw_bls24315/pairing_test.go b/std/algebra/native/sw_bls24315/pairing_test.go
index 6f4ed72d97..c406556e9a 100644
--- a/std/algebra/native/sw_bls24315/pairing_test.go
+++ b/std/algebra/native/sw_bls24315/pairing_test.go
@@ -182,12 +182,12 @@ func TestDoublePairingFixedBLS315(t *testing.T) {
 
 }
 
-type pairingCheckBLS315 struct {
+type pairingCheck struct {
 	P1, P2 G1Affine
 	Q1, Q2 G2Affine
 }
 
-func (circuit *pairingCheckBLS315) Define(api frontend.API) error {
+func (circuit *pairingCheck) Define(api frontend.API) error {
 
 	err := PairingCheck(api, []G1Affine{circuit.P1, circuit.P2}, []G2Affine{circuit.Q1, circuit.Q2})
 
@@ -198,19 +198,48 @@ func (circuit *pairingCheckBLS315) Define(api frontend.API) error {
 	return nil
 }
 
-func TestPairingCheckBLS315(t *testing.T) {
+func TestPairingCheck(t *testing.T) {
 
 	// pairing test data
 	P, Q := pairingCheckData()
-	witness := pairingCheckBLS315{
+	witness := pairingCheck{
 		P1: NewG1Affine(P[0]),
 		P2: NewG1Affine(P[1]),
 		Q1: NewG2Affine(Q[0]),
 		Q2: NewG2Affine(Q[1]),
 	}
 	assert := test.NewAssert(t)
-	assert.CheckCircuit(&pairingCheckBLS315{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
+	assert.CheckCircuit(&pairingCheck{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
+}
+
+type doublePairingCheck struct {
+	P1, P2 G1Affine
+	Q1, Q2 G2Affine
+}
+
+func (circuit *doublePairingCheck) Define(api frontend.API) error {
 
+	err := DoublePairingCheck(api, [2]G1Affine{circuit.P1, circuit.P2}, [2]G2Affine{circuit.Q1, circuit.Q2})
+
+	if err != nil {
+		return fmt.Errorf("pair: %w", err)
+	}
+
+	return nil
+}
+
+func TestDoublePairingCheck(t *testing.T) {
+
+	// pairing test data
+	P, Q := pairingCheckData()
+	witness := doublePairingCheck{
+		P1: NewG1Affine(P[0]),
+		P2: NewG1Affine(P[1]),
+		Q1: NewG2Affine(Q[0]),
+		Q2: NewG2Affine(Q[1]),
+	}
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&doublePairingCheck{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
 }
 
 // utils
diff --git a/std/commitments/kzg/verifier.go b/std/commitments/kzg/verifier.go
index bf23203a18..ee4c295d46 100644
--- a/std/commitments/kzg/verifier.go
+++ b/std/commitments/kzg/verifier.go
@@ -436,9 +436,9 @@ func (v *Verifier[FR, G1El, G2El, GTEl]) CheckOpeningProof(commitment Commitment
 	totalG1 = v.curve.Add(totalG1, commitmentNeg)
 
 	// e([f(a)-f(α)-a*H(α)]G₁], G₂).e([H(α)]G₁, [α]G₂) == 1
-	if err := v.pairing.PairingCheck(
-		[]*G1El{totalG1, &proof.Quotient},
-		[]*G2El{&vk.G2[0], &vk.G2[1]},
+	if err := v.pairing.DoublePairingCheck(
+		[2]*G1El{totalG1, &proof.Quotient},
+		[2]*G2El{&vk.G2[0], &vk.G2[1]},
 	); err != nil {
 		return fmt.Errorf("pairing check: %w", err)
 	}
@@ -580,9 +580,9 @@ func (v *Verifier[FR, G1El, G2El, GTEl]) BatchVerifyMultiPoints(digests []Commit
 	}
 
 	// pairing check
-	err = v.pairing.PairingCheck(
-		[]*G1El{foldedDigest, foldedQuotients},
-		[]*G2El{&vk.G2[0], &vk.G2[1]},
+	err = v.pairing.DoublePairingCheck(
+		[2]*G1El{foldedDigest, foldedQuotients},
+		[2]*G2El{&vk.G2[0], &vk.G2[1]},
 	)
 	if err != nil {
 		return fmt.Errorf("pairingcheck: %w", err)
diff --git a/std/commitments/pedersen/verifier.go b/std/commitments/pedersen/verifier.go
index 136ad0a6bc..cfa5e84f85 100644
--- a/std/commitments/pedersen/verifier.go
+++ b/std/commitments/pedersen/verifier.go
@@ -63,7 +63,7 @@ func (v *Verifier[FR, G1El, G2El, GtEl]) AssertCommitment(commitment Commitment[
 		v.pairing.AssertIsOnG1(&knowledgeProof.G1El)
 	}
 
-	v.pairing.PairingCheck([]*G1El{&commitment.G1El, &knowledgeProof.G1El}, []*G2El{&vk.G, &vk.GRootSigmaNeg})
+	v.pairing.DoublePairingCheck([2]*G1El{&commitment.G1El, &knowledgeProof.G1El}, [2]*G2El{&vk.G, &vk.GRootSigmaNeg})
 	return nil
 }
 
diff --git a/std/recursion/groth16/verifier.go b/std/recursion/groth16/verifier.go
index dfa452174b..c8cdbed213 100644
--- a/std/recursion/groth16/verifier.go
+++ b/std/recursion/groth16/verifier.go
@@ -163,9 +163,11 @@ func PlaceholderProof[G1El algebra.G1ElementT, G2El algebra.G2ElementT](ccs cons
 // witness creation use the method [ValueOfVerifyingKey] and for stub
 // placeholder use [PlaceholderVerifyingKey].
 type VerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl algebra.GtElementT] struct {
-	E                            GtEl
-	G1                           struct{ K []G1El }
-	G2                           struct{ GammaNeg, DeltaNeg G2El }
+	G1 struct {
+		Alpha G1El
+		K     []G1El
+	}
+	G2                           struct{ GammaNeg, Beta, DeltaNeg G2El }
 	CommitmentKey                pedersen.VerifyingKey[G2El]
 	PublicAndCommitmentCommitted [][]int
 }
@@ -179,7 +181,10 @@ func PlaceholderVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, G
 	commitmentWires := commitments.CommitmentIndexes()
 
 	return VerifyingKey[G1El, G2El, GtEl]{
-		G1: struct{ K []G1El }{
+		G1: struct {
+			Alpha G1El
+			K     []G1El
+		}{
 			K: make([]G1El, ccs.GetNbPublicVariables()+len(ccs.GetCommitments().(constraint.Groth16Commitments))),
 		},
 		PublicAndCommitmentCommitted: commitments.GetPublicAndCommitmentCommitted(commitmentWires, ccs.GetNbPublicVariables()),
@@ -197,12 +202,10 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		if !ok {
 			return ret, fmt.Errorf("expected bn254.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bn254.Pair([]bn254.G1Affine{tVk.G1.Alpha}, []bn254.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bn254.NewGTEl(e)
+		var alpha bn254.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bn254.NewG1Affine(alpha)
+		s.G2.Beta = sw_bn254.NewG2Affine(tVk.G2.Beta)
 		s.G1.K = make([]sw_bn254.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bn254.NewG1Affine(tVk.G1.K[i])
@@ -212,6 +215,7 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bn254.NewG2Affine(deltaNeg)
 		s.G2.GammaNeg = sw_bn254.NewG2Affine(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKey[sw_bn254.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -219,14 +223,12 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 	case *VerifyingKey[sw_bls12377.G1Affine, sw_bls12377.G2Affine, sw_bls12377.GT]:
 		tVk, ok := vk.(*groth16backend_bls12377.VerifyingKey)
 		if !ok {
-			return ret, fmt.Errorf("expected bn254.VerifyingKey, got %T", vk)
+			return ret, fmt.Errorf("expected bls12377.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bls12377.Pair([]bls12377.G1Affine{tVk.G1.Alpha}, []bls12377.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bls12377.NewGTEl(e)
+		var alpha bls12377.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bls12377.NewG1Affine(alpha)
+		s.G2.Beta = sw_bls12377.NewG2Affine(tVk.G2.Beta)
 		s.G1.K = make([]sw_bls12377.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bls12377.NewG1Affine(tVk.G1.K[i])
@@ -236,6 +238,7 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bls12377.NewG2Affine(deltaNeg)
 		s.G2.GammaNeg = sw_bls12377.NewG2Affine(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKey[sw_bls12377.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -245,12 +248,10 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		if !ok {
 			return ret, fmt.Errorf("expected bls12381.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bls12381.Pair([]bls12381.G1Affine{tVk.G1.Alpha}, []bls12381.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bls12381.NewGTEl(e)
+		var alpha bls12381.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bls12381.NewG1Affine(alpha)
+		s.G2.Beta = sw_bls12381.NewG2Affine(tVk.G2.Beta)
 		s.G1.K = make([]sw_bls12381.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bls12381.NewG1Affine(tVk.G1.K[i])
@@ -260,6 +261,7 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bls12381.NewG2Affine(deltaNeg)
 		s.G2.GammaNeg = sw_bls12381.NewG2Affine(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKey[sw_bls12381.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -267,14 +269,12 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 	case *VerifyingKey[sw_bls24315.G1Affine, sw_bls24315.G2Affine, sw_bls24315.GT]:
 		tVk, ok := vk.(*groth16backend_bls24315.VerifyingKey)
 		if !ok {
-			return ret, fmt.Errorf("expected bls12381.VerifyingKey, got %T", vk)
-		}
-		// compute E
-		e, err := bls24315.Pair([]bls24315.G1Affine{tVk.G1.Alpha}, []bls24315.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
+			return ret, fmt.Errorf("expected bls24315.VerifyingKey, got %T", vk)
 		}
-		s.E = sw_bls24315.NewGTEl(e)
+		var alpha bls24315.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bls24315.NewG1Affine(alpha)
+		s.G2.Beta = sw_bls24315.NewG2Affine(tVk.G2.Beta)
 		s.G1.K = make([]sw_bls24315.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bls24315.NewG1Affine(tVk.G1.K[i])
@@ -284,6 +284,7 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bls24315.NewG2Affine(deltaNeg)
 		s.G2.GammaNeg = sw_bls24315.NewG2Affine(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKey[sw_bls24315.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -293,12 +294,10 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		if !ok {
 			return ret, fmt.Errorf("expected bw6761.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bw6761.Pair([]bw6761.G1Affine{tVk.G1.Alpha}, []bw6761.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bw6761.NewGTEl(e)
+		var alpha bw6761.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bw6761.NewG1Affine(alpha)
+		s.G2.Beta = sw_bw6761.NewG2Affine(tVk.G2.Beta)
 		s.G1.K = make([]sw_bw6761.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bw6761.NewG1Affine(tVk.G1.K[i])
@@ -308,6 +307,7 @@ func ValueOfVerifyingKey[G1El algebra.G1ElementT, G2El algebra.G2ElementT, GtEl
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bw6761.NewG2Affine(deltaNeg)
 		s.G2.GammaNeg = sw_bw6761.NewG2Affine(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKey[sw_bw6761.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -329,12 +329,10 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		if !ok {
 			return ret, fmt.Errorf("expected bn254.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bn254.Pair([]bn254.G1Affine{tVk.G1.Alpha}, []bn254.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bn254.NewGTEl(e)
+		var alpha bn254.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bn254.NewG1Affine(alpha)
+		s.G2.Beta = sw_bn254.NewG2AffineFixed(tVk.G2.Beta)
 		s.G1.K = make([]sw_bn254.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bn254.NewG1Affine(tVk.G1.K[i])
@@ -344,6 +342,7 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bn254.NewG2AffineFixed(deltaNeg)
 		s.G2.GammaNeg = sw_bn254.NewG2AffineFixed(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKeyFixed[sw_bn254.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -351,14 +350,12 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 	case *VerifyingKey[sw_bls12377.G1Affine, sw_bls12377.G2Affine, sw_bls12377.GT]:
 		tVk, ok := vk.(*groth16backend_bls12377.VerifyingKey)
 		if !ok {
-			return ret, fmt.Errorf("expected bn254.VerifyingKey, got %T", vk)
-		}
-		// compute E
-		e, err := bls12377.Pair([]bls12377.G1Affine{tVk.G1.Alpha}, []bls12377.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
+			return ret, fmt.Errorf("expected bls12377.VerifyingKey, got %T", vk)
 		}
-		s.E = sw_bls12377.NewGTEl(e)
+		var alpha bls12377.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bls12377.NewG1Affine(alpha)
+		s.G2.Beta = sw_bls12377.NewG2AffineFixed(tVk.G2.Beta)
 		s.G1.K = make([]sw_bls12377.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bls12377.NewG1Affine(tVk.G1.K[i])
@@ -368,6 +365,7 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bls12377.NewG2AffineFixed(deltaNeg)
 		s.G2.GammaNeg = sw_bls12377.NewG2AffineFixed(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKeyFixed[sw_bls12377.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -377,12 +375,10 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		if !ok {
 			return ret, fmt.Errorf("expected bls12381.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bls12381.Pair([]bls12381.G1Affine{tVk.G1.Alpha}, []bls12381.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bls12381.NewGTEl(e)
+		var alpha bls12381.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bls12381.NewG1Affine(alpha)
+		s.G2.Beta = sw_bls12381.NewG2AffineFixed(tVk.G2.Beta)
 		s.G1.K = make([]sw_bls12381.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bls12381.NewG1Affine(tVk.G1.K[i])
@@ -392,6 +388,7 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bls12381.NewG2AffineFixed(deltaNeg)
 		s.G2.GammaNeg = sw_bls12381.NewG2AffineFixed(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKeyFixed[sw_bls12381.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -399,14 +396,12 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 	case *VerifyingKey[sw_bls24315.G1Affine, sw_bls24315.G2Affine, sw_bls24315.GT]:
 		tVk, ok := vk.(*groth16backend_bls24315.VerifyingKey)
 		if !ok {
-			return ret, fmt.Errorf("expected bls12381.VerifyingKey, got %T", vk)
+			return ret, fmt.Errorf("expected bls24315.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bls24315.Pair([]bls24315.G1Affine{tVk.G1.Alpha}, []bls24315.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bls24315.NewGTEl(e)
+		var alpha bls24315.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bls24315.NewG1Affine(alpha)
+		s.G2.Beta = sw_bls24315.NewG2AffineFixed(tVk.G2.Beta)
 		s.G1.K = make([]sw_bls24315.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bls24315.NewG1Affine(tVk.G1.K[i])
@@ -416,6 +411,7 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bls24315.NewG2AffineFixed(deltaNeg)
 		s.G2.GammaNeg = sw_bls24315.NewG2AffineFixed(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKeyFixed[sw_bls24315.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -425,12 +421,10 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		if !ok {
 			return ret, fmt.Errorf("expected bw6761.VerifyingKey, got %T", vk)
 		}
-		// compute E
-		e, err := bw6761.Pair([]bw6761.G1Affine{tVk.G1.Alpha}, []bw6761.G2Affine{tVk.G2.Beta})
-		if err != nil {
-			return ret, fmt.Errorf("precompute pairing: %w", err)
-		}
-		s.E = sw_bw6761.NewGTEl(e)
+		var alpha bw6761.G1Affine
+		alpha.Neg(&tVk.G1.Alpha)
+		s.G1.Alpha = sw_bw6761.NewG1Affine(alpha)
+		s.G2.Beta = sw_bw6761.NewG2AffineFixed(tVk.G2.Beta)
 		s.G1.K = make([]sw_bw6761.G1Affine, len(tVk.G1.K))
 		for i := range s.G1.K {
 			s.G1.K[i] = sw_bw6761.NewG1Affine(tVk.G1.K[i])
@@ -440,6 +434,7 @@ func ValueOfVerifyingKeyFixed[G1El algebra.G1ElementT, G2El algebra.G2ElementT,
 		gammaNeg.Neg(&tVk.G2.Gamma)
 		s.G2.DeltaNeg = sw_bw6761.NewG2AffineFixed(deltaNeg)
 		s.G2.GammaNeg = sw_bw6761.NewG2AffineFixed(gammaNeg)
+		var err error
 		s.CommitmentKey, err = pedersen.ValueOfVerifyingKeyFixed[sw_bw6761.G2Affine](&tVk.CommitmentKey)
 		if err != nil {
 			return ret, fmt.Errorf("commitment key: %w", err)
@@ -638,10 +633,6 @@ func (v *Verifier[FR, G1El, G2El, GtEl]) AssertProof(vk VerifyingKey[G1El, G2El,
 		v.pairing.AssertIsOnG1(&proof.Krs)
 		v.pairing.AssertIsOnG2(&proof.Bs)
 	}
-	pairing, err := v.pairing.Pair([]*G1El{kSum, &proof.Krs, &proof.Ar}, []*G2El{&vk.G2.GammaNeg, &vk.G2.DeltaNeg, &proof.Bs})
-	if err != nil {
-		return fmt.Errorf("pairing: %w", err)
-	}
-	v.pairing.AssertIsEqual(pairing, &vk.E)
+	v.pairing.PairingCheck([]*G1El{kSum, &proof.Krs, &proof.Ar, &vk.G1.Alpha}, []*G2El{&vk.G2.GammaNeg, &vk.G2.DeltaNeg, &proof.Bs, &vk.G2.Beta})
 	return nil
 }