func evaluateHDomain(poly []fr.Element, domainH *fft.Domain) []fr.Element {
	res := make([]fr.Element, domainH.Cardinality)
	utils.Parallelize(len(poly), func(start, end int) {
		for i := start; i < end; i++ {
			res[i].Mul(&poly[i], &domainH.CosetTable[0][i])
		}
	}, runtime.NumCPU()/2)
	domainH.FFT(res, fft.DIF, 0)
	return res
}
evaluateHDomain evaluates poly (canonical form) of degree m<n where n=domainH.Cardinality
on the odd coset of (Z/2nZ)/(Z/nZ).

Puts the result in res of size n.
Warning: result is in bit reversed order, we do a bit reverse operation only once in computeH
package plonk
import (
	"crypto/sha256"
	"math/big"
	"math/bits"
	"runtime"
	"sync"
	"github.com/consensys/gnark-crypto/ecc/bn254/fr"
	curve "github.com/consensys/gnark-crypto/ecc/bn254"
	"github.com/consensys/gnark-crypto/ecc/bn254/fr/polynomial"
	"github.com/consensys/gnark-crypto/ecc/bn254/fr/kzg"
	"github.com/consensys/gnark-crypto/ecc/bn254/fr/fft"
	bn254witness "github.com/consensys/gnark/internal/backend/bn254/witness"
	"github.com/consensys/gnark/internal/backend/bn254/cs"
	"github.com/consensys/gnark-crypto/fiat-shamir"
	"github.com/consensys/gnark/backend"
	"github.com/consensys/gnark/internal/utils"
)
type Proof struct {
	LRO [3]kzg.Digest
	Z kzg.Digest
	H [3]kzg.Digest
	BatchedProof kzg.BatchOpeningProof
	ZShiftedOpening kzg.OpeningProof
}
func Prove(spr *cs.SparseR1CS, pk *ProvingKey, fullWitness bn254witness.Witness, opt backend.ProverOption) (*Proof, error) {
	hFunc := sha256.New()
	fs := fiatshamir.NewTranscript(hFunc, "gamma", "alpha", "zeta")
	proof := &Proof{}
	var solution []fr.Element
	var err error
	if solution, err = spr.Solve(fullWitness, opt); err != nil {
		if !opt.Force {
			return nil, err
		} else {
			var r fr.Element
			_, _ = r.SetRandom()
			for i := spr.NbPublicVariables + spr.NbSecretVariables; i < len(solution); i++ {
				solution[i] = r
				r.Double(&r)
			}
		}
	}
	ll, lr, lo := computeLRO(spr, pk, solution)
	bcl, bcr, bco, err := computeBlindedLRO(ll, lr, lo, &pk.DomainNum)
	if err != nil {
		return nil, err
	}
	if err := commitToLRO(bcl, bcr, bco, proof, pk.Vk.KZGSRS); err != nil {
		return nil, err
	}
	gamma, err := deriveRandomness(&fs, "gamma", &proof.LRO[0], &proof.LRO[1], &proof.LRO[2])
	if err != nil {
		return nil, err
	}
	var bz polynomial.Polynomial
	chZ := make(chan error, 1)
	var alpha fr.Element
	go func() {
		var err error
		bz, err = computeBlindedZ(ll, lr, lo, pk, gamma)
		if err != nil {
			chZ <- err
			close(chZ)
			return
		}
		if proof.Z, err = kzg.Commit(bz, pk.Vk.KZGSRS, runtime.NumCPU()*2); err != nil {
			chZ <- err
			close(chZ)
			return
		}
		alpha, err = deriveRandomness(&fs, "alpha", &proof.Z)
		chZ <- err
		close(chZ)
	}()
	var evalBL, evalBR, evalBO, evalBZ polynomial.Polynomial
	chEvalBL := make(chan struct{}, 1)
	chEvalBR := make(chan struct{}, 1)
	chEvalBO := make(chan struct{}, 1)
	go func() {
		evalBL = evaluateHDomain(bcl, &pk.DomainH)
		close(chEvalBL)
	}()
	go func() {
		evalBR = evaluateHDomain(bcr, &pk.DomainH)
		close(chEvalBR)
	}()
	go func() {
		evalBO = evaluateHDomain(bco, &pk.DomainH)
		close(chEvalBO)
	}()
	var constraintsInd, constraintsOrdering polynomial.Polynomial
	chConstraintInd := make(chan struct{}, 1)
	go func() {
		qk := make(polynomial.Polynomial, pk.DomainNum.Cardinality)
		copy(qk, fullWitness[:spr.NbPublicVariables])
		copy(qk[spr.NbPublicVariables:], pk.LQk[spr.NbPublicVariables:])
		pk.DomainNum.FFTInverse(qk, fft.DIF, 0)
		fft.BitReverse(qk)
		<-chEvalBL
		<-chEvalBR
		<-chEvalBO
		constraintsInd = evalConstraints(pk, evalBL, evalBR, evalBO, qk)
		close(chConstraintInd)
	}()
	chConstraintOrdering := make(chan error, 1)
	go func() {
		if err := <-chZ; err != nil {
			chConstraintOrdering <- err
			return
		}
		evalBZ = evaluateHDomain(bz, &pk.DomainH)
		<-chEvalBL
		<-chEvalBR
		<-chEvalBO
		constraintsOrdering = evalConstraintOrdering(pk, evalBZ, evalBL, evalBR, evalBO, gamma)
		chConstraintOrdering <- nil
		close(chConstraintOrdering)
	}()
	if err := <-chConstraintOrdering; err != nil {
		return nil, err
	}
	<-chConstraintInd
	h1, h2, h3 := computeH(pk, constraintsInd, constraintsOrdering, evalBZ, alpha)
	if err := commitToH(h1, h2, h3, proof, pk.Vk.KZGSRS); err != nil {
		return nil, err
	}
	zeta, err := deriveRandomness(&fs, "zeta", &proof.H[0], &proof.H[1], &proof.H[2])
	if err != nil {
		return nil, err
	}
	var blzeta, brzeta, bozeta fr.Element
	var wgZetaEvals sync.WaitGroup
	wgZetaEvals.Add(3)
	go func() {
		blzeta = bcl.Eval(&zeta)
		wgZetaEvals.Done()
	}()
	go func() {
		brzeta = bcr.Eval(&zeta)
		wgZetaEvals.Done()
	}()
	go func() {
		bozeta = bco.Eval(&zeta)
		wgZetaEvals.Done()
	}()
	var zetaShifted fr.Element
	zetaShifted.Mul(&zeta, &pk.Vk.Generator)
	proof.ZShiftedOpening, err = kzg.Open(
		bz,
		&zetaShifted,
		&pk.DomainH,
		pk.Vk.KZGSRS,
	)
	if err != nil {
		return nil, err
	}
	bzuzeta := proof.ZShiftedOpening.ClaimedValue
	var (
		linearizedPolynomial       polynomial.Polynomial
		linearizedPolynomialDigest curve.G1Affine
		errLPoly                   error
	)
	chLpoly := make(chan struct{}, 1)
	go func() {
		wgZetaEvals.Wait()
		linearizedPolynomial = computeLinearizedPolynomial(
			blzeta,
			brzeta,
			bozeta,
			alpha,
			gamma,
			zeta,
			bzuzeta,
			bz,
			pk,
		)
		linearizedPolynomialDigest, errLPoly = kzg.Commit(linearizedPolynomial, pk.Vk.KZGSRS)
		close(chLpoly)
	}()
	var bZetaPowerm, bSize big.Int
	bSize.SetUint64(pk.DomainNum.Cardinality + 2)  
	var zetaPowerm fr.Element
	zetaPowerm.Exp(zeta, &bSize)
	zetaPowerm.ToBigIntRegular(&bZetaPowerm)
	foldedHDigest := proof.H[2]
	foldedHDigest.ScalarMultiplication(&foldedHDigest, &bZetaPowerm)
	foldedHDigest.Add(&foldedHDigest, &proof.H[1])                    
	foldedHDigest.ScalarMultiplication(&foldedHDigest, &bZetaPowerm)  
	foldedHDigest.Add(&foldedHDigest, &proof.H[0])                    
	foldedH := h3
	utils.Parallelize(len(foldedH), func(start, end int) {
		for i := start; i < end; i++ {
			foldedH[i].Mul(&foldedH[i], &zetaPowerm)  
			foldedH[i].Add(&foldedH[i], &h2[i])       
			foldedH[i].Mul(&foldedH[i], &zetaPowerm)  
			foldedH[i].Add(&foldedH[i], &h1[i])       
		}
	})
	<-chLpoly
	if errLPoly != nil {
		return nil, errLPoly
	}
	proof.BatchedProof, err = kzg.BatchOpenSinglePoint(
		[]polynomial.Polynomial{
			foldedH,
			linearizedPolynomial,
			bcl,
			bcr,
			bco,
			pk.CS1,
			pk.CS2,
		},
		[]kzg.Digest{
			foldedHDigest,
			linearizedPolynomialDigest,
			proof.LRO[0],
			proof.LRO[1],
			proof.LRO[2],
			pk.Vk.S[0],
			pk.Vk.S[1],
		},
		&zeta,
		hFunc,
		&pk.DomainH,
		pk.Vk.KZGSRS,
	)
	if err != nil {
		return nil, err
	}
	return proof, nil
}
func commitToLRO(bcl, bcr, bco polynomial.Polynomial, proof *Proof, srs *kzg.SRS) error {
	n := runtime.NumCPU() / 2
	var err0, err1, err2 error
	chCommit0 := make(chan struct{}, 1)
	chCommit1 := make(chan struct{}, 1)
	go func() {
		proof.LRO[0], err0 = kzg.Commit(bcl, srs, n)
		close(chCommit0)
	}()
	go func() {
		proof.LRO[1], err1 = kzg.Commit(bcr, srs, n)
		close(chCommit1)
	}()
	if proof.LRO[2], err2 = kzg.Commit(bco, srs, n); err2 != nil {
		return err2
	}
	<-chCommit0
	<-chCommit1
	if err0 != nil {
		return err0
	}
	return err1
}
func commitToH(h1, h2, h3 polynomial.Polynomial, proof *Proof, srs *kzg.SRS) error {
	n := runtime.NumCPU() / 2
	var err0, err1, err2 error
	chCommit0 := make(chan struct{}, 1)
	chCommit1 := make(chan struct{}, 1)
	go func() {
		proof.H[0], err0 = kzg.Commit(h1, srs, n)
		close(chCommit0)
	}()
	go func() {
		proof.H[1], err1 = kzg.Commit(h2, srs, n)
		close(chCommit1)
	}()
	if proof.H[2], err2 = kzg.Commit(h3, srs, n); err2 != nil {
		return err2
	}
	<-chCommit0
	<-chCommit1
	if err0 != nil {
		return err0
	}
	return err1
}
func computeBlindedLRO(ll, lr, lo polynomial.Polynomial, domain *fft.Domain) (bcl, bcr, bco polynomial.Polynomial, err error) {
	cl := make(polynomial.Polynomial, domain.Cardinality, domain.Cardinality+2)
	cr := make(polynomial.Polynomial, domain.Cardinality, domain.Cardinality+2)
	co := make(polynomial.Polynomial, domain.Cardinality, domain.Cardinality+2)
	chDone := make(chan error, 2)
	go func() {
		var err error
		copy(cl, ll)
		domain.FFTInverse(cl, fft.DIF, 0)
		fft.BitReverse(cl)
		bcl, err = blindPoly(cl, domain.Cardinality, 1)
		chDone <- err
	}()
	go func() {
		var err error
		copy(cr, lr)
		domain.FFTInverse(cr, fft.DIF, 0)
		fft.BitReverse(cr)
		bcr, err = blindPoly(cr, domain.Cardinality, 1)
		chDone <- err
	}()
	copy(co, lo)
	domain.FFTInverse(co, fft.DIF, 0)
	fft.BitReverse(co)
	if bco, err = blindPoly(co, domain.Cardinality, 1); err != nil {
		return
	}
	err = <-chDone
	if err != nil {
		return
	}
	err = <-chDone
	return
}
func blindPoly(cp polynomial.Polynomial, rou, bo uint64) (polynomial.Polynomial, error) {
	totalDegree := rou + bo
	res := cp[:totalDegree+1]
	blindingPoly := make(polynomial.Polynomial, bo+1)
	for i := uint64(0); i < bo+1; i++ {
		if _, err := blindingPoly[i].SetRandom(); err != nil {
			return nil, err
		}
	}
	for i := uint64(0); i < bo+1; i++ {
		res[i].Sub(&res[i], &blindingPoly[i])
		res[rou+i].Add(&res[rou+i], &blindingPoly[i])
	}
	return res, nil
}
func computeLRO(spr *cs.SparseR1CS, pk *ProvingKey, solution []fr.Element) (polynomial.Polynomial, polynomial.Polynomial, polynomial.Polynomial) {
	s := int(pk.DomainNum.Cardinality)
	var l, r, o polynomial.Polynomial
	l = make([]fr.Element, s)
	r = make([]fr.Element, s)
	o = make([]fr.Element, s)
	s0 := solution[0]
	for i := 0; i < spr.NbPublicVariables; i++ {  
		l[i] = solution[i]
		r[i] = s0
		o[i] = s0
	}
	offset := spr.NbPublicVariables
	for i := 0; i < len(spr.Constraints); i++ {  
		l[offset+i] = solution[spr.Constraints[i].L.VariableID()]
		r[offset+i] = solution[spr.Constraints[i].R.VariableID()]
		o[offset+i] = solution[spr.Constraints[i].O.VariableID()]
	}
	offset += len(spr.Constraints)
	for i := 0; i < s-offset; i++ {  
		l[offset+i] = s0
		r[offset+i] = s0
		o[offset+i] = s0
	}
	return l, r, o
}
func computeBlindedZ(l, r, o polynomial.Polynomial, pk *ProvingKey, gamma fr.Element) (polynomial.Polynomial, error) {
	z := make(polynomial.Polynomial, pk.DomainNum.Cardinality, pk.DomainNum.Cardinality+3)
	nbElmts := int(pk.DomainNum.Cardinality)
	gInv := make(polynomial.Polynomial, pk.DomainNum.Cardinality)
	z[0].SetOne()
	gInv[0].SetOne()
	utils.Parallelize(nbElmts-1, func(start, end int) {
		var f [3]fr.Element
		var g [3]fr.Element
		var u [3]fr.Element
		u[0].Exp(pk.DomainNum.Generator, new(big.Int).SetInt64(int64(start)))
		u[1].Mul(&u[0], &pk.Vk.Shifter[0])
		u[2].Mul(&u[0], &pk.Vk.Shifter[1])
		for i := start; i < end; i++ {
			f[0].Add(&l[i], &u[0]).Add(&f[0], &gamma)  
			f[1].Add(&r[i], &u[1]).Add(&f[1], &gamma)  
			f[2].Add(&o[i], &u[2]).Add(&f[2], &gamma)  
			g[0].Add(&l[i], &pk.LS1[i]).Add(&g[0], &gamma)  
			g[1].Add(&r[i], &pk.LS2[i]).Add(&g[1], &gamma)  
			g[2].Add(&o[i], &pk.LS3[i]).Add(&g[2], &gamma)  
			f[0].Mul(&f[0], &f[1]).Mul(&f[0], &f[2])  
			g[0].Mul(&g[0], &g[1]).Mul(&g[0], &g[2])  
			gInv[i+1] = g[0]
			z[i+1] = f[0]
			u[0].Mul(&u[0], &pk.DomainNum.Generator)  
			u[1].Mul(&u[1], &pk.DomainNum.Generator)  
			u[2].Mul(&u[2], &pk.DomainNum.Generator)  
		}
	})
	gInv = fr.BatchInvert(gInv)
	for i := 1; i < nbElmts; i++ {
		z[i].Mul(&z[i], &z[i-1]).
			Mul(&z[i], &gInv[i])
	}
	pk.DomainNum.FFTInverse(z, fft.DIF, 0)
	fft.BitReverse(z)
	return blindPoly(z, pk.DomainNum.Cardinality, 2)
}
func evalConstraints(pk *ProvingKey, evalL, evalR, evalO, qk []fr.Element) []fr.Element {
	var evalQl, evalQr, evalQm, evalQo, evalQk polynomial.Polynomial
	var wg sync.WaitGroup
	wg.Add(4)
	go func() {
		evalQl = evaluateHDomain(pk.Ql, &pk.DomainH)
		wg.Done()
	}()
	go func() {
		evalQr = evaluateHDomain(pk.Qr, &pk.DomainH)
		wg.Done()
	}()
	go func() {
		evalQm = evaluateHDomain(pk.Qm, &pk.DomainH)
		wg.Done()
	}()
	go func() {
		evalQo = evaluateHDomain(pk.Qo, &pk.DomainH)
		wg.Done()
	}()
	evalQk = evaluateHDomain(qk, &pk.DomainH)
	wg.Wait()
	utils.Parallelize(len(evalQk), func(start, end int) {
		var t0, t1 fr.Element
		for i := start; i < end; i++ {
			t1.Mul(&evalQm[i], &evalR[i])  
			t1.Add(&t1, &evalQl[i])        
			t1.Mul(&t1, &evalL[i])         
			t0.Mul(&evalQr[i], &evalR[i])
			t0.Add(&t0, &t1)  
			t1.Mul(&evalQo[i], &evalO[i])
			t0.Add(&t0, &t1)                
			evalQk[i].Add(&t0, &evalQk[i])  
		}
	})
	return evalQk
}
func evalIDCosets(pk *ProvingKey) (id polynomial.Polynomial) {
	id = make([]fr.Element, pk.DomainH.Cardinality)
	utils.Parallelize(int(pk.DomainH.Cardinality), func(start, end int) {
		var acc fr.Element
		acc.Exp(pk.DomainH.Generator, new(big.Int).SetInt64(int64(start)))
		for i := start; i < end; i++ {
			id[i].Mul(&acc, &pk.DomainH.FinerGenerator)
			acc.Mul(&acc, &pk.DomainH.Generator)
		}
	})
	return id
}
func evalConstraintOrdering(pk *ProvingKey, evalZ, evalL, evalR, evalO polynomial.Polynomial, gamma fr.Element) polynomial.Polynomial {
	evalID := evalIDCosets(pk)
	var wg sync.WaitGroup
	wg.Add(2)
	var evalS1, evalS2, evalS3 polynomial.Polynomial
	go func() {
		evalS1 = evaluateHDomain(pk.CS1, &pk.DomainH)
		wg.Done()
	}()
	go func() {
		evalS2 = evaluateHDomain(pk.CS2, &pk.DomainH)
		wg.Done()
	}()
	evalS3 = evaluateHDomain(pk.CS3, &pk.DomainH)
	wg.Wait()
	res := evalS1  
	s := uint64(len(evalZ))
	nn := uint64(64 - bits.TrailingZeros64(uint64(s)))
	toShift := pk.DomainH.Cardinality / pk.DomainNum.Cardinality
	utils.Parallelize(int(pk.DomainH.Cardinality), func(start, end int) {
		var f [3]fr.Element
		var g [3]fr.Element
		var eID fr.Element
		for i := start; i < end; i++ {
			irev := bits.Reverse64(uint64(i)) >> nn
			eID = evalID[irev]
			shiftedZ := bits.Reverse64(uint64((irev+toShift)%s)) >> nn
			f[0].Add(&eID, &evalL[i]).Add(&f[0], &gamma)  
			f[1].Mul(&eID, &pk.Vk.Shifter[0])
			f[2].Mul(&eID, &pk.Vk.Shifter[1])
			f[1].Add(&f[1], &evalR[i]).Add(&f[1], &gamma)  
			f[2].Add(&f[2], &evalO[i]).Add(&f[2], &gamma)  
			g[0].Add(&evalL[i], &evalS1[i]).Add(&g[0], &gamma)  
			g[1].Add(&evalR[i], &evalS2[i]).Add(&g[1], &gamma)  
			g[2].Add(&evalO[i], &evalS3[i]).Add(&g[2], &gamma)  
			f[0].Mul(&f[0], &f[1]).
				Mul(&f[0], &f[2]).
				Mul(&f[0], &evalZ[i])  
			g[0].Mul(&g[0], &g[1]).
				Mul(&g[0], &g[2]).
				Mul(&g[0], &evalZ[shiftedZ])  
			res[i].Sub(&g[0], &f[0])
		}
	})
	return res
}