From a697775c4f26fabf63d0f9b7a67f6833d0dbd0ae Mon Sep 17 00:00:00 2001 From: Ben Edgington Date: Fri, 26 Feb 2021 14:38:04 +0000 Subject: [PATCH] Implement zero polynomial - passes tests, but could be tidier --- src/Makefile | 4 +- src/c_kzg_util.c | 14 ++ src/c_kzg_util.h | 1 + src/debug_util.c | 11 ++ src/debug_util.h | 1 + src/test_util.c | 11 ++ src/test_util.h | 1 + src/utility.c | 11 +- src/zero_poly.c | 248 ++++++++++++++++++++++++++++++++++ src/zero_poly.h | 32 +++++ src/zero_poly_test.c | 307 +++++++++++++++++++++++++++++++++++++++++++ 11 files changed, 638 insertions(+), 3 deletions(-) create mode 100644 src/zero_poly.c create mode 100644 src/zero_poly.h create mode 100644 src/zero_poly_test.c diff --git a/src/Makefile b/src/Makefile index 27b99f0..fc3bb5c 100644 --- a/src/Makefile +++ b/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 = diff --git a/src/c_kzg_util.c b/src/c_kzg_util.c index 35a7346..771f0f4 100644 --- a/src/c_kzg_util.c +++ b/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. * diff --git a/src/c_kzg_util.h b/src/c_kzg_util.h index 14b4411..d057cfb 100644 --- a/src/c_kzg_util.h +++ b/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); diff --git a/src/debug_util.c b/src/debug_util.c index 8c13d31..74d72e3 100644 --- a/src/debug_util.c +++ b/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 // diff --git a/src/debug_util.h b/src/debug_util.h index 85ec244..ba4bd71 100644 --- a/src/debug_util.h +++ b/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]); diff --git a/src/test_util.c b/src/test_util.c index ff29532..0e78007 100644 --- a/src/test_util.c +++ b/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) {} diff --git a/src/test_util.h b/src/test_util.h index 1b3e1b1..f280f77 100644 --- a/src/test_util.h +++ b/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); diff --git a/src/utility.c b/src/utility.c index 83bfda8..03c090a 100644 --- a/src/utility.c +++ b/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; diff --git a/src/zero_poly.c b/src/zero_poly.c new file mode 100644 index 0000000..83a211f --- /dev/null +++ b/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; +} diff --git a/src/zero_poly.h b/src/zero_poly.h new file mode 100644 index 0000000..ffc0d85 --- /dev/null +++ b/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); diff --git a/src/zero_poly_test.c b/src/zero_poly_test.c new file mode 100644 index 0000000..27993ee --- /dev/null +++ b/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 */ +};