diff --git a/.github/scripts/install_llvm.sh b/.github/scripts/install_llvm.sh new file mode 100644 index 0000000..c00db83 --- /dev/null +++ b/.github/scripts/install_llvm.sh @@ -0,0 +1,4 @@ +#!/bin/bash + +sudo apt update +sudo apt install -y clang libc++-dev libclang-rt-dev lld diff --git a/.github/scripts/setup.sh b/.github/scripts/setup.sh new file mode 100644 index 0000000..ac23e92 --- /dev/null +++ b/.github/scripts/setup.sh @@ -0,0 +1,4 @@ +#!/bin/bash + +sudo apt update +sudo apt install -y ninja-build build-essential cmake diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 4f389d2..65e22ca 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,31 +1,55 @@ name: Test - -on: - push: - branches: - - main +run-name: Test +on: pull_request jobs: - ctest: + gcc: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v5 + - name: Set up + run: bash .github/scripts/setup.sh + - name: Install gcc + run: sudo apt install -y gcc + - name: Build tests + run: | + cmake --preset test_gcc -B build/debug/gcc + cd ./build/debug/gcc && ninja + - name: Run tests + run: | + ./build/debug/gcc/test/mvPolynomial_test + + clang: runs-on: ubuntu-latest - permissions: - contents: write steps: - - uses: actions/checkout@v4 - - name: Set up vcpkg + - uses: actions/checkout@v5 + - name: Set up + run: bash .github/scripts/setup.sh + - name: Install llvm + run: bash .github/scripts/install_llvm.sh + - name: Build tests run: | - git clone --depth 1 -b 2025.08.27 https://github.com/microsoft/vcpkg.git vcpkg - - name: Install build tools + cmake --preset test_clang -B build/debug/clang + cd ./build/debug/clang && ninja + - name: Run tests run: | - sudo apt update - sudo apt install -y ninja-build build-essential cmake - echo "$(pwd)/vcpkg" >> $GITHUB_PATH + ./build/debug/clang/test/mvPolynomial_test + + clang_memory: + # libstdc++ and libc++ have uninitialized variables, + # so I skip this job temporarily. + if: "false" + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v5 + - name: Set up + run: bash .github/scripts/setup.sh + - name: Install llvm + run: bash .github/scripts/install_llvm.sh - name: Build tests run: | - ./vcpkg/bootstrap-vcpkg.sh - cmake --preset CI_test - cd ./build/debug && ninja - env: - VCPKG_ROOT: ${{ github.workspace }}/vcpkg + cmake --preset test_clang_memory -B build/debug/clang/memory + cd ./build/debug/clang/memory && ninja - name: Run tests - run: ctest --test-dir ./build/debug/test + run: | + ./build/debug/clang/memory/test/mvPolynomial_test diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index 24582cd..0000000 --- a/.gitmodules +++ /dev/null @@ -1,3 +0,0 @@ -[submodule "platanus"] - path = platanus - url = https://github.com/sukeya/platanus.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b61a1a..f36d95b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,36 +7,122 @@ set(CMAKE_CXX_EXTENSIONS OFF) set(CMAKE_OPTIMIZE_DEPENDENCIES ON) # Library -# Boost -find_package(Boost CONFIG REQUIRED) -list(APPEND extra_headers ${Boost_INCLUDE_DIRS}) +include(FetchContent) +if (MVPOLYNOMIAL_BUILD_TESTS) + FetchContent_Declare( + eigen3 + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG 549bf8c75b6aae071cde2f28aa48f16ee3ae60b0 # 5.0.0 + OVERRIDE_FIND_PACKAGE + ) + # The reason why I use CACHE INTERNAL is written in the following link: + # https://cmake.org/cmake/help/latest/module/FetchContent.html#overriding-where-to-find-cmakelists-txt + set(EIGEN_BUILD_TESTING OFF CACHE INTERNAL "") + set(EIGEN_BUILD_CMAKE_PACKAGE ON CACHE INTERNAL "") + FetchContent_MakeAvailable(eigen3) -# Eigen -find_package(Eigen3 CONFIG REQUIRED NO_MODULE) -list(APPEND extra_libs Eigen3::Eigen) + FetchContent_Declare( + fmt + GIT_REPOSITORY https://github.com/fmtlib/fmt.git + GIT_TAG e424e3f2e607da02742f73db84873b8084fc714c # 12.0.0 + OVERRIDE_FIND_PACKAGE + ) + FetchContent_MakeAvailable(fmt) +else() + find_package(Eigen3 5.0.0 CONFIG REQUIRED) + find_package(fmt 12.0.0 CONFIG REQUIRED) +endif() -# fmt -find_package(fmt CONFIG REQUIRED) +list(APPEND extra_libs Eigen3::Eigen) list(APPEND extra_libs fmt::fmt) -add_library(mvPolynomial_compiler_flags INTERFACE) -target_compile_features(mvPolynomial_compiler_flags INTERFACE cxx_std_20) +# To avoid add_subdirectory, +FetchContent_Declare( + platanus + GIT_REPOSITORY https://github.com/sukeya/platanus.git + GIT_TAG 23a9ba9cd6bfe38b01fd7e87e28ae42dd8b2744f # origin/main +) +FetchContent_MakeAvailable(platanus) +list(APPEND extra_libs platanus) + add_library(${PROJECT_NAME} INTERFACE) +target_compile_features(${PROJECT_NAME} INTERFACE cxx_std_20) target_include_directories( ${PROJECT_NAME} INTERFACE - ${extra_headers} ${CMAKE_CURRENT_SOURCE_DIR}/include ) target_link_libraries( ${PROJECT_NAME} INTERFACE ${extra_libs} - mvPolynomial_compiler_flags ) -if (BUILD_TESTS) - include(CTest) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +set(gcc "$") +set(clang "$") +set(msvc "$") + +set(gcc_like $ #include +#include namespace mvPolynomial { /** @@ -12,33 +13,24 @@ namespace mvPolynomial { * \tparam IntType the type of elements of indices. * \tparam D the dimension of indices. */ -template -class IndexComparer { +template +class IndexComparer final { public: static_assert(D > 0, "IndexComparer: the dimension must be positive."); - using Index = IndexType; + using Int = Int_; + using Index = IndexType; + + static constexpr int dim = D; /** - * \brief If the elems of lhd are equal to that of rhd until i th time and lhd[i] is greater than - * rhd[i], return true: otherwise, false. + * \brief If lhd[i] == rhd[i] for i = 0, ..., N - 1 and lhd[N] > rhd[N], + * return greater: if lhd[N] < rhd[N], return less: otherwise, return equal. * \param[in] lhd an index * \param[in] rhd an index */ - constexpr bool operator()(const Index& lhd, const Index& rhd) const noexcept { - for (std::size_t i = 0; i != lhd.size(); ++i) { - auto comp = lhd[i] <=> rhd[i]; - if (comp > 0) { - return true; - } else if (comp < 0) { - return false; - } - } - return false; - } - - constexpr std::strong_ordering get_ordering(const Index& lhd, const Index& rhd) const noexcept { - for (std::size_t i = 0; i != lhd.size(); ++i) { + static constexpr std::strong_ordering Compare(const Index& lhd, const Index& rhd) { + for (std::size_t i = 0; i < lhd.size(); ++i) { auto comp = lhd[i] <=> rhd[i]; if (comp > 0) { return std::strong_ordering::greater; @@ -48,14 +40,23 @@ class IndexComparer { } return std::strong_ordering::equal; } -}; -template -class IndexComparer { - public: - using Index = IntType; - - constexpr bool operator()(Index lhd, Index rhd) const noexcept { return lhd > rhd; } + /** + * \brief If lhd[i] == rhd[i] for i = 0, ..., N - 1 and lhd[N] > rhd[N], + * return true: otherwise, false. + * \param[in] lhd an index + * \param[in] rhd an index + */ + bool operator()(const Index& lhd, const Index& rhd) const { + for (std::size_t i = 0; i < lhd.size(); ++i) { + if (lhd[i] > rhd[i]) { + return true; + } else if (lhd[i] < rhd[i]) { + return false; + } + } + return false; + } }; } // namespace mvPolynomial diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index 816e95c..5803a87 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -3,55 +3,97 @@ #include "mvPolynomial/type.hpp" #include "mvPolynomial/index_comparer.hpp" -#include "mvPolynomial/polynomial.hpp" #include #include -#include -#include -#include - -#include "boost/container/flat_map.hpp" -#include "boost/container/new_allocator.hpp" -#include "boost/tuple/tuple.hpp" -#include "boost/iterator/zip_iterator.hpp" +#include + #include "Eigen/Core" #include "fmt/core.h" +#include "platanus/btree_map.hpp" namespace mvPolynomial { +namespace { +void CheckAxis(int dim, int axis) { + if (axis < 0 || axis >= dim) { + throw std::runtime_error( + fmt::format("CheckAxis: Given axis {} must be in [0, {}).", axis, dim) + ); + } +} + +template +auto OfImpl(Iterator begin, Iterator end, int dim, int axis, const Coord& x) { + using Index = std::remove_cvref_t; + + CheckAxis(dim, axis); + + if (axis == dim - 1) { + auto last_coeff = begin->second; + auto last_index = Index{begin->first}; + for (auto it = std::next(begin); it != end; ++it) { + const auto& [next_index, next_coeff] = *it; + last_coeff *= std::pow(x[axis], last_index[axis] - next_index[axis]); + last_coeff += next_coeff; + last_index = next_index; + } + last_coeff *= x.pow(last_index.template cast()).prod(); + return last_coeff; + } else { + auto sum = typename Iterator::value_type::second_type(0); + while (true) { + const auto& [first_index, first_coeff] = *begin; + auto partition_point = std::partition_point( + begin, + end, + [axis, &first_index](const typename Iterator::value_type& pair) { + return pair.first[axis] == first_index[axis]; + } + ); + const auto& [p_index, p_coeff] = *partition_point; + sum += OfImpl(begin, partition_point, dim, axis + 1, x); + if (partition_point == end) { + // The calculation ends. + break; + } + begin = partition_point; + } + return sum; + } +} +} // namespace + template < std::signed_integral IntType, std::floating_point R, int D, - class Comparer = IndexComparer, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -class MVPolynomial { + class Allocator = std::allocator, R>>> +class MVPolynomial final { public: static_assert(D > 0, "MVPolynomial: the dimension must be greater than 0."); static constexpr int dim = D; - using index_type = typename AllocatorOrContainer::value_type::first_type; + // This setting is too strict, so I expect users to set tolerance. + static R tolerance; + + using index_type = IndexType; using coord_type = CoordType; private: - using IndexContainer = boost::container::flat_map; + using Comparer = IndexComparer; + using IndexContainer = platanus::btree_map; public: - using key_type = IndexContainer::key_type; - using mapped_type = IndexContainer::mapped_type; - using value_type = IndexContainer::value_type; - using movable_value_type = IndexContainer::movable_value_type; + using key_type = IndexContainer::key_type; + using value_type = IndexContainer::value_type; + using mapped_type = IndexContainer::mapped_type; + using coeff_type = mapped_type; using key_compare = IndexContainer::key_compare; using value_compare = IndexContainer::value_compare; - using sequence_type = IndexContainer::sequence_type; - - using allocator_type = IndexContainer::allocator_type; - using allocator_traits_type = IndexContainer::allocator_traits_type; - using stored_allocator_type = IndexContainer::stored_allocator_type; + using allocator_type = IndexContainer::allocator_type; using pointer = IndexContainer::pointer; using const_pointer = IndexContainer::const_pointer; @@ -67,16 +109,24 @@ class MVPolynomial { using reverse_iterator = IndexContainer::reverse_iterator; using const_reverse_iterator = IndexContainer::const_reverse_iterator; - explicit MVPolynomial(const allocator_type& allocator) : index2value_(allocator) { - CheckSelfIndexes(); - } + MVPolynomial(const MVPolynomial& other) = default; + MVPolynomial& operator=(const MVPolynomial& other) = default; + MVPolynomial(MVPolynomial&& other) = default; + MVPolynomial& operator=(MVPolynomial&& other) = default; + ~MVPolynomial() = default; - explicit MVPolynomial(const Comparer& comparer) : index2value_(comparer) { CheckSelfIndexes(); } + MVPolynomial() + : index2value_({ + {index_type::Zero(), 0} + }) {} - explicit MVPolynomial(const Comparer& comparer, const allocator_type& allocator) - : index2value_(comparer, allocator) { - CheckSelfIndexes(); - } + explicit MVPolynomial(const allocator_type& allocator) + : index2value_( + { + {index_type::Zero(), 0} + }, + allocator + ) {} template MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e) { @@ -89,146 +139,28 @@ class MVPolynomial { CheckSelfIndexes(); } - template - explicit MVPolynomial(InputIterator s, InputIterator e, const Comparer& c) - : index2value_(s, e, c) { - CheckSelfIndexes(); - } - - template - explicit MVPolynomial( - InputIterator s, InputIterator e, const Comparer& c, const allocator_type& a - ) - : index2value_(s, e, c, a) { - CheckSelfIndexes(); - } - - template - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, InputIterator s, InputIterator e - ) - : index2value_(o, s, e) { - CheckSelfIndexes(); - } - - template - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, - InputIterator s, - InputIterator e, - const Comparer& c - ) - : index2value_(o, s, e, c) { - CheckSelfIndexes(); - } - - template - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, - InputIterator s, - InputIterator e, - const Comparer& c, - const allocator_type& a - ) - : index2value_(o, s, e, c, a) { - CheckSelfIndexes(); - } - - template - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, - InputIterator s, - InputIterator e, - const allocator_type& a - ) - : index2value_(o, s, e, a) { - CheckSelfIndexes(); - } - - explicit MVPolynomial(std::initializer_list l) : index2value_(l) { - CheckSelfIndexes(); - } - - explicit MVPolynomial(std::initializer_list l, const allocator_type& a) + MVPolynomial(std::initializer_list l, const allocator_type& a = allocator_type{}) : index2value_(l, a) { CheckSelfIndexes(); } - explicit MVPolynomial(std::initializer_list l, const Comparer& c) - : index2value_(l, c) { - CheckSelfIndexes(); - } - - explicit MVPolynomial( - std::initializer_list l, const Comparer& c, const allocator_type& a - ) - : index2value_(l, c, a) { - CheckSelfIndexes(); - } - - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, std::initializer_list l - ) - : index2value_(o, l) { - CheckSelfIndexes(); - } - - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, - std::initializer_list l, - const Comparer& c - ) - : index2value_(o, l, c) { - CheckSelfIndexes(); - } - - explicit MVPolynomial( - boost::container::ordered_unique_range_t o, - std::initializer_list l, - const Comparer& c, - const allocator_type& a - ) - : index2value_(o, l, c, a) { - CheckSelfIndexes(); - } - explicit MVPolynomial(const MVPolynomial& m, const allocator_type& a) - : index2value_(m.index2value_, a) { - CheckSelfIndexes(); - } + : index2value_(m.index2value_, a) {} explicit MVPolynomial(MVPolynomial&& m, const allocator_type& a) - : index2value_(std::move(m.index2value_), a) { - CheckSelfIndexes(); - } + : index2value_(std::move(m.index2value_), a) {} MVPolynomial& operator=(std::initializer_list l) { - index2value_ = l; + index2value_.insert(l.begin(), l.end()); CheckSelfIndexes(); return *this; } explicit MVPolynomial(mapped_type r) { index2value_.at(index_type::Zero()) = r; } - MVPolynomial() = default; - MVPolynomial(const MVPolynomial& other) = default; - MVPolynomial& operator=(const MVPolynomial& other) = default; - MVPolynomial(MVPolynomial&& other) = default; - MVPolynomial& operator=(MVPolynomial&& other) = default; - virtual ~MVPolynomial() = default; - - static void CheckAxis(std::size_t axis) { - if (axis < 0 || axis >= dim) { - throw std::runtime_error( - fmt::format("CheckAxis: Given axis {} must be in [0, {}).", axis, dim) - ); - } - } - allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); } - auto get_stored_allocator() noexcept { return index2value_.get_stored_allocator(); } - auto get_stored_allocator() const noexcept { return index2value_.get_stored_allocator(); } + Comparer key_comp() const noexcept { return index2value_.key_comp(); } iterator begin() noexcept { return index2value_.begin(); } const_iterator begin() const noexcept { return index2value_.begin(); } @@ -248,18 +180,11 @@ class MVPolynomial { const_reverse_iterator crbegin() const noexcept { return index2value_.crbegin(); } const_reverse_iterator crend() const noexcept { return index2value_.crend(); } - bool empty() const noexcept { return index2value_.empty(); } - + bool empty() const noexcept { return index2value_.empty(); } size_type size() const noexcept { return index2value_.size(); } - size_type max_size() const noexcept { return index2value_.max_size(); } - size_type capacity() const noexcept { return index2value_.capacity(); } - void reserve(size_type size) { return index2value_.reserve(size); } - - void shrink_to_fit() { return index2value_.shrink_to_fit(); } - mapped_type& operator[](const key_type& index) { return index2value_[index]; } mapped_type& operator[](key_type&& index) { return index2value_[index]; } @@ -270,141 +195,30 @@ class MVPolynomial { return const_cast(this)->operator[](index); } - template - std::pair insert_or_assign(const key_type& i, M&& m) { - CheckIndex(i); - return index2value_.insert_or_assign(i, std::move(m)); - } - - template - std::pair insert_or_assign(key_type&& i, M&& m) { - CheckIndex(i); - return index2value_.insert_or_assign(std::move(i), std::move(m)); - } - - template - iterator insert_or_assign(const_iterator ci, const key_type& i, M&& m) { - CheckIndex(i); - return index2value_.insert_or_assign(ci, i, std::move(m)); - } - - template - iterator insert_or_assign(const_iterator ci, key_type&& i, M&& m) { - CheckIndex(i); - return index2value_.insert_or_assign(ci, std::move(i), std::move(m)); - } - - iterator nth(size_type size) noexcept { return index2value_.nth(size); } - const_iterator nth(size_type size) const noexcept { return index2value_.nth(size); } - - size_type index_of(iterator i) noexcept { return index2value_.index_of(i); } - size_type index_of(const_iterator ci) const noexcept { return index2value_.index_of(ci); } - mapped_type& at(const key_type& i) { return index2value_.at(i); } const mapped_type& at(const key_type& i) const { return index2value_.at(i); } - template - std::pair emplace(Args&&... args) { - return CheckIndexOfPairOfIterAndIsInserted(index2value_.emplace(std::move(args...))); - } - - template - iterator emplace_hint(const_iterator ci, Args&&... args) { - auto iter = index2value_.emplace_hint(ci, std::move(args)...); - CheckIndex(iter->first); - return iter; - } - - template - std::pair try_emplace(const key_type& i, Args&&... args) { - return CheckIndexOfPairOfIterAndIsInserted(index2value_.try_emplace(i, std::move(args...))); - } - - template - iterator try_emplace(const_iterator ci, const key_type& i, Args&&... args) { - auto iter = index2value_.try_emplace(ci, i, std::move(args...)); - CheckIndex(iter->first); - return iter; - } - - template - std::pair try_emplace(key_type&& i, Args&&... args) { - return CheckIndexOfPairOfIterAndIsInserted( - index2value_.try_emplace(std::move(i), std::move(args...)) - ); - } - - template - iterator try_emplace(const_iterator ci, key_type&& i, Args&&... args) { - auto iter = index2value_.try_emplace(ci, std::move(i), std::move(args...)); - CheckIndex(iter->first); - return iter; - } - - std::pair insert(const value_type& i_and_v) { - return CheckIndexOfPairOfIterAndIsInserted(index2value_.insert(i_and_v)); - } - - std::pair insert(value_type&& i_and_v) { - return CheckIndexOfPairOfIterAndIsInserted(index2value_.insert(std::move(i_and_v))); - } - - template - std::pair insert(Pair&& p) { - return CheckIndexOfPairOfIterAndIsInserted(index2value_.insert(std::move(p))); - } - - iterator insert(const_iterator ci, const value_type& i_and_v) { - auto iter = index2value_.insert(ci, i_and_v); - CheckIndex(iter->first); - return iter; - } + void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); } - iterator insert(const_iterator ci, value_type&& i_and_v) { - auto iter = index2value_.insert(ci, std::move(i_and_v)); - CheckIndex(iter->first); - return iter; - } + // I didn't want to add insert and erase, but for efficiency, I did. + std::pair insert(const value_type& x) { return index2value_.insert(x); } + std::pair insert(value_type&& x) { return index2value_.insert(std::move(x)); } - template - iterator insert(const_iterator ci, Pair&& p) { - auto iter = index2value_.insert(ci, std::move(p)); - CheckIndex(iter->first); - return iter; - } + iterator insert(iterator hint, const value_type& x) { return index2value_.insert(hint, x); } + iterator insert(iterator hint, value_type&& x) { return index2value_.insert(hint, std::move(x)); } template - void insert(InputIterator s, InputIterator e) { - index2value_.insert(s, e); - CheckSelfIndexes(); + void insert(InputIterator b, InputIterator e) { + index2value_.insert(b, e); } - template - void insert(boost::container::ordered_unique_range_t o, InputIterator s, InputIterator e) { - index2value_.insert(o, s, e); - CheckSelfIndexes(); + void insert(std::initializer_list list) { + index2value_.insert(list.begin(), list.end()); } - void insert(std::initializer_list l) { - index2value_.insert(l); - CheckSelfIndexes(); - } - - void insert(boost::container::ordered_unique_range_t o, std::initializer_list l) { - index2value_.insert(o, l); - CheckSelfIndexes(); - } - - iterator erase(const_iterator ci) { return index2value_.erase(ci); } - size_type erase(const value_type& i_and_v) { return index2value_.erase(i_and_v); } - iterator erase(const_iterator s, const_iterator e) { return index2value_.erase(s, e); } - - void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); } - - void clear() noexcept { index2value_.clear(); } - - key_compare key_comp() const { return index2value_.key_comp(); } - value_compare value_comp() const { return index2value_.value_comp(); } + size_type erase(const key_type& key) { return index2value_.erase(key); } + iterator erase(const iterator& iter) { return index2value_.erase(iter); } + size_type erase(const iterator& b, const iterator& e) { return index2value_.erase(b, e); } iterator find(const key_type& i) { return index2value_.find(i); } const_iterator find(const key_type& i) const { return index2value_.find(i); } @@ -417,12 +231,6 @@ class MVPolynomial { return index2value_.find(i); } - size_type count(const key_type& i) const { return index2value_.count(i); } - template - size_type count(const K& i) const { - return index2value_.count(i); - } - bool contains(const key_type& i) const { return index2value_.contains(i); } template bool contains(const K& i) const { @@ -451,38 +259,7 @@ class MVPolynomial { return index2value_.upper_bound(i); } - std::pair equal_range(const key_type& i) { - return index2value_.equal_range(i); - } - - std::pair equal_range(const key_type& i) const { - return index2value_.equal_range(i); - } - - template - std::pair equal_range(const K& i) { - return index2value_.equal_range(i); - } - - template - std::pair equal_range(const K& i) const { - return index2value_.equal_range(i); - } - - sequence_type extract_sequence() { return index2value_.extract_sequence(); } - - void adopt_sequence(sequence_type&& seq) { index2value_.adopt_sequence(std::move(seq)); } - void adopt_sequence(boost::container::ordered_unique_range_t o, sequence_type&& seq) { - index2value_.adopt_sequence(o, std::move(seq)); - } - - const sequence_type& sequence() const noexcept { return index2value_.sequence(); } - - reference front() { return *(index2value_.begin()); } - const_reference front() const { return *(index2value_.cbegin()); } - - reference back() { return *(index2value_.rbegin()); } - const_reference back() const { return *(index2value_.crbegin()); } + R operator()(const coord_type& x) const { return OfImpl(cbegin(), cend(), dim, 0, x); } MVPolynomial operator+() const { return *this; } @@ -504,713 +281,140 @@ class MVPolynomial { } // friend functions - // TODO consider tolerance. friend bool operator==(const MVPolynomial& l, const MVPolynomial& r) { - return l.index2value_ == r.index2value_; + if (l.size() != r.size()) { + return false; + } + auto l_it = l.cbegin(); + auto r_it = r.cbegin(); + for (std::size_t i = 0; i != l.size(); ++i) { + const auto& [l_idx, l_coeff] = *l_it; + const auto& [r_idx, r_coeff] = *r_it; + if ((l_idx != r_idx).any()) { + return false; + } + if (std::abs(l_coeff - r_coeff) >= tolerance) { + return false; + } + ++l_it; + ++r_it; + } + return true; } - // TODO consider tolerance. friend bool operator!=(const MVPolynomial& l, const MVPolynomial& r) { return !(l == r); } friend MVPolynomial operator+(const MVPolynomial& l, const MVPolynomial& r) { - auto comparer = l.key_comp(); - auto sum = MVPolynomial(comparer, l.get_allocator()); - sum.reserve(l.size() + r.size()); - auto l_it = l.begin(); - auto r_it = r.begin(); - // Like Merge sort algorithm, insert or sum data to sum. - while (l_it != l.end() && r_it != r.end()) { - const auto& [l_idx, l_v] = *l_it; - const auto& [r_idx, r_v] = *r_it; - if (comparer(l_idx, r_idx)) { - sum.emplace_hint(sum.end(), *l_it); - ++l_it; - } else if (comparer(r_idx, l_idx)) { - sum.emplace_hint(sum.end(), *r_it); - ++r_it; + auto sum = MVPolynomial(l); + for (const auto& [idx, coeff] : r) { + if (sum.contains(idx)) { + sum[idx] += coeff; } else { - sum.emplace_hint(sum.end(), l_idx, l_v + r_v); - ++l_it; - ++r_it; + sum[idx] = coeff; } } - // If r or l iterator don't equals its end, insert the extra to sum. - while (l_it != l.end()) { - sum.emplace_hint(sum.end(), *l_it); - ++l_it; - } - while (r_it != r.end()) { - sum.emplace_hint(sum.end(), *r_it); - ++r_it; - } return sum; } friend MVPolynomial operator-(const MVPolynomial& l, const MVPolynomial& r) { - auto comparer = l.key_comp(); - auto sub = MVPolynomial(comparer, l.get_allocator()); - sub.reserve(l.size() + r.size()); - auto l_it = l.begin(); - auto r_it = r.begin(); - // Like Merge sort algorithm, insert or subtract data to sub. - while (l_it != l.end() && r_it != r.end()) { - const auto& [l_idx, l_v] = *l_it; - const auto& [r_idx, r_v] = *r_it; - if (comparer(l_idx, r_idx)) { - sub.emplace_hint(sub.end(), *l_it); - ++l_it; - } else if (comparer(r_idx, l_idx)) { - sub.emplace_hint(sub.end(), r_idx, -r_v); - ++r_it; + auto sub = MVPolynomial(l); + for (const auto& [idx, coeff] : r) { + if (sub.contains(idx)) { + sub[idx] -= coeff; } else { - sub.emplace_hint(sub.end(), l_idx, l_v - r_v); - ++l_it; - ++r_it; + sub[idx] = -coeff; } } - // If r or l iterator don't equals its end, insert the extra to sub. - while (l_it != l.end()) { - sub.emplace_hint(sub.end(), *l_it); - ++l_it; - } - while (r_it != r.end()) { - const auto& [r_idx, r_v] = *r_it; - sub.emplace_hint(sub.end(), r_idx, -r_v); - ++r_it; - } return sub; } friend MVPolynomial operator*(const MVPolynomial& l, const MVPolynomial& r) { auto comparer = l.key_comp(); - auto mul = std::vector(); - mul.reserve(l.size() * r.size()); + + auto mul = MVPolynomial(l.get_allocator()); // Calculate all product of each l's term and r's term. for (const auto& l_p : l) { const auto& [l_idx, l_v] = l_p; for (const auto& r_p : r) { const auto& [r_idx, r_v] = r_p; - mul.emplace_back(l_idx + r_idx, l_v * r_v); + const auto idx = l_idx + r_idx; + const auto v = l_v * r_v; + if (mul.contains(idx)) { + mul[idx] += v; + } else { + mul[idx] = v; + } } } - std::sort(mul.begin(), mul.end(), [&comparer](const value_type& l, const value_type& r) { - return comparer(l.first, r.first); - }); - auto rm_begin = - std::unique(mul.begin(), mul.end(), [](const value_type& l, const value_type& r) { - return (l.first == r.first).all(); - }); - // For simplicity, I make a MVPolynomial now. - auto mp = MVPolynomial( - boost::container::ordered_unique_range_t(), - mul.begin(), - rm_begin, - comparer, - l.get_allocator() - ); - // Add duplicated elements to mp. - for (auto it = rm_begin; it != mul.end(); ++it) { - mp[it->first] += it->second; - } - return mp; + return mul; } - friend void swap(MVPolynomial& l, MVPolynomial& r) { swap(l.index2value_, r.index2value_); } - private: - void CheckIndex(const key_type& index) const { - if ((index < index_type::Zero()).any()) { - auto err_msg_stream = std::stringstream(); - for (auto i = 0; i != index.size() - 1; ++i) { - err_msg_stream << i << ", "; - } - err_msg_stream << index[index.size() - 1]; - - throw std::runtime_error( - fmt::format("Each element of the index ({}) must be non-negative.", err_msg_stream.str()) - ); - } - } - void CheckSelfIndexes() const { - for (const auto& index_and_value : index2value_) { - const auto& [index, value] = index_and_value; - CheckIndex(index); + // The last index is the lowest index of all index, + // so I only have to check if each of its elements is non-negative. + if ((std::prev(index2value_.end())->first < 0).any()) { + throw std::invalid_argument(fmt::format("Negative index not supported!")); } } - std::pair CheckIndexOfPairOfIterAndIsInserted( - const std::pair& iter_and_is_inserted - ) { - const auto& [iter, is_inserted] = iter_and_is_inserted; - if (is_inserted) { - CheckIndex(iter->first); - } - return iter_and_is_inserted; - } - - IndexContainer index2value_{ - {index_type::Zero(), 0} - }; + IndexContainer index2value_; }; -template < - std::signed_integral IntType, - class R, - int D, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -using DefaultMVPolynomial = - MVPolynomial, AllocatorOrContainer>; - -template < - std::signed_integral IntType, - class R, - int Dim, - class Comparer, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -auto D(const MVPolynomial& p, std::size_t axis) { - using MP = MVPolynomial; - - MP::CheckAxis(axis); - - auto new_index2value_seq = typename MP::sequence_type(p.get_allocator()); - new_index2value_seq.reserve(p.size()); - for (auto index_and_value : p) { - auto [index, value] = index_and_value; - if (index[axis] == 0) { - continue; - } else { - value *= index[axis]--; - new_index2value_seq.emplace_back(index, value); - } - } - const auto& comparer = p.key_comp(); - // Sort indexes in order not to violate the order by comparer. - std::sort( - new_index2value_seq.begin(), - new_index2value_seq.end(), - [&comparer](const MP::value_type& l, const MP::value_type& r) { - return comparer(l.first, r.first); - } - ); - return MP( - boost::container::ordered_unique_range_t(), - new_index2value_seq.begin(), - new_index2value_seq.end(), - comparer, - p.get_allocator() - ); -} +template +R MVPolynomial::tolerance = std::numeric_limits::min_exponent10; -template < - std::signed_integral IntType, - class R, - int Dim, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -auto D(const DefaultMVPolynomial& p, std::size_t axis) { - using MP = DefaultMVPolynomial; +template +auto D(const MVPolynomial& p, int axis) { + using MP = MVPolynomial; + using Index = typename MP::index_type; + using IndexAndCoeff = typename MP::value_type; - MP::CheckAxis(axis); + CheckAxis(MP::dim, axis); - auto new_index2value_seq = typename MP::sequence_type(p.get_allocator()); - new_index2value_seq.reserve(p.size()); + auto dp = MP{p.get_allocator()}; auto p_it = p.begin(); while (p_it != p.end()) { - auto [index, value] = *p_it; + auto value = p_it->second; + auto index = Index{p_it->first}; if (index[axis] == 0) { - if (axis == MP::dim - 1) { - ++p_it; - } auto d_end_it = p.end(); - for (std::size_t ith_axis = 0; ith_axis <= axis; ++ith_axis) { - d_end_it = std::partition_point( - p_it, - d_end_it, - [ith_axis, &index](const typename MP::value_type& v) { - return v.first[ith_axis] == index[ith_axis]; - } - ); + for (int ith_axis = 0; ith_axis <= axis; ++ith_axis) { + d_end_it = std::partition_point(p_it, d_end_it, [ith_axis, &index](const IndexAndCoeff& v) { + return v.first[ith_axis] == index[ith_axis]; + }); } // Skip indexes which axis-th element is zero. p_it = d_end_it; } else { value *= index[axis]--; - new_index2value_seq.emplace_back(index, value); + dp[index] = value; ++p_it; } } - return MP( - boost::container::ordered_unique_range_t(), - new_index2value_seq.begin(), - new_index2value_seq.end(), - p.key_comp(), - p.get_allocator() - ); + return dp; } -template < - std::signed_integral IntType, - class R, - int D, - class Comparer = IndexComparer, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -auto Integrate(MVPolynomial&& p, std::size_t axis) { - using MP = MVPolynomial; - - MP::CheckAxis(axis); - - auto index2value = p.extract_sequence(); - for (auto& index_and_value : index2value) { - auto& [index, value] = index_and_value; - value /= ++index[axis]; - } - const auto& comparer = p.key_comp(); - std::sort( - index2value.begin(), - index2value.end(), - [&comparer](const typename MP::value_type& l, const typename MP::value_type& r) { - return comparer(l.first, r.first); - } - ); - p.adopt_sequence(std::move(index2value)); - return std::move(p); -} - -template < - std::signed_integral IntType, - class R, - int D, - class Comparer = IndexComparer, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -auto Integrate( - const MVPolynomial& p, std::size_t axis -) { - return Integrate(MVPolynomial(p), axis); -} +template +auto Integrate(MVPolynomial&& p, int axis) { + using MP = MVPolynomial; -template < - std::signed_integral IntType, - class R, - int D, - class Comparer, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -auto Of( - const MVPolynomial& p, - const typename MVPolynomial::coord_type& x -) { - using MP = MVPolynomial; - typename MP::mapped_type sum = 0; - for (const auto& index_and_value : p) { - const auto& [index, value] = index_and_value; - sum += value * (x.array().pow(index.template cast())).prod(); - } - return sum; -} + CheckAxis(D, axis); -template -auto OfImpl(Iterator begin, Iterator end, int dim, std::size_t axis, const Coord& x) { - assert(axis >= 0 && axis < dim); - if (axis == dim - 1) { - auto [last_index, last_coeff] = *begin; - for (auto it = std::next(begin); it != end; ++it) { - const auto& [next_index, next_coeff] = *it; - last_coeff *= std::pow(x[axis], last_index[axis] - next_index[axis]); - last_coeff += next_coeff; - last_index = next_index; - } - last_coeff *= x.pow(last_index.template cast()).prod(); - return last_coeff; - } else { - auto sum = typename Iterator::value_type::second_type(0); - while (true) { - const auto& [first_index, first_coeff] = *begin; - auto partition_point = std::partition_point( - begin, - end, - [axis, &first_index](const typename Iterator::value_type& pair) { - return pair.first[axis] == first_index[axis]; - } - ); - const auto& [p_index, p_coeff] = *partition_point; - sum += OfImpl(begin, partition_point, dim, axis + 1, x); - if (partition_point == end) { - // The calculation ends. - break; - } - begin = partition_point; - } - return sum; + for (auto& index_and_value : p) { + auto& value = index_and_value.second; + auto& index = const_cast(index_and_value.first); + value /= ++index[axis]; } + return std::move(p); } -template < - std::signed_integral IntType, - class R, - int D, - class AllocatorOrContainer = - boost::container::new_allocator, R>>> -auto Of( - const DefaultMVPolynomial& p, - const typename DefaultMVPolynomial::coord_type& x -) { - using MP = DefaultMVPolynomial; - return OfImpl(p.cbegin(), p.cend(), MP::dim, 0, x); -} - -template