Implement DAS extension

README.md

@@ -11,6 +11,7 @@ Done so far:
   - Polynomial multi commitment and verification
   - [FK20](https://github.com/khovratovich/Kate/blob/master/Kate_amortized.pdf) single proof method (normal, and optimised for data availability)
   - FK20 multi proof method (normal, and optimised for data availability)
+  - Polynomial extension for data availability sampling
 
 ## Install
 

src/Makefile

@@ -1,7 +1,7 @@
 TESTS = bls12_381_test c_kzg_util_test fft_common_test fft_fr_test fft_g1_test \
 	fk20_proofs_test kzg_proofs_test poly_test
 BENCH = fft_fr_bench fft_g1_bench
-LIB_SRC = bls12_381.c c_kzg_util.c fft_common.c fft_fr.c fft_g1.c fk20_proofs.c kzg_proofs.c poly.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
 LIB_OBJ = $(LIB_SRC:.c=.o)
 
 CFLAGS =
@@ -24,12 +24,12 @@ libckzg.a: $(LIB_OBJ) Makefile
 	clang -Wall $(CFLAGS) -o $@ $*_test.c debug_util.o test_util.o libckzg.a -L../lib -lblst
 
 # Benchmarks
-%_bench: CFLAGS += -O3
-%_bench: %_bench.c bench_util.o $(LIB_OBJ) Makefile
-	clang -Wall $(CFLAGS) -o $@ $@.c bench_util.o $(LIB_OBJ) -L../lib -lblst
+%_bench: CFLAGS += -O
+%_bench: %_bench.c bench_util.o test_util.o $(LIB_OBJ) Makefile
+	clang -Wall $(CFLAGS) -o $@ $@.c bench_util.o test_util.o $(LIB_OBJ) -L../lib -lblst
 	./$@
 
-lib: CFLAGS += -O3
+lib: CFLAGS += -O
 lib: clean libckzg.a
 
 profilelib: CFLAGS += -fprofile-instr-generate -fcoverage-mapping

src/bench_util.c

@@ -14,33 +14,8 @@
  * limitations under the License.
  */
 
-#include <stdlib.h> // rand()
 #include "bench_util.h"
 
 unsigned long tdiff(timespec_t start, timespec_t end) {
     return (end.tv_sec - start.tv_sec) * NANO + (end.tv_nsec - start.tv_nsec);
 }
-
-uint64_t rand_uint64() {
-    uint64_t a = (uint64_t)rand();
-    uint64_t b = (uint64_t)rand();
-    return a << 32 | b;
-}
-
-fr_t rand_fr() {
-    fr_t ret;
-    uint64_t a[4];
-    a[0] = rand_uint64();
-    a[1] = rand_uint64();
-    a[2] = rand_uint64();
-    a[3] = rand_uint64();
-    fr_from_uint64s(&ret, a);
-    return ret;
-}
-
-g1_t rand_g1() {
-    g1_t ret;
-    fr_t random = rand_fr();
-    g1_mul(&ret, &g1_generator, &random);
-    return ret;
-}

src/bench_util.h

@@ -23,6 +23,3 @@ typedef struct timespec timespec_t;
 #define NSEC 1
 
 unsigned long tdiff(timespec_t start, timespec_t end);
-uint64_t rand_uint64();
-fr_t rand_fr();
-g1_t rand_g1();

src/das_extension.c

@@ -0,0 +1,111 @@
+/*
+ * 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 das_extension.c
+ *
+ * Perform polynomial extension for data availability sampling.
+ */
+
+#include "das_extension.h"
+
+/**
+ * Recursive implementation of #das_fft_extension.
+ *
+ * @param[in, out] ab     Input: values of the even indices. Output: values of the odd indices (in-place)
+ * @param[in]      n      The length of @p ab
+ * @param[in]      stride The step length through the roots of unity
+ * @param[in]      fs     The FFT settings previously initialised with #new_fft_settings
+ */
+static void das_fft_extension_stride(fr_t *ab, uint64_t n, uint64_t stride, const FFTSettings *fs) {
+
+    if (n < 2) return;
+
+    if (n == 2) {
+        fr_t x, y, tmp;
+        fr_add(&x, &ab[0], &ab[1]);
+        fr_sub(&y, &ab[0], &ab[1]);
+        fr_mul(&tmp, &y, &fs->expanded_roots_of_unity[stride]);
+        fr_add(&ab[0], &x, &tmp);
+        fr_sub(&ab[1], &x, &tmp);
+    } else {
+        uint64_t half = n, halfhalf = half / 2;
+        fr_t *ab_half_0s = ab;
+        fr_t *ab_half_1s = ab + halfhalf;
+
+        // Modify ab_half_* in-place, rather than allocating L0 and L1 arrays.
+        // L0[i] = (((a_half0 + a_half1) % modulus) * inv2) % modulus
+        // R0[i] = (((a_half0 - L0[i]) % modulus) * inverse_domain[i * 2]) % modulus
+        for (uint64_t i = 0; i < halfhalf; i++) {
+            fr_t tmp1, tmp2;
+            fr_t *a_half_0 = ab_half_0s + i;
+            fr_t *a_half_1 = ab_half_1s + i;
+            fr_add(&tmp1, a_half_0, a_half_1);
+            fr_sub(&tmp2, a_half_0, a_half_1);
+            fr_mul(a_half_1, &tmp2, &fs->reverse_roots_of_unity[i * 2 * stride]);
+            *a_half_0 = tmp1;
+        }
+
+        // Recurse
+        das_fft_extension_stride(ab_half_0s, halfhalf, stride * 2, fs);
+        das_fft_extension_stride(ab_half_1s, halfhalf, stride * 2, fs);
+
+        // The odd deduced outputs are written to the output array already, but then updated in-place
+        // L1 = b[:halfHalf]
+        // R1 = b[halfHalf:]
+
+        for (uint64_t i = 0; i < halfhalf; i++) {
+            fr_t y_times_root;
+            fr_t x = ab_half_0s[i];
+            fr_t y = ab_half_1s[i];
+            fr_mul(&y_times_root, &y, &fs->expanded_roots_of_unity[(1 + 2 * i) * stride]);
+            // write outputs in place, avoid unnecessary list allocations
+            fr_add(&ab_half_0s[i], &x, &y_times_root);
+            fr_sub(&ab_half_1s[i], &x, &y_times_root);
+        }
+    }
+}
+
+/**
+ * Perform polynomial extension for data availability sampling.
+ *
+ * The input is the even-numbered indices, which is replaced by the odd indices required to make the right half of the
+ * coefficients of the inverse FFT of the combined indices zero.
+ *
+ * @remark The input (even index) values are replace by the output (odd index) values.
+ *
+ * @param[in, out] vals Input: values of the even indices. Output: values of the odd indices (in place)
+ * @param[in]      n    The length of @p vals
+ * @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 das_fft_extension(fr_t *vals, uint64_t n, const FFTSettings *fs) {
+    fr_t invlen;
+
+    CHECK(n * 2 <= fs->max_width);
+    CHECK(n >= 2);
+
+    das_fft_extension_stride(vals, n, 1, fs);
+
+    fr_from_uint64(&invlen, n);
+    fr_inv(&invlen, &invlen);
+    for (uint64_t i = 0; i < n; i++) {
+        fr_mul(&vals[i], &vals[i], &invlen);
+    }
+
+    return C_KZG_OK;
+}

src/das_extension.h

@@ -0,0 +1,26 @@
+/*
+ * 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 das_extension.h
+ *
+ * Perform polynomial extension for data availability sampling.
+ */
+
+#include "c_kzg.h"
+#include "fft_common.h"
+
+C_KZG_RET das_fft_extension(fr_t *vals, uint64_t n, const FFTSettings *fs);

src/das_extension_test.c

@@ -0,0 +1,107 @@
+/*
+ * 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 "fft_fr.h"
+#include "das_extension.h"
+
+void das_extension_test_known(void) {
+    FFTSettings fs;
+    uint64_t half;
+    fr_t *data;
+    // The expected output (from go-kzg):
+    const uint64_t expected_u[8][4] = {
+        {0x740e9eb4cef44b7fL, 0xb66dfe6438d1316aL, 0xc62c567f7b9c5e97L, 0x5a4f75d3eb7951c4L},
+        {0x52617bb5c60d8fabL, 0x0e225ff8d6c658d9L, 0x73ce7b30718ed2aeL, 0x0d94457acc0cb888L},
+        {0x8bf0614a310bb489L, 0xb14ca59c512a2a94L, 0xdfbbb4ba1802749cL, 0x199e317f3e242b82L},
+        {0x98cbae3b8f66f769L, 0x69198f8ae9ed978dL, 0x360058041982845aL, 0x3171960f86bb881dL},
+        {0x8bf0614a310bb489L, 0xb14ca59c512a2a94L, 0xdfbbb4ba1802749cL, 0x199e317f3e242b82L},
+        {0x52617bb5c60d8fabL, 0x0e225ff8d6c658d9L, 0x73ce7b30718ed2aeL, 0x0d94457acc0cb888L},
+        {0x740e9eb4cef44b7fL, 0xb66dfe6438d1316aL, 0xc62c567f7b9c5e97L, 0x5a4f75d3eb7951c4L},
+        {0xc2735a57e47de950L, 0xa665548b548a12beL, 0x3040233ff907b7f0L, 0x2753864e0ac8841bL}};
+
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4));
+    half = fs.max_width / 2;
+
+    TEST_CHECK(C_KZG_OK == new_fr(&data, half));
+    for (uint64_t i = 0; i < half; i++) {
+        fr_from_uint64(data + i, i);
+    }
+
+    TEST_CHECK(C_KZG_OK == das_fft_extension(data, half, &fs));
+
+    // Check against the expected values
+    for (uint64_t i = 0; i < 8; i++) {
+        fr_t expected;
+        fr_from_uint64s(&expected, expected_u[i]);
+        TEST_CHECK(fr_equal(&expected, data + i));
+    }
+
+    free(data);
+    free_fft_settings(&fs);
+}
+
+// Caution: uses random data
+void das_extension_test_random(void) {
+    FFTSettings fs;
+    fr_t *even_data, *odd_data, *data, *coeffs;
+    for (int scale = 4; scale < 10; scale++) {
+        TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, scale));
+        TEST_CHECK(C_KZG_OK == new_fr(&even_data, fs.max_width / 2));
+        TEST_CHECK(C_KZG_OK == new_fr(&odd_data, fs.max_width / 2));
+        TEST_CHECK(C_KZG_OK == new_fr(&data, fs.max_width));
+        TEST_CHECK(C_KZG_OK == new_fr(&coeffs, fs.max_width));
+
+        for (int rep = 0; rep < 4; rep++) {
+
+            // Make random input data, and save a copy of it
+            for (int i = 0; i < fs.max_width / 2; i++) {
+                even_data[i] = rand_fr();
+                odd_data[i] = even_data[i];
+            }
+
+            // Extend the odd data
+            TEST_CHECK(C_KZG_OK == das_fft_extension(odd_data, fs.max_width / 2, &fs));
+
+            // Reconstruct the data
+            for (int i = 0; i < fs.max_width; i += 2) {
+                data[i] = even_data[i / 2];
+                data[i + 1] = odd_data[i / 2];
+            }
+            TEST_CHECK(C_KZG_OK == fft_fr(coeffs, data, true, fs.max_width, &fs));
+
+            // Second half of the coefficients should be all zeros
+            for (int i = fs.max_width / 2; i < fs.max_width; i++) {
+                TEST_CHECK(fr_is_zero(&coeffs[i]));
+            }
+        }
+
+        free(even_data);
+        free(odd_data);
+        free(data);
+        free(coeffs);
+        free_fft_settings(&fs);
+    }
+}
+
+TEST_LIST = {
+    {"DAS_EXTENSION_TEST", title},
+    {"das_extension_test_known", das_extension_test_known},
+    {"das_extension_test_random", das_extension_test_random},
+    {NULL, NULL} /* zero record marks the end of the list */
+};

