Move FFT functions to fft.{c,h} files

bindings/node.js/Makefile

@@ -36,12 +36,14 @@ build: install clean
 	@cp ../../src/common.c deps/c-kzg
 	@cp ../../src/eip4844.c deps/c-kzg
 	@cp ../../src/eip7594.c deps/c-kzg
+	@cp ../../src/fft.c deps/c-kzg
 	@cp ../../src/setup.c deps/c-kzg
 	@# Copy header files
 	@cp ../../src/ckzg.h deps/c-kzg
 	@cp ../../src/common.h deps/c-kzg
 	@cp ../../src/eip4844.h deps/c-kzg
 	@cp ../../src/eip7594.h deps/c-kzg
+	@cp ../../src/fft.h deps/c-kzg
 	@cp ../../src/setup.h deps/c-kzg
 	@# Copy trusted setup
 	@cp ../../src/trusted_setup.txt deps/c-kzg

src/Makefile

@@ -40,7 +40,8 @@ BLST_BUILDSCRIPT_FLAGS = -D__BLST_PORTABLE__
 LIBS = $(BLST_LIBRARY)
 
 # List of individual implementation files.
-IMPL_FILES = common.c eip4844.c eip7594.c setup.c
+IMPL_FILES = common.c debug.c eip4844.c eip7594.c fft.c setup.c
+HEADER_FILES = $(patsubst %.c, %.h, $(IMPL_FILES))
 
 ###############################################################################
 # Core
@@ -162,7 +163,7 @@ analyze: $(IMPL_FILES)
 
 .PHONY: format
 format:
-	@clang-format -i --sort-includes common.* eip4844.* eip7594.* setup.* debug.* ckzg.* tests.c
+	@clang-format -i --sort-includes $(IMPL_FILES) $(HEADER_FILES) tests.c
 
 .PHONY: clean
 clean:

src/ckzg.c

@@ -17,4 +17,5 @@
 #include "common.c"
 #include "eip4844.c"
 #include "eip7594.c"
+#include "fft.c"
 #include "setup.c"

src/common.c

@@ -608,96 +608,3 @@ C_KZG_RET g1_lincomb_fast(g1_t *out, const g1_t *p, const fr_t *coeffs, size_t l
     c_kzg_free(scalars);
     return ret;
 }
-
-/**
- * Fast Fourier Transform.
- *
- * Recursively divide and conquer.
- *
- * @param[out] out          The results (length `n`)
- * @param[in]  in           The input data (length `n * stride`)
- * @param[in]  stride       The input data stride
- * @param[in]  roots        Roots of unity (length `n * roots_stride`)
- * @param[in]  roots_stride The stride interval among the roots of unity
- * @param[in]  n            Length of the FFT, must be a power of two
- */
-static void g1_fft_fast(
-    g1_t *out, const g1_t *in, uint64_t stride, const fr_t *roots, uint64_t roots_stride, uint64_t n
-) {
-    g1_t y_times_root;
-    uint64_t half = n / 2;
-    if (half > 0) { /* Tunable parameter */
-        g1_fft_fast(out, in, stride * 2, roots, roots_stride * 2, half);
-        g1_fft_fast(out + half, in + stride, stride * 2, roots, roots_stride * 2, half);
-        for (uint64_t i = 0; i < half; i++) {
-            /* If the point is infinity, we can skip the calculation */
-            if (blst_p1_is_inf(&out[i + half])) {
-                out[i + half] = out[i];
-            } else {
-                /* If the scalar is one, we can skip the multiplication */
-                if (fr_is_one(&roots[i * roots_stride])) {
-                    y_times_root = out[i + half];
-                } else {
-                    g1_mul(&y_times_root, &out[i + half], &roots[i * roots_stride]);
-                }
-                g1_sub(&out[i + half], &out[i], &y_times_root);
-                blst_p1_add_or_double(&out[i], &out[i], &y_times_root);
-            }
-        }
-    } else {
-        *out = *in;
-    }
-}
-
-/**
- * The entry point for forward FFT over G1 points.
- *
- * @param[out]  out The results (array of length n)
- * @param[in]   in  The input data (array of length n)
- * @param[in]   n   Length of the arrays
- * @param[in]   s   The trusted setup
- *
- * @remark The array lengths must be a power of two.
- * @remark Use g1_ifft() for inverse transformation.
- */
-C_KZG_RET g1_fft(g1_t *out, const g1_t *in, size_t n, const KZGSettings *s) {
-    /* Ensure the length is valid */
-    if (n > s->max_width || !is_power_of_two(n)) {
-        return C_KZG_BADARGS;
-    }
-
-    uint64_t stride = s->max_width / n;
-    g1_fft_fast(out, in, 1, s->expanded_roots_of_unity, stride, n);
-
-    return C_KZG_OK;
-}
-
-/**
- * The entry point for inverse FFT over G1 points.
- *
- * @param[out]  out The results (array of length n)
- * @param[in]   in  The input data (array of length n)
- * @param[in]   n   Length of the arrays
- * @param[in]   s   The trusted setup
- *
- * @remark The array lengths must be a power of two.
- * @remark Use g1_fft() for forward transformation.
- */
-C_KZG_RET g1_ifft(g1_t *out, const g1_t *in, size_t n, const KZGSettings *s) {
-    /* Ensure the length is valid */
-    if (n > s->max_width || !is_power_of_two(n)) {
-        return C_KZG_BADARGS;
-    }
-
-    uint64_t stride = s->max_width / n;
-    g1_fft_fast(out, in, 1, s->reverse_roots_of_unity, stride, n);
-
-    fr_t inv_len;
-    fr_from_uint64(&inv_len, n);
-    blst_fr_eucl_inverse(&inv_len, &inv_len);
-    for (uint64_t i = 0; i < n; i++) {
-        g1_mul(&out[i], &out[i], &inv_len);
-    }
-
-    return C_KZG_OK;
-}

