diff --git a/.github/workflows/publish-results.yml b/.github/workflows/publish-results.yml new file mode 100644 index 0000000..4b7b9f3 --- /dev/null +++ b/.github/workflows/publish-results.yml @@ -0,0 +1,13 @@ +# Needed to publish test results in fork +name: Testing Callback + +on: + workflow_run: + workflows: ["PR Testing"] + types: + - completed + +jobs: + call-reusable-workflow: + name: Call Reusable Testing Callback Workflow + uses: NilFoundation/ci-cd/.github/workflows/reusable-crypto3-publish-result.yml@v1 diff --git a/.github/workflows/pull-request.yml b/.github/workflows/pull-request.yml new file mode 100644 index 0000000..0639277 --- /dev/null +++ b/.github/workflows/pull-request.yml @@ -0,0 +1,24 @@ +name: PR Testing + +on: + pull_request: + types: + - opened + - synchronize + +jobs: + handle-syncwith: + name: Call Reusable SyncWith Handler + uses: NilFoundation/ci-cd/.github/workflows/reusable-handle-syncwith.yml@v1.1.2 + with: + ci-cd-ref: 'v1.1.2' + secrets: inherit + + matrix-test: + name: Call Reusable Crypto3 Testing + needs: + - handle-syncwith + uses: NilFoundation/ci-cd/.github/workflows/reusable-crypto3-testing.yml@v1.1.2 + with: + submodules-refs: ${{ needs.handle-syncwith.outputs.prs-refs }} + secrets: inherit diff --git a/include/nil/crypto3/detail/inline_variable.hpp b/include/nil/crypto3/detail/inline_variable.hpp deleted file mode 100644 index 1c36127..0000000 --- a/include/nil/crypto3/detail/inline_variable.hpp +++ /dev/null @@ -1,60 +0,0 @@ -//---------------------------------------------------------------------------// -// Copyright (c) 2018-2020 Mikhail Komarov -// -// MIT License -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -//---------------------------------------------------------------------------// - -#ifndef CRYPTO3_INLINE_VARIABLE_HPP -#define CRYPTO3_INLINE_VARIABLE_HPP - -#define CRYPTO3_CXX_STD_14 201402L -#define CRYPTO3_CXX_STD_17 201703L - -#if defined(_MSVC_LANG) && _MSVC_LANG > __cplusplus // Older clangs define _MSVC_LANG < __cplusplus -#define CRYPTO3_CXX_VER _MSVC_LANG -#else -#define CRYPTO3_CXX_VER __cplusplus -#endif - -#ifndef CRYPTO3_CXX17_INLINE_VARIABLES -#ifdef __cpp_inline_variables -#define CRYPTO3_CXX17_INLINE_VARIABLES __cpp_inline_variables -#else -#define CRYPTO3_CXX17_INLINE_VARIABLES (CRYPTO3_CXX_VER >= CRYPTO3_CXX_STD_17) -#endif -#endif - -#ifdef CRYPTO3_CXX17_INLINE_VARIABLES -#define CRYPTO3_INLINE_VARIABLE(TYPE, NAME, VALUE) \ - constexpr static inline TYPE NAME() { \ - return TYPE VALUE; \ - } -#else -#define CRYPTO3_INLINE_VARIABLE(TYPE, NAME, VALUE) \ - struct NAME { \ - inline TYPE const &operator()() const { \ - static TYPE const v VALUE; \ - return v; \ - } \ - }; -#endif - -#endif // CRYPTO3_INLINE_VARIABLE_HPP diff --git a/include/nil/crypto3/detail/make_array.hpp b/include/nil/crypto3/detail/make_array.hpp deleted file mode 100644 index cd01123..0000000 --- a/include/nil/crypto3/detail/make_array.hpp +++ /dev/null @@ -1,75 +0,0 @@ -//---------------------------------------------------------------------------// -// Copyright (c) 2018-2020 Mikhail Komarov -// -// MIT License -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -//---------------------------------------------------------------------------// - -#ifndef CRYPTO3_MAKE_ARRAY_HPP -#define CRYPTO3_MAKE_ARRAY_HPP - -#include -#include - -namespace nil { - namespace crypto3 { - namespace detail { - template - struct indices { - using next = indices; - }; - template - struct build_indices { - using type = typename build_indices::type::next; - }; - template<> - struct build_indices<0> { - using type = indices<>; - }; - template - using BuildIndices = typename build_indices::type; - - template - using ValueType = typename std::iterator_traits::value_type; - - // internal overload with indices tag - - template, sizeof...(I)>> - - Array make_array(InputIterator first, indices) { - return Array {{(void(I), *first++)...}}; - } - } // namespace detail - - // externally visible interface - template - std::array, N> make_array(RandomAccessIterator first, - RandomAccessIterator last) { - // last is not relevant if we're assuming the size is N - // I'll assert it is correct anyway - assert(last - first == N); - return make_array(first, detail::BuildIndices {}); - } - } // namespace crypto3 -} // namespace nil - -#endif // CRYPTO3_MAKE_ARRAY_HPP diff --git a/include/nil/crypto3/detail/make_uint_t.hpp b/include/nil/crypto3/detail/make_uint_t.hpp deleted file mode 100644 index 8bebc31..0000000 --- a/include/nil/crypto3/detail/make_uint_t.hpp +++ /dev/null @@ -1,63 +0,0 @@ -//---------------------------------------------------------------------------// -// Copyright (c) 2018-2020 Mikhail Komarov -// -// MIT License -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -//---------------------------------------------------------------------------/// - -#ifndef CRYPTO3_MAKE_UINT_T_HPP -#define CRYPTO3_MAKE_UINT_T_HPP - -#include - -#include - -namespace nil { - namespace crypto3 { - namespace detail { - template - static inline typename boost::uint_t::exact extract_uint_t(Integer v, std::size_t position) { - return static_cast::exact>(v >> - (((~position) & (sizeof(Integer) - 1)) << 3)); - } - - template - static inline typename boost::uint_t::exact make_uint_t(const std::initializer_list &args) { - typedef typename std::initializer_list::value_type value_type; - typename boost::uint_t::exact result = 0; - - - for (const value_type &itr : args) { - result = static_cast::exact>( - (result << std::numeric_limits::digits) | itr); - } - - return result; - } - - template - static inline typename boost::uint_t::exact make_uint_t(Args... args) { - return make_uint_t>::type>({args...}); - } - } // namespace detail - } // namespace crypto3 -} // namespace nil - -#endif // CRYPTO3_MAKE_UINT_T_HPP diff --git a/include/nil/crypto3/detail/state_adder.hpp b/include/nil/crypto3/detail/state_adder.hpp deleted file mode 100644 index bb69c58..0000000 --- a/include/nil/crypto3/detail/state_adder.hpp +++ /dev/null @@ -1,46 +0,0 @@ -//---------------------------------------------------------------------------// -// Copyright (c) 2018-2020 Mikhail Komarov -// -// MIT License -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -//---------------------------------------------------------------------------// - -#ifndef CRYPTO3_DETAIL_STATE_ADDER_HPP -#define CRYPTO3_DETAIL_STATE_ADDER_HPP - -namespace nil { - namespace crypto3 { - namespace detail { - struct state_adder { - template - void operator()(T &s1, T const &s2) { - typedef typename T::size_type size_type; - size_type n = (s2.size() < s1.size() ? s2.size() : s1.size()); - for (typename T::size_type i = 0; i < n; ++i) { - s1[i] += s2[i]; - } - } - }; - - } // namespace detail - } // namespace crypto3 -} // namespace nil - -#endif // CRYPTO3_BLOCK_DETAIL_STATE_ADDER_HPP diff --git a/include/nil/crypto3/detail/unbounded_shift.hpp b/include/nil/crypto3/detail/unbounded_shift.hpp deleted file mode 100644 index 832cec4..0000000 --- a/include/nil/crypto3/detail/unbounded_shift.hpp +++ /dev/null @@ -1,76 +0,0 @@ -//---------------------------------------------------------------------------// -// Copyright (c) 2018-2020 Mikhail Komarov -// -// MIT License -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -//---------------------------------------------------------------------------// - -#ifndef CRYPTO3_DETAIL_UNBOUNDED_SHIFT_HPP -#define CRYPTO3_DETAIL_UNBOUNDED_SHIFT_HPP - -#include - -namespace nil { - namespace crypto3 { - namespace detail { - - template - struct unbounded_shifter { - static T shl(T x) { - return unbounded_shifter::shl(T(x << 1)); - } - - static T shr(T x) { - return unbounded_shifter::shr(T(x >> 1)); - } - }; - - template - struct unbounded_shifter<0, T> { - static T shl(T x) { - return x; - } - - static T shr(T x) { - return x; - } - }; - - template - T unbounded_shl(T x) { - return unbounded_shifter::shl(x); - } - - template - T unbounded_shr(T x) { - return unbounded_shifter::shr(x); - } - - template - T low_bits(T x) { - T highmask = unbounded_shl(~T()); - return T(x & ~highmask); - } - - } // namespace detail - } // namespace crypto3 -} // namespace nil - -#endif // CRYPTO3_BLOCK_DETAIL_UNBOUNDED_SHIFT_HPP diff --git a/include/nil/crypto3/math/algorithms/make_evaluation_domain.hpp b/include/nil/crypto3/math/algorithms/make_evaluation_domain.hpp index 491860c..44c16ee 100644 --- a/include/nil/crypto3/math/algorithms/make_evaluation_domain.hpp +++ b/include/nil/crypto3/math/algorithms/make_evaluation_domain.hpp @@ -49,6 +49,7 @@ namespace nil { */ template std::shared_ptr> make_evaluation_domain(std::size_t m) { + typedef std::shared_ptr> result_type; const std::size_t big = 1ul << (std::size_t(std::ceil(std::log2(m))) - 1); diff --git a/include/nil/crypto3/math/detail/field_utils.hpp b/include/nil/crypto3/math/detail/field_utils.hpp index a217950..102694b 100644 --- a/include/nil/crypto3/math/detail/field_utils.hpp +++ b/include/nil/crypto3/math/detail/field_utils.hpp @@ -40,7 +40,7 @@ namespace nil { using namespace nil::crypto3::algebra; - std::size_t bitreverse(std::size_t n, const std::size_t l) { + inline std::size_t bitreverse(std::size_t n, const std::size_t l) { std::size_t r = 0; for (std::size_t k = 0; k < l; ++k) { r = (r << 1) | (n & 1); diff --git a/include/nil/crypto3/math/domains/arithmetic_sequence_domain.hpp b/include/nil/crypto3/math/domains/arithmetic_sequence_domain.hpp index 56161af..71c253e 100755 --- a/include/nil/crypto3/math/domains/arithmetic_sequence_domain.hpp +++ b/include/nil/crypto3/math/domains/arithmetic_sequence_domain.hpp @@ -82,7 +82,7 @@ namespace nil { precomputation_sentinel = false; } - void fft(std::vector &a) { + void fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -116,7 +116,7 @@ namespace nil { } } - void inverse_fft(std::vector &a) { + void inverse_fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -151,7 +151,7 @@ namespace nil { newton_to_monomial_basis(a, subproduct_tree, this->m); } - std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) { + std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) override { /* Compute Lagrange polynomial of size m, with m+1 points (x_0, y_0), ... ,(x_m, y_m) */ /* Evaluate for x = t */ /* Return coeffs for each l_j(x) = (l / l_i[j]) * w[j] */ @@ -202,11 +202,11 @@ namespace nil { } std::vector evaluate_all_lagrange_polynomials(const typename std::vector::const_iterator &t_powers_begin, - const typename std::vector::const_iterator &t_powers_end) { - if(std::distance(t_powers_begin, t_powers_end) < this->m) { + const typename std::vector::const_iterator &t_powers_end) override { + if(std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) { throw std::invalid_argument("arithmetic_sequence_radix2: expected std::distance(t_powers_begin, t_powers_end) >= this->m"); } - + /* Compute Lagrange polynomial of size m, with m+1 points (x_0, y_0), ... ,(x_m, y_m) */ /* Evaluate for x = t */ /* Return coeffs for each l_j(x) = (l / l_i[j]) * w[j] */ @@ -263,20 +263,26 @@ namespace nil { for(std::size_t j = 0; j < l[i].size(); ++j) { result[i] = result[i] + t_powers_begin[j] * l[i][j]; - } + } result[i] = result[i] * w[i]; } return result; } - field_value_type get_domain_element(const std::size_t idx) { + // This one is not the unity root actually, but it's ok for our purposes. + const field_value_type& get_unity_root() override { + return arithmetic_generator; + } + + field_value_type get_domain_element(const std::size_t idx) override { if (!this->precomputation_sentinel) do_precomputation(); return this->arithmetic_sequence[idx]; } - field_value_type compute_vanishing_polynomial(const field_value_type &t) { + + field_value_type compute_vanishing_polynomial(const field_value_type &t) override { if (!this->precomputation_sentinel) do_precomputation(); @@ -287,7 +293,8 @@ namespace nil { } return Z; } - polynomial get_vanishing_polynomial() { + + polynomial get_vanishing_polynomial() override { if (!precomputation_sentinel) do_precomputation(); @@ -297,7 +304,8 @@ namespace nil { } return z; } - void add_poly_z(const field_value_type &coeff, std::vector &H) { + + void add_poly_z(const field_value_type &coeff, std::vector &H) override { if (H.size() != this->m + 1) throw std::invalid_argument("arithmetic: expected H.size() == this->m+1"); @@ -321,7 +329,8 @@ namespace nil { H[i] += (x[i] * coeff); } } - void divide_by_z_on_coset(std::vector &P) { + + void divide_by_z_on_coset(std::vector &P) override { const field_value_type coset = this->arithmetic_generator; /* coset in arithmetic sequence? */ const field_value_type Z_inverse_at_coset = this->compute_vanishing_polynomial(coset).inversed(); for (std::size_t i = 0; i < this->m; ++i) { diff --git a/include/nil/crypto3/math/domains/basic_radix2_domain.hpp b/include/nil/crypto3/math/domains/basic_radix2_domain.hpp index 916512a..07158b0 100755 --- a/include/nil/crypto3/math/domains/basic_radix2_domain.hpp +++ b/include/nil/crypto3/math/domains/basic_radix2_domain.hpp @@ -48,13 +48,22 @@ namespace nil { class basic_radix2_domain : public evaluation_domain { typedef typename FieldType::value_type field_value_type; typedef ValueType value_type; + typedef std::pair, std::vector> cache_type; + std::unique_ptr fft_cache; + void create_fft_cache() { + fft_cache = std::make_unique(std::vector(), std::vector()); + detail::create_fft_cache(this->m, omega, fft_cache->first); + detail::create_fft_cache(this->m, omega.inversed(), fft_cache->second); + } public: typedef FieldType field_type; - field_value_type omega; + const field_value_type omega; - basic_radix2_domain(const std::size_t m) : evaluation_domain(m) { + basic_radix2_domain(const std::size_t m) + : evaluation_domain(m), + omega(unity_root(m)) { if (m <= 1) throw std::invalid_argument("basic_radix2(): expected m > 1"); @@ -64,11 +73,9 @@ namespace nil { throw std::invalid_argument( "basic_radix2(): expected logm <= fields::arithmetic_params::s"); } - - omega = unity_root(m); } - void fft(std::vector &a) { + void fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -77,10 +84,13 @@ namespace nil { } } - detail::basic_radix2_fft(a, omega); + if (fft_cache == nullptr) { + create_fft_cache(); + } + detail::basic_radix2_fft_cached(a, fft_cache->first); } - void inverse_fft(std::vector &a) { + void inverse_fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -89,7 +99,10 @@ namespace nil { } } - detail::basic_radix2_fft(a, omega.inversed()); + if (fft_cache == nullptr) { + create_fft_cache(); + } + detail::basic_radix2_fft_cached(a, fft_cache->second); const field_value_type sconst = field_value_type(a.size()).inversed(); for (std::size_t i = 0; i < a.size(); ++i) { @@ -97,13 +110,13 @@ namespace nil { } } - std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) { + std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) override { return detail::basic_radix2_evaluate_all_lagrange_polynomials(this->m, t); } std::vector evaluate_all_lagrange_polynomials(const typename std::vector::const_iterator &t_powers_begin, - const typename std::vector::const_iterator &t_powers_end) { - if (std::distance(t_powers_begin, t_powers_end) < this->m) { + const typename std::vector::const_iterator &t_powers_end) override { + if (std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) { throw std::invalid_argument("basic_radix2: expected std::distance(t_powers_begin, t_powers_end) >= this->m"); } std::vector tmp(t_powers_begin, t_powers_begin + this->m); @@ -111,22 +124,26 @@ namespace nil { return tmp; } - field_value_type get_domain_element(const std::size_t idx) { + const field_value_type& get_unity_root() override { + return omega; + } + + field_value_type get_domain_element(const std::size_t idx) override { return omega.pow(idx); } - field_value_type compute_vanishing_polynomial(const field_value_type &t) { + field_value_type compute_vanishing_polynomial(const field_value_type &t) override { return (t.pow(this->m)) - field_value_type::one(); } - polynomial get_vanishing_polynomial() { + polynomial get_vanishing_polynomial() override { polynomial z(this->m + 1, field_value_type::zero()); z[this->m] = field_value_type::one(); z[0] = -field_value_type::one(); return z; } - void add_poly_z(const field_value_type &coeff, std::vector &H) { + void add_poly_z(const field_value_type &coeff, std::vector &H) override { if (H.size() != this->m + 1) throw std::invalid_argument("basic_radix2: expected H.size() == this->m+1"); @@ -134,7 +151,7 @@ namespace nil { H[0] -= coeff; } - void divide_by_z_on_coset(std::vector &P) { + void divide_by_z_on_coset(std::vector &P) override { const field_value_type coset = fields::arithmetic_params::multiplicative_generator; const field_value_type Z_inverse_at_coset = this->compute_vanishing_polynomial(coset).inversed(); for (std::size_t i = 0; i < this->m; ++i) { diff --git a/include/nil/crypto3/math/domains/detail/basic_radix2_domain_aux.hpp b/include/nil/crypto3/math/domains/detail/basic_radix2_domain_aux.hpp index 3cee717..85dad4f 100644 --- a/include/nil/crypto3/math/domains/detail/basic_radix2_domain_aux.hpp +++ b/include/nil/crypto3/math/domains/detail/basic_radix2_domain_aux.hpp @@ -1,6 +1,7 @@ //---------------------------------------------------------------------------// // Copyright (c) 2020-2021 Mikhail Komarov // Copyright (c) 2020-2021 Nikita Kaskov +// Copyright (c) 2024 Dmitrii Tabalin // // MIT License // @@ -27,6 +28,7 @@ #define CRYPTO3_MATH_BASIC_RADIX2_DOMAIN_AUX_HPP #include +#include #include #include @@ -39,17 +41,32 @@ namespace nil { namespace math { namespace detail { + /* + * Building caches for fft operations + */ + template + void create_fft_cache( + const std::size_t size, + const typename FieldType::value_type &omega, + std::vector &cache) { + typedef typename FieldType::value_type value_type; + cache.resize(size); + cache[0] = value_type::one(); + for (std::size_t i = 1; i < size; ++i) { + cache[i] = cache[i - 1] * omega; + } + } + /* * Below we make use of pseudocode from [CLRS 2n Ed, pp. 864]. * Also, note that it's the caller's responsibility to multiply by 1/N. */ template - void basic_radix2_fft(Range &a, const typename FieldType::value_type &omega) { + void basic_radix2_fft_cached(Range &a, const std::vector &omega_cache) { typedef typename std::iterator_traits()))>::value_type value_type; - typedef typename FieldType::value_type field_value_type; BOOST_STATIC_ASSERT(algebra::is_field::value); - + // It now supports curve elements too, should probably some other assertion about the field type and value type // BOOST_STATIC_ASSERT(std::is_same::value); @@ -64,23 +81,36 @@ namespace nil { std::swap(a[k], a[rk]); } - std::size_t m = 1; // invariant: m = 2^{s-1} - for (std::size_t s = 1; s <= logn; ++s) { + // invariant: m = 2^{s-1} + value_type t; + for (std::size_t s = 1, m = 1, inc = n / 2; s <= logn; ++s, m <<= 1, inc >>= 1) { // w_m is 2^s-th root of unity now - const field_value_type w_m = omega.pow(n / (2 * m)); - - asm volatile("/* pre-inner */"); for (std::size_t k = 0; k < n; k += 2 * m) { - field_value_type w = field_value_type::one(); - for (std::size_t j = 0; j < m; ++j) { - const value_type t = a[k + j + m] * w; - a[k + j + m] = a[k + j] - t; - a[k + j] = a[k + j] + t; - w *= w_m; + for (std::size_t j = 0, idx = 0; j < m; ++j, idx += inc) { + t = a[k + j + m]; + t *= omega_cache[idx].data; + a[k + j + m] = a[k + j]; + a[k + j + m] -= t; + a[k + j] += t; } } - asm volatile("/* post-inner */"); - m *= 2; + } + } + + /** + * Note that it's the caller's responsibility to multiply by 1/N. + */ + template + void basic_radix2_fft( + Range &a, const typename FieldType::value_type &omega, + std::shared_ptr> omega_cache = nullptr) { + + if (omega_cache == nullptr) { + std::vector omega_powers; + create_fft_cache(a.size(), omega, omega_powers); + basic_radix2_fft_cached(a, omega_powers); + } else { + basic_radix2_fft_cached(a, *omega_cache); } } diff --git a/include/nil/crypto3/math/domains/evaluation_domain.hpp b/include/nil/crypto3/math/domains/evaluation_domain.hpp index 78a18bf..6851afd 100644 --- a/include/nil/crypto3/math/domains/evaluation_domain.hpp +++ b/include/nil/crypto3/math/domains/evaluation_domain.hpp @@ -47,16 +47,8 @@ namespace nil { public: typedef FieldType field_type; - field_value_type root; - field_value_type root_inverse; - field_value_type domain; - field_value_type domain_inverse; - field_value_type generator; - field_value_type generator_inverse; - - std::size_t m; - std::size_t log2_size; - std::size_t generator_size; + const std::size_t m; + const std::size_t log2_size; /** * Construct an evaluation domain S of size m, if possible. @@ -69,6 +61,16 @@ namespace nil { return m; } + /* + * Virtual destructor. + */ + virtual ~evaluation_domain() {}; + + /** + * Get the unity root. + */ + virtual const field_value_type& get_unity_root() = 0; + /** * Get the idx-th element in S. */ @@ -104,8 +106,9 @@ namespace nil { * The output is a vector (b_{0},...,b_{m-1}) * where b_{i} is the evaluation of L_{i,S}(z) at z = t. */ - virtual std::vector evaluate_all_lagrange_polynomials(const typename std::vector::const_iterator &t_powers_begin, - const typename std::vector::const_iterator &t_powers_end) = 0; + virtual std::vector evaluate_all_lagrange_polynomials( + const typename std::vector::const_iterator &t_powers_begin, + const typename std::vector::const_iterator &t_powers_end) = 0; /** * Evaluate the vanishing polynomial of S at the field element t. @@ -128,10 +131,7 @@ namespace nil { virtual void divide_by_z_on_coset(std::vector &P) = 0; bool operator==(const evaluation_domain &rhs) const { - return root == rhs.root && root_inverse == rhs.root_inverse && domain == rhs.domain && - domain_inverse == rhs.domain_inverse && generator == rhs.generator && - generator_inverse == rhs.generator_inverse && m == rhs.m && log2_size == rhs.log2_size && - generator_size == rhs.generator_size; + return m == rhs.m && log2_size == rhs.log2_size; } }; } // namespace math diff --git a/include/nil/crypto3/math/domains/extended_radix2_domain.hpp b/include/nil/crypto3/math/domains/extended_radix2_domain.hpp index 172c977..57d4ca9 100755 --- a/include/nil/crypto3/math/domains/extended_radix2_domain.hpp +++ b/include/nil/crypto3/math/domains/extended_radix2_domain.hpp @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -46,15 +47,28 @@ namespace nil { class extended_radix2_domain : public evaluation_domain { typedef typename FieldType::value_type field_value_type; typedef ValueType value_type; + typedef std::pair, std::vector> cache_type; + std::unique_ptr fft_cache; + + void create_fft_cache() { + fft_cache = std::make_unique(std::vector(), + std::vector()); + detail::create_fft_cache(small_m, omega, fft_cache->first); + detail::create_fft_cache(small_m, omega.inversed(), fft_cache->second); + } public: typedef FieldType field_type; - std::size_t small_m; - field_value_type omega; - field_value_type shift; + const std::size_t small_m; + const field_value_type omega; + const field_value_type shift; - extended_radix2_domain(const std::size_t m) : evaluation_domain(m) { + extended_radix2_domain(const std::size_t m) + : evaluation_domain(m), + small_m(m / 2), + omega(unity_root(small_m)), + shift(detail::coset_shift()) { if (m <= 1) throw std::invalid_argument("extended_radix2(): expected m > 1"); @@ -64,15 +78,9 @@ namespace nil { throw std::invalid_argument( "extended_radix2(): expected logm == fields::arithmetic_params::s + 1"); } - - small_m = m / 2; - - omega = unity_root(small_m); - - shift = detail::coset_shift(); } - void fft(std::vector &a) { + void fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -94,8 +102,11 @@ namespace nil { shift_i *= shift; } - detail::basic_radix2_fft(a0, omega); - detail::basic_radix2_fft(a1, omega); + if (fft_cache == nullptr) { + create_fft_cache(); + } + detail::basic_radix2_fft_cached(a0, fft_cache->first); + detail::basic_radix2_fft_cached(a1, fft_cache->first); for (std::size_t i = 0; i < small_m; ++i) { a[i] = a0[i]; @@ -103,7 +114,7 @@ namespace nil { } } - void inverse_fft(std::vector &a) { + void inverse_fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -116,9 +127,11 @@ namespace nil { std::vector a0(a.begin(), a.begin() + small_m); std::vector a1(a.begin() + small_m, a.end()); - const field_value_type omega_inverse = omega.inversed(); - detail::basic_radix2_fft(a0, omega_inverse); - detail::basic_radix2_fft(a1, omega_inverse); + if (fft_cache == nullptr) { + create_fft_cache(); + } + detail::basic_radix2_fft_cached(a0, fft_cache->second); + detail::basic_radix2_fft_cached(a1, fft_cache->second); const field_value_type shift_to_small_m = shift.pow(small_m); const field_value_type sconst = (field_value_type(small_m) * (field_value_type::one() - shift_to_small_m)).inversed(); @@ -134,7 +147,7 @@ namespace nil { } } - std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) { + std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) override { const std::vector T0 = detail::basic_radix2_evaluate_all_lagrange_polynomials(small_m, t); const std::vector T1 = @@ -157,18 +170,18 @@ namespace nil { } std::vector evaluate_all_lagrange_polynomials(const typename std::vector::const_iterator &t_powers_begin, - const typename std::vector::const_iterator &t_powers_end) { - if(std::distance(t_powers_begin, t_powers_end) < this->m) { + const typename std::vector::const_iterator &t_powers_end) override { + if(std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) { throw std::invalid_argument("extended_radix2: expected std::distance(t_powers_begin, t_powers_end) >= this->m"); } - + basic_radix2_domain basic_domain(small_m); - + std::vector T0 = basic_domain.evaluate_all_lagrange_polynomials(t_powers_begin, t_powers_end); std::vector T0_times_t_to_small_m = basic_domain.evaluate_all_lagrange_polynomials(t_powers_begin + small_m, t_powers_end); - + field_value_type shift_inverse = shift.inversed(); std::vector shift_inv_t_powers(small_m); std::vector shift_inv_t_powers_times_t_to_small_m(small_m); @@ -198,7 +211,11 @@ namespace nil { return result; } - field_value_type get_domain_element(const std::size_t idx) { + const field_value_type& get_unity_root() override { + return omega; + } + + field_value_type get_domain_element(const std::size_t idx) override { if (idx < small_m) { return omega.pow(idx); } else { @@ -206,11 +223,11 @@ namespace nil { } } - field_value_type compute_vanishing_polynomial(const field_value_type &t) { + field_value_type compute_vanishing_polynomial(const field_value_type &t) override { return (t.pow(small_m) - field_value_type::one()) * (t.pow(small_m) - shift.pow(small_m)); } - polynomial get_vanishing_polynomial() { + polynomial get_vanishing_polynomial() override { polynomial z(2*small_m + 1, field_value_type::zero()); field_value_type shift_to_small_m = shift.pow(small_m); z[2*small_m] = field_value_type::one(); @@ -219,7 +236,7 @@ namespace nil { return z; } - void add_poly_z(const field_value_type &coeff, std::vector &H) { + void add_poly_z(const field_value_type &coeff, std::vector &H) override { // if (H.size() != this->m + 1) // throw std::invalid_argument("extended_radix2: expected H.size() == this->m+1"); @@ -230,7 +247,7 @@ namespace nil { H[0] += coeff * shift_to_small_m; } - void divide_by_z_on_coset(std::vector &P) { + void divide_by_z_on_coset(std::vector &P) override { const field_value_type coset = fields::arithmetic_params::multiplicative_generator; const field_value_type coset_to_small_m = coset.pow(small_m); diff --git a/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp b/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp index e801219..f50f332 100755 --- a/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp +++ b/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp @@ -53,8 +53,11 @@ namespace nil { bool precomputation_sentinel; std::vector geometric_sequence; std::vector geometric_triangular_sequence; + field_value_type geometric_generator; void do_precomputation() { + geometric_generator = field_value_type(fields::arithmetic_params::geometric_generator); + geometric_sequence = std::vector(this->m, field_value_type::zero()); geometric_sequence[0] = field_value_type::one(); @@ -62,8 +65,7 @@ namespace nil { geometric_triangular_sequence[0] = field_value_type::one(); for (std::size_t i = 1; i < this->m; i++) { - geometric_sequence[i] = - geometric_sequence[i - 1] * fields::arithmetic_params::geometric_generator; + geometric_sequence[i] = geometric_sequence[i - 1] * geometric_generator; geometric_triangular_sequence[i] = geometric_triangular_sequence[i - 1] * geometric_sequence[i - 1]; } @@ -86,7 +88,7 @@ namespace nil { precomputation_sentinel = false; } - void fft(std::vector &a) { + void fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -120,7 +122,8 @@ namespace nil { a[i] = a[i] * T[i].inversed(); } } - void inverse_fft(std::vector &a) { + + void inverse_fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type::zero()); @@ -159,7 +162,8 @@ namespace nil { newton_to_monomial_basis_geometric(a, geometric_sequence, geometric_triangular_sequence, this->m); } - std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) { + + std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) override { /* Compute Lagrange polynomial of size m, with m+1 points (x_0, y_0), ... ,(x_m, y_m) */ /* Evaluate for x = t */ /* Return coeffs for each l_j(x) = (l / l_i[j]) * w[j] */ @@ -220,11 +224,11 @@ namespace nil { } std::vector evaluate_all_lagrange_polynomials(const typename std::vector::const_iterator &t_powers_begin, - const typename std::vector::const_iterator &t_powers_end) { - if(std::distance(t_powers_begin, t_powers_end) < this->m) { + const typename std::vector::const_iterator &t_powers_end) override { + if(std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) { throw std::invalid_argument("geometric_sequence_radix2: expected std::distance(t_powers_begin, t_powers_end) >= this->m"); } - + /* Compute Lagrange polynomial of size m, with m+1 points (x_0, y_0), ... ,(x_m, y_m) */ /* Evaluate for x = t */ /* Return coeffs for each l_j(x) = (l / l_i[j]) * w[j] */ @@ -280,18 +284,18 @@ namespace nil { } std::vector result(this->m, value_type::zero()); - + for(std::size_t j = 0; j < l[0].size(); ++j) { result[0] = result[0] + t_powers_begin[j] * l[0][j]; } result[0] = result[0] * g_i[0]; for (std::size_t i = 1; i < this->m; i++) { g_i[i] = g_i[i - 1] * g[this->m - i] * -g[i].inversed() * geometric_sequence[i]; - + for(std::size_t j = 0; j < l[i].size(); ++j) { result[i] = result[i] + t_powers_begin[j] * l[i][j]; } - + result[i] = result[i] * (r_i * g_i[i]); r_i *= r; } @@ -299,13 +303,19 @@ namespace nil { return result; } - field_value_type get_domain_element(const std::size_t idx) { + // This one is not the unity root actually, but it's ok for our purposes. + const field_value_type& get_unity_root() override { + return geometric_generator; + } + + field_value_type get_domain_element(const std::size_t idx) override { if (!precomputation_sentinel) do_precomputation(); return this->geometric_sequence[idx]; } - field_value_type compute_vanishing_polynomial(const field_value_type &t) { + + field_value_type compute_vanishing_polynomial(const field_value_type &t) override { if (!precomputation_sentinel) do_precomputation(); @@ -317,7 +327,8 @@ namespace nil { } return Z; } - polynomial get_vanishing_polynomial() { + + polynomial get_vanishing_polynomial() override { if (!precomputation_sentinel) do_precomputation(); @@ -327,7 +338,8 @@ namespace nil { } return z; } - void add_poly_z(const field_value_type &coeff, std::vector &H) { + + void add_poly_z(const field_value_type &coeff, std::vector &H) override { if (H.size() != this->m + 1) throw std::invalid_argument("geometric: expected H.size() == this->m+1"); @@ -351,7 +363,8 @@ namespace nil { H[i] += (x[i] * coeff); } } - void divide_by_z_on_coset(std::vector &P) { + + void divide_by_z_on_coset(std::vector &P) override { const field_value_type coset = field_value_type( fields::arithmetic_params::multiplicative_generator); /* coset in geometric sequence? */ diff --git a/include/nil/crypto3/math/domains/step_radix2_domain.hpp b/include/nil/crypto3/math/domains/step_radix2_domain.hpp index e469f39..7dce339 100755 --- a/include/nil/crypto3/math/domains/step_radix2_domain.hpp +++ b/include/nil/crypto3/math/domains/step_radix2_domain.hpp @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -46,33 +47,44 @@ namespace nil { class step_radix2_domain : public evaluation_domain { typedef typename FieldType::value_type field_value_type; typedef ValueType value_type; - + typedef std::pair, std::vector> cache_type; + + std::unique_ptr small_fft_cache, big_fft_cache; + + void create_fft_cache() { + small_fft_cache = std::make_unique( + std::make_pair(std::vector(), std::vector())); + big_fft_cache = std::make_unique( + std::make_pair(std::vector(), std::vector())); + detail::create_fft_cache(big_m, big_omega, big_fft_cache->first); + detail::create_fft_cache(big_m, big_omega.inversed(), big_fft_cache->second); + detail::create_fft_cache(small_m, small_omega, small_fft_cache->first); + detail::create_fft_cache(small_m, small_omega.inversed(), small_fft_cache->second); + } public: typedef FieldType field_type; - std::size_t big_m; - std::size_t small_m; - field_value_type omega; - field_value_type big_omega; - field_value_type small_omega; - - step_radix2_domain(const std::size_t m) : evaluation_domain(m) { + const std::size_t big_m; + const std::size_t small_m; + const field_value_type omega; + const field_value_type big_omega; + const field_value_type small_omega; + + step_radix2_domain(const std::size_t m) + : evaluation_domain(m), + big_m(1ul << (static_cast(std::ceil(std::log2(m))) - 1)), + small_m(m - big_m), + omega(unity_root(1ul << static_cast(std::ceil(std::log2(m))))), + big_omega(omega.squared()), + small_omega(unity_root(small_m)) { if (m <= 1) throw std::invalid_argument("step_radix2(): expected m > 1"); - big_m = 1ul << (static_cast(std::ceil(std::log2(m))) - 1); - small_m = m - big_m; - if (small_m != 1ul << static_cast(std::ceil(std::log2(small_m)))) throw std::invalid_argument("step_radix2(): expected small_m == 1ul<(1ul << static_cast(std::ceil(std::log2(m)))); - - big_omega = omega.squared(); - small_omega = unity_root(small_m); } - void fft(std::vector &a) { + void fft(std::vector &a) override { if (a.size() != this->m) { if (a.size() < this->m) { a.resize(this->m, value_type()); @@ -100,8 +112,11 @@ namespace nil { } } - detail::basic_radix2_fft(c, omega.squared()); - detail::basic_radix2_fft(e, unity_root(small_m)); + if (small_fft_cache == nullptr) { + create_fft_cache(); + } + detail::basic_radix2_fft_cached(c, big_fft_cache->first); + detail::basic_radix2_fft_cached(e, small_fft_cache->first); for (std::size_t i = 0; i < big_m; ++i) { a[i] = c[i]; @@ -111,15 +126,18 @@ namespace nil { a[i + big_m] = e[i]; } } - void inverse_fft(std::vector &a) { + void inverse_fft(std::vector &a) override { if (a.size() != this->m) throw std::invalid_argument("step_radix2: expected a.size() == this->m"); std::vector U0(a.begin(), a.begin() + big_m); std::vector U1(a.begin() + big_m, a.end()); - detail::basic_radix2_fft(U0, omega.squared().inversed()); - detail::basic_radix2_fft(U1, unity_root(small_m).inversed()); + if (small_fft_cache == nullptr) { + create_fft_cache(); + } + detail::basic_radix2_fft_cached(U0, big_fft_cache->second); + detail::basic_radix2_fft_cached(U1, small_fft_cache->second); const field_value_type U0_size_inv = field_value_type(big_m).inversed(); for (std::size_t i = 0; i < big_m; ++i) { @@ -170,7 +188,7 @@ namespace nil { } } - std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) { + std::vector evaluate_all_lagrange_polynomials(const field_value_type &t) override { std::vector inner_big = detail::basic_radix2_evaluate_all_lagrange_polynomials(big_m, t); std::vector inner_small = @@ -199,17 +217,17 @@ namespace nil { } std::vector evaluate_all_lagrange_polynomials(const typename std::vector::const_iterator &t_powers_begin, - const typename std::vector::const_iterator &t_powers_end) { - if(std::distance(t_powers_begin, t_powers_end) < this-> m) { + const typename std::vector::const_iterator &t_powers_end) override { + if(std::size_t(std::distance(t_powers_begin, t_powers_end)) < this->m) { throw std::invalid_argument("extended_radix2: expected std::distance(t_powers_begin, t_powers_end) >= this->m"); } - + basic_radix2_domain basic_domain_big(big_m); std::vector inner_big = basic_domain_big.evaluate_all_lagrange_polynomials(t_powers_begin, t_powers_end); std::vector inner_big_times_t_to_small_m = basic_domain_big.evaluate_all_lagrange_polynomials(t_powers_begin + small_m, t_powers_end); - + basic_radix2_domain basic_domain_small(small_m); std::vector omega_inverse_t_powers(small_m); std::vector omega_inverse_t_powers_times_t_to_big_m(small_m); @@ -226,7 +244,7 @@ namespace nil { basic_domain_small.evaluate_all_lagrange_polynomials(omega_inverse_t_powers_times_t_to_big_m.cbegin(), omega_inverse_t_powers_times_t_to_big_m.cend()); std::vector result(this->m, value_type::zero()); - + const field_value_type omega_to_small_m = omega.pow(small_m); const field_value_type big_omega_to_small_m = big_omega.pow(small_m); field_value_type elt = field_value_type::one(); @@ -244,7 +262,11 @@ namespace nil { return result; } - field_value_type get_domain_element(const std::size_t idx) { + const field_value_type& get_unity_root() override { + return omega; + } + + field_value_type get_domain_element(const std::size_t idx) override { if (idx < big_m) { return big_omega.pow(idx); } else { @@ -252,11 +274,11 @@ namespace nil { } } - field_value_type compute_vanishing_polynomial(const field_value_type &t) { + field_value_type compute_vanishing_polynomial(const field_value_type &t) override { return (t.pow(big_m) - field_value_type::one()) * (t.pow(small_m) - omega.pow(small_m)); } - polynomial get_vanishing_polynomial() { + polynomial get_vanishing_polynomial() override { polynomial z(big_m + small_m + 1, field_value_type::zero()); field_value_type omega_to_small_m = omega.pow(small_m); z[big_m + small_m] = field_value_type::one(); @@ -267,7 +289,7 @@ namespace nil { return z; } - void add_poly_z(const field_value_type &coeff, std::vector &H) { + void add_poly_z(const field_value_type &coeff, std::vector &H) override { // if (H.size() != this->m + 1) // throw std::invalid_argument("step_radix2: expected H.size() == this->m+1"); @@ -278,7 +300,7 @@ namespace nil { H[small_m] -= coeff; H[0] += coeff * omega_to_small_m; } - void divide_by_z_on_coset(std::vector &P) { + void divide_by_z_on_coset(std::vector &P) override { // (c^{2^k}-1) * (c^{2^r} * w^{2^{r+1}*i) - w^{2^r}) const field_value_type coset = fields::arithmetic_params::multiplicative_generator; diff --git a/include/nil/crypto3/math/polynomial/basic_operations.hpp b/include/nil/crypto3/math/polynomial/basic_operations.hpp index b0272bb..3747964 100755 --- a/include/nil/crypto3/math/polynomial/basic_operations.hpp +++ b/include/nil/crypto3/math/polynomial/basic_operations.hpp @@ -268,7 +268,7 @@ namespace nil { const field_value_type sconst = field_value_type(n).inversed(); - for (std::size_t i = 0; i < n; ++i) { + for(std::size_t i = 0; i < n; ++i) { c[i] = c[i] * sconst; } @@ -281,11 +281,7 @@ namespace nil { * [Bostan, Lecerf, & Schost, 2003. Tellegen's Principle in Practice, on page 39]. */ template - AlgebraicRange - transpose_multiplication(const std::size_t &n, const AlgebraicRange &a, const FieldRange &c) { - typedef - typename std::iterator_traits()))>::value_type value_type; + AlgebraicRange transpose_multiplication(const std::size_t &n, const AlgebraicRange &a, const FieldRange &c) { const std::size_t m = a.size(); // if (c.size() - 1 > m + n) diff --git a/include/nil/crypto3/math/polynomial/basis_change.hpp b/include/nil/crypto3/math/polynomial/basis_change.hpp index ec3df8a..12d7bae 100755 --- a/include/nil/crypto3/math/polynomial/basis_change.hpp +++ b/include/nil/crypto3/math/polynomial/basis_change.hpp @@ -157,7 +157,6 @@ namespace nil { size_t n) { typedef typename Range::value_type value_type; - typedef typename FieldType::value_type field_vale_type; std::size_t m = log2(n); // if (T.size() != m + 1u) diff --git a/include/nil/crypto3/math/polynomial/evaluate.hpp b/include/nil/crypto3/math/polynomial/evaluate.hpp index 92f8a89..6d6a3c3 100755 --- a/include/nil/crypto3/math/polynomial/evaluate.hpp +++ b/include/nil/crypto3/math/polynomial/evaluate.hpp @@ -47,7 +47,7 @@ namespace nil { template inline FieldValueType evaluate_polynomial(ContiguousIterator first, ContiguousIterator last, const FieldValueType &t, std::size_t m) { - BOOST_ASSERT(std::distance(first, last) == m); + BOOST_ASSERT(std::size_t(std::distance(first, last)) == m); return boost::math::tools::evaluate_polynomial(&*first, t, m); } @@ -77,7 +77,7 @@ namespace nil { BOOST_STATIC_ASSERT(std::is_same::value); - if (m != std::distance(first, last)) { + if (m != std::size_t(std::distance(first, last))) { throw std::invalid_argument("expected m == domain.size()"); } if (idx >= m) { diff --git a/include/nil/crypto3/math/polynomial/polynomial_dfs.hpp b/include/nil/crypto3/math/polynomial/polynomial_dfs.hpp index 3fd4a78..9767203 100644 --- a/include/nil/crypto3/math/polynomial/polynomial_dfs.hpp +++ b/include/nil/crypto3/math/polynomial/polynomial_dfs.hpp @@ -28,11 +28,15 @@ #define CRYPTO3_MATH_POLYNOMIAL_POLYNOM_DFT_HPP #include +#include +#include #include #include #include +#include #include +#include #include namespace nil { @@ -86,15 +90,17 @@ namespace nil { template polynomial_dfs(size_t d, InputIterator first, InputIterator last) : val(first, last), _d(d) { - BOOST_ASSERT_MSG(std::distance(first, last) == detail::power_of_two(std::distance(first, last)), - "DFS optimal polynomial size must be a power of two"); + BOOST_ASSERT_MSG( + std::size_t(std::distance(first, last)) == detail::power_of_two(std::distance(first, last)), + "DFS optimal polynomial size must be a power of two"); } template polynomial_dfs(size_t d, InputIterator first, InputIterator last, const allocator_type& a) : val(first, last, a), _d(d) { - BOOST_ASSERT_MSG(std::distance(first, last) == detail::power_of_two(std::distance(first, last)), - "DFS optimal polynomial size must be a power of two"); + BOOST_ASSERT_MSG( + std::size_t(std::distance(first, last)) == detail::power_of_two(std::distance(first, last)), + "DFS optimal polynomial size must be a power of two"); } ~polynomial_dfs() = default; @@ -121,7 +127,7 @@ namespace nil { , _d(x._d) { } - polynomial_dfs(polynomial_dfs&& x, const allocator_type& a) + polynomial_dfs(polynomial_dfs&& x, const allocator_type& a) : val(std::move(x.val), a) , _d(x._d) { } @@ -148,21 +154,6 @@ namespace nil { return *this; } - // polynomial_dfs& operator=(const container_type& x) { - // val = x; - // return *this; - // } - // - // polynomial_dfs& operator=(container_type&& x) { - // val = x; - // return *this; - // } - - // polynomial_dfs& operator=(std::initializer_list il) { - // val.assign(il.begin(), il.end()); - // return *this; - // } - bool operator==(const polynomial_dfs& rhs) const { return val == rhs.val && _d == rhs._d; } @@ -170,20 +161,6 @@ namespace nil { return !(rhs == *this && _d == rhs._d); } - // template - // typename std::iterator_traits::reference assign(InputIterator first, - // InputIterator last) { - // return val.assign(first, last); - // } - // - // void assign(size_type n, const_reference u) { - // return val.assign(n, u); - // } - // - // void assign(std::initializer_list il) { - // assign(il.begin(), il.end()); - // } - allocator_type get_allocator() const BOOST_NOEXCEPT { return this->val.__alloc(); } @@ -350,7 +327,9 @@ namespace nil { val.clear(); } - void resize(size_type _sz) { + void resize(size_type _sz, + std::shared_ptr> old_domain = nullptr, + std::shared_ptr> new_domain = nullptr) { if (this->size() == _sz) { return; @@ -358,13 +337,23 @@ namespace nil { BOOST_ASSERT_MSG(_sz >= _d, "Resizing DFS polynomial to a size less than degree is prohibited: can't restore the polynomial in the future."); if (this->degree() == 0) { // Here we cannot write this->val.resize(_sz, this->val[0]), it will segfault. - auto value = this->val[0]; + auto value = this->val[0]; this->val.resize(_sz, value); } else { typedef typename value_type::field_type FieldType; - make_evaluation_domain(this->size())->inverse_fft(this->val); + if (old_domain == nullptr) { + old_domain = make_evaluation_domain(this->size()); + } else { + BOOST_ASSERT_MSG(old_domain->size() == this->size(), "Old domain size is not equal to the polynomial size"); + } + old_domain->inverse_fft(this->val); this->val.resize(_sz, FieldValueType::zero()); - make_evaluation_domain(_sz)->fft(this->val); + if (new_domain == nullptr) { + new_domain = make_evaluation_domain(_sz); + } else { + BOOST_ASSERT_MSG(new_domain->size() == _sz, "New domain size is not equal to the polynomial size"); + } + new_domain->fft(this->val); } } @@ -393,7 +382,7 @@ namespace nil { } return true; } - + /** * Returns true if polynomial is a one polynomial. */ @@ -406,11 +395,11 @@ namespace nil { } inline static polynomial_dfs zero() { - return polynomial_dfs(); + return polynomial_dfs(); } inline static polynomial_dfs one() { - return polynomial_dfs(0, size_type(1), value_type(1)); + return polynomial_dfs(0, size_type(1), value_type(1)); } /** @@ -429,24 +418,12 @@ namespace nil { */ polynomial_dfs operator+(const polynomial_dfs& other) const { polynomial_dfs result = *this; - if (other.size() > this->size()) { - result.resize(other.size()); - } - result._d = std::max(this->_d, other._d); - if (this->size() > other.size()) { - polynomial_dfs tmp(other); - tmp.resize(this->size()); - std::transform(tmp.begin(), tmp.end(), result.begin(), result.begin(), - std::plus()); - return result; - } - std::transform(other.begin(), other.end(), result.begin(), result.begin(), - std::plus()); + result += other; return result; } /** - * Computes the standard polynomial addition, polynomial A + polynomial B, + * Computes the standard polynomial addition, polynomial A + polynomial B, * and stores result in polynomial A. */ polynomial_dfs& operator+=(const polynomial_dfs& other) { @@ -466,16 +443,16 @@ namespace nil { } /** - * Computes polynomial A + constant c, + * Computes polynomial A + constant c, * and stores result in polynomial A. */ polynomial_dfs& operator+=(const FieldValueType& c) { for( auto it = this->begin(); it!=this->end(); it++) *it += c; return *this; } - + /** - * Computes polynomial A - constant c, + * Computes polynomial A - constant c, * and stores result in polynomial A. */ polynomial_dfs operator-() const { @@ -485,29 +462,17 @@ namespace nil { } /** - * Computes the standard polynomial subtraction, polynomial A - polynomial B, + * Computes the standard polynomial subtraction, polynomial A - polynomial B, * and stores result in polynomial C. */ polynomial_dfs operator-(const polynomial_dfs& other) const { polynomial_dfs result = *this; - if (other.size() > this->size()) { - result.resize(other.size()); - } - result._d = std::max(_d, other._d); - if (this->size() > other.size()) { - polynomial_dfs tmp(other); - tmp.resize(this->size()); - std::transform(result.begin(), result.end(), tmp.begin(), result.begin(), - std::minus()); - return result; - } - std::transform(result.begin(), result.end(), other.begin(), result.begin(), - std::minus()); + result -= other; return result; } /** - * Computes the standard polynomial subtraction, polynomial A - polynomial B, + * Computes the standard polynomial subtraction, polynomial A - polynomial B, * and stores result in polynomial A. */ polynomial_dfs& operator-=(const polynomial_dfs& other) { @@ -527,7 +492,7 @@ namespace nil { } /** - * Computes tpolynomial A - constant c + * Computes tpolynomial A - constant c * and stores result in polynomial A. */ polynomial_dfs& operator-=(const FieldValueType& c) { @@ -536,43 +501,38 @@ namespace nil { } /** - * Perform the multiplication of two polynomials, polynomial A * polynomial B, + * Perform the multiplication of two polynomials, polynomial A * polynomial B, * and stores result in polynomial C. */ polynomial_dfs operator*(const polynomial_dfs& other) const { polynomial_dfs result = *this; + result *= other; - size_t polynomial_s = - detail::power_of_two(std::max({this->size(), other.size(), this->degree() + other.degree() + 1})); - - if (result.size() < polynomial_s) { - result.resize(polynomial_s); - } - // Change the degree only here, after a possible resize, otherwise we have a polynomial - // with a high degree but small size, which sometimes segfaults. - result._d = this->degree() + other.degree(); - - if (other.size() < polynomial_s) { - polynomial_dfs tmp(other); - tmp.resize(polynomial_s); - std::transform(tmp.begin(), tmp.end(), result.begin(), result.begin(), std::multiplies()); - return result; - } - std::transform(other.begin(), other.end(), result.begin(), result.begin(), std::multiplies()); - return result; } /** - * Perform the multiplication of two polynomials, polynomial A * polynomial B, + * Perform the multiplication of two polynomials, polynomial A * polynomial B, * and stores result in polynomial A. */ polynomial_dfs& operator*=(const polynomial_dfs& other) { - size_t polynomial_s = + return cached_multiplication(other); + } + + /** + * Performs multiplication of two polynomials, but with domain caches + */ + polynomial_dfs& cached_multiplication( + const polynomial_dfs& other, + std::shared_ptr> domain = nullptr, + std::shared_ptr> other_domain = nullptr, + std::shared_ptr> new_domain = nullptr) { + + const size_t polynomial_s = detail::power_of_two(std::max({this->size(), other.size(), this->degree() + other.degree() + 1})); if (this->size() < polynomial_s) { - this->resize(polynomial_s); + this->resize(polynomial_s, domain, new_domain); } // Change the degree only here, after a possible resize, otherwise we have a polynomial @@ -581,7 +541,7 @@ namespace nil { if (other.size() < polynomial_s) { polynomial_dfs tmp(other); - tmp.resize(polynomial_s); + tmp.resize(polynomial_s, other_domain, new_domain); std::transform(tmp.begin(), tmp.end(), this->begin(), this->begin(), std::multiplies()); return *this; @@ -589,16 +549,16 @@ namespace nil { std::transform(this->begin(), this->end(), other.begin(), this->begin(), std::multiplies()); return *this; } - + /** - * Perform the multiplication of two polynomials, polynomial A * constant alpha, + * Perform the multiplication of two polynomials, polynomial A * constant alpha, * and stores result in polynomial A. */ polynomial_dfs& operator*=(const FieldValueType& alpha) { for( auto it = this->begin(); it!=this->end(); it++) *it *= alpha; return *this; } - + /** * Perform the standard Euclidean Division algorithm. * Input: Polynomial A, Polynomial B, where A / B @@ -650,18 +610,23 @@ namespace nil { detail::basic_radix2_fft(val, omega); } - std::vector coefficients() const { + std::vector coefficients( + std::shared_ptr> domain = nullptr) const { typedef typename value_type::field_type FieldType; value_type omega = unity_root(this->size()); std::vector tmp(this->begin(), this->end()); - detail::basic_radix2_fft(tmp, omega.inversed()); + if (domain == nullptr) { + detail::basic_radix2_fft(tmp, omega.inversed()); + const value_type sconst = value_type(this->size()).inversed(); + std::transform(tmp.begin(), + tmp.end(), + tmp.begin(), + std::bind(std::multiplies(), sconst, std::placeholders::_1)); + } else { + domain->inverse_fft(tmp); + } - const value_type sconst = value_type(this->size()).inversed(); - std::transform(tmp.begin(), - tmp.end(), - tmp.begin(), - std::bind(std::multiplies(), sconst, std::placeholders::_1)); size_t r_size = tmp.size(); while (r_size > 1 && tmp[r_size - 1] == FieldValueType::zero()) { --r_size; @@ -677,13 +642,13 @@ namespace nil { polynomial_dfs power_of_2 = *this; size_t expected_size = detail::power_of_two( - std::max({this->size(), this->degree() * power + 1})); + std::max({this->size(), this->degree() * power + 1})); power_of_2.resize(expected_size); polynomial_dfs result(0, expected_size, FieldValueType::one()); while (power) { if (power % 2 == 1) { result *= power_of_2; - } + } power /= 2; if (power == 0) break; @@ -715,7 +680,7 @@ namespace nil { } return result; } - + template, typename = typename std::enable_if::value>::type> @@ -796,6 +761,69 @@ namespace nil { return os; } + template + static inline polynomial_dfs polynomial_product( + std::vector> multipliers) { + // Pre-create all the domains. We could do this on-the-go, but we want this function to be more + // parallelization-friendly. This single-threaded version may look a bit complicated, + // but it's now very similar to what we have in parallel code. + std::unordered_map>> domain_cache; + + std::size_t min_domain_size = std::numeric_limits::max(); + std::size_t max_domain_size = 0; + std::size_t total_degree = 0; + for (std::size_t i = 0; i < multipliers.size(); i++) { + min_domain_size = std::min(min_domain_size, multipliers[i].size()); + max_domain_size = std::max(max_domain_size, multipliers[i].size()); + total_degree += multipliers[i].degree(); + } + max_domain_size = std::max(max_domain_size, detail::power_of_two(total_degree + 1)); + + std::vector needed_domain_sizes; + for (std::size_t i = min_domain_size; i <= max_domain_size; i *= 2) { + needed_domain_sizes.push_back(i); + // On the next line I want to create the tree structure, then create the evaluation domains. + // This way filling of this structure can be done in parallel. + domain_cache[i] = nullptr; + } + + // This loop will run in parallel. + for (const auto& domain_size: needed_domain_sizes) { + domain_cache[domain_size] = make_evaluation_domain(domain_size); + } + + for (std::size_t stride = 1; stride < multipliers.size(); stride <<= 1) { + const std::size_t double_stride = stride << 1; + // This loop will run in parallel. + std::size_t max_i = (multipliers.size() - stride) / double_stride; + if ((multipliers.size() - stride) % double_stride != 0) + max_i++; + + for(std::size_t i = 0; i < max_i; i++) { + std::size_t index1 = i * double_stride; + std::size_t index2 = index1 + stride; + + const std::size_t current_domain_size = multipliers[index1].size(); + const std::size_t next_domain_size = multipliers[index2].size(); + const std::size_t new_domain_size = + detail::power_of_two(std::max( + {current_domain_size, + next_domain_size, + multipliers[index1].degree() + multipliers[index2].degree() + 1})); + + multipliers[index1].cached_multiplication( + multipliers[index2], + domain_cache[current_domain_size], + domain_cache[next_domain_size], + domain_cache[new_domain_size]); + + // Free the memory we are not going to use anymore. + multipliers[index2] = polynomial_dfs(); + } + } + return multipliers[0]; + } + } // namespace math } // namespace crypto3 } // namespace nil diff --git a/include/nil/crypto3/math/polynomial/polynomial_dfs_view.hpp b/include/nil/crypto3/math/polynomial/polynomial_dfs_view.hpp index 2d45629..77b1411 100644 --- a/include/nil/crypto3/math/polynomial/polynomial_dfs_view.hpp +++ b/include/nil/crypto3/math/polynomial/polynomial_dfs_view.hpp @@ -333,8 +333,6 @@ namespace nil { FieldValueType evaluate(const FieldValueType& value) const { - typedef typename value_type::field_type FieldType; - std::vector tmp = this->coefficients(); FieldValueType result = 0; auto end = tmp.end(); diff --git a/include/nil/crypto3/math/polynomial/shift.hpp b/include/nil/crypto3/math/polynomial/shift.hpp index 095e9a6..1eec617 100644 --- a/include/nil/crypto3/math/polynomial/shift.hpp +++ b/include/nil/crypto3/math/polynomial/shift.hpp @@ -38,7 +38,7 @@ namespace nil { const FieldValueType &x) { polynomial f_shifted(f); FieldValueType x_power = x; - for (int i = 1; i < f.size(); i++) { + for (std::size_t i = 1; i < f.size(); i++) { f_shifted[i] = f_shifted[i] * x_power; x_power *= x; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e8e021c..dc6bf7b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,19 +20,31 @@ cm_test_link_libraries(${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME} ${Boost_LIBRARIES}) macro(define_math_test name) - cm_test(NAME math_${name}_test SOURCES ${name}.cpp) + set(test_name "math_${name}_test") - target_include_directories(math_${name}_test PRIVATE + set(additional_args "") + if(ENABLE_JUNIT_TEST_OUTPUT) + set(TEST_RESULTS_DIR "${CMAKE_CURRENT_BINARY_DIR}/junit_results") + set(TEST_LOGS_DIR "${TEST_RESULTS_DIR}/logs") + file(MAKE_DIRECTORY ${TEST_LOGS_DIR}) + + set(additional_args "--log_format=JUNIT" + "--log_sink=${TEST_LOGS_DIR}/${test_name}.xml") + endif() + + cm_test(NAME ${test_name} SOURCES ${name}.cpp ARGS ${additional_args}) + + target_include_directories(${test_name} PRIVATE "$" "$" ${Boost_INCLUDE_DIRS}) - set_target_properties(math_${name}_test PROPERTIES CXX_STANDARD 17) + set_target_properties(${test_name} PROPERTIES CXX_STANDARD 17) get_target_property(target_type Boost::unit_test_framework TYPE) if(target_type STREQUAL "SHARED_LIB") - target_compile_definitions(math_${name}_test PRIVATE BOOST_TEST_DYN_LINK) + target_compile_definitions(${test_name} PRIVATE BOOST_TEST_DYN_LINK) elseif(target_type STREQUAL "STATIC_LIB") endif() @@ -41,13 +53,14 @@ endmacro() set(TESTS_NAMES "evaluation_domain" "expression" - "kronecker_substitution" + # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI "polynomial_arithmetic" "polynomial" "polynomial_view" "polynomial_dfs" "polynomial_dfs_view" - "lagrange_interpolation") + "lagrange_interpolation" + "basic_radix2_domain") foreach(TEST_NAME ${TESTS_NAMES}) define_math_test(${TEST_NAME}) diff --git a/test/basic_radix2_domain.cpp b/test/basic_radix2_domain.cpp new file mode 100644 index 0000000..00205da --- /dev/null +++ b/test/basic_radix2_domain.cpp @@ -0,0 +1,152 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2024 Dmitrii Tabalin +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE basic_radix2_domain_test + +#include +#include + +#include +#include +#include + + +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace nil::crypto3::algebra; +using namespace nil::crypto3::math; + +typedef fields::bls12_fr<381> FieldType; + +BOOST_AUTO_TEST_SUITE(basic_radix2_domain_test_suit) + +BOOST_AUTO_TEST_CASE(basic_radix2_domain_benchmark, *boost::unit_test::disabled()) { + using value_type = FieldType::value_type; + const std::size_t fft_count = 5; + const std::array fft_sizes = {1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20}; + std::array, fft_count> test_data; + std::chrono::time_point gen_start(std::chrono::high_resolution_clock::now()); + for (std::size_t i = 0; i < fft_count; ++i) { + test_data[i].resize(fft_sizes[i]); + for (std::size_t j = 0; j < fft_sizes[i]; ++j) { + test_data[i][j] = nil::crypto3::algebra::random_element(); + } + } + std::cout << "Generation: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - gen_start) + .count() + << " ms" << std::endl; + + // manually calculate the power, saving all the intermediate powers + std::chrono::time_point cache_start(std::chrono::high_resolution_clock::now()); + std::vector>> omega_powers(fft_count); + for (std::size_t i = 0; i < fft_count; i++) { + omega_powers[i].reset(new std::vector); + omega_powers[i]->resize(fft_sizes[i]); + (*omega_powers[i])[0] = unity_root(fft_sizes[i]); + for (std::size_t j = 1; j < fft_sizes[i]; j++) { + (*omega_powers[i])[j] = (*omega_powers[i])[j - 1] * (*omega_powers[i])[0]; + } + } + std::cout << "Cache: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - cache_start) + .count() + << " ms" << std::endl; + + std::chrono::time_point start_fft(std::chrono::high_resolution_clock::now()); + for (std::size_t i = 0; i < fft_count; ++i) { + nil::crypto3::math::detail::basic_radix2_fft( + test_data[i], + unity_root(fft_sizes[i])); //omega_powers[i]); + } + + std::cout << "Uncached FFT: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - start_fft) + .count() + << " ms" << std::endl; + + std::chrono::time_point start_cached(std::chrono::high_resolution_clock::now()); + for (std::size_t i = 0; i < fft_count; ++i) { + nil::crypto3::math::detail::basic_radix2_fft( + test_data[i], + unity_root(fft_sizes[i]), + omega_powers[i]); + } + std::cout << "Cached FFT: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - start_cached + ).count() + << " ms" << std::endl; +} + +BOOST_AUTO_TEST_CASE(fft_vs_multiplication_benchmark) { + using value_type = FieldType::value_type; + const std::size_t fft_size = 1 << 16; + std::vector test_data(fft_size); + std::chrono::time_point gen_start(std::chrono::high_resolution_clock::now()); + for (std::size_t i = 0; i < fft_size; ++i) { + test_data[i] = nil::crypto3::algebra::random_element(); + } + std::vector duped_data(test_data); + value_type random_mult = nil::crypto3::algebra::random_element(); + std::cout << "Generation: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - gen_start) + .count() + << " ms" << std::endl; + + std::chrono::time_point start_fft(std::chrono::high_resolution_clock::now()); + nil::crypto3::math::detail::basic_radix2_fft( + test_data, + unity_root(fft_size)); + std::cout << "FFT: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - start_fft) + .count() + << " ms" << std::endl; + + std::chrono::time_point start_mult(std::chrono::high_resolution_clock::now()); + + for (std::size_t i = 0; i < fft_size; ++i) { + duped_data[i] *= random_mult; + } + std::cout << "Multiplication: " + << std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - start_mult) + .count() + << " ms" << std::endl; +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/test/evaluation_domain.cpp b/test/evaluation_domain.cpp index 258a640..c635885 100755 --- a/test/evaluation_domain.cpp +++ b/test/evaluation_domain.cpp @@ -45,6 +45,11 @@ #include #include +#include +#include + +#include + #include #include #include @@ -73,9 +78,7 @@ void test_fft() { const std::size_t m = 4; std::vector f = {2, 5, 3, 8}; - std::shared_ptr> domain; - - domain = make_evaluation_domain(m); + std::shared_ptr> domain = make_evaluation_domain(m); std::vector a(f); @@ -101,15 +104,44 @@ void test_fft() { std::cout << "is_arithmetic_sequence_domain = " << detail::is_arithmetic_sequence_domain(m) << std::endl; } +template +void test_fft_performance(std::string field_name, size_t poly_size, size_t domain_size) { + typedef typename FieldType::value_type value_type; + + std::vector f; + + for (std::size_t i = 0; i < poly_size; ++i) { + f.push_back(nil::crypto3::algebra::random_element()); + } + + std::shared_ptr> domain = make_evaluation_domain(domain_size); + + std::vector a(f); + + std::chrono::time_point start(std::chrono::high_resolution_clock::now()); + + size_t SAMPLES = 20; + + for (std::size_t i = 0; i < SAMPLES; ++i) { + a.resize(poly_size); + domain->fft(a); + } + + auto elapsed = std::chrono::duration_cast( + std::chrono::high_resolution_clock::now() - start); + + size_t time = elapsed.count() / SAMPLES / 1000 / 1000; + std::cout << "FFT on " << field_name << " of poly size " << poly_size << " and domain size " << domain_size << " took: " << std::fixed << std::setprecision(3) + << time / 1000 << " sec " << time % 1000 << " ms" << std::endl; +} + template void test_inverse_fft_of_fft() { typedef typename FieldType::value_type value_type; const std::size_t m = 4; std::vector f = {2, 5, 3, 8}; - std::shared_ptr> domain; - - domain = make_evaluation_domain(m); + std::shared_ptr> domain = make_evaluation_domain(m); std::vector a(f); domain->fft(a); @@ -202,7 +234,7 @@ void test_fft_curve_elements() { typedef typename FieldType::value_type field_value_type; std::size_t m = 4; - + // Make sure the results are reproducible. std::srand(0); std::vector f(m); @@ -211,7 +243,7 @@ void test_fft_curve_elements() { for(std::size_t i = 0; i < m; ++i) { g[i] = value_type::one() * f[i]; } - + std::shared_ptr> domain; domain.reset(new EvaluationDomainType(m)); @@ -224,7 +256,7 @@ void test_fft_curve_elements() { curve_element_domain->fft(g); BOOST_CHECK_EQUAL(f.size(), g.size()); - + for(std::size_t i = 0; i < f.size(); ++i) { BOOST_CHECK(f[i] * value_type::one() == g[i]); } @@ -247,7 +279,7 @@ void test_inverse_fft_curve_elements() { for(std::size_t i = 0; i < m; ++i) { g[i] = value_type::one() * f[i]; } - + std::shared_ptr> domain; domain.reset(new EvaluationDomainType(m)); @@ -260,7 +292,7 @@ void test_inverse_fft_curve_elements() { curve_element_domain->inverse_fft(g); BOOST_CHECK_EQUAL(f.size(), g.size()); - + for(std::size_t i = 0; i < f.size(); ++i) { BOOST_CHECK(f[i] * value_type::one() == g[i]); } @@ -280,7 +312,7 @@ void test_lagrange_coefficients_from_powers(std::size_t m) { for(std::size_t i = 1; i < m; ++i) { t_powers[i] = t_powers[i-1] * t; } - + std::shared_ptr> domain; domain.reset(new EvaluationDomainType(m)); @@ -310,10 +342,10 @@ void test_lagrange_coefficients_curve_elements(std::size_t m) { for(std::size_t i = 1; i < m; ++i) { t_powers[i] = t_powers[i-1] * t; } - + std::shared_ptr> domain; domain.reset(new EvaluationDomainType(m)); - + std::shared_ptr> curve_element_domain; curve_element_domain.reset(new GroupEvaluationDomainType(m)); @@ -336,7 +368,7 @@ void test_get_vanishing_polynomial(std::size_t m) { // Make sure the results are reproducible. std::srand(0); field_value_type t = std::rand(); - + std::shared_ptr> domain; domain.reset(new EvaluationDomainType(m)); @@ -353,38 +385,49 @@ BOOST_AUTO_TEST_SUITE(fft_evaluation_domain_test_suite) BOOST_AUTO_TEST_CASE(fft) { test_fft>(); test_fft>(); + test_fft(); +} + +BOOST_AUTO_TEST_CASE(fft_perf_test, *boost::unit_test::disabled()) { + for (int i = 15; i <= 22; ++i) { + test_fft_performance>("BLS12<381> Scalar", 1 << i, 1 << i); + } } BOOST_AUTO_TEST_CASE(inverse_fft_to_fft) { test_inverse_fft_of_fft>(); test_inverse_fft_of_fft>(); + test_inverse_fft_of_fft(); } BOOST_AUTO_TEST_CASE(inverse_coset_ftt_to_coset_fft) { test_inverse_coset_ftt_of_coset_fft>(); test_inverse_coset_ftt_of_coset_fft>(); + test_inverse_coset_ftt_of_coset_fft(); } BOOST_AUTO_TEST_CASE(lagrange_coefficients) { test_lagrange_coefficients>(); test_lagrange_coefficients>(); + test_lagrange_coefficients(); } BOOST_AUTO_TEST_CASE(compute_z) { test_compute_z>(); test_compute_z>(); + test_compute_z(); } BOOST_AUTO_TEST_CASE(curve_elements_fft) { typedef curves::bls12<381>::scalar_field_type field_type; typedef curves::bls12<381>::g1_type<> group_type; using group_value_type = group_type::value_type; - + test_fft_curve_elements, basic_radix2_domain>(); - // not applicable for any m < 100 for this field + // not applicable for any m < 100 for this field // test_fft_curve_elements, @@ -407,7 +450,7 @@ BOOST_AUTO_TEST_CASE(curve_elements_inverse_fft) { typedef curves::bls12<381>::scalar_field_type field_type; typedef curves::bls12<381>::g1_type<> group_type; using group_value_type = group_type::value_type; - + test_inverse_fft_curve_elements, @@ -433,7 +476,7 @@ BOOST_AUTO_TEST_CASE(curve_elements_inverse_fft) { BOOST_AUTO_TEST_CASE(lagrange_coefficients_from_powers) { typedef curves::bls12<381>::scalar_field_type field_type; - + test_lagrange_coefficients_from_powers>(4); // not applicable for any m < 100 for this field, testing with base field instead @@ -451,12 +494,12 @@ BOOST_AUTO_TEST_CASE(curve_elements_lagrange_coefficients) { typedef curves::bls12<381>::scalar_field_type field_type; typedef curves::bls12<381>::g1_type<> group_type; using group_value_type = group_type::value_type; - + test_lagrange_coefficients_curve_elements, basic_radix2_domain>(4); - // not applicable for any m < 100 for this field + // not applicable for any m < 100 for this field // test_lagrange_coefficients_curve_elements, @@ -477,7 +520,7 @@ BOOST_AUTO_TEST_CASE(curve_elements_lagrange_coefficients) { BOOST_AUTO_TEST_CASE(get_vanishing_polynomial) { typedef curves::bls12<381>::scalar_field_type field_type; - + test_get_vanishing_polynomial>(4); // not applicable for any m < 100 for this field, testing with base field instead diff --git a/test/lagrange_interpolation.cpp b/test/lagrange_interpolation.cpp index 306b89c..8e66e15 100644 --- a/test/lagrange_interpolation.cpp +++ b/test/lagrange_interpolation.cpp @@ -69,8 +69,7 @@ BOOST_AUTO_TEST_CASE(polynomial_lagrange_interpolation_manual_test) { BOOST_AUTO_TEST_CASE(polynomial_lagrange_interpolation_random_test) { using field_type = fields::bls12_fr<381>; - using integral_type = typename field_type::integral_type; - auto one = field_type::value_type::one(); + auto one = field_type::value_type::one(); std::size_t n = std::rand() % 50 + 1; std::vector p_coeffs(2 * n); diff --git a/test/polynomial_dfs.cpp b/test/polynomial_dfs.cpp index 13f01cf..bde6aba 100644 --- a/test/polynomial_dfs.cpp +++ b/test/polynomial_dfs.cpp @@ -34,6 +34,7 @@ #include #include +#include #include #include @@ -1024,7 +1025,7 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_division) { } BOOST_AUTO_TEST_CASE(polynomial_dfs_shift) { - + polynomial_dfs a = { 0, {0x01, @@ -1039,7 +1040,7 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_shift) { std::vector a_coefficients = a.coefficients(); - polynomial a_normal = + polynomial a_normal = polynomial(a_coefficients); polynomial_dfs a_ans; @@ -1134,10 +1135,10 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_sub_constant) { 0x2292a9419fc78dfa76e936e03b4374c51f1b1702fa1c61018f51109bb9ab2d1b_cppui253, 0x69fd599ad882439cb1024001d88240000000bfffffffffff5_cppui253, 0x225563a322dac633d1c77fc14960780266d5b167b8a62ccc942322dcfdb216f7_cppui253}}; - + polynomial_dfs c2_res = { 5, - {0x73eda753299d7d483339d80809a1d80553bda402fffe5bfefffffffefffffffe_cppui253, + {0x73eda753299d7d483339d80809a1d80553bda402fffe5bfefffffffefffffffe_cppui253, 0x2292a9419fc78dfa76e936e03b4374c51f1b1702fa1c61018f51109bb9ab2d31_cppui253, 0x69fd599ad882439cb1024001d88240000000c00000000000b_cppui253, 0x225563a322dac633d1c77fc14960780266d5b167b8a62ccc942322dcfdb2170d_cppui253, @@ -1180,7 +1181,7 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_mul_constant) { 0x4aa11aeeffeaf5820679a035dc419bec892909fe768e4dbce45b05c05b05aaf_cppui253, 0x4eeb4b3446e9b3c3c79e9cedc89af1d7017a5b7f4dcd61879d800d4191e44d8a_cppui253 }}; - + polynomial_dfs c1 = a * c; polynomial_dfs c2 = c * a; a *= c; @@ -1210,7 +1211,7 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_pow_eq_test) { polynomial_dfs res = a; for (int i = 1; i < 7; ++i) res *= a; - + BOOST_CHECK_EQUAL(res, a.pow(7)); } @@ -1258,7 +1259,7 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_evaluate_after_resize_and_shift_test) { typename FieldType::value_type point = 0x10_cppui253; polynomial_dfs large_poly = small_poly; - small_poly = polynomial_shift(small_poly, -1, 8); + small_poly = polynomial_shift(small_poly, -1, 8); for (size_t new_size : {16, 32, 64, 128, 256, 512}) { large_poly.resize(new_size); auto large_poly_shifted = polynomial_shift(large_poly, -1, 8); @@ -1293,4 +1294,58 @@ BOOST_AUTO_TEST_CASE(polynomial_dfs_zero_one_test) { BOOST_CHECK((small_poly - one * small_poly).is_zero()); } +BOOST_AUTO_TEST_CASE(polynomial_dfs_multiplication_perf_test, *boost::unit_test::disabled()) { + size_t size = 131072 * 4; + + polynomial_dfs poly = { + size / 128, size, nil::crypto3::algebra::random_element()}; + + std::vector> poly4(64, poly); + + auto start = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < poly4.size(); ++i) { + for (int j = 0; j < 32; ++j) + poly4[i] *= poly; + } + + for (int i = 1; i < poly4.size(); ++i) { + BOOST_CHECK(poly4[i] == poly4[0]); + } + + // Record the end time + auto end = std::chrono::high_resolution_clock::now(); + + // Calculate the duration + auto duration = std::chrono::duration_cast(end - start); + + std::cout << "Multiplication time: " << duration.count() << " microseconds." << std::endl; +} + +BOOST_AUTO_TEST_CASE(polynomial_dfs_resize_perf_test, *boost::unit_test::disabled()) { + std::vector values; + size_t size = 131072 * 16; + for (int i = 0; i < size; i++) { + values.push_back(nil::crypto3::algebra::random_element()); + } + + polynomial_dfs poly = { + size - 1, values}; + + auto start = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < 10; ++i) { + auto poly2 = poly; + poly2.resize(8 * size); + + BOOST_CHECK(poly2.size() == 8 * size); + } + + // Record the end time + auto end = std::chrono::high_resolution_clock::now(); + + // Calculate the duration + auto duration = std::chrono::duration_cast(end - start); + + std::cout << "Resize time: " << duration.count() << " microseconds." << std::endl; +} + BOOST_AUTO_TEST_SUITE_END()