diff --git a/README.md b/README.md
index c43caee473..e222cd5411 100644
--- a/README.md
+++ b/README.md
@@ -116,6 +116,32 @@ func main() {
 
 ```
 
+### GPU Support
+
+#### Icicle Library
+
+The following schemes and curves support experimental use of Ingomyama's Icicle GPU library for low level zk-SNARK primitives such as MSM, NTT, and polynomial operations:
+
+- [x] [Groth16](https://eprint.iacr.org/2016/260)
+
+instantiated with the following curve(s)
+
+- [x] BN254
+
+To use GPUs, add the `icicle` buildtag to your build/run commands, e.g. `go run -tags=icicle main.go`.
+
+You can then toggle on or off icicle acceleration by providing the `WithIcicleAcceleration` backend ProverOption:
+
+```go
+    // toggle on
+    proofIci, err := groth16.Prove(ccs, pk, secretWitness, backend.WithIcicleAcceleration())
+    
+    // toggle off
+    proof, err := groth16.Prove(ccs, pk, secretWitness)
+```
+
+For more information about prerequisites see the [Icicle repo](https://github.com/ingonyama-zk/icicle).
+
 ## Citing
 
 If you use `gnark` in your research a citation would be appreciated.
diff --git a/backend/backend.go b/backend/backend.go
index aab30f2df1..50f0ab15e5 100644
--- a/backend/backend.go
+++ b/backend/backend.go
@@ -62,6 +62,7 @@ type ProverConfig struct {
 	HashToFieldFn  hash.Hash
 	ChallengeHash  hash.Hash
 	KZGFoldingHash hash.Hash
+	Accelerator    string
 }
 
 // NewProverConfig returns a default ProverConfig with given prover options opts
@@ -122,6 +123,19 @@ func WithProverKZGFoldingHashFunction(hFunc hash.Hash) ProverOption {
 	}
 }
 
+// WithIcicleAcceleration requests to use [ICICLE] GPU proving backend for the
+// prover. This option requires that the program is compiled with `icicle` build
+// tag and the ICICLE dependencies are properly installed. See [ICICLE] for
+// installation description.
+//
+// [ICICLE]: https://github.com/ingonyama-zk/icicle
+func WithIcicleAcceleration() ProverOption {
+	return func(pc *ProverConfig) error {
+		pc.Accelerator = "icicle"
+		return nil
+	}
+}
+
 // VerifierOption defines option for altering the behavior of the verifier. See
 // the descriptions of functions returning instances of this type for
 // implemented options.
diff --git a/backend/groth16/bls12-377/prove.go b/backend/groth16/bls12-377/prove.go
index e715e5e7c9..d1ee7d1103 100644
--- a/backend/groth16/bls12-377/prove.go
+++ b/backend/groth16/bls12-377/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/bls12-381/prove.go b/backend/groth16/bls12-381/prove.go
index e79094f0bf..c96394908d 100644
--- a/backend/groth16/bls12-381/prove.go
+++ b/backend/groth16/bls12-381/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/bls24-315/prove.go b/backend/groth16/bls24-315/prove.go
index 8d2d5a3599..88b8218711 100644
--- a/backend/groth16/bls24-315/prove.go
+++ b/backend/groth16/bls24-315/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/bls24-317/prove.go b/backend/groth16/bls24-317/prove.go
index 41aa6c80b0..de4230d50c 100644
--- a/backend/groth16/bls24-317/prove.go
+++ b/backend/groth16/bls24-317/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/bn254/icicle/doc.go b/backend/groth16/bn254/icicle/doc.go
new file mode 100644
index 0000000000..a77a7fbde9
--- /dev/null
+++ b/backend/groth16/bn254/icicle/doc.go
@@ -0,0 +1,2 @@
+// Package icicle_bn254 implements ICICLE acceleration for BN254 Groth16 backend.
+package icicle_bn254
diff --git a/backend/groth16/bn254/icicle/icicle.go b/backend/groth16/bn254/icicle/icicle.go
new file mode 100644
index 0000000000..5b1b235d33
--- /dev/null
+++ b/backend/groth16/bn254/icicle/icicle.go
@@ -0,0 +1,513 @@
+//go:build icicle
+
+package icicle_bn254
+
+import (
+	"fmt"
+	"math/big"
+	"math/bits"
+	"time"
+	"unsafe"
+
+	curve "github.com/consensys/gnark-crypto/ecc/bn254"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fp"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/hash_to_field"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/pedersen"
+	"github.com/consensys/gnark/backend"
+	groth16_bn254 "github.com/consensys/gnark/backend/groth16/bn254"
+	"github.com/consensys/gnark/backend/groth16/internal"
+	"github.com/consensys/gnark/backend/witness"
+	"github.com/consensys/gnark/constraint"
+	cs "github.com/consensys/gnark/constraint/bn254"
+	"github.com/consensys/gnark/constraint/solver"
+	"github.com/consensys/gnark/internal/utils"
+	"github.com/consensys/gnark/logger"
+	iciclegnark "github.com/ingonyama-zk/iciclegnark/curves/bn254"
+)
+
+const HasIcicle = true
+
+func (pk *ProvingKey) setupDevicePointers() error {
+	if pk.deviceInfo != nil {
+		return nil
+	}
+	pk.deviceInfo = &deviceInfo{}
+	n := int(pk.Domain.Cardinality)
+	sizeBytes := n * fr.Bytes
+
+	/*************************  Start Domain Device Setup  ***************************/
+	copyCosetInvDone := make(chan unsafe.Pointer, 1)
+	copyCosetDone := make(chan unsafe.Pointer, 1)
+	copyDenDone := make(chan unsafe.Pointer, 1)
+	/*************************     CosetTableInv      ***************************/
+	go iciclegnark.CopyToDevice(pk.Domain.CosetTableInv, sizeBytes, copyCosetInvDone)
+
+	/*************************     CosetTable      ***************************/
+	go iciclegnark.CopyToDevice(pk.Domain.CosetTable, sizeBytes, copyCosetDone)
+
+	/*************************     Den      ***************************/
+	var denI, oneI fr.Element
+	oneI.SetOne()
+	denI.Exp(pk.Domain.FrMultiplicativeGen, big.NewInt(int64(pk.Domain.Cardinality)))
+	denI.Sub(&denI, &oneI).Inverse(&denI)
+
+	log2SizeFloor := bits.Len(uint(n)) - 1
+	denIcicleArr := []fr.Element{denI}
+	for i := 0; i < log2SizeFloor; i++ {
+		denIcicleArr = append(denIcicleArr, denIcicleArr...)
+	}
+	pow2Remainder := n - 1<<log2SizeFloor
+	for i := 0; i < pow2Remainder; i++ {
+		denIcicleArr = append(denIcicleArr, denI)
+	}
+
+	go iciclegnark.CopyToDevice(denIcicleArr, sizeBytes, copyDenDone)
+
+	/*************************     Twiddles and Twiddles Inv    ***************************/
+	twiddlesInv_d_gen, twddles_err := iciclegnark.GenerateTwiddleFactors(n, true)
+	if twddles_err != nil {
+		return twddles_err
+	}
+
+	twiddles_d_gen, twddles_err := iciclegnark.GenerateTwiddleFactors(n, false)
+	if twddles_err != nil {
+		return twddles_err
+	}
+
+	/*************************  End Domain Device Setup  ***************************/
+	pk.DomainDevice.Twiddles = twiddles_d_gen
+	pk.DomainDevice.TwiddlesInv = twiddlesInv_d_gen
+
+	pk.DomainDevice.CosetTableInv = <-copyCosetInvDone
+	pk.DomainDevice.CosetTable = <-copyCosetDone
+	pk.DenDevice = <-copyDenDone
+
+	/*************************  Start G1 Device Setup  ***************************/
+	/*************************     A      ***************************/
+	pointsBytesA := len(pk.G1.A) * fp.Bytes * 2
+	copyADone := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyPointsToDevice(pk.G1.A, pointsBytesA, copyADone) // Make a function for points
+
+	/*************************     B      ***************************/
+	pointsBytesB := len(pk.G1.B) * fp.Bytes * 2
+	copyBDone := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyPointsToDevice(pk.G1.B, pointsBytesB, copyBDone) // Make a function for points
+
+	/*************************     K      ***************************/
+	var pointsNoInfinity []curve.G1Affine
+	for i, gnarkPoint := range pk.G1.K {
+		if gnarkPoint.IsInfinity() {
+			pk.InfinityPointIndicesK = append(pk.InfinityPointIndicesK, i)
+		} else {
+			pointsNoInfinity = append(pointsNoInfinity, gnarkPoint)
+		}
+	}
+
+	pointsBytesK := len(pointsNoInfinity) * fp.Bytes * 2
+	copyKDone := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyPointsToDevice(pointsNoInfinity, pointsBytesK, copyKDone) // Make a function for points
+
+	/*************************     Z      ***************************/
+	pointsBytesZ := len(pk.G1.Z) * fp.Bytes * 2
+	copyZDone := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyPointsToDevice(pk.G1.Z, pointsBytesZ, copyZDone) // Make a function for points
+
+	/*************************  End G1 Device Setup  ***************************/
+	pk.G1Device.A = <-copyADone
+	pk.G1Device.B = <-copyBDone
+	pk.G1Device.K = <-copyKDone
+	pk.G1Device.Z = <-copyZDone
+
+	/*************************  Start G2 Device Setup  ***************************/
+	pointsBytesB2 := len(pk.G2.B) * fp.Bytes * 4
+	copyG2BDone := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyG2PointsToDevice(pk.G2.B, pointsBytesB2, copyG2BDone) // Make a function for points
+	pk.G2Device.B = <-copyG2BDone
+
+	/*************************  End G2 Device Setup  ***************************/
+	return nil
+}
+
+// Prove generates the proof of knowledge of a r1cs with full witness (secret + public part).
+func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...backend.ProverOption) (*groth16_bn254.Proof, error) {
+	opt, err := backend.NewProverConfig(opts...)
+	if err != nil {
+		return nil, fmt.Errorf("new prover config: %w", err)
+	}
+	if opt.HashToFieldFn == nil {
+		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
+	}
+	if opt.Accelerator != "icicle" {
+		return groth16_bn254.Prove(r1cs, &pk.ProvingKey, fullWitness, opts...)
+	}
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "icicle").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	if pk.deviceInfo == nil {
+		log.Debug().Msg("precomputing proving key in GPU")
+		if err := pk.setupDevicePointers(); err != nil {
+			return nil, fmt.Errorf("setup device pointers: %w", err)
+		}
+	}
+
+	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
+
+	proof := &groth16_bn254.Proof{Commitments: make([]curve.G1Affine, len(commitmentInfo))}
+
+	solverOpts := opt.SolverOpts[:len(opt.SolverOpts):len(opt.SolverOpts)]
+
+	privateCommittedValues := make([][]fr.Element, len(commitmentInfo))
+	for i := range commitmentInfo {
+		solverOpts = append(solverOpts, solver.OverrideHint(commitmentInfo[i].HintID, func(i int) solver.Hint {
+			return func(_ *big.Int, in []*big.Int, out []*big.Int) error {
+				privateCommittedValues[i] = make([]fr.Element, len(commitmentInfo[i].PrivateCommitted))
+				hashed := in[:len(commitmentInfo[i].PublicAndCommitmentCommitted)]
+				committed := in[len(hashed):]
+				for j, inJ := range committed {
+					privateCommittedValues[i][j].SetBigInt(inJ)
+				}
+
+				var err error
+				if proof.Commitments[i], err = pk.CommitmentKeys[i].Commit(privateCommittedValues[i]); err != nil {
+					return err
+				}
+
+				opt.HashToFieldFn.Write(constraint.SerializeCommitment(proof.Commitments[i].Marshal(), hashed, (fr.Bits-1)/8+1))
+				hashBts := opt.HashToFieldFn.Sum(nil)
+				opt.HashToFieldFn.Reset()
+				nbBuf := fr.Bytes
+				if opt.HashToFieldFn.Size() < fr.Bytes {
+					nbBuf = opt.HashToFieldFn.Size()
+				}
+				var res fr.Element
+				res.SetBytes(hashBts[:nbBuf])
+				res.BigInt(out[0])
+				return err
+			}
+		}(i)))
+	}
+
+	if r1cs.GkrInfo.Is() {
+		var gkrData cs.GkrSolvingData
+		solverOpts = append(solverOpts,
+			solver.OverrideHint(r1cs.GkrInfo.SolveHintID, cs.GkrSolveHint(r1cs.GkrInfo, &gkrData)),
+			solver.OverrideHint(r1cs.GkrInfo.ProveHintID, cs.GkrProveHint(r1cs.GkrInfo.HashName, &gkrData)))
+	}
+
+	_solution, err := r1cs.Solve(fullWitness, solverOpts...)
+	if err != nil {
+		return nil, err
+	}
+
+	solution := _solution.(*cs.R1CSSolution)
+	wireValues := []fr.Element(solution.W)
+
+	start := time.Now()
+
+	commitmentsSerialized := make([]byte, fr.Bytes*len(commitmentInfo))
+	for i := range commitmentInfo {
+		copy(commitmentsSerialized[fr.Bytes*i:], wireValues[commitmentInfo[i].CommitmentIndex].Marshal())
+	}
+
+	if proof.CommitmentPok, err = pedersen.BatchProve(pk.CommitmentKeys, privateCommittedValues, commitmentsSerialized); err != nil {
+		return nil, err
+	}
+
+	// H (witness reduction / FFT part)
+	var h unsafe.Pointer
+	chHDone := make(chan struct{}, 1)
+	go func() {
+		h = computeH(solution.A, solution.B, solution.C, pk)
+		solution.A = nil
+		solution.B = nil
+		solution.C = nil
+		chHDone <- struct{}{}
+	}()
+
+	// we need to copy and filter the wireValues for each multi exp
+	// as pk.G1.A, pk.G1.B and pk.G2.B may have (a significant) number of point at infinity
+	var wireValuesADevice, wireValuesBDevice iciclegnark.OnDeviceData
+	chWireValuesA, chWireValuesB := make(chan struct{}, 1), make(chan struct{}, 1)
+
+	go func() {
+		wireValuesA := make([]fr.Element, len(wireValues)-int(pk.NbInfinityA))
+		for i, j := 0, 0; j < len(wireValuesA); i++ {
+			if pk.InfinityA[i] {
+				continue
+			}
+			wireValuesA[j] = wireValues[i]
+			j++
+		}
+		wireValuesASize := len(wireValuesA)
+		scalarBytes := wireValuesASize * fr.Bytes
+
+		// Copy scalars to the device and retain ptr to them
+		copyDone := make(chan unsafe.Pointer, 1)
+		iciclegnark.CopyToDevice(wireValuesA, scalarBytes, copyDone)
+		wireValuesADevicePtr := <-copyDone
+
+		wireValuesADevice = iciclegnark.OnDeviceData{
+			P:    wireValuesADevicePtr,
+			Size: wireValuesASize,
+		}
+
+		close(chWireValuesA)
+	}()
+	go func() {
+		wireValuesB := make([]fr.Element, len(wireValues)-int(pk.NbInfinityB))
+		for i, j := 0, 0; j < len(wireValuesB); i++ {
+			if pk.InfinityB[i] {
+				continue
+			}
+			wireValuesB[j] = wireValues[i]
+			j++
+		}
+		wireValuesBSize := len(wireValuesB)
+		scalarBytes := wireValuesBSize * fr.Bytes
+
+		// Copy scalars to the device and retain ptr to them
+		copyDone := make(chan unsafe.Pointer, 1)
+		iciclegnark.CopyToDevice(wireValuesB, scalarBytes, copyDone)
+		wireValuesBDevicePtr := <-copyDone
+
+		wireValuesBDevice = iciclegnark.OnDeviceData{
+			P:    wireValuesBDevicePtr,
+			Size: wireValuesBSize,
+		}
+
+		close(chWireValuesB)
+	}()
+
+	// sample random r and s
+	var r, s big.Int
+	var _r, _s, _kr fr.Element
+	if _, err := _r.SetRandom(); err != nil {
+		return nil, err
+	}
+	if _, err := _s.SetRandom(); err != nil {
+		return nil, err
+	}
+	_kr.Mul(&_r, &_s).Neg(&_kr)
+
+	_r.BigInt(&r)
+	_s.BigInt(&s)
+
+	// computes r[δ], s[δ], kr[δ]
+	deltas := curve.BatchScalarMultiplicationG1(&pk.G1.Delta, []fr.Element{_r, _s, _kr})
+
+	var bs1, ar curve.G1Jac
+
+	computeBS1 := func() error {
+		<-chWireValuesB
+
+		if bs1, _, err = iciclegnark.MsmOnDevice(wireValuesBDevice.P, pk.G1Device.B, wireValuesBDevice.Size, true); err != nil {
+			return err
+		}
+
+		bs1.AddMixed(&pk.G1.Beta)
+		bs1.AddMixed(&deltas[1])
+
+		return nil
+	}
+
+	computeAR1 := func() error {
+		<-chWireValuesA
+
+		if ar, _, err = iciclegnark.MsmOnDevice(wireValuesADevice.P, pk.G1Device.A, wireValuesADevice.Size, true); err != nil {
+			return err
+		}
+
+		ar.AddMixed(&pk.G1.Alpha)
+		ar.AddMixed(&deltas[0])
+		proof.Ar.FromJacobian(&ar)
+
+		return nil
+	}
+
+	computeKRS := func() error {
+		var krs, krs2, p1 curve.G1Jac
+		sizeH := int(pk.Domain.Cardinality - 1) // comes from the fact the deg(H)=(n-1)+(n-1)-n=n-2
+
+		// check for small circuits as iciclegnark doesn't handle zero sizes well
+		if len(pk.G1.Z) > 0 {
+			if krs2, _, err = iciclegnark.MsmOnDevice(h, pk.G1Device.Z, sizeH, true); err != nil {
+				return err
+			}
+		}
+
+		// filter the wire values if needed
+		// TODO Perf @Tabaie worst memory allocation offender
+		toRemove := commitmentInfo.GetPrivateCommitted()
+		toRemove = append(toRemove, commitmentInfo.CommitmentIndexes())
+		scalars := filterHeap(wireValues[r1cs.GetNbPublicVariables():], r1cs.GetNbPublicVariables(), internal.ConcatAll(toRemove...))
+
+		// filter zero/infinity points since icicle doesn't handle them
+		// See https://github.com/ingonyama-zk/icicle/issues/169 for more info
+		for _, indexToRemove := range pk.InfinityPointIndicesK {
+			scalars = append(scalars[:indexToRemove], scalars[indexToRemove+1:]...)
+		}
+
+		scalarBytes := len(scalars) * fr.Bytes
+
+		copyDone := make(chan unsafe.Pointer, 1)
+		iciclegnark.CopyToDevice(scalars, scalarBytes, copyDone)
+		scalars_d := <-copyDone
+
+		krs, _, err = iciclegnark.MsmOnDevice(scalars_d, pk.G1Device.K, len(scalars), true)
+		iciclegnark.FreeDevicePointer(scalars_d)
+
+		if err != nil {
+			return err
+		}
+
+		krs.AddMixed(&deltas[2])
+
+		krs.AddAssign(&krs2)
+
+		p1.ScalarMultiplication(&ar, &s)
+		krs.AddAssign(&p1)
+
+		p1.ScalarMultiplication(&bs1, &r)
+		krs.AddAssign(&p1)
+
+		proof.Krs.FromJacobian(&krs)
+
+		return nil
+	}
+
+	computeBS2 := func() error {
+		// Bs2 (1 multi exp G2 - size = len(wires))
+		var Bs, deltaS curve.G2Jac
+
+		<-chWireValuesB
+		if Bs, _, err = iciclegnark.MsmG2OnDevice(wireValuesBDevice.P, pk.G2Device.B, wireValuesBDevice.Size, true); err != nil {
+			return err
+		}
+
+		deltaS.FromAffine(&pk.G2.Delta)
+		deltaS.ScalarMultiplication(&deltaS, &s)
+		Bs.AddAssign(&deltaS)
+		Bs.AddMixed(&pk.G2.Beta)
+
+		proof.Bs.FromJacobian(&Bs)
+		return nil
+	}
+
+	// wait for FFT to end
+	<-chHDone
+
+	// schedule our proof part computations
+	if err := computeAR1(); err != nil {
+		return nil, err
+	}
+	if err := computeBS1(); err != nil {
+		return nil, err
+	}
+	if err := computeKRS(); err != nil {
+		return nil, err
+	}
+	if err := computeBS2(); err != nil {
+		return nil, err
+	}
+
+	log.Debug().Dur("took", time.Since(start)).Msg("prover done")
+
+	// free device/GPU memory that is not needed for future proofs (scalars/hpoly)
+	go func() {
+		iciclegnark.FreeDevicePointer(wireValuesADevice.P)
+		iciclegnark.FreeDevicePointer(wireValuesBDevice.P)
+		iciclegnark.FreeDevicePointer(h)
+	}()
+
+	return proof, nil
+}
+
+// if len(toRemove) == 0, returns slice
+// else, returns a new slice without the indexes in toRemove. The first value in the slice is taken as indexes as sliceFirstIndex
+// this assumes len(slice) > len(toRemove)
+// filterHeap modifies toRemove
+func filterHeap(slice []fr.Element, sliceFirstIndex int, toRemove []int) (r []fr.Element) {
+
+	if len(toRemove) == 0 {
+		return slice
+	}
+
+	heap := utils.IntHeap(toRemove)
+	heap.Heapify()
+
+	r = make([]fr.Element, 0, len(slice))
+
+	// note: we can optimize that for the likely case where len(slice) >>> len(toRemove)
+	for i := 0; i < len(slice); i++ {
+		if len(heap) > 0 && i+sliceFirstIndex == heap[0] {
+			for len(heap) > 0 && i+sliceFirstIndex == heap[0] {
+				heap.Pop()
+			}
+			continue
+		}
+		r = append(r, slice[i])
+	}
+
+	return
+}
+
+func computeH(a, b, c []fr.Element, pk *ProvingKey) unsafe.Pointer {
+	// H part of Krs
+	// Compute H (hz=ab-c, where z=-2 on ker X^n+1 (z(x)=x^n-1))
+	// 	1 - _a = ifft(a), _b = ifft(b), _c = ifft(c)
+	// 	2 - ca = fft_coset(_a), ba = fft_coset(_b), cc = fft_coset(_c)
+	// 	3 - h = ifft_coset(ca o cb - cc)
+
+	n := len(a)
+
+	// add padding to ensure input length is domain cardinality
+	padding := make([]fr.Element, int(pk.Domain.Cardinality)-n)
+	a = append(a, padding...)
+	b = append(b, padding...)
+	c = append(c, padding...)
+	n = len(a)
+
+	sizeBytes := n * fr.Bytes
+
+	/*********** Copy a,b,c to Device Start ************/
+	// Individual channels are necessary to know which device pointers
+	// point to which vector
+	copyADone := make(chan unsafe.Pointer, 1)
+	copyBDone := make(chan unsafe.Pointer, 1)
+	copyCDone := make(chan unsafe.Pointer, 1)
+
+	go iciclegnark.CopyToDevice(a, sizeBytes, copyADone)
+	go iciclegnark.CopyToDevice(b, sizeBytes, copyBDone)
+	go iciclegnark.CopyToDevice(c, sizeBytes, copyCDone)
+
+	a_device := <-copyADone
+	b_device := <-copyBDone
+	c_device := <-copyCDone
+	/*********** Copy a,b,c to Device End ************/
+
+	computeInttNttDone := make(chan error, 1)
+	computeInttNttOnDevice := func(devicePointer unsafe.Pointer) {
+		a_intt_d := iciclegnark.INttOnDevice(devicePointer, pk.DomainDevice.TwiddlesInv, nil, n, sizeBytes, false)
+		iciclegnark.NttOnDevice(devicePointer, a_intt_d, pk.DomainDevice.Twiddles, pk.DomainDevice.CosetTable, n, n, sizeBytes, true)
+		computeInttNttDone <- nil
+		iciclegnark.FreeDevicePointer(a_intt_d)
+	}
+
+	go computeInttNttOnDevice(a_device)
+	go computeInttNttOnDevice(b_device)
+	go computeInttNttOnDevice(c_device)
+	_, _, _ = <-computeInttNttDone, <-computeInttNttDone, <-computeInttNttDone
+
+	iciclegnark.PolyOps(a_device, b_device, c_device, pk.DenDevice, n)
+
+	h := iciclegnark.INttOnDevice(a_device, pk.DomainDevice.TwiddlesInv, pk.DomainDevice.CosetTableInv, n, sizeBytes, true)
+
+	go func() {
+		iciclegnark.FreeDevicePointer(a_device)
+		iciclegnark.FreeDevicePointer(b_device)
+		iciclegnark.FreeDevicePointer(c_device)
+	}()
+
+	iciclegnark.ReverseScalars(h, n)
+
+	return h
+}
diff --git a/backend/groth16/bn254/icicle/marshal_test.go b/backend/groth16/bn254/icicle/marshal_test.go
new file mode 100644
index 0000000000..75c5a2b57e
--- /dev/null
+++ b/backend/groth16/bn254/icicle/marshal_test.go
@@ -0,0 +1,67 @@
+package icicle_bn254_test
+
+import (
+	"bytes"
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	"github.com/consensys/gnark/backend/groth16"
+	groth16_bn254 "github.com/consensys/gnark/backend/groth16/bn254"
+	icicle_bn254 "github.com/consensys/gnark/backend/groth16/bn254/icicle"
+	cs_bn254 "github.com/consensys/gnark/constraint/bn254"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/frontend/cs/r1cs"
+	"github.com/consensys/gnark/test"
+)
+
+type circuit struct {
+	A, B frontend.Variable `gnark:",public"`
+	Res  frontend.Variable
+}
+
+func (c *circuit) Define(api frontend.API) error {
+	api.AssertIsEqual(api.Mul(c.A, c.B), c.Res)
+	return nil
+}
+
+func TestMarshal(t *testing.T) {
+	assert := test.NewAssert(t)
+	ccs, err := frontend.Compile(ecc.BN254.ScalarField(), r1cs.NewBuilder, &circuit{})
+	assert.NoError(err)
+	tCcs := ccs.(*cs_bn254.R1CS)
+	nativePK := groth16_bn254.ProvingKey{}
+	nativeVK := groth16_bn254.VerifyingKey{}
+	err = groth16_bn254.Setup(tCcs, &nativePK, &nativeVK)
+	assert.NoError(err)
+
+	pk := groth16.NewProvingKey(ecc.BN254)
+	buf := new(bytes.Buffer)
+	_, err = nativePK.WriteTo(buf)
+	assert.NoError(err)
+	_, err = pk.ReadFrom(buf)
+	assert.NoError(err)
+	if pk.IsDifferent(&nativePK) {
+		t.Error("marshal output difference")
+	}
+}
+
+func TestMarshal2(t *testing.T) {
+	assert := test.NewAssert(t)
+	ccs, err := frontend.Compile(ecc.BN254.ScalarField(), r1cs.NewBuilder, &circuit{})
+	assert.NoError(err)
+	tCcs := ccs.(*cs_bn254.R1CS)
+	iciPK := icicle_bn254.ProvingKey{}
+	iciVK := groth16_bn254.VerifyingKey{}
+	err = icicle_bn254.Setup(tCcs, &iciPK, &iciVK)
+	assert.NoError(err)
+
+	nativePK := groth16_bn254.ProvingKey{}
+	buf := new(bytes.Buffer)
+	_, err = iciPK.WriteTo(buf)
+	assert.NoError(err)
+	_, err = nativePK.ReadFrom(buf)
+	assert.NoError(err)
+	if iciPK.IsDifferent(&nativePK) {
+		t.Error("marshal output difference")
+	}
+}
diff --git a/backend/groth16/bn254/icicle/noicicle.go b/backend/groth16/bn254/icicle/noicicle.go
new file mode 100644
index 0000000000..87703339ce
--- /dev/null
+++ b/backend/groth16/bn254/icicle/noicicle.go
@@ -0,0 +1,18 @@
+//go:build !icicle
+
+package icicle_bn254
+
+import (
+	"fmt"
+
+	"github.com/consensys/gnark/backend"
+	groth16_bn254 "github.com/consensys/gnark/backend/groth16/bn254"
+	"github.com/consensys/gnark/backend/witness"
+	cs "github.com/consensys/gnark/constraint/bn254"
+)
+
+const HasIcicle = false
+
+func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...backend.ProverOption) (*groth16_bn254.Proof, error) {
+	return nil, fmt.Errorf("icicle backend requested but program compiled without 'icicle' build tag")
+}
diff --git a/backend/groth16/bn254/icicle/provingkey.go b/backend/groth16/bn254/icicle/provingkey.go
new file mode 100644
index 0000000000..146a794255
--- /dev/null
+++ b/backend/groth16/bn254/icicle/provingkey.go
@@ -0,0 +1,36 @@
+package icicle_bn254
+
+import (
+	"unsafe"
+
+	groth16_bn254 "github.com/consensys/gnark/backend/groth16/bn254"
+	cs "github.com/consensys/gnark/constraint/bn254"
+)
+
+type deviceInfo struct {
+	G1Device struct {
+		A, B, K, Z unsafe.Pointer
+	}
+	DomainDevice struct {
+		Twiddles, TwiddlesInv     unsafe.Pointer
+		CosetTable, CosetTableInv unsafe.Pointer
+	}
+	G2Device struct {
+		B unsafe.Pointer
+	}
+	DenDevice             unsafe.Pointer
+	InfinityPointIndicesK []int
+}
+
+type ProvingKey struct {
+	groth16_bn254.ProvingKey
+	*deviceInfo
+}
+
+func Setup(r1cs *cs.R1CS, pk *ProvingKey, vk *groth16_bn254.VerifyingKey) error {
+	return groth16_bn254.Setup(r1cs, &pk.ProvingKey, vk)
+}
+
+func DummySetup(r1cs *cs.R1CS, pk *ProvingKey) error {
+	return groth16_bn254.DummySetup(r1cs, &pk.ProvingKey)
+}
diff --git a/backend/groth16/bn254/prove.go b/backend/groth16/bn254/prove.go
index 7bb46ef26a..b20b6c7cae 100644
--- a/backend/groth16/bn254/prove.go
+++ b/backend/groth16/bn254/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/bw6-633/prove.go b/backend/groth16/bw6-633/prove.go
index 70c8087b96..6d0ee420f0 100644
--- a/backend/groth16/bw6-633/prove.go
+++ b/backend/groth16/bw6-633/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/bw6-761/prove.go b/backend/groth16/bw6-761/prove.go
index 814c4057ad..203d1be898 100644
--- a/backend/groth16/bw6-761/prove.go
+++ b/backend/groth16/bw6-761/prove.go
@@ -67,7 +67,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)
 
diff --git a/backend/groth16/groth16.go b/backend/groth16/groth16.go
index 823e7c3f4b..25b60cebca 100644
--- a/backend/groth16/groth16.go
+++ b/backend/groth16/groth16.go
@@ -49,6 +49,7 @@ import (
 	groth16_bls24315 "github.com/consensys/gnark/backend/groth16/bls24-315"
 	groth16_bls24317 "github.com/consensys/gnark/backend/groth16/bls24-317"
 	groth16_bn254 "github.com/consensys/gnark/backend/groth16/bn254"
+	icicle_bn254 "github.com/consensys/gnark/backend/groth16/bn254/icicle"
 	groth16_bw6633 "github.com/consensys/gnark/backend/groth16/bw6-633"
 	groth16_bw6761 "github.com/consensys/gnark/backend/groth16/bw6-761"
 )
@@ -167,7 +168,6 @@ func Verify(proof Proof, vk VerifyingKey, publicWitness witness.Witness, opts ..
 //	 will produce an invalid proof
 //		internally, the solution vector to the R1CS will be filled with random values which may impact benchmarking
 func Prove(r1cs constraint.ConstraintSystem, pk ProvingKey, fullWitness witness.Witness, opts ...backend.ProverOption) (Proof, error) {
-
 	switch _r1cs := r1cs.(type) {
 	case *cs_bls12377.R1CS:
 		return groth16_bls12377.Prove(_r1cs, pk.(*groth16_bls12377.ProvingKey), fullWitness, opts...)
@@ -176,6 +176,9 @@ func Prove(r1cs constraint.ConstraintSystem, pk ProvingKey, fullWitness witness.
 		return groth16_bls12381.Prove(_r1cs, pk.(*groth16_bls12381.ProvingKey), fullWitness, opts...)
 
 	case *cs_bn254.R1CS:
+		if icicle_bn254.HasIcicle {
+			return icicle_bn254.Prove(_r1cs, pk.(*icicle_bn254.ProvingKey), fullWitness, opts...)
+		}
 		return groth16_bn254.Prove(_r1cs, pk.(*groth16_bn254.ProvingKey), fullWitness, opts...)
 
 	case *cs_bw6761.R1CS:
@@ -221,8 +224,15 @@ func Setup(r1cs constraint.ConstraintSystem) (ProvingKey, VerifyingKey, error) {
 		}
 		return &pk, &vk, nil
 	case *cs_bn254.R1CS:
-		var pk groth16_bn254.ProvingKey
 		var vk groth16_bn254.VerifyingKey
+		if icicle_bn254.HasIcicle {
+			var pk icicle_bn254.ProvingKey
+			if err := icicle_bn254.Setup(_r1cs, &pk, &vk); err != nil {
+				return nil, nil, err
+			}
+			return &pk, &vk, nil
+		}
+		var pk groth16_bn254.ProvingKey
 		if err := groth16_bn254.Setup(_r1cs, &pk, &vk); err != nil {
 			return nil, nil, err
 		}
@@ -277,6 +287,13 @@ func DummySetup(r1cs constraint.ConstraintSystem) (ProvingKey, error) {
 		}
 		return &pk, nil
 	case *cs_bn254.R1CS:
+		if icicle_bn254.HasIcicle {
+			var pk icicle_bn254.ProvingKey
+			if err := icicle_bn254.DummySetup(_r1cs, &pk); err != nil {
+				return nil, err
+			}
+			return &pk, nil
+		}
 		var pk groth16_bn254.ProvingKey
 		if err := groth16_bn254.DummySetup(_r1cs, &pk); err != nil {
 			return nil, err
@@ -318,6 +335,9 @@ func NewProvingKey(curveID ecc.ID) ProvingKey {
 	switch curveID {
 	case ecc.BN254:
 		pk = &groth16_bn254.ProvingKey{}
+		if icicle_bn254.HasIcicle {
+			pk = &icicle_bn254.ProvingKey{}
+		}
 	case ecc.BLS12_377:
 		pk = &groth16_bls12377.ProvingKey{}
 	case ecc.BLS12_381:
diff --git a/backend/plonk/bn254/icicle/doc.go b/backend/plonk/bn254/icicle/doc.go
new file mode 100644
index 0000000000..07bcd4568e
--- /dev/null
+++ b/backend/plonk/bn254/icicle/doc.go
@@ -0,0 +1,2 @@
+// Package icicle_bn254 implements ICICLE acceleration for BN254 Plonk backend.
+package icicle_bn254
diff --git a/backend/plonk/bn254/icicle/icicle.go b/backend/plonk/bn254/icicle/icicle.go
new file mode 100644
index 0000000000..93740aec2b
--- /dev/null
+++ b/backend/plonk/bn254/icicle/icicle.go
@@ -0,0 +1,1698 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by gnark DO NOT EDIT
+
+package icicle_bn254 
+
+import (
+	"context"
+	"unsafe"
+	"errors"
+	"fmt"
+	"hash"
+	"math/big"
+	"math/bits"
+	"runtime"
+	"sync"
+	"time"
+
+	"golang.org/x/sync/errgroup"
+
+	//"github.com/consensys/gnark-crypto/ecc"
+
+	curve "github.com/consensys/gnark-crypto/ecc/bn254"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/fp"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/fft"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/hash_to_field"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/iop"
+
+	kzg "github.com/consensys/gnark-crypto/ecc/bn254/kzg"
+	fiatshamir "github.com/consensys/gnark-crypto/fiat-shamir"
+	"github.com/consensys/gnark/backend"
+	"github.com/consensys/gnark/backend/witness"
+
+	"github.com/consensys/gnark/constraint"
+	cs "github.com/consensys/gnark/constraint/bn254"
+	"github.com/consensys/gnark/constraint/solver"
+	"github.com/consensys/gnark/internal/utils"
+	"github.com/consensys/gnark/logger"
+
+	iciclegnark "github.com/ingonyama-zk/iciclegnark/curves/bn254"
+)
+
+const HasIcicle = true
+
+const (
+	id_L int = iota
+	id_R
+	id_O
+	id_Z
+	id_ZS
+	id_Ql
+	id_Qr
+	id_Qm
+	id_Qo
+	id_Qk
+	id_S1
+	id_S2
+	id_S3
+	id_ID
+	id_LOne
+	id_Qci // [ .. , Qc_i, Pi_i, ...]
+)
+
+// blinding factors
+const (
+	id_Bl int = iota
+	id_Br
+	id_Bo
+	id_Bz
+	nb_blinding_polynomials
+)
+
+// blinding orders (-1 to deactivate)
+const (
+	order_blinding_L = 1
+	order_blinding_R = 1
+	order_blinding_O = 1
+	order_blinding_Z = 2
+)
+
+type Proof struct {
+
+	// Commitments to the solution vectors
+	LRO [3]kzg.Digest
+
+	// Commitment to Z, the permutation polynomial
+	Z kzg.Digest
+
+	// Commitments to h1, h2, h3 such that h = h1 + Xh2 + X**2h3 is the quotient polynomial
+	H [3]kzg.Digest
+
+	Bsb22Commitments []kzg.Digest
+
+	// Batch opening proof of h1 + zeta*h2 + zeta**2h3, linearizedPolynomial, l, r, o, s1, s2, qCPrime
+	BatchedProof kzg.BatchOpeningProof
+
+	// Opening proof of Z at zeta*mu
+	ZShiftedOpening kzg.OpeningProof
+}
+
+func (pk *ProvingKey) setupDevicePointers() error {
+	log := logger.Logger().With().Str("position", "start").Logger()
+	log.Info().Msg("setupDevicePointers")
+
+	start := time.Now()
+
+	if pk.deviceInfo != nil {
+		return nil
+	}
+
+	// TODO is [0] the correct part of the array
+	pk.deviceInfo = &deviceInfo{}
+	n := int(pk.Domain[0].Cardinality)
+	sizeBytes := n * fr.Bytes
+
+	/*************************  Start Domain Device Setup  ***************************/
+	copyCosetInvDone := make(chan unsafe.Pointer, 1)
+	copyCosetDone := make(chan unsafe.Pointer, 1)
+	copyDenDone := make(chan unsafe.Pointer, 1)
+
+	/*************************     CosetTableInv      ***************************/
+	go iciclegnark.CopyToDevice(pk.Domain[0].CosetTableInv, sizeBytes, copyCosetInvDone)
+
+	/*************************     CosetTable      ***************************/
+	go iciclegnark.CopyToDevice(pk.Domain[0].CosetTable, sizeBytes, copyCosetDone)
+
+	/*************************     Den      ***************************/
+	var denI, oneI fr.Element
+	oneI.SetOne()
+	denI.Exp(pk.Domain[0].FrMultiplicativeGen, big.NewInt(int64(pk.Domain[0].Cardinality)))
+	denI.Sub(&denI, &oneI).Inverse(&denI)
+
+	log2SizeFloor := bits.Len(uint(n)) - 1
+	denIcicleArr := []fr.Element{denI}
+	for i := 0; i < log2SizeFloor; i++ {
+		denIcicleArr = append(denIcicleArr, denIcicleArr...)
+	}
+	pow2Remainder := n - 1<<log2SizeFloor
+	for i := 0; i < pow2Remainder; i++ {
+		denIcicleArr = append(denIcicleArr, denI)
+	}
+
+	go iciclegnark.CopyToDevice(denIcicleArr, sizeBytes, copyDenDone)
+
+	/*************************     Twiddles and Twiddles Inv    ***************************/
+	log = logger.Logger().With().Str("position", "start").Logger()
+	log.Info().Msg("Generating Twiddle Factors")
+
+	twiddlesInv_d_gen, twddles_err := iciclegnark.GenerateTwiddleFactors(n, true)
+	if twddles_err != nil {
+		return twddles_err
+	}
+
+	twiddles_d_gen, twddles_err := iciclegnark.GenerateTwiddleFactors(n, false)
+	if twddles_err != nil {
+		return twddles_err
+	}
+
+	/*************************  End Domain Device Setup  ***************************/
+	pk.deviceInfo.DomainDevice.Twiddles = twiddles_d_gen
+	pk.deviceInfo.DomainDevice.TwiddlesInv = twiddlesInv_d_gen
+
+	pk.deviceInfo.DomainDevice.CosetTableInv = <-copyCosetInvDone
+	pk.deviceInfo.DomainDevice.CosetTable = <-copyCosetDone
+	pk.deviceInfo.DenDevice = <-copyDenDone
+
+	/*************************  G1 Device Setup ***************************/
+	pointsBytesG1 := len(pk.Kzg.G1) * fp.Bytes * 2
+	log.Info().Int("size", len(pk.Kzg.G1)).Int("size bytes G1", pointsBytesG1).Msg("G1 Device Setup")
+	copyG1Done := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyPointsToDevice(pk.Kzg.G1, pointsBytesG1, copyG1Done) // Make a function for points
+
+	pointsBytesLagrangeG1 := len(pk.KzgLagrange.G1) * fp.Bytes * 2
+	log.Info().Int("size", len(pk.KzgLagrange.G1)).Int("size bytes G1Lagrange", pointsBytesLagrangeG1).Msg("G1Lagrange Device Setup")
+	copyLagrangeG1Done := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyPointsToDevice(pk.KzgLagrange.G1, pointsBytesLagrangeG1, copyLagrangeG1Done) // Make a function for points
+
+	pk.deviceInfo.G1Device.G1 = <-copyG1Done
+	pk.deviceInfo.G1Device.G1Lagrange = <-copyLagrangeG1Done
+
+	log.Debug().Dur("took", time.Since(start)).Msg("Device Setup Complete")
+	return nil
+}
+
+func Prove(spr *cs.SparseR1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...backend.ProverOption) (*Proof, error) {
+
+	log := logger.Logger().With().
+		Str("curve", spr.CurveID().String()).
+		Int("nbConstraints", spr.GetNbConstraints()).
+		Str("backend", "plonk").Logger()
+
+	// parse the options
+	opt, err := backend.NewProverConfig(opts...)
+	if err != nil {
+		return nil, fmt.Errorf("get prover options: %w", err)
+	}
+
+	if pk.deviceInfo == nil {
+		log.Debug().Msg("precomputing proving key in GPU")
+		if err := pk.setupDevicePointers(); err != nil {
+			return nil, fmt.Errorf("setup device pointers: %w", err)
+		}
+	}
+
+	start := time.Now()
+
+	// init instance
+	g, ctx := errgroup.WithContext(context.Background())
+	instance, err := newInstance(ctx, spr, pk, fullWitness, &opt)
+	if err != nil {
+		return nil, fmt.Errorf("new instance: %w", err)
+	}
+
+	// solve constraints
+	g.Go(instance.solveConstraints)
+
+	// compute numerator data
+	g.Go(instance.initComputeNumerator)
+
+	// complete qk
+	g.Go(instance.completeQk)
+
+	// init blinding polynomials
+	g.Go(instance.initBlindingPolynomials)
+
+	// derive gamma, beta (copy constraint)
+	g.Go(instance.deriveGammaAndBeta)
+
+	// compute accumulating ratio for the copy constraint
+	g.Go(instance.buildRatioCopyConstraint)
+
+	// compute h
+	g.Go(instance.evaluateConstraints)
+
+	// open Z (blinded) at ωζ (proof.ZShiftedOpening)
+	g.Go(instance.openZ)
+
+	// fold the commitment to H ([H₀] + ζᵐ⁺²*[H₁] + ζ²⁽ᵐ⁺²⁾[H₂])
+	g.Go(instance.foldH)
+
+	// linearized polynomial
+	g.Go(instance.computeLinearizedPolynomial)
+
+	// Batch opening
+	g.Go(instance.batchOpening)
+
+	if err := g.Wait(); err != nil {
+		return nil, err
+	}
+
+	log.Debug().Dur("took", time.Since(start)).Msg("prover done")
+	return instance.proof, nil
+}
+
+// represents a Prover instance
+type instance struct {
+	ctx context.Context
+
+	pk    *ProvingKey
+	proof *Proof
+	spr   *cs.SparseR1CS
+	opt   *backend.ProverConfig
+
+	fs             fiatshamir.Transcript
+	kzgFoldingHash hash.Hash // for KZG folding
+	htfFunc        hash.Hash // hash to field function
+
+	// polynomials
+	x        []*iop.Polynomial // x stores tracks the polynomial we need
+	bp       []*iop.Polynomial // blinding polynomials
+	h        *iop.Polynomial   // h is the quotient polynomial
+	blindedZ []fr.Element      // blindedZ is the blinded version of Z
+
+	foldedH       []fr.Element // foldedH is the folded version of H
+	foldedHDigest kzg.Digest   // foldedHDigest is the kzg commitment of foldedH
+
+	linearizedPolynomial       []fr.Element
+	linearizedPolynomialDigest kzg.Digest
+
+	fullWitness witness.Witness
+
+	// bsb22 commitment stuff
+	commitmentInfo constraint.PlonkCommitments
+	commitmentVal  []fr.Element
+	cCommitments   []*iop.Polynomial
+
+	// challenges
+	gamma, beta, alpha, zeta fr.Element
+
+	// compute numerator data
+	cres, twiddles0, cosetTableRev, twiddlesRev []fr.Element
+
+	// channel to wait for the steps
+	chLRO,
+	chQk,
+	chbp,
+	chZ,
+	chH,
+	chRestoreLRO,
+	chZOpening,
+	chLinearizedPolynomial,
+	chFoldedH,
+	chNumeratorInit,
+	chGammaBeta chan struct{}
+}
+
+func newInstance(ctx context.Context, spr *cs.SparseR1CS, pk *ProvingKey, fullWitness witness.Witness, opts *backend.ProverConfig) (*instance, error) {
+	if opts.HashToFieldFn == nil {
+		opts.HashToFieldFn = hash_to_field.New([]byte("BSB22-Plonk"))
+	}
+	s := instance{
+		ctx:                    ctx,
+		pk:                     pk,
+		proof:                  &Proof{},
+		spr:                    spr,
+		opt:                    opts,
+		fullWitness:            fullWitness,
+		bp:                     make([]*iop.Polynomial, nb_blinding_polynomials),
+		fs:                     fiatshamir.NewTranscript(opts.ChallengeHash, "gamma", "beta", "alpha", "zeta"),
+		kzgFoldingHash:         opts.KZGFoldingHash,
+		htfFunc:                opts.HashToFieldFn,
+		chLRO:                  make(chan struct{}, 1),
+		chQk:                   make(chan struct{}, 1),
+		chbp:                   make(chan struct{}, 1),
+		chGammaBeta:            make(chan struct{}, 1),
+		chZ:                    make(chan struct{}, 1),
+		chH:                    make(chan struct{}, 1),
+		chZOpening:             make(chan struct{}, 1),
+		chLinearizedPolynomial: make(chan struct{}, 1),
+		chFoldedH:              make(chan struct{}, 1),
+		chRestoreLRO:           make(chan struct{}, 1),
+		chNumeratorInit:        make(chan struct{}, 1),
+	}
+	s.initBSB22Commitments()
+	s.setupGKRHints()
+	s.x = make([]*iop.Polynomial, id_Qci+2*len(s.commitmentInfo))
+
+	return &s, nil
+}
+
+func (s *instance) initComputeNumerator() error {
+	n := s.pk.Domain[0].Cardinality
+	s.cres = make([]fr.Element, s.pk.Domain[1].Cardinality)
+	s.twiddles0 = make([]fr.Element, n)
+	if n == 1 {
+		// edge case
+		s.twiddles0[0].SetOne()
+	} else {
+		copy(s.twiddles0, s.pk.Domain[0].Twiddles[0])
+		for i := len(s.pk.Domain[0].Twiddles[0]); i < len(s.twiddles0); i++ {
+			s.twiddles0[i].Mul(&s.twiddles0[i-1], &s.twiddles0[1])
+		}
+	}
+
+	cosetTable := s.pk.Domain[0].CosetTable
+	twiddles := s.pk.Domain[1].Twiddles[0][:n]
+
+	s.cosetTableRev = make([]fr.Element, len(cosetTable))
+	copy(s.cosetTableRev, cosetTable)
+	fft.BitReverse(s.cosetTableRev)
+
+	s.twiddlesRev = make([]fr.Element, len(twiddles))
+	copy(s.twiddlesRev, twiddles)
+	fft.BitReverse(s.twiddlesRev)
+
+	close(s.chNumeratorInit)
+
+	return nil
+}
+
+func (s *instance) initBlindingPolynomials() error {
+	s.bp[id_Bl] = getRandomPolynomial(order_blinding_L)
+	s.bp[id_Br] = getRandomPolynomial(order_blinding_R)
+	s.bp[id_Bo] = getRandomPolynomial(order_blinding_O)
+	s.bp[id_Bz] = getRandomPolynomial(order_blinding_Z)
+	close(s.chbp)
+	return nil
+}
+
+func (s *instance) initBSB22Commitments() {
+	s.commitmentInfo = s.spr.CommitmentInfo.(constraint.PlonkCommitments)
+	s.commitmentVal = make([]fr.Element, len(s.commitmentInfo)) // TODO @Tabaie get rid of this
+	s.cCommitments = make([]*iop.Polynomial, len(s.commitmentInfo))
+	s.proof.Bsb22Commitments = make([]kzg.Digest, len(s.commitmentInfo))
+
+	// override the hint for the commitment constraints
+	for i := range s.commitmentInfo {
+		s.opt.SolverOpts = append(s.opt.SolverOpts,
+			solver.OverrideHint(s.commitmentInfo[i].HintID, s.bsb22Hint(i)))
+	}
+}
+
+// Computing and verifying Bsb22 multi-commits explained in https://hackmd.io/x8KsadW3RRyX7YTCFJIkHg
+func (s *instance) bsb22Hint(commDepth int) solver.Hint {
+	return func(_ *big.Int, ins, outs []*big.Int) error {
+		var err error
+
+		res := &s.commitmentVal[commDepth]
+
+		commitmentInfo := s.spr.CommitmentInfo.(constraint.PlonkCommitments)[commDepth]
+		committedValues := make([]fr.Element, s.pk.Domain[0].Cardinality)
+		offset := s.spr.GetNbPublicVariables()
+		for i := range ins {
+			committedValues[offset+commitmentInfo.Committed[i]].SetBigInt(ins[i])
+		}
+		if _, err = committedValues[offset+commitmentInfo.CommitmentIndex].SetRandom(); err != nil { // Commitment injection constraint has qcp = 0. Safe to use for blinding.
+			return err
+		}
+		if _, err = committedValues[offset+s.spr.GetNbConstraints()-1].SetRandom(); err != nil { // Last constraint has qcp = 0. Safe to use for blinding
+			return err
+		}
+		s.cCommitments[commDepth] = iop.NewPolynomial(&committedValues, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
+		if s.proof.Bsb22Commitments[commDepth], err = kzg.Commit(s.cCommitments[commDepth].Coefficients(), s.pk.KzgLagrange); err != nil {
+			return err
+		}
+		s.cCommitments[commDepth].ToCanonical(&s.pk.Domain[0]).ToRegular()
+
+		s.htfFunc.Write(s.proof.Bsb22Commitments[commDepth].Marshal())
+		hashBts := s.htfFunc.Sum(nil)
+		s.htfFunc.Reset()
+		nbBuf := fr.Bytes
+		if s.htfFunc.Size() < fr.Bytes {
+			nbBuf = s.htfFunc.Size()
+		}
+		res.SetBytes(hashBts[:nbBuf]) // TODO @Tabaie use CommitmentIndex for this; create a new variable CommitmentConstraintIndex for other uses
+		res.BigInt(outs[0])
+
+		return nil
+	}
+}
+
+func (s *instance) setupGKRHints() {
+	if s.spr.GkrInfo.Is() {
+		var gkrData cs.GkrSolvingData
+		s.opt.SolverOpts = append(s.opt.SolverOpts,
+			solver.OverrideHint(s.spr.GkrInfo.SolveHintID, cs.GkrSolveHint(s.spr.GkrInfo, &gkrData)),
+			solver.OverrideHint(s.spr.GkrInfo.ProveHintID, cs.GkrProveHint(s.spr.GkrInfo.HashName, &gkrData)))
+	}
+}
+
+// solveConstraints computes the evaluation of the polynomials L, R, O
+// and sets x[id_L], x[id_R], x[id_O] in canonical form
+func (s *instance) solveConstraints() error {
+	_solution, err := s.spr.Solve(s.fullWitness, s.opt.SolverOpts...)
+	if err != nil {
+		return err
+	}
+	solution := _solution.(*cs.SparseR1CSSolution)
+	evaluationLDomainSmall := []fr.Element(solution.L)
+	evaluationRDomainSmall := []fr.Element(solution.R)
+	evaluationODomainSmall := []fr.Element(solution.O)
+	var wg sync.WaitGroup
+	wg.Add(2)
+	go func() {
+		s.x[id_L] = iop.NewPolynomial(&evaluationLDomainSmall, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
+		wg.Done()
+	}()
+	go func() {
+		s.x[id_R] = iop.NewPolynomial(&evaluationRDomainSmall, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
+		wg.Done()
+	}()
+
+	s.x[id_O] = iop.NewPolynomial(&evaluationODomainSmall, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
+
+	wg.Wait()
+
+	// commit to l, r, o and add blinding factors
+	if err := s.commitToLRO(); err != nil {
+		return err
+	}
+	close(s.chLRO)
+	return nil
+}
+
+func (s *instance) completeQk() error {
+	log := logger.Logger()
+	start := time.Now()
+
+	qk := s.pk.trace.Qk.Clone().ToLagrange(&s.pk.Domain[0]).ToRegular()
+	qkCoeffs := qk.Coefficients()
+
+	wWitness, ok := s.fullWitness.Vector().(fr.Vector)
+	if !ok {
+		return witness.ErrInvalidWitness
+	}
+
+	copy(qkCoeffs, wWitness[:len(s.spr.Public)])
+
+	// wait for solver to be done
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chLRO:
+	}
+
+	for i := range s.commitmentInfo {
+		qkCoeffs[s.spr.GetNbPublicVariables()+s.commitmentInfo[i].CommitmentIndex] = s.commitmentVal[i]
+	}
+
+	s.x[id_Qk] = qk
+	close(s.chQk)
+
+	log.Debug().Dur("took", time.Since(start)).Msg("Complete Qk:")
+	return nil
+}
+
+func (s *instance) commitToLRO() error {
+	// wait for blinding polynomials to be initialized or context to be done
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chbp:
+	}
+
+	g := new(errgroup.Group)
+
+	g.Go(func() (err error) {
+		s.proof.LRO[0], err = s.commitToPolyAndBlinding(s.x[id_L], s.bp[id_Bl])
+		return
+	})
+	g.Wait()
+
+	g.Go(func() (err error) {
+		s.proof.LRO[1], err = s.commitToPolyAndBlinding(s.x[id_R], s.bp[id_Br])
+		return
+	})
+	g.Wait()
+
+	g.Go(func() (err error) {
+		s.proof.LRO[2], err = s.commitToPolyAndBlinding(s.x[id_O], s.bp[id_Bo])
+		return
+	})
+	g.Wait()
+
+	return g.Wait()
+}
+
+// deriveGammaAndBeta (copy constraint)
+func (s *instance) deriveGammaAndBeta() error {
+	wWitness, ok := s.fullWitness.Vector().(fr.Vector)
+	if !ok {
+		return witness.ErrInvalidWitness
+	}
+
+	if err := bindPublicData(&s.fs, "gamma", s.pk.Vk, wWitness[:len(s.spr.Public)]); err != nil {
+		return err
+	}
+
+	// wait for LRO to be committed
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chLRO:
+	}
+
+	gamma, err := deriveRandomness(&s.fs, "gamma", &s.proof.LRO[0], &s.proof.LRO[1], &s.proof.LRO[2])
+	if err != nil {
+		return err
+	}
+
+	bbeta, err := s.fs.ComputeChallenge("beta")
+	if err != nil {
+		return err
+	}
+	s.gamma = gamma
+	s.beta.SetBytes(bbeta)
+
+	close(s.chGammaBeta)
+
+	return nil
+}
+
+// commitToPolyAndBlinding computes the KZG commitment of a polynomial p
+// in Lagrange form (large degree)
+// and add the contribution of a blinding polynomial b (small degree)
+// /!\ The polynomial p is supposed to be in Lagrange form.
+func (s *instance) commitToPolyAndBlinding(p, b *iop.Polynomial) (commit curve.G1Affine, err error) {
+	log := logger.Logger()
+	start := time.Now()
+
+	commit, err = kzg.OnDeviceCommit(p.Coefficients(), s.pk.deviceInfo.G1Device.G1Lagrange)
+
+	// we add in the blinding contribution
+	n := int(s.pk.Domain[0].Cardinality)
+	cb := deviceCommitBlindingFactor(n, b, s.pk)
+
+	commit.Add(&commit, &cb)
+
+	log.Debug().Dur("took", time.Since(start)).Int("size", p.Size()).Msg("MSM (commitToPolyAndBlinding):")
+	return
+}
+
+func (s *instance) deriveAlpha() (err error) {
+	alphaDeps := make([]*curve.G1Affine, len(s.proof.Bsb22Commitments)+1)
+	for i := range s.proof.Bsb22Commitments {
+		alphaDeps[i] = &s.proof.Bsb22Commitments[i]
+	}
+	alphaDeps[len(alphaDeps)-1] = &s.proof.Z
+	s.alpha, err = deriveRandomness(&s.fs, "alpha", alphaDeps...)
+	return err
+}
+
+func (s *instance) deriveZeta() (err error) {
+	s.zeta, err = deriveRandomness(&s.fs, "zeta", &s.proof.H[0], &s.proof.H[1], &s.proof.H[2])
+	return
+}
+
+// evaluateConstraints computes H
+func (s *instance) evaluateConstraints() (err error) {
+	log := logger.Logger()
+	start := time.Now()
+
+	// clone polys from the proving key.
+	s.x[id_Ql] = s.pk.trace.Ql.Clone()
+	s.x[id_Qr] = s.pk.trace.Qr.Clone()
+	s.x[id_Qm] = s.pk.trace.Qm.Clone()
+	s.x[id_Qo] = s.pk.trace.Qo.Clone()
+	s.x[id_S1] = s.pk.trace.S1.Clone()
+	s.x[id_S2] = s.pk.trace.S2.Clone()
+	s.x[id_S3] = s.pk.trace.S3.Clone()
+
+	for i := 0; i < len(s.commitmentInfo); i++ {
+		s.x[id_Qci+2*i] = s.pk.trace.Qcp[i].Clone()
+	}
+
+	n := s.pk.Domain[0].Cardinality
+	lone := make([]fr.Element, n)
+	lone[0].SetOne()
+
+	// wait for solver to be done
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chLRO:
+	}
+
+	for i := 0; i < len(s.commitmentInfo); i++ {
+		s.x[id_Qci+2*i+1] = s.cCommitments[i].Clone()
+	}
+
+	// wait for Z to be committed or context done
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chZ:
+	}
+
+	// derive alpha
+	if err = s.deriveAlpha(); err != nil {
+		return err
+	}
+
+	// TODO complete waste of memory find another way to do that
+	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()
+	if err != nil {
+		return err
+	}
+
+	s.h, err = divideByXMinusOne(numerator, [2]*fft.Domain{&s.pk.Domain[0], &s.pk.Domain[1]})
+	if err != nil {
+		return err
+	}
+
+	// commit to h
+	if err := commitToQuotient(s.h1(), s.h2(), s.h3(), s.proof, s.pk); err != nil {
+		return err
+	}
+
+	if err := s.deriveZeta(); err != nil {
+		return err
+	}
+
+	// wait for clean up tasks to be done
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chRestoreLRO:
+	}
+
+	close(s.chH)
+
+	log.Debug().Dur("took", time.Since(start)).Msg("evaluateConstraints:")
+	return nil
+}
+
+func (s *instance) buildRatioCopyConstraint() (err error) {
+	log := logger.Logger()
+	// wait for gamma and beta to be derived (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chGammaBeta:
+	}
+
+	// TODO @gbotrel having iop.BuildRatioCopyConstraint return something
+	// with capacity = len() + 4 would avoid extra alloc / copy during openZ
+	start := time.Now()
+	s.x[id_Z], err = iop.BuildRatioCopyConstraint(
+		[]*iop.Polynomial{
+			s.x[id_L],
+			s.x[id_R],
+			s.x[id_O],
+		},
+		s.pk.trace.S,
+		s.beta,
+		s.gamma,
+		iop.Form{Basis: iop.Lagrange, Layout: iop.Regular},
+		&s.pk.Domain[0],
+	)
+	if err != nil {
+		return err
+	}
+	log.Debug().Dur("took", time.Since(start)).Msg("FFT (BuildRatioCopyConstraint):")
+
+	// commit to the blinded version of z
+	s.proof.Z, err = s.commitToPolyAndBlinding(s.x[id_Z], s.bp[id_Bz])
+	log.Debug().Dur("took", time.Since(start)).Msg("BuildRatioCopyConstraints:")
+
+	close(s.chZ)
+
+	return
+}
+
+// open Z (blinded) at ωζ
+func (s *instance) openZ() (err error) {
+	log := logger.Logger()
+
+	// wait for H to be committed and zeta to be derived (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chH:
+	}
+	var zetaShifted fr.Element
+	zetaShifted.Mul(&s.zeta, &s.pk.Vk.Generator)
+	s.blindedZ = getBlindedCoefficients(s.x[id_Z], s.bp[id_Bz])
+
+	// open z at zeta
+	start := time.Now()
+
+	s.proof.ZShiftedOpening, err = kzg.OnDeviceOpen(s.blindedZ, zetaShifted, s.pk.deviceInfo.G1Device.G1)
+
+	if err != nil {
+		return err
+	}
+	log.Debug().Dur("took", time.Since(start)).Msg("MSM (Open Z)")
+
+	close(s.chZOpening)
+	return nil
+}
+
+func (s *instance) h1() []fr.Element {
+	h1 := s.h.Coefficients()[:s.pk.Domain[0].Cardinality+2]
+	return h1
+}
+
+func (s *instance) h2() []fr.Element {
+	h2 := s.h.Coefficients()[s.pk.Domain[0].Cardinality+2 : 2*(s.pk.Domain[0].Cardinality+2)]
+	return h2
+}
+
+func (s *instance) h3() []fr.Element {
+	h3 := s.h.Coefficients()[2*(s.pk.Domain[0].Cardinality+2) : 3*(s.pk.Domain[0].Cardinality+2)]
+	return h3
+}
+
+// fold the commitment to H ([H₀] + ζᵐ⁺²*[H₁] + ζ²⁽ᵐ⁺²⁾[H₂])
+func (s *instance) foldH() error {
+	// wait for H to be committed and zeta to be derived (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chH:
+	}
+	var n big.Int
+	n.SetUint64(s.pk.Domain[0].Cardinality + 2)
+
+	var zetaPowerNplusTwo fr.Element
+	zetaPowerNplusTwo.Exp(s.zeta, &n)
+	zetaPowerNplusTwo.BigInt(&n)
+
+	s.foldedHDigest.ScalarMultiplication(&s.proof.H[2], &n)
+	s.foldedHDigest.Add(&s.foldedHDigest, &s.proof.H[1])       // ζᵐ⁺²*Comm(h3)
+	s.foldedHDigest.ScalarMultiplication(&s.foldedHDigest, &n) // ζ²⁽ᵐ⁺²⁾*Comm(h3) + ζᵐ⁺²*Comm(h2)
+	s.foldedHDigest.Add(&s.foldedHDigest, &s.proof.H[0])
+
+	// fold H (H₀ + ζᵐ⁺²*H₁ + ζ²⁽ᵐ⁺²⁾H₂))
+	h1 := s.h1()
+	h2 := s.h2()
+	s.foldedH = s.h3()
+
+	for i := 0; i < int(s.pk.Domain[0].Cardinality)+2; i++ {
+		s.foldedH[i].
+			Mul(&s.foldedH[i], &zetaPowerNplusTwo).
+			Add(&s.foldedH[i], &h2[i]).
+			Mul(&s.foldedH[i], &zetaPowerNplusTwo).
+			Add(&s.foldedH[i], &h1[i])
+	}
+
+	close(s.chFoldedH)
+
+	return nil
+}
+
+func (s *instance) computeLinearizedPolynomial() error {
+	log := logger.Logger()
+
+	// wait for H to be committed and zeta to be derived (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chH:
+	}
+
+	qcpzeta := make([]fr.Element, len(s.commitmentInfo))
+	var blzeta, brzeta, bozeta fr.Element
+	var wg sync.WaitGroup
+	wg.Add(3 + len(s.commitmentInfo))
+
+	for i := 0; i < len(s.commitmentInfo); i++ {
+		go func(i int) {
+			qcpzeta[i] = s.pk.trace.Qcp[i].Evaluate(s.zeta)
+			wg.Done()
+		}(i)
+	}
+
+	go func() {
+		blzeta = evaluateBlinded(s.x[id_L], s.bp[id_Bl], s.zeta)
+		wg.Done()
+	}()
+
+	go func() {
+		brzeta = evaluateBlinded(s.x[id_R], s.bp[id_Br], s.zeta)
+		wg.Done()
+	}()
+
+	go func() {
+		bozeta = evaluateBlinded(s.x[id_O], s.bp[id_Bo], s.zeta)
+		wg.Done()
+	}()
+
+	// wait for Z to be opened at zeta (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chZOpening:
+	}
+	bzuzeta := s.proof.ZShiftedOpening.ClaimedValue
+
+	wg.Wait()
+
+	s.linearizedPolynomial = computeLinearizedPolynomial(
+		blzeta,
+		brzeta,
+		bozeta,
+		s.alpha,
+		s.beta,
+		s.gamma,
+		s.zeta,
+		bzuzeta,
+		qcpzeta,
+		s.blindedZ,
+		coefficients(s.cCommitments),
+		s.pk,
+	)
+
+	var err error
+	timeCommit := time.Now()
+	s.linearizedPolynomialDigest, err = kzg.OnDeviceCommit(s.linearizedPolynomial, s.pk.deviceInfo.G1Device.G1, runtime.NumCPU()*2)
+	if err != nil {
+		return err
+	}
+	log.Debug().Dur("took", time.Since(timeCommit)).Int("size", len(s.linearizedPolynomial)).Msg("MSM (linearizedPolynomial):")
+	close(s.chLinearizedPolynomial)
+	return nil
+}
+
+func (s *instance) batchOpening() error {
+	log := logger.Logger()
+	polysQcp := coefficients(s.pk.trace.Qcp)
+	polysToOpen := make([][]fr.Element, 7+len(polysQcp))
+	copy(polysToOpen[7:], polysQcp)
+
+	// wait for LRO to be committed (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chLRO:
+	}
+
+	// wait for foldedH to be computed (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chFoldedH:
+	}
+
+	// wait for linearizedPolynomial to be computed (or ctx.Done())
+	select {
+	case <-s.ctx.Done():
+		return errContextDone
+	case <-s.chLinearizedPolynomial:
+	}
+
+	polysToOpen[0] = s.foldedH
+	polysToOpen[1] = s.linearizedPolynomial
+	polysToOpen[2] = getBlindedCoefficients(s.x[id_L], s.bp[id_Bl])
+	polysToOpen[3] = getBlindedCoefficients(s.x[id_R], s.bp[id_Br])
+	polysToOpen[4] = getBlindedCoefficients(s.x[id_O], s.bp[id_Bo])
+	polysToOpen[5] = s.pk.trace.S1.Coefficients()
+	polysToOpen[6] = s.pk.trace.S2.Coefficients()
+
+	digestsToOpen := make([]curve.G1Affine, len(s.pk.Vk.Qcp)+7)
+	copy(digestsToOpen[7:], s.pk.Vk.Qcp)
+
+	digestsToOpen[0] = s.foldedHDigest
+	digestsToOpen[1] = s.linearizedPolynomialDigest
+	digestsToOpen[2] = s.proof.LRO[0]
+	digestsToOpen[3] = s.proof.LRO[1]
+	digestsToOpen[4] = s.proof.LRO[2]
+	digestsToOpen[5] = s.pk.Vk.S[0]
+	digestsToOpen[6] = s.pk.Vk.S[1]
+
+	var err error
+	start := time.Now()
+	s.proof.BatchedProof, err = kzg.OnDeviceBatchOpenSinglePoint(
+		polysToOpen,
+		digestsToOpen,
+		s.zeta,
+		s.kzgFoldingHash,
+		s.pk.Kzg,
+		s.pk.deviceInfo.G1Device.G1,
+		s.proof.ZShiftedOpening.ClaimedValue.Marshal(),
+	)
+
+	log.Debug().Dur("took", time.Since(start)).Msg("MSM (batchOpeningSinglePoint):")
+
+	return err
+}
+
+// evaluate the full set of constraints, all polynomials in x are back in
+// canonical regular form at the end
+func (s *instance) computeNumerator() (*iop.Polynomial, error) {
+	log := logger.Logger()
+	start := time.Now()
+
+	n := s.pk.Domain[0].Cardinality
+
+	nbBsbGates := (len(s.x) - id_Qci + 1) >> 1
+
+	gateConstraint := func(u ...fr.Element) fr.Element {
+
+		var ic, tmp fr.Element
+
+		ic.Mul(&u[id_Ql], &u[id_L])
+		tmp.Mul(&u[id_Qr], &u[id_R])
+		ic.Add(&ic, &tmp)
+		tmp.Mul(&u[id_Qm], &u[id_L]).Mul(&tmp, &u[id_R])
+		ic.Add(&ic, &tmp)
+		tmp.Mul(&u[id_Qo], &u[id_O])
+		ic.Add(&ic, &tmp).Add(&ic, &u[id_Qk])
+		for i := 0; i < nbBsbGates; i++ {
+			tmp.Mul(&u[id_Qci+2*i], &u[id_Qci+2*i+1])
+			ic.Add(&ic, &tmp)
+		}
+
+		return ic
+	}
+
+	var cs, css fr.Element
+	cs.Set(&s.pk.Domain[1].FrMultiplicativeGen)
+	css.Square(&cs)
+
+	orderingConstraint := func(u ...fr.Element) fr.Element {
+		gamma := s.gamma
+
+		var a, b, c, r, l fr.Element
+
+		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)
+		r.Mul(&a, &b).Mul(&r, &c).Mul(&r, &u[id_Z])
+
+		a.Add(&u[id_S1], &u[id_L]).Add(&a, &gamma)
+		b.Add(&u[id_S2], &u[id_R]).Add(&b, &gamma)
+		c.Add(&u[id_S3], &u[id_O]).Add(&c, &gamma)
+		l.Mul(&a, &b).Mul(&l, &c).Mul(&l, &u[id_ZS])
+
+		l.Sub(&l, &r)
+
+		return l
+	}
+
+	ratioLocalConstraint := func(u ...fr.Element) fr.Element {
+
+		var res fr.Element
+		res.SetOne()
+		res.Sub(&u[id_Z], &res).Mul(&res, &u[id_LOne])
+
+		return res
+	}
+
+	rho := int(s.pk.Domain[1].Cardinality / n)
+	shifters := make([]fr.Element, rho)
+	shifters[0].Set(&s.pk.Domain[1].FrMultiplicativeGen)
+	for i := 1; i < rho; i++ {
+		shifters[i].Set(&s.pk.Domain[1].Generator)
+	}
+
+	// stores the current coset shifter
+	var coset fr.Element
+	coset.SetOne()
+
+	var tmp, one fr.Element
+	one.SetOne()
+	bn := big.NewInt(int64(n))
+
+	// wait for init go routine
+	<-s.chNumeratorInit
+
+
+	// Copy the twiddles slice to the device
+	sizeBytes := len(s.x[0].Coefficients()) * fr.Bytes
+
+	copyTDone := make(chan unsafe.Pointer, 1)
+	go iciclegnark.CopyToDevice(s.pk.Domain[1].Twiddles[0][:n], sizeBytes, copyTDone)
+	twiddlesPtr := <-copyTDone
+
+	// init the result polynomial & buffer
+	cres := s.cres
+	buf := make([]fr.Element, n)
+	var wgBuf sync.WaitGroup
+	
+	allConstraints := func(i 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)
+		u[id_S3].Mul(&u[id_S3], &s.beta)
+
+		// blind L, R, O, Z, ZS
+		var y fr.Element
+		y = s.bp[id_Bl].Evaluate(s.twiddles0[i])
+		u[id_L].Add(&u[id_L], &y)
+		y = s.bp[id_Br].Evaluate(s.twiddles0[i])
+		u[id_R].Add(&u[id_R], &y)
+		y = s.bp[id_Bo].Evaluate(s.twiddles0[i])
+		u[id_O].Add(&u[id_O], &y)
+		y = s.bp[id_Bz].Evaluate(s.twiddles0[i])
+		u[id_Z].Add(&u[id_Z], &y)
+
+		// ZS is shifted by 1; need to get correct twiddle
+		y = s.bp[id_Bz].Evaluate(s.twiddles0[(i+1)%int(n)])
+		u[id_ZS].Add(&u[id_ZS], &y)
+
+		a := gateConstraint(u...)
+		b := orderingConstraint(u...)
+		c := ratioLocalConstraint(u...)
+		c.Mul(&c, &s.alpha).Add(&c, &b).Mul(&c, &s.alpha).Add(&c, &a)
+		return c
+	}
+
+	// select the correct scaling vector to scale by shifter[i]
+	selectScalingVector := func(i int, l iop.Layout) unsafe.Pointer {
+		var w unsafe.Pointer
+		if i == 0 {
+			if l == iop.Regular {
+				w = s.pk.deviceInfo.DomainDevice.CosetTable
+			} else {
+				w = s.pk.deviceInfo.DomainDevice.CosetTableInv
+			}
+		} else {
+			if l == iop.Regular {
+				w = twiddlesPtr
+			} else {
+				w = s.pk.deviceInfo.DomainDevice.TwiddlesInv
+			}
+		}
+		return w
+	}
+
+	// pre-computed to compute the bit reverse index
+	// of the result polynomial
+	m := uint64(s.pk.Domain[1].Cardinality)
+	mm := uint64(64 - bits.TrailingZeros64(m))
+
+	// Set up device pointers to polynomial coefficients
+	var devicePointers []unsafe.Pointer
+	for i := 0; i < len(s.x); i++ {
+		sizeBytes := s.x[i].Size() * fr.Bytes
+
+		copyADone := make(chan unsafe.Pointer, 1)
+		go iciclegnark.CopyToDevice(s.x[i].Coefficients(), sizeBytes, copyADone)
+		a_device := <-copyADone
+
+		devicePointers = append(devicePointers, a_device)
+	}
+
+	for i := 0; i < rho; i++ {
+
+		coset.Mul(&coset, &shifters[i])
+		tmp.Exp(coset, bn).Sub(&tmp, &one)
+
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		for _, q := range s.bp {
+			cq := q.Coefficients()
+			acc := tmp
+			for j := 0; j < len(cq); j++ {
+				cq[j].Mul(&cq[j], &acc)
+				acc.Mul(&acc, &shifters[i])
+			}
+		}
+
+		batchTime := time.Now()
+		batchApply(s.x, func(p *iop.Polynomial, pn int) {
+			nbTasks := calculateNbTasks(len(s.x)-1) * 2
+
+			n := p.Size()
+			sizeBytes := p.Size() * fr.Bytes
+
+			// scale by shifter[i]
+			w_device := selectScalingVector(i, p.Layout)
+
+			// Initialize channels
+			computeInttNttDone := make(chan error, 1)
+
+			// INTT and NTT on device
+			computeInttNttOnDevice := func(scaleVecPtr, devicePointer unsafe.Pointer) {
+				a_intt_d := iciclegnark.INttOnDevice(devicePointer, s.pk.deviceInfo.DomainDevice.TwiddlesInv, nil, n, sizeBytes, false)
+
+				// Multiply by shifter[i]
+				iciclegnark.VecMulOnDevice(a_intt_d, scaleVecPtr, n)
+				iciclegnark.NttOnDevice(devicePointer, a_intt_d, s.pk.deviceInfo.DomainDevice.Twiddles, nil, n, n, sizeBytes, false)
+				iciclegnark.MontConvOnDevice(devicePointer, n, true)
+
+				// Copy results back to host
+				res := iciclegnark.CopyScalarsToHost(devicePointer, n, sizeBytes)
+
+				// Set the coefficients of the polynomial
+				cp := p.Coefficients()
+				utils.Parallelize(len(cp), func(start, end int) {
+					for j := start; j < end; j++ {
+						cp[j].Set(&res[j])
+					}
+				}, nbTasks)
+
+				// Convert back to Monticarlo form
+				iciclegnark.MontConvOnDevice(devicePointer, n, false)
+
+				computeInttNttDone <- nil
+				iciclegnark.FreeDevicePointer(a_intt_d)
+			}
+
+			// NTT on device
+			computeNttDone := make(chan error, 1)
+			computeNttOnDevice := func(scaleVecPtr, devicePointer unsafe.Pointer) {
+				// Multiply by shifter[i]
+				iciclegnark.VecMulOnDevice(devicePointer, scaleVecPtr, n)
+				iciclegnark.NttOnDevice(devicePointer, devicePointer, s.pk.deviceInfo.DomainDevice.Twiddles, nil, n, n, sizeBytes, false)
+				iciclegnark.MontConvOnDevice(devicePointer, n, true)
+
+				// Copy back to host
+				res := iciclegnark.CopyScalarsToHost(devicePointer, n, sizeBytes)
+				s.x[pn] = iop.NewPolynomial(&res, iop.Form{Basis: iop.Lagrange, Layout: iop.Regular})
+
+				// Convert back to Monticarlo form
+				iciclegnark.MontConvOnDevice(devicePointer, n, false)
+
+				computeNttDone <- nil
+			}
+
+			// Run INTT and NTT on device
+			if p.Basis == iop.Lagrange {
+				go computeInttNttOnDevice(w_device, devicePointers[pn])
+				_ = <-computeInttNttDone
+			} else {
+				go computeNttOnDevice(w_device, devicePointers[pn])
+				_ = <- computeNttDone
+			}
+
+		})
+		log.Debug().Dur("took", time.Since(batchTime)).Msg("FFT (batchApply)")
+		 
+		wgBuf.Wait()
+		if _, err := iop.Evaluate(
+			allConstraints,
+			buf,
+			iop.Form{Basis: iop.Lagrange, Layout: iop.Regular},
+			s.x...,
+		); err != nil {
+			return nil, err
+		}
+		wgBuf.Add(1)
+		go func(i int) {
+			for j := 0; j < int(n); j++ {
+				// we build the polynomial in bit reverse order
+				cres[bits.Reverse64(uint64(rho*j+i))>>mm] = buf[j]
+			}
+			wgBuf.Done()
+		}(i)
+
+		tmp.Inverse(&tmp)
+		// bl <- bl *( (s*ωⁱ)ⁿ-1 )s
+		for _, q := range s.bp {
+			cq := q.Coefficients()
+			for j := 0; j < len(cq); j++ {
+				cq[j].Mul(&cq[j], &tmp)
+			}
+		}
+	}
+
+	// scale everything back
+	go func() {
+		batchTime := time.Now()
+		for i := id_ZS; i < len(s.x); i++ {
+			s.x[i] = nil
+		}
+
+		var cs fr.Element
+		cs.Set(&shifters[0])
+		for i := 1; i < len(shifters); i++ {
+			cs.Mul(&cs, &shifters[i])
+		}
+		cs.Inverse(&cs)
+
+		// send the scale back factor to the GPU
+		var acc fr.Element
+		acc.SetOne()
+		accList := make([]fr.Element, s.x[0].Size())
+		for j := 0; j < s.x[0].Size(); j++ {
+			accList[j].Set(&acc)
+			acc.Mul(&acc, &cs)
+		}
+
+		sizeBytes := s.x[0].Size() * fr.Bytes
+		copyWDone := make(chan unsafe.Pointer, 1)
+		go iciclegnark.CopyToDevice(accList, sizeBytes, copyWDone)
+		w_device := <-copyWDone
+
+		batchApply(s.x[:id_ZS], func(p *iop.Polynomial, pn int) {
+			// ON Device
+			n := p.Size()
+			sizeBytes := p.Size() * fr.Bytes
+
+			// Initialize channels
+			computeInttNttDone := make(chan error, 1)
+			computeInttNttOnDevice := func(scaleVecPtr, devicePointer unsafe.Pointer) {
+				a_intt_d := iciclegnark.INttOnDevice(devicePointer, s.pk.deviceInfo.DomainDevice.TwiddlesInv, nil, n, sizeBytes, false)
+
+				iciclegnark.VecMulOnDevice(a_intt_d, scaleVecPtr, n)
+				iciclegnark.MontConvOnDevice(a_intt_d, n, true)
+
+				res := iciclegnark.CopyScalarsToHost(a_intt_d, n, sizeBytes)
+				s.x[pn] = iop.NewPolynomial(&res, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
+
+				// Convert back to Monticarlo form
+				iciclegnark.MontConvOnDevice(a_intt_d, n, false)
+
+				computeInttNttDone <- nil
+				iciclegnark.FreeDevicePointer(a_intt_d)
+			}
+			// Run computeInttNttOnDevice on device
+			go computeInttNttOnDevice(w_device, devicePointers[pn])
+			_ = <-computeInttNttDone
+
+		})
+		for _, q := range s.bp {
+			scalePowers(q, cs)
+		}
+
+		log.Debug().Dur("took", time.Since(batchTime)).Msg("FFT (Scale back batchApply):")
+		close(s.chRestoreLRO)
+	}()
+
+	// ensure all the goroutines are done
+	wgBuf.Wait()
+	
+	res := iop.NewPolynomial(&cres, iop.Form{Basis: iop.LagrangeCoset, Layout: iop.BitReverse})
+	log.Debug().Dur("took", time.Since(start)).Msg("computeNumerator")
+	return res, nil
+
+}
+
+func calculateNbTasks(n int) int {
+	nbAvailableCPU := runtime.NumCPU() - n
+	if nbAvailableCPU < 0 {
+		nbAvailableCPU = 1
+	}
+	nbTasks := 1 + (nbAvailableCPU / n)
+	return nbTasks
+}
+
+// batchApply executes fn on all polynomials in x except x[id_ZS] in parallel.
+func batchApply(x []*iop.Polynomial, fn func(*iop.Polynomial, int)) {
+	//var wg sync.WaitGroup
+	for i := 0; i < len(x); i++ {
+		if i == id_ZS {
+			continue
+		}
+			fn(x[i], i)
+		}
+}
+
+// p <- <p, (1, w, .., wⁿ) >
+// p is supposed to be in canonical form
+func scalePowers(p *iop.Polynomial, w fr.Element) {
+	var acc fr.Element
+	acc.SetOne()
+	cp := p.Coefficients()
+	for i := 0; i < p.Size(); i++ {
+		cp[i].Mul(&cp[i], &acc)
+		acc.Mul(&acc, &w)
+	}
+}
+
+func evaluateBlinded(p, bp *iop.Polynomial, zeta fr.Element) fr.Element {
+	// Get the size of the polynomial
+	n := big.NewInt(int64(p.Size()))
+
+	var pEvaluatedAtZeta fr.Element
+
+	// Evaluate the polynomial and blinded polynomial at zeta
+	chP := make(chan struct{}, 1)
+	go func() {
+		pEvaluatedAtZeta = p.Evaluate(zeta)
+		close(chP)
+	}()
+
+	bpEvaluatedAtZeta := bp.Evaluate(zeta)
+
+	// Multiply the evaluated blinded polynomial by tempElement
+	var t fr.Element
+	one := fr.One()
+	t.Exp(zeta, n).Sub(&t, &one)
+	bpEvaluatedAtZeta.Mul(&bpEvaluatedAtZeta, &t)
+
+	// Add the evaluated polynomial and the evaluated blinded polynomial
+	<-chP
+	pEvaluatedAtZeta.Add(&pEvaluatedAtZeta, &bpEvaluatedAtZeta)
+
+	// Return the result
+	return pEvaluatedAtZeta
+}
+
+// /!\ modifies p's underlying array of coefficients, in particular the size changes
+func getBlindedCoefficients(p, bp *iop.Polynomial) []fr.Element {
+	cp := p.Coefficients()
+	cbp := bp.Coefficients()
+	cp = append(cp, cbp...)
+	for i := 0; i < len(cbp); i++ {
+		cp[i].Sub(&cp[i], &cbp[i])
+	}
+	return cp
+}
+
+func deviceCommitBlindingFactor(n int, b *iop.Polynomial, pk *ProvingKey) curve.G1Affine {
+	// scalars
+	cp := b.Coefficients()
+	np := b.Size()
+
+	// get slice of points from unsafe pointer
+	resPtr := unsafe.Pointer(uintptr(unsafe.Pointer(pk.deviceInfo.G1Device.G1)) + uintptr(n)*unsafe.Sizeof(curve.G1Affine{}))
+
+	// Initialize channels
+	copyCpDone := make(chan unsafe.Pointer, 1)
+	cpDeviceData := make(chan iciclegnark.OnDeviceData, 1)
+
+	// Start asynchronous routine
+	go func() {
+
+		// tbd mul * 2
+		sizeBytesScalars := np * fr.Bytes
+
+		// Perform copy operation
+		iciclegnark.CopyToDevice(cp, sizeBytesScalars, copyCpDone)
+
+		// Receive result once copy operation is done
+		cpDevice := <-copyCpDone
+
+		// Create OnDeviceData
+		cpDeviceValue := iciclegnark.OnDeviceData{
+			P:    cpDevice,
+			Size: sizeBytesScalars,
+		}
+
+		// Send OnDeviceData to respective channel
+		cpDeviceData <- cpDeviceValue
+
+		// Close channels
+		close(copyCpDone)
+		close(cpDeviceData)
+	}()
+
+	// Wait for copy operation to finish
+	cpDeviceValue := <-cpDeviceData
+
+	var wg sync.WaitGroup
+
+	// Calculate(lo) commitment on Device
+	wg.Add(1)
+	tmpChan := make(chan curve.G1Affine, 1)
+	go func() {
+		defer wg.Done()
+		tmpVal, _, err := iciclegnark.MsmOnDevice(cpDeviceValue.P, pk.deviceInfo.G1Device.G1, np, true)
+		if err != nil {
+			fmt.Print("Error", err)
+		}
+		var tmp curve.G1Affine
+		tmp.FromJacobian(&tmpVal)
+		tmpChan <- tmp
+	}()
+	wg.Wait()
+
+	tmpAffinePoint := <-tmpChan
+
+	// Calculate(hi) commitment on Device
+	wg.Add(1)
+	resChan := make(chan curve.G1Affine, 1)
+	go func() {
+		defer wg.Done()
+
+		resVal, _, err := iciclegnark.MsmOnDevice(cpDeviceValue.P, resPtr, np, true)
+		if err != nil {
+			fmt.Print("Error", err)
+		}
+		var res curve.G1Affine
+		res.FromJacobian(&resVal)
+
+		resChan <- res
+	}()
+	wg.Wait()
+	resAffinePoint := <-resChan
+
+	// Sub(lo, hi) to get the final commitment
+	resAffinePoint.Sub(&resAffinePoint, &tmpAffinePoint)
+
+	// Free device memory
+	go func() {
+		iciclegnark.FreeDevicePointer(unsafe.Pointer(&cpDeviceValue))
+	}()
+
+	return resAffinePoint
+}
+
+// return a random polynomial of degree n, if n==-1 cancel the blinding
+func getRandomPolynomial(n int) *iop.Polynomial {
+	var a []fr.Element
+	if n == -1 {
+		a := make([]fr.Element, 1)
+		a[0].SetZero()
+	} else {
+		a = make([]fr.Element, n+1)
+		for i := 0; i <= n; i++ {
+			a[i].SetRandom()
+		}
+	}
+	res := iop.NewPolynomial(&a, iop.Form{
+		Basis: iop.Canonical, Layout: iop.Regular})
+	return res
+}
+
+func coefficients(p []*iop.Polynomial) [][]fr.Element {
+	res := make([][]fr.Element, len(p))
+	for i, pI := range p {
+		res[i] = pI.Coefficients()
+	}
+	return res
+}
+
+func commitToQuotient(h1, h2, h3 []fr.Element, proof *Proof, pk *ProvingKey) error {
+	log := logger.Logger()
+	start := time.Now()
+
+	g := new(errgroup.Group)
+
+	g.Go(func() (err error) {
+		start := time.Now()
+		proof.H[0], err = kzg.OnDeviceCommit(h1, pk.deviceInfo.G1Device.G1)
+		log.Debug().Dur("took", time.Since(start)).Int("size", len(h3)).Msg("MSM (commitToQuotient):")
+		return
+	})
+	g.Wait()
+
+	g.Go(func() (err error) {
+		start := time.Now()
+		proof.H[1], err = kzg.OnDeviceCommit(h2, pk.deviceInfo.G1Device.G1)
+		log.Debug().Dur("took", time.Since(start)).Int("size", len(h3)).Msg("MSM (commitToQuotient):")
+		return
+	})
+	g.Wait()
+
+	g.Go(func() (err error) {
+		start := time.Now()
+		proof.H[2], err = kzg.OnDeviceCommit(h3, pk.deviceInfo.G1Device.G1)
+		log.Debug().Dur("took", time.Since(start)).Int("size", len(h3)).Msg("MSM (commitToQuotient):")
+		return
+	})
+	g.Wait()
+
+	log.Debug().Dur("took", time.Since(start)).Msg("commitToQuotient")
+	return g.Wait()
+}
+
+// divideByXMinusOne
+// The input must be in LagrangeCoset.
+// The result is in Canonical Regular. (in place using a)
+func divideByXMinusOne(a *iop.Polynomial, domains [2]*fft.Domain) (*iop.Polynomial, error) {
+	log := logger.Logger()
+	start := time.Now()
+	// check that the basis is LagrangeCoset
+	if a.Basis != iop.LagrangeCoset || a.Layout != iop.BitReverse {
+		return nil, errors.New("invalid form")
+	}
+
+	// prepare the evaluations of x^n-1 on the big domain's coset
+	xnMinusOneInverseLagrangeCoset := evaluateXnMinusOneDomainBigCoset(domains)
+	rho := int(domains[1].Cardinality / domains[0].Cardinality)
+
+	r := a.Coefficients()
+	n := uint64(len(r))
+	nn := uint64(64 - bits.TrailingZeros64(n))
+
+	utils.Parallelize(len(r), func(start, end int) {
+		for i := start; i < end; i++ {
+			iRev := bits.Reverse64(uint64(i)) >> nn
+			r[i].Mul(&r[i], &xnMinusOneInverseLagrangeCoset[int(iRev)%rho])
+		}
+	})
+
+	// since a is in bit reverse order, ToRegular shouldn't do anything
+	a.ToCanonical(domains[1]).ToRegular()
+
+	log.Debug().Dur("took", time.Since(start)).Msg("FFT (divideByXMinusOne):")
+	return a, nil
+
+}
+
+// evaluateXnMinusOneDomainBigCoset evaluates Xᵐ-1 on DomainBig coset
+func evaluateXnMinusOneDomainBigCoset(domains [2]*fft.Domain) []fr.Element {
+
+	rho := domains[1].Cardinality / domains[0].Cardinality
+
+	res := make([]fr.Element, rho)
+
+	expo := big.NewInt(int64(domains[0].Cardinality))
+	res[0].Exp(domains[1].FrMultiplicativeGen, expo)
+
+	var t fr.Element
+	t.Exp(domains[1].Generator, big.NewInt(int64(domains[0].Cardinality)))
+
+	one := fr.One()
+
+	for i := 1; i < int(rho); i++ {
+		res[i].Mul(&res[i-1], &t)
+		res[i-1].Sub(&res[i-1], &one)
+	}
+	res[len(res)-1].Sub(&res[len(res)-1], &one)
+
+	res = fr.BatchInvert(res)
+
+	return res
+}
+
+// computeLinearizedPolynomial computes the linearized polynomial in canonical basis.
+// The purpose is to commit and open all in one ql, qr, qm, qo, qk.
+// * lZeta, rZeta, oZeta are the evaluation of l, r, o at zeta
+// * z is the permutation polynomial, zu is Z(μX), the shifted version of Z
+// * pk is the proving key: the linearized polynomial is a linear combination of ql, qr, qm, qo, qk.
+//
+// The Linearized polynomial is:
+//
+// α²*L₁(ζ)*Z(X)
+// + α*( (l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*Z(μζ)*s3(X) - Z(X)*(l(ζ)+β*id1(ζ)+γ)*(r(ζ)+β*id2(ζ)+γ)*(o(ζ)+β*id3(ζ)+γ))
+// + l(ζ)*Ql(X) + l(ζ)r(ζ)*Qm(X) + r(ζ)*Qr(X) + o(ζ)*Qo(X) + Qk(X)
+func computeLinearizedPolynomial(lZeta, rZeta, oZeta, alpha, beta, gamma, zeta, zu fr.Element, qcpZeta, blindedZCanonical []fr.Element, pi2Canonical [][]fr.Element, pk *ProvingKey) []fr.Element {
+
+	// first part: individual constraints
+	var rl fr.Element
+	rl.Mul(&rZeta, &lZeta)
+
+	// second part:
+	// Z(μζ)(l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*β*s3(X)-Z(X)(l(ζ)+β*id1(ζ)+γ)*(r(ζ)+β*id2(ζ)+γ)*(o(ζ)+β*id3(ζ)+γ)
+	var s1, s2 fr.Element
+	chS1 := make(chan struct{}, 1)
+	go func() {
+		s1 = pk.trace.S1.Evaluate(zeta)                      // s1(ζ)
+		s1.Mul(&s1, &beta).Add(&s1, &lZeta).Add(&s1, &gamma) // (l(ζ)+β*s1(ζ)+γ)
+		close(chS1)
+	}()
+	// ps2 := iop.NewPolynomial(&pk.S2Canonical, iop.Form{Basis: iop.Canonical, Layout: iop.Regular})
+	tmp := pk.trace.S2.Evaluate(zeta)                        // s2(ζ)
+	tmp.Mul(&tmp, &beta).Add(&tmp, &rZeta).Add(&tmp, &gamma) // (r(ζ)+β*s2(ζ)+γ)
+	<-chS1
+	s1.Mul(&s1, &tmp).Mul(&s1, &zu).Mul(&s1, &beta) // (l(ζ)+β*s1(β)+γ)*(r(ζ)+β*s2(β)+γ)*β*Z(μζ)
+
+	var uzeta, uuzeta fr.Element
+	uzeta.Mul(&zeta, &pk.Vk.CosetShift)
+	uuzeta.Mul(&uzeta, &pk.Vk.CosetShift)
+
+	s2.Mul(&beta, &zeta).Add(&s2, &lZeta).Add(&s2, &gamma)      // (l(ζ)+β*ζ+γ)
+	tmp.Mul(&beta, &uzeta).Add(&tmp, &rZeta).Add(&tmp, &gamma)  // (r(ζ)+β*u*ζ+γ)
+	s2.Mul(&s2, &tmp)                                           // (l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)
+	tmp.Mul(&beta, &uuzeta).Add(&tmp, &oZeta).Add(&tmp, &gamma) // (o(ζ)+β*u²*ζ+γ)
+	s2.Mul(&s2, &tmp)                                           // (l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ)
+	s2.Neg(&s2)                                                 // -(l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ)
+
+	// third part L₁(ζ)*α²*Z
+	var lagrangeZeta, one, den, frNbElmt fr.Element
+	one.SetOne()
+	nbElmt := int64(pk.Domain[0].Cardinality)
+	lagrangeZeta.Set(&zeta).
+		Exp(lagrangeZeta, big.NewInt(nbElmt)).
+		Sub(&lagrangeZeta, &one)
+	frNbElmt.SetUint64(uint64(nbElmt))
+	den.Sub(&zeta, &one).
+		Inverse(&den)
+	lagrangeZeta.Mul(&lagrangeZeta, &den). // L₁ = (ζⁿ⁻¹)/(ζ-1)
+						Mul(&lagrangeZeta, &alpha).
+						Mul(&lagrangeZeta, &alpha).
+						Mul(&lagrangeZeta, &pk.Domain[0].CardinalityInv) // (1/n)*α²*L₁(ζ)
+
+	s3canonical := pk.trace.S3.Coefficients()
+
+	utils.Parallelize(len(blindedZCanonical), func(start, end int) {
+
+		cql := pk.trace.Ql.Coefficients()
+		cqr := pk.trace.Qr.Coefficients()
+		cqm := pk.trace.Qm.Coefficients()
+		cqo := pk.trace.Qo.Coefficients()
+		cqk := pk.trace.Qk.Coefficients()
+
+		var t, t0, t1 fr.Element
+
+		for i := start; i < end; i++ {
+
+			t.Mul(&blindedZCanonical[i], &s2) // -Z(X)*(l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ)
+
+			if i < len(s3canonical) {
+
+				t0.Mul(&s3canonical[i], &s1) // (l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*Z(μζ)*β*s3(X)
+
+				t.Add(&t, &t0)
+			}
+
+			t.Mul(&t, &alpha) // α*( (l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*Z(μζ)*s3(X) - Z(X)*(l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ))
+
+			if i < len(cqm) {
+
+				t1.Mul(&cqm[i], &rl) // linPol = linPol + l(ζ)r(ζ)*Qm(X)
+
+				t0.Mul(&cql[i], &lZeta)
+				t0.Add(&t0, &t1)
+
+				t.Add(&t, &t0) // linPol = linPol + l(ζ)*Ql(X)
+
+				t0.Mul(&cqr[i], &rZeta)
+				t.Add(&t, &t0) // linPol = linPol + r(ζ)*Qr(X)
+
+				t0.Mul(&cqo[i], &oZeta)
+				t0.Add(&t0, &cqk[i])
+
+				t.Add(&t, &t0) // linPol = linPol + o(ζ)*Qo(X) + Qk(X)
+
+				for j := range qcpZeta {
+					t0.Mul(&pi2Canonical[j][i], &qcpZeta[j])
+					t.Add(&t, &t0)
+				}
+			}
+
+			t0.Mul(&blindedZCanonical[i], &lagrangeZeta)
+			blindedZCanonical[i].Add(&t, &t0) // finish the computation
+		}
+	})
+	return blindedZCanonical
+}
+
+var errContextDone = errors.New("context done")
diff --git a/backend/plonk/bn254/icicle/marshal.go b/backend/plonk/bn254/icicle/marshal.go
new file mode 100644
index 0000000000..8912b29d1e
--- /dev/null
+++ b/backend/plonk/bn254/icicle/marshal.go
@@ -0,0 +1,406 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by gnark DO NOT EDIT
+
+package icicle_bn254
+
+import (
+	curve "github.com/consensys/gnark-crypto/ecc/bn254"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+
+	"errors"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/iop"
+	"github.com/consensys/gnark-crypto/ecc/bn254/kzg"
+	"io"
+)
+
+// WriteRawTo writes binary encoding of Proof to w without point compression
+func (proof *Proof) WriteRawTo(w io.Writer) (int64, error) {
+	return proof.writeTo(w, curve.RawEncoding())
+}
+
+// WriteTo writes binary encoding of Proof to w with point compression
+func (proof *Proof) WriteTo(w io.Writer) (int64, error) {
+	return proof.writeTo(w)
+}
+
+func (proof *Proof) writeTo(w io.Writer, options ...func(*curve.Encoder)) (int64, error) {
+	enc := curve.NewEncoder(w, options...)
+
+	toEncode := []interface{}{
+		&proof.LRO[0],
+		&proof.LRO[1],
+		&proof.LRO[2],
+		&proof.Z,
+		&proof.H[0],
+		&proof.H[1],
+		&proof.H[2],
+		&proof.BatchedProof.H,
+		proof.BatchedProof.ClaimedValues,
+		&proof.ZShiftedOpening.H,
+		&proof.ZShiftedOpening.ClaimedValue,
+		proof.Bsb22Commitments,
+	}
+
+	for _, v := range toEncode {
+		if err := enc.Encode(v); err != nil {
+			return enc.BytesWritten(), err
+		}
+	}
+
+	return enc.BytesWritten(), nil
+}
+
+// ReadFrom reads binary representation of Proof from r
+func (proof *Proof) ReadFrom(r io.Reader) (int64, error) {
+	dec := curve.NewDecoder(r)
+	toDecode := []interface{}{
+		&proof.LRO[0],
+		&proof.LRO[1],
+		&proof.LRO[2],
+		&proof.Z,
+		&proof.H[0],
+		&proof.H[1],
+		&proof.H[2],
+		&proof.BatchedProof.H,
+		&proof.BatchedProof.ClaimedValues,
+		&proof.ZShiftedOpening.H,
+		&proof.ZShiftedOpening.ClaimedValue,
+		&proof.Bsb22Commitments,
+	}
+
+	for _, v := range toDecode {
+		if err := dec.Decode(v); err != nil {
+			return dec.BytesRead(), err
+		}
+	}
+
+	if proof.Bsb22Commitments == nil {
+		proof.Bsb22Commitments = []kzg.Digest{}
+	}
+
+	return dec.BytesRead(), nil
+}
+
+// WriteTo writes binary encoding of ProvingKey to w
+func (pk *ProvingKey) WriteTo(w io.Writer) (n int64, err error) {
+	return pk.writeTo(w, true)
+}
+
+// WriteRawTo writes binary encoding of ProvingKey to w without point compression
+func (pk *ProvingKey) WriteRawTo(w io.Writer) (n int64, err error) {
+	return pk.writeTo(w, false)
+}
+
+func (pk *ProvingKey) writeTo(w io.Writer, withCompression bool) (n int64, err error) {
+	// encode the verifying key
+	if withCompression {
+		n, err = pk.Vk.WriteTo(w)
+	} else {
+		n, err = pk.Vk.WriteRawTo(w)
+	}
+	if err != nil {
+		return
+	}
+
+	// fft domains
+	n2, err := pk.Domain[0].WriteTo(w)
+	if err != nil {
+		return
+	}
+	n += n2
+
+	n2, err = pk.Domain[1].WriteTo(w)
+	if err != nil {
+		return
+	}
+	n += n2
+
+	// KZG key
+	if withCompression {
+		n2, err = pk.Kzg.WriteTo(w)
+	} else {
+		n2, err = pk.Kzg.WriteRawTo(w)
+	}
+	if err != nil {
+		return
+	}
+	n += n2
+	if withCompression {
+		n2, err = pk.KzgLagrange.WriteTo(w)
+	} else {
+		n2, err = pk.KzgLagrange.WriteRawTo(w)
+	}
+	if err != nil {
+		return
+	}
+	n += n2
+
+	
+	// sanity check len(Permutation) == 3*int(pk.Domain[0].Cardinality)
+	if len(pk.trace.S) != (3 * int(pk.Domain[0].Cardinality)) {
+		return n, errors.New("invalid permutation size, expected 3*domain cardinality")
+	}
+
+	enc := curve.NewEncoder(w)
+	// note: type Polynomial, which is handled by default binary.Write(...) op and doesn't
+	// encode the size (nor does it convert from Montgomery to Regular form)
+	// so we explicitly transmit []fr.Element
+	toEncode := []interface{}{
+		pk.trace.Ql.Coefficients(),
+		pk.trace.Qr.Coefficients(),
+		pk.trace.Qm.Coefficients(),
+		pk.trace.Qo.Coefficients(),
+		pk.trace.Qk.Coefficients(),
+		coefficients(pk.trace.Qcp),
+		pk.trace.S1.Coefficients(),
+		pk.trace.S2.Coefficients(),
+		pk.trace.S3.Coefficients(),
+		pk.trace.S,
+	}
+
+	for _, v := range toEncode {
+		if err := enc.Encode(v); err != nil {
+			return n + enc.BytesWritten(), err
+		}
+	}
+
+	return n + enc.BytesWritten(), nil
+}
+
+// ReadFrom reads from binary representation in r into ProvingKey
+func (pk *ProvingKey) ReadFrom(r io.Reader) (int64, error) {
+	return pk.readFrom(r, true)
+}
+
+// UnsafeReadFrom reads from binary representation in r into ProvingKey without subgroup checks
+func (pk *ProvingKey) UnsafeReadFrom(r io.Reader) (int64, error) {
+	return pk.readFrom(r, false)
+}
+
+func (pk *ProvingKey) readFrom(r io.Reader, withSubgroupChecks bool) (int64, error) {
+	pk.Vk = &VerifyingKey{}
+	n, err := pk.Vk.ReadFrom(r)
+	if err != nil {
+		return n, err
+	}
+
+	n2, err, chDomain0 := pk.Domain[0].AsyncReadFrom(r)
+	n += n2
+	if err != nil {
+		return n, err
+	}
+
+	n2, err, chDomain1 := pk.Domain[1].AsyncReadFrom(r)
+	n += n2
+	if err != nil {
+		return n, err
+	}
+
+	if withSubgroupChecks {
+		n2, err = pk.Kzg.ReadFrom(r)
+	} else {
+		n2, err = pk.Kzg.UnsafeReadFrom(r)
+	}
+	n += n2
+	if err != nil {
+		return n, err
+	}
+	if withSubgroupChecks {
+		n2, err = pk.KzgLagrange.ReadFrom(r)
+	} else {
+		n2, err = pk.KzgLagrange.UnsafeReadFrom(r)
+	}
+	n += n2
+	if err != nil {
+		return n, err
+	}
+
+	pk.trace.S = make([]int64, 3*pk.Domain[0].Cardinality)
+
+	dec := curve.NewDecoder(r)
+
+	var ql, qr, qm, qo, qk, s1, s2, s3 []fr.Element
+	var qcp [][]fr.Element
+
+	// TODO @gbotrel: this is a bit ugly, we should probably refactor this.
+	// The order of the variables is important, as it matches the order in which they are
+	// encoded in the WriteTo(...) method.
+
+	// Note: instead of calling dec.Decode(...) for each of the above variables,
+	// we call AsyncReadFrom when possible which allows to consume bytes from the reader
+	// and perform the decoding in parallel
+
+	type v struct {
+		data  *fr.Vector
+		chErr chan error
+	}
+
+	vectors := make([]v, 8)
+	vectors[0] = v{data: (*fr.Vector)(&ql)}
+	vectors[1] = v{data: (*fr.Vector)(&qr)}
+	vectors[2] = v{data: (*fr.Vector)(&qm)}
+	vectors[3] = v{data: (*fr.Vector)(&qo)}
+	vectors[4] = v{data: (*fr.Vector)(&qk)}
+	vectors[5] = v{data: (*fr.Vector)(&s1)}
+	vectors[6] = v{data: (*fr.Vector)(&s2)}
+	vectors[7] = v{data: (*fr.Vector)(&s3)}
+
+	// read ql, qr, qm, qo, qk
+	for i := 0; i < 5; i++ {
+		n2, err, ch := vectors[i].data.AsyncReadFrom(r)
+		n += n2
+		if err != nil {
+			return n, err
+		}
+		vectors[i].chErr = ch
+	}
+
+	// read qcp
+	if err := dec.Decode(&qcp); err != nil {
+		return n + dec.BytesRead(), err
+	}
+
+	// read lqk, s1, s2, s3
+	for i := 5; i < 8; i++ {
+		n2, err, ch := vectors[i].data.AsyncReadFrom(r)
+		n += n2
+		if err != nil {
+			return n, err
+		}
+		vectors[i].chErr = ch
+	}
+
+	// read pk.Trace.S
+	if err := dec.Decode(&pk.trace.S); err != nil {
+		return n + dec.BytesRead(), err
+	}
+
+	// wait for all AsyncReadFrom(...) to complete
+	for i := range vectors {
+		if err := <-vectors[i].chErr; err != nil {
+			return n, err
+		}
+	}
+
+	canReg := iop.Form{Basis: iop.Canonical, Layout: iop.Regular}
+	pk.trace.Ql = iop.NewPolynomial(&ql, canReg)
+	pk.trace.Qr = iop.NewPolynomial(&qr, canReg)
+	pk.trace.Qm = iop.NewPolynomial(&qm, canReg)
+	pk.trace.Qo = iop.NewPolynomial(&qo, canReg)
+	pk.trace.Qk = iop.NewPolynomial(&qk, canReg)
+	pk.trace.S1 = iop.NewPolynomial(&s1, canReg)
+	pk.trace.S2 = iop.NewPolynomial(&s2, canReg)
+	pk.trace.S3 = iop.NewPolynomial(&s3, canReg)
+
+	pk.trace.Qcp = make([]*iop.Polynomial, len(qcp))
+	for i := range qcp {
+		pk.trace.Qcp[i] = iop.NewPolynomial(&qcp[i], canReg)
+	}
+
+	// wait for FFT to be precomputed
+	<-chDomain0
+	<-chDomain1
+
+	return n + dec.BytesRead(), nil
+
+}
+
+// WriteTo writes binary encoding of VerifyingKey to w
+func (vk *VerifyingKey) WriteTo(w io.Writer) (n int64, err error) {
+	return vk.writeTo(w)
+}
+
+// WriteRawTo writes binary encoding of VerifyingKey to w without point compression
+func (vk *VerifyingKey) WriteRawTo(w io.Writer) (int64, error) {
+	return vk.writeTo(w, curve.RawEncoding())
+}
+
+func (vk *VerifyingKey) writeTo(w io.Writer, options ...func(*curve.Encoder)) (n int64, err error) {
+	enc := curve.NewEncoder(w)
+
+	toEncode := []interface{}{
+		vk.Size,
+		&vk.SizeInv,
+		&vk.Generator,
+		vk.NbPublicVariables,
+		&vk.CosetShift,
+		&vk.S[0],
+		&vk.S[1],
+		&vk.S[2],
+		&vk.Ql,
+		&vk.Qr,
+		&vk.Qm,
+		&vk.Qo,
+		&vk.Qk,
+		vk.Qcp,
+		&vk.Kzg.G1,
+		&vk.Kzg.G2[0],
+		&vk.Kzg.G2[1],
+		vk.CommitmentConstraintIndexes,
+	}
+
+	for _, v := range toEncode {
+		if err := enc.Encode(v); err != nil {
+			return enc.BytesWritten(), err
+		}
+	}
+
+	return enc.BytesWritten(), nil
+}
+
+// UnsafeReadFrom reads from binary representation in r into VerifyingKey.
+// Current implementation is a passthrough to ReadFrom
+func (vk *VerifyingKey) UnsafeReadFrom(r io.Reader) (int64, error) {
+	return vk.ReadFrom(r)
+}
+
+// ReadFrom reads from binary representation in r into VerifyingKey
+func (vk *VerifyingKey) ReadFrom(r io.Reader) (int64, error) {
+	dec := curve.NewDecoder(r)
+	toDecode := []interface{}{
+		&vk.Size,
+		&vk.SizeInv,
+		&vk.Generator,
+		&vk.NbPublicVariables,
+		&vk.CosetShift,
+		&vk.S[0],
+		&vk.S[1],
+		&vk.S[2],
+		&vk.Ql,
+		&vk.Qr,
+		&vk.Qm,
+		&vk.Qo,
+		&vk.Qk,
+		&vk.Qcp,
+		&vk.Kzg.G1,
+		&vk.Kzg.G2[0],
+		&vk.Kzg.G2[1],
+		&vk.CommitmentConstraintIndexes,
+	}
+
+	for _, v := range toDecode {
+		if err := dec.Decode(v); err != nil {
+			return dec.BytesRead(), err
+		}
+	}
+
+	if vk.Qcp == nil {
+		vk.Qcp = []kzg.Digest{}
+	}
+
+	return dec.BytesRead(), nil
+}
diff --git a/backend/plonk/bn254/icicle/marshal_test.go b/backend/plonk/bn254/icicle/marshal_test.go
new file mode 100644
index 0000000000..349ab38887
--- /dev/null
+++ b/backend/plonk/bn254/icicle/marshal_test.go
@@ -0,0 +1,175 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by gnark DO NOT EDIT
+
+package icicle_bn254
+
+import (
+	curve "github.com/consensys/gnark-crypto/ecc/bn254"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/fft"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/iop"
+	"github.com/consensys/gnark/io"
+	"math/big"
+	"math/rand"
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestProofSerialization(t *testing.T) {
+	// create a  proof
+	var proof Proof
+	proof.randomize()
+
+	assert.NoError(t, io.RoundTripCheck(&proof, func() interface{} { return new(Proof) }))
+}
+
+func TestProvingKeySerialization(t *testing.T) {
+	// random pk
+	var pk ProvingKey
+	pk.randomize()
+
+	assert.NoError(t, io.RoundTripCheck(&pk, func() interface{} { return new(ProvingKey) }))
+}
+
+func TestVerifyingKeySerialization(t *testing.T) {
+	// create a random vk
+	var vk VerifyingKey
+	vk.randomize()
+
+	assert.NoError(t, io.RoundTripCheck(&vk, func() interface{} { return new(VerifyingKey) }))
+}
+
+func (pk *ProvingKey) randomize() {
+
+	var vk VerifyingKey
+	vk.randomize()
+	pk.Vk = &vk
+	pk.Domain[0] = *fft.NewDomain(32)
+	pk.Domain[1] = *fft.NewDomain(4 * 32)
+
+	pk.Kzg.G1 = make([]curve.G1Affine, 32)
+	pk.KzgLagrange.G1 = make([]curve.G1Affine, 32)
+	for i := range pk.Kzg.G1 {
+		pk.Kzg.G1[i] = randomG1Point()
+		pk.KzgLagrange.G1[i] = randomG1Point()
+	}
+
+	n := int(pk.Domain[0].Cardinality)
+	ql := randomScalars(n)
+	qr := randomScalars(n)
+	qm := randomScalars(n)
+	qo := randomScalars(n)
+	qk := randomScalars(n)
+	s1 := randomScalars(n)
+	s2 := randomScalars(n)
+	s3 := randomScalars(n)
+
+	canReg := iop.Form{Basis: iop.Canonical, Layout: iop.Regular}
+	pk.trace.Ql = iop.NewPolynomial(&ql, canReg)
+	pk.trace.Qr = iop.NewPolynomial(&qr, canReg)
+	pk.trace.Qm = iop.NewPolynomial(&qm, canReg)
+	pk.trace.Qo = iop.NewPolynomial(&qo, canReg)
+	pk.trace.Qk = iop.NewPolynomial(&qk, canReg)
+	pk.trace.S1 = iop.NewPolynomial(&s1, canReg)
+	pk.trace.S2 = iop.NewPolynomial(&s2, canReg)
+	pk.trace.S3 = iop.NewPolynomial(&s3, canReg)
+
+	pk.trace.Qcp = make([]*iop.Polynomial, rand.Intn(4)) //#nosec G404 weak rng is fine here
+	for i := range pk.trace.Qcp {
+		qcp := randomScalars(rand.Intn(n / 4)) //#nosec G404 weak rng is fine here
+		pk.trace.Qcp[i] = iop.NewPolynomial(&qcp, canReg)
+	}
+
+	pk.trace.S = make([]int64, 3*pk.Domain[0].Cardinality)
+	pk.trace.S[0] = -12
+	pk.trace.S[len(pk.trace.S)-1] = 8888
+
+}
+
+func (vk *VerifyingKey) randomize() {
+	vk.Size = rand.Uint64() //#nosec G404 weak rng is fine here
+	vk.SizeInv.SetRandom()
+	vk.Generator.SetRandom()
+	vk.NbPublicVariables = rand.Uint64()                     //#nosec G404 weak rng is fine here
+	vk.CommitmentConstraintIndexes = []uint64{rand.Uint64()} //#nosec G404 weak rng is fine here
+	vk.CosetShift.SetRandom()
+
+	vk.S[0] = randomG1Point()
+	vk.S[1] = randomG1Point()
+	vk.S[2] = randomG1Point()
+
+	vk.Kzg.G1 = randomG1Point()
+	vk.Kzg.G2[0] = randomG2Point()
+	vk.Kzg.G2[1] = randomG2Point()
+
+	vk.Ql = randomG1Point()
+	vk.Qr = randomG1Point()
+	vk.Qm = randomG1Point()
+	vk.Qo = randomG1Point()
+	vk.Qk = randomG1Point()
+	vk.Qcp = randomG1Points(rand.Intn(4)) //#nosec G404 weak rng is fine here
+}
+
+func (proof *Proof) randomize() {
+	proof.LRO[0] = randomG1Point()
+	proof.LRO[1] = randomG1Point()
+	proof.LRO[2] = randomG1Point()
+	proof.Z = randomG1Point()
+	proof.H[0] = randomG1Point()
+	proof.H[1] = randomG1Point()
+	proof.H[2] = randomG1Point()
+	proof.BatchedProof.H = randomG1Point()
+	proof.BatchedProof.ClaimedValues = randomScalars(2)
+	proof.ZShiftedOpening.H = randomG1Point()
+	proof.ZShiftedOpening.ClaimedValue.SetRandom()
+	proof.Bsb22Commitments = randomG1Points(rand.Intn(4)) //#nosec G404 weak rng is fine here
+}
+
+func randomG2Point() curve.G2Affine {
+	_, _, _, r := curve.Generators()
+	r.ScalarMultiplication(&r, big.NewInt(int64(rand.Uint64()))) //#nosec G404 weak rng is fine here
+	return r
+}
+
+func randomG1Point() curve.G1Affine {
+	_, _, r, _ := curve.Generators()
+	r.ScalarMultiplication(&r, big.NewInt(int64(rand.Uint64()))) //#nosec G404 weak rng is fine here
+	return r
+}
+
+func randomG1Points(n int) []curve.G1Affine {
+	res := make([]curve.G1Affine, n)
+	for i := range res {
+		res[i] = randomG1Point()
+	}
+	return res
+}
+
+func randomScalars(n int) []fr.Element {
+	v := make([]fr.Element, n)
+	one := fr.One()
+	for i := 0; i < len(v); i++ {
+		if i == 0 {
+			v[i].SetRandom()
+		} else {
+			v[i].Add(&v[i-1], &one)
+		}
+	}
+	return v
+}
diff --git a/backend/plonk/bn254/icicle/setup.go b/backend/plonk/bn254/icicle/setup.go
new file mode 100644
index 0000000000..c46ba58794
--- /dev/null
+++ b/backend/plonk/bn254/icicle/setup.go
@@ -0,0 +1,410 @@
+package icicle_bn254
+
+import (
+	"errors"
+	"fmt"
+	"unsafe"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/fft"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/iop"
+	"github.com/consensys/gnark-crypto/ecc/bn254/kzg"
+	"github.com/consensys/gnark/backend/plonk/internal"
+	"github.com/consensys/gnark/constraint"
+	cs "github.com/consensys/gnark/constraint/bn254"
+)
+
+// VerifyingKey stores the data needed to verify a proof:
+// * The commitment scheme
+// * Commitments of ql prepended with as many ones as there are public inputs
+// * Commitments of qr, qm, qo, qk prepended with as many zeroes as there are public inputs
+// * Commitments to S1, S2, S3
+type VerifyingKey struct {
+	// Size circuit
+	Size              uint64
+	SizeInv           fr.Element
+	Generator         fr.Element
+	NbPublicVariables uint64
+
+	// Commitment scheme that is used for an instantiation of PLONK
+	Kzg kzg.VerifyingKey
+
+	// cosetShift generator of the coset on the small domain
+	CosetShift fr.Element
+
+	// S commitments to S1, S2, S3
+	S [3]kzg.Digest
+
+	// Commitments to ql, qr, qm, qo, qcp prepended with as many zeroes (ones for l) as there are public inputs.
+	// In particular Qk is not complete.
+	Ql, Qr, Qm, Qo, Qk kzg.Digest
+	Qcp                []kzg.Digest
+
+	CommitmentConstraintIndexes []uint64
+}
+
+// Trace stores a plonk trace as columns
+type Trace struct {
+	// Constants describing a plonk circuit. The first entries
+	// of LQk (whose index correspond to the public inputs) are set to 0, and are to be
+	// completed by the prover. At those indices i (so from 0 to nb_public_variables), LQl[i]=-1
+	// so the first nb_public_variables constraints look like this:
+	// -1*Wire[i] + 0* + 0 . It is zero when the constant coefficient is replaced by Wire[i].
+	Ql, Qr, Qm, Qo, Qk *iop.Polynomial
+	Qcp                []*iop.Polynomial
+
+	// Polynomials representing the splitted permutation. The full permutation's support is 3*N where N=nb wires.
+	// The set of interpolation is <g> of size N, so to represent the permutation S we let S acts on the
+	// set A=(<g>, u*<g>, u^{2}*<g>) of size 3*N, where u is outside <g> (its use is to shift the set <g>).
+	// We obtain a permutation of A, A'. We split A' in 3 (A'_{1}, A'_{2}, A'_{3}), and S1, S2, S3 are
+	// respectively the interpolation of A'_{1}, A'_{2}, A'_{3} on <g>.
+	S1, S2, S3 *iop.Polynomial
+
+	// S full permutation, i -> S[i]
+	S []int64
+}
+
+// ProvingKey stores the data needed to generate a proof:
+// * the commitment scheme
+// * ql, prepended with as many ones as they are public inputs
+// * qr, qm, qo prepended with as many zeroes as there are public inputs.
+// * qk, prepended with as many zeroes as public inputs, to be completed by the prover
+// with the list of public inputs.
+// * sigma_1, sigma_2, sigma_3 in both basis
+// * the copy constraint permutation
+type ProvingKey struct {
+	// stores ql, qr, qm, qo, qk (-> to be completed by the prover)
+	// and s1, s2, s3. They are set in canonical basis before generating the proof, they will be used
+	// for computing the opening proofs (hence the canonical form). The canonical version
+	// of qk incomplete is used in the linearisation polynomial.
+	// The polynomials in trace are in canonical basis.
+	trace Trace
+
+	Kzg, KzgLagrange kzg.ProvingKey
+
+	// Verifying Key is embedded into the proving key (needed by Prove)
+	Vk *VerifyingKey
+
+	// Domains used for the FFTs.
+	// Domain[0] = small Domain
+	// Domain[1] = big Domain
+	Domain [2]fft.Domain
+
+	deviceInfo *deviceInfo
+}
+
+type deviceInfo struct {
+	G1Device struct {
+		G1         unsafe.Pointer
+		G1Lagrange unsafe.Pointer
+	}
+	DomainDevice struct {
+		Twiddles, TwiddlesInv     unsafe.Pointer
+		CosetTable, CosetTableInv unsafe.Pointer
+	}
+	DenDevice             unsafe.Pointer
+	InfinityPointIndicesK []int
+}
+
+// TODO modify the signature to receive the SRS in Lagrange form (optional argument ?)
+func Setup(spr *cs.SparseR1CS, kzgSrs kzg.SRS) (*ProvingKey, *VerifyingKey, error) {
+
+	var pk ProvingKey
+	var vk VerifyingKey
+	pk.Vk = &vk
+	vk.CommitmentConstraintIndexes = internal.IntSliceToUint64Slice(spr.CommitmentInfo.CommitmentIndexes())
+
+	pk.deviceInfo = nil
+
+	// step 0: set the fft domains
+	pk.initDomains(spr)
+	if pk.Domain[0].Cardinality < 2 {
+		return nil, nil, fmt.Errorf("circuit has only %d constraints; unsupported by the current implementation", spr.GetNbConstraints())
+	}
+
+	// step 1: set the verifying key
+	pk.Vk.CosetShift.Set(&pk.Domain[0].FrMultiplicativeGen)
+	vk.Size = pk.Domain[0].Cardinality
+	vk.SizeInv.SetUint64(vk.Size).Inverse(&vk.SizeInv)
+	vk.Generator.Set(&pk.Domain[0].Generator)
+	vk.NbPublicVariables = uint64(len(spr.Public))
+	if len(kzgSrs.Pk.G1) < int(vk.Size)+3 { // + 3 for the kzg.Open of blinded poly
+		return nil, nil, errors.New("kzg srs is too small")
+	}
+	pk.Kzg.G1 = kzgSrs.Pk.G1[:int(vk.Size)+3]
+	var err error
+	pk.KzgLagrange.G1, err = kzg.ToLagrangeG1(kzgSrs.Pk.G1[:int(vk.Size)])
+	if err != nil {
+		return nil, nil, err
+	}
+	vk.Kzg = kzgSrs.Vk
+
+	// step 2: ql, qr, qm, qo, qk, qcp in Lagrange Basis
+	BuildTrace(spr, &pk.trace)
+
+	// step 3: build the permutation and build the polynomials S1, S2, S3 to encode the permutation.
+	// Note: at this stage, the permutation takes in account the placeholders
+	nbVariables := spr.NbInternalVariables + len(spr.Public) + len(spr.Secret)
+	buildPermutation(spr, &pk.trace, nbVariables)
+	s := computePermutationPolynomials(&pk.trace, &pk.Domain[0])
+	pk.trace.S1 = s[0]
+	pk.trace.S2 = s[1]
+	pk.trace.S3 = s[2]
+
+	// step 4: commit to s1, s2, s3, ql, qr, qm, qo, and (the incomplete version of) qk.
+	// All the above polynomials are expressed in canonical basis afterwards. This is why
+	// we save lqk before, because the prover needs to complete it in Lagrange form, and
+	// then express it on the Lagrange coset basis.
+	if err = commitTrace(&pk.trace, &pk); err != nil {
+		return nil, nil, err
+	}
+
+	return &pk, &vk, nil
+}
+
+// NbPublicWitness returns the expected public witness size (number of field elements)
+func (vk *VerifyingKey) NbPublicWitness() int {
+	return int(vk.NbPublicVariables)
+}
+
+// VerifyingKey returns pk.Vk
+func (pk *ProvingKey) VerifyingKey() interface{} {
+	return pk.Vk
+}
+
+// BuildTrace fills the constant columns ql, qr, qm, qo, qk from the sparser1cs.
+// Size is the size of the system that is nb_constraints+nb_public_variables
+func BuildTrace(spr *cs.SparseR1CS, pt *Trace) {
+
+	nbConstraints := spr.GetNbConstraints()
+	sizeSystem := uint64(nbConstraints + len(spr.Public))
+	size := ecc.NextPowerOfTwo(sizeSystem)
+	commitmentInfo := spr.CommitmentInfo.(constraint.PlonkCommitments)
+
+	ql := make([]fr.Element, size)
+	qr := make([]fr.Element, size)
+	qm := make([]fr.Element, size)
+	qo := make([]fr.Element, size)
+	qk := make([]fr.Element, size)
+	qcp := make([][]fr.Element, len(commitmentInfo))
+
+	for i := 0; i < len(spr.Public); i++ { // placeholders (-PUB_INPUT_i + qk_i = 0) TODO should return error if size is inconsistent
+		ql[i].SetOne().Neg(&ql[i])
+		qr[i].SetZero()
+		qm[i].SetZero()
+		qo[i].SetZero()
+		qk[i].SetZero() // → to be completed by the prover
+	}
+	offset := len(spr.Public)
+
+	j := 0
+	it := spr.GetSparseR1CIterator()
+	for c := it.Next(); c != nil; c = it.Next() {
+		ql[offset+j].Set(&spr.Coefficients[c.QL])
+		qr[offset+j].Set(&spr.Coefficients[c.QR])
+		qm[offset+j].Set(&spr.Coefficients[c.QM])
+		qo[offset+j].Set(&spr.Coefficients[c.QO])
+		qk[offset+j].Set(&spr.Coefficients[c.QC])
+		j++
+	}
+
+	lagReg := iop.Form{Basis: iop.Lagrange, Layout: iop.Regular}
+
+	pt.Ql = iop.NewPolynomial(&ql, lagReg)
+	pt.Qr = iop.NewPolynomial(&qr, lagReg)
+	pt.Qm = iop.NewPolynomial(&qm, lagReg)
+	pt.Qo = iop.NewPolynomial(&qo, lagReg)
+	pt.Qk = iop.NewPolynomial(&qk, lagReg)
+	pt.Qcp = make([]*iop.Polynomial, len(qcp))
+
+	for i := range commitmentInfo {
+		qcp[i] = make([]fr.Element, size)
+		for _, committed := range commitmentInfo[i].Committed {
+			qcp[i][offset+committed].SetOne()
+		}
+		pt.Qcp[i] = iop.NewPolynomial(&qcp[i], lagReg)
+	}
+}
+
+// commitTrace commits to every polynomial in the trace, and put
+// the commitments int the verifying key.
+func commitTrace(trace *Trace, pk *ProvingKey) error {
+
+	trace.Ql.ToCanonical(&pk.Domain[0]).ToRegular()
+	trace.Qr.ToCanonical(&pk.Domain[0]).ToRegular()
+	trace.Qm.ToCanonical(&pk.Domain[0]).ToRegular()
+	trace.Qo.ToCanonical(&pk.Domain[0]).ToRegular()
+	trace.Qk.ToCanonical(&pk.Domain[0]).ToRegular() // -> qk is not complete
+	trace.S1.ToCanonical(&pk.Domain[0]).ToRegular()
+	trace.S2.ToCanonical(&pk.Domain[0]).ToRegular()
+	trace.S3.ToCanonical(&pk.Domain[0]).ToRegular()
+
+	var err error
+	pk.Vk.Qcp = make([]kzg.Digest, len(trace.Qcp))
+	for i := range trace.Qcp {
+		trace.Qcp[i].ToCanonical(&pk.Domain[0]).ToRegular()
+		if pk.Vk.Qcp[i], err = kzg.Commit(pk.trace.Qcp[i].Coefficients(), pk.Kzg); err != nil {
+			return err
+		}
+	}
+	if pk.Vk.Ql, err = kzg.Commit(pk.trace.Ql.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.Qr, err = kzg.Commit(pk.trace.Qr.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.Qm, err = kzg.Commit(pk.trace.Qm.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.Qo, err = kzg.Commit(pk.trace.Qo.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.Qk, err = kzg.Commit(pk.trace.Qk.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.S[0], err = kzg.Commit(pk.trace.S1.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.S[1], err = kzg.Commit(pk.trace.S2.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	if pk.Vk.S[2], err = kzg.Commit(pk.trace.S3.Coefficients(), pk.Kzg); err != nil {
+		return err
+	}
+	return nil
+}
+
+func (pk *ProvingKey) initDomains(spr *cs.SparseR1CS) {
+
+	nbConstraints := spr.GetNbConstraints()
+	sizeSystem := uint64(nbConstraints + len(spr.Public)) // len(spr.Public) is for the placeholder constraints
+	pk.Domain[0] = *fft.NewDomain(sizeSystem)
+
+	// h, the quotient polynomial is of degree 3(n+1)+2, so it's in a 3(n+2) dim vector space,
+	// the domain is the next power of 2 superior to 3(n+2). 4*domainNum is enough in all cases
+	// except when n<6.
+	if sizeSystem < 6 {
+		pk.Domain[1] = *fft.NewDomain(8 * sizeSystem)
+	} else {
+		pk.Domain[1] = *fft.NewDomain(4 * sizeSystem)
+	}
+
+}
+
+// buildPermutation builds the Permutation associated with a circuit.
+//
+// The permutation s is composed of cycles of maximum length such that
+//
+//	s. (l∥r∥o) = (l∥r∥o)
+//
+// , where l∥r∥o is the concatenation of the indices of l, r, o in
+// ql.l+qr.r+qm.l.r+qo.O+k = 0.
+//
+// The permutation is encoded as a slice s of size 3*size(l), where the
+// i-th entry of l∥r∥o is sent to the s[i]-th entry, so it acts on a tab
+// like this: for i in tab: tab[i] = tab[permutation[i]]
+func buildPermutation(spr *cs.SparseR1CS, pt *Trace, nbVariables int) {
+
+	// nbVariables := spr.NbInternalVariables + len(spr.Public) + len(spr.Secret)
+	sizeSolution := len(pt.Ql.Coefficients())
+	sizePermutation := 3 * sizeSolution
+
+	// init permutation
+	permutation := make([]int64, sizePermutation)
+	for i := 0; i < len(permutation); i++ {
+		permutation[i] = -1
+	}
+
+	// init LRO position -> variable_ID
+	lro := make([]int, sizePermutation) // position -> variable_ID
+	for i := 0; i < len(spr.Public); i++ {
+		lro[i] = i // IDs of LRO associated to placeholders (only L needs to be taken care of)
+	}
+
+	offset := len(spr.Public)
+
+	j := 0
+	it := spr.GetSparseR1CIterator()
+	for c := it.Next(); c != nil; c = it.Next() {
+		lro[offset+j] = int(c.XA)
+		lro[sizeSolution+offset+j] = int(c.XB)
+		lro[2*sizeSolution+offset+j] = int(c.XC)
+
+		j++
+	}
+
+	// init cycle:
+	// map ID -> last position the ID was seen
+	cycle := make([]int64, nbVariables)
+	for i := 0; i < len(cycle); i++ {
+		cycle[i] = -1
+	}
+
+	for i := 0; i < len(lro); i++ {
+		if cycle[lro[i]] != -1 {
+			// if != -1, it means we already encountered this value
+			// so we need to set the corresponding permutation index.
+			permutation[i] = cycle[lro[i]]
+		}
+		cycle[lro[i]] = int64(i)
+	}
+
+	// complete the Permutation by filling the first IDs encountered
+	for i := 0; i < sizePermutation; i++ {
+		if permutation[i] == -1 {
+			permutation[i] = cycle[lro[i]]
+		}
+	}
+
+	pt.S = permutation
+}
+
+// computePermutationPolynomials computes the LDE (Lagrange basis) of the permutation.
+// We let the permutation act on <g> || u<g> || u^{2}<g>, split the result in 3 parts,
+// and interpolate each of the 3 parts on <g>.
+func computePermutationPolynomials(pt *Trace, domain *fft.Domain) [3]*iop.Polynomial {
+
+	nbElmts := int(domain.Cardinality)
+
+	var res [3]*iop.Polynomial
+
+	// Lagrange form of ID
+	evaluationIDSmallDomain := getSupportPermutation(domain)
+
+	// Lagrange form of S1, S2, S3
+	s1Canonical := make([]fr.Element, nbElmts)
+	s2Canonical := make([]fr.Element, nbElmts)
+	s3Canonical := make([]fr.Element, nbElmts)
+	for i := 0; i < nbElmts; i++ {
+		s1Canonical[i].Set(&evaluationIDSmallDomain[pt.S[i]])
+		s2Canonical[i].Set(&evaluationIDSmallDomain[pt.S[nbElmts+i]])
+		s3Canonical[i].Set(&evaluationIDSmallDomain[pt.S[2*nbElmts+i]])
+	}
+
+	lagReg := iop.Form{Basis: iop.Lagrange, Layout: iop.Regular}
+	res[0] = iop.NewPolynomial(&s1Canonical, lagReg)
+	res[1] = iop.NewPolynomial(&s2Canonical, lagReg)
+	res[2] = iop.NewPolynomial(&s3Canonical, lagReg)
+
+	return res
+}
+
+// getSupportPermutation returns the support on which the permutation acts, it is
+// <g> || u<g> || u^{2}<g>
+func getSupportPermutation(domain *fft.Domain) []fr.Element {
+
+	res := make([]fr.Element, 3*domain.Cardinality)
+
+	res[0].SetOne()
+	res[domain.Cardinality].Set(&domain.FrMultiplicativeGen)
+	res[2*domain.Cardinality].Square(&domain.FrMultiplicativeGen)
+
+	for i := uint64(1); i < domain.Cardinality; i++ {
+		res[i].Mul(&res[i-1], &domain.Generator)
+		res[domain.Cardinality+i].Mul(&res[domain.Cardinality+i-1], &domain.Generator)
+		res[2*domain.Cardinality+i].Mul(&res[2*domain.Cardinality+i-1], &domain.Generator)
+	}
+
+	return res
+}
diff --git a/backend/plonk/bn254/icicle/solidity.go b/backend/plonk/bn254/icicle/solidity.go
new file mode 100644
index 0000000000..22b7f36454
--- /dev/null
+++ b/backend/plonk/bn254/icicle/solidity.go
@@ -0,0 +1,1420 @@
+package icicle_bn254
+
+const tmplSolidityVerifier = `// SPDX-License-Identifier: Apache-2.0
+
+// Copyright 2023 Consensys Software Inc.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by gnark DO NOT EDIT
+
+pragma solidity ^0.8.19;
+
+contract PlonkVerifier {
+
+  uint256 private constant R_MOD = 21888242871839275222246405745257275088548364400416034343698204186575808495617;
+  uint256 private constant P_MOD = 21888242871839275222246405745257275088696311157297823662689037894645226208583;
+  {{ range $index, $element := .Kzg.G2 }}
+  uint256 private constant G2_SRS_{{ $index }}_X_0 = {{ (fpstr $element.X.A1) }};
+  uint256 private constant G2_SRS_{{ $index }}_X_1 = {{ (fpstr $element.X.A0) }};
+  uint256 private constant G2_SRS_{{ $index }}_Y_0 = {{ (fpstr $element.Y.A1) }};
+  uint256 private constant G2_SRS_{{ $index }}_Y_1 = {{ (fpstr $element.Y.A0) }};
+  {{ end }}
+  uint256 private constant G1_SRS_X = {{ fpstr .Kzg.G1.X }};
+  uint256 private constant G1_SRS_Y = {{ fpstr .Kzg.G1.Y }};
+
+  // ----------------------- vk ---------------------
+  uint256 private constant VK_NB_PUBLIC_INPUTS = {{ .NbPublicVariables }};
+  uint256 private constant VK_DOMAIN_SIZE = {{ .Size }};
+  uint256 private constant VK_INV_DOMAIN_SIZE = {{ (frstr .SizeInv) }};
+  uint256 private constant VK_OMEGA = {{ (frstr .Generator) }};
+  uint256 private constant VK_QL_COM_X = {{ (fpstr .Ql.X) }};
+  uint256 private constant VK_QL_COM_Y = {{ (fpstr .Ql.Y) }};
+  uint256 private constant VK_QR_COM_X = {{ (fpstr .Qr.X) }};
+  uint256 private constant VK_QR_COM_Y = {{ (fpstr .Qr.Y) }};
+  uint256 private constant VK_QM_COM_X = {{ (fpstr .Qm.X) }};
+  uint256 private constant VK_QM_COM_Y = {{ (fpstr .Qm.Y) }};
+  uint256 private constant VK_QO_COM_X = {{ (fpstr .Qo.X) }};
+  uint256 private constant VK_QO_COM_Y = {{ (fpstr .Qo.Y) }};
+  uint256 private constant VK_QK_COM_X = {{ (fpstr .Qk.X) }};
+  uint256 private constant VK_QK_COM_Y = {{ (fpstr .Qk.Y) }};
+  {{ range $index, $element := .S }}
+  uint256 private constant VK_S{{ inc $index }}_COM_X = {{ (fpstr $element.X) }};
+  uint256 private constant VK_S{{ inc $index }}_COM_Y = {{ (fpstr $element.Y) }};
+  {{ end }}
+  uint256 private constant VK_COSET_SHIFT = 5;
+  
+  {{ range $index, $element := .Qcp}}
+  uint256 private constant VK_QCP_{{ $index }}_X = {{ (fpstr $element.X) }};
+  uint256 private constant VK_QCP_{{ $index }}_Y = {{ (fpstr $element.Y) }};
+  {{ end }}
+  
+  {{ range $index, $element := .CommitmentConstraintIndexes -}}
+  uint256 private constant VK_INDEX_COMMIT_API{{ $index }} = {{ $element }};
+  {{ end -}}
+  uint256 private constant VK_NB_CUSTOM_GATES = {{ len .CommitmentConstraintIndexes }};
+
+  // ------------------------------------------------
+
+  // offset proof
+  uint256 private constant PROOF_L_COM_X = 0x00;
+  uint256 private constant PROOF_L_COM_Y = 0x20;
+  uint256 private constant PROOF_R_COM_X = 0x40;
+  uint256 private constant PROOF_R_COM_Y = 0x60;
+  uint256 private constant PROOF_O_COM_X = 0x80;
+  uint256 private constant PROOF_O_COM_Y = 0xa0;
+
+  // h = h_0 + x^{n+2}h_1 + x^{2(n+2)}h_2
+  uint256 private constant PROOF_H_0_X = 0xc0;
+  uint256 private constant PROOF_H_0_Y = 0xe0;
+  uint256 private constant PROOF_H_1_X = 0x100;
+  uint256 private constant PROOF_H_1_Y = 0x120;
+  uint256 private constant PROOF_H_2_X = 0x140;
+  uint256 private constant PROOF_H_2_Y = 0x160;
+
+  // wire values at zeta
+  uint256 private constant PROOF_L_AT_ZETA = 0x180;
+  uint256 private constant PROOF_R_AT_ZETA = 0x1a0;
+  uint256 private constant PROOF_O_AT_ZETA = 0x1c0;
+
+  //uint256[STATE_WIDTH-1] permutation_polynomials_at_zeta; // Sσ1(zeta),Sσ2(zeta)
+  uint256 private constant PROOF_S1_AT_ZETA = 0x1e0; // Sσ1(zeta)
+  uint256 private constant PROOF_S2_AT_ZETA = 0x200; // Sσ2(zeta)
+
+  //Bn254.G1Point grand_product_commitment;                 // [z(x)]
+  uint256 private constant PROOF_GRAND_PRODUCT_COMMITMENT_X = 0x220;
+  uint256 private constant PROOF_GRAND_PRODUCT_COMMITMENT_Y = 0x240;
+
+  uint256 private constant PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA = 0x260; // z(w*zeta)
+  uint256 private constant PROOF_QUOTIENT_POLYNOMIAL_AT_ZETA = 0x280; // t(zeta)
+  uint256 private constant PROOF_LINEARISED_POLYNOMIAL_AT_ZETA = 0x2a0; // r(zeta)
+
+  // Folded proof for the opening of H, linearised poly, l, r, o, s_1, s_2, qcp
+  uint256 private constant PROOF_BATCH_OPENING_AT_ZETA_X = 0x2c0; // [Wzeta]
+  uint256 private constant PROOF_BATCH_OPENING_AT_ZETA_Y = 0x2e0;
+
+  uint256 private constant PROOF_OPENING_AT_ZETA_OMEGA_X = 0x300;
+  uint256 private constant PROOF_OPENING_AT_ZETA_OMEGA_Y = 0x320;
+
+  uint256 private constant PROOF_OPENING_QCP_AT_ZETA = 0x340;
+  uint256 private constant PROOF_COMMITMENTS_WIRES_CUSTOM_GATES = {{ hex (add 832 (mul (len .CommitmentConstraintIndexes) 32 ) )}};
+
+  // -> next part of proof is
+  // [ openings_selector_commits || commitments_wires_commit_api]
+
+  // -------- offset state
+
+  // challenges to check the claimed quotient
+  uint256 private constant STATE_ALPHA = 0x00;
+  uint256 private constant STATE_BETA = 0x20;
+  uint256 private constant STATE_GAMMA = 0x40;
+  uint256 private constant STATE_ZETA = 0x60;
+
+  // reusable value
+  uint256 private constant STATE_ALPHA_SQUARE_LAGRANGE_0 = 0x80;
+
+  // commitment to H
+  uint256 private constant STATE_FOLDED_H_X = 0xa0;
+  uint256 private constant STATE_FOLDED_H_Y = 0xc0;
+
+  // commitment to the linearised polynomial
+  uint256 private constant STATE_LINEARISED_POLYNOMIAL_X = 0xe0;
+  uint256 private constant STATE_LINEARISED_POLYNOMIAL_Y = 0x100;
+
+  // Folded proof for the opening of H, linearised poly, l, r, o, s_1, s_2, qcp
+  uint256 private constant STATE_FOLDED_CLAIMED_VALUES = 0x120;
+
+  // folded digests of H, linearised poly, l, r, o, s_1, s_2, qcp
+  uint256 private constant STATE_FOLDED_DIGESTS_X = 0x140;
+  uint256 private constant STATE_FOLDED_DIGESTS_Y = 0x160;
+
+  uint256 private constant STATE_PI = 0x180;
+
+  uint256 private constant STATE_ZETA_POWER_N_MINUS_ONE = 0x1a0;
+
+  uint256 private constant STATE_GAMMA_KZG = 0x1c0;
+
+  uint256 private constant STATE_SUCCESS = 0x1e0;
+  uint256 private constant STATE_CHECK_VAR = 0x200; // /!\ this slot is used for debugging only
+
+  uint256 private constant STATE_LAST_MEM = 0x220;
+
+  // -------- errors
+  uint256 private constant ERROR_STRING_ID = 0x08c379a000000000000000000000000000000000000000000000000000000000; // selector for function Error(string)
+
+  {{ if (gt (len .CommitmentConstraintIndexes) 0 )}}
+  // -------- utils (for hash_fr)
+	uint256 private constant HASH_FR_BB = 340282366920938463463374607431768211456; // 2**128
+	uint256 private constant HASH_FR_ZERO_UINT256 = 0;
+
+	uint8 private constant HASH_FR_LEN_IN_BYTES = 48;
+	uint8 private constant HASH_FR_SIZE_DOMAIN = 11;
+	uint8 private constant HASH_FR_ONE = 1;
+	uint8 private constant HASH_FR_TWO = 2;
+  {{ end }}
+  
+  /// Verify a Plonk proof.
+  /// Reverts if the proof or the public inputs are malformed.
+  /// @param proof serialised plonk proof (using gnark's MarshalSolidity)
+  /// @param public_inputs (must be reduced)
+  /// @return success true if the proof passes false otherwise
+  function Verify(bytes calldata proof, uint256[] calldata public_inputs) 
+  public view returns(bool success) {
+
+    assembly {
+
+      let mem := mload(0x40)
+      let freeMem := add(mem, STATE_LAST_MEM)
+
+      // sanity checks
+      check_number_of_public_inputs(public_inputs.length)
+      check_inputs_size(public_inputs.length, public_inputs.offset)
+      check_proof_size(proof.length)
+      check_proof_openings_size(proof.offset)
+
+      // compute the challenges
+      let prev_challenge_non_reduced
+      prev_challenge_non_reduced := derive_gamma(proof.offset, public_inputs.length, public_inputs.offset)
+      prev_challenge_non_reduced := derive_beta(prev_challenge_non_reduced)
+      prev_challenge_non_reduced := derive_alpha(proof.offset, prev_challenge_non_reduced)
+      derive_zeta(proof.offset, prev_challenge_non_reduced)
+
+      // evaluation of Z=Xⁿ-1 at ζ, we save this value
+      let zeta := mload(add(mem, STATE_ZETA))
+      let zeta_power_n_minus_one := addmod(pow(zeta, VK_DOMAIN_SIZE, freeMem), sub(R_MOD, 1), R_MOD)
+      mstore(add(mem, STATE_ZETA_POWER_N_MINUS_ONE), zeta_power_n_minus_one)
+
+      // public inputs contribution
+      let l_pi := sum_pi_wo_api_commit(public_inputs.offset, public_inputs.length, freeMem)
+      {{ if (gt (len .CommitmentConstraintIndexes) 0 ) -}}
+      let l_wocommit := sum_pi_commit(proof.offset, public_inputs.length, freeMem)
+      l_pi := addmod(l_wocommit, l_pi, R_MOD)
+      {{ end -}}
+      mstore(add(mem, STATE_PI), l_pi)
+
+      compute_alpha_square_lagrange_0()
+      verify_quotient_poly_eval_at_zeta(proof.offset)
+      fold_h(proof.offset)
+      compute_commitment_linearised_polynomial(proof.offset)
+      compute_gamma_kzg(proof.offset)
+      fold_state(proof.offset)
+      batch_verify_multi_points(proof.offset)
+
+      success := mload(add(mem, STATE_SUCCESS))
+
+      // Beginning errors -------------------------------------------------
+
+      function error_nb_public_inputs() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0x1d)
+        mstore(add(ptError, 0x44), "wrong number of public inputs")
+        revert(ptError, 0x64)
+      }
+
+      /// Called when an operation on Bn254 fails
+      /// @dev for instance when calling EcMul on a point not on Bn254.
+      function error_ec_op() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0x12)
+        mstore(add(ptError, 0x44), "error ec operation")
+        revert(ptError, 0x64)
+      }
+
+      /// Called when one of the public inputs is not reduced.
+      function error_inputs_size() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0x18)
+        mstore(add(ptError, 0x44), "inputs are bigger than r")
+        revert(ptError, 0x64)
+      }
+
+      /// Called when the size proof is not as expected
+      /// @dev to avoid overflow attack for instance
+      function error_proof_size() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0x10)
+        mstore(add(ptError, 0x44), "wrong proof size")
+        revert(ptError, 0x64)
+      }
+
+      /// Called when one the openings is bigger than r
+      /// The openings are the claimed evalutions of a polynomial
+      /// in a Kzg proof.
+      function error_proof_openings_size() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0x16)
+        mstore(add(ptError, 0x44), "openings bigger than r")
+        revert(ptError, 0x64)
+      }
+
+      function error_verify() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0xc)
+        mstore(add(ptError, 0x44), "error verify")
+        revert(ptError, 0x64)
+      }
+
+      function error_random_generation() {
+        let ptError := mload(0x40)
+        mstore(ptError, ERROR_STRING_ID) // selector for function Error(string)
+        mstore(add(ptError, 0x4), 0x20)
+        mstore(add(ptError, 0x24), 0x14)
+        mstore(add(ptError, 0x44), "error random gen kzg")
+        revert(ptError, 0x64)
+      }
+      // end errors -------------------------------------------------
+
+      // Beginning checks -------------------------------------------------
+      
+      /// @param s actual number of public inputs
+      function check_number_of_public_inputs(s) {
+        if iszero(eq(s, VK_NB_PUBLIC_INPUTS)) {
+          error_nb_public_inputs()
+        }
+      }
+    
+      /// Checks that the public inputs are < R_MOD.
+      /// @param s number of public inputs
+      /// @param p pointer to the public inputs array
+      function check_inputs_size(s, p) {
+        let input_checks := 1
+        for {let i} lt(i, s) {i:=add(i,1)}
+        {
+          input_checks := and(input_checks,lt(calldataload(p), R_MOD))
+          p := add(p, 0x20)
+        }
+        if iszero(input_checks) {
+          error_inputs_size()
+        }
+      }
+
+      /// Checks if the proof is of the correct size
+      /// @param actual_proof_size size of the proof (not the expected size)
+      function check_proof_size(actual_proof_size) {
+        let expected_proof_size := add(0x340, mul(VK_NB_CUSTOM_GATES,0x60))
+        if iszero(eq(actual_proof_size, expected_proof_size)) {
+         error_proof_size() 
+        }
+      }
+    
+      /// Checks if the multiple openings of the polynomials are < R_MOD.
+      /// @param aproof pointer to the beginning of the proof
+      /// @dev the 'a' prepending proof is to have a local name
+      function check_proof_openings_size(aproof) {
+  
+        let openings_check := 1
+      
+        // linearised polynomial at zeta
+        let p := add(aproof, PROOF_LINEARISED_POLYNOMIAL_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+
+        // quotient polynomial at zeta
+        p := add(aproof, PROOF_QUOTIENT_POLYNOMIAL_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+        
+        // PROOF_L_AT_ZETA
+        p := add(aproof, PROOF_L_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+
+        // PROOF_R_AT_ZETA
+        p := add(aproof, PROOF_R_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+
+        // PROOF_O_AT_ZETA
+        p := add(aproof, PROOF_O_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+
+        // PROOF_S1_AT_ZETA
+        p := add(aproof, PROOF_S1_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+        
+        // PROOF_S2_AT_ZETA
+        p := add(aproof, PROOF_S2_AT_ZETA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+
+        // PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA
+        p := add(aproof, PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA)
+        openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+
+        // PROOF_OPENING_QCP_AT_ZETA
+        
+        p := add(aproof, PROOF_OPENING_QCP_AT_ZETA)
+        for {let i:=0} lt(i, VK_NB_CUSTOM_GATES) {i:=add(i,1)}
+        {
+          openings_check := and(openings_check, lt(calldataload(p), R_MOD))
+          p := add(p, 0x20)
+        }
+      
+        if iszero(openings_check) {
+          error_proof_openings_size()
+        }
+
+      }
+      // end checks -------------------------------------------------
+
+      // Beginning challenges -------------------------------------------------
+
+      /// Derive gamma as Sha256(<transcript>)
+      /// @param aproof pointer to the proof
+      /// @param nb_pi number of public inputs
+      /// @param pi pointer to the array of public inputs
+      /// @return the challenge gamma, not reduced
+      /// @notice The transcript is the concatenation (in this order) of:
+      /// * the word "gamma" in ascii, equal to [0x67,0x61,0x6d, 0x6d, 0x61] and encoded as a uint256.
+      /// * the commitments to the permutation polynomials S1, S2, S3, where we concatenate the coordinates of those points
+      /// * the commitments of Ql, Qr, Qm, Qo, Qk
+      /// * the public inputs
+      /// * the commitments of the wires related to the custom gates (commitments_wires_commit_api)
+      /// * commitments to L, R, O (proof_<l,r,o>_com_<x,y>)
+      /// The data described above is written starting at mPtr. "gamma" lies on 5 bytes,
+      /// and is encoded as a uint256 number n. In basis b = 256, the number looks like this
+      /// [0 0 0 .. 0x67 0x61 0x6d, 0x6d, 0x61]. The first non zero entry is at position 27=0x1b
+      /// Gamma reduced (the actual challenge) is stored at add(state, state_gamma)
+      function derive_gamma(aproof, nb_pi, pi)->gamma_not_reduced {
+        
+        let state := mload(0x40)
+        let mPtr := add(state, STATE_LAST_MEM)
+
+        // gamma
+        // gamma in ascii is [0x67,0x61,0x6d, 0x6d, 0x61]
+        // (same for alpha, beta, zeta)
+        mstore(mPtr, 0x67616d6d61) // "gamma"
+
+        mstore(add(mPtr, 0x20), VK_S1_COM_X)
+        mstore(add(mPtr, 0x40), VK_S1_COM_Y)
+        mstore(add(mPtr, 0x60), VK_S2_COM_X)
+        mstore(add(mPtr, 0x80), VK_S2_COM_Y)
+        mstore(add(mPtr, 0xa0), VK_S3_COM_X)
+        mstore(add(mPtr, 0xc0), VK_S3_COM_Y)
+        mstore(add(mPtr, 0xe0), VK_QL_COM_X)
+        mstore(add(mPtr, 0x100), VK_QL_COM_Y)
+        mstore(add(mPtr, 0x120), VK_QR_COM_X)
+        mstore(add(mPtr, 0x140), VK_QR_COM_Y)
+        mstore(add(mPtr, 0x160), VK_QM_COM_X)
+        mstore(add(mPtr, 0x180), VK_QM_COM_Y)
+        mstore(add(mPtr, 0x1a0), VK_QO_COM_X)
+        mstore(add(mPtr, 0x1c0), VK_QO_COM_Y)
+        mstore(add(mPtr, 0x1e0), VK_QK_COM_X)
+        mstore(add(mPtr, 0x200), VK_QK_COM_Y)
+        {{ range $index, $element := .CommitmentConstraintIndexes}}
+        mstore(add(mPtr, {{ hex (add 544 (mul $index 64)) }}), VK_QCP_{{ $index }}_X)
+        mstore(add(mPtr, {{ hex (add 576 (mul $index 64)) }}), VK_QCP_{{ $index }}_Y)
+        {{ end }}
+        // public inputs
+        let _mPtr := add(mPtr, {{ hex (add (mul (len .CommitmentConstraintIndexes) 64) 544) }})
+        let size_pi_in_bytes := mul(nb_pi, 0x20)
+        calldatacopy(_mPtr, pi, size_pi_in_bytes)
+        _mPtr := add(_mPtr, size_pi_in_bytes)
+
+        // commitments to l, r, o
+        let size_commitments_lro_in_bytes := 0xc0
+        calldatacopy(_mPtr, aproof, size_commitments_lro_in_bytes)
+        _mPtr := add(_mPtr, size_commitments_lro_in_bytes)
+
+        // total size is :
+        // sizegamma(=0x5) + 11*64(=0x2c0)
+        // + nb_public_inputs*0x20
+        // + nb_custom gates*0x40
+        let size := add(0x2c5, size_pi_in_bytes)
+        {{ if (gt (len .CommitmentConstraintIndexes) 0 )}}
+        size := add(size, mul(VK_NB_CUSTOM_GATES, 0x40))
+        {{ end -}}
+        let l_success := staticcall(gas(), 0x2, add(mPtr, 0x1b), size, mPtr, 0x20) //0x1b -> 000.."gamma"
+        if iszero(l_success) {
+          error_verify()
+        }
+        gamma_not_reduced := mload(mPtr)
+        mstore(add(state, STATE_GAMMA), mod(gamma_not_reduced, R_MOD))
+      }
+
+      /// derive beta as Sha256<transcript>
+      /// @param gamma_not_reduced the previous challenge (gamma) not reduced
+      /// @return beta_not_reduced the next challenge, beta, not reduced
+      /// @notice the transcript consists of the previous challenge only.
+      /// The reduced version of beta is stored at add(state, state_beta)
+      function derive_beta(gamma_not_reduced)->beta_not_reduced{
+        
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+
+        // beta
+        mstore(mPtr, 0x62657461) // "beta"
+        mstore(add(mPtr, 0x20), gamma_not_reduced)
+        let l_success := staticcall(gas(), 0x2, add(mPtr, 0x1c), 0x24, mPtr, 0x20) //0x1b -> 000.."gamma"
+        if iszero(l_success) {
+          error_verify()
+        }
+        beta_not_reduced := mload(mPtr)
+        mstore(add(state, STATE_BETA), mod(beta_not_reduced, R_MOD))
+      }
+
+      /// derive alpha as sha256<transcript>
+      /// @param aproof pointer to the proof object
+      /// @param beta_not_reduced the previous challenge (beta) not reduced
+      /// @return alpha_not_reduced the next challenge, alpha, not reduced
+      /// @notice the transcript consists of the previous challenge (beta)
+      /// not reduced, the commitments to the wires associated to the QCP_i,
+      /// and the commitment to the grand product polynomial 
+      function derive_alpha(aproof, beta_not_reduced)->alpha_not_reduced {
+        
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+        let full_size := 0x65 // size("alpha") + 0x20 (previous challenge)
+
+        // alpha
+        mstore(mPtr, 0x616C706861) // "alpha"
+        let _mPtr := add(mPtr, 0x20)
+        mstore(_mPtr, beta_not_reduced)
+        _mPtr := add(_mPtr, 0x20)
+        {{ if (gt (len .CommitmentConstraintIndexes) 0 )}}
+        // Bsb22Commitments
+        let proof_bsb_commitments := add(aproof, PROOF_COMMITMENTS_WIRES_CUSTOM_GATES)
+        let size_bsb_commitments := mul(0x40, VK_NB_CUSTOM_GATES)
+        calldatacopy(_mPtr, proof_bsb_commitments, size_bsb_commitments)
+        _mPtr := add(_mPtr, size_bsb_commitments)
+        full_size := add(full_size, size_bsb_commitments)
+        {{ end }}
+        // [Z], the commitment to the grand product polynomial
+        calldatacopy(_mPtr, add(aproof, PROOF_GRAND_PRODUCT_COMMITMENT_X), 0x40)
+        let l_success := staticcall(gas(), 0x2, add(mPtr, 0x1b), full_size, mPtr, 0x20)
+        if iszero(l_success) {
+          error_verify()
+        }
+
+        alpha_not_reduced := mload(mPtr)
+        mstore(add(state, STATE_ALPHA), mod(alpha_not_reduced, R_MOD))
+      }
+
+      /// derive zeta as sha256<transcript>
+      /// @param aproof pointer to the proof object
+      /// @param alpha_not_reduced the previous challenge (alpha) not reduced
+      /// The transcript consists of the previous challenge and the commitment to
+      /// the quotient polynomial h.
+      function derive_zeta(aproof, alpha_not_reduced) {
+        
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+
+        // zeta
+        mstore(mPtr, 0x7a657461) // "zeta"
+        mstore(add(mPtr, 0x20), alpha_not_reduced)
+        calldatacopy(add(mPtr, 0x40), add(aproof, PROOF_H_0_X), 0xc0)
+        let l_success := staticcall(gas(), 0x2, add(mPtr, 0x1c), 0xe4, mPtr, 0x20)
+        if iszero(l_success) {
+          error_verify()
+        }
+        let zeta_not_reduced := mload(mPtr)
+        mstore(add(state, STATE_ZETA), mod(zeta_not_reduced, R_MOD))
+      }
+      // END challenges -------------------------------------------------
+
+      // BEGINNING compute_pi -------------------------------------------------
+
+      /// sum_pi_wo_api_commit computes the public inputs contributions,
+      /// except for the public inputs coming from the custom gate
+      /// @param ins pointer to the public inputs
+      /// @param n number of public inputs
+      /// @param mPtr free memory
+      /// @return pi_wo_commit public inputs contribution (except the public inputs coming from the custom gate)
+      function sum_pi_wo_api_commit(ins, n, mPtr)->pi_wo_commit {
+        
+        let state := mload(0x40)
+        let z := mload(add(state, STATE_ZETA))
+        let zpnmo := mload(add(state, STATE_ZETA_POWER_N_MINUS_ONE))
+
+        let li := mPtr
+        batch_compute_lagranges_at_z(z, zpnmo, n, li)
+
+        let tmp := 0
+        for {let i:=0} lt(i,n) {i:=add(i,1)}
+        {
+          tmp := mulmod(mload(li), calldataload(ins), R_MOD)
+          pi_wo_commit := addmod(pi_wo_commit, tmp, R_MOD)
+          li := add(li, 0x20)
+          ins := add(ins, 0x20)
+        }
+        
+      }
+
+      /// batch_compute_lagranges_at_z computes [L_0(z), .., L_{n-1}(z)]
+      /// @param z point at which the Lagranges are evaluated
+      /// @param zpnmo ζⁿ-1
+      /// @param n number of public inputs (number of Lagranges to compute)
+      /// @param mPtr pointer to which the results are stored
+      function batch_compute_lagranges_at_z(z, zpnmo, n, mPtr) {
+
+        let zn := mulmod(zpnmo, VK_INV_DOMAIN_SIZE, R_MOD) // 1/n * (ζⁿ - 1)
+        
+        let _w := 1
+        let _mPtr := mPtr
+        for {let i:=0} lt(i,n) {i:=add(i,1)}
+        {
+          mstore(_mPtr, addmod(z,sub(R_MOD, _w), R_MOD))
+          _w := mulmod(_w, VK_OMEGA, R_MOD)
+          _mPtr := add(_mPtr, 0x20)
+        }
+        batch_invert(mPtr, n, _mPtr)
+        _mPtr := mPtr
+        _w := 1
+        for {let i:=0} lt(i,n) {i:=add(i,1)}
+        {
+          mstore(_mPtr, mulmod(mulmod(mload(_mPtr), zn , R_MOD), _w, R_MOD))
+          _mPtr := add(_mPtr, 0x20)
+          _w := mulmod(_w, VK_OMEGA, R_MOD)
+        }
+      } 
+
+      /// @notice Montgomery trick for batch inversion mod R_MOD
+      /// @param ins pointer to the data to batch invert
+      /// @param number of elements to batch invert
+      /// @param mPtr free memory
+      function batch_invert(ins, nb_ins, mPtr) {
+        mstore(mPtr, 1)
+        let offset := 0
+        for {let i:=0} lt(i, nb_ins) {i:=add(i,1)}
+        {
+          let prev := mload(add(mPtr, offset))
+          let cur := mload(add(ins, offset))
+          cur := mulmod(prev, cur, R_MOD)
+          offset := add(offset, 0x20)
+          mstore(add(mPtr, offset), cur)
+        }
+        ins := add(ins, sub(offset, 0x20))
+        mPtr := add(mPtr, offset)
+        let inv := pow(mload(mPtr), sub(R_MOD,2), add(mPtr, 0x20))
+        for {let i:=0} lt(i, nb_ins) {i:=add(i,1)}
+        {
+          mPtr := sub(mPtr, 0x20)
+          let tmp := mload(ins)
+          let cur := mulmod(inv, mload(mPtr), R_MOD)
+          mstore(ins, cur)
+          inv := mulmod(inv, tmp, R_MOD)
+          ins := sub(ins, 0x20)
+        }
+      }
+
+      {{ if (gt (len .CommitmentConstraintIndexes) 0 )}}
+      /// Public inputs (the ones coming from the custom gate) contribution
+      /// @param aproof pointer to the proof
+      /// @param nb_public_inputs number of public inputs
+      /// @param mPtr pointer to free memory
+      /// @return pi_commit custom gate public inputs contribution
+      function sum_pi_commit(aproof, nb_public_inputs, mPtr)->pi_commit {
+
+        let state := mload(0x40)
+        let z := mload(add(state, STATE_ZETA))
+        let zpnmo := mload(add(state, STATE_ZETA_POWER_N_MINUS_ONE))
+
+        let p := add(aproof, PROOF_COMMITMENTS_WIRES_CUSTOM_GATES)
+
+        let h_fr, ith_lagrange
+       
+        {{ range $index, $element := .CommitmentConstraintIndexes}}
+        h_fr := hash_fr(calldataload(p), calldataload(add(p, 0x20)), mPtr)
+        ith_lagrange := compute_ith_lagrange_at_z(z, zpnmo, add(nb_public_inputs, VK_INDEX_COMMIT_API{{ $index }}), mPtr)
+        pi_commit := addmod(pi_commit, mulmod(h_fr, ith_lagrange, R_MOD), R_MOD)
+        p := add(p, 0x40)
+        {{ end }}
+
+      }
+
+      /// Computes L_i(zeta) =  ωⁱ/n * (ζⁿ-1)/(ζ-ωⁱ) where:
+      /// @param z zeta
+      /// @param zpmno ζⁿ-1
+      /// @param i i-th lagrange
+      /// @param mPtr free memory
+      /// @return res = ωⁱ/n * (ζⁿ-1)/(ζ-ωⁱ) 
+      function compute_ith_lagrange_at_z(z, zpnmo, i, mPtr)->res {
+
+        let w := pow(VK_OMEGA, i, mPtr) // w**i
+        i := addmod(z, sub(R_MOD, w), R_MOD) // z-w**i
+        w := mulmod(w, VK_INV_DOMAIN_SIZE, R_MOD) // w**i/n
+        i := pow(i, sub(R_MOD,2), mPtr) // (z-w**i)**-1
+        w := mulmod(w, i, R_MOD) // w**i/n*(z-w)**-1
+        res := mulmod(w, zpnmo, R_MOD)
+      
+      }
+
+      /// @dev https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-06#section-5.2
+      /// @param x x coordinate of a point on Bn254(𝔽_p)
+      /// @param y y coordinate of a point on Bn254(𝔽_p)
+      /// @param mPtr free memory
+      /// @return res an element mod R_MOD
+      function hash_fr(x, y, mPtr)->res {
+
+        // [0x00, .. , 0x00 || x, y, || 0, 48, 0, dst, HASH_FR_SIZE_DOMAIN]
+        // <-  64 bytes  ->  <-64b -> <-       1 bytes each     ->
+
+        // [0x00, .., 0x00] 64 bytes of zero
+        mstore(mPtr, HASH_FR_ZERO_UINT256)
+        mstore(add(mPtr, 0x20), HASH_FR_ZERO_UINT256)
+    
+        // msg =  x || y , both on 32 bytes
+        mstore(add(mPtr, 0x40), x)
+        mstore(add(mPtr, 0x60), y)
+
+        // 0 || 48 || 0 all on 1 byte
+        mstore8(add(mPtr, 0x80), 0)
+        mstore8(add(mPtr, 0x81), HASH_FR_LEN_IN_BYTES)
+        mstore8(add(mPtr, 0x82), 0)
+
+        // "BSB22-Plonk" = [42, 53, 42, 32, 32, 2d, 50, 6c, 6f, 6e, 6b,]
+        mstore8(add(mPtr, 0x83), 0x42)
+        mstore8(add(mPtr, 0x84), 0x53)
+        mstore8(add(mPtr, 0x85), 0x42)
+        mstore8(add(mPtr, 0x86), 0x32)
+        mstore8(add(mPtr, 0x87), 0x32)
+        mstore8(add(mPtr, 0x88), 0x2d)
+        mstore8(add(mPtr, 0x89), 0x50)
+        mstore8(add(mPtr, 0x8a), 0x6c)
+        mstore8(add(mPtr, 0x8b), 0x6f)
+        mstore8(add(mPtr, 0x8c), 0x6e)
+        mstore8(add(mPtr, 0x8d), 0x6b)
+
+        // size domain
+        mstore8(add(mPtr, 0x8e), HASH_FR_SIZE_DOMAIN)
+
+        let l_success := staticcall(gas(), 0x2, mPtr, 0x8f, mPtr, 0x20)
+        if iszero(l_success) {
+          error_verify()
+        }
+
+        let b0 := mload(mPtr)
+
+        // [b0         || one || dst || HASH_FR_SIZE_DOMAIN]
+        // <-64bytes ->  <-    1 byte each      ->
+        mstore8(add(mPtr, 0x20), HASH_FR_ONE) // 1
+        
+        mstore8(add(mPtr, 0x21), 0x42) // dst
+        mstore8(add(mPtr, 0x22), 0x53)
+        mstore8(add(mPtr, 0x23), 0x42)
+        mstore8(add(mPtr, 0x24), 0x32)
+        mstore8(add(mPtr, 0x25), 0x32)
+        mstore8(add(mPtr, 0x26), 0x2d)
+        mstore8(add(mPtr, 0x27), 0x50)
+        mstore8(add(mPtr, 0x28), 0x6c)
+        mstore8(add(mPtr, 0x29), 0x6f)
+        mstore8(add(mPtr, 0x2a), 0x6e)
+        mstore8(add(mPtr, 0x2b), 0x6b)
+
+        mstore8(add(mPtr, 0x2c), HASH_FR_SIZE_DOMAIN) // size domain
+        l_success := staticcall(gas(), 0x2, mPtr, 0x2d, mPtr, 0x20)
+        if iszero(l_success) {
+          error_verify()
+        }
+
+        // b1 is located at mPtr. We store b2 at add(mPtr, 0x20)
+
+        // [b0^b1      || two || dst || HASH_FR_SIZE_DOMAIN]
+        // <-64bytes ->  <-    1 byte each      ->
+        mstore(add(mPtr, 0x20), xor(mload(mPtr), b0))
+        mstore8(add(mPtr, 0x40), HASH_FR_TWO)
+
+        mstore8(add(mPtr, 0x41), 0x42) // dst
+        mstore8(add(mPtr, 0x42), 0x53)
+        mstore8(add(mPtr, 0x43), 0x42)
+        mstore8(add(mPtr, 0x44), 0x32)
+        mstore8(add(mPtr, 0x45), 0x32)
+        mstore8(add(mPtr, 0x46), 0x2d)
+        mstore8(add(mPtr, 0x47), 0x50)
+        mstore8(add(mPtr, 0x48), 0x6c)
+        mstore8(add(mPtr, 0x49), 0x6f)
+        mstore8(add(mPtr, 0x4a), 0x6e)
+        mstore8(add(mPtr, 0x4b), 0x6b)
+
+        mstore8(add(mPtr, 0x4c), HASH_FR_SIZE_DOMAIN) // size domain
+
+        let offset := add(mPtr, 0x20)
+        l_success := staticcall(gas(), 0x2, offset, 0x2d, offset, 0x20)
+        if iszero(l_success) {
+          error_verify()
+        }
+
+        // at this point we have mPtr = [ b1 || b2] where b1 is on 32byes and b2 in 16bytes.
+        // we interpret it as a big integer mod r in big endian (similar to regular decimal notation)
+        // the result is then 2**(8*16)*mPtr[32:] + mPtr[32:48]
+        res := mulmod(mload(mPtr), HASH_FR_BB, R_MOD) // <- res = 2**128 * mPtr[:32]
+        let b1 := shr(128, mload(add(mPtr, 0x20))) // b1 <- [0, 0, .., 0 ||  b2[:16] ]
+        res := addmod(res, b1, R_MOD)
+
+      }
+      {{ end }}
+      // END compute_pi -------------------------------------------------
+
+      /// @notice compute α² * 1/n * (ζ{n}-1)/(ζ - 1) where
+      /// *  α = challenge derived in derive_gamma_beta_alpha_zeta
+      /// * n = vk_domain_size
+      /// * ω = vk_omega (generator of the multiplicative cyclic group of order n in (ℤ/rℤ)*)
+      /// * ζ = zeta (challenge derived with Fiat Shamir)
+      function compute_alpha_square_lagrange_0() {   
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+
+        let res := mload(add(state, STATE_ZETA_POWER_N_MINUS_ONE))
+        let den := addmod(mload(add(state, STATE_ZETA)), sub(R_MOD, 1), R_MOD)
+        den := pow(den, sub(R_MOD, 2), mPtr)
+        den := mulmod(den, VK_INV_DOMAIN_SIZE, R_MOD)
+        res := mulmod(den, res, R_MOD)
+
+        let l_alpha := mload(add(state, STATE_ALPHA))
+        res := mulmod(res, l_alpha, R_MOD)
+        res := mulmod(res, l_alpha, R_MOD)
+        mstore(add(state, STATE_ALPHA_SQUARE_LAGRANGE_0), res)
+      }
+
+      /// @notice follows alg. p.13 of https://eprint.iacr.org/2019/953.pdf
+      /// with t₁ = t₂ = 1, and the proofs are ([digest] + [quotient] +purported evaluation):
+      /// * [state_folded_state_digests], [proof_batch_opening_at_zeta_x], state_folded_evals
+      /// * [proof_grand_product_commitment], [proof_opening_at_zeta_omega_x], [proof_grand_product_at_zeta_omega]
+      /// @param aproof pointer to the proof
+      function batch_verify_multi_points(aproof) {
+        let state := mload(0x40)
+        let mPtr := add(state, STATE_LAST_MEM)
+
+        // derive a random number. As there is no random generator, we
+        // do an FS like challenge derivation, depending on both digests and
+        // ζ to ensure that the prover cannot control the random numger.
+        // Note: adding the other point ζω is not needed, as ω is known beforehand.
+        mstore(mPtr, mload(add(state, STATE_FOLDED_DIGESTS_X)))
+        mstore(add(mPtr, 0x20), mload(add(state, STATE_FOLDED_DIGESTS_Y)))
+        mstore(add(mPtr, 0x40), calldataload(add(aproof, PROOF_BATCH_OPENING_AT_ZETA_X)))
+        mstore(add(mPtr, 0x60), calldataload(add(aproof, PROOF_BATCH_OPENING_AT_ZETA_Y)))
+        mstore(add(mPtr, 0x80), calldataload(add(aproof, PROOF_GRAND_PRODUCT_COMMITMENT_X)))
+        mstore(add(mPtr, 0xa0), calldataload(add(aproof, PROOF_GRAND_PRODUCT_COMMITMENT_Y)))
+        mstore(add(mPtr, 0xc0), calldataload(add(aproof, PROOF_OPENING_AT_ZETA_OMEGA_X)))
+        mstore(add(mPtr, 0xe0), calldataload(add(aproof, PROOF_OPENING_AT_ZETA_OMEGA_Y)))
+        mstore(add(mPtr, 0x100), mload(add(state, STATE_ZETA)))
+        mstore(add(mPtr, 0x120), mload(add(state, STATE_GAMMA_KZG)))
+        let random := staticcall(gas(), 0x2, mPtr, 0x140, mPtr, 0x20)
+        if iszero(random){
+          error_random_generation()
+        }
+        random := mod(mload(mPtr), R_MOD) // use the same variable as we are one variable away from getting stack-too-deep error...
+
+        let folded_quotients := mPtr
+        mPtr := add(folded_quotients, 0x40)
+        mstore(folded_quotients, calldataload(add(aproof, PROOF_BATCH_OPENING_AT_ZETA_X)))
+        mstore(add(folded_quotients, 0x20), calldataload(add(aproof, PROOF_BATCH_OPENING_AT_ZETA_Y)))
+        point_acc_mul_calldata(folded_quotients, add(aproof, PROOF_OPENING_AT_ZETA_OMEGA_X), random, mPtr)
+
+        let folded_digests := add(state, STATE_FOLDED_DIGESTS_X)
+        point_acc_mul_calldata(folded_digests, add(aproof, PROOF_GRAND_PRODUCT_COMMITMENT_X), random, mPtr)
+
+        let folded_evals := add(state, STATE_FOLDED_CLAIMED_VALUES)
+        fr_acc_mul_calldata(folded_evals, add(aproof, PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA), random)
+
+        let folded_evals_commit := mPtr
+        mPtr := add(folded_evals_commit, 0x40)
+        mstore(folded_evals_commit, G1_SRS_X)
+        mstore(add(folded_evals_commit, 0x20), G1_SRS_Y)
+        mstore(add(folded_evals_commit, 0x40), mload(folded_evals))
+        let check_staticcall := staticcall(gas(), 7, folded_evals_commit, 0x60, folded_evals_commit, 0x40)
+        if iszero(check_staticcall) {
+          error_verify()
+        }
+
+        let folded_evals_commit_y := add(folded_evals_commit, 0x20)
+        mstore(folded_evals_commit_y, sub(P_MOD, mload(folded_evals_commit_y)))
+        point_add(folded_digests, folded_digests, folded_evals_commit, mPtr)
+
+        let folded_points_quotients := mPtr
+        mPtr := add(mPtr, 0x40)
+        point_mul_calldata(
+          folded_points_quotients,
+          add(aproof, PROOF_BATCH_OPENING_AT_ZETA_X),
+          mload(add(state, STATE_ZETA)),
+          mPtr
+        )
+        let zeta_omega := mulmod(mload(add(state, STATE_ZETA)), VK_OMEGA, R_MOD)
+        random := mulmod(random, zeta_omega, R_MOD)
+        point_acc_mul_calldata(folded_points_quotients, add(aproof, PROOF_OPENING_AT_ZETA_OMEGA_X), random, mPtr)
+
+        point_add(folded_digests, folded_digests, folded_points_quotients, mPtr)
+
+        let folded_quotients_y := add(folded_quotients, 0x20)
+        mstore(folded_quotients_y, sub(P_MOD, mload(folded_quotients_y)))
+
+        mstore(mPtr, mload(folded_digests))
+        mstore(add(mPtr, 0x20), mload(add(folded_digests, 0x20)))
+        mstore(add(mPtr, 0x40), G2_SRS_0_X_0) // the 4 lines are the canonical G2 point on BN254
+        mstore(add(mPtr, 0x60), G2_SRS_0_X_1)
+        mstore(add(mPtr, 0x80), G2_SRS_0_Y_0)
+        mstore(add(mPtr, 0xa0), G2_SRS_0_Y_1)
+        mstore(add(mPtr, 0xc0), mload(folded_quotients))
+        mstore(add(mPtr, 0xe0), mload(add(folded_quotients, 0x20)))
+        mstore(add(mPtr, 0x100), G2_SRS_1_X_0)
+        mstore(add(mPtr, 0x120), G2_SRS_1_X_1)
+        mstore(add(mPtr, 0x140), G2_SRS_1_Y_0)
+        mstore(add(mPtr, 0x160), G2_SRS_1_Y_1)
+        check_pairing_kzg(mPtr)
+      }
+
+      /// @notice check_pairing_kzg checks the result of the final pairing product of the batched
+      /// kzg verification. The purpose of this function is to avoid exhausting the stack
+      /// in the function batch_verify_multi_points.
+      /// @param mPtr pointer storing the tuple of pairs
+      function check_pairing_kzg(mPtr) {
+        let state := mload(0x40)
+
+        // TODO test the staticcall using the method from audit_4-5
+        let l_success := staticcall(gas(), 8, mPtr, 0x180, 0x00, 0x20)
+        let res_pairing := mload(0x00)
+        let s_success := mload(add(state, STATE_SUCCESS))
+        res_pairing := and(and(res_pairing, l_success), s_success)
+        mstore(add(state, STATE_SUCCESS), res_pairing)
+      }
+
+      /// @notice Fold the opening proofs at ζ:
+      /// * at state+state_folded_digest we store: [H] + γ[Linearised_polynomial]+γ²[L] + γ³[R] + γ⁴[O] + γ⁵[S₁] +γ⁶[S₂] + ∑ᵢγ⁶⁺ⁱ[Pi_{i}]
+      /// * at state+state_folded_claimed_values we store: H(ζ) + γLinearised_polynomial(ζ)+γ²L(ζ) + γ³R(ζ)+ γ⁴O(ζ) + γ⁵S₁(ζ) +γ⁶S₂(ζ) + ∑ᵢγ⁶⁺ⁱPi_{i}(ζ)
+      /// @param aproof pointer to the proof
+      /// acc_gamma stores the γⁱ
+      function fold_state(aproof) {
+
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+        let mPtr20 := add(mPtr, 0x20)
+        let mPtr40 := add(mPtr, 0x40)
+
+        let l_gamma_kzg := mload(add(state, STATE_GAMMA_KZG))
+        let acc_gamma := l_gamma_kzg
+        let state_folded_digests := add(state, STATE_FOLDED_DIGESTS_X)
+
+        mstore(add(state, STATE_FOLDED_DIGESTS_X), mload(add(state, STATE_FOLDED_H_X)))
+        mstore(add(state, STATE_FOLDED_DIGESTS_Y), mload(add(state, STATE_FOLDED_H_Y)))
+        mstore(add(state, STATE_FOLDED_CLAIMED_VALUES), calldataload(add(aproof, PROOF_QUOTIENT_POLYNOMIAL_AT_ZETA)))
+
+        point_acc_mul(state_folded_digests, add(state, STATE_LINEARISED_POLYNOMIAL_X), acc_gamma, mPtr)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), add(aproof, PROOF_LINEARISED_POLYNOMIAL_AT_ZETA), acc_gamma)
+
+        acc_gamma := mulmod(acc_gamma, l_gamma_kzg, R_MOD)
+        point_acc_mul_calldata(add(state, STATE_FOLDED_DIGESTS_X), add(aproof, PROOF_L_COM_X), acc_gamma, mPtr)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), add(aproof, PROOF_L_AT_ZETA), acc_gamma)
+
+        acc_gamma := mulmod(acc_gamma, l_gamma_kzg, R_MOD)
+        point_acc_mul_calldata(state_folded_digests, add(aproof, PROOF_R_COM_X), acc_gamma, mPtr)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), add(aproof, PROOF_R_AT_ZETA), acc_gamma)
+
+        acc_gamma := mulmod(acc_gamma, l_gamma_kzg, R_MOD)
+        point_acc_mul_calldata(state_folded_digests, add(aproof, PROOF_O_COM_X), acc_gamma, mPtr)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), add(aproof, PROOF_O_AT_ZETA), acc_gamma)
+
+        acc_gamma := mulmod(acc_gamma, l_gamma_kzg, R_MOD)
+        mstore(mPtr, VK_S1_COM_X)
+        mstore(mPtr20, VK_S1_COM_Y)
+        point_acc_mul(state_folded_digests, mPtr, acc_gamma, mPtr40)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), add(aproof, PROOF_S1_AT_ZETA), acc_gamma)
+
+        acc_gamma := mulmod(acc_gamma, l_gamma_kzg, R_MOD)
+        mstore(mPtr, VK_S2_COM_X)
+        mstore(mPtr20, VK_S2_COM_Y)
+        point_acc_mul(state_folded_digests, mPtr, acc_gamma, mPtr40)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), add(aproof, PROOF_S2_AT_ZETA), acc_gamma)
+
+        {{- if (gt (len .CommitmentConstraintIndexes) 0 ) }}
+        let poscaz := add(aproof, PROOF_OPENING_QCP_AT_ZETA)
+        {{ end -}}
+
+        {{ range $index, $element := .CommitmentConstraintIndexes }}
+        acc_gamma := mulmod(acc_gamma, l_gamma_kzg, R_MOD)
+        mstore(mPtr, VK_QCP_{{ $index }}_X)
+        mstore(mPtr20, VK_QCP_{{ $index }}_Y)
+        point_acc_mul(state_folded_digests, mPtr, acc_gamma, mPtr40)
+        fr_acc_mul_calldata(add(state, STATE_FOLDED_CLAIMED_VALUES), poscaz, acc_gamma)
+        poscaz := add(poscaz, 0x20)
+        {{ end }}
+
+      }
+
+      /// @notice generate the challenge (using Fiat Shamir) to fold the opening proofs
+      /// at ζ.
+      /// The process for deriving γ is the same as in derive_gamma but this time the inputs are
+      /// in this order (the [] means it's a commitment):
+      /// * ζ
+      /// * [H] ( = H₁ + ζᵐ⁺²*H₂ + ζ²⁽ᵐ⁺²⁾*H₃ )
+      /// * [Linearised polynomial]
+      /// * [L], [R], [O]
+      /// * [S₁] [S₂]
+      /// * [Pi_{i}] (wires associated to custom gates)
+      /// Then there are the purported evaluations of the previous committed polynomials:
+      /// * H(ζ)
+      /// * Linearised_polynomial(ζ)
+      /// * L(ζ), R(ζ), O(ζ), S₁(ζ), S₂(ζ)
+      /// * Pi_{i}(ζ)
+      /// * Z(ζω)
+      /// @param aproof pointer to the proof
+      function compute_gamma_kzg(aproof) {
+
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+        mstore(mPtr, 0x67616d6d61) // "gamma"
+        mstore(add(mPtr, 0x20), mload(add(state, STATE_ZETA)))
+        mstore(add(mPtr,0x40), mload(add(state, STATE_FOLDED_H_X)))
+        mstore(add(mPtr,0x60), mload(add(state, STATE_FOLDED_H_Y)))
+        mstore(add(mPtr,0x80), mload(add(state, STATE_LINEARISED_POLYNOMIAL_X)))
+        mstore(add(mPtr,0xa0), mload(add(state, STATE_LINEARISED_POLYNOMIAL_Y)))
+        calldatacopy(add(mPtr, 0xc0), add(aproof, PROOF_L_COM_X), 0xc0)
+        mstore(add(mPtr,0x180), VK_S1_COM_X)
+        mstore(add(mPtr,0x1a0), VK_S1_COM_Y)
+        mstore(add(mPtr,0x1c0), VK_S2_COM_X)
+        mstore(add(mPtr,0x1e0), VK_S2_COM_Y)
+        
+        let offset := 0x200
+        {{ range $index, $element := .CommitmentConstraintIndexes }}
+        mstore(add(mPtr,offset), VK_QCP_{{ $index }}_X)
+        mstore(add(mPtr,add(offset, 0x20)), VK_QCP_{{ $index }}_Y)
+        offset := add(offset, 0x40)
+        {{ end }}
+
+        mstore(add(mPtr, offset), calldataload(add(aproof, PROOF_QUOTIENT_POLYNOMIAL_AT_ZETA)))
+        mstore(add(mPtr, add(offset, 0x20)), calldataload(add(aproof, PROOF_LINEARISED_POLYNOMIAL_AT_ZETA)))
+        mstore(add(mPtr, add(offset, 0x40)), calldataload(add(aproof, PROOF_L_AT_ZETA)))
+        mstore(add(mPtr, add(offset, 0x60)), calldataload(add(aproof, PROOF_R_AT_ZETA)))
+        mstore(add(mPtr, add(offset, 0x80)), calldataload(add(aproof, PROOF_O_AT_ZETA)))
+        mstore(add(mPtr, add(offset, 0xa0)), calldataload(add(aproof, PROOF_S1_AT_ZETA)))
+        mstore(add(mPtr, add(offset, 0xc0)), calldataload(add(aproof, PROOF_S2_AT_ZETA)))
+
+        let _mPtr := add(mPtr, add(offset, 0xe0))
+        {{ if (gt (len .CommitmentConstraintIndexes) 0 )}}
+        let _poscaz := add(aproof, PROOF_OPENING_QCP_AT_ZETA)
+        for {let i:=0} lt(i, VK_NB_CUSTOM_GATES) {i:=add(i,1)}
+        {
+          mstore(_mPtr, calldataload(_poscaz))
+          _poscaz := add(_poscaz, 0x20)
+          _mPtr := add(_mPtr, 0x20)
+        }
+        {{ end }}
+
+        mstore(_mPtr, calldataload(add(aproof, PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA)))
+
+        let start_input := 0x1b // 00.."gamma"
+        let size_input := add(0x17, mul(VK_NB_CUSTOM_GATES,3)) // number of 32bytes elmts = 0x17 (zeta+2*7+7 for the digests+openings) + 2*VK_NB_CUSTOM_GATES (for the commitments of the selectors) + VK_NB_CUSTOM_GATES (for the openings of the selectors)
+        size_input := add(0x5, mul(size_input, 0x20)) // size in bytes: 15*32 bytes + 5 bytes for gamma
+        let check_staticcall := staticcall(gas(), 0x2, add(mPtr,start_input), size_input, add(state, STATE_GAMMA_KZG), 0x20)
+        if iszero(check_staticcall) {
+          error_verify()
+        }
+        mstore(add(state, STATE_GAMMA_KZG), mod(mload(add(state, STATE_GAMMA_KZG)), R_MOD))
+      }
+
+      function compute_commitment_linearised_polynomial_ec(aproof, s1, s2) {
+        let state := mload(0x40)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+
+        mstore(mPtr, VK_QL_COM_X)
+        mstore(add(mPtr, 0x20), VK_QL_COM_Y)
+        point_mul(
+          add(state, STATE_LINEARISED_POLYNOMIAL_X),
+          mPtr,
+          calldataload(add(aproof, PROOF_L_AT_ZETA)),
+          add(mPtr, 0x40)
+        )
+
+        mstore(mPtr, VK_QR_COM_X)
+        mstore(add(mPtr, 0x20), VK_QR_COM_Y)
+        point_acc_mul(
+          add(state, STATE_LINEARISED_POLYNOMIAL_X),
+          mPtr,
+          calldataload(add(aproof, PROOF_R_AT_ZETA)),
+          add(mPtr, 0x40)
+        )
+
+        let rl := mulmod(calldataload(add(aproof, PROOF_L_AT_ZETA)), calldataload(add(aproof, PROOF_R_AT_ZETA)), R_MOD)
+        mstore(mPtr, VK_QM_COM_X)
+        mstore(add(mPtr, 0x20), VK_QM_COM_Y)
+        point_acc_mul(add(state, STATE_LINEARISED_POLYNOMIAL_X), mPtr, rl, add(mPtr, 0x40))
+
+        mstore(mPtr, VK_QO_COM_X)
+        mstore(add(mPtr, 0x20), VK_QO_COM_Y)
+        point_acc_mul(
+          add(state, STATE_LINEARISED_POLYNOMIAL_X),
+          mPtr,
+          calldataload(add(aproof, PROOF_O_AT_ZETA)),
+          add(mPtr, 0x40)
+        )
+
+        mstore(mPtr, VK_QK_COM_X)
+        mstore(add(mPtr, 0x20), VK_QK_COM_Y)
+        point_add(
+          add(state, STATE_LINEARISED_POLYNOMIAL_X),
+          add(state, STATE_LINEARISED_POLYNOMIAL_X),
+          mPtr,
+          add(mPtr, 0x40)
+        )
+
+        let commits_api_at_zeta := add(aproof, PROOF_OPENING_QCP_AT_ZETA)
+        let commits_api := add(aproof, PROOF_COMMITMENTS_WIRES_CUSTOM_GATES)
+        for {
+          let i := 0
+        } lt(i, VK_NB_CUSTOM_GATES) {
+          i := add(i, 1)
+        } {
+          mstore(mPtr, calldataload(commits_api))
+          mstore(add(mPtr, 0x20), calldataload(add(commits_api, 0x20)))
+          point_acc_mul(
+            add(state, STATE_LINEARISED_POLYNOMIAL_X),
+            mPtr,
+            calldataload(commits_api_at_zeta),
+            add(mPtr, 0x40)
+          )
+          commits_api_at_zeta := add(commits_api_at_zeta, 0x20)
+          commits_api := add(commits_api, 0x40)
+        }
+
+        mstore(mPtr, VK_S3_COM_X)
+        mstore(add(mPtr, 0x20), VK_S3_COM_Y)
+        point_acc_mul(add(state, STATE_LINEARISED_POLYNOMIAL_X), mPtr, s1, add(mPtr, 0x40))
+
+        mstore(mPtr, calldataload(add(aproof, PROOF_GRAND_PRODUCT_COMMITMENT_X)))
+        mstore(add(mPtr, 0x20), calldataload(add(aproof, PROOF_GRAND_PRODUCT_COMMITMENT_Y)))
+        point_acc_mul(add(state, STATE_LINEARISED_POLYNOMIAL_X), mPtr, s2, add(mPtr, 0x40))
+      }
+
+      /// @notice Compute the commitment to the linearized polynomial equal to
+      ///	L(ζ)[Qₗ]+r(ζ)[Qᵣ]+R(ζ)L(ζ)[Qₘ]+O(ζ)[Qₒ]+[Qₖ]+Σᵢqc'ᵢ(ζ)[BsbCommitmentᵢ] +
+      ///	α*( Z(μζ)(L(ζ)+β*S₁(ζ)+γ)*(R(ζ)+β*S₂(ζ)+γ)[S₃]-[Z](L(ζ)+β*id_{1}(ζ)+γ)*(R(ζ)+β*id_{2(ζ)+γ)*(O(ζ)+β*id_{3}(ζ)+γ) ) +
+      ///	α²*L₁(ζ)[Z]
+      /// where
+      /// * id_1 = id, id_2 = vk_coset_shift*id, id_3 = vk_coset_shift^{2}*id
+      /// * the [] means that it's a commitment (i.e. a point on Bn254(F_p))
+      /// @param aproof pointer to the proof
+      function compute_commitment_linearised_polynomial(aproof) {
+        let state := mload(0x40)
+        let l_beta := mload(add(state, STATE_BETA))
+        let l_gamma := mload(add(state, STATE_GAMMA))
+        let l_zeta := mload(add(state, STATE_ZETA))
+        let l_alpha := mload(add(state, STATE_ALPHA))
+
+        let u := mulmod(calldataload(add(aproof, PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA)), l_beta, R_MOD)
+        let v := mulmod(l_beta, calldataload(add(aproof, PROOF_S1_AT_ZETA)), R_MOD)
+        v := addmod(v, calldataload(add(aproof, PROOF_L_AT_ZETA)), R_MOD)
+        v := addmod(v, l_gamma, R_MOD)
+
+        let w := mulmod(l_beta, calldataload(add(aproof, PROOF_S2_AT_ZETA)), R_MOD)
+        w := addmod(w, calldataload(add(aproof, PROOF_R_AT_ZETA)), R_MOD)
+        w := addmod(w, l_gamma, R_MOD)
+
+        let s1 := mulmod(u, v, R_MOD)
+        s1 := mulmod(s1, w, R_MOD)
+        s1 := mulmod(s1, l_alpha, R_MOD)
+
+        let coset_square := mulmod(VK_COSET_SHIFT, VK_COSET_SHIFT, R_MOD)
+        let betazeta := mulmod(l_beta, l_zeta, R_MOD)
+        u := addmod(betazeta, calldataload(add(aproof, PROOF_L_AT_ZETA)), R_MOD)
+        u := addmod(u, l_gamma, R_MOD)
+
+        v := mulmod(betazeta, VK_COSET_SHIFT, R_MOD)
+        v := addmod(v, calldataload(add(aproof, PROOF_R_AT_ZETA)), R_MOD)
+        v := addmod(v, l_gamma, R_MOD)
+
+        w := mulmod(betazeta, coset_square, R_MOD)
+        w := addmod(w, calldataload(add(aproof, PROOF_O_AT_ZETA)), R_MOD)
+        w := addmod(w, l_gamma, R_MOD)
+
+        let s2 := mulmod(u, v, R_MOD)
+        s2 := mulmod(s2, w, R_MOD)
+        s2 := sub(R_MOD, s2)
+        s2 := mulmod(s2, l_alpha, R_MOD)
+        s2 := addmod(s2, mload(add(state, STATE_ALPHA_SQUARE_LAGRANGE_0)), R_MOD)
+
+        // at this stage:
+        // * s₁ = α*Z(μζ)(l(ζ)+β*s₁(ζ)+γ)*(r(ζ)+β*s₂(ζ)+γ)*β
+        // * s₂ = -α*(l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ) + α²*L₁(ζ)
+
+        compute_commitment_linearised_polynomial_ec(aproof, s1, s2)
+      }
+
+      /// @notice compute H₁ + ζᵐ⁺²*H₂ + ζ²⁽ᵐ⁺²⁾*H₃ and store the result at
+      /// state + state_folded_h
+      /// @param aproof pointer to the proof
+      function fold_h(aproof) {
+        let state := mload(0x40)
+        let n_plus_two := add(VK_DOMAIN_SIZE, 2)
+        let mPtr := add(mload(0x40), STATE_LAST_MEM)
+        let zeta_power_n_plus_two := pow(mload(add(state, STATE_ZETA)), n_plus_two, mPtr)
+        point_mul_calldata(add(state, STATE_FOLDED_H_X), add(aproof, PROOF_H_2_X), zeta_power_n_plus_two, mPtr)
+        point_add_calldata(add(state, STATE_FOLDED_H_X), add(state, STATE_FOLDED_H_X), add(aproof, PROOF_H_1_X), mPtr)
+        point_mul(add(state, STATE_FOLDED_H_X), add(state, STATE_FOLDED_H_X), zeta_power_n_plus_two, mPtr)
+        point_add_calldata(add(state, STATE_FOLDED_H_X), add(state, STATE_FOLDED_H_X), add(aproof, PROOF_H_0_X), mPtr)
+      }
+
+      /// @notice check that
+      ///	L(ζ)Qₗ(ζ)+r(ζ)Qᵣ(ζ)+R(ζ)L(ζ)Qₘ(ζ)+O(ζ)Qₒ(ζ)+Qₖ(ζ)+Σᵢqc'ᵢ(ζ)BsbCommitmentᵢ(ζ) +
+      ///  α*( Z(μζ)(l(ζ)+β*s₁(ζ)+γ)*(r(ζ)+β*s₂(ζ)+γ)*β*s₃(X)-Z(X)(l(ζ)+β*id_1(ζ)+γ)*(r(ζ)+β*id_2(ζ)+γ)*(o(ζ)+β*id_3(ζ)+γ) ) )
+      /// + α²*L₁(ζ) = 
+      /// (ζⁿ-1)H(ζ)
+      /// @param aproof pointer to the proof
+      function verify_quotient_poly_eval_at_zeta(aproof) {
+        let state := mload(0x40)
+
+        // (l(ζ)+β*s1(ζ)+γ)
+        let s1 := add(mload(0x40), STATE_LAST_MEM)
+        mstore(s1, mulmod(calldataload(add(aproof, PROOF_S1_AT_ZETA)), mload(add(state, STATE_BETA)), R_MOD))
+        mstore(s1, addmod(mload(s1), mload(add(state, STATE_GAMMA)), R_MOD))
+        mstore(s1, addmod(mload(s1), calldataload(add(aproof, PROOF_L_AT_ZETA)), R_MOD))
+
+        // (r(ζ)+β*s2(ζ)+γ)
+        let s2 := add(s1, 0x20)
+        mstore(s2, mulmod(calldataload(add(aproof, PROOF_S2_AT_ZETA)), mload(add(state, STATE_BETA)), R_MOD))
+        mstore(s2, addmod(mload(s2), mload(add(state, STATE_GAMMA)), R_MOD))
+        mstore(s2, addmod(mload(s2), calldataload(add(aproof, PROOF_R_AT_ZETA)), R_MOD))
+        // _s2 := mload(s2)
+
+        // (o(ζ)+γ)
+        let o := add(s1, 0x40)
+        mstore(o, addmod(calldataload(add(aproof, PROOF_O_AT_ZETA)), mload(add(state, STATE_GAMMA)), R_MOD))
+
+        //  α*(Z(μζ))*(l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*(o(ζ)+γ)
+        mstore(s1, mulmod(mload(s1), mload(s2), R_MOD))
+        mstore(s1, mulmod(mload(s1), mload(o), R_MOD))
+        mstore(s1, mulmod(mload(s1), mload(add(state, STATE_ALPHA)), R_MOD))
+        mstore(s1, mulmod(mload(s1), calldataload(add(aproof, PROOF_GRAND_PRODUCT_AT_ZETA_OMEGA)), R_MOD))
+
+        let computed_quotient := add(s1, 0x60)
+
+        // linearizedpolynomial + pi(zeta)
+        mstore(computed_quotient,addmod(calldataload(add(aproof, PROOF_LINEARISED_POLYNOMIAL_AT_ZETA)), mload(add(state, STATE_PI)), R_MOD))
+        mstore(computed_quotient, addmod(mload(computed_quotient), mload(s1), R_MOD))
+        mstore(computed_quotient,addmod(mload(computed_quotient), sub(R_MOD, mload(add(state, STATE_ALPHA_SQUARE_LAGRANGE_0))), R_MOD))
+        mstore(s2,mulmod(calldataload(add(aproof, PROOF_QUOTIENT_POLYNOMIAL_AT_ZETA)),mload(add(state, STATE_ZETA_POWER_N_MINUS_ONE)),R_MOD))
+
+        mstore(add(state, STATE_SUCCESS), eq(mload(computed_quotient), mload(s2)))
+      }
+
+      // BEGINNING utils math functions -------------------------------------------------
+      
+      /// @param dst pointer storing the result
+      /// @param p pointer to the first point
+      /// @param q pointer to the second point
+      /// @param mPtr pointer to free memory
+      function point_add(dst, p, q, mPtr) {
+        let state := mload(0x40)
+        mstore(mPtr, mload(p))
+        mstore(add(mPtr, 0x20), mload(add(p, 0x20)))
+        mstore(add(mPtr, 0x40), mload(q))
+        mstore(add(mPtr, 0x60), mload(add(q, 0x20)))
+        let l_success := staticcall(gas(),6,mPtr,0x80,dst,0x40)
+        if iszero(l_success) {
+          error_ec_op()
+        }
+      }
+
+      /// @param dst pointer storing the result
+      /// @param p pointer to the first point (calldata)
+      /// @param q pointer to the second point (calladata)
+      /// @param mPtr pointer to free memory
+      function point_add_calldata(dst, p, q, mPtr) {
+        let state := mload(0x40)
+        mstore(mPtr, mload(p))
+        mstore(add(mPtr, 0x20), mload(add(p, 0x20)))
+        mstore(add(mPtr, 0x40), calldataload(q))
+        mstore(add(mPtr, 0x60), calldataload(add(q, 0x20)))
+        let l_success := staticcall(gas(), 6, mPtr, 0x80, dst, 0x40)
+        if iszero(l_success) {
+          error_ec_op()
+        }
+      }
+
+      /// @parma dst pointer storing the result
+      /// @param src pointer to a point on Bn254(𝔽_p)
+      /// @param s scalar
+      /// @param mPtr free memory
+      function point_mul(dst,src,s, mPtr) {
+        let state := mload(0x40)
+        mstore(mPtr,mload(src))
+        mstore(add(mPtr,0x20),mload(add(src,0x20)))
+        mstore(add(mPtr,0x40),s)
+        let l_success := staticcall(gas(),7,mPtr,0x60,dst,0x40)
+        if iszero(l_success) {
+          error_ec_op()
+        }
+      }
+
+      /// @parma dst pointer storing the result
+      /// @param src pointer to a point on Bn254(𝔽_p) on calldata
+      /// @param s scalar
+      /// @param mPtr free memory
+      function point_mul_calldata(dst, src, s, mPtr) {
+        let state := mload(0x40)
+        mstore(mPtr, calldataload(src))
+        mstore(add(mPtr, 0x20), calldataload(add(src, 0x20)))
+        mstore(add(mPtr, 0x40), s)
+        let l_success := staticcall(gas(), 7, mPtr, 0x60, dst, 0x40)
+        if iszero(l_success) {
+          error_ec_op()
+        }
+      }
+
+      /// @notice dst <- dst + [s]src (Elliptic curve)
+      /// @param dst pointer accumulator point storing the result
+      /// @param src pointer to the point to multiply and add
+      /// @param s scalar
+      /// @param mPtr free memory
+      function point_acc_mul(dst,src,s, mPtr) {
+        let state := mload(0x40)
+        mstore(mPtr,mload(src))
+        mstore(add(mPtr,0x20),mload(add(src,0x20)))
+        mstore(add(mPtr,0x40),s)
+        let l_success := staticcall(gas(),7,mPtr,0x60,mPtr,0x40)
+        mstore(add(mPtr,0x40),mload(dst))
+        mstore(add(mPtr,0x60),mload(add(dst,0x20)))
+        l_success := and(l_success, staticcall(gas(),6,mPtr,0x80,dst, 0x40))
+        if iszero(l_success) {
+          error_ec_op()
+        }
+      }
+
+      /// @notice dst <- dst + [s]src (Elliptic curve)
+      /// @param dst pointer accumulator point storing the result
+      /// @param src pointer to the point to multiply and add (on calldata)
+      /// @param s scalar
+      /// @mPtr free memory
+      function point_acc_mul_calldata(dst, src, s, mPtr) {
+        let state := mload(0x40)
+        mstore(mPtr, calldataload(src))
+        mstore(add(mPtr, 0x20), calldataload(add(src, 0x20)))
+        mstore(add(mPtr, 0x40), s)
+        let l_success := staticcall(gas(), 7, mPtr, 0x60, mPtr, 0x40)
+        mstore(add(mPtr, 0x40), mload(dst))
+        mstore(add(mPtr, 0x60), mload(add(dst, 0x20)))
+        l_success := and(l_success, staticcall(gas(), 6, mPtr, 0x80, dst, 0x40))
+        if iszero(l_success) {
+          error_ec_op()
+        }
+      }
+
+      /// @notice dst <- dst + src*s (Fr) dst,src are addresses, s is a value
+      /// @param dst pointer storing the result
+      /// @param src pointer to the scalar to multiply and add (on calldata)
+      /// @param s scalar
+      function fr_acc_mul_calldata(dst, src, s) {
+        let tmp :=  mulmod(calldataload(src), s, R_MOD)
+        mstore(dst, addmod(mload(dst), tmp, R_MOD))
+      }
+
+      /// @param x element to exponentiate
+      /// @param e exponent
+      /// @param mPtr free memory
+      /// @return res x ** e mod r
+      function pow(x, e, mPtr)->res {
+        mstore(mPtr, 0x20)
+        mstore(add(mPtr, 0x20), 0x20)
+        mstore(add(mPtr, 0x40), 0x20)
+        mstore(add(mPtr, 0x60), x)
+        mstore(add(mPtr, 0x80), e)
+        mstore(add(mPtr, 0xa0), R_MOD)
+        let check_staticcall := staticcall(gas(),0x05,mPtr,0xc0,mPtr,0x20)
+        if eq(check_staticcall, 0) {
+          error_verify()
+        }
+        res := mload(mPtr)
+      }
+    }
+  }
+}
+`
+
+// MarshalSolidity converts a proof to a byte array that can be used in a
+// Solidity contract.
+func (proof *Proof) MarshalSolidity() []byte {
+
+	res := make([]byte, 0, 1024)
+
+	// uint256 l_com_x;
+	// uint256 l_com_y;
+	// uint256 r_com_x;
+	// uint256 r_com_y;
+	// uint256 o_com_x;
+	// uint256 o_com_y;
+	var tmp64 [64]byte
+	for i := 0; i < 3; i++ {
+		tmp64 = proof.LRO[i].RawBytes()
+		res = append(res, tmp64[:]...)
+	}
+
+	// uint256 h_0_x;
+	// uint256 h_0_y;
+	// uint256 h_1_x;
+	// uint256 h_1_y;
+	// uint256 h_2_x;
+	// uint256 h_2_y;
+	for i := 0; i < 3; i++ {
+		tmp64 = proof.H[i].RawBytes()
+		res = append(res, tmp64[:]...)
+	}
+	var tmp32 [32]byte
+
+	// uint256 l_at_zeta;
+	// uint256 r_at_zeta;
+	// uint256 o_at_zeta;
+	// uint256 s1_at_zeta;
+	// uint256 s2_at_zeta;
+	for i := 2; i < 7; i++ {
+		tmp32 = proof.BatchedProof.ClaimedValues[i].Bytes()
+		res = append(res, tmp32[:]...)
+	}
+
+	// uint256 grand_product_commitment_x;
+	// uint256 grand_product_commitment_y;
+	tmp64 = proof.Z.RawBytes()
+	res = append(res, tmp64[:]...)
+
+	// uint256 grand_product_at_zeta_omega;
+	tmp32 = proof.ZShiftedOpening.ClaimedValue.Bytes()
+	res = append(res, tmp32[:]...)
+
+	// uint256 quotient_polynomial_at_zeta;
+	// uint256 linearization_polynomial_at_zeta;
+	tmp32 = proof.BatchedProof.ClaimedValues[0].Bytes()
+	res = append(res, tmp32[:]...)
+	tmp32 = proof.BatchedProof.ClaimedValues[1].Bytes()
+	res = append(res, tmp32[:]...)
+
+	// uint256 opening_at_zeta_proof_x;
+	// uint256 opening_at_zeta_proof_y;
+	tmp64 = proof.BatchedProof.H.RawBytes()
+	res = append(res, tmp64[:]...)
+
+	// uint256 opening_at_zeta_omega_proof_x;
+	// uint256 opening_at_zeta_omega_proof_y;
+	tmp64 = proof.ZShiftedOpening.H.RawBytes()
+	res = append(res, tmp64[:]...)
+
+	// uint256[] selector_commit_api_at_zeta;
+	// uint256[] wire_committed_commitments;
+	if len(proof.Bsb22Commitments) > 0 {
+		for i := 0; i < len(proof.Bsb22Commitments); i++ {
+			tmp32 = proof.BatchedProof.ClaimedValues[7+i].Bytes()
+			res = append(res, tmp32[:]...)
+		}
+
+		for _, bc := range proof.Bsb22Commitments {
+			tmp64 = bc.RawBytes()
+			res = append(res, tmp64[:]...)
+		}
+	}
+
+	return res
+}
diff --git a/backend/plonk/bn254/icicle/unmarshal.go b/backend/plonk/bn254/icicle/unmarshal.go
new file mode 100644
index 0000000000..9edd6a4730
--- /dev/null
+++ b/backend/plonk/bn254/icicle/unmarshal.go
@@ -0,0 +1,88 @@
+package icicle_bn254
+
+import (
+	"github.com/consensys/gnark-crypto/ecc/bn254"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+)
+
+func UnmarshalSolidity(s []byte, nbCommits int) Proof {
+
+	var proof Proof
+	offset := 0
+	point_size := 64
+	fr_size := 32
+	proof.BatchedProof.ClaimedValues = make([]fr.Element, 7+nbCommits)
+	proof.Bsb22Commitments = make([]bn254.G1Affine, nbCommits)
+
+	// uint256 l_com_x;
+	// uint256 l_com_y;
+	// uint256 r_com_x;
+	// uint256 r_com_y;
+	// uint256 o_com_x;
+	// uint256 o_com_y;
+	for i := 0; i < 3; i++ {
+		proof.LRO[i].Unmarshal(s[offset : offset+point_size])
+		offset += point_size
+	}
+
+	// uint256 h_0_x;
+	// uint256 h_0_y;
+	// uint256 h_1_x;
+	// uint256 h_1_y;
+	// uint256 h_2_x;
+	// uint256 h_2_y;
+	for i := 0; i < 3; i++ {
+		proof.H[i].Unmarshal(s[offset : offset+point_size])
+		offset += point_size
+	}
+
+	// uint256 l_at_zeta;
+	// uint256 r_at_zeta;
+	// uint256 o_at_zeta;
+	// uint256 s1_at_zeta;
+	// uint256 s2_at_zeta;
+	for i := 2; i < 7; i++ {
+		proof.BatchedProof.ClaimedValues[i].SetBytes(s[offset : offset+fr_size])
+		offset += fr_size
+	}
+
+	// uint256 grand_product_commitment_x;
+	// uint256 grand_product_commitment_y;
+	proof.Z.Unmarshal(s[offset : offset+point_size])
+	offset += point_size
+
+	// uint256 grand_product_at_zeta_omega;
+	proof.ZShiftedOpening.ClaimedValue.SetBytes(s[offset : offset+fr_size])
+	offset += fr_size
+
+	// uint256 quotient_polynomial_at_zeta;
+	// uint256 linearization_polynomial_at_zeta;
+	proof.BatchedProof.ClaimedValues[0].SetBytes(s[offset : offset+fr_size])
+	offset += fr_size
+	proof.BatchedProof.ClaimedValues[1].SetBytes(s[offset : offset+fr_size])
+	offset += fr_size
+
+	// uint256 opening_at_zeta_proof_x;
+	// uint256 opening_at_zeta_proof_y;
+	proof.BatchedProof.H.Unmarshal(s[offset : offset+point_size])
+	offset += point_size
+
+	// uint256 opening_at_zeta_omega_proof_x;
+	// uint256 opening_at_zeta_omega_proof_y;
+	proof.ZShiftedOpening.H.Unmarshal(s[offset : offset+point_size])
+	offset += point_size
+
+	// uint256[] selector_commit_api_at_zeta;
+	// uint256[] wire_committed_commitments;
+	for i := 0; i < nbCommits; i++ {
+		proof.BatchedProof.ClaimedValues[7+i].SetBytes(s[offset : offset+fr_size])
+		offset += fr_size
+	}
+
+	for i := 0; i < nbCommits; i++ {
+		proof.Bsb22Commitments[i].Unmarshal(s[offset : offset+point_size])
+		offset += point_size
+	}
+
+	return proof
+}
diff --git a/backend/plonk/bn254/icicle/verify.go b/backend/plonk/bn254/icicle/verify.go
new file mode 100644
index 0000000000..1b7902282e
--- /dev/null
+++ b/backend/plonk/bn254/icicle/verify.go
@@ -0,0 +1,401 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by gnark DO NOT EDIT
+
+package icicle_bn254
+
+import (
+	"errors"
+	"fmt"
+	"io"
+	"math/big"
+	"text/template"
+	"time"
+
+	"github.com/consensys/gnark-crypto/ecc"
+
+	curve "github.com/consensys/gnark-crypto/ecc/bn254"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/fp"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
+	"github.com/consensys/gnark-crypto/ecc/bn254/fr/hash_to_field"
+
+	"github.com/consensys/gnark-crypto/ecc/bn254/kzg"
+	fiatshamir "github.com/consensys/gnark-crypto/fiat-shamir"
+	"github.com/consensys/gnark/backend"
+	"github.com/consensys/gnark/logger"
+)
+
+var (
+	errWrongClaimedQuotient = errors.New("claimed quotient is not as expected")
+)
+
+func Verify(proof *Proof, vk *VerifyingKey, publicWitness fr.Vector, opts ...backend.VerifierOption) error {
+	log := logger.Logger().With().Str("curve", "bn254").Str("backend", "plonk").Logger()
+	start := time.Now()
+	cfg, err := backend.NewVerifierConfig(opts...)
+	if err != nil {
+		return fmt.Errorf("create backend config: %w", err)
+	}
+
+	if len(proof.Bsb22Commitments) != len(vk.Qcp) {
+		return errors.New("BSB22 Commitment number mismatch")
+	}
+
+	// transcript to derive the challenge
+	fs := fiatshamir.NewTranscript(cfg.ChallengeHash, "gamma", "beta", "alpha", "zeta")
+
+	// The first challenge is derived using the public data: the commitments to the permutation,
+	// the coefficients of the circuit, and the public inputs.
+	// derive gamma from the Comm(blinded cl), Comm(blinded cr), Comm(blinded co)
+	if err := bindPublicData(&fs, "gamma", vk, publicWitness); err != nil {
+		return err
+	}
+	gamma, err := deriveRandomness(&fs, "gamma", &proof.LRO[0], &proof.LRO[1], &proof.LRO[2])
+	if err != nil {
+		return err
+	}
+
+	// derive beta from Comm(l), Comm(r), Comm(o)
+	beta, err := deriveRandomness(&fs, "beta")
+	if err != nil {
+		return err
+	}
+
+	// derive alpha from Comm(l), Comm(r), Comm(o), Com(Z), Bsb22Commitments
+	alphaDeps := make([]*curve.G1Affine, len(proof.Bsb22Commitments)+1)
+	for i := range proof.Bsb22Commitments {
+		alphaDeps[i] = &proof.Bsb22Commitments[i]
+	}
+	alphaDeps[len(alphaDeps)-1] = &proof.Z
+	alpha, err := deriveRandomness(&fs, "alpha", alphaDeps...)
+	if err != nil {
+		return err
+	}
+
+	// derive zeta, the point of evaluation
+	zeta, err := deriveRandomness(&fs, "zeta", &proof.H[0], &proof.H[1], &proof.H[2])
+	if err != nil {
+		return err
+	}
+
+	// evaluation of Z=Xⁿ⁻¹ at ζ
+	var zetaPowerM, zzeta fr.Element
+	var bExpo big.Int
+	one := fr.One()
+	bExpo.SetUint64(vk.Size)
+	zetaPowerM.Exp(zeta, &bExpo)
+	zzeta.Sub(&zetaPowerM, &one)
+
+	// compute PI = ∑_{i<n} Lᵢ*wᵢ
+	// TODO use batch inversion
+	var pi, lagrangeOne fr.Element
+	{
+		var den, xiLi fr.Element
+		lagrange := zzeta // ζⁿ⁻¹
+		wPowI := fr.One()
+		den.Sub(&zeta, &wPowI)
+		lagrange.Div(&lagrange, &den).Mul(&lagrange, &vk.SizeInv) // (1/n)(ζⁿ-1)/(ζ-1)
+		lagrangeOne.Set(&lagrange)                                // save it for later
+		for i := 0; i < len(publicWitness); i++ {
+
+			xiLi.Mul(&lagrange, &publicWitness[i])
+			pi.Add(&pi, &xiLi)
+
+			// use Lᵢ₊₁ = w×Lᵢ(ζ-wⁱ)/(ζ-wⁱ⁺¹)
+			if i+1 != len(publicWitness) {
+				lagrange.Mul(&lagrange, &vk.Generator).
+					Mul(&lagrange, &den)
+				wPowI.Mul(&wPowI, &vk.Generator)
+				den.Sub(&zeta, &wPowI)
+				lagrange.Div(&lagrange, &den)
+			}
+		}
+
+		if cfg.HashToFieldFn == nil {
+			cfg.HashToFieldFn = hash_to_field.New([]byte("BSB22-Plonk"))
+		}
+		var hashBts []byte
+		var hashedCmt fr.Element
+		nbBuf := fr.Bytes
+		if cfg.HashToFieldFn.Size() < fr.Bytes {
+			nbBuf = cfg.HashToFieldFn.Size()
+		}
+		for i := range vk.CommitmentConstraintIndexes {
+			cfg.HashToFieldFn.Write(proof.Bsb22Commitments[i].Marshal())
+			hashBts = cfg.HashToFieldFn.Sum(hashBts[0:])
+			cfg.HashToFieldFn.Reset()
+			hashedCmt.SetBytes(hashBts[:nbBuf])
+
+			// Computing L_{CommitmentIndex}
+
+			wPowI.Exp(vk.Generator, big.NewInt(int64(vk.NbPublicVariables)+int64(vk.CommitmentConstraintIndexes[i])))
+			den.Sub(&zeta, &wPowI) // ζ-wⁱ
+
+			lagrange.SetOne().
+				Sub(&zeta, &lagrange).       // ζ-1
+				Mul(&lagrange, &wPowI).      // wⁱ(ζ-1)
+				Div(&lagrange, &den).        // wⁱ(ζ-1)/(ζ-wⁱ)
+				Mul(&lagrange, &lagrangeOne) // wⁱ/n (ζⁿ-1)/(ζ-wⁱ)
+
+			xiLi.Mul(&lagrange, &hashedCmt)
+			pi.Add(&pi, &xiLi)
+		}
+	}
+
+	// linearizedpolynomial + pi(ζ) + α*(Z(μζ))*(l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*(o(ζ)+γ) - α²*L₁(ζ)
+	var _s1, _s2, _o, alphaSquareLagrange fr.Element
+
+	zu := proof.ZShiftedOpening.ClaimedValue
+
+	claimedQuotient := proof.BatchedProof.ClaimedValues[0]
+	linearizedPolynomialZeta := proof.BatchedProof.ClaimedValues[1]
+	l := proof.BatchedProof.ClaimedValues[2]
+	r := proof.BatchedProof.ClaimedValues[3]
+	o := proof.BatchedProof.ClaimedValues[4]
+	s1 := proof.BatchedProof.ClaimedValues[5]
+	s2 := proof.BatchedProof.ClaimedValues[6]
+
+	_s1.Mul(&s1, &beta).Add(&_s1, &l).Add(&_s1, &gamma) // (l(ζ)+β*s1(ζ)+γ)
+	_s2.Mul(&s2, &beta).Add(&_s2, &r).Add(&_s2, &gamma) // (r(ζ)+β*s2(ζ)+γ)
+	_o.Add(&o, &gamma)                                  // (o(ζ)+γ)
+
+	_s1.Mul(&_s1, &_s2).
+		Mul(&_s1, &_o).
+		Mul(&_s1, &alpha).
+		Mul(&_s1, &zu) //  α*(Z(μζ))*(l(ζ)+β*s1(ζ)+γ)*(r(ζ)+β*s2(ζ)+γ)*(o(ζ)+γ)
+
+	alphaSquareLagrange.Mul(&lagrangeOne, &alpha).
+		Mul(&alphaSquareLagrange, &alpha) // α²*L₁(ζ)
+
+	linearizedPolynomialZeta.
+		Add(&linearizedPolynomialZeta, &pi).                 // linearizedpolynomial + pi(zeta)
+		Add(&linearizedPolynomialZeta, &_s1).                // linearizedpolynomial+pi(zeta)+α*(Z(μζ))*(l(ζ)+s1(ζ)+γ)*(r(ζ)+s2(ζ)+γ)*(o(ζ)+γ)
+		Sub(&linearizedPolynomialZeta, &alphaSquareLagrange) // linearizedpolynomial+pi(zeta)+α*(Z(μζ))*(l(ζ)+s1(ζ)+γ)*(r(ζ)+s2(ζ)+γ)*(o(ζ)+γ)-α²*L₁(ζ)
+
+	// Compute H(ζ) using the previous result: H(ζ) = prev_result/(ζⁿ-1)
+	var zetaPowerMMinusOne fr.Element
+	zetaPowerMMinusOne.Sub(&zetaPowerM, &one)
+	linearizedPolynomialZeta.Div(&linearizedPolynomialZeta, &zetaPowerMMinusOne)
+
+	// check that H(ζ) is as claimed
+	if !claimedQuotient.Equal(&linearizedPolynomialZeta) {
+		return errWrongClaimedQuotient
+	}
+
+	// compute the folded commitment to H: Comm(h₁) + ζᵐ⁺²*Comm(h₂) + ζ²⁽ᵐ⁺²⁾*Comm(h₃)
+	mPlusTwo := big.NewInt(int64(vk.Size) + 2)
+	var zetaMPlusTwo fr.Element
+	zetaMPlusTwo.Exp(zeta, mPlusTwo)
+	var zetaMPlusTwoBigInt big.Int
+	zetaMPlusTwo.BigInt(&zetaMPlusTwoBigInt)
+	foldedH := proof.H[2]
+	foldedH.ScalarMultiplication(&foldedH, &zetaMPlusTwoBigInt)
+	foldedH.Add(&foldedH, &proof.H[1])
+	foldedH.ScalarMultiplication(&foldedH, &zetaMPlusTwoBigInt)
+	foldedH.Add(&foldedH, &proof.H[0])
+
+	// Compute the commitment to the linearized polynomial
+	// linearizedPolynomialDigest =
+	// 		l(ζ)*ql+r(ζ)*qr+r(ζ)l(ζ)*qm+o(ζ)*qo+qk+Σᵢqc'ᵢ(ζ)*BsbCommitmentᵢ +
+	// 		α*( Z(μζ)(l(ζ)+β*s₁(ζ)+γ)*(r(ζ)+β*s₂(ζ)+γ)*s₃(X)-Z(X)(l(ζ)+β*id_1(ζ)+γ)*(r(ζ)+β*id_2(ζ)+γ)*(o(ζ)+β*id_3(ζ)+γ) ) +
+	// 		α²*L₁(ζ)*Z
+	// first part: individual constraints
+	var rl fr.Element
+	rl.Mul(&l, &r)
+
+	var linearizedPolynomialDigest curve.G1Affine
+
+	// second part: α*( Z(μζ)(l(ζ)+β*s₁(ζ)+γ)*(r(ζ)+β*s₂(ζ)+γ)*β*s₃(X)-Z(X)(l(ζ)+β*id_1(ζ)+γ)*(r(ζ)+β*id_2(ζ)+γ)*(o(ζ)+β*id_3(ζ)+γ) ) )
+
+	var u, v, w, cosetsquare fr.Element
+	u.Mul(&zu, &beta)
+	v.Mul(&beta, &s1).Add(&v, &l).Add(&v, &gamma)
+	w.Mul(&beta, &s2).Add(&w, &r).Add(&w, &gamma)
+	_s1.Mul(&u, &v).Mul(&_s1, &w).Mul(&_s1, &alpha) // α*Z(μζ)(l(ζ)+β*s₁(ζ)+γ)*(r(ζ)+β*s₂(ζ)+γ)*β
+
+	cosetsquare.Square(&vk.CosetShift)
+	u.Mul(&beta, &zeta).Add(&u, &l).Add(&u, &gamma)                         // (l(ζ)+β*ζ+γ)
+	v.Mul(&beta, &zeta).Mul(&v, &vk.CosetShift).Add(&v, &r).Add(&v, &gamma) // (r(ζ)+β*μ*ζ+γ)
+	w.Mul(&beta, &zeta).Mul(&w, &cosetsquare).Add(&w, &o).Add(&w, &gamma)   // (o(ζ)+β*μ²*ζ+γ)
+	_s2.Mul(&u, &v).Mul(&_s2, &w).Neg(&_s2)                                 // -(l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ)
+
+	// note since third part =  α²*L₁(ζ)*Z
+	_s2.Mul(&_s2, &alpha).Add(&_s2, &alphaSquareLagrange) // -α*(l(ζ)+β*ζ+γ)*(r(ζ)+β*u*ζ+γ)*(o(ζ)+β*u²*ζ+γ) + α²*L₁(ζ)
+
+	points := append(proof.Bsb22Commitments,
+		vk.Ql, vk.Qr, vk.Qm, vk.Qo, vk.Qk, // first part
+		vk.S[2], proof.Z, // second & third part
+	)
+
+	qC := make([]fr.Element, len(proof.Bsb22Commitments))
+	copy(qC, proof.BatchedProof.ClaimedValues[7:])
+	scalars := append(qC,
+		l, r, rl, o, one, /* TODO Perf @Tabaie Consider just adding Qk instead */ // first part
+		_s1, _s2, // second & third part
+	)
+	if _, err := linearizedPolynomialDigest.MultiExp(points, scalars, ecc.MultiExpConfig{}); err != nil {
+		return err
+	}
+
+	// Fold the first proof
+	digestsToFold := make([]curve.G1Affine, len(vk.Qcp)+7)
+	copy(digestsToFold[7:], vk.Qcp)
+	digestsToFold[0] = foldedH
+	digestsToFold[1] = linearizedPolynomialDigest
+	digestsToFold[2] = proof.LRO[0]
+	digestsToFold[3] = proof.LRO[1]
+	digestsToFold[4] = proof.LRO[2]
+	digestsToFold[5] = vk.S[0]
+	digestsToFold[6] = vk.S[1]
+	foldedProof, foldedDigest, err := kzg.FoldProof(
+		digestsToFold,
+		&proof.BatchedProof,
+		zeta,
+		cfg.KZGFoldingHash,
+		zu.Marshal(),
+	)
+	if err != nil {
+		return err
+	}
+
+	// Batch verify
+	var shiftedZeta fr.Element
+	shiftedZeta.Mul(&zeta, &vk.Generator)
+	err = kzg.BatchVerifyMultiPoints([]kzg.Digest{
+		foldedDigest,
+		proof.Z,
+	},
+		[]kzg.OpeningProof{
+			foldedProof,
+			proof.ZShiftedOpening,
+		},
+		[]fr.Element{
+			zeta,
+			shiftedZeta,
+		},
+		vk.Kzg,
+	)
+
+	log.Debug().Dur("took", time.Since(start)).Msg("verifier done")
+
+	return err
+}
+
+func bindPublicData(fs *fiatshamir.Transcript, challenge string, vk *VerifyingKey, publicInputs []fr.Element) error {
+
+	// permutation
+	if err := fs.Bind(challenge, vk.S[0].Marshal()); err != nil {
+		return err
+	}
+	if err := fs.Bind(challenge, vk.S[1].Marshal()); err != nil {
+		return err
+	}
+	if err := fs.Bind(challenge, vk.S[2].Marshal()); err != nil {
+		return err
+	}
+
+	// coefficients
+	if err := fs.Bind(challenge, vk.Ql.Marshal()); err != nil {
+		return err
+	}
+	if err := fs.Bind(challenge, vk.Qr.Marshal()); err != nil {
+		return err
+	}
+	if err := fs.Bind(challenge, vk.Qm.Marshal()); err != nil {
+		return err
+	}
+	if err := fs.Bind(challenge, vk.Qo.Marshal()); err != nil {
+		return err
+	}
+	if err := fs.Bind(challenge, vk.Qk.Marshal()); err != nil {
+		return err
+	}
+	for i := range vk.Qcp {
+		if err := fs.Bind(challenge, vk.Qcp[i].Marshal()); err != nil {
+			return err
+		}
+	}
+
+	// public inputs
+	for i := 0; i < len(publicInputs); i++ {
+		if err := fs.Bind(challenge, publicInputs[i].Marshal()); err != nil {
+			return err
+		}
+	}
+
+	return nil
+
+}
+
+func deriveRandomness(fs *fiatshamir.Transcript, challenge string, points ...*curve.G1Affine) (fr.Element, error) {
+
+	var buf [curve.SizeOfG1AffineUncompressed]byte
+	var r fr.Element
+
+	for _, p := range points {
+		buf = p.RawBytes()
+		if err := fs.Bind(challenge, buf[:]); err != nil {
+			return r, err
+		}
+	}
+
+	b, err := fs.ComputeChallenge(challenge)
+	if err != nil {
+		return r, err
+	}
+	r.SetBytes(b)
+	return r, nil
+}
+
+// ExportSolidity exports the verifying key to a solidity smart contract.
+//
+// See https://github.com/ConsenSys/gnark-tests for example usage.
+//
+// Code has not been audited and is provided as-is, we make no guarantees or warranties to its safety and reliability.
+func (vk *VerifyingKey) ExportSolidity(w io.Writer) error {
+	funcMap := template.FuncMap{
+		"hex": func(i int) string {
+			return fmt.Sprintf("0x%x", i)
+		},
+		"mul": func(a, b int) int {
+			return a * b
+		},
+		"inc": func(i int) int {
+			return i + 1
+		},
+		"frstr": func(x fr.Element) string {
+			// we use big.Int to always get a positive string.
+			// not the most efficient hack, but it works better for .sol generation.
+			bv := new(big.Int)
+			x.BigInt(bv)
+			return bv.String()
+		},
+		"fpstr": func(x fp.Element) string {
+			bv := new(big.Int)
+			x.BigInt(bv)
+			return bv.String()
+		},
+		"add": func(i, j int) int {
+			return i + j
+		},
+	}
+
+	t, err := template.New("t").Funcs(funcMap).Parse(tmplSolidityVerifier)
+	if err != nil {
+		return err
+	}
+	return t.Execute(w, vk)
+}
diff --git a/backend/plonk/bn254/prove.go b/backend/plonk/bn254/prove.go
index 1a237418d6..f117c88e95 100644
--- a/backend/plonk/bn254/prove.go
+++ b/backend/plonk/bn254/prove.go
@@ -113,7 +113,7 @@ type Proof struct {
 }
 
 func Prove(spr *cs.SparseR1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...backend.ProverOption) (*Proof, error) {
-
+	
 	log := logger.Logger().With().
 		Str("curve", spr.CurveID().String()).
 		Int("nbConstraints", spr.GetNbConstraints()).
diff --git a/backend/plonk/plonk.go b/backend/plonk/plonk.go
index 2fd11b18aa..3c5cceac7b 100644
--- a/backend/plonk/plonk.go
+++ b/backend/plonk/plonk.go
@@ -41,6 +41,7 @@ import (
 	plonk_bls24315 "github.com/consensys/gnark/backend/plonk/bls24-315"
 	plonk_bls24317 "github.com/consensys/gnark/backend/plonk/bls24-317"
 	plonk_bn254 "github.com/consensys/gnark/backend/plonk/bn254"
+	icicle_bn254 "github.com/consensys/gnark/backend/plonk/bn254/icicle"
 	plonk_bw6633 "github.com/consensys/gnark/backend/plonk/bw6-633"
 	plonk_bw6761 "github.com/consensys/gnark/backend/plonk/bw6-761"
 
@@ -100,6 +101,9 @@ func Setup(ccs constraint.ConstraintSystem, kzgSrs kzg.SRS) (ProvingKey, Verifyi
 
 	switch tccs := ccs.(type) {
 	case *cs_bn254.SparseR1CS:
+		if icicle_bn254.HasIcicle {
+			return icicle_bn254.Setup(tccs, *kzgSrs.(*kzg_bn254.SRS))
+		}
 		return plonk_bn254.Setup(tccs, *kzgSrs.(*kzg_bn254.SRS))
 	case *cs_bls12381.SparseR1CS:
 		return plonk_bls12381.Setup(tccs, *kzgSrs.(*kzg_bls12381.SRS))
@@ -129,6 +133,9 @@ func Prove(ccs constraint.ConstraintSystem, pk ProvingKey, fullWitness witness.W
 
 	switch tccs := ccs.(type) {
 	case *cs_bn254.SparseR1CS:
+		if icicle_bn254.HasIcicle {
+			return icicle_bn254.Prove(tccs, pk.(*icicle_bn254.ProvingKey), fullWitness, opts...)
+		}
 		return plonk_bn254.Prove(tccs, pk.(*plonk_bn254.ProvingKey), fullWitness, opts...)
 
 	case *cs_bls12381.SparseR1CS:
@@ -159,6 +166,13 @@ func Verify(proof Proof, vk VerifyingKey, publicWitness witness.Witness, opts ..
 
 	switch _proof := proof.(type) {
 
+	case *icicle_bn254.Proof:
+		w, ok := publicWitness.Vector().(fr_bn254.Vector)
+		if !ok {
+			return witness.ErrInvalidWitness
+		}
+		return icicle_bn254.Verify(_proof, vk.(*icicle_bn254.VerifyingKey), w, opts...)
+
 	case *plonk_bn254.Proof:
 		w, ok := publicWitness.Vector().(fr_bn254.Vector)
 		if !ok {
diff --git a/go.mod b/go.mod
index 181f81b1f2..27ed2d46eb 100644
--- a/go.mod
+++ b/go.mod
@@ -6,31 +6,31 @@ require (
 	github.com/bits-and-blooms/bitset v1.8.0
 	github.com/blang/semver/v4 v4.0.0
 	github.com/consensys/bavard v0.1.13
-	github.com/consensys/gnark-crypto v0.12.2-0.20231023220848-538dff926c15
+	github.com/consensys/gnark-crypto v0.12.2-0.20231208203441-d4eab6ddd2af
 	github.com/fxamacker/cbor/v2 v2.5.0
 	github.com/google/go-cmp v0.5.9
 	github.com/google/pprof v0.0.0-20230817174616-7a8ec2ada47b
+	//github.com/ingonyama-zk/iciclegnark v0.1.0.
+	github.com/ingonyama-zk/iciclegnark v0.1.2-0.20240207160517-9558f9ad3baf
 	github.com/leanovate/gopter v0.2.9
 	github.com/rs/zerolog v1.30.0
 	github.com/stretchr/testify v1.8.4
-	golang.org/x/crypto v0.12.0
+	golang.org/x/crypto v0.17.0
 	golang.org/x/exp v0.0.0-20230817173708-d852ddb80c63
 	golang.org/x/sync v0.3.0
 )
 
-require (
-	github.com/kr/text v0.2.0 // indirect
-	github.com/rogpeppe/go-internal v1.11.0 // indirect
-)
-
 require (
 	github.com/davecgh/go-spew v1.1.1 // indirect
+	github.com/ingonyama-zk/icicle v0.1.0 // indirect
 	github.com/mattn/go-colorable v0.1.13 // indirect
 	github.com/mattn/go-isatty v0.0.19 // indirect
 	github.com/mmcloughlin/addchain v0.4.0 // indirect
 	github.com/pmezard/go-difflib v1.0.0 // indirect
 	github.com/x448/float16 v0.8.4 // indirect
-	golang.org/x/sys v0.11.0 // indirect
+	golang.org/x/sys v0.15.0 // indirect
 	gopkg.in/yaml.v3 v3.0.1 // indirect
 	rsc.io/tmplfunc v0.0.3 // indirect
 )
+
+replace github.com/consensys/gnark-crypto v0.12.2-0.20231208203441-d4eab6ddd2af => github.com/ingonyama-zk/gnark-crypto v0.0.0-20240207133741-7e078bb23cf3
diff --git a/go.sum b/go.sum
index 93b1eff8c0..b4778a24e6 100644
--- a/go.sum
+++ b/go.sum
@@ -4,12 +4,7 @@ github.com/blang/semver/v4 v4.0.0 h1:1PFHFE6yCCTv8C1TeyNNarDzntLi7wMI5i/pzqYIsAM
 github.com/blang/semver/v4 v4.0.0/go.mod h1:IbckMUScFkM3pff0VJDNKRiT6TG/YpiHIM2yvyW5YoQ=
 github.com/consensys/bavard v0.1.13 h1:oLhMLOFGTLdlda/kma4VOJazblc7IM5y5QPd2A/YjhQ=
 github.com/consensys/bavard v0.1.13/go.mod h1:9ItSMtA/dXMAiL7BG6bqW2m3NdSEObYWoH223nGHukI=
-github.com/consensys/gnark-crypto v0.12.2-0.20231017134050-6652a8b98254 h1:21iGClQpXhJ0PA6lRBpRbpmJFtfdV7R9QMCsQJ3A9mA=
-github.com/consensys/gnark-crypto v0.12.2-0.20231017134050-6652a8b98254/go.mod h1:v2Gy7L/4ZRosZ7Ivs+9SfUDr0f5UlG+EM5t7MPHiLuY=
-github.com/consensys/gnark-crypto v0.12.2-0.20231023220848-538dff926c15 h1:fu5ienFKWWqrfMPbWnhw4zfIFZW3pzVIbv3KtASymbU=
-github.com/consensys/gnark-crypto v0.12.2-0.20231023220848-538dff926c15/go.mod h1:v2Gy7L/4ZRosZ7Ivs+9SfUDr0f5UlG+EM5t7MPHiLuY=
 github.com/coreos/go-systemd/v22 v22.5.0/go.mod h1:Y58oyj3AT4RCenI/lSvhwexgC+NSVTIJ3seZv2GcEnc=
-github.com/creack/pty v1.1.9/go.mod h1:oKZEueFk5CKHvIhNR5MUki03XCEU+Q6VDXinZuGJ33E=
 github.com/davecgh/go-spew v1.1.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c=
 github.com/davecgh/go-spew v1.1.1/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=
 github.com/fxamacker/cbor/v2 v2.5.0 h1:oHsG0V/Q6E/wqTS2O1Cozzsy69nqCiguo5Q1a1ADivE=
@@ -20,9 +15,16 @@ github.com/google/go-cmp v0.5.9/go.mod h1:17dUlkBOakJ0+DkrSSNjCkIjxS6bF9zb3elmeN
 github.com/google/pprof v0.0.0-20230817174616-7a8ec2ada47b h1:h9U78+dx9a4BKdQkBBos92HalKpaGKHrp+3Uo6yTodo=
 github.com/google/pprof v0.0.0-20230817174616-7a8ec2ada47b/go.mod h1:czg5+yv1E0ZGTi6S6vVK1mke0fV+FaUhNGcd6VRS9Ik=
 github.com/google/subcommands v1.2.0/go.mod h1:ZjhPrFU+Olkh9WazFPsl27BQ4UPiG37m3yTrtFlrHVk=
+github.com/ingonyama-zk/gnark-crypto v0.0.0-20240207133741-7e078bb23cf3 h1:N5oouLARz+PIZD9FWr2U/UvNf6IMIvqch9tDlC4t3c8=
+github.com/ingonyama-zk/gnark-crypto v0.0.0-20240207133741-7e078bb23cf3/go.mod h1:Rkbv3haumAdndgfVlSAxACC2p0661ly3oAhTbQBknCY=
+github.com/ingonyama-zk/icicle v0.1.0 h1:9zbHaYv8/4g3HWRabBCpeH+64U8GJ99K1qeqE2jO6LM=
+github.com/ingonyama-zk/icicle v0.1.0/go.mod h1:kAK8/EoN7fUEmakzgZIYdWy1a2rBnpCaZLqSHwZWxEk=
+github.com/ingonyama-zk/iciclegnark v0.1.2-0.20240131141109-5f8923e3fbd5 h1:LZeh9IVCrZ8RiHZy1JuC+kKlTESUENxRrdBHwVBDFZ4=
+github.com/ingonyama-zk/iciclegnark v0.1.2-0.20240131141109-5f8923e3fbd5/go.mod h1:g17CDuMfNBiN4hhZ4aA0rGF24Abv5GBFHJqE7aLxaZQ=
+github.com/ingonyama-zk/iciclegnark v0.1.2-0.20240207160517-9558f9ad3baf h1:Z08V0nMJMwHoa6c4ASysBEZc1UU6GzRxG7XulNUGqHw=
+github.com/ingonyama-zk/iciclegnark v0.1.2-0.20240207160517-9558f9ad3baf/go.mod h1:g17CDuMfNBiN4hhZ4aA0rGF24Abv5GBFHJqE7aLxaZQ=
 github.com/kr/pretty v0.3.1 h1:flRD4NNwYAUpkphVc1HcthR4KEIFJ65n8Mw5qdRn3LE=
 github.com/kr/text v0.2.0 h1:5Nx0Ya0ZqY2ygV366QzturHI13Jq95ApcVaJBhpS+AY=
-github.com/kr/text v0.2.0/go.mod h1:eLer722TekiGuMkidMxC/pM04lWEeraHUUmBw8l2grE=
 github.com/leanovate/gopter v0.2.9 h1:fQjYxZaynp97ozCzfOyOuAGOU4aU/z37zf/tOujFk7c=
 github.com/leanovate/gopter v0.2.9/go.mod h1:U2L/78B+KVFIx2VmW6onHJQzXtFb+p5y3y2Sh+Jxxv8=
 github.com/mattn/go-colorable v0.1.12/go.mod h1:u5H1YNBxpqRaxsYJYSkiCWKzEfiAb1Gb520KVy5xxl4=
@@ -38,8 +40,7 @@ github.com/mmcloughlin/profile v0.1.1/go.mod h1:IhHD7q1ooxgwTgjxQYkACGA77oFTDdFV
 github.com/pkg/errors v0.9.1/go.mod h1:bwawxfHBFNV+L2hUp1rHADufV3IMtnDRdf1r5NINEl0=
 github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM=
 github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4=
-github.com/rogpeppe/go-internal v1.11.0 h1:cWPaGQEPrBb5/AsnsZesgZZ9yb1OQ+GOISoDNXVBh4M=
-github.com/rogpeppe/go-internal v1.11.0/go.mod h1:ddIwULY96R17DhadqLgMfk9H9tvdUzkipdSkR5nkCZA=
+github.com/rogpeppe/go-internal v1.12.0 h1:exVL4IDcn6na9z1rAb56Vxr+CgyK3nn3O+epU5NdKM8=
 github.com/rs/xid v1.5.0/go.mod h1:trrq9SKmegXys3aeAKXMUTdJsYXVwGY3RLcfgqegfbg=
 github.com/rs/zerolog v1.30.0 h1:SymVODrcRsaRaSInD9yQtKbtWqwsfoPcRff/oRXLj4c=
 github.com/rs/zerolog v1.30.0/go.mod h1:/tk+P47gFdPXq4QYjvCmT5/Gsug2nagsFWBWhAiSi1w=
@@ -47,8 +48,8 @@ github.com/stretchr/testify v1.8.4 h1:CcVxjf3Q8PM0mHUKJCdn+eZZtm5yQwehR5yeSVQQcU
 github.com/stretchr/testify v1.8.4/go.mod h1:sz/lmYIOXD/1dqDmKjjqLyZ2RngseejIcXlSw2iwfAo=
 github.com/x448/float16 v0.8.4 h1:qLwI1I70+NjRFUR3zs1JPUCgaCXSh3SW62uAKT1mSBM=
 github.com/x448/float16 v0.8.4/go.mod h1:14CWIYCyZA/cWjXOioeEpHeN/83MdbZDRQHoFcYsOfg=
-golang.org/x/crypto v0.12.0 h1:tFM/ta59kqch6LlvYnPa0yx5a83cL2nHflFhYKvv9Yk=
-golang.org/x/crypto v0.12.0/go.mod h1:NF0Gs7EO5K4qLn+Ylc+fih8BSTeIjAP05siRnAh98yw=
+golang.org/x/crypto v0.17.0 h1:r8bRNjWL3GshPW3gkd+RpvzWrZAwPS49OmTGZ/uhM4k=
+golang.org/x/crypto v0.17.0/go.mod h1:gCAAfMLgwOJRpTjQ2zCCt2OcSfYMTeZVSRtQlPC7Nq4=
 golang.org/x/exp v0.0.0-20230817173708-d852ddb80c63 h1:m64FZMko/V45gv0bNmrNYoDEq8U5YUhetc9cBWKS1TQ=
 golang.org/x/exp v0.0.0-20230817173708-d852ddb80c63/go.mod h1:0v4NqG35kSWCMzLaMeX+IQrlSnVE/bqGSyC2cz/9Le8=
 golang.org/x/sync v0.3.0 h1:ftCYgMx6zT/asHUrPw8BLLscYtGznsLAnjq5RH9P66E=
@@ -57,8 +58,8 @@ golang.org/x/sys v0.0.0-20210630005230-0f9fa26af87c/go.mod h1:oPkhp1MJrh7nUepCBc
 golang.org/x/sys v0.0.0-20210927094055-39ccf1dd6fa6/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg=
 golang.org/x/sys v0.0.0-20220811171246-fbc7d0a398ab/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg=
 golang.org/x/sys v0.6.0/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg=
-golang.org/x/sys v0.11.0 h1:eG7RXZHdqOJ1i+0lgLgCpSXAp6M3LYlAo6osgSi0xOM=
-golang.org/x/sys v0.11.0/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg=
+golang.org/x/sys v0.15.0 h1:h48lPFYpsTvQJZF4EKyI4aLHaev3CxivZmv7yZig9pc=
+golang.org/x/sys v0.15.0/go.mod h1:/VUhepiaJMQUp4+oa/7Zr1D23ma6VTLIYjOOTFZPUcA=
 gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0=
 gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c h1:Hei/4ADfdWqJk1ZMxUNpqntNwaWcugrBjAiHlqqRiVk=
 gopkg.in/yaml.v2 v2.4.0 h1:D8xgwECY7CYvx+Y2n4sBz93Jn9JRvxdiyyo8CTfuKaY=
diff --git a/internal/generator/backend/template/zkpschemes/groth16/groth16.prove.go.tmpl b/internal/generator/backend/template/zkpschemes/groth16/groth16.prove.go.tmpl
index a495057267..635f71c8cf 100644
--- a/internal/generator/backend/template/zkpschemes/groth16/groth16.prove.go.tmpl
+++ b/internal/generator/backend/template/zkpschemes/groth16/groth16.prove.go.tmpl
@@ -51,7 +51,7 @@ func Prove(r1cs *cs.R1CS, pk *ProvingKey, fullWitness witness.Witness, opts ...b
 		opt.HashToFieldFn = hash_to_field.New([]byte(constraint.CommitmentDst))
 	}
 
-	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
+	log := logger.Logger().With().Str("curve", r1cs.CurveID().String()).Str("acceleration", "none").Int("nbConstraints", r1cs.GetNbConstraints()).Str("backend", "groth16").Logger()
 
 	commitmentInfo := r1cs.CommitmentInfo.(constraint.Groth16Commitments)