Polynomial refactoring

[mratsim]

The Polynomial primitives need to be refactored.

• EIP-7594 (EIP-7594: Data Availability Sampling #341) will require a lot of new polynomial primitives:
  https://github.com/ethereum/consensus-specs/blob/fe8db03/specs/_features/eip7594/polynomial-commitments-sampling.md

  Polynomials in coefficient form
  polynomial_eval_to_coeff
  add_polynomialcoeff
  neg_polynomialcoeff
  multiply_polynomialcoeff
  divide_polynomialcoeff
  shift_polynomialcoeff
  interpolate_polynomialcoeff
  vanishing_polynomialcoeff
  evaluate_polynomialcoeff

• As will multilinear polynomial commitment schemes Multilinear Polynomial Commitments #299

The first step is unifying the Lagrange / Barycentric form polynomial division that are implemented twice:

• for KZG:

  constantine/constantine/math/polynomials/polynomials.nim

  Lines 157 to 221 in 4e964ff

  	func differenceQuotientEvalInDomain*[N: static int, Field](
  	r: var PolynomialEval[N, Field],
  	poly: PolynomialEval[N, Field],
  	zIndex: uint32,
  	invRootsMinusZ: array[N, Field],
  	domain: PolyDomainEval[N, Field],
  	isBitReversedDomain: static bool) =
  	## Compute r(x) = (p(x) - p(z)) / (x - z)
  	##
  	## for z = ωⁱ a power of a root of unity
  	##
  	## Input:
  	## - poly: p(x) a polynomial in evaluation form as an array of p(ωⁱ)
  	## - rootsOfUnity: ωⁱ
  	## - invRootsMinusZ: 1/(ωⁱ-z)
  	## - zIndex: the index of the root of unity power that matches z = ωⁱᵈˣ
  	static:
  	# For powers of 2: x mod N == x and (N-1)
  	doAssert N.isPowerOf2_vartime()
  	r.evals[zIndex].setZero()
  	for i in 0'u32 ..< N:
  	if i == zIndex:
  	# https://dankradfeist.de/ethereum/2021/06/18/pcs-multiproofs.html
  	# section "Dividing when one of the points is zero".
  	continue
  	# qᵢ = (p(ωⁱ) - p(z))/(ωⁱ-z)
  	var qi {.noinit.}: Field
  	qi.diff(poly.evals[i], poly.evals[zIndex])
  	r.evals[i].prod(qi, invRootsMinusZ[i])
  	# q'ᵢ = -qᵢ * ωⁱ/z
  	# q'idx = ∑ q'ᵢ
  	var ri {.noinit.}: Field
  	ri.neg(r.evals[i]) # -qᵢ
  	when isBitReversedDomain:
  	const logN = log2_vartime(uint32 N)
  	let invZidx = N - reverseBits(uint32 zIndex, logN)
  	let canonI = reverseBits(i, logN)
  	let idx = reverseBits((canonI + invZidx) and (N-1), logN)
  	ri *= domain.rootsOfUnity[idx] # -qᵢ * ωⁱ/z (explanation at the bottom)
  	else:
  	ri *= domain.rootsOfUnity[(i+N-zIndex) and (N-1)] # -qᵢ * ωⁱ/z (explanation at the bottom)
  	r.evals[zIndex] += ri # r[zIndex] = ∑ -qᵢ * ωⁱ/z
  	# * 1/z computation detail
  	# from ωⁿ = 1 and z = ωⁱᵈˣ
  	# hence ωⁿ⁻ⁱᵈˣ = 1/z
  	# However our z may be in bit-reversal permuted
  	#
  	# * We want ωⁱ/z which translate to ωⁱ*ωⁿ⁻ⁱᵈˣ hence ωⁱ⁺ⁿ⁻ⁱᵈˣ
  	# with the roots of unity being a cyclic group of order N so we compute i+N-zIndex (mod N)
  	#
  	# However some protocols use bit-reversal permutation (brp) to store the ωⁱ
  	# Hence retrieving the data requires roots[brp((brp(i)-n-brp(idx)) mod n)] for those (note: n = brp(n))
  	#
  	# For Ethereum:
  	# A 254~255-bit multiplication takes 11ns / 38 cycles (Fr[BLS12-381]),
  	# A brp with n = 2¹² = 4096 (for EIP4844) takes about 6ns
  	# We could also cache either ωⁿ⁻ⁱ or a map i' = brp(n - brp(i))
  	# in non-brp order but cache misses are expensive
  	# and brp can benefits from instruction-level parallelism

• for IPA:

  constantine/constantine/eth_verkle_ipa/barycentric_form.nim

  Lines 149 to 178 in 4e964ff

  	func divisionOnDomain*(res: var array[VerkleDomain,Fr[Banderwagon]], precomp: PrecomputedWeights, index: var int, f: openArray[Fr[Banderwagon]]) =
  	## Computes f(x) - f(x_i) / x - x_i using the barycentric weights, where x_i is an element in the
  	var is_negative = true
  	var y = f[index]
  	for i in 0 ..< VerkleDomain:
  	if i != index:
  	var denominator = i - index
  	var absDenominator {.noInit.}: int
  	absDenominator.absIntChecker(denominator)
  	if absDenominator > 0:
  	is_negative = false
  	var denominatorInv {.noInit.}: Fr[Banderwagon]
  	denominatorInv.getInvertedElement(precomp, absDenominator, is_negative)
  	res[i].diff(f[i], y)
  	res[i] *= denominatorInv
  	var weight_ratios {.noInit.}: Fr[Banderwagon]
  	var dummy {.noInit.} : int
  	dummy = i
  	weight_ratios.getWeightRatios(precomp, index, dummy)
  	var tmp {.noInit.}: Fr[Banderwagon]
  	tmp.prod(weight_ratios, res[i])
  	res[index] -= tmp