Compute KZG proof directly from evaluation form

Previously we were using IFFT to go from evaluation form to coefficient form before computing the KZG proof in the
traditional way.

We can instead do the KZG proof division directly in evaluation form by using point-by-point division:
   https://dankradfeist.de/ethereum/2021/06/18/pcs-multiproofs.html#dividing-when-one-of-the-points-is-zero

Two fine points:
- We need to check that the denominator of the division is not zero. It's very unlikely because `z` comes out of
Fiat-Shamir, but we still add an assert to make sure.

- We need to shift our main polynomial by `p(z)` as in the traditional KZG proof formula. We can avoid this
in coefficient form because with long division we can extract the quotient polynomial and ignore the remainder.
However that's not possible when doing point-by-point division in evaluation form and hence we need to ensure
that there is no remainder by shifting by `p(z)`.

specs/eip4844/polynomial-commitments.md

@@ -150,18 +150,24 @@ def verify_kzg_proof(polynomial_kzg: KZGCommitment,
 #### `compute_kzg_proof`
 
 ```python
-def compute_kzg_proof(polynomial: Sequence[BLSFieldElement], x: BLSFieldElement) -> KZGProof:
+def compute_kzg_proof(polynomial: Sequence[BLSFieldElement], z: BLSFieldElement) -> KZGProof:
+    """Compute KZG proof at point `z` with `polynomial` being in evaluation form"""
+
     # To avoid SSZ overflow/underflow, convert element into int
     polynomial = [int(i) for i in polynomial]
+    z = int(z)
+
+    # Shift our polynomial first (in evaluation form we can't handle the division remainder)
+    y = evaluate_polynomial_in_evaluation_form(polynomial, z)
+    polynomial_shifted = [(p - int(y)) % BLS_MODULUS for p in polynomial]
 
-    # Convert `polynomial` to coefficient form
-    assert pow(ROOTS_OF_UNITY[1], len(polynomial), BLS_MODULUS) == 1
-    fft_output = fft(polynomial, ROOTS_OF_UNITY)
-    inv_length = pow(len(polynomial), BLS_MODULUS - 2, BLS_MODULUS)
-    polynomial_in_coefficient_form = [fft_output[-i] * inv_length % BLS_MODULUS for i in range(len(fft_output))]
+    # Make sure we won't divide by zero during division
+    assert z not in ROOTS_OF_UNITY
+    denominator_poly = [(x - z) % BLS_MODULUS for x in ROOTS_OF_UNITY]
 
-    quotient_polynomial = div_polys(polynomial_in_coefficient_form, [-int(x), 1])
-    return KZGProof(lincomb(KZG_SETUP_G1[:len(quotient_polynomial)], quotient_polynomial))
+    # Calculate quotient polynomial by doing point-by-point division
+    quotient_polynomial = [div(a, b) for a, b in zip(polynomial_shifted, denominator_poly)]
+    return KZGProof(lincomb(KZG_SETUP_LAGRANGE[:len(quotient_polynomial)], quotient_polynomial))
 ```
 
 ### Polynomials
@@ -187,43 +193,3 @@ def evaluate_polynomial_in_evaluation_form(polynomial: Sequence[BLSFieldElement]
     return result
 ```
 
-#### `fft`
-
-```python
-def fft(vals, domain):
-    """
-    FFT for ``BLSFieldElement``.
-    """
-    if len(vals) == 1:
-        return vals
-    L = fft(vals[::2], domain[::2])
-    R = fft(vals[1::2], domain[::2])
-    result = [0] * len(vals)
-    for i, (x, y) in enumerate(zip(L, R)):
-        y_times_root = y * domain[i] % BLS_MODULUS
-        result[i] = x + y_times_root % BLS_MODULUS
-        result[i + len(L)] = x + (BLS_MODULUS - y_times_root) % BLS_MODULUS
-    return result
-```
-
-#### `div_polys`
-
-```python
-def div_polys(a: Sequence[int], b: Sequence[int]) -> Sequence[BLSFieldElement]:
-    """
-    Long polynomial division for two polynomials in coefficient form.
-    """
-    a = a.copy()  # avoid side effects
-    result = []
-    a_position = len(a) - 1
-    b_position = len(b) - 1
-    diff = a_position - b_position
-    while diff >= 0:
-        quotient = div(a[a_position], b[b_position])
-        result.insert(0, quotient)
-        for i in range(b_position, -1, -1):
-            a[diff + i] -= b[i] * quotient
-        a_position -= 1
-        diff -= 1
-    return [x % BLS_MODULUS for x in result]
-```