1package fri
  2
  3import (
  4	"fmt"
  5	"math"
  6	"math/big"
  7	"math/bits"
  8
  9	"github.com/consensys/gnark-crypto/field/goldilocks"
 10	"github.com/consensys/gnark/frontend"
 11	gl "github.com/succinctlabs/gnark-plonky2-verifier/goldilocks"
 12	"github.com/succinctlabs/gnark-plonky2-verifier/poseidon"
 13	"github.com/succinctlabs/gnark-plonky2-verifier/types"
 14)
 15
 16type Chip struct {
 17	api               frontend.API     `gnark:"-"`
 18	gl                gl.GoldilocksApi `gnark:"-"`
 19	poseidonBN254Chip *poseidon.BN254Chip
 20	friParams         *types.FriParams `gnark:"-"`
 21}
 22
 23func NewChip(
 24	api frontend.API,
 25	friParams *types.FriParams,
 26) *Chip {
 27	poseidonBN254Chip := poseidon.NewBN254Chip(api)
 28	return &Chip{
 29		api:               api,
 30		poseidonBN254Chip: poseidonBN254Chip,
 31		friParams:         friParams,
 32		gl:                *gl.NewGoldilocksApi(api),
 33	}
 34}
 35
 36func (f *Chip) assertLeadingZeros(powWitness gl.GoldilocksVariable, friConfig types.FriConfig) {
 37	// Asserts that powWitness'es big-endian bit representation has at least `leading_zeros` leading zeros.
 38	// Note that this is assuming that the Goldilocks field is being used.  Specfically that the
 39	// field is 64 bits long
 40	maxPowWitness := uint64(math.Pow(2, float64(64-friConfig.ProofOfWorkBits))) - 1
 41	reducedPowWitness := f.gl.Reduce(powWitness)
 42	f.api.AssertIsLessOrEqual(reducedPowWitness.Limb, frontend.Variable(maxPowWitness))
 43}
 44
 45func (f *Chip) fromOpeningsAndAlpha(
 46	openings *Openings,
 47	alpha gl.QuadraticExtensionVariable,
 48) []gl.QuadraticExtensionVariable {
 49	// One reduced opening for all openings evaluated at point Zeta.
 50	// Another one for all openings evaluated at point Zeta * Omega (which is only PlonkZsNext polynomial)
 51
 52	reducedOpenings := make([]gl.QuadraticExtensionVariable, 0, 2)
 53	for _, batch := range openings.Batches {
 54		reducedOpenings = append(reducedOpenings, f.gl.ReduceWithPowers(batch.Values, alpha))
 55	}
 56
 57	return reducedOpenings
 58}
 59
 60func (f *Chip) verifyMerkleProofToCapWithCapIndex(
 61	leafData []gl.GoldilocksVariable,
 62	leafIndexBits []frontend.Variable,
 63	capIndexBits []frontend.Variable,
 64	merkleCap types.FriMerkleCap,
 65	proof *types.FriMerkleProof,
 66) {
 67	currentDigest := f.poseidonBN254Chip.HashOrNoop(leafData)
 68	for i, sibling := range proof.Siblings {
 69		bit := leafIndexBits[i]
 70		// TODO: Don't need to do two hashes by using a trick that the plonky2 verifier circuit does
 71		// https://github.com/mir-protocol/plonky2/blob/973624f12d2d12d74422b3ea051358b9eaacb050/plonky2/src/gates/poseidon.rs#L298
 72		leftHash := f.poseidonBN254Chip.TwoToOne(sibling, currentDigest)
 73		rightHash := f.poseidonBN254Chip.TwoToOne(currentDigest, sibling)
 74		currentDigest = f.api.Select(bit, leftHash, rightHash)
 75	}
 76
 77	// We assume that the cap_height is 4.  Create two levels of the Lookup2 circuit
 78	if len(capIndexBits) != 4 || len(merkleCap) != 16 {
 79		errorMsg, _ := fmt.Printf(
 80			"capIndexBits length should be 4 and the merkleCap length should be 16.  Actual values (capIndexBits: %d, merkleCap: %d)\n",
 81			len(capIndexBits),
 82			len(merkleCap),
 83		)
 84		panic(errorMsg)
 85	}
 86
 87	const NUM_LEAF_LOOKUPS = 4
 88	var leafLookups [NUM_LEAF_LOOKUPS]poseidon.BN254HashOut
 89	// First create the "leaf" lookup2 circuits
 90	// The will use the least significant bits of the capIndexBits array
 91	for i := 0; i < NUM_LEAF_LOOKUPS; i++ {
 92		leafLookups[i] = f.api.Lookup2(
 93			capIndexBits[0], capIndexBits[1],
 94			merkleCap[i*NUM_LEAF_LOOKUPS], merkleCap[i*NUM_LEAF_LOOKUPS+1], merkleCap[i*NUM_LEAF_LOOKUPS+2], merkleCap[i*NUM_LEAF_LOOKUPS+3],
 95		)
 96	}
 97
 98	// Use the most 2 significant bits of the capIndexBits array for the "root" lookup
 99	merkleCapEntry := f.api.Lookup2(capIndexBits[2], capIndexBits[3], leafLookups[0], leafLookups[1], leafLookups[2], leafLookups[3])
100	f.api.AssertIsEqual(currentDigest, merkleCapEntry)
101}
102
103func (f *Chip) verifyInitialProof(xIndexBits []frontend.Variable, proof *types.FriInitialTreeProof, initialMerkleCaps []types.FriMerkleCap, capIndexBits []frontend.Variable) {
104	if len(proof.EvalsProofs) != len(initialMerkleCaps) {
105		panic("length of eval proofs in fri proof should equal length of initial merkle caps")
106	}
107
108	for i := 0; i < len(initialMerkleCaps); i++ {
109		evals := proof.EvalsProofs[i].Elements
110		merkleProof := proof.EvalsProofs[i].MerkleProof
111		cap := initialMerkleCaps[i]
112		f.verifyMerkleProofToCapWithCapIndex(evals, xIndexBits, capIndexBits, cap, &merkleProof)
113	}
114}
115
116// / We decompose FRI query indices into bits without verifying that the decomposition given by
117// / the prover is the canonical one. In particular, if `x_index < 2^field_bits - p`, then the
118// / prover could supply the binary encoding of either `x_index` or `x_index + p`, since they are
119// / congruent mod `p`. However, this only occurs with probability
120// /     p_ambiguous = (2^field_bits - p) / p
121// / which is small for the field that we use in practice.
122// /
123// / In particular, the soundness error of one FRI query is roughly the codeword rate, which
124// / is much larger than this ambiguous-element probability given any reasonable parameters.
125// / Thus ambiguous elements contribute a negligible amount to soundness error.
126// /
127// / Here we compare the probabilities as a sanity check, to verify the claim above.
128func (f *Chip) assertNoncanonicalIndicesOK() {
129	numAmbiguousElems := uint64(math.MaxUint64) - goldilocks.Modulus().Uint64() + 1
130	queryError := f.friParams.Config.Rate()
131	pAmbiguous := float64(numAmbiguousElems) / float64(goldilocks.Modulus().Uint64())
132
133	// TODO:  Check that pAmbiguous value is the same as the one in plonky2 verifier
134	if pAmbiguous >= queryError*1e-5 {
135		panic("A non-negligible portion of field elements are in the range that permits non-canonical encodings. Need to do more analysis or enforce canonical encodings.")
136	}
137}
138
139func (f *Chip) expFromBitsConstBase(
140	base goldilocks.Element,
141	exponentBits []frontend.Variable,
142) gl.GoldilocksVariable {
143	product := gl.One()
144	for i, bit := range exponentBits {
145		// If the bit is on, we multiply product by base^pow.
146		// We can arithmetize this as:
147		//     product *= 1 + bit (base^pow - 1)
148		//     product = (base^pow - 1) product bit + product
149		pow := int64(1 << i)
150		basePow := goldilocks.NewElement(0)
151		basePow.Exp(base, big.NewInt(pow))
152		basePowVariable := gl.NewVariable(basePow.Uint64() - 1)
153		product = f.gl.Add(
154			f.gl.Mul(
155				f.gl.Mul(
156					basePowVariable,
157					product,
158				),
159				gl.NewVariable(bit),
160			),
161			product,
162		)
163	}
164	return product
165}
166
167func (f *Chip) calculateSubgroupX(
168	xIndexBits []frontend.Variable,
169	nLog uint64,
170) gl.GoldilocksVariable {
171	// Compute x from its index
172	// `subgroup_x` is `subgroup[x_index]`, i.e., the actual field element in the domain.
173	// TODO - Make these as global values
174	g := gl.NewVariable(gl.MULTIPLICATIVE_GROUP_GENERATOR.Uint64())
175	base := gl.PrimitiveRootOfUnity(nLog)
176
177	// Create a reverse list of xIndexBits
178	xIndexBitsRev := make([]frontend.Variable, 0)
179	for i := len(xIndexBits) - 1; i >= 0; i-- {
180		xIndexBitsRev = append(xIndexBitsRev, xIndexBits[i])
181	}
182
183	product := f.expFromBitsConstBase(base, xIndexBitsRev)
184
185	return f.gl.Mul(g, product)
186}
187
188func (f *Chip) friCombineInitial(
189	instance InstanceInfo,
190	proof types.FriInitialTreeProof,
191	friAlpha gl.QuadraticExtensionVariable,
192	subgroupX_QE gl.QuadraticExtensionVariable,
193	precomputedReducedEval []gl.QuadraticExtensionVariable,
194) gl.QuadraticExtensionVariable {
195	sum := gl.ZeroExtension()
196
197	if len(instance.Batches) != len(precomputedReducedEval) {
198		panic("len(openings) != len(precomputedReducedEval)")
199	}
200
201	for i := 0; i < len(instance.Batches); i++ {
202		batch := instance.Batches[i]
203		reducedOpenings := precomputedReducedEval[i]
204
205		point := batch.Point
206		evals := make([]gl.QuadraticExtensionVariable, 0)
207		for _, polynomial := range batch.Polynomials {
208			evals = append(
209				evals,
210				gl.QuadraticExtensionVariable{
211					proof.EvalsProofs[polynomial.OracleIndex].Elements[polynomial.PolynomialInfo],
212					gl.Zero(),
213				},
214			)
215		}
216
217		reducedEvals := f.gl.ReduceWithPowers(evals, friAlpha)
218		numerator := f.gl.SubExtensionNoReduce(reducedEvals, reducedOpenings)
219		denominator := f.gl.SubExtension(subgroupX_QE, point)
220		sum = f.gl.MulExtension(f.gl.ExpExtension(friAlpha, uint64(len(evals))), sum)
221		sum = f.gl.MulAddExtension(
222			numerator,
223			f.gl.InverseExtension(denominator),
224			sum,
225		)
226	}
227
228	return sum
229}
230
231func (f *Chip) finalPolyEval(finalPoly types.PolynomialCoeffs, point gl.QuadraticExtensionVariable) gl.QuadraticExtensionVariable {
232	ret := gl.ZeroExtension()
233	for i := len(finalPoly.Coeffs) - 1; i >= 0; i-- {
234		ret = f.gl.MulAddExtension(ret, point, finalPoly.Coeffs[i])
235	}
236	return ret
237}
238
239func (f *Chip) interpolate(
240	x gl.QuadraticExtensionVariable,
241	xPoints []gl.QuadraticExtensionVariable,
242	yPoints []gl.QuadraticExtensionVariable,
243	barycentricWeights []gl.QuadraticExtensionVariable,
244) gl.QuadraticExtensionVariable {
245	if len(xPoints) != len(yPoints) || len(xPoints) != len(barycentricWeights) {
246		panic("length of xPoints, yPoints, and barycentricWeights are inconsistent")
247	}
248
249	lX := gl.OneExtension()
250	for i := 0; i < len(xPoints); i++ {
251		lX = f.gl.SubMulExtension(x, xPoints[i], lX)
252	}
253
254	sum := gl.ZeroExtension()
255	for i := 0; i < len(xPoints); i++ {
256		sum = f.gl.AddExtension(
257			f.gl.MulExtension(
258				f.gl.DivExtension(
259					barycentricWeights[i],
260					f.gl.SubExtension(
261						x,
262						xPoints[i],
263					),
264				),
265				yPoints[i],
266			),
267			sum,
268		)
269	}
270
271	interpolation := f.gl.MulExtension(lX, sum)
272
273	returnField := interpolation
274	// Now check if x is already within the xPoints
275	for i := 0; i < len(xPoints); i++ {
276		returnField = f.gl.Lookup(
277			f.gl.IsZero(f.gl.SubExtension(x, xPoints[i])),
278			returnField,
279			yPoints[i],
280		)
281	}
282
283	return returnField
284}
285
286func (f *Chip) computeEvaluation(
287	x gl.GoldilocksVariable,
288	xIndexWithinCosetBits []frontend.Variable,
289	arityBits uint64,
290	evals []gl.QuadraticExtensionVariable,
291	beta gl.QuadraticExtensionVariable,
292) gl.QuadraticExtensionVariable {
293	arity := 1 << arityBits
294	if (len(evals)) != arity {
295		panic("len(evals) ! arity")
296	}
297	if arityBits > 8 {
298		panic("currently assuming that arityBits is <= 8")
299	}
300
301	g := gl.PrimitiveRootOfUnity(arityBits)
302	gInv := goldilocks.NewElement(0)
303	gInv.Exp(g, big.NewInt(int64(arity-1)))
304
305	// The evaluation vector needs to be reordered first.  Permute the evals array such that each
306	// element's new index is the bit reverse of it's original index.
307	// TODO:  Optimization - Since the size of the evals array should be constant (e.g. 2^arityBits),
308	//        we can just hard code the permutation.
309	permutedEvals := make([]gl.QuadraticExtensionVariable, len(evals))
310	for i := uint8(0); i < uint8(len(evals)); i++ {
311		newIndex := bits.Reverse8(i) >> arityBits
312		permutedEvals[newIndex] = evals[i]
313	}
314
315	// Want `g^(arity - rev_x_index_within_coset)` as in the out-of-circuit version. Compute it
316	// as `(g^-1)^rev_x_index_within_coset`.
317	revXIndexWithinCosetBits := make([]frontend.Variable, len(xIndexWithinCosetBits))
318	for i := 0; i < len(xIndexWithinCosetBits); i++ {
319		revXIndexWithinCosetBits[len(xIndexWithinCosetBits)-1-i] = xIndexWithinCosetBits[i]
320	}
321	start := f.expFromBitsConstBase(gInv, revXIndexWithinCosetBits)
322	cosetStart := f.gl.Mul(start, x)
323
324	xPoints := make([]gl.QuadraticExtensionVariable, len(evals))
325	yPoints := permutedEvals
326
327	// TODO: Make g_F a constant
328	g_F := gl.NewVariable(g.Uint64()).ToQuadraticExtension()
329	xPoints[0] = gl.QuadraticExtensionVariable{cosetStart, gl.Zero()}
330	for i := 1; i < len(evals); i++ {
331		xPoints[i] = f.gl.MulExtension(xPoints[i-1], g_F)
332	}
333
334	// TODO:  This is n^2.  Is there a way to do this better?
335	// Compute the barycentric weights
336	barycentricWeights := make([]gl.QuadraticExtensionVariable, len(xPoints))
337	for i := 0; i < len(xPoints); i++ {
338		barycentricWeights[i] = gl.OneExtension()
339		for j := 0; j < len(xPoints); j++ {
340			if i != j {
341				barycentricWeights[i] = f.gl.SubMulExtension(
342					xPoints[i],
343					xPoints[j],
344					barycentricWeights[i],
345				)
346			}
347		}
348		// Take the inverse of the barycentric weights
349		// TODO: Can provide a witness to this value
350		barycentricWeights[i] = f.gl.InverseExtension(barycentricWeights[i])
351	}
352
353	return f.interpolate(beta, xPoints, yPoints, barycentricWeights)
354}
355
356func (f *Chip) verifyQueryRound(
357	instance InstanceInfo,
358	challenges *types.FriChallenges,
359	precomputedReducedEval []gl.QuadraticExtensionVariable,
360	initialMerkleCaps []types.FriMerkleCap,
361	proof *types.FriProof,
362	xIndex gl.GoldilocksVariable,
363	n uint64,
364	nLog uint64,
365	roundProof *types.FriQueryRound,
366) {
367	f.assertNoncanonicalIndicesOK()
368	xIndex = f.gl.Reduce(xIndex)
369	xIndexBits := f.api.ToBinary(xIndex.Limb, 64)[0 : f.friParams.DegreeBits+f.friParams.Config.RateBits]
370	capIndexBits := xIndexBits[len(xIndexBits)-int(f.friParams.Config.CapHeight):]
371
372	f.verifyInitialProof(xIndexBits, &roundProof.InitialTreesProof, initialMerkleCaps, capIndexBits)
373
374	subgroupX := f.calculateSubgroupX(
375		xIndexBits,
376		nLog,
377	)
378
379	subgroupX_QE := subgroupX.ToQuadraticExtension()
380
381	oldEval := f.friCombineInitial(
382		instance,
383		roundProof.InitialTreesProof,
384		challenges.FriAlpha,
385		subgroupX_QE,
386		precomputedReducedEval,
387	)
388
389	for i, arityBits := range f.friParams.ReductionArityBits {
390		evals := roundProof.Steps[i].Evals
391
392		cosetIndexBits := xIndexBits[arityBits:]
393		xIndexWithinCosetBits := xIndexBits[:arityBits]
394
395		// Assumes that the arity bits will be 4.  That means that the range of
396		// xIndexWithCoset is [0,2^4-1].  This is based on plonky2's circuit recursive
397		// config:  https://github.com/mir-protocol/plonky2/blob/main/plonky2/src/plonk/circuit_data.rs#L63
398		// Will use a two levels tree of 4-selector gadgets.
399		if arityBits != 4 {
400			panic("assuming arity bits is 4")
401		}
402
403		const NUM_LEAF_LOOKUPS = 4
404		var leafLookups [NUM_LEAF_LOOKUPS]gl.QuadraticExtensionVariable
405		// First create the "leaf" lookup2 circuits
406		// The will use the least significant bits of the xIndexWithCosetBits array
407		for i := 0; i < NUM_LEAF_LOOKUPS; i++ {
408			leafLookups[i] = f.gl.Lookup2(
409				xIndexWithinCosetBits[0],
410				xIndexWithinCosetBits[1],
411				evals[i*NUM_LEAF_LOOKUPS],
412				evals[i*NUM_LEAF_LOOKUPS+1],
413				evals[i*NUM_LEAF_LOOKUPS+2],
414				evals[i*NUM_LEAF_LOOKUPS+3],
415			)
416		}
417
418		// Use the most 2 significant bits of the xIndexWithCosetBits array for the "root" lookup
419		newEval := f.gl.Lookup2(
420			xIndexWithinCosetBits[2],
421			xIndexWithinCosetBits[3],
422			leafLookups[0],
423			leafLookups[1],
424			leafLookups[2],
425			leafLookups[3],
426		)
427
428		f.gl.AssertIsEqual(newEval[0], oldEval[0])
429		f.gl.AssertIsEqual(newEval[1], oldEval[1])
430
431		oldEval = f.computeEvaluation(
432			subgroupX,
433			xIndexWithinCosetBits,
434			arityBits,
435			evals,
436			challenges.FriBetas[i],
437		)
438
439		// Convert evals (array of QE) to fields by taking their 0th degree coefficients
440		fieldEvals := make([]gl.GoldilocksVariable, 0, 2*len(evals))
441		for j := 0; j < len(evals); j++ {
442			fieldEvals = append(fieldEvals, evals[j][0])
443			fieldEvals = append(fieldEvals, evals[j][1])
444		}
445		f.verifyMerkleProofToCapWithCapIndex(
446			fieldEvals,
447			cosetIndexBits,
448			capIndexBits,
449			proof.CommitPhaseMerkleCaps[i],
450			&roundProof.Steps[i].MerkleProof,
451		)
452
453		// Update the point x to x^arity.
454		for j := uint64(0); j < arityBits; j++ {
455			subgroupX = f.gl.Mul(subgroupX, subgroupX)
456		}
457
458		xIndexBits = cosetIndexBits
459	}
460
461	subgroupX_QE = subgroupX.ToQuadraticExtension()
462	finalPolyEval := f.finalPolyEval(proof.FinalPoly, subgroupX_QE)
463
464	f.gl.AssertIsEqual(oldEval[0], finalPolyEval[0])
465	f.gl.AssertIsEqual(oldEval[1], finalPolyEval[1])
466}
467
468func (f *Chip) VerifyFriProof(
469	instance InstanceInfo,
470	openings Openings,
471	friChallenges *types.FriChallenges,
472	initialMerkleCaps []types.FriMerkleCap,
473	friProof *types.FriProof,
474) {
475	// TODO:  Check fri config
476	/* if let Some(max_arity_bits) = params.max_arity_bits() {
477		self.check_recursion_config::<C>(max_arity_bits);
478	}
479
480	debug_assert_eq!(
481		params.final_poly_len(),
482		proof.final_poly.len(),
483		"Final polynomial has wrong degree."
484	); */
485
486	// Check POW
487
488	f.assertLeadingZeros(friChallenges.FriPowResponse, f.friParams.Config)
489
490	precomputedReducedEvals := f.fromOpeningsAndAlpha(&openings, friChallenges.FriAlpha)
491
492	// Size of the LDE domain.
493	nLog := f.friParams.DegreeBits + f.friParams.Config.RateBits
494	n := uint64(math.Pow(2, float64(nLog)))
495
496	if len(friChallenges.FriQueryIndices) != len(friProof.QueryRoundProofs) {
497		panic(fmt.Sprintf(
498			"Number of query indices (%d) should equal number of query round proofs (%d)",
499			len(friChallenges.FriQueryIndices),
500			len(friProof.QueryRoundProofs),
501		))
502	}
503
504	for idx, xIndex := range friChallenges.FriQueryIndices {
505		roundProof := friProof.QueryRoundProofs[idx]
506
507		f.verifyQueryRound(
508			instance,
509			friChallenges,
510			precomputedReducedEvals,
511			initialMerkleCaps,
512			friProof,
513			xIndex,
514			n,
515			nLog,
516			&roundProof,
517		)
518	}
519}