src/common.h

@@ -256,8 +256,6 @@ void g1_mul(g1_t *out, const g1_t *a, const fr_t *b);
 bool pairings_verify(const g1_t *a1, const g2_t *a2, const g1_t *b1, const g2_t *b2);
 void g1_lincomb_naive(g1_t *out, const g1_t *p, const fr_t *coeffs, uint64_t len);
 C_KZG_RET g1_lincomb_fast(g1_t *out, const g1_t *p, const fr_t *coeffs, size_t len);
-C_KZG_RET g1_fft(g1_t *out, const g1_t *in, size_t n, const KZGSettings *s);
-C_KZG_RET g1_ifft(g1_t *out, const g1_t *in, size_t n, const KZGSettings *s);
 
 #ifdef __cplusplus
 }

src/eip7594.c

@@ -18,6 +18,7 @@
 
 #include "blst.h"
 #include "common.h"
+#include "fft.h"
 
 #include <inttypes.h>
 #include <stdio.h>
@@ -107,92 +108,6 @@ static uint32_t reverse_bits_limited(uint32_t n, uint32_t value) {
     return reverse_bits(value) >> unused_bit_len;
 }
 
-////////////////////////////////////////////////////////////////////////////////////////////////////
-// Fast Fourier Transform
-////////////////////////////////////////////////////////////////////////////////////////////////////
-
-/**
- * Fast Fourier Transform.
- *
- * Recursively divide and conquer.
- *
- * @param[out]  out             The results (length `n`)
- * @param[in]   in              The input data (length `n * stride`)
- * @param[in]   stride          The input data stride
- * @param[in]   roots           Roots of unity (length `n * roots_stride`)
- * @param[in]   roots_stride    The stride interval among the roots of unity
- * @param[in]   n               Length of the FFT, must be a power of two
- */
-static void fr_fft_fast(
-    fr_t *out, const fr_t *in, size_t stride, const fr_t *roots, size_t roots_stride, size_t n
-) {
-    size_t half = n / 2;
-    if (half > 0) {
-        fr_t y_times_root;
-        fr_fft_fast(out, in, stride * 2, roots, roots_stride * 2, half);
-        fr_fft_fast(out + half, in + stride, stride * 2, roots, roots_stride * 2, half);
-        for (size_t i = 0; i < half; i++) {
-            blst_fr_mul(&y_times_root, &out[i + half], &roots[i * roots_stride]);
-            blst_fr_sub(&out[i + half], &out[i], &y_times_root);
-            blst_fr_add(&out[i], &out[i], &y_times_root);
-        }
-    } else {
-        *out = *in;
-    }
-}
-
-/**
- * The entry point for forward FFT over field elements.
- *
- * @param[out]  out The results (array of length n)
- * @param[in]   in  The input data (array of length n)
- * @param[in]   n   Length of the arrays
- * @param[in]   s   The trusted setup
- *
- * @remark The array lengths must be a power of two.
- * @remark Use fr_ifft() for inverse transformation.
- */
-static C_KZG_RET fr_fft(fr_t *out, const fr_t *in, size_t n, const KZGSettings *s) {
-    /* Ensure the length is valid */
-    if (n > s->max_width || !is_power_of_two(n)) {
-        return C_KZG_BADARGS;
-    }
-
-    size_t stride = s->max_width / n;
-    fr_fft_fast(out, in, 1, s->expanded_roots_of_unity, stride, n);
-
-    return C_KZG_OK;
-}
-
-/**
- * The entry point for inverse FFT over field elements.
- *
- * @param[out]  out The results (array of length n)
- * @param[in]   in  The input data (array of length n)
- * @param[in]   n   Length of the arrays
- * @param[in]   s   The trusted setup
- *
- * @remark The array lengths must be a power of two.
- * @remark Use fr_fft() for forward transformation.
- */
-static C_KZG_RET fr_ifft(fr_t *out, const fr_t *in, size_t n, const KZGSettings *s) {
-    /* Ensure the length is valid */
-    if (n > s->max_width || !is_power_of_two(n)) {
-        return C_KZG_BADARGS;
-    }
-
-    size_t stride = s->max_width / n;
-    fr_fft_fast(out, in, 1, s->reverse_roots_of_unity, stride, n);
-
-    fr_t inv_len;
-    fr_from_uint64(&inv_len, n);
-    blst_fr_inverse(&inv_len, &inv_len);
-    for (size_t i = 0; i < n; i++) {
-        blst_fr_mul(&out[i], &out[i], &inv_len);
-    }
-    return C_KZG_OK;
-}
-
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Vanishing Polynomial
 ////////////////////////////////////////////////////////////////////////////////////////////////////