From 31ece42c198855a50d257afd266965c4984914d7 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Thu, 25 Dec 2025 19:55:18 +0900 Subject: [PATCH 1/8] Add clear func. --- include/mvPolynomial/mvPolynomial.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index 14a55dd..6e24d84 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -198,6 +198,8 @@ class MVPolynomial final { void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); } + void clear() { index2value_.clear(); } + // 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)); } @@ -375,8 +377,7 @@ class MVPolynomial final { auto comparer = l.key_comp(); auto mul = MVPolynomial(l.get_allocator()); - // Clear - mul.erase(mul.begin()); + mul.clear(); // 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; From 789970d2cefc5c1b2adf7a21b46bb39d01b8de45 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Thu, 25 Dec 2025 19:55:49 +0900 Subject: [PATCH 2/8] Explicitly use an allocator. --- include/mvPolynomial/mvPolynomial.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index 6e24d84..a017c0c 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -264,7 +264,7 @@ class MVPolynomial final { MVPolynomial operator+() const { return *this; } MVPolynomial operator-() const& { - auto m = MVPolynomial(*this); + auto m = MVPolynomial(*this, get_allocator()); for (auto& i_and_v : m) { auto& [_, v] = i_and_v; v = -v; @@ -344,7 +344,7 @@ class MVPolynomial final { friend bool operator!=(const MVPolynomial& l, const MVPolynomial& r) { return !(l == r); } friend MVPolynomial operator+(const MVPolynomial& l, const MVPolynomial& r) { - return MVPolynomial(l) + r; + return MVPolynomial(l, l.get_allocator()) + r; } friend MVPolynomial operator+(MVPolynomial&& l, const MVPolynomial& r) { @@ -359,7 +359,7 @@ class MVPolynomial final { friend MVPolynomial operator+(MVPolynomial&& l, MVPolynomial&& r) { return std::move(l) + r; } friend MVPolynomial operator-(const MVPolynomial& l, const MVPolynomial& r) { - return MVPolynomial(l) - r; + return MVPolynomial(l, l.get_allocator()) - r; } friend MVPolynomial operator-(MVPolynomial&& l, const MVPolynomial& r) { @@ -458,7 +458,7 @@ auto Integrate(MVPolynomial&& p, int axis) { template auto Integrate(const MVPolynomial& p, int axis) { - return Integrate(MVPolynomial(p), axis); + return Integrate(MVPolynomial(p, p.get_allocator()), axis); } } // namespace mvPolynomial From 981182a93f6e0470470b4d152189ae21f90ee809 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Thu, 25 Dec 2025 19:57:41 +0900 Subject: [PATCH 3/8] Add compound assignment ops for a scalar. --- include/mvPolynomial/mvPolynomial.hpp | 30 +++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index a017c0c..9c93e3c 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -280,6 +280,16 @@ class MVPolynomial final { return std::move(*this); } + MVPolynomial& operator+=(mapped_type r) { + Index idx = Index::Zero(); + if (contains(idx)) { + (*this)[idx] += r; + } else { + (*this)[idx] = r; + } + return *this; + } + MVPolynomial& operator+=(const MVPolynomial& r) { for (const auto& [idx, coeff] : r) { if (contains(idx)) { @@ -291,6 +301,16 @@ class MVPolynomial final { return *this; } + MVPolynomial& operator-=(mapped_type r) { + Index idx = Index::Zero(); + if (contains(idx)) { + (*this)[idx] -= r; + } else { + (*this)[idx] = r; + } + return *this; + } + MVPolynomial& operator-=(const MVPolynomial& r) { for (const auto& [idx, coeff] : r) { if (contains(idx)) { @@ -302,6 +322,16 @@ class MVPolynomial final { return *this; } + MVPolynomial& operator*=(mapped_type r) { + Index idx = Index::Zero(); + if (contains(idx)) { + (*this)[idx] *= r; + } else { + (*this)[idx] = r; + } + return *this; + } + MVPolynomial& operator*=(const MVPolynomial& r) { if (r.size() == 1) { const auto& [r_index, r_coeff] = *(r.begin()); From fc94df2d782af6df71a9caf902279d48905ecdc5 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Thu, 25 Dec 2025 23:18:23 +0900 Subject: [PATCH 4/8] Add operator()(const MVPolynomial&). --- include/mvPolynomial/mvPolynomial.hpp | 82 +++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index 9c93e3c..b1abc9e 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -261,6 +261,15 @@ class MVPolynomial final { R operator()(const coord_type& x) const { return details::OfImpl(crbegin(), crend(), dim, 0, x); } + MVPolynomial operator()(const MVPolynomial& x, int axis) const { + CheckAxis(dim, axis); + + auto composed_mvp = MVPolynomial{get_allocator()}; + composed_mvp.clear(); + OfImpl(composed_mvp, begin(), end(), 0, axis, x); + return composed_mvp; + } + MVPolynomial operator+() const { return *this; } MVPolynomial operator-() const& { @@ -435,6 +444,79 @@ class MVPolynomial final { } } + void OfImpl(MVPolynomial& composed_mvp, const_iterator begin, const_iterator end, int i, int axis, const MVPolynomial& x) const { + using Index = std::remove_cvref_t; + + CheckAxis(dim, i); + + if (i == axis) { + if (begin->first[i] == 0) { + for (auto it = begin; it != end; ++it) { + composed_mvp.insert(*it); + } + return; + } + if (axis == dim - 1) { + auto last_mvp = MVPolynomial{get_allocator()}; + last_mvp += 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_mvp *= x.pow(last_index[i] - next_index[i]); + last_mvp += next_coeff; + last_index = next_index; + } + last_coeff *= MVPolynomial{ + {{last_index, 1}}, get_allocator() + }; + composed_mvp += last_coeff; + return; + } + auto mvp = MVPolynomial{get_allocator()}; + 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]; + } + ); + for (auto it = begin; it != partition_point; ++it) { + auto index = Index{it->first}; + // set the index of axis to 0 for multiplying powed `x` all at once. + index[axis] = 0; + mvp.insert(std::make_pair(index, it->second)); + } + composed_mvp += mvp * x.pow(first_index[axis]); + + // The calculation ends. + if (partition_point == end) { + break; + } + begin = partition_point; + mvp.clear(); + } + } else { + while (true) { + const auto& [first_index, first_coeff] = *begin; + auto partition_point = std::partition_point( + begin, + end, + [i, &first_index](const typename Iterator::value_type& pair) { + return pair.first[i] == first_index[i]; + } + ); + OfImpl(composed_mvp, begin, partition_point, i + 1, axis, x); + if (partition_point == end) { + // The calculation ends. + break; + } + begin = partition_point; + } + } + } + IndexContainer index2value_; }; From 09496900cc93b7e4102bc3f5087d1b4da3ff0a22 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Mon, 29 Dec 2025 16:44:41 +0900 Subject: [PATCH 5/8] Add `pow` member func. --- include/mvPolynomial/mvPolynomial.hpp | 54 ++++++++++++++++++++++++--- test/mvPolynomial_test.cpp | 22 +++++++++++ 2 files changed, 71 insertions(+), 5 deletions(-) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index b1abc9e..cd5a0cd 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -5,7 +5,10 @@ #include "mvPolynomial/index_comparer.hpp" #include +#include +#include #include +#include #include #include "Eigen/Core" @@ -259,6 +262,42 @@ class MVPolynomial final { return index2value_.upper_bound(i); } + MVPolynomial pow(int exp) const { + if (exp < 0) { + throw std::invalid_argument("Given exp must be positive."); + } + switch (exp) { + case 0: + return MVPolynomial{1, get_allocator()}; + case 1: + return *this; + case 2: + return (*this) * (*this); + default: + auto max_pow2_under_exp = std::bit_floor(static_cast(exp)); + auto max_bit_width = std::bit_width(max_pow2_under_exp); + auto cache = std::vector(max_bit_width - 1, get_allocator()); + cache.at(0) = (*this) * (*this); + for (int i = 2; i < max_bit_width; ++i) { + cache.at(i - 1) = cache.at(i - 2) * cache.at(i - 2); + } + + assert(!cache.empty()); + auto powed_mvp = std::move(cache.back()); + exp -= max_pow2_under_exp; + while (exp > 1) { + auto max_pow2_under_exp = std::bit_floor(static_cast(exp)); + auto max_bit_width = std::bit_width(max_pow2_under_exp); + powed_mvp *= cache[max_bit_width - 2]; + exp -= max_pow2_under_exp; + } + if (exp == 1) { + powed_mvp *= *this; + } + return powed_mvp; + } + } + R operator()(const coord_type& x) const { return details::OfImpl(crbegin(), crend(), dim, 0, x); } MVPolynomial operator()(const MVPolynomial& x, int axis) const { @@ -444,7 +483,14 @@ class MVPolynomial final { } } - void OfImpl(MVPolynomial& composed_mvp, const_iterator begin, const_iterator end, int i, int axis, const MVPolynomial& x) const { + void OfImpl( + MVPolynomial& composed_mvp, + const_iterator begin, + const_iterator end, + int i, + int axis, + const MVPolynomial& x + ) const { using Index = std::remove_cvref_t; CheckAxis(dim, i); @@ -466,16 +512,14 @@ class MVPolynomial final { last_mvp += next_coeff; last_index = next_index; } - last_coeff *= MVPolynomial{ - {{last_index, 1}}, get_allocator() - }; + last_coeff *= MVPolynomial{{{last_index, 1}}, get_allocator()}; composed_mvp += last_coeff; return; } auto mvp = MVPolynomial{get_allocator()}; while (true) { const auto& [first_index, first_coeff] = *begin; - auto partition_point = std::partition_point( + auto partition_point = std::partition_point( begin, end, [axis, &first_index](const typename Iterator::value_type& pair) { diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp index 0aa3175..3806204 100644 --- a/test/mvPolynomial_test.cpp +++ b/test/mvPolynomial_test.cpp @@ -65,6 +65,28 @@ TEST_CASE("constructor", "[mvPolynomial]") { } } +TEST_CASE("pow", "[mvPolynomial]") { + auto m = MP2({ + {{0, 0}, 1}, + {{1, 0}, 2}, + {{0, 1}, 3}, + }); + + auto m3 = m * m * m; + + SECTION("0") { REQUIRE(m.pow(0) == MP2(1)); } + + SECTION("1") { REQUIRE(m.pow(1) == m); } + + SECTION("2") { REQUIRE(m.pow(2) == m * m); } + + SECTION("3") { REQUIRE(m.pow(3) == m * m * m); } + + SECTION("7") { REQUIRE(m.pow(7) == m3 * m3 * m); } + + SECTION("15") { REQUIRE(m.pow(15) == m3 * m3 * m3 * m3 * m3); } +} + TEST_CASE("operator()", "[mvPolynomial]") { auto m = MP2({ {{0, 0}, 1}, From 760d997a4c63adb9fc05b092444fe564dbb8f069 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Mon, 29 Dec 2025 16:54:06 +0900 Subject: [PATCH 6/8] Fix bugs. --- include/mvPolynomial/mvPolynomial.hpp | 44 +++++++++++++++------------ 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index cd5a0cd..fe12eb6 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -301,11 +301,11 @@ class MVPolynomial final { R operator()(const coord_type& x) const { return details::OfImpl(crbegin(), crend(), dim, 0, x); } MVPolynomial operator()(const MVPolynomial& x, int axis) const { - CheckAxis(dim, axis); + details::CheckAxis(dim, axis); auto composed_mvp = MVPolynomial{get_allocator()}; composed_mvp.clear(); - OfImpl(composed_mvp, begin(), end(), 0, axis, x); + OfImpl(composed_mvp, rbegin(), rend(), 0, axis, x); return composed_mvp; } @@ -329,7 +329,7 @@ class MVPolynomial final { } MVPolynomial& operator+=(mapped_type r) { - Index idx = Index::Zero(); + auto idx = index_type::Zero(); if (contains(idx)) { (*this)[idx] += r; } else { @@ -350,7 +350,7 @@ class MVPolynomial final { } MVPolynomial& operator-=(mapped_type r) { - Index idx = Index::Zero(); + auto idx = index_type::Zero(); if (contains(idx)) { (*this)[idx] -= r; } else { @@ -371,7 +371,7 @@ class MVPolynomial final { } MVPolynomial& operator*=(mapped_type r) { - Index idx = Index::Zero(); + auto idx = index_type::Zero(); if (contains(idx)) { (*this)[idx] *= r; } else { @@ -484,21 +484,25 @@ class MVPolynomial final { } void OfImpl( - MVPolynomial& composed_mvp, - const_iterator begin, - const_iterator end, - int i, - int axis, - const MVPolynomial& x + MVPolynomial& composed_mvp, + const_reverse_iterator begin, + const_reverse_iterator end, + int i, + int axis, + const MVPolynomial& x ) const { - using Index = std::remove_cvref_t; + using Index = std::remove_cvref_t; - CheckAxis(dim, i); + details::CheckAxis(dim, i); if (i == axis) { if (begin->first[i] == 0) { for (auto it = begin; it != end; ++it) { - composed_mvp.insert(*it); + if (composed_mvp.contains(it->first)) { + composed_mvp.at(it->first) += it->second; + } else { + composed_mvp.insert(*it); + } } return; } @@ -512,17 +516,18 @@ class MVPolynomial final { last_mvp += next_coeff; last_index = next_index; } - last_coeff *= MVPolynomial{{{last_index, 1}}, get_allocator()}; - composed_mvp += last_coeff; + last_mvp *= MVPolynomial{{{last_index, 1}}, get_allocator()}; + composed_mvp += last_mvp; return; } auto mvp = MVPolynomial{get_allocator()}; + mvp.clear(); 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) { + [axis, &first_index](const typename const_reverse_iterator::value_type& pair) { return pair.first[axis] == first_index[axis]; } ); @@ -547,7 +552,7 @@ class MVPolynomial final { auto partition_point = std::partition_point( begin, end, - [i, &first_index](const typename Iterator::value_type& pair) { + [i, &first_index](const typename const_reverse_iterator::value_type& pair) { return pair.first[i] == first_index[i]; } ); @@ -565,7 +570,8 @@ class MVPolynomial final { }; template -R MVPolynomial::tolerance = std::numeric_limits::min_exponent10; +R MVPolynomial::tolerance = + std::ldexp(std::numeric_limits::epsilon(), std::numeric_limits::min_exponent); template auto D(const MVPolynomial& p, int axis) { From 061a076d11eb2d3bc0b35cb5ec7e574dd17d909c Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Mon, 29 Dec 2025 16:55:58 +0900 Subject: [PATCH 7/8] Add allocator arg to the ctr from scalar. --- include/mvPolynomial/mvPolynomial.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index fe12eb6..2c993f1 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -11,7 +11,6 @@ #include #include -#include "Eigen/Core" #include "fmt/core.h" #include "platanus/btree_map.hpp" @@ -157,7 +156,7 @@ class MVPolynomial final { return *this; } - MVPolynomial(mapped_type r) { index2value_.at(index_type::Zero()) = r; } + MVPolynomial(mapped_type r, const allocator_type& a = allocator_type{}) : index2value_({{index_type::Zero(), r}}, a) {} allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); } @@ -479,7 +478,7 @@ class MVPolynomial final { // The first index is the lowest index of all index, // so I only have to check if each of its elements is non-negative. if ((index2value_.begin()->first < 0).any()) { - throw std::invalid_argument(fmt::format("Negative index not supported!")); + throw std::invalid_argument("Negative index not supported!"); } } @@ -507,8 +506,7 @@ class MVPolynomial final { return; } if (axis == dim - 1) { - auto last_mvp = MVPolynomial{get_allocator()}; - last_mvp += begin->second; + auto last_mvp = MVPolynomial{begin->second, get_allocator()}; auto last_index = Index{begin->first}; for (auto it = std::next(begin); it != end; ++it) { const auto& [next_index, next_coeff] = *it; From 15f83ee96f80f1284794eda7ce57a6caa8f85b95 Mon Sep 17 00:00:00 2001 From: Yuya Asano <64895419+sukeya@users.noreply.github.com> Date: Mon, 29 Dec 2025 16:57:06 +0900 Subject: [PATCH 8/8] Add tests for assigning mvPoly. --- test/mvPolynomial_test.cpp | 98 ++++++++++++++++++++++++++++++++++---- 1 file changed, 88 insertions(+), 10 deletions(-) diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp index 3806204..4db4ad0 100644 --- a/test/mvPolynomial_test.cpp +++ b/test/mvPolynomial_test.cpp @@ -88,17 +88,95 @@ TEST_CASE("pow", "[mvPolynomial]") { } TEST_CASE("operator()", "[mvPolynomial]") { - auto m = MP2({ - {{0, 0}, 1}, - {{1, 0}, 2}, - {{0, 1}, 3}, - {{1, 1}, 4}, - {{2, 0}, 5}, - {{0, 2}, 6}, - }); + SECTION("point") { + auto m = MP2({ + {{0, 0}, 1}, + {{1, 0}, 2}, + {{0, 1}, 3}, + {{1, 1}, 4}, + {{2, 0}, 5}, + {{0, 2}, 6}, + }); + + REQUIRE(m(Eigen::Vector2d::Zero()) == 1); + REQUIRE(m(Eigen::Vector2d({2, 3})) == 112); + } - REQUIRE(m(Eigen::Vector2d::Zero()) == 1); - REQUIRE(m(Eigen::Vector2d({2, 3})) == 112); + SECTION("mvpolynomial_1d") { + auto m = MP3({ + {{0, 0, 0}, 1}, + {{1, 0, 0}, 2}, + {{0, 1, 0}, 3}, + {{0, 0, 1}, 4}, + }); + auto x = MP3({ + {{0, 0, 0}, 1}, + {{1, 0, 0}, 2}, + {{0, 1, 0}, 3}, + }); + REQUIRE( + m(x, 2) + == MP3({ + {{0, 0, 0}, 5}, + {{1, 0, 0}, 10}, + {{0, 1, 0}, 15}, + }) + ); + auto y = MP3({ + {{0, 0, 0}, 1}, + {{1, 0, 0}, 2}, + {{0, 0, 1}, 4}, + }); + REQUIRE( + m(y, 1) + == MP3({ + {{0, 0, 0}, 4}, + {{1, 0, 0}, 8}, + {{0, 0, 1}, 16}, + }) + ); + } + + SECTION("mvpolynomial_2d") { + auto m = MP3({ + {{0, 0, 0}, 1}, + {{2, 0, 0}, 2}, + {{0, 2, 0}, 3}, + {{0, 0, 2}, 4}, + }); + auto x = MP3({ + {{0, 0, 0}, 1}, + {{1, 0, 0}, 2}, + {{0, 1, 0}, 3}, + }); + REQUIRE( + m(x, 2) + == MP3({ + {{0, 0, 0}, 5}, + {{1, 0, 0}, 16}, + {{0, 1, 0}, 24}, + {{1, 1, 0}, 48}, + {{2, 0, 0}, 18}, + {{0, 2, 0}, 39}, + }) + ); + auto y = MP3({ + {{0, 0, 0}, 1}, + {{1, 0, 0}, 2}, + {{0, 0, 1}, 4}, + }); + REQUIRE( + m(y, 1) + == MP3({ + {{0, 0, 0}, 4}, + {{1, 0, 0}, 12}, + {{0, 0, 1}, 24}, + {{1, 0, 1}, 48}, + {{2, 0, 0}, 14}, + {{0, 0, 2}, 52}, + }) + ); + } } TEST_CASE("D", "[mvPolynomial]") {