From 3a407c10927bda2b2ed04cc04757b4ba2bf81c32 Mon Sep 17 00:00:00 2001
From: ThomasPiellard <thomas.piellard@consensys.net>
Date: Fri, 17 Jan 2025 11:14:14 +0100
Subject: [PATCH] Feat/plonk memory optim (#1395)

---
 backend/plonk/bls12-377/prove.go              | 116 +++++++++++-----
 backend/plonk/bls12-381/prove.go              | 116 +++++++++++-----
 backend/plonk/bls24-315/prove.go              | 116 +++++++++++-----
 backend/plonk/bls24-317/prove.go              | 116 +++++++++++-----
 backend/plonk/bn254/prove.go                  | 116 +++++++++++-----
 backend/plonk/bw6-633/prove.go                | 116 +++++++++++-----
 backend/plonk/bw6-761/prove.go                | 116 +++++++++++-----
 .../zkpschemes/plonk/plonk.prove.go.tmpl      | 124 ++++++++++++------
 8 files changed, 644 insertions(+), 292 deletions(-)

diff --git a/backend/plonk/bls12-377/prove.go b/backend/plonk/bls12-377/prove.go
index 9d465eac1f..94dddc9c8b 100644
--- a/backend/plonk/bls12-377/prove.go
+++ b/backend/plonk/bls12-377/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/backend/plonk/bls12-381/prove.go b/backend/plonk/bls12-381/prove.go
index 936ae723d5..7792f9dd67 100644
--- a/backend/plonk/bls12-381/prove.go
+++ b/backend/plonk/bls12-381/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/backend/plonk/bls24-315/prove.go b/backend/plonk/bls24-315/prove.go
index d25ba73087..5b1653c4c8 100644
--- a/backend/plonk/bls24-315/prove.go
+++ b/backend/plonk/bls24-315/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/backend/plonk/bls24-317/prove.go b/backend/plonk/bls24-317/prove.go
index ebf1b860f4..8d81908bef 100644
--- a/backend/plonk/bls24-317/prove.go
+++ b/backend/plonk/bls24-317/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/backend/plonk/bn254/prove.go b/backend/plonk/bn254/prove.go
index fa9255139d..9b52a4c716 100644
--- a/backend/plonk/bn254/prove.go
+++ b/backend/plonk/bn254/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/backend/plonk/bw6-633/prove.go b/backend/plonk/bw6-633/prove.go
index 5d7f2883c6..e3079d805e 100644
--- a/backend/plonk/bw6-633/prove.go
+++ b/backend/plonk/bw6-633/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/backend/plonk/bw6-761/prove.go b/backend/plonk/bw6-761/prove.go
index 23b3fae45e..cc44cac2a4 100644
--- a/backend/plonk/bw6-761/prove.go
+++ b/backend/plonk/bw6-761/prove.go
@@ -55,8 +55,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -174,6 +172,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -377,6 +376,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -516,8 +523,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -808,14 +813,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -828,11 +847,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -844,14 +864,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -862,7 +874,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -870,22 +883,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -901,15 +914,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -957,6 +981,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -974,20 +999,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1022,6 +1046,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i < len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec) - 1; i > 0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
diff --git a/internal/generator/backend/template/zkpschemes/plonk/plonk.prove.go.tmpl b/internal/generator/backend/template/zkpschemes/plonk/plonk.prove.go.tmpl
index ce258ef5c8..f3cf53a4f9 100644
--- a/internal/generator/backend/template/zkpschemes/plonk/plonk.prove.go.tmpl
+++ b/internal/generator/backend/template/zkpschemes/plonk/plonk.prove.go.tmpl
@@ -43,8 +43,6 @@ const (
 	id_S1
 	id_S2
 	id_S3
-	id_ID
-	id_LOne
 	id_Qci // [ .. , Qc_i, Pi_i, ...]
 )
 