src/fft_fr_bench.c

@@ -19,6 +19,7 @@
 #include <assert.h> // assert()
 #include <unistd.h> // EXIT_SUCCESS/FAILURE
 #include "bench_util.h"
+#include "test_util.h"
 #include "fft_fr.h"
 
 // Run the benchmark for `max_seconds` and return the time per iteration in nanoseconds.

src/fft_g1_bench.c

@@ -19,6 +19,7 @@
 #include <assert.h> // assert()
 #include <unistd.h> // EXIT_SUCCESS/FAILURE
 #include "bench_util.h"
+#include "test_util.h"
 #include "fft_g1.h"
 
 // Run the benchmark for `max_seconds` and return the time per iteration in nanoseconds.

src/test_util.c

@@ -14,6 +14,7 @@
  * limitations under the License.
  */
 
+#include <stdlib.h> // rand()
 #include "test_util.h"
 
 void generate_trusted_setup(g1_t *s1, g2_t *s2, const scalar_t *secret, const uint64_t n) {
@@ -28,5 +29,29 @@ void generate_trusted_setup(g1_t *s1, g2_t *s2, const scalar_t *secret, const ui
     }
 }
 
+uint64_t rand_uint64() {
+    uint64_t a = (uint64_t)rand();
+    uint64_t b = (uint64_t)rand();
+    return a << 32 | b;
+}
+
+fr_t rand_fr() {
+    fr_t ret;
+    uint64_t a[4];
+    a[0] = rand_uint64();
+    a[1] = rand_uint64();
+    a[2] = rand_uint64();
+    a[3] = rand_uint64();
+    fr_from_uint64s(&ret, a);
+    return ret;
+}
+
+g1_t rand_g1() {
+    g1_t ret;
+    fr_t random = rand_fr();
+    g1_mul(&ret, &g1_generator, &random);
+    return ret;
+}
+
 // Dummy function used to get the test-suite to print a title
 void title(void) {}

src/test_util.h

@@ -22,4 +22,7 @@ static const scalar_t secret = {0xa4, 0x73, 0x31, 0x95, 0x28, 0xc8, 0xb6, 0xea,
                                 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
 
 void generate_trusted_setup(g1_t *s1, g2_t *s2, const scalar_t *secret, const uint64_t n);
+uint64_t rand_uint64();
+fr_t rand_fr();
+g1_t rand_g1();
 void title(void);