diff --git a/src/Makefile b/src/Makefile index 87e91c6..6c64ba0 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,10 +1,14 @@ -tests = fft_fr_test +tests = fft_fr_test fft_util_test .PRECIOUS: %.o %.o: %.c %.h c-kzg.h clang -Wall -c $*.c +fft_fr_test: fft_fr.o fft_fr_test.c test_util.o fft_util.o + clang -Wall -o $@ $@.c test_util.o fft_fr.o fft_util.o -L../lib -lblst + ./$@ + %_test: %.o %_test.c test_util.o clang -Wall -o $@ $@.c test_util.o $*.o -L../lib -lblst ./$@ diff --git a/src/c-kzg.h b/src/c-kzg.h index 6a66d97..bbc4978 100644 --- a/src/c-kzg.h +++ b/src/c-kzg.h @@ -15,4 +15,5 @@ */ #include +#include #include "../inc/blst.h" diff --git a/src/fft_fr.c b/src/fft_fr.c index 8b1ef0b..439049b 100644 --- a/src/fft_fr.c +++ b/src/fft_fr.c @@ -16,62 +16,6 @@ #include "fft_fr.h" -bool is_one(const blst_fr *fr_p) { - uint64_t a[4]; - blst_uint64_from_fr(a, fr_p); - return a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 0; -} - -bool is_power_of_two(uint64_t n) { - return (n & (n - 1)) == 0; -} - -void fr_from_uint64(blst_fr *a, uint64_t n) { - uint64_t vals[] = {n, 0, 0, 0}; - blst_fr_from_uint64(a, vals); -} - -// Returns an array of powers of the root of unity -// Allocates space for the array that needs to be freed later -blst_fr *expand_root_of_unity(blst_fr *root_of_unity, uint64_t width) { - blst_fr *roots = malloc((width + 1) * sizeof(blst_fr)); - roots[0] = one; - roots[1] = *root_of_unity; - - for (int i = 2; !is_one(&roots[i - 1]); i++) { - assert(i <= width); - blst_fr_mul(&roots[i], &roots[i - 1], root_of_unity); - } - assert(is_one(&roots[width])); - - return roots; -} - -// Return a reversed copy of the list of Fr provided -// `width` is one less than the length of `r` -// Allocates space for the array that needs to be freed later -blst_fr *reverse(blst_fr *r, uint64_t width) { - blst_fr *rr = malloc((width + 1) * sizeof(blst_fr)); - for (int i = 0; i <= width; i++) { - rr[i] = r[width - i]; - } - return rr; -} - -FFTSettings new_fft_settings(unsigned int max_scale) { - FFTSettings s; - s.max_width = (uint64_t)1 << max_scale; - blst_fr_from_uint64(&s.root_of_unity, scale2_root_of_unity[max_scale]); - s.expanded_roots_of_unity = expand_root_of_unity(&s.root_of_unity, s.max_width); - s.reverse_roots_of_unity = reverse(s.expanded_roots_of_unity, s.max_width); - return s; -} - -void free_fft_settings(FFTSettings *s) { - free(s->expanded_roots_of_unity); - free(s->reverse_roots_of_unity); -} - // Slow Fourier Transform (simple, good for small sizes) void slow_ft(blst_fr *out, blst_fr *in, uint64_t offset, uint64_t stride, blst_fr *roots, uint64_t roots_stride, uint64_t l) { blst_fr v, last, tmp; diff --git a/src/fft_fr.h b/src/fft_fr.h index 4b1a9f2..1006bca 100644 --- a/src/fft_fr.h +++ b/src/fft_fr.h @@ -14,70 +14,9 @@ * limitations under the License. */ -#include -#include -#include #include "c-kzg.h" +#include "fft_util.h" -// This is 1 in Blst's `blst_fr` limb representation. Crazy but true. -static const blst_fr one = - {0x00000001fffffffeL, 0x5884b7fa00034802L, 0x998c4fefecbc4ff5L, 0x1824b159acc5056fL}; - -// MODULUS = 52435875175126190479447740508185965837690552500527637822603658699938581184513 -// PRIMITIVE_ROOT = 5 -// [pow(PRIMITIVE_ROOT, (MODULUS - 1) // (2**i), MODULUS) for i in range(32)] -// -// These are not in `blst_fr` limb format and must be converted via `blst_fr_from_uint64()` -static const uint64_t scale2_root_of_unity[][4] = - { - {0x0000000000000001L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L}, - {0xffffffff00000000L, 0x53bda402fffe5bfeL, 0x3339d80809a1d805L, 0x73eda753299d7d48L}, - {0x0001000000000000L, 0xec03000276030000L, 0x8d51ccce760304d0L, 0x0000000000000000L}, - {0x8dd702cb688bc087L, 0xa032824078eaa4feL, 0xa733b23a98ca5b22L, 0x3f96405d25a31660L}, - {0x95df4b360411fe73L, 0x1ef4e672ebc1e1bbL, 0x6e92a9c460afca4aL, 0x4f2c596e753e4fccL}, - {0x9733f3a6ba60eaa6L, 0xbd7fdf9c77487ae7L, 0xd84f8612c8b6cc00L, 0x476fa2fb6162ffabL}, - {0xd2fc5e69ac5db47fL, 0xa12a70a615d0b8e4L, 0x293b1d67bc8de5d9L, 0x0e4840ac57f86f5eL}, - {0xb750da4cab28e208L, 0x501dff643be95635L, 0x8cbe2437f0b4b276L, 0x07d0c802a94a946eL}, - {0x2cabadec2fe322b8L, 0x752c84f315412560L, 0x32a732ae1a3b0aefL, 0x2e95da59a33dcbf2L}, - {0x33811ea1fe0c65f4L, 0x15c1ad4c687f28a2L, 0xecfbede342dee7f4L, 0x1bb466679a5d88b1L}, - {0xd58a5af42d010ff9L, 0x79efd6b0570bf109L, 0x3ed6d55a6350721dL, 0x2f27b09858f43cefL}, - {0x74a1f6718c130477L, 0xa534af14b61e0abeL, 0xeb674a1a620890d7L, 0x43527a8bca252472L}, - {0x450d9f977ea8ee05L, 0x565af17137d56fc0L, 0xe155cb4893f9e9acL, 0x110cebd0c8e9101bL}, - {0x23c9159959a0be92L, 0x87d188ce7a027759L, 0x70491431cab3c3ccL, 0x0ac00eb8b3f7f8daL}, - {0x13e96ade69583404L, 0x82c057275306243dL, 0x77e48bf529ca9f2aL, 0x50646ac81fe19595L}, - {0xe6a354dda97eccd4L, 0x39929d2e88fbbc57L, 0xa22ba63dd6e7b1c8L, 0x42c22911f5f07f43L}, - {0x137b458acfc35f7aL, 0x0caba63a29c01b06L, 0x0409ee987a02402cL, 0x6709c6cd56aa725bL}, - {0x10251f7d8831e03eL, 0x77d85a937ff858ecL, 0xebe905bd4fb9ac5cL, 0x05deb333f8727901L}, - {0xbf87b689b9009408L, 0x4f730e7ddd3ccc96L, 0xfd7f05ba4610300cL, 0x5ef5e8db0b8ac903L}, - {0x6499688417cd0c14L, 0xa672867368812f7fL, 0x2e1d9a1922cc3253L, 0x3a689e83aa0a1d80L}, - {0x20b53cbe41144deaL, 0x870c46fac2f0fcbdL, 0x556c35f6537d6971L, 0x3436287f5f686d91L}, - {0x007e082a436ba2e7L, 0x67c6630f9116e877L, 0x36f8f165fb4460f7L, 0x6eee34d57e7046e0L}, - {0xc5b670eea53a56d1L, 0x127d1f4253037d7bL, 0x57d4257ea722c2e2L, 0x03ae26a333cbd838L}, - {0x1e91484876504cf8L, 0x55bbbf1eb63edd02L, 0xbcdafec84e55aa02L, 0x5145c4cd2dc0beb0L}, - {0x5b90153a1ab70e2cL, 0x8deffa3175fb0ab8L, 0xc553ae2346900c95L, 0x1d31dcdc6bd3118cL}, - {0x801c894c59a2e8ebL, 0xbc535c5ce12fc974L, 0x95508d2747d39803L, 0x16d9d3cdac5d094fL}, - {0x810fa372cca1d8beL, 0xc67b8c2882e0bfa7L, 0xdbb4edf0e2d35bc2L, 0x712d15805087c995L}, - {0xeb162203fd88f133L, 0xac96c38ff010ea74L, 0x4307987fe64cfc70L, 0x350fe98d37b7a114L}, - {0xaba2f51842f2a254L, 0x4d7f3c3aa71efc0cL, 0x97ae418dd274a80aL, 0x2967385d5e3e7682L}, - {0x75c55c7b575a0b79L, 0x3ba4a15774a7ded1L, 0xc3974d73a04fccf3L, 0x705aba4f4a939684L}, - {0x8409a9ea14ebb608L, 0xfad0084e66bac611L, 0x04287254811c1dfbL, 0x086d072b23b30c29L}, - {0xb427c9b367e4756aL, 0xc7537fb902ebc38dL, 0x51de21becd6a205fL, 0x6064ab727923597dL} - }; - -typedef struct { - uint64_t max_width; - blst_fr root_of_unity; - blst_fr *expanded_roots_of_unity; - blst_fr *reverse_roots_of_unity; -} FFTSettings; - -bool is_one(const blst_fr *fr_p); -bool is_power_of_two(uint64_t n); -void fr_from_uint64(blst_fr *a, uint64_t n); -blst_fr *expand_root_of_unity(blst_fr *root_of_unity, uint64_t width); -blst_fr *reverse(blst_fr *r, uint64_t width); -FFTSettings new_fft_settings(unsigned int max_scale); -void free_fft_settings(FFTSettings *s); void slow_ft(blst_fr *out, blst_fr *in, uint64_t offset, uint64_t stride, blst_fr *roots, uint64_t roots_stride, uint64_t l); void fast_ft(blst_fr *out, blst_fr *in, uint64_t offset, uint64_t stride, blst_fr *roots, uint64_t roots_stride, uint64_t l); void fft_helper(blst_fr *out, blst_fr *in, uint64_t offset, uint64_t stride, blst_fr *roots, uint64_t roots_stride, uint64_t l); diff --git a/src/fft_fr_test.c b/src/fft_fr_test.c index 49ffd28..8f30782 100644 --- a/src/fft_fr_test.c +++ b/src/fft_fr_test.c @@ -18,8 +18,6 @@ #include "test_util.h" #include "fft_fr.h" -#define NUM_ROOTS 32 - const uint64_t inv_fft_expected[][4] = { {0x7fffffff80000008L, 0xa9ded2017fff2dffL, 0x199cec0404d0ec02L, 0x39f6d3a994cebea4L}, @@ -40,117 +38,6 @@ const uint64_t inv_fft_expected[][4] = {0xd8cda22e753b3117L, 0x880454ec72238f55L, 0xcaf6199fc14a5353L, 0x197df7c2f05866d4L} }; -void is_one_works(void) { - TEST_CHECK(true == is_one(&one)); -} - -void roots_of_unity_is_the_expected_size(void) { - TEST_CHECK(NUM_ROOTS == - sizeof(scale2_root_of_unity) / sizeof(scale2_root_of_unity[0])); -} - -void roots_of_unity_are_plausible(void) { - blst_fr r; - for (unsigned int i = 0; i < NUM_ROOTS; i++) { - blst_fr_from_uint64(&r, scale2_root_of_unity[i]); - for (unsigned int j = 0; j < i; j++) { - blst_fr_sqr(&r, &r); - } - TEST_CHECK(true == is_one(&r)); - TEST_MSG("Root %d failed", i); - } -} - -void reverse_works(void) { - int n = 24; - blst_fr arr[n + 1]; - blst_fr *rev, diff; - - // Initialise - increasing values - arr[0] = one; - for (int i = 1; i <= n; i++) { - blst_fr_add(arr + i, arr + i - 1, &one); - } - - // Reverse - rev = reverse(arr, n); - - // Verify - decreasing values - for (int i = 0; i < n; i++) { - blst_fr_sub(&diff, rev + i, rev + i + 1); - TEST_CHECK(true == is_one(&diff)); - } - TEST_CHECK(true == is_one(rev + n)); - - free(rev); -} - -void expand_roots_is_plausible(void) { - // Just test one (largeish) value of scale - unsigned int scale = 20; - unsigned int width = 1 << scale; - blst_fr root, *expanded, prod; - - // Initialise - blst_fr_from_uint64(&root, scale2_root_of_unity[scale]); - expanded = expand_root_of_unity(&root, width); - - // Verify - each pair should multiply to one - TEST_CHECK(true == is_one(expanded + 0)); - TEST_CHECK(true == is_one(expanded + width)); - for (unsigned int i = 1; i <= width / 2; i++) { - blst_fr_mul(&prod, expanded + i, expanded + width - i); - TEST_CHECK(true == is_one(&prod)); - } - - free(expanded); -} - -void new_fft_settings_is_plausible(void) { - // Just test one (largeish) value of scale - unsigned int scale = 21; - unsigned int width = 1 << scale; - blst_fr prod; - FFTSettings s = new_fft_settings(scale); - - // Verify - each pair should multiply to one - for (unsigned int i = 1; i <= width; i++) { - blst_fr_mul(&prod, s.expanded_roots_of_unity + i, s.reverse_roots_of_unity + i); - TEST_CHECK(true == is_one(&prod)); - } - - free_fft_settings(&s); -} - -void is_power_of_two_works(void) { - // All actual powers of two - for (int i = 0; i <=63; i++) { - TEST_CHECK(true == is_power_of_two((uint64_t)1 << i)); - TEST_MSG("Case %d", i); - } - - // This is a bit weird - TEST_CHECK(true == is_power_of_two(0)); - - // Not powers of two - TEST_CHECK(false == is_power_of_two(123)); - TEST_CHECK(false == is_power_of_two(1234567)); -} - -void fr_from_uint64_works(void) { - blst_fr a; - fr_from_uint64(&a, 1); - TEST_CHECK(true == is_one(&a)); -} - -void fr_equal_works(void) { - blst_fr a, b; - blst_fr_from_uint64(&a, scale2_root_of_unity[15]); - blst_fr_from_uint64(&b, scale2_root_of_unity[16]); - TEST_CHECK(true == fr_equal(&a, &a)); - TEST_CHECK(false == fr_equal(&a, &b)); -} - void compare_sft_fft(void) { // Initialise: ascending values of i (could be anything), and arbitrary size unsigned int size = 8; @@ -220,15 +107,6 @@ void inverse_fft(void) { TEST_LIST = { - {"is_one_works", is_one_works }, - {"roots_of_unity_is_the_expected_size", roots_of_unity_is_the_expected_size}, - {"roots_of_unity_are_plausible", roots_of_unity_are_plausible}, - {"reverse_works", reverse_works}, - {"expand_roots_is_plausible", expand_roots_is_plausible}, - {"new_fft_settings_is_plausible", new_fft_settings_is_plausible}, - {"is_power_of_two_works", is_power_of_two_works}, - {"fr_from_uint64_works", fr_from_uint64_works}, - {"fr_equal_works", fr_equal_works}, {"compare_sft_fft", compare_sft_fft}, {"roundtrip_fft", roundtrip_fft}, {"inverse_fft", inverse_fft}, diff --git a/src/fft_util.c b/src/fft_util.c new file mode 100644 index 0000000..b926f38 --- /dev/null +++ b/src/fft_util.c @@ -0,0 +1,73 @@ +/* + * 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 "fft_util.h" + +bool is_one(const blst_fr *fr_p) { + uint64_t a[4]; + blst_uint64_from_fr(a, fr_p); + return a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 0; +} + +bool is_power_of_two(uint64_t n) { + return (n & (n - 1)) == 0; +} + +void fr_from_uint64(blst_fr *a, uint64_t n) { + uint64_t vals[] = {n, 0, 0, 0}; + blst_fr_from_uint64(a, vals); +} + +// Returns an array of powers of the root of unity +// Allocates space for the array that needs to be freed later +blst_fr *expand_root_of_unity(blst_fr *root_of_unity, uint64_t width) { + blst_fr *roots = malloc((width + 1) * sizeof(blst_fr)); + roots[0] = one; + roots[1] = *root_of_unity; + + for (int i = 2; !is_one(&roots[i - 1]); i++) { + assert(i <= width); + blst_fr_mul(&roots[i], &roots[i - 1], root_of_unity); + } + assert(is_one(&roots[width])); + + return roots; +} + +// Return a reversed copy of the list of Fr provided +// `width` is one less than the length of `r` +// Allocates space for the array that needs to be freed later +blst_fr *reverse(blst_fr *r, uint64_t width) { + blst_fr *rr = malloc((width + 1) * sizeof(blst_fr)); + for (int i = 0; i <= width; i++) { + rr[i] = r[width - i]; + } + return rr; +} + +FFTSettings new_fft_settings(unsigned int max_scale) { + FFTSettings s; + s.max_width = (uint64_t)1 << max_scale; + blst_fr_from_uint64(&s.root_of_unity, scale2_root_of_unity[max_scale]); + s.expanded_roots_of_unity = expand_root_of_unity(&s.root_of_unity, s.max_width); + s.reverse_roots_of_unity = reverse(s.expanded_roots_of_unity, s.max_width); + return s; +} + +void free_fft_settings(FFTSettings *s) { + free(s->expanded_roots_of_unity); + free(s->reverse_roots_of_unity); +} diff --git a/src/fft_util.h b/src/fft_util.h new file mode 100644 index 0000000..82baed4 --- /dev/null +++ b/src/fft_util.h @@ -0,0 +1,78 @@ +/* + * 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 +#include "c-kzg.h" + +// This is 1 in Blst's `blst_fr` limb representation. Crazy but true. +static const blst_fr one = + {0x00000001fffffffeL, 0x5884b7fa00034802L, 0x998c4fefecbc4ff5L, 0x1824b159acc5056fL}; + +// MODULUS = 52435875175126190479447740508185965837690552500527637822603658699938581184513 +// PRIMITIVE_ROOT = 5 +// [pow(PRIMITIVE_ROOT, (MODULUS - 1) // (2**i), MODULUS) for i in range(32)] +// +// These are not in `blst_fr` limb format and must be converted via `blst_fr_from_uint64()` +static const uint64_t scale2_root_of_unity[][4] = + { + {0x0000000000000001L, 0x0000000000000000L, 0x0000000000000000L, 0x0000000000000000L}, + {0xffffffff00000000L, 0x53bda402fffe5bfeL, 0x3339d80809a1d805L, 0x73eda753299d7d48L}, + {0x0001000000000000L, 0xec03000276030000L, 0x8d51ccce760304d0L, 0x0000000000000000L}, + {0x8dd702cb688bc087L, 0xa032824078eaa4feL, 0xa733b23a98ca5b22L, 0x3f96405d25a31660L}, + {0x95df4b360411fe73L, 0x1ef4e672ebc1e1bbL, 0x6e92a9c460afca4aL, 0x4f2c596e753e4fccL}, + {0x9733f3a6ba60eaa6L, 0xbd7fdf9c77487ae7L, 0xd84f8612c8b6cc00L, 0x476fa2fb6162ffabL}, + {0xd2fc5e69ac5db47fL, 0xa12a70a615d0b8e4L, 0x293b1d67bc8de5d9L, 0x0e4840ac57f86f5eL}, + {0xb750da4cab28e208L, 0x501dff643be95635L, 0x8cbe2437f0b4b276L, 0x07d0c802a94a946eL}, + {0x2cabadec2fe322b8L, 0x752c84f315412560L, 0x32a732ae1a3b0aefL, 0x2e95da59a33dcbf2L}, + {0x33811ea1fe0c65f4L, 0x15c1ad4c687f28a2L, 0xecfbede342dee7f4L, 0x1bb466679a5d88b1L}, + {0xd58a5af42d010ff9L, 0x79efd6b0570bf109L, 0x3ed6d55a6350721dL, 0x2f27b09858f43cefL}, + {0x74a1f6718c130477L, 0xa534af14b61e0abeL, 0xeb674a1a620890d7L, 0x43527a8bca252472L}, + {0x450d9f977ea8ee05L, 0x565af17137d56fc0L, 0xe155cb4893f9e9acL, 0x110cebd0c8e9101bL}, + {0x23c9159959a0be92L, 0x87d188ce7a027759L, 0x70491431cab3c3ccL, 0x0ac00eb8b3f7f8daL}, + {0x13e96ade69583404L, 0x82c057275306243dL, 0x77e48bf529ca9f2aL, 0x50646ac81fe19595L}, + {0xe6a354dda97eccd4L, 0x39929d2e88fbbc57L, 0xa22ba63dd6e7b1c8L, 0x42c22911f5f07f43L}, + {0x137b458acfc35f7aL, 0x0caba63a29c01b06L, 0x0409ee987a02402cL, 0x6709c6cd56aa725bL}, + {0x10251f7d8831e03eL, 0x77d85a937ff858ecL, 0xebe905bd4fb9ac5cL, 0x05deb333f8727901L}, + {0xbf87b689b9009408L, 0x4f730e7ddd3ccc96L, 0xfd7f05ba4610300cL, 0x5ef5e8db0b8ac903L}, + {0x6499688417cd0c14L, 0xa672867368812f7fL, 0x2e1d9a1922cc3253L, 0x3a689e83aa0a1d80L}, + {0x20b53cbe41144deaL, 0x870c46fac2f0fcbdL, 0x556c35f6537d6971L, 0x3436287f5f686d91L}, + {0x007e082a436ba2e7L, 0x67c6630f9116e877L, 0x36f8f165fb4460f7L, 0x6eee34d57e7046e0L}, + {0xc5b670eea53a56d1L, 0x127d1f4253037d7bL, 0x57d4257ea722c2e2L, 0x03ae26a333cbd838L}, + {0x1e91484876504cf8L, 0x55bbbf1eb63edd02L, 0xbcdafec84e55aa02L, 0x5145c4cd2dc0beb0L}, + {0x5b90153a1ab70e2cL, 0x8deffa3175fb0ab8L, 0xc553ae2346900c95L, 0x1d31dcdc6bd3118cL}, + {0x801c894c59a2e8ebL, 0xbc535c5ce12fc974L, 0x95508d2747d39803L, 0x16d9d3cdac5d094fL}, + {0x810fa372cca1d8beL, 0xc67b8c2882e0bfa7L, 0xdbb4edf0e2d35bc2L, 0x712d15805087c995L}, + {0xeb162203fd88f133L, 0xac96c38ff010ea74L, 0x4307987fe64cfc70L, 0x350fe98d37b7a114L}, + {0xaba2f51842f2a254L, 0x4d7f3c3aa71efc0cL, 0x97ae418dd274a80aL, 0x2967385d5e3e7682L}, + {0x75c55c7b575a0b79L, 0x3ba4a15774a7ded1L, 0xc3974d73a04fccf3L, 0x705aba4f4a939684L}, + {0x8409a9ea14ebb608L, 0xfad0084e66bac611L, 0x04287254811c1dfbL, 0x086d072b23b30c29L}, + {0xb427c9b367e4756aL, 0xc7537fb902ebc38dL, 0x51de21becd6a205fL, 0x6064ab727923597dL} + }; + +typedef struct { + uint64_t max_width; + blst_fr root_of_unity; + blst_fr *expanded_roots_of_unity; + blst_fr *reverse_roots_of_unity; +} FFTSettings; + +bool is_one(const blst_fr *fr_p); +bool is_power_of_two(uint64_t n); +void fr_from_uint64(blst_fr *a, uint64_t n); +blst_fr *expand_root_of_unity(blst_fr *root_of_unity, uint64_t width); +blst_fr *reverse(blst_fr *r, uint64_t width); +FFTSettings new_fft_settings(unsigned int max_scale); +void free_fft_settings(FFTSettings *s); diff --git a/src/fft_util_test.c b/src/fft_util_test.c new file mode 100644 index 0000000..684c0ff --- /dev/null +++ b/src/fft_util_test.c @@ -0,0 +1,146 @@ +/* + * 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 "test_util.h" +#include "fft_util.h" + +#define NUM_ROOTS 32 + +void is_one_works(void) { + TEST_CHECK(true == is_one(&one)); +} + +void roots_of_unity_is_the_expected_size(void) { + TEST_CHECK(NUM_ROOTS == + sizeof(scale2_root_of_unity) / sizeof(scale2_root_of_unity[0])); +} + +void roots_of_unity_are_plausible(void) { + blst_fr r; + for (unsigned int i = 0; i < NUM_ROOTS; i++) { + blst_fr_from_uint64(&r, scale2_root_of_unity[i]); + for (unsigned int j = 0; j < i; j++) { + blst_fr_sqr(&r, &r); + } + TEST_CHECK(true == is_one(&r)); + TEST_MSG("Root %d failed", i); + } +} + +void reverse_works(void) { + int n = 24; + blst_fr arr[n + 1]; + blst_fr *rev, diff; + + // Initialise - increasing values + arr[0] = one; + for (int i = 1; i <= n; i++) { + blst_fr_add(arr + i, arr + i - 1, &one); + } + + // Reverse + rev = reverse(arr, n); + + // Verify - decreasing values + for (int i = 0; i < n; i++) { + blst_fr_sub(&diff, rev + i, rev + i + 1); + TEST_CHECK(true == is_one(&diff)); + } + TEST_CHECK(true == is_one(rev + n)); + + free(rev); +} + +void expand_roots_is_plausible(void) { + // Just test one (largeish) value of scale + unsigned int scale = 20; + unsigned int width = 1 << scale; + blst_fr root, *expanded, prod; + + // Initialise + blst_fr_from_uint64(&root, scale2_root_of_unity[scale]); + expanded = expand_root_of_unity(&root, width); + + // Verify - each pair should multiply to one + TEST_CHECK(true == is_one(expanded + 0)); + TEST_CHECK(true == is_one(expanded + width)); + for (unsigned int i = 1; i <= width / 2; i++) { + blst_fr_mul(&prod, expanded + i, expanded + width - i); + TEST_CHECK(true == is_one(&prod)); + } + + free(expanded); +} + +void new_fft_settings_is_plausible(void) { + // Just test one (largeish) value of scale + unsigned int scale = 21; + unsigned int width = 1 << scale; + blst_fr prod; + FFTSettings s = new_fft_settings(scale); + + // Verify - each pair should multiply to one + for (unsigned int i = 1; i <= width; i++) { + blst_fr_mul(&prod, s.expanded_roots_of_unity + i, s.reverse_roots_of_unity + i); + TEST_CHECK(true == is_one(&prod)); + } + + free_fft_settings(&s); +} + +void is_power_of_two_works(void) { + // All actual powers of two + for (int i = 0; i <=63; i++) { + TEST_CHECK(true == is_power_of_two((uint64_t)1 << i)); + TEST_MSG("Case %d", i); + } + + // This is a bit weird + TEST_CHECK(true == is_power_of_two(0)); + + // Not powers of two + TEST_CHECK(false == is_power_of_two(123)); + TEST_CHECK(false == is_power_of_two(1234567)); +} + +void fr_from_uint64_works(void) { + blst_fr a; + fr_from_uint64(&a, 1); + TEST_CHECK(true == is_one(&a)); +} + +void fr_equal_works(void) { + blst_fr a, b; + blst_fr_from_uint64(&a, scale2_root_of_unity[15]); + blst_fr_from_uint64(&b, scale2_root_of_unity[16]); + TEST_CHECK(true == fr_equal(&a, &a)); + TEST_CHECK(false == fr_equal(&a, &b)); +} + +TEST_LIST = + { + {"is_one_works", is_one_works }, + {"roots_of_unity_is_the_expected_size", roots_of_unity_is_the_expected_size}, + {"roots_of_unity_are_plausible", roots_of_unity_are_plausible}, + {"reverse_works", reverse_works}, + {"expand_roots_is_plausible", expand_roots_is_plausible}, + {"new_fft_settings_is_plausible", new_fft_settings_is_plausible}, + {"is_power_of_two_works", is_power_of_two_works}, + {"fr_from_uint64_works", fr_from_uint64_works}, + {"fr_equal_works", fr_equal_works}, + { NULL, NULL } /* zero record marks the end of the list */ + };