Implement zero polynomial - passes tests, but could be tidier

src/Makefile

@@ -1,7 +1,7 @@
 TESTS = bls12_381_test das_extension_test c_kzg_util_test fft_common_test fft_fr_test fft_g1_test \
-	fk20_proofs_test kzg_proofs_test poly_test utility_test
+	fk20_proofs_test kzg_proofs_test poly_test utility_test zero_poly_test
 BENCH = fft_fr_bench fft_g1_bench fft_fr_in_place_0_bench fft_fr_in_place_1_bench
-LIB_SRC = bls12_381.c c_kzg_util.c das_extension.c fft_common.c fft_fr.c fft_g1.c fk20_proofs.c kzg_proofs.c poly.c utility.c
+LIB_SRC = bls12_381.c c_kzg_util.c das_extension.c fft_common.c fft_fr.c fft_g1.c fk20_proofs.c kzg_proofs.c poly.c utility.c zero_poly.c
 LIB_OBJ = $(LIB_SRC:.c=.o)
 
 CFLAGS =

src/c_kzg_util.c

@@ -39,6 +39,20 @@ C_KZG_RET c_kzg_malloc(void **x, size_t n) {
     return C_KZG_OK;
 }
 
+/**
+ * Allocate memory for an array of uint64_t.
+ *
+ * @remark Free the space later using `free()`.
+ *
+ * @param[out] x Pointer to the allocated space
+ * @param[in]  n The number of utin64_t to be allocated
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_MALLOC  Memory allocation failed
+ */
+C_KZG_RET new_uint64_array(uint64_t **x, size_t n) {
+    return c_kzg_malloc((void **)x, n * sizeof **x);
+}
+
 /**
  * Allocate memory for an array of field elements.
  *

src/c_kzg_util.h

@@ -20,6 +20,7 @@
 #include "c_kzg.h"
 
 C_KZG_RET c_kzg_malloc(void **p, size_t n);
+C_KZG_RET new_uint64_array(uint64_t **x, size_t n);
 C_KZG_RET new_fr_array(fr_t **x, size_t n);
 C_KZG_RET new_fr_array_2(fr_t ***x, size_t n);
 C_KZG_RET new_g1_array(g1_t **x, size_t n);

src/debug_util.c

@@ -47,6 +47,17 @@ void print_fr(const fr_t *a) {
     print_bytes_as_hex_le(b.b, 0, 32);
 }
 
+// Print a vector of Frs
+void print_frs(const char *s, const fr_t *x, uint64_t n) {
+    printf("\n----\n");
+    for (uint64_t i = 0; i < n; i++) {
+        printf("%s %lu: ", s, i);
+        print_fr(x + i);
+        printf("\n");
+    }
+    printf("----\n");
+}
+
 //
 // Fp Utilities
 //

src/debug_util.h

@@ -19,6 +19,7 @@
 
 // Fr utilities
 void print_fr(const fr_t *a);
+void print_frs(const char *s, const fr_t *x, uint64_t n);
 
 // G1 and G2 utilities
 void print_g1_bytes(byte p1[96]);

src/test_util.c

@@ -71,5 +71,16 @@ g1_t rand_g1() {
     return ret;
 }
 
+void shuffle(uint64_t *a, uint64_t n) {
+    uint64_t i = n, j, tmp;
+    while (i > 0) {
+        j = rand_uint64() % i;
+        i--;
+        tmp = a[j];
+        a[j] = a[i];
+        a[i] = tmp;
+    }
+}
+
 // Dummy function used to get the test-suite to print a title
 void title(void) {}

src/test_util.h

@@ -25,4 +25,5 @@ void generate_trusted_setup(g1_t *s1, g2_t *s2, const scalar_t *secret, const ui
 uint64_t rand_uint64();
 fr_t rand_fr();
 g1_t rand_g1();
+void shuffle(uint64_t *a, uint64_t n);
 void title(void);

src/utility.c

@@ -58,7 +58,16 @@ int log2_pow2(uint32_t n) {
     return r;
 }
 
-// Adapted from https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
+/**
+ * Return the next highest power of two.
+ *
+ * If @p v is already a power of two, it is returned as-is.
+ *
+ * Adapted from https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
+ *
+ * @param[in] v A 64-bit unsigned integer <= 2^31
+ * @return      The lowest power of two equal or larger than @p v
+ */
 uint64_t next_power_of_two(uint64_t v) {
     v--;
     v |= v >> 1;

src/zero_poly.c

@@ -0,0 +1,248 @@
+/*
+ * Copyright 2021 Benjamin Edgington
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ *  @file zero_poly.c
+ *
+ *  Methods for constructing zero polynomials and reconstructing polynomials from partial evaluations on a subgroup.
+ */
+
+#include "zero_poly.h"
+#include "fft_fr.h"
+#include "c_kzg_util.h"
+#include "utility.h"
+
+#include "debug_util.h"
+
+/**
+ * Process a slice of missing indices into a leaf.
+ *
+ * @param[out] dst         The resulting leaf, length @p len_dst
+ * @param[in]  len_dst     Length of the output leaf, @p dst
+ * @param[in]  indices     Array of missing indices of length @p len_indices
+ * @param[in]  len_indices Length of the missing indices array, @p indices
+ * @param[in]  stride      Stride length through the powers of the root of unity
+ * @param[in]  fs          The FFT settings previously initialised with #new_fft_settings
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_BADARGS Invalid parameters were supplied
+ */
+C_KZG_RET do_zero_poly_mul_leaf(fr_t *dst, uint64_t len_dst, const uint64_t *indices, uint64_t len_indices,
+                                uint64_t stride, const FFTSettings *fs) {
+
+    CHECK(len_dst >= len_indices + 1);
+
+    for (uint64_t i = 0; i < len_indices; i++) {
+        fr_t neg_di;
+        fr_negate(&neg_di, &fs->expanded_roots_of_unity[indices[i] * stride]);
+        dst[i] = neg_di;
+        if (i > 0) {
+            fr_add(&dst[i], &dst[i], &dst[i - 1]);
+            for (uint64_t j = i - 1; j > 0; j--) {
+                fr_mul(&dst[j], &dst[j], &neg_di);
+                fr_add(&dst[j], &dst[j], &dst[j - 1]);
+            }
+            fr_mul(&dst[0], &dst[0], &neg_di);
+        }
+    }
+
+    dst[len_indices] = fr_one;
+
+    for (uint64_t i = len_indices + 1; i < len_dst; i++) {
+        dst[i] = fr_zero;
+    }
+
+    return C_KZG_OK;
+}
+
+/**
+ * Copy @p p to @p out, padding to length @p n with zeros.
+ *
+ * @param[out] out     A copy of @p p padded to length @p n with zeros
+ * @param[in]  out_len The length of the desired output data, @p out
+ * @param[in]  p       The data to be copied and padded, length @p p_len
+ * @param[in]  p_len   The length of the data to be copied and padded
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_BADARGS Invalid parameters were supplied
+ */
+C_KZG_RET pad_p(fr_t *out, uint64_t out_len, const fr_t *p, uint64_t p_len) {
+    CHECK(out_len >= p_len);
+    for (uint64_t i = 0; i < p_len; i++) {
+        out[i] = p[i];
+    }
+    for (uint64_t i = p_len; i < out_len; i++) {
+        out[i] = fr_zero;
+    }
+    return C_KZG_OK;
+}
+
+/**
+ * Calculate the convolution of the input polynomials.
+ *
+ * Pad the polynomials in @p ps, perform FFTs, multiply the results together, and apply an inverse FFT to the result.
+ *
+ * @param[out] dst         The result of the convolution
+ * @param[in]  len_dst     Length of the output, a power of two
+ * @param      scratch     Scratch space of size at least 3 times the output size
+ * @param[in]  len_scratch Length of @p scratch, at least 3 x @p len_dst
+ * @param[in]  ps          Array of polynomial coefficients ps[@p len_ps][@p len_p]
+ * @param[in]  len_ps      The number of polynomials
+ * @param[in]  len_p       Array of lengths of each polynomial, size @p len_ps
+ * @param[in]  fs          The FFT settings previously initialised with #new_fft_settings
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_BADARGS Invalid parameters were supplied
+ * @retval C_CZK_ERROR   An internal error occurred
+ *
+ * @todo Check if we can make `ps` a proper 2d array rather than an array of pointers to arrays.
+ */
+C_KZG_RET reduce_leaves(fr_t *dst, uint64_t len_dst, fr_t *scratch, uint64_t len_scratch, blst_fr **ps, uint64_t len_ps,
+                        const uint64_t *len_p, const FFTSettings *fs) {
+    CHECK(is_power_of_two(len_dst));
+    CHECK(len_ps > 0);
+    CHECK(len_ps * len_p[0] <= len_dst); // needed, or just len_p <= len_dst
+    CHECK(len_scratch >= 3 * len_dst);
+
+    // Split `scratch` up into three equally sized working arrays
+    fr_t *p_padded = scratch;
+    fr_t *mul_eval_ps = scratch + len_dst;
+    fr_t *p_eval = scratch + 2 * len_dst;
+
+    // Do the last leaf first: it may be shorter than the others, and the padding can remain in place for the rest.
+    TRY(pad_p(p_padded, len_dst, ps[len_ps - 1], len_p[len_ps - 1]));
+    TRY(fft_fr(mul_eval_ps, p_padded, false, len_dst, fs));
+
+    for (uint64_t i = 0; i < len_ps - 1; i++) {
+        TRY(pad_p(p_padded, len_p[i], ps[i], len_p[i]));
+        TRY(fft_fr(p_eval, p_padded, false, len_dst, fs));
+        for (uint64_t j = 0; j < len_dst; j++) {
+            fr_mul(&mul_eval_ps[j], &mul_eval_ps[j], &p_eval[j]);
+        }
+    }
+
+    TRY(fft_fr(dst, mul_eval_ps, true, len_dst, fs));
+
+    return C_KZG_OK;
+}
+
+/**
+ * Calculate a polynomial that evaluates to zero for powers of roots of one that correspond to missing indices.
+ *
+ * Also calculates the FFT.
+ *
+ * @param[out] zero_eval Array length @p length (TODO: description)
+ * @param[out] zero_poly Array length @p length (TODO: description)
+ * @param[out] zero_poly_len The length of the resulting @p zero_poly
+ * @param[in]  length Length of the output arrays
+ * @param[in]  missing_indices Array length @p len_missing (TODO: description)
+ * @param[in]  len_missing Length of @p missing_indices
+ * @param[in]  fs     The FFT settings previously initialised with #new_fft_settings
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_BADARGS Invalid parameters were supplied
+ * @retval C_CZK_ERROR   An internal error occurred
+ * @retval C_CZK_MALLOC  Memory allocation failed
+ */
+C_KZG_RET zero_polynomial_via_multiplication(fr_t *zero_eval, fr_t *zero_poly, uint64_t *zero_poly_len, uint64_t length,
+                                             const uint64_t *missing_indices, uint64_t len_missing,
+                                             const FFTSettings *fs) {
+    if (len_missing == 0) {
+        return C_KZG_OK;
+    }
+    CHECK(length <= fs->max_width);
+    CHECK(is_power_of_two(length));
+
+    uint64_t domain_stride = fs->max_width / length;
+    uint64_t per_leaf_poly = 64;
+    uint64_t per_leaf = per_leaf_poly - 1;
+    uint64_t leaf_count = (len_missing + per_leaf - 1) / per_leaf;
+    uint64_t n = next_power_of_two(leaf_count * per_leaf_poly);
+
+    if (len_missing <= per_leaf) {
+        TRY(do_zero_poly_mul_leaf(zero_poly, length, missing_indices, len_missing, domain_stride, fs));
+        TRY(fft_fr(zero_eval, zero_poly, false, length, fs));
+        *zero_poly_len = length;
+    } else {
+        CHECK(n <= length);
+
+        // Work space for reducing the leaves - we could use `zero_poly` for this if it were always big enough.
+        fr_t *work = zero_poly;
+
+        // Build the leaves.
+
+        // Just allocate pointers here since we're re-using `work` for the leaf processing
+        // Combining leaves can be done mostly in-place, using a scratchpad.
+        fr_t **leaves, *scratch, *reduced;
+        uint64_t *leaf_lengths;
+        TRY(new_fr_array_2(&leaves, leaf_count));
+        TRY(new_uint64_array(&leaf_lengths, leaf_count));
+        uint64_t offset = 0, out_offset = 0, max = len_missing;
+        for (int i = 0; i < leaf_count; i++) {
+            uint64_t end = offset + per_leaf;
+            if (end > max) end = max;
+            leaves[i] = &work[out_offset];
+            leaf_lengths[i] = per_leaf_poly;
+            TRY(do_zero_poly_mul_leaf(leaves[i], leaf_lengths[i], &missing_indices[offset], end - offset, domain_stride,
+                                      fs));
+            offset += per_leaf;
+            out_offset += per_leaf_poly;
+        }
+
+        // Now reduce all the leaves to a single poly
+
+        int reduction_factor = 4; // must be a power of 2
+        TRY(new_fr_array(&scratch, n * 3));
+        // All the leaves are the same length, except possibly the last leaf, but that's ok.
+        while (leaf_count > 1) {
+            uint64_t reduced_count = (leaf_count + reduction_factor - 1) / reduction_factor;
+            uint64_t leaf_size = leaf_lengths[0];
+            for (uint64_t i = 0; i < reduced_count; i++) {
+                uint64_t start = i * reduction_factor;
+                uint64_t end = start + reduction_factor;
+                // E.g. if we *started* with 2 leaves, we won't have more than that since it is already a power
+                // of 2. If we had 3, it would have been rounded up anyway. So just pick the end
+                uint64_t out_end = end * leaf_size;
+                if (out_end > n) {
+                    out_end = n;
+                }
+                reduced = work + start * leaf_size;
+                uint64_t reduced_len = out_end - start * leaf_size;
+                if (end > leaf_count) {
+                    end = leaf_count;
+                }
+                uint64_t leaves_slice_len = end - start;
+                if (end > start + 1) {
+                    TRY(reduce_leaves(reduced, reduced_len, scratch, n * 3, &leaves[start], leaves_slice_len,
+                                      &leaf_lengths[start], fs));
+                    leaf_lengths[i] = reduced_len;
+                } else {
+                    leaf_lengths[i] = leaf_lengths[start];
+                }
+                leaves[i] = reduced;
+            }
+            leaf_count = reduced_count;
+        }
+
+        *zero_poly_len = leaf_lengths[0];
+        for (uint64_t i = 0; i < length; i++) {
+            zero_poly[i] = i < *zero_poly_len ? leaves[0][i] : fr_zero;
+        }
+        TRY(fft_fr(zero_eval, zero_poly, false, length, fs));
+
+        free(leaves);
+        free(leaf_lengths);
+        free(scratch);
+    }
+
+    return C_KZG_OK;
+}

src/zero_poly.h

@@ -0,0 +1,32 @@
+/*
+ * Copyright 2021 Benjamin Edgington
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ *  @file zero_poly.h
+ *
+ *  Methods for constructing zero polynomials and reconstructing polynomials from partial evaluations on a subgroup
+ */
+
+#include "c_kzg.h"
+#include "fft_common.h"
+
+C_KZG_RET do_zero_poly_mul_leaf(fr_t *dst, uint64_t len_dst, const uint64_t *indices, uint64_t len_indices,
+                                uint64_t stride, const FFTSettings *fs);
+C_KZG_RET reduce_leaves(fr_t *dst, uint64_t len_dst, fr_t *scratch, uint64_t len_scratch, blst_fr **ps, uint64_t len_ps,
+                        const uint64_t *len_p, const FFTSettings *fs);
+C_KZG_RET zero_polynomial_via_multiplication(fr_t *zero_eval, fr_t *zero_poly, uint64_t *zero_poly_len, uint64_t length,
+                                             const uint64_t *missing_indices, uint64_t len_missing,
+                                             const FFTSettings *fs);

src/zero_poly_test.c

@@ -0,0 +1,307 @@
+/*
+ * Copyright 2021 Benjamin Edgington
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "../inc/acutest.h"
+#include "c_kzg_util.h"
+#include "test_util.h"
+#include "zero_poly.h"
+#include "poly.h"
+#include "fft_fr.h"
+#include "debug_util.h"
+
+bool exists[16] = {true,  false, false, true, false, true, true,  false,
+                   false, false, true,  true, false, true, false, true};
+
+uint64_t expected_eval_u64[16][4] = {
+    {0xf675fcb368535efaL, 0xe702bee472f5a74cL, 0xb2f500c4418d44d8L, 0x204089b477319517L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x2be1bf25823353ecL, 0xe98177cae115131bL, 0xe0de4495f16788fbL, 0x37e5487beb15a91eL},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0xa8fd50243ec6f6daL, 0xb5863f0c04559733L, 0xbb55a8d735b8ceafL, 0x15856a55a6ba245bL},
+    {0x40d8d622337027e7L, 0xd0c41e3defe394e5L, 0x25d1a6848cfbe861L, 0x6615977f56ab9ad1L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x19b6d37343ac8596L, 0x9ac16b5b3f0c39eaL, 0x1938f2cc6f656899L, 0x2bc6a69eab7ebeadL},
+    {0x75ceddca83d9b1e4L, 0x69917e9ccac289bcL, 0x7564f74fd58cc97aL, 0x7215036c8f20939fL},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0088e6ba87233593L, 0xcc4a412d77455e7eL, 0x06ce406c147ada85L, 0x44275d7e26f9392cL},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x05ced2791378da2bL, 0xd16275df7a713f92L, 0x0cd24cf43668722dL, 0x22635b695b0fd198L}};
+
+uint64_t expected_poly_u64[16][4] = {
+    {0x6a20b4c8fbee018eL, 0x34c8bd90143c7a43L, 0xc4a72e43a8f20dbbL, 0x24c14de4b45f2d7bL},
+    {0xba227dc25dab47c2L, 0xfa1cdd366cf44de2L, 0x2920a9a04dd15d06L, 0x0174305e712df7baL},
+    {0xa3c8b170d759d6c4L, 0x846e2f5bfc241b81L, 0x1e4c5e807b5793eeL, 0x0758eca45c6dec8aL},
+    {0x2c280194f3795affL, 0x55035b9ba568dd4fL, 0x91dda79960525b60L, 0x3fbfd2edd4a105f3L},
+    {0x537cca635e26d630L, 0xaed6c42a88801d8fL, 0x41b2fdf16c422f7dL, 0x1d45a831fe3bf66eL},
+    {0x037b0169fc698ffdL, 0xe982a4842fc849f0L, 0xdd398294c762e031L, 0x4092c5b8416d2c2fL},
+    {0x19d4cdbb82bb00fbL, 0x5f31525ea0987c51L, 0xe80dcdb499dca94aL, 0x3aae0972562d375fL},
+    {0x91757d97669b7cc0L, 0x8e9c261ef753ba83L, 0x747c849bb4e1e1d8L, 0x02472328ddfa1df6L},
+    {0x0000000000000001L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L},
+    {0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L}};
+
+void test_reduce_leaves(void) {
+    FFTSettings fs;
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4));
+    fr_t from_tree_reduction[16], from_direct[9], scratch[48];
+
+    // Via reduce_leaves
+
+    fr_t *leaves[4];
+    fr_t leaf0[3], leaf1[3], leaf2[3], leaf3[3];
+    leaves[0] = leaf0;
+    leaves[1] = leaf1;
+    leaves[2] = leaf2;
+    leaves[3] = leaf3;
+    uint64_t leaf_lengths[] = {3, 3, 3, 3};
+    const uint64_t leaf_indices[4][2] = {{1, 3}, {7, 8}, {9, 10}, {12, 13}};
+    for (int i = 0; i < 4; i++) {
+        TEST_CHECK(C_KZG_OK == do_zero_poly_mul_leaf(leaves[i], 3, leaf_indices[i], 2, 1, &fs));
+    }
+    TEST_CHECK(C_KZG_OK == reduce_leaves(from_tree_reduction, 16, scratch, 48, leaves, 4, leaf_lengths, &fs));
+
+    // Direct
+    uint64_t indices[] = {1, 3, 7, 8, 9, 10, 12, 13};
+    TEST_CHECK(C_KZG_OK == do_zero_poly_mul_leaf(from_direct, 9, indices, 8, 1, &fs));
+
+    // Compare
+    for (int i = 0; i < 9; i++) {
+        TEST_CHECK(fr_equal(&from_tree_reduction[i], &from_direct[i]));
+    }
+
+    free_fft_settings(&fs);
+}
+
+void reduce_leaves_random(void) {
+    for (int scale = 5; scale < 13; scale++) {
+        for (int ii = 1; ii <= 7; ii++) {
+            float missing_ratio = 0.1 * ii;
+
+            FFTSettings fs;
+            new_fft_settings(&fs, scale);
+            uint64_t point_count = fs.max_width;
+            uint64_t missing_count = point_count * missing_ratio;
+
+            uint64_t *missing;
+            TEST_CHECK(C_KZG_OK == new_uint64_array(&missing, point_count));
+            for (uint64_t i = 0; i < point_count; i++) {
+                missing[i] = i;
+            }
+            shuffle(missing, point_count);
+
+            // Build the leaves
+            fr_t **leaves;
+            const int points_per_leaf = 63;
+            uint64_t indices[points_per_leaf];
+            uint64_t leaf_count = (missing_count + points_per_leaf - 1) / points_per_leaf;
+            uint64_t *leaf_lengths;
+            TEST_CHECK(C_KZG_OK == new_uint64_array(&leaf_lengths, leaf_count));
+            TEST_CHECK(C_KZG_OK == new_fr_array_2(&leaves, leaf_count));
+            for (uint64_t i = 0; i < leaf_count; i++) {
+                uint64_t start = i * points_per_leaf;
+                uint64_t end = start + points_per_leaf;
+                if (end > missing_count) end = missing_count;
+                uint64_t leaf_size = end - start;
+                TEST_CHECK(C_KZG_OK == new_fr_array(&leaves[i], leaf_size + 1));
+                for (int j = 0; j < leaf_size; j++) {
+                    indices[j] = missing[i * points_per_leaf + j];
+                }
+                leaf_lengths[i] = leaf_size + 1;
+                TEST_CHECK(C_KZG_OK == do_zero_poly_mul_leaf(leaves[i], leaf_lengths[i], indices, leaf_size, 1, &fs));
+            }
+
+            // From tree reduction
+            fr_t *from_tree_reduction, *scratch;
+            TEST_CHECK(C_KZG_OK == new_fr_array(&from_tree_reduction, point_count));
+            TEST_CHECK(C_KZG_OK == new_fr_array(&scratch, point_count * 3));
+            TEST_CHECK(C_KZG_OK == reduce_leaves(from_tree_reduction, point_count, scratch, point_count * 3, leaves,
+                                                 leaf_count, leaf_lengths, &fs));
+
+            // From direct
+            fr_t *from_direct;
+            TEST_CHECK(C_KZG_OK == new_fr_array(&from_direct, missing_count + 1));
+            TEST_CHECK(C_KZG_OK == do_zero_poly_mul_leaf(from_direct, missing_count + 1, missing, missing_count,
+                                                         fs.max_width / point_count, &fs));
+
+            for (uint64_t i = 0; i < missing_count + 1; i++) {
+                TEST_CHECK(fr_equal(&from_tree_reduction[i], &from_direct[i]));
+            }
+
+            free(from_tree_reduction);
+            free(from_direct);
+            free(scratch);
+            for (uint64_t i = 0; i < leaf_count; i++) {
+                free(leaves[i]);
+            }
+            free(leaves);
+            free(leaf_lengths);
+            free(missing);
+            free_fft_settings(&fs);
+        }
+    }
+}
+
+void check_test_data(void) {
+    FFTSettings fs;
+    poly expected_eval, expected_poly, tmp_poly;
+    new_poly(&expected_eval, 16);
+    new_poly(&expected_poly, 16);
+    new_poly(&tmp_poly, 16);
+    new_fft_settings(&fs, 4);
+
+    for (int i = 0; i < 16; i++) {
+        fr_from_uint64s(&expected_eval.coeffs[i], expected_eval_u64[i]);
+        fr_from_uint64s(&expected_poly.coeffs[i], expected_poly_u64[i]);
+    }
+
+    // Polynomial evalutes to zero at the expected places
+    for (int i = 0; i < 16; i++) {
+        if (!exists[i]) {
+            fr_t tmp;
+            eval_poly(&tmp, &expected_poly, &fs.expanded_roots_of_unity[i]);
+            TEST_CHECK(fr_is_zero(&tmp));
+            TEST_MSG("Failed for i = %d", i);
+        }
+    }
+
+    // This is a curiosity
+    for (int i = 1; i < 8; i++) {
+        fr_t tmp;
+        eval_poly(&tmp, &expected_eval, &fs.expanded_roots_of_unity[i]);
+        TEST_CHECK(fr_is_zero(&tmp));
+        TEST_MSG("Failed for i = %d", i);
+    }
+
+    // The eval poly is the FFT of the zero poly
+    TEST_CHECK(C_KZG_OK == fft_fr(tmp_poly.coeffs, expected_eval.coeffs, true, tmp_poly.length, &fs));
+    for (int i = 0; i < 16; i++) {
+        TEST_CHECK(fr_equal(&tmp_poly.coeffs[i], &expected_poly.coeffs[i]));
+        TEST_MSG("Failed for i = %d", i);
+    }
+
+    free_poly(&expected_poly);
+    free_poly(&expected_eval);
+    free_poly(&tmp_poly);
+    free_fft_settings(&fs);
+}
+
+void zero_poly_known(void) {
+    FFTSettings fs;
+    poly expected_eval, expected_poly, zero_eval, zero_poly;
+    uint64_t missing[16];
+    uint64_t len_missing = 0;
+    uint64_t zero_poly_len;
+    new_poly(&expected_eval, 16);
+    new_poly(&expected_poly, 16);
+    new_poly(&zero_eval, 16);
+    new_poly(&zero_poly, 16);
+    new_fft_settings(&fs, 4);
+
+    for (int i = 0; i < 16; i++) {
+        fr_from_uint64s(&expected_eval.coeffs[i], expected_eval_u64[i]);
+        fr_from_uint64s(&expected_poly.coeffs[i], expected_poly_u64[i]);
+        if (!exists[i]) {
+            missing[len_missing++] = i;
+        }
+    }
+
+    TEST_CHECK(C_KZG_OK == zero_polynomial_via_multiplication(zero_eval.coeffs, zero_poly.coeffs, &zero_poly_len,
+                                                              zero_eval.length, missing, len_missing, &fs));
+
+    TEST_CHECK(expected_poly.length == zero_poly_len);
+    TEST_MSG("Expected %lu, got %lu", expected_poly.length, zero_poly_len);
+
+    for (int i = 0; i < expected_eval.length; i++) {
+        TEST_CHECK(fr_equal(&expected_eval.coeffs[i], &zero_eval.coeffs[i]));
+        TEST_CHECK(fr_equal(&expected_poly.coeffs[i], &zero_poly.coeffs[i]));
+    }
+
+    free_poly(&expected_poly);
+    free_poly(&expected_eval);
+    free_poly(&zero_poly);
+    free_poly(&zero_eval);
+    free_fft_settings(&fs);
+}
+
+void zero_poly_random(void) {
+    for (int its = 0; its < 8; its++) {
+        for (int scale = 3; scale < 13; scale++) {
+            FFTSettings fs;
+            TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, scale));
+
+            uint64_t missing[fs.max_width];
+            int len_missing = 0;
+
+            for (int i = 0; i < fs.max_width; i++) {
+                if (rand() % 2) {
+                    missing[len_missing++] = i;
+                }
+            }
+
+            fr_t *zero_eval, *zero_poly;
+            uint64_t zero_poly_len;
+            TEST_CHECK(C_KZG_OK == new_fr_array(&zero_eval, fs.max_width));
+            TEST_CHECK(C_KZG_OK == new_fr_array(&zero_poly, fs.max_width));
+            TEST_CHECK(C_KZG_OK == zero_polynomial_via_multiplication(zero_eval, zero_poly, &zero_poly_len,
+                                                                      fs.max_width, missing, len_missing, &fs));
+
+            poly p;
+            p.length = zero_poly_len;
+            p.coeffs = zero_poly;
+            int ret = 0;
+            for (int i = 0; i < len_missing; i++) {
+                fr_t out;
+                eval_poly(&out, &p, &fs.expanded_roots_of_unity[missing[i]]);
+                ret = TEST_CHECK(fr_is_zero(&out));
+                TEST_MSG("Failed for missing[%d] = %lu", i, missing[i]);
+            }
+            TEST_MSG("Failed for scale %d", scale);
+
+            fr_t *zero_eval_fft;
+            TEST_CHECK(C_KZG_OK == new_fr_array(&zero_eval_fft, fs.max_width));
+            TEST_CHECK(C_KZG_OK == fft_fr(zero_eval_fft, zero_eval, true, fs.max_width, &fs));
+            for (uint64_t i = 0; i < zero_poly_len; i++) {
+                TEST_CHECK(fr_equal(&zero_poly[i], &zero_eval_fft[i]));
+            }
+            for (uint64_t i = zero_poly_len; i < fs.max_width; i++) {
+                TEST_CHECK(fr_is_zero(&zero_eval_fft[i]));
+            }
+
+            free(zero_poly);
+            free(zero_eval);
+            free(zero_eval_fft);
+            free_fft_settings(&fs);
+        }
+    }
+}
+
+TEST_LIST = {
+    {"ZERO_POLY_TEST", title},
+    {"test_reduce_leaves", test_reduce_leaves},
+    {"check_test_data", check_test_data},
+    {"reduce_leaves_random", reduce_leaves_random},
+    {"zero_poly_known", zero_poly_known},
+    {"zero_poly_random", zero_poly_random},
+    {NULL, NULL} /* zero record marks the end of the list */
+};