Implement Composition Polynomial

[varunthakore]

This PR implements the composition polynomial library.
https://github.com/starkware-libs/stone-prover/tree/main/src/starkware/composition_polynomial

[rot256]

Is this different from self.constraints.len()?

[varunthakore]

its the same, removed n_constraints.

[rot256]

Consider introducing a struct for the argument: that way they are named and it is easy to see what a ConstraintFunction takes e.g.

struct ConstraintEval {
  neighbors: &[F]
  periodic_columns: &[F]
  random_coefficients: $[F]
  point: &[F]
  gen_powers: &[F]
  precomp_domains: &[F]
}

[rot256]

I thought the random_coefficients where used to combine each constraint? Does each constraint need to take the random coefficients or can we compute the combination here?

[rot256]

Could this just be a method on CompositionPolynomial<F>?

e.g. it takes random_coefficients which is really what defines a CompositionPolynomial from an AIR. So if we move it there we would not need to take e.g. random_coefficients

[rot256]

Could this take ownership of the Air?

i.e.

fn into_composition_polynomial(
   self,
   trace_generator,
   random_coefficients
) -> CompositionPolynomial<F>

That way we avoid the Box and Clone?

[rot256]

nit: assert_eq

[rot256]

assert!(values.is_power_of_two())

[rot256]

What's the difference between start_point (as used in get_coset) and offset? Both should be the first evaluation point of the first entry of the first copy of the periodic table?

[rot256]

All these are restricted to powers of two?

[rot256]

Nvm: the LDE manager is used to interpolate exactly this polynomial.

Why do we need the lde_manager?

We have values and assuming the length of the table divides the group order the polynomial for the column has a succinct description (and can be evaluated efficiently). So can't we just evaluate it directly?

For instance, the polynomial f(X) which evaluates to:

f(1)   = y0
f(w)   = y1
f(w^2) = y0
f(w^3) = y1
...

i.e. alternating. For a domain of length $2^k$, is simply:

$$ f(X) = a + b * X^{2^k / 2} $$

Where $a, b$:

$$ f(1) = a + b = y0 $$

$$ f(w) = a + b * w^{2^k / 2} = a - b = y1 $$

So $a = (y0 + y1) / 2$ and $b = (y0 - y1) / 2$. A similar deal applies when you have a "starting_point", except now the polynomial is composed with * starting_point.

Here is a Sage script if you want to play around this:

F = GF(2^251 + 17 * 2^192 + 1)
k = 192

g = F.multiplicative_generator()
g = g^(g.multiplicative_order() / 2^k)

order = 2^8
w = g^(g.multiplicative_order() / order)

print(f"order: {order}")
print(g.multiplicative_order())

table = [
    1, 2, 1, 2
]

ys = table * (order // len(table))
xs = [w^i for i in range(order)]

R.<X> = PolynomialRing(F)

f = R.lagrange_polynomial(zip(xs, ys))

print(f)

Otherwise, let's hop on a call.