@@ -162,6 +160,7 @@ type instance struct {
 	blindedZ                  []fr.Element      // blindedZ is the blinded version of Z
 	quotientShardsRandomizers [2]fr.Element     // random elements for blinding the shards of the quotient
 
+	precomputedDenominators    []fr.Element // stores the denominators of the Lagrange polynomials
 	linearizedPolynomial       []fr.Element
 	linearizedPolynomialDigest kzg.Digest
 
@@ -365,6 +364,14 @@ func (s *instance) completeQk() error {
 	return nil
 }
 
+// computeLagrangeOneOnCoset computes 1/n (x**n-1)/(x-1) on coset*ωⁱ
+func (s *instance) computeLagrangeOneOnCoset(cosetExpMinusOne fr.Element, index int) fr.Element {
+	var res fr.Element
+	res.Mul(&cosetExpMinusOne, &s.domain0.CardinalityInv).
+		Mul(&res, &s.precomputedDenominators[index])
+	return res
+}
+
 func (s *instance) commitToLRO() error {
 	// wait for blinding polynomials to be initialized or context to be done
 	select {
@@ -504,8 +511,6 @@ func (s *instance) computeQuotient() (err error) {
 	identity := make([]fr.Element, n)
 	identity[1].Set(&s.beta)
 
-	s.x[id_ID] = iop.NewPolynomial(&identity, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
-	s.x[id_LOne] = iop.NewPolynomial(&lone, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
 	s.x[id_ZS] = s.x[id_Z].ShallowClone().Shift(1)
 
 	numerator, err := s.computeNumerator()
@@ -796,14 +801,28 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	cs.Set(&s.domain1.FrMultiplicativeGen)
 	css.Square(&cs)
 
-	orderingConstraint := func(u ...fr.Element) fr.Element {
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	// cosetExponentiatedToNMinusOne stores <coset>^n-1
+	var cosetExponentiatedToNMinusOne, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	orderingConstraint := func(index int, u ...fr.Element) fr.Element {
+
 		gamma := s.gamma
 
-		var a, b, c, r, l fr.Element
+		// ordering constraint
+		var a, b, c, r, l, id fr.Element
+
+		// evaluation of ID at coset*ωⁱ where i:=index
+		id.Mul(&twiddles0[index], &coset).Mul(&id, &s.beta)
 
-		a.Add(&gamma, &u[id_L]).Add(&a, &u[id_ID])
-		b.Mul(&u[id_ID], &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
-		c.Mul(&u[id_ID], &css).Add(&c, &u[id_O]).Add(&c, &gamma)
+		a.Add(&gamma, &u[id_L]).Add(&a, &id)
+		b.Mul(&id, &cs).Add(&b, &u[id_R]).Add(&b, &gamma)
+		c.Mul(&id, &css).Add(&c, &u[id_O]).Add(&c, &gamma)
 		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
 
 		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
@@ -816,11 +835,12 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		return l
 	}
 
-	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
-
-		var res fr.Element
+	localConstraint := func(index int, u ...fr.Element) fr.Element {
+		// local constraint
+		var res, lone fr.Element
+		lone = s.computeLagrangeOneOnCoset(cosetExponentiatedToNMinusOne, index)
 		res.SetOne()
-		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+		res.Sub(&u[id_Z], &res).Mul(&res, &lone)
 
 		return res
 	}
@@ -832,14 +852,6 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		shifters[i].Set(&s.domain1.Generator)
 	}
 
-	// stores the current coset shifter
-	var coset fr.Element
-	coset.SetOne()
-
-	var tmp, one fr.Element
-	one.SetOne()
-	bn := big.NewInt(int64(n))
-
 	cosetTable, err := s.domain0.CosetTable()
 	if err != nil {
 		return nil, err
@@ -850,7 +862,8 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	buf := make([]fr.Element, n)
 	var wgBuf sync.WaitGroup
 
-	allConstraints := func(i int, u ...fr.Element) fr.Element {
+	allConstraints := func(index int, u ...fr.Element) fr.Element {
+
 		// scale S1, S2, S3 by β
 		u[id_S1].Mul(&u[id_S1], &s.beta)
 		u[id_S2].Mul(&u[id_S2], &s.beta)
@@ -858,22 +871,22 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 		// blind L, R, O, Z, ZS
 		var y fr.Element
-		y = s.bp[id_Bl].Evaluate(twiddles0[i])
+		y = s.bp[id_Bl].Evaluate(twiddles0[index])
 		u[id_L].Add(&u[id_L], &y)
-		y = s.bp[id_Br].Evaluate(twiddles0[i])
+		y = s.bp[id_Br].Evaluate(twiddles0[index])
 		u[id_R].Add(&u[id_R], &y)
-		y = s.bp[id_Bo].Evaluate(twiddles0[i])
+		y = s.bp[id_Bo].Evaluate(twiddles0[index])
 		u[id_O].Add(&u[id_O], &y)
-		y = s.bp[id_Bz].Evaluate(twiddles0[i])
+		y = s.bp[id_Bz].Evaluate(twiddles0[index])
 		u[id_Z].Add(&u[id_Z], &y)
 
 		// ZS is shifted by 1; need to get correct twiddle
-		y = s.bp[id_Bz].Evaluate(twiddles0[(i+1)%int(n)])
+		y = s.bp[id_Bz].Evaluate(twiddles0[(index+1)%int(n)])
 		u[id_ZS].Add(&u[id_ZS], &y)
 
 		a := gateConstraint(u...)
-		b := orderingConstraint(u...)
-		c := ratioLocalConstraint(u...)
+		b := orderingConstraint(index, u...)
+		c := localConstraint(index, u...)
 		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
 		return c
 	}
@@ -889,15 +902,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 	m := uint64(s.domain1.Cardinality)
 	mm := uint64(64 - bits.TrailingZeros64(m))
 
+	s.precomputedDenominators = make([]fr.Element, s.domain0.Cardinality)
+	bufBatchInvert := make([]fr.Element, s.domain0.Cardinality)
+
 	for i := 0; i < rho; i++ {
 
 		coset.Mul(&coset, &shifters[i])
-		tmp.Exp(coset, bn).Sub(&tmp, &one)
+		cosetExponentiatedToNMinusOne.Exp(coset, bn).
+			Sub(&cosetExponentiatedToNMinusOne, &one)
+
+		for j := 0; j < int(s.domain0.Cardinality); j++ {
+			s.precomputedDenominators[j].
+				Mul(&coset, &twiddles0[j]).
+				Sub(&s.precomputedDenominators[j], &one)
+		}
+		batchInvert(s.precomputedDenominators, bufBatchInvert)
 
 		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
 		for _, q := range s.bp {
 			cq := q.Coefficients()
-			acc := tmp
+			acc := cosetExponentiatedToNMinusOne
 			for j := 0; j < len(cq); j++ {
 				cq[j].Mul(&cq[j], &acc)
 				acc.Mul(&acc, &shifters[i])
@@ -945,6 +969,7 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 		})
 
 		wgBuf.Wait()
+
 		if _, err := iop.Evaluate(
 			allConstraints,
 			buf,
@@ -962,20 +987,19 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 			wgBuf.Done()
 		}(i)
 
-		tmp.Inverse(&tmp)
-		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		cosetExponentiatedToNMinusOne.
+			Inverse(&cosetExponentiatedToNMinusOne)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )**-1
 		for _, q := range s.bp {
 			cq := q.Coefficients()
 			for j := 0; j < len(cq); j++ {
-				cq[j].Mul(&cq[j], &tmp)
+				cq[j].Mul(&cq[j], &cosetExponentiatedToNMinusOne)
 			}
 		}
 	}
 
 	// scale everything back
 	go func() {
-		s.x[id_ID] = nil
-		s.x[id_LOne] = nil
 		s.x[id_ZS] = nil
 		s.x[id_Qk] = nil
 
@@ -1010,6 +1034,26 @@ func (s *instance) computeNumerator() (*iop.Polynomial, error) {
 
 }
 
+// batchInvert modifies in place vec, with vec[i]<-vec[i]^{-1}, using
+// the Montgomery batch inversion trick. We don't use gnark-crypto's batchInvert
+// because we want to use a buffer preallocated, to avoid wasting memory.
+// /!\ it doesn't check that all vec's inputs or non zero, it is ensured by the size
+// of the field /!\
+func batchInvert(vec, buf []fr.Element) {
+	// local function only, vec and buf are of the same size
+	copy(buf, vec)
+	for i := 1; i<len(vec); i++ {
+		vec[i].Mul(&vec[i], &vec[i-1])
+	}
+	acc := vec[len(vec)-1]
+	acc.Inverse(&acc)
+	for i := len(vec)-1; i>0; i-- {
+		vec[i].Mul(&acc, &vec[i-1])
+		acc.Mul(&acc, &buf[i])
+	}
+	vec[0].Set(&acc)
+}
+
 func calculateNbTasks(n int) int {
 	nbAvailableCPU := runtime.NumCPU() - n
 	if nbAvailableCPU < 0 {
@@ -1263,11 +1307,11 @@ func (s *instance) innerComputeLinearizedPoly(lZeta, rZeta, oZeta, alpha, beta,
 	one.SetOne()
 	nbElmt := int64(s.domain0.Cardinality)
 	alphaSquareLagrangeZero.Set(&zeta).Exp(alphaSquareLagrangeZero, big.NewInt(nbElmt)) // ζⁿ
-	zetaNPlusTwo.Mul(&alphaSquareLagrangeZero, &zeta).Mul(&zetaNPlusTwo, &zeta)        // ζⁿ⁺²
+	zetaNPlusTwo.Mul(&alphaSquareLagrangeZero, &zeta).Mul(&zetaNPlusTwo, &zeta)         // ζⁿ⁺²
 	alphaSquareLagrangeZero.Sub(&alphaSquareLagrangeZero, &one)                         // ζⁿ - 1
-	zhZeta.Set(&alphaSquareLagrangeZero)                                               // Z_h(ζ) = ζⁿ - 1
+	zhZeta.Set(&alphaSquareLagrangeZero)                                                // Z_h(ζ) = ζⁿ - 1
 	frNbElmt.SetUint64(uint64(nbElmt))
-	den.Sub(&zeta, &one).Inverse(&den)                         // 1/(ζ-1)
+	den.Sub(&zeta, &one).Inverse(&den)                           // 1/(ζ-1)
 	alphaSquareLagrangeZero.Mul(&alphaSquareLagrangeZero, &den). // L₁ = (ζⁿ - 1)/(ζ-1)
 									Mul(&alphaSquareLagrangeZero, &alpha).
 									Mul(&alphaSquareLagrangeZero, &alpha).
@@ -1319,7 +1363,7 @@ func (s *instance) innerComputeLinearizedPoly(lZeta, rZeta, oZeta, alpha, beta,
 			}
 
 			t0.Mul(&blindedZCanonical[i], &alphaSquareLagrangeZero) // α²L₁(ζ)Z(X)
-			blindedZCanonical[i].Add(&t, &t0)                      // linPol += α²L₁(ζ)Z(X)
+			blindedZCanonical[i].Add(&t, &t0)                       // linPol += α²L₁(ζ)Z(X)
 
 			// if statistical zeroknowledge is deactivated, len(h1)=len(h2)=len(h3)=len(blindedZ)-1.
 			// Else len(h1)=len(h2)=len(blindedZCanonical)=len(h3)+1