From 392d8d7773e987ee2eb1d15a7fd57d94197350eb Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sat, 16 May 2026 15:36:17 +0900
Subject: [PATCH 01/41] Remove an unused member func.

---
 include/mvPolynomial/index_comparer.hpp | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/include/mvPolynomial/index_comparer.hpp b/include/mvPolynomial/index_comparer.hpp
index 78ecf96..2cd1a1a 100644
--- a/include/mvPolynomial/index_comparer.hpp
+++ b/include/mvPolynomial/index_comparer.hpp
@@ -23,24 +23,6 @@ class IndexComparer final {
 
   static constexpr int dim = D;
 
-  /**
-   * \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
-   */
-  static constexpr std::strong_ordering Compare(const Index& lhd, const Index& rhd) {
-    for (int i = 0; i < lhd.size(); ++i) {
-      auto comp = lhd[i] <=> rhd[i];
-      if (comp < 0) {
-        return std::strong_ordering::greater;
-      } else if (comp > 0) {
-        return std::strong_ordering::less;
-      }
-    }
-    return std::strong_ordering::equal;
-  }
-
   /**
    * \brief If lhd[i] == rhd[i] for i = 0, ..., N - 1 and lhd[N] < rhd[N],
    * return true: otherwise, false.

From 88ef36dca5be00161af79431b3ee979a43a0b350 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sat, 16 May 2026 15:37:27 +0900
Subject: [PATCH 02/41] Add clang-format to pre-commit setting.

---
 .github/workflows/pre_commit.yml | 2 +-
 .pre-commit-config.yaml          | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/pre_commit.yml b/.github/workflows/pre_commit.yml
index aa92e2f..51212ab 100644
--- a/.github/workflows/pre_commit.yml
+++ b/.github/workflows/pre_commit.yml
@@ -16,6 +16,6 @@ jobs:
       - name: Set up pre-commit
         run: |
           python -m pip install --upgrade pip
-          pip install pre-commit
+          pip install pre-commit clang-format
       - name: Run pre-commit for all files
         run: pre-commit run --all-files
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index c3a3919..53a3b9f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -8,7 +8,7 @@ repos:
     -   id: check-yaml
     -   id: check-added-large-files
 -   repo: https://github.com/pre-commit/mirrors-clang-format
-    rev: v21.1.0
+    rev: v22.1.5
     hooks:
     -   id: clang-format
         types_or: [c++]

From d963694d59977186e5b1057266f95d9eae8a7c6c Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sat, 16 May 2026 23:25:18 +0900
Subject: [PATCH 03/41] Fix format.

---
 include/mvPolynomial/mvPolynomial.hpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index c574bba..a8cb53a 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -160,8 +160,7 @@ class MVPolynomial final {
   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(); }
-
-  Comparer key_comp() const noexcept { return index2value_.key_comp(); }
+  key_compare    key_comp() const noexcept { return index2value_.key_comp(); }
 
   iterator       begin() noexcept { return index2value_.begin(); }
   const_iterator begin() const noexcept { return index2value_.begin(); }
@@ -181,7 +180,8 @@ class MVPolynomial final {
   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(); }

From edb48c7b0568b301f6af2638853c27b2ebb09085 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sat, 16 May 2026 23:25:35 +0900
Subject: [PATCH 04/41] Simplyfy Integrate.

---
 include/mvPolynomial/mvPolynomial.hpp | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index a8cb53a..87f9263 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -596,7 +596,7 @@ auto D(const MVPolynomial<IntType, R, Dim, Allocator>& p, int axis) {
 }
 
 template <std::signed_integral IntType, std::floating_point R, int D, class Allocator>
-auto Integrate(MVPolynomial<IntType, R, D, Allocator>&& p, int axis) {
+auto Integrate(MVPolynomial<IntType, R, D, Allocator> p, int axis) {
   using MP = MVPolynomial<IntType, R, D, Allocator>;
 
   details::CheckAxis(D, axis);
@@ -609,11 +609,6 @@ auto Integrate(MVPolynomial<IntType, R, D, Allocator>&& p, int axis) {
   return std::move(p);
 }
 
-template <std::signed_integral IntType, std::floating_point R, int D, class Allocator>
-auto Integrate(const MVPolynomial<IntType, R, D, Allocator>& p, int axis) {
-  return Integrate(MVPolynomial<IntType, R, D, Allocator>(p, p.get_allocator()), axis);
-}
-
 }  // namespace mvPolynomial
 
 #endif

From cd53f864c74c88b8135a97720980eeaaea0fc984 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 00:07:44 +0900
Subject: [PATCH 05/41] Simplify the impl of operator().

---
 include/mvPolynomial/mvPolynomial.hpp | 150 ++++----------------------
 1 file changed, 22 insertions(+), 128 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 87f9263..55b2e25 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -23,45 +23,6 @@ inline void CheckAxis(int dim, int axis) {
     );
   }
 }
-
-template <class Iterator, class Coord>
-auto OfImpl(Iterator begin, Iterator end, int dim, int axis, const Coord& x) {
-  using Index = std::remove_cvref_t<typename Iterator::value_type::first_type>;
-
-  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<typename Coord::Scalar>()).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];
-          }
-      );
-      sum += OfImpl(begin, partition_point, dim, axis + 1, x);
-      if (partition_point == end) {
-        // The calculation ends.
-        break;
-      }
-      begin = partition_point;
-    }
-    return sum;
-  }
-}
 }  // namespace details
 
 template <
@@ -291,14 +252,29 @@ class MVPolynomial final {
     }
   }
 
-  R operator()(const coord_type& x) const { return details::OfImpl(crbegin(), crend(), dim, 0, x); }
+  R operator()(const coord_type& x) const {
+    auto sum = R(0.0);
+    for (const auto& index_and_coeff : index2value_) {
+      auto        coeff = index_and_coeff.second;
+      const auto& index = index_and_coeff.first;
+      sum += coeff * (x.array().pow(index.template cast<double>())).prod();
+    }
+    return sum;
+  }
 
-  MVPolynomial operator()(const MVPolynomial& x, int axis) const {
+  MVPolynomial operator()(const MVPolynomial& mvp, int axis) const {
     details::CheckAxis(dim, axis);
 
     auto composed_mvp = MVPolynomial{get_allocator()};
-    composed_mvp.clear();
-    OfImpl(composed_mvp, rbegin(), rend(), 0, axis, x);
+    for (const auto& index_and_coeff : index2value_) {
+      auto        coeff = index_and_coeff.second;
+      const auto& index = index_and_coeff.first;
+
+      auto index_with_axis_0  = index;
+      index_with_axis_0[axis] = 0;
+      composed_mvp +=
+          MVPolynomial{{{index_with_axis_0, coeff}}, get_allocator()} * mvp.pow(index[axis]);
+    }
     return composed_mvp;
   }
 
@@ -468,91 +444,9 @@ class MVPolynomial final {
 
  private:
   void CheckSelfIndexes() const {
-    // 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("Negative index not supported!");
-    }
-  }
-
-  void OfImpl(
-      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<typename const_iterator::value_type::first_type>;
-
-    details::CheckAxis(dim, i);
-
-    if (i == axis) {
-      if (begin->first[i] == 0) {
-        for (auto it = begin; it != end; ++it) {
-          if (composed_mvp.contains(it->first)) {
-            composed_mvp.at(it->first) += it->second;
-          } else {
-            composed_mvp.insert(*it);
-          }
-        }
-        return;
-      }
-      if (axis == dim - 1) {
-        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;
-          last_mvp *= x.pow(last_index[i] - next_index[i]);
-          last_mvp += next_coeff;
-          last_index = next_index;
-        }
-        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 const_reverse_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 const_reverse_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;
+    for (const auto& index_and_coeff : index2value_) {
+      if ((index_and_coeff.first < 0).any()) {
+        throw std::invalid_argument("Negative index not supported!");
       }
     }
   }

From 7a7daf094292768caa5db563a401e57c24b451a5 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 00:09:09 +0900
Subject: [PATCH 06/41] Remove deleting funcs.

---
 include/mvPolynomial/mvPolynomial.hpp | 24 +-----------------------
 1 file changed, 1 insertion(+), 23 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 55b2e25..5d5c844 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -155,28 +155,6 @@ 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<iterator, bool> insert(const value_type& x) { return index2value_.insert(x); }
-  std::pair<iterator, bool> insert(value_type&& x) { return index2value_.insert(std::move(x)); }
-
-  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 <typename InputIterator>
-  void insert(InputIterator b, InputIterator e) {
-    index2value_.insert(b, e);
-  }
-
-  void insert(std::initializer_list<value_type> list) {
-    index2value_.insert(list.begin(), list.end());
-  }
-
-  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); }
   template <typename K>
@@ -423,7 +401,7 @@ class MVPolynomial final {
     auto comparer = l.key_comp();
 
     auto mul = MVPolynomial(l.get_allocator());
-    mul.clear();
+    mul.index2value_.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 22175af9aa8a418a57e7f8ae048031ef19702012 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 00:09:45 +0900
Subject: [PATCH 07/41] Remove unused alias type.

---
 include/mvPolynomial/mvPolynomial.hpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 5d5c844..1c0cff5 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -43,8 +43,7 @@ class MVPolynomial final {
   using coord_type = CoordType<R, dim>;
 
  private:
-  using Comparer       = IndexComparer<IntType, D>;
-  using IndexContainer = platanus::btree_map<index_type, R, Comparer, Allocator>;
+  using IndexContainer = platanus::btree_map<index_type, R, IndexComparer<IntType, D>, Allocator>;
 
  public:
   using key_type    = IndexContainer::key_type;

From c5c801b3d5e68dd9c894e1bc1f12769087cabecd Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:10:30 +0900
Subject: [PATCH 08/41] Simplify .clang-format.

---
 .clang-format | 219 ++------------------------------------------------
 1 file changed, 7 insertions(+), 212 deletions(-)

diff --git a/.clang-format b/.clang-format
index 42591e6..496efca 100644
--- a/.clang-format
+++ b/.clang-format
@@ -1,225 +1,20 @@
 ---
-Language:        Cpp
-# BasedOnStyle:  Google
-AccessModifierOffset: -1
+Language: Cpp
+BasedOnStyle: Google
 AlignAfterOpenBracket: BlockIndent
 AlignArrayOfStructures: Right
 AlignConsecutiveMacros: Consecutive
 AlignConsecutiveAssignments: Consecutive
 AlignConsecutiveBitFields: Consecutive
 AlignConsecutiveDeclarations: Consecutive
-AlignEscapedNewlines: Left
-AlignOperands:   Align
-AlignTrailingComments: true
 AllowAllArgumentsOnNextLine: false
-AllowAllParametersOfDeclarationOnNextLine: true
-AllowShortEnumsOnASingleLine: true
-AllowShortBlocksOnASingleLine: Never
-AllowShortCaseLabelsOnASingleLine: false
-AllowShortFunctionsOnASingleLine: All
-AllowShortLambdasOnASingleLine: All
-AllowShortIfStatementsOnASingleLine: WithoutElse
-AllowShortLoopsOnASingleLine: true
-AlwaysBreakAfterDefinitionReturnType: None
-AlwaysBreakAfterReturnType: None
-AlwaysBreakBeforeMultilineStrings: true
-AlwaysBreakTemplateDeclarations: Yes
-AttributeMacros:
-  - __capability
 BinPackArguments: false
 BinPackParameters: false
-BraceWrapping:
-  AfterCaseLabel:  false
-  AfterClass:      false
-  AfterControlStatement: Never
-  AfterEnum:       false
-  AfterFunction:   false
-  AfterNamespace:  false
-  AfterObjCDeclaration: false
-  AfterStruct:     false
-  AfterUnion:      false
-  AfterExternBlock: false
-  BeforeCatch:     false
-  BeforeElse:      false
-  BeforeLambdaBody: false
-  BeforeWhile:     false
-  IndentBraces:    false
-  SplitEmptyFunction: true
-  SplitEmptyRecord: true
-  SplitEmptyNamespace: true
-BreakBeforeBinaryOperators: NonAssignment
-BreakBeforeConceptDeclarations: true
-BreakBeforeBraces: Attach
-BreakBeforeInheritanceComma: false
-BreakInheritanceList: BeforeColon
-BreakBeforeTernaryOperators: true
-BreakConstructorInitializersBeforeComma: false
-BreakConstructorInitializers: BeforeColon
-BreakAfterJavaFieldAnnotations: false
-BreakStringLiterals: true
-ColumnLimit:     100
-CommentPragmas:  '^ IWYU pragma:'
-QualifierAlignment: Leave
-CompactNamespaces: false
-ConstructorInitializerIndentWidth: 4
-ContinuationIndentWidth: 4
-Cpp11BracedListStyle: true
-DeriveLineEnding: true
-DerivePointerAlignment: false
-DisableFormat:   false
-EmptyLineAfterAccessModifier: Never
-EmptyLineBeforeAccessModifier: LogicalBlock
-ExperimentalAutoDetectBinPacking: false
+ColumnLimit: 120
+IndentRequires: false
+#PackArguments.BinPackArgumentsStyle: OnePerLine
 PackConstructorInitializers: NextLine
-BasedOnStyle:    ''
-ConstructorInitializerAllOnOneLineOrOnePerLine: false
-AllowAllConstructorInitializersOnNextLine: true
-FixNamespaceComments: true
-ForEachMacros:
-  - foreach
-  - Q_FOREACH
-  - BOOST_FOREACH
-IfMacros:
-  - KJ_IF_MAYBE
-IncludeBlocks:   Regroup
-IncludeCategories:
-  - Regex:           '^<ext/.*\.h>'
-    Priority:        2
-    SortPriority:    0
-    CaseSensitive:   false
-  - Regex:           '^<.*\.h>'
-    Priority:        1
-    SortPriority:    0
-    CaseSensitive:   false
-  - Regex:           '^<.*'
-    Priority:        2
-    SortPriority:    0
-    CaseSensitive:   false
-  - Regex:           '.*'
-    Priority:        3
-    SortPriority:    0
-    CaseSensitive:   false
-IncludeIsMainRegex: '([-_](test|unittest))?$'
-IncludeIsMainSourceRegex: ''
-IndentAccessModifiers: false
-IndentCaseLabels: true
-IndentCaseBlocks: false
-IndentGotoLabels: true
-IndentPPDirectives: None
-IndentExternBlock: AfterExternBlock
-IndentRequires:  false
-IndentWidth:     2
-IndentWrappedFunctionNames: false
-InsertTrailingCommas: None
-JavaScriptQuotes: Leave
-JavaScriptWrapImports: true
-KeepEmptyLinesAtTheStartOfBlocks: false
-LambdaBodyIndentation: Signature
-MacroBlockBegin: ''
-MacroBlockEnd:   ''
-MaxEmptyLinesToKeep: 1
-NamespaceIndentation: None
-ObjCBinPackProtocolList: Never
-ObjCBlockIndentWidth: 2
-ObjCBreakBeforeNestedBlockParam: true
-ObjCSpaceAfterProperty: false
-ObjCSpaceBeforeProtocolList: true
-PenaltyBreakAssignment: 2
-PenaltyBreakBeforeFirstCallParameter: 1
-PenaltyBreakComment: 300
-PenaltyBreakFirstLessLess: 120
-PenaltyBreakOpenParenthesis: 0
-PenaltyBreakString: 1000
-PenaltyBreakTemplateDeclaration: 10
-PenaltyExcessCharacter: 1000000
-PenaltyReturnTypeOnItsOwnLine: 200
-PenaltyIndentedWhitespace: 0
+#PackParameters.BinPackParametersStyle: OnePerLine
 PointerAlignment: Left
-PPIndentWidth:   -1
-RawStringFormats:
-  - Language:        Cpp
-    Delimiters:
-      - cc
-      - CC
-      - cpp
-      - Cpp
-      - CPP
-      - 'c++'
-      - 'C++'
-    CanonicalDelimiter: ''
-    BasedOnStyle:    google
-  - Language:        TextProto
-    Delimiters:
-      - pb
-      - PB
-      - proto
-      - PROTO
-    EnclosingFunctions:
-      - EqualsProto
-      - EquivToProto
-      - PARSE_PARTIAL_TEXT_PROTO
-      - PARSE_TEST_PROTO
-      - PARSE_TEXT_PROTO
-      - ParseTextOrDie
-      - ParseTextProtoOrDie
-      - ParseTestProto
-      - ParsePartialTestProto
-    CanonicalDelimiter: pb
-    BasedOnStyle:    google
-ReferenceAlignment: Pointer
-ReflowComments:  true
-RemoveBracesLLVM: false
-SeparateDefinitionBlocks: Leave
-ShortNamespaceLines: 1
-SortIncludes:    Never
-SortJavaStaticImport: Before
-SortUsingDeclarations: true
-SpaceAfterCStyleCast: false
-SpaceAfterLogicalNot: false
-SpaceAfterTemplateKeyword: true
-SpaceBeforeAssignmentOperators: true
-SpaceBeforeCaseColon: false
-SpaceBeforeCpp11BracedList: false
-SpaceBeforeCtorInitializerColon: true
-SpaceBeforeInheritanceColon: true
-SpaceBeforeParens: ControlStatements
-SpaceBeforeParensOptions:
-  AfterControlStatements: true
-  AfterForeachMacros: true
-  AfterFunctionDefinitionName: false
-  AfterFunctionDeclarationName: false
-  AfterIfMacros:   true
-  AfterOverloadedOperator: false
-  BeforeNonEmptyParentheses: false
-SpaceAroundPointerQualifiers: Default
-SpaceBeforeRangeBasedForLoopColon: true
-SpaceInEmptyBlock: false
-SpaceInEmptyParentheses: false
-SpacesBeforeTrailingComments: 2
-SpacesInAngles:  Never
-SpacesInConditionalStatement: false
-SpacesInContainerLiterals: true
-SpacesInCStyleCastParentheses: false
-SpacesInLineCommentPrefix:
-  Minimum:         1
-  Maximum:         -1
-SpacesInParentheses: false
-SpacesInSquareBrackets: false
-SpaceBeforeSquareBrackets: false
-BitFieldColonSpacing: Both
-Standard:        Auto
-StatementAttributeLikeMacros:
-  - Q_EMIT
-StatementMacros:
-  - Q_UNUSED
-  - QT_REQUIRE_VERSION
-TabWidth:        8
-UseCRLF:         false
-UseTab:          Never
-WhitespaceSensitiveMacros:
-  - STRINGIZE
-  - PP_STRINGIZE
-  - BOOST_PP_STRINGIZE
-  - NS_SWIFT_NAME
-  - CF_SWIFT_NAME
+SortIncludes: Never
 ...

From 9f0467517f01f999bb9122d1bdbcc832f1383f4b Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:11:00 +0900
Subject: [PATCH 09/41] Format.

---
 include/mvPolynomial/mvPolynomial.hpp | 24 ++++----------
 test/mvPolynomial_test.cpp            | 48 ++++++++++++---------------
 2 files changed, 29 insertions(+), 43 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 1c0cff5..758c73f 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -18,9 +18,7 @@ namespace mvPolynomial {
 namespace details {
 inline 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)
-    );
+    throw std::runtime_error(fmt::format("CheckAxis: Given axis {} must be in [0, {}).", axis, dim));
   }
 }
 }  // namespace details
@@ -95,20 +93,17 @@ class MVPolynomial final {
   }
 
   template <typename InputIterator>
-  MVPolynomial(InputIterator s, InputIterator e, const allocator_type& allocator)
-      : index2value_(s, e, allocator) {
+  MVPolynomial(InputIterator s, InputIterator e, const allocator_type& allocator) : index2value_(s, e, allocator) {
     CheckSelfIndexes();
   }
 
-  MVPolynomial(std::initializer_list<value_type> l, const allocator_type& a = allocator_type{})
-      : index2value_(l, a) {
+  MVPolynomial(std::initializer_list<value_type> l, const allocator_type& a = allocator_type{}) : index2value_(l, a) {
     CheckSelfIndexes();
   }
 
   MVPolynomial(const MVPolynomial& m, const allocator_type& a) : index2value_(m.index2value_, a) {}
 
-  MVPolynomial(MVPolynomial&& m, const allocator_type& a)
-      : index2value_(std::move(m.index2value_), a) {}
+  MVPolynomial(MVPolynomial&& m, const allocator_type& a) : index2value_(std::move(m.index2value_), a) {}
 
   MVPolynomial& operator=(std::initializer_list<value_type> l) {
     index2value_.clear();
@@ -249,8 +244,7 @@ class MVPolynomial final {
 
       auto index_with_axis_0  = index;
       index_with_axis_0[axis] = 0;
-      composed_mvp +=
-          MVPolynomial{{{index_with_axis_0, coeff}}, get_allocator()} * mvp.pow(index[axis]);
+      composed_mvp += MVPolynomial{{{index_with_axis_0, coeff}}, get_allocator()} * mvp.pow(index[axis]);
     }
     return composed_mvp;
   }
@@ -375,9 +369,7 @@ class MVPolynomial final {
     return std::move(l);
   }
 
-  friend MVPolynomial operator+(const MVPolynomial& l, MVPolynomial&& r) {
-    return std::move(r) + l;
-  }
+  friend MVPolynomial operator+(const MVPolynomial& l, MVPolynomial&& r) { return std::move(r) + l; }
 
   friend MVPolynomial operator+(MVPolynomial&& l, MVPolynomial&& r) { return std::move(l) + r; }
 
@@ -390,9 +382,7 @@ class MVPolynomial final {
     return std::move(l);
   }
 
-  friend MVPolynomial operator-(const MVPolynomial& l, MVPolynomial&& r) {
-    return -std::move(r) + l;
-  }
+  friend MVPolynomial operator-(const MVPolynomial& l, MVPolynomial&& r) { return -std::move(r) + l; }
 
   friend MVPolynomial operator-(MVPolynomial&& l, MVPolynomial&& r) { return std::move(l) - r; }
 
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 4db4ad0..88c40d3 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -115,11 +115,10 @@ TEST_CASE("operator()", "[mvPolynomial]") {
         {{0, 1, 0}, 3},
     });
     REQUIRE(
-        m(x, 2)
-        == MP3({
-            {{0, 0, 0},  5},
-            {{1, 0, 0}, 10},
-            {{0, 1, 0}, 15},
+        m(x, 2) == MP3({
+                       {{0, 0, 0},  5},
+                       {{1, 0, 0}, 10},
+                       {{0, 1, 0}, 15},
     })
     );
     auto y = MP3({
@@ -128,11 +127,10 @@ TEST_CASE("operator()", "[mvPolynomial]") {
         {{0, 0, 1}, 4},
     });
     REQUIRE(
-        m(y, 1)
-        == MP3({
-            {{0, 0, 0},  4},
-            {{1, 0, 0},  8},
-            {{0, 0, 1}, 16},
+        m(y, 1) == MP3({
+                       {{0, 0, 0},  4},
+                       {{1, 0, 0},  8},
+                       {{0, 0, 1}, 16},
     })
     );
   }
@@ -150,14 +148,13 @@ TEST_CASE("operator()", "[mvPolynomial]") {
         {{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},
+        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({
@@ -166,14 +163,13 @@ TEST_CASE("operator()", "[mvPolynomial]") {
         {{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},
+        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},
     })
     );
   }

From 99f68153bc2125c41168673b8839581a6a2822ba Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:15:31 +0900
Subject: [PATCH 10/41] Fix operator*=(scalar).

---
 include/mvPolynomial/mvPolynomial.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 758c73f..b550c78 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -314,7 +314,7 @@ class MVPolynomial final {
     auto idx = index_type::Zero();
     for (auto& index_and_coeff : *this) {
       auto& coeff = index_and_coeff.second;
-      coeff *= 2;
+      coeff *= r;
     }
     return *this;
   }

From eef1856df63690ba47ea0d61eee7df57236ef190 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:16:54 +0900
Subject: [PATCH 11/41] Fix operator-=(scalar).

---
 include/mvPolynomial/mvPolynomial.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index b550c78..10331c7 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -294,7 +294,7 @@ class MVPolynomial final {
     if (contains(idx)) {
       (*this)[idx] -= r;
     } else {
-      (*this)[idx] = r;
+      (*this)[idx] = -r;
     }
     return *this;
   }

From c50667ebca361239d6fb72c200b22e2b45ada727 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:19:16 +0900
Subject: [PATCH 12/41] Fix test cmake setting.

---
 test/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 449adcc..8f54824 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -9,7 +9,7 @@ FetchContent_MakeAvailable(Catch2)
 add_executable(mvPolynomial_test mvPolynomial_test.cpp)
 target_compile_options(
   mvPolynomial_test
-  INTERFACE
+  PRIVATE
     # Common options for gcc and clang
     "$<${gcc_like}:-Wall;-Wextra;-Werror>"
     "$<${msvc}:/Wall;/WX>"

From d7de6db9ad4c15d44445f19b777659f86c760989 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:20:55 +0900
Subject: [PATCH 13/41] Fix README.md.

---
 README.md | 14 --------------
 1 file changed, 14 deletions(-)

diff --git a/README.md b/README.md
index 1aef24c..61f4b02 100644
--- a/README.md
+++ b/README.md
@@ -33,17 +33,3 @@ The examples exist in "test" directory.
 
 ## Multi-variable polynomial
 A class `MVPolynomial` implements multi-variable polynomials.
-Its interface is the same as Boost's flat_map except that it lacks `merge` member function.
-So, please see the document of Boost's flat_map.
-
-A class `ExactOf` calculates f of x as less multiplication as possible, but uses more memories.
-
-## Polynomial product
-A class `PolynomialProduct` implements a product of polynomials which have different variable each other (example: (1 + x + x^2) * (y + 6 * y^3 + y^10) * (1 + z)).
-Its interface is the same as std::array.
-So, please see the document of std::array.
-
-## Polynomial
-A class `Polynomial` implements one-variable polynomials.
-Its interface is also the same as Boost's flat_map except that it lacks `merge` member function.
-So, please see the document of Boost's flat_map.

From 682a37aaef35affe42ce7e128eebfd0335c3eb7f Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 01:22:44 +0900
Subject: [PATCH 14/41] Add using std header.

---
 include/mvPolynomial/mvPolynomial.hpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 10331c7..9f7e7ca 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -6,10 +6,13 @@
 
 #include <algorithm>
 #include <bit>
+#include <cassert>
 #include <cmath>
 #include <iterator>
+#include <limits>
 #include <stdexcept>
 #include <type_traits>
+#include <vector>
 
 #include "fmt/core.h"
 #include "platanus/btree_map.hpp"

From fba80e946fcdce67e3bea9b732a625b38a40ec74 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 07:58:22 +0900
Subject: [PATCH 15/41] Fix.

---
 CMakeLists.txt                        | 2 +-
 include/mvPolynomial/mvPolynomial.hpp | 2 --
 test/mvPolynomial_test.cpp            | 4 ++--
 3 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index fe3e4be..45e274f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -63,7 +63,7 @@ set(gcc "$<COMPILE_LANG_AND_ID:CXX,GNU>")
 set(clang "$<COMPILE_LANG_AND_ID:CXX,Clang>")
 set(msvc "$<COMPILE_LANG_AND_ID:CXX,MSVC>")
 
-set(gcc_like $<OR:${gcc},${clang})
+set(gcc_like $<OR:${gcc},${clang}>)
 
 if (MVPOLYNOMIAL_BUILD_TESTS)
   enable_testing()
diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 9f7e7ca..9d56a57 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -390,8 +390,6 @@ class MVPolynomial final {
   friend MVPolynomial operator-(MVPolynomial&& l, MVPolynomial&& r) { return std::move(l) - r; }
 
   friend MVPolynomial operator*(const MVPolynomial& l, const MVPolynomial& r) {
-    auto comparer = l.key_comp();
-
     auto mul = MVPolynomial(l.get_allocator());
     mul.index2value_.clear();
     // Calculate all product of each l's term and r's term.
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 88c40d3..a4f53de 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -16,7 +16,7 @@ TEST_CASE("constructor", "[mvPolynomial]") {
     auto m   = MP2();
 
     REQUIRE(m.size() == ans.size());
-    for (auto i = 0; i < ans.size(); ++i) {
+    for (size_t i = 0; i < ans.size(); ++i) {
       REQUIRE(m[ans[i].first] == ans[i].second);
     }
   }
@@ -343,7 +343,7 @@ TEST_CASE("sub", "[mvPolynomial]") {
   auto sub = l - r;
 
   REQUIRE(sub.size() == ans.size());
-  for (auto i = 0; i < ans.size(); ++i) {
+  for (size_t i = 0; i < ans.size(); ++i) {
     REQUIRE(sub[ans[i].first] == ans[i].second);
   }
 }

From daf486387348228ce9fd1dce02e785cc7f83d8cf Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 08:01:37 +0900
Subject: [PATCH 16/41] Take measure to be collapsed format.

---
 include/mvPolynomial/mvPolynomial.hpp | 21 ++++++++-------------
 1 file changed, 8 insertions(+), 13 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 9d56a57..3f1aafb 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -77,18 +77,11 @@ class MVPolynomial final {
   MVPolynomial& operator=(MVPolynomial&& other)      = default;
   ~MVPolynomial()                                    = default;
 
-  MVPolynomial()
-      : index2value_({
-            {index_type::Zero(), 0}
-  }) {}
-
-  explicit MVPolynomial(const allocator_type& allocator)
-      : index2value_(
-            {
-                {index_type::Zero(), 0}
-  },
-            allocator
-        ) {}
+  MVPolynomial() : index2value_() { index2value_[index_type::Zero()] = R{0}; }
+
+  explicit MVPolynomial(const allocator_type& allocator) : index2value_(allocator) {
+    index2value_[index_type::Zero()] = R{0};
+  }
 
   template <typename InputIterator>
   MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e) {
@@ -115,7 +108,9 @@ class MVPolynomial final {
     return *this;
   }
 
-  MVPolynomial(mapped_type r, const allocator_type& a = allocator_type{}) : index2value_({{index_type::Zero(), r}}, a) {}
+  MVPolynomial(mapped_type r, const allocator_type& a = allocator_type{}) : index2value_(a) {
+    index2value_[index_type::Zero()] = r;
+  }
 
   allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); }
   key_compare    key_comp() const noexcept { return index2value_.key_comp(); }

From 833b7c1793582a2934bdbdc1d7c4edb739d74b1b Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 08:02:35 +0900
Subject: [PATCH 17/41] Remove empty mem func becuase mvPoly must not be empty.

---
 include/mvPolynomial/mvPolynomial.hpp | 2 --
 1 file changed, 2 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 3f1aafb..fc0d0d0 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -133,8 +133,6 @@ class MVPolynomial final {
   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(); }
-
   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(); }

From 6882a3901e08a59d46866bab1f7a7414c6d23815 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 08:51:47 +0900
Subject: [PATCH 18/41] Remove const_cast.

---
 include/mvPolynomial/mvPolynomial.hpp | 32 +++++++++++------
 test/mvPolynomial_test.cpp            | 52 +++++++++++++++++++++++++++
 2 files changed, 74 insertions(+), 10 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index fc0d0d0..081a427 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -318,12 +318,18 @@ class MVPolynomial final {
   MVPolynomial& operator*=(const MVPolynomial& r) {
     if (r.size() == 1) {
       const auto& [r_index, r_coeff] = *(r.begin());
-      for (auto& index_and_coeff : *this) {
-        auto& index = const_cast<index_type&>(index_and_coeff.first);
-        auto& coeff = index_and_coeff.second;
-        index += r_index;
-        coeff *= r_coeff;
+      auto result                    = MVPolynomial(get_allocator());
+      result.index2value_.clear();
+      for (const auto& [index, coeff] : *this) {
+        const auto new_index = index + r_index;
+        const auto new_coeff = coeff * r_coeff;
+        if (result.contains(new_index)) {
+          result[new_index] += new_coeff;
+        } else {
+          result[new_index] = new_coeff;
+        }
       }
+      swap(result);
     } else {
       *this = *this * r;
     }
@@ -456,12 +462,18 @@ auto Integrate(MVPolynomial<IntType, R, D, Allocator> p, int axis) {
 
   details::CheckAxis(D, axis);
 
-  for (auto& index_and_value : p) {
-    auto& value = index_and_value.second;
-    auto& index = const_cast<typename MP::index_type&>(index_and_value.first);
-    value /= ++index[axis];
+  auto result = MP(std::initializer_list<typename MP::value_type>{}, p.get_allocator());
+  for (const auto& [index, value] : p) {
+    auto new_index = index;
+    ++new_index[axis];
+    const auto new_value = value / new_index[axis];
+    if (result.contains(new_index)) {
+      result[new_index] += new_value;
+    } else {
+      result[new_index] = new_value;
+    }
   }
-  return std::move(p);
+  return result;
 }
 
 }  // namespace mvPolynomial
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index a4f53de..b5410dd 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -257,6 +257,34 @@ TEST_CASE("integral", "[mvPolynomial]") {
       REQUIRE(sm[ans[i].first] == ans[i].second);
     }
   }
+
+  SECTION("constant term preserves searchable ordering") {
+    auto c  = MP2({
+        {{0, 0}, 9},
+    });
+    auto sc = mvPolynomial::Integrate(c, 1);
+
+    REQUIRE(sc.size() == 1);
+    REQUIRE(sc.contains(Eigen::Array2i({0, 1})));
+    REQUIRE(sc.find(Eigen::Array2i({0, 1})) != sc.end());
+    REQUIRE(sc.lower_bound(Eigen::Array2i({0, 1})) != sc.end());
+    REQUIRE((sc.lower_bound(Eigen::Array2i({0, 1}))->first == Eigen::Array2i({0, 1})).all());
+    REQUIRE(sc.at(Eigen::Array2i({0, 1})) == 9);
+  }
+
+  SECTION("higher degree term keeps lookup behavior") {
+    auto high = MP2({
+        {{3, 2}, 8},
+    });
+    auto shi  = mvPolynomial::Integrate(high, 0);
+
+    REQUIRE(shi.size() == 1);
+    REQUIRE(shi.contains(Eigen::Array2i({4, 2})));
+    REQUIRE(shi.find(Eigen::Array2i({4, 2})) != shi.end());
+    REQUIRE(shi.lower_bound(Eigen::Array2i({4, 2})) != shi.end());
+    REQUIRE((shi.lower_bound(Eigen::Array2i({4, 2}))->first == Eigen::Array2i({4, 2})).all());
+    REQUIRE(shi.at(Eigen::Array2i({4, 2})) == 2);
+  }
 }
 
 TEST_CASE("multiply", "[mvPolynomial]") {
@@ -290,6 +318,30 @@ TEST_CASE("multiply", "[mvPolynomial]") {
   for (size_t i = 0; i < ans.size(); ++i) {
     REQUIRE(prod[ans[i].first] == ans[i].second);
   }
+
+  SECTION("multiply assign by monomial matches binary multiply") {
+    auto lhs      = MP2({
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+    auto monomial = MP2({
+        {{2, 1}, 5},
+    });
+    auto expected = MP2({
+        {{3, 1}, 10},
+        {{2, 2}, 15},
+    });
+
+    auto inplace = lhs;
+    inplace *= monomial;
+
+    REQUIRE(inplace == expected);
+    REQUIRE(inplace == lhs * monomial);
+    REQUIRE(inplace.contains(Eigen::Array2i({3, 1})));
+    REQUIRE(inplace.find(Eigen::Array2i({2, 2})) != inplace.end());
+    REQUIRE(inplace.lower_bound(Eigen::Array2i({2, 2})) != inplace.end());
+    REQUIRE((inplace.lower_bound(Eigen::Array2i({2, 2}))->first == Eigen::Array2i({2, 2})).all());
+  }
 }
 
 TEST_CASE("sum", "[mvPolynomial]") {

From 87adb5335e6bf051799b5fb0886011310b9feab3 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 09:51:29 +0900
Subject: [PATCH 19/41] Add abs and relative tolerance.

---
 include/mvPolynomial/mvPolynomial.hpp | 58 +++++++++++++++++++++++----
 1 file changed, 51 insertions(+), 7 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 081a427..7b7e2c7 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -8,8 +8,10 @@
 #include <bit>
 #include <cassert>
 #include <cmath>
+#include <cstddef>
 #include <iterator>
 #include <limits>
+#include <ranges>
 #include <stdexcept>
 #include <type_traits>
 #include <vector>
@@ -38,7 +40,8 @@ class MVPolynomial final {
   static constexpr int dim = D;
 
   // This setting is too strict, so I expect users to set tolerance.
-  static R tolerance;
+  inline static R rel_tolerance = std::numeric_limits<R>::epsilon();
+  inline static R abs_tolerance = std::numeric_limits<R>::min();
 
   using index_type = IndexType<IntType, D>;
   using coord_type = CoordType<R, dim>;
@@ -220,6 +223,31 @@ class MVPolynomial final {
     }
   }
 
+  void DeleteZeroCoeffTerm() {
+    std::vector<size_t> removed_term_indexes;
+    size_t              i                     = 0;
+    bool                should_set_const_zero = false;
+    for (const auto& index_and_coeff : index2value_) {
+      auto        coeff = index_and_coeff.second;
+      const auto& index = index_and_coeff.first;
+      if ((index == index_type::Zero()).all()) {
+        if (coeff < abs_tolerance) {
+          should_set_const_zero = true;
+        }
+      }
+      if (std::abs(coeff) < abs_tolerance) {
+        removed_term_indexes.push_back(i);
+      }
+      ++i;
+    }
+    if (should_set_const_zero) {
+      index2value_[index_type::Zero()] = R(0.0);
+    }
+    for (auto removed_index : removed_term_indexes | std::views::reverse) {
+      index2value_.erase(std::next(index2value_.begin(), removed_index));
+    }
+  }
+
   R operator()(const coord_type& x) const {
     auto sum = R(0.0);
     for (const auto& index_and_coeff : index2value_) {
@@ -271,6 +299,7 @@ class MVPolynomial final {
     } else {
       (*this)[idx] = r;
     }
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -282,6 +311,7 @@ class MVPolynomial final {
         (*this)[idx] = coeff;
       }
     }
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -292,6 +322,7 @@ class MVPolynomial final {
     } else {
       (*this)[idx] = -r;
     }
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -303,6 +334,7 @@ class MVPolynomial final {
         (*this)[idx] = -coeff;
       }
     }
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -312,6 +344,7 @@ class MVPolynomial final {
       auto& coeff = index_and_coeff.second;
       coeff *= r;
     }
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -333,6 +366,7 @@ class MVPolynomial final {
     } else {
       *this = *this * r;
     }
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -351,8 +385,15 @@ class MVPolynomial final {
       if ((l_idx != r_idx).any()) {
         return false;
       }
-      if (std::abs(l_coeff - r_coeff) >= tolerance) {
-        return false;
+      auto abs_diff = std::abs(l_coeff - r_coeff);
+      if (abs_diff < 1) {
+        if (!(abs_diff < abs_tolerance)) {
+          return false;
+        }
+      } else {
+        if (!(abs_diff < rel_tolerance * std::max(std::abs(l_coeff), std::abs(r_coeff)))) {
+          return false;
+        }
       }
       ++l_it;
       ++r_it;
@@ -406,6 +447,8 @@ class MVPolynomial final {
       }
     }
 
+    mul.DeleteZeroCoeffTerm();
+
     return mul;
   }
 
@@ -421,10 +464,6 @@ class MVPolynomial final {
   IndexContainer index2value_;
 };
 
-template <std::signed_integral IntType, std::floating_point R, int Dim, class Allocator>
-R MVPolynomial<IntType, R, Dim, Allocator>::tolerance =
-    std::ldexp(std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::min_exponent);
-
 template <std::signed_integral IntType, std::floating_point R, int Dim, class Allocator>
 auto D(const MVPolynomial<IntType, R, Dim, Allocator>& p, int axis) {
   using MP            = MVPolynomial<IntType, R, Dim, Allocator>;
@@ -453,6 +492,8 @@ auto D(const MVPolynomial<IntType, R, Dim, Allocator>& p, int axis) {
       ++p_it;
     }
   }
+  dp.DeleteZeroCoeffTerm();
+
   return dp;
 }
 
@@ -473,6 +514,9 @@ auto Integrate(MVPolynomial<IntType, R, D, Allocator> p, int axis) {
       result[new_index] = new_value;
     }
   }
+
+  result.DeleteZeroCoeffTerm();
+
   return result;
 }
 

From 75c7f4270b8257faaed3c5bde76ab5e756a46cc2 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 09:55:19 +0900
Subject: [PATCH 20/41] Check index at operator[] and fix format.

---
 include/mvPolynomial/mvPolynomial.hpp | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 7b7e2c7..407c5fe 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -140,8 +140,14 @@ class MVPolynomial final {
   size_type max_size() const noexcept { return index2value_.max_size(); }
   size_type capacity() const noexcept { return index2value_.capacity(); }
 
-  mapped_type& operator[](const key_type& index) { return index2value_[index]; }
-  mapped_type& operator[](key_type&& index) { return index2value_[index]; }
+  mapped_type& operator[](const key_type& index) {
+    CheckIndexIncludingNegative(index);
+    return index2value_[index];
+  }
+  mapped_type& operator[](key_type&& index) {
+    CheckIndexIncludingNegative(index);
+    return index2value_[index];
+  }
 
   mapped_type&       at(const key_type& i) { return index2value_.at(i); }
   const mapped_type& at(const key_type& i) const { return index2value_.at(i); }
@@ -225,9 +231,8 @@ class MVPolynomial final {
 
   void DeleteZeroCoeffTerm() {
     std::vector<size_t> removed_term_indexes;
-    size_t              i                     = 0;
     bool                should_set_const_zero = false;
-    for (const auto& index_and_coeff : index2value_) {
+    for (size_t i = 0; const auto& index_and_coeff : index2value_) {
       auto        coeff = index_and_coeff.second;
       const auto& index = index_and_coeff.first;
       if ((index == index_type::Zero()).all()) {
@@ -453,11 +458,15 @@ class MVPolynomial final {
   }
 
  private:
+  void CheckIndexIncludingNegative(const index_type& index) const {
+    if ((index < 0).any()) {
+      throw std::invalid_argument("Negative index not supported!");
+    }
+  }
+
   void CheckSelfIndexes() const {
     for (const auto& index_and_coeff : index2value_) {
-      if ((index_and_coeff.first < 0).any()) {
-        throw std::invalid_argument("Negative index not supported!");
-      }
+      CheckIndexIncludingNegative(index_and_coeff.first);
     }
   }
 

From 3de42f1b04be9e1c33027112c37ed33adaab8493 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 10:06:58 +0900
Subject: [PATCH 21/41] Add tests and fix how to delete zero coeff.

---
 include/mvPolynomial/mvPolynomial.hpp | 17 +++------
 test/mvPolynomial_test.cpp            | 50 +++++++++++++++++++++++++++
 2 files changed, 54 insertions(+), 13 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 407c5fe..5e8ba7c 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -231,26 +231,18 @@ class MVPolynomial final {
 
   void DeleteZeroCoeffTerm() {
     std::vector<size_t> removed_term_indexes;
-    bool                should_set_const_zero = false;
     for (size_t i = 0; const auto& index_and_coeff : index2value_) {
-      auto        coeff = index_and_coeff.second;
-      const auto& index = index_and_coeff.first;
-      if ((index == index_type::Zero()).all()) {
-        if (coeff < abs_tolerance) {
-          should_set_const_zero = true;
-        }
-      }
-      if (std::abs(coeff) < abs_tolerance) {
+      if (std::abs(index_and_coeff.second) < abs_tolerance) {
         removed_term_indexes.push_back(i);
       }
       ++i;
     }
-    if (should_set_const_zero) {
-      index2value_[index_type::Zero()] = R(0.0);
-    }
     for (auto removed_index : removed_term_indexes | std::views::reverse) {
       index2value_.erase(std::next(index2value_.begin(), removed_index));
     }
+    if (index2value_.empty()) {
+      index2value_[index_type::Zero()] = R{0.0};
+    }
   }
 
   R operator()(const coord_type& x) const {
@@ -344,7 +336,6 @@ class MVPolynomial final {
   }
 
   MVPolynomial& operator*=(mapped_type r) {
-    auto idx = index_type::Zero();
     for (auto& index_and_coeff : *this) {
       auto& coeff = index_and_coeff.second;
       coeff *= r;
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index b5410dd..167dbb4 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -1,6 +1,7 @@
 #include "mvPolynomial/mvPolynomial.hpp"
 
 #include <cstddef>
+#include <stdexcept>
 #include <vector>
 
 #include <catch2/catch_test_macros.hpp>
@@ -63,6 +64,55 @@ TEST_CASE("constructor", "[mvPolynomial]") {
       REQUIRE(m.at(ans[i].first) == ans[i].second);
     }
   }
+
+  SECTION("negative index is rejected") {
+    REQUIRE_THROWS_AS(
+        MP2({
+            {{-1, 0}, 1},
+    }),
+        std::invalid_argument
+    );
+  }
+}
+
+TEST_CASE("invariants", "[mvPolynomial]") {
+  SECTION("subtracting self stays canonical zero polynomial") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+
+    p -= p;
+
+    REQUIRE(p.size() == 1);
+    REQUIRE(p == MP2());
+    REQUIRE(p.contains(Eigen::Array2i::Zero()));
+    REQUIRE(p.at(Eigen::Array2i::Zero()) == 0);
+  }
+
+  SECTION("multiplying by zero stays canonical zero polynomial") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+
+    p *= 0.0;
+
+    REQUIRE(p.size() == 1);
+    REQUIRE(p == MP2());
+    REQUIRE(p.contains(Eigen::Array2i::Zero()));
+    REQUIRE(p.at(Eigen::Array2i::Zero()) == 0);
+  }
+
+  SECTION("operator[] rejects negative index for lvalue and rvalue keys") {
+    auto p       = MP2();
+    auto bad_key = Eigen::Array2i({-1, 0});
+
+    REQUIRE_THROWS_AS(p[bad_key] = 1.0, std::invalid_argument);
+    REQUIRE_THROWS_AS(p[Eigen::Array2i({0, -1})] = 2.0, std::invalid_argument);
+  }
 }
 
 TEST_CASE("pow", "[mvPolynomial]") {

From be04a96ed99b5e35d04909015f0adb179738d62f Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 13:41:56 +0900
Subject: [PATCH 22/41] Optimize operator()(coord_type).

---
 include/mvPolynomial/mvPolynomial.hpp | 26 ++++++++++++++++++++------
 1 file changed, 20 insertions(+), 6 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 5e8ba7c..cd5aef5 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -11,6 +11,7 @@
 #include <cstddef>
 #include <iterator>
 #include <limits>
+#include <numeric>
 #include <ranges>
 #include <stdexcept>
 #include <type_traits>
@@ -246,13 +247,26 @@ class MVPolynomial final {
   }
 
   R operator()(const coord_type& x) const {
-    auto sum = R(0.0);
-    for (const auto& index_and_coeff : index2value_) {
-      auto        coeff = index_and_coeff.second;
-      const auto& index = index_and_coeff.first;
-      sum += coeff * (x.array().pow(index.template cast<double>())).prod();
+    auto           partial_sum = rbegin()->second;
+    std::vector<R> partial_sums;
+    for (auto it = rbegin(); it != std::prev(rend()); ++it) {
+      const auto& index = it->first;
+
+      auto        next_it    = std::next(it);
+      auto        next_coeff = next_it->second;
+      const auto& next_index = next_it->first;
+
+      index_type index_diff = index - next_index;
+      if ((index_diff > 0).all()) {
+        partial_sum = next_coeff + partial_sum * (x.array().pow(index_diff.template cast<R>())).prod();
+      } else {
+        partial_sum *= (x.array().pow(index.template cast<R>())).prod();
+        partial_sums.push_back(partial_sum);
+        partial_sum = next_coeff;
+      }
     }
-    return sum;
+    // Because partial_sum is constant, so
+    return std::reduce(partial_sums.cbegin(), partial_sums.cend()) + partial_sum;
   }
 
   MVPolynomial operator()(const MVPolynomial& mvp, int axis) const {

From 1167380a3911994cb2b312d05c459d5dbd4af1bd Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 18:28:22 +0900
Subject: [PATCH 23/41] Fix a bug.

---
 include/mvPolynomial/mvPolynomial.hpp |  7 +++++--
 test/mvPolynomial_test.cpp            | 13 +++++++++++++
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index cd5aef5..40f02fc 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -247,8 +247,8 @@ class MVPolynomial final {
   }
 
   R operator()(const coord_type& x) const {
-    auto           partial_sum = rbegin()->second;
     std::vector<R> partial_sums;
+    auto           partial_sum = rbegin()->second;
     for (auto it = rbegin(); it != std::prev(rend()); ++it) {
       const auto& index = it->first;
 
@@ -265,7 +265,10 @@ class MVPolynomial final {
         partial_sum = next_coeff;
       }
     }
-    // Because partial_sum is constant, so
+    {
+      const auto& first_index = begin()->first;
+      partial_sum *= (x.array().pow(first_index.template cast<R>())).prod();
+    }
     return std::reduce(partial_sums.cbegin(), partial_sums.cend()) + partial_sum;
   }
 
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 167dbb4..2a15d95 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -152,6 +152,19 @@ TEST_CASE("operator()", "[mvPolynomial]") {
     REQUIRE(m(Eigen::Vector2d({2, 3})) == 112);
   }
 
+  SECTION("point without constant") {
+    auto m = MP2({
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+        {{1, 1}, 4},
+        {{2, 0}, 5},
+        {{0, 2}, 6},
+    });
+
+    REQUIRE(m(Eigen::Vector2d::Zero()) == 0);
+    REQUIRE(m(Eigen::Vector2d({2, 3})) == 111);
+  }
+
   SECTION("mvpolynomial_1d") {
     auto m = MP3({
         {{0, 0, 0}, 1},

From 27b2070b67ea4ce7f254276b1a617667734849cd Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 18:55:47 +0900
Subject: [PATCH 24/41] Optimize operator()(MVPolynomial).

---
 include/mvPolynomial/mvPolynomial.hpp | 40 +++++++++++++++++++++------
 test/mvPolynomial_test.cpp            | 20 ++++++++++++++
 2 files changed, 52 insertions(+), 8 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 40f02fc..ab3345b 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -275,16 +275,40 @@ class MVPolynomial final {
   MVPolynomial operator()(const MVPolynomial& mvp, int axis) const {
     details::CheckAxis(dim, axis);
 
-    auto composed_mvp = MVPolynomial{get_allocator()};
-    for (const auto& index_and_coeff : index2value_) {
-      auto        coeff = index_and_coeff.second;
-      const auto& index = index_and_coeff.first;
+    auto partial_sums = std::vector<MVPolynomial>{get_allocator()};
+    auto partial_sum  = MVPolynomial{rbegin()->second, get_allocator()};
+    for (auto it = rbegin(); it != std::prev(rend()); ++it) {
+      const auto& index = it->first;
+
+      auto        next_it    = std::next(it);
+      auto        next_coeff = next_it->second;
+      const auto& next_index = next_it->first;
+
+      index_type index_diff              = index - next_index;
+      index_type index_diff_without_axis = index_diff;
+      index_diff_without_axis[axis]      = 0;
 
-      auto index_with_axis_0  = index;
-      index_with_axis_0[axis] = 0;
-      composed_mvp += MVPolynomial{{{index_with_axis_0, coeff}}, get_allocator()} * mvp.pow(index[axis]);
+      if ((index_diff > 0).all()) {
+        auto tmp_mvp = MVPolynomial{{{index_diff_without_axis, 1}}, get_allocator()};
+        partial_sum  = next_coeff + partial_sum * (mvp.pow(index[axis])) * tmp_mvp;
+      } else {
+        index_type index_without_axis = index;
+        index_without_axis[axis]      = 0;
+        auto tmp_mvp                  = MVPolynomial{{{index_without_axis, 1}}, get_allocator()};
+
+        partial_sum *= mvp.pow(index[axis]) * tmp_mvp;
+        partial_sums.push_back(std::move(partial_sum));
+        partial_sum = MVPolynomial{next_coeff, get_allocator()};
+      }
     }
-    return composed_mvp;
+    {
+      const auto& first_index              = begin()->first;
+      index_type  first_index_without_axis = first_index;
+      first_index_without_axis[axis]       = 0;
+      auto tmp_mvp                         = MVPolynomial{{{first_index_without_axis, 1}}, get_allocator()};
+      partial_sum *= mvp.pow(first_index[axis]) * tmp_mvp;
+    }
+    return std::reduce(partial_sums.cbegin(), partial_sums.cend()) + partial_sum;
   }
 
   MVPolynomial operator+() const { return *this; }
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 2a15d95..99f95f0 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -198,6 +198,26 @@ TEST_CASE("operator()", "[mvPolynomial]") {
     );
   }
 
+  SECTION("mvpolynomial_1d_without_constant") {
+    auto m = MP3({
+        {{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},  4},
+                       {{1, 0, 0}, 10},
+                       {{0, 1, 0}, 15},
+    })
+    );
+  }
+
   SECTION("mvpolynomial_2d") {
     auto m = MP3({
         {{0, 0, 0}, 1},

From 423e9ec4462b4b4ad8667ee08ce27553cbdc9239 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 19:02:12 +0900
Subject: [PATCH 25/41] Fix to pass allocator.

---
 include/mvPolynomial/mvPolynomial.hpp | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index ab3345b..8134b67 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -81,14 +81,14 @@ class MVPolynomial final {
   MVPolynomial& operator=(MVPolynomial&& other)      = default;
   ~MVPolynomial()                                    = default;
 
-  MVPolynomial() : index2value_() { index2value_[index_type::Zero()] = R{0}; }
+  MVPolynomial() : index2value_(allocator_type{}) { index2value_[index_type::Zero()] = R{0}; }
 
   explicit MVPolynomial(const allocator_type& allocator) : index2value_(allocator) {
     index2value_[index_type::Zero()] = R{0};
   }
 
   template <typename InputIterator>
-  MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e) {
+  MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e, allocator_type{}) {
     CheckSelfIndexes();
   }
 
@@ -231,7 +231,7 @@ class MVPolynomial final {
   }
 
   void DeleteZeroCoeffTerm() {
-    std::vector<size_t> removed_term_indexes;
+    auto removed_term_indexes = std::vector<size_t>{get_allocator()};
     for (size_t i = 0; const auto& index_and_coeff : index2value_) {
       if (std::abs(index_and_coeff.second) < abs_tolerance) {
         removed_term_indexes.push_back(i);
@@ -247,8 +247,8 @@ class MVPolynomial final {
   }
 
   R operator()(const coord_type& x) const {
-    std::vector<R> partial_sums;
-    auto           partial_sum = rbegin()->second;
+    auto partial_sums = std::vector<R>{get_allocator()};
+    auto partial_sum  = rbegin()->second;
     for (auto it = rbegin(); it != std::prev(rend()); ++it) {
       const auto& index = it->first;
 

From c7da6dec88d7f110893afb2ddcaa3e0d0447eb18 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 21:25:00 +0900
Subject: [PATCH 26/41] Fix bugs destructing normalization of MVPoly.

---
 include/mvPolynomial/mvPolynomial.hpp | 59 ++++++++++++++-------------
 test/mvPolynomial_test.cpp            | 34 +++++++--------
 2 files changed, 47 insertions(+), 46 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 8134b67..1fd6ca7 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -90,15 +90,18 @@ class MVPolynomial final {
   template <typename InputIterator>
   MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e, allocator_type{}) {
     CheckSelfIndexes();
+    DeleteZeroCoeffTerm();
   }
 
   template <typename InputIterator>
   MVPolynomial(InputIterator s, InputIterator e, const allocator_type& allocator) : index2value_(s, e, allocator) {
     CheckSelfIndexes();
+    DeleteZeroCoeffTerm();
   }
 
   MVPolynomial(std::initializer_list<value_type> l, const allocator_type& a = allocator_type{}) : index2value_(l, a) {
     CheckSelfIndexes();
+    DeleteZeroCoeffTerm();
   }
 
   MVPolynomial(const MVPolynomial& m, const allocator_type& a) : index2value_(m.index2value_, a) {}
@@ -109,11 +112,13 @@ class MVPolynomial final {
     index2value_.clear();
     index2value_.insert(l.begin(), l.end());
     CheckSelfIndexes();
+    DeleteZeroCoeffTerm();
     return *this;
   }
 
   MVPolynomial(mapped_type r, const allocator_type& a = allocator_type{}) : index2value_(a) {
     index2value_[index_type::Zero()] = r;
+    DeleteZeroCoeffTerm();
   }
 
   allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); }
@@ -141,18 +146,14 @@ class MVPolynomial final {
   size_type max_size() const noexcept { return index2value_.max_size(); }
   size_type capacity() const noexcept { return index2value_.capacity(); }
 
-  mapped_type& operator[](const key_type& index) {
-    CheckIndexIncludingNegative(index);
-    return index2value_[index];
-  }
-  mapped_type& operator[](key_type&& index) {
+  const mapped_type& get(const key_type& index) const { return index2value_.at(index); }
+
+  void set(const key_type& index, R coeff) {
     CheckIndexIncludingNegative(index);
-    return index2value_[index];
+    index2value_[index] = coeff;
+    DeleteZeroCoeffTerm();
   }
 
-  mapped_type&       at(const key_type& i) { return index2value_.at(i); }
-  const mapped_type& at(const key_type& i) const { return index2value_.at(i); }
-
   void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); }
 
   iterator       find(const key_type& i) { return index2value_.find(i); }
@@ -208,7 +209,7 @@ class MVPolynomial final {
       default:
         auto max_pow2_under_exp = std::bit_floor(static_cast<unsigned int>(exp));
         auto max_bit_width      = std::bit_width(max_pow2_under_exp);
-        auto cache              = std::vector<MVPolynomial>(max_bit_width - 1, get_allocator());
+        auto cache              = std::vector<MVPolynomial>(max_bit_width - 1);
         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);
@@ -231,7 +232,7 @@ class MVPolynomial final {
   }
 
   void DeleteZeroCoeffTerm() {
-    auto removed_term_indexes = std::vector<size_t>{get_allocator()};
+    auto removed_term_indexes = std::vector<size_t>{};
     for (size_t i = 0; const auto& index_and_coeff : index2value_) {
       if (std::abs(index_and_coeff.second) < abs_tolerance) {
         removed_term_indexes.push_back(i);
@@ -247,7 +248,7 @@ class MVPolynomial final {
   }
 
   R operator()(const coord_type& x) const {
-    auto partial_sums = std::vector<R>{get_allocator()};
+    auto partial_sums = std::vector<R>{};
     auto partial_sum  = rbegin()->second;
     for (auto it = rbegin(); it != std::prev(rend()); ++it) {
       const auto& index = it->first;
@@ -275,7 +276,7 @@ class MVPolynomial final {
   MVPolynomial operator()(const MVPolynomial& mvp, int axis) const {
     details::CheckAxis(dim, axis);
 
-    auto partial_sums = std::vector<MVPolynomial>{get_allocator()};
+    auto partial_sums = std::vector<MVPolynomial>{};
     auto partial_sum  = MVPolynomial{rbegin()->second, get_allocator()};
     for (auto it = rbegin(); it != std::prev(rend()); ++it) {
       const auto& index = it->first;
@@ -333,9 +334,9 @@ class MVPolynomial final {
   MVPolynomial& operator+=(mapped_type r) {
     auto idx = index_type::Zero();
     if (contains(idx)) {
-      (*this)[idx] += r;
+      index2value_[idx] += r;
     } else {
-      (*this)[idx] = r;
+      index2value_[idx] = r;
     }
     DeleteZeroCoeffTerm();
     return *this;
@@ -344,9 +345,9 @@ class MVPolynomial final {
   MVPolynomial& operator+=(const MVPolynomial& r) {
     for (const auto& [idx, coeff] : r) {
       if (contains(idx)) {
-        (*this)[idx] += coeff;
+        index2value_[idx] += coeff;
       } else {
-        (*this)[idx] = coeff;
+        index2value_[idx] = coeff;
       }
     }
     DeleteZeroCoeffTerm();
@@ -356,9 +357,9 @@ class MVPolynomial final {
   MVPolynomial& operator-=(mapped_type r) {
     auto idx = index_type::Zero();
     if (contains(idx)) {
-      (*this)[idx] -= r;
+      index2value_[idx] -= r;
     } else {
-      (*this)[idx] = -r;
+      index2value_[idx] = -r;
     }
     DeleteZeroCoeffTerm();
     return *this;
@@ -367,9 +368,9 @@ class MVPolynomial final {
   MVPolynomial& operator-=(const MVPolynomial& r) {
     for (const auto& [idx, coeff] : r) {
       if (contains(idx)) {
-        (*this)[idx] -= coeff;
+        index2value_[idx] -= coeff;
       } else {
-        (*this)[idx] = -coeff;
+        index2value_[idx] = -coeff;
       }
     }
     DeleteZeroCoeffTerm();
@@ -394,9 +395,9 @@ class MVPolynomial final {
         const auto new_index = index + r_index;
         const auto new_coeff = coeff * r_coeff;
         if (result.contains(new_index)) {
-          result[new_index] += new_coeff;
+          result.index2value_[new_index] += new_coeff;
         } else {
-          result[new_index] = new_coeff;
+          result.index2value_[new_index] = new_coeff;
         }
       }
       swap(result);
@@ -477,9 +478,9 @@ class MVPolynomial final {
         const auto idx           = l_idx + r_idx;
         const auto v             = l_v * r_v;
         if (mul.contains(idx)) {
-          mul[idx] += v;
+          mul.index2value_[idx] += v;
         } else {
-          mul[idx] = v;
+          mul.index2value_[idx] = v;
         }
       }
     }
@@ -529,7 +530,7 @@ auto D(const MVPolynomial<IntType, R, Dim, Allocator>& p, int axis) {
       p_it = d_end_it;
     } else {
       value *= index[axis]--;
-      dp[index] = value;
+      dp.set(index, value);
       ++p_it;
     }
   }
@@ -544,15 +545,15 @@ auto Integrate(MVPolynomial<IntType, R, D, Allocator> p, int axis) {
 
   details::CheckAxis(D, axis);
 
-  auto result = MP(std::initializer_list<typename MP::value_type>{}, p.get_allocator());
+  auto result = MP(p.get_allocator());
   for (const auto& [index, value] : p) {
     auto new_index = index;
     ++new_index[axis];
     const auto new_value = value / new_index[axis];
     if (result.contains(new_index)) {
-      result[new_index] += new_value;
+      result.set(new_index, result.get(new_index) + new_value);
     } else {
-      result[new_index] = new_value;
+      result.set(new_index, new_value);
     }
   }
 
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 99f95f0..a2ae23b 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -18,7 +18,7 @@ TEST_CASE("constructor", "[mvPolynomial]") {
 
     REQUIRE(m.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(m[ans[i].first] == ans[i].second);
+      REQUIRE(m.get(ans[i].first) == ans[i].second);
     }
   }
 
@@ -36,7 +36,7 @@ TEST_CASE("constructor", "[mvPolynomial]") {
 
     REQUIRE(m.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(m.at(ans[i].first) == ans[i].second);
+      REQUIRE(m.get(ans[i].first) == ans[i].second);
     }
   }
 
@@ -61,7 +61,7 @@ TEST_CASE("constructor", "[mvPolynomial]") {
 
     REQUIRE(m.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(m.at(ans[i].first) == ans[i].second);
+      REQUIRE(m.get(ans[i].first) == ans[i].second);
     }
   }
 
@@ -88,7 +88,7 @@ TEST_CASE("invariants", "[mvPolynomial]") {
     REQUIRE(p.size() == 1);
     REQUIRE(p == MP2());
     REQUIRE(p.contains(Eigen::Array2i::Zero()));
-    REQUIRE(p.at(Eigen::Array2i::Zero()) == 0);
+    REQUIRE(p.get(Eigen::Array2i::Zero()) == 0);
   }
 
   SECTION("multiplying by zero stays canonical zero polynomial") {
@@ -103,15 +103,15 @@ TEST_CASE("invariants", "[mvPolynomial]") {
     REQUIRE(p.size() == 1);
     REQUIRE(p == MP2());
     REQUIRE(p.contains(Eigen::Array2i::Zero()));
-    REQUIRE(p.at(Eigen::Array2i::Zero()) == 0);
+    REQUIRE(p.get(Eigen::Array2i::Zero()) == 0);
   }
 
-  SECTION("operator[] rejects negative index for lvalue and rvalue keys") {
+  SECTION("set rejects negative index for lvalue and rvalue keys") {
     auto p       = MP2();
     auto bad_key = Eigen::Array2i({-1, 0});
 
-    REQUIRE_THROWS_AS(p[bad_key] = 1.0, std::invalid_argument);
-    REQUIRE_THROWS_AS(p[Eigen::Array2i({0, -1})] = 2.0, std::invalid_argument);
+    REQUIRE_THROWS_AS(p.set(bad_key, 1.0), std::invalid_argument);
+    REQUIRE_THROWS_AS(p.set(Eigen::Array2i({0, -1}), 2.0), std::invalid_argument);
   }
 }
 
@@ -278,7 +278,7 @@ TEST_CASE("D", "[mvPolynomial]") {
 
     REQUIRE(dm0.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(dm0[ans[i].first] == ans[i].second);
+      REQUIRE(dm0.get(ans[i].first) == ans[i].second);
     }
   }
 
@@ -292,7 +292,7 @@ TEST_CASE("D", "[mvPolynomial]") {
 
     REQUIRE(dm1.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(dm1[ans[i].first] == ans[i].second);
+      REQUIRE(dm1.get(ans[i].first) == ans[i].second);
     }
   }
 }
@@ -320,7 +320,7 @@ TEST_CASE("integral", "[mvPolynomial]") {
 
     REQUIRE(sm.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(sm[ans[i].first] == ans[i].second);
+      REQUIRE(sm.get(ans[i].first) == ans[i].second);
     }
   }
 
@@ -337,7 +337,7 @@ TEST_CASE("integral", "[mvPolynomial]") {
 
     REQUIRE(sm.size() == ans.size());
     for (size_t i = 0; i < ans.size(); ++i) {
-      REQUIRE(sm[ans[i].first] == ans[i].second);
+      REQUIRE(sm.get(ans[i].first) == ans[i].second);
     }
   }
 
@@ -352,7 +352,7 @@ TEST_CASE("integral", "[mvPolynomial]") {
     REQUIRE(sc.find(Eigen::Array2i({0, 1})) != sc.end());
     REQUIRE(sc.lower_bound(Eigen::Array2i({0, 1})) != sc.end());
     REQUIRE((sc.lower_bound(Eigen::Array2i({0, 1}))->first == Eigen::Array2i({0, 1})).all());
-    REQUIRE(sc.at(Eigen::Array2i({0, 1})) == 9);
+    REQUIRE(sc.get(Eigen::Array2i({0, 1})) == 9);
   }
 
   SECTION("higher degree term keeps lookup behavior") {
@@ -366,7 +366,7 @@ TEST_CASE("integral", "[mvPolynomial]") {
     REQUIRE(shi.find(Eigen::Array2i({4, 2})) != shi.end());
     REQUIRE(shi.lower_bound(Eigen::Array2i({4, 2})) != shi.end());
     REQUIRE((shi.lower_bound(Eigen::Array2i({4, 2}))->first == Eigen::Array2i({4, 2})).all());
-    REQUIRE(shi.at(Eigen::Array2i({4, 2})) == 2);
+    REQUIRE(shi.get(Eigen::Array2i({4, 2})) == 2);
   }
 }
 
@@ -399,7 +399,7 @@ TEST_CASE("multiply", "[mvPolynomial]") {
 
   REQUIRE(prod.size() == ans.size());
   for (size_t i = 0; i < ans.size(); ++i) {
-    REQUIRE(prod[ans[i].first] == ans[i].second);
+    REQUIRE(prod.get(ans[i].first) == ans[i].second);
   }
 
   SECTION("multiply assign by monomial matches binary multiply") {
@@ -452,7 +452,7 @@ TEST_CASE("sum", "[mvPolynomial]") {
 
   REQUIRE(sum.size() == ans.size());
   for (size_t i = 0; i < ans.size(); ++i) {
-    REQUIRE(sum[ans[i].first] == ans[i].second);
+    REQUIRE(sum.get(ans[i].first) == ans[i].second);
   }
 }
 
@@ -479,6 +479,6 @@ TEST_CASE("sub", "[mvPolynomial]") {
 
   REQUIRE(sub.size() == ans.size());
   for (size_t i = 0; i < ans.size(); ++i) {
-    REQUIRE(sub[ans[i].first] == ans[i].second);
+    REQUIRE(sub.get(ans[i].first) == ans[i].second);
   }
 }

From f65c4ac41ebf2afe6e0e93636b067234e0013000 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 21:29:21 +0900
Subject: [PATCH 27/41] Add tests.

---
 test/mvPolynomial_test.cpp | 49 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)

diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index a2ae23b..afc5670 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -40,6 +40,17 @@ TEST_CASE("constructor", "[mvPolynomial]") {
     }
   }
 
+  SECTION("empty range is canonical zero polynomial") {
+    auto empty = std::vector<std::pair<Eigen::Array2i, double>>{};
+    auto m     = MP2(empty.begin(), empty.end());
+
+    REQUIRE(m.size() == 1);
+    REQUIRE(m == MP2());
+    REQUIRE(m.contains(Eigen::Array2i::Zero()));
+    REQUIRE(m.get(Eigen::Array2i::Zero()) == 0);
+    REQUIRE(m(Eigen::Vector2d::Zero()) == 0);
+  }
+
   SECTION("initializer_list") {
     auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
         {{0, 0}, 1},
@@ -65,6 +76,30 @@ TEST_CASE("constructor", "[mvPolynomial]") {
     }
   }
 
+  SECTION("empty initializer_list is canonical zero polynomial") {
+    auto m = MP2({});
+
+    REQUIRE(m.size() == 1);
+    REQUIRE(m == MP2());
+    REQUIRE(m.contains(Eigen::Array2i::Zero()));
+    REQUIRE(m.get(Eigen::Array2i::Zero()) == 0);
+    REQUIRE(m(Eigen::Vector2d({2, 3})) == 0);
+  }
+
+  SECTION("zero coefficient terms are normalized away at construction") {
+    auto m = MP2({
+        {{1, 0}, 0},
+        {{0, 1}, 0},
+        {{0, 0}, 0},
+    });
+
+    REQUIRE(m.size() == 1);
+    REQUIRE(m == MP2());
+    REQUIRE(m.contains(Eigen::Array2i::Zero()));
+    REQUIRE(m.get(Eigen::Array2i::Zero()) == 0);
+    REQUIRE(m(Eigen::Vector2d({2, 3})) == 0);
+  }
+
   SECTION("negative index is rejected") {
     REQUIRE_THROWS_AS(
         MP2({
@@ -113,6 +148,20 @@ TEST_CASE("invariants", "[mvPolynomial]") {
     REQUIRE_THROWS_AS(p.set(bad_key, 1.0), std::invalid_argument);
     REQUIRE_THROWS_AS(p.set(Eigen::Array2i({0, -1}), 2.0), std::invalid_argument);
   }
+
+  SECTION("setting a coefficient to zero removes the term") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+    });
+
+    p.set(Eigen::Array2i({1, 0}), 0.0);
+
+    REQUIRE(p.size() == 1);
+    REQUIRE(p == MP2(1));
+    REQUIRE_FALSE(p.contains(Eigen::Array2i({1, 0})));
+    REQUIRE(p.get(Eigen::Array2i::Zero()) == 1);
+  }
 }
 
 TEST_CASE("pow", "[mvPolynomial]") {

From 06d6f4e64b3be490e2cc055ee491fec2dda3f2f6 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 21:34:16 +0900
Subject: [PATCH 28/41] Fix bugs breaking 0 coeff normalization.

---
 include/mvPolynomial/mvPolynomial.hpp | 27 +++------------------------
 1 file changed, 3 insertions(+), 24 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 1fd6ca7..b90418f 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -124,16 +124,10 @@ class MVPolynomial final {
   allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); }
   key_compare    key_comp() const noexcept { return index2value_.key_comp(); }
 
-  iterator       begin() noexcept { return index2value_.begin(); }
   const_iterator begin() const noexcept { return index2value_.begin(); }
-
-  iterator       end() noexcept { return index2value_.end(); }
   const_iterator end() const noexcept { return index2value_.end(); }
 
-  reverse_iterator       rbegin() noexcept { return index2value_.rbegin(); }
   const_reverse_iterator rbegin() const noexcept { return index2value_.rbegin(); }
-
-  reverse_iterator       rend() noexcept { return index2value_.rend(); }
   const_reverse_iterator rend() const noexcept { return index2value_.rend(); }
 
   const_iterator cbegin() const noexcept { return index2value_.cbegin(); }
@@ -156,13 +150,8 @@ class MVPolynomial final {
 
   void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); }
 
-  iterator       find(const key_type& i) { return index2value_.find(i); }
   const_iterator find(const key_type& i) const { return index2value_.find(i); }
   template <typename K>
-  iterator find(const K& i) {
-    return index2value_.find(i);
-  }
-  template <typename K>
   const_iterator find(const K& i) const {
     return index2value_.find(i);
   }
@@ -173,24 +162,14 @@ class MVPolynomial final {
     return index2value_.contains(i);
   }
 
-  iterator       lower_bound(const key_type& i) { return index2value_.lower_bound(i); }
   const_iterator lower_bound(const key_type& i) const { return index2value_.lower_bound(i); }
   template <typename K>
-  iterator lower_bound(const K& i) {
-    return index2value_.lower_bound(i);
-  }
-  template <typename K>
   const_iterator lower_bound(const K& i) const {
     return index2value_.lower_bound(i);
   }
 
-  iterator       upper_bound(const key_type& i) { return index2value_.upper_bound(i); }
   const_iterator upper_bound(const key_type& i) const { return index2value_.upper_bound(i); }
   template <typename K>
-  iterator upper_bound(const K& i) {
-    return index2value_.upper_bound(i);
-  }
-  template <typename K>
   const_iterator upper_bound(const K& i) const {
     return index2value_.upper_bound(i);
   }
@@ -316,7 +295,7 @@ class MVPolynomial final {
 
   MVPolynomial operator-() const& {
     auto m = MVPolynomial(*this, get_allocator());
-    for (auto& i_and_v : m) {
+    for (auto& i_and_v : m.index2value_) {
       auto& [_, v] = i_and_v;
       v            = -v;
     }
@@ -324,7 +303,7 @@ class MVPolynomial final {
   }
 
   MVPolynomial operator-() && {
-    for (auto& i_and_v : *this) {
+    for (auto& i_and_v : index2value_) {
       auto& [_, v] = i_and_v;
       v            = -v;
     }
@@ -378,7 +357,7 @@ class MVPolynomial final {
   }
 
   MVPolynomial& operator*=(mapped_type r) {
-    for (auto& index_and_coeff : *this) {
+    for (auto& index_and_coeff : index2value_) {
       auto& coeff = index_and_coeff.second;
       coeff *= r;
     }

From 72e976da3b3ce99558618ef0734b380eee0b437b Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 21:38:49 +0900
Subject: [PATCH 29/41] Unify index access to index2value_.

---
 include/mvPolynomial/mvPolynomial.hpp | 63 +++++++++++----------------
 1 file changed, 25 insertions(+), 38 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index b90418f..77d79a6 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -81,10 +81,10 @@ class MVPolynomial final {
   MVPolynomial& operator=(MVPolynomial&& other)      = default;
   ~MVPolynomial()                                    = default;
 
-  MVPolynomial() : index2value_(allocator_type{}) { index2value_[index_type::Zero()] = R{0}; }
+  MVPolynomial() : index2value_(allocator_type{}) { AssignCoeffWithoutNormalization(index_type::Zero(), R{0}); }
 
   explicit MVPolynomial(const allocator_type& allocator) : index2value_(allocator) {
-    index2value_[index_type::Zero()] = R{0};
+    AssignCoeffWithoutNormalization(index_type::Zero(), R{0});
   }
 
   template <typename InputIterator>
@@ -117,7 +117,7 @@ class MVPolynomial final {
   }
 
   MVPolynomial(mapped_type r, const allocator_type& a = allocator_type{}) : index2value_(a) {
-    index2value_[index_type::Zero()] = r;
+    AssignCoeffWithoutNormalization(index_type::Zero(), r);
     DeleteZeroCoeffTerm();
   }
 
@@ -143,8 +143,7 @@ class MVPolynomial final {
   const mapped_type& get(const key_type& index) const { return index2value_.at(index); }
 
   void set(const key_type& index, R coeff) {
-    CheckIndexIncludingNegative(index);
-    index2value_[index] = coeff;
+    AssignCoeffWithoutNormalization(index, coeff);
     DeleteZeroCoeffTerm();
   }
 
@@ -222,7 +221,7 @@ class MVPolynomial final {
       index2value_.erase(std::next(index2value_.begin(), removed_index));
     }
     if (index2value_.empty()) {
-      index2value_[index_type::Zero()] = R{0.0};
+      AssignCoeffWithoutNormalization(index_type::Zero(), R{0});
     }
   }
 
@@ -311,46 +310,28 @@ class MVPolynomial final {
   }
 
   MVPolynomial& operator+=(mapped_type r) {
-    auto idx = index_type::Zero();
-    if (contains(idx)) {
-      index2value_[idx] += r;
-    } else {
-      index2value_[idx] = r;
-    }
+    AddCoeffWithoutNormalization(index_type::Zero(), r);
     DeleteZeroCoeffTerm();
     return *this;
   }
 
   MVPolynomial& operator+=(const MVPolynomial& r) {
     for (const auto& [idx, coeff] : r) {
-      if (contains(idx)) {
-        index2value_[idx] += coeff;
-      } else {
-        index2value_[idx] = coeff;
-      }
+      AddCoeffWithoutNormalization(idx, coeff);
     }
     DeleteZeroCoeffTerm();
     return *this;
   }
 
   MVPolynomial& operator-=(mapped_type r) {
-    auto idx = index_type::Zero();
-    if (contains(idx)) {
-      index2value_[idx] -= r;
-    } else {
-      index2value_[idx] = -r;
-    }
+    AddCoeffWithoutNormalization(index_type::Zero(), -r);
     DeleteZeroCoeffTerm();
     return *this;
   }
 
   MVPolynomial& operator-=(const MVPolynomial& r) {
     for (const auto& [idx, coeff] : r) {
-      if (contains(idx)) {
-        index2value_[idx] -= coeff;
-      } else {
-        index2value_[idx] = -coeff;
-      }
+      AddCoeffWithoutNormalization(idx, -coeff);
     }
     DeleteZeroCoeffTerm();
     return *this;
@@ -373,11 +354,7 @@ class MVPolynomial final {
       for (const auto& [index, coeff] : *this) {
         const auto new_index = index + r_index;
         const auto new_coeff = coeff * r_coeff;
-        if (result.contains(new_index)) {
-          result.index2value_[new_index] += new_coeff;
-        } else {
-          result.index2value_[new_index] = new_coeff;
-        }
+        result.AddCoeffWithoutNormalization(new_index, new_coeff);
       }
       swap(result);
     } else {
@@ -456,11 +433,7 @@ class MVPolynomial final {
         const auto& [r_idx, r_v] = r_p;
         const auto idx           = l_idx + r_idx;
         const auto v             = l_v * r_v;
-        if (mul.contains(idx)) {
-          mul.index2value_[idx] += v;
-        } else {
-          mul.index2value_[idx] = v;
-        }
+        mul.AddCoeffWithoutNormalization(idx, v);
       }
     }
 
@@ -470,6 +443,20 @@ class MVPolynomial final {
   }
 
  private:
+  void AssignCoeffWithoutNormalization(const key_type& index, mapped_type coeff) {
+    CheckIndexIncludingNegative(index);
+    index2value_[index] = coeff;
+  }
+
+  void AddCoeffWithoutNormalization(const key_type& index, mapped_type delta) {
+    CheckIndexIncludingNegative(index);
+    if (auto it = index2value_.find(index); it != index2value_.end()) {
+      it->second += delta;
+    } else {
+      index2value_[index] = delta;
+    }
+  }
+
   void CheckIndexIncludingNegative(const index_type& index) const {
     if ((index < 0).any()) {
       throw std::invalid_argument("Negative index not supported!");

From c845508e6424096905135eb21cb351c23cdd4352 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 21:46:07 +0900
Subject: [PATCH 30/41] Unify allocator type.

---
 include/mvPolynomial/mvPolynomial.hpp | 22 +++++--
 test/mvPolynomial_test.cpp            | 92 +++++++++++++++++++++++++++
 2 files changed, 108 insertions(+), 6 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 77d79a6..d1e19bc 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -11,6 +11,7 @@
 #include <cstddef>
 #include <iterator>
 #include <limits>
+#include <memory>
 #include <numeric>
 #include <ranges>
 #include <stdexcept>
@@ -75,6 +76,14 @@ class MVPolynomial final {
   using reverse_iterator       = IndexContainer::reverse_iterator;
   using const_reverse_iterator = IndexContainer::const_reverse_iterator;
 
+ private:
+  template <typename T>
+  using rebound_allocator_type = typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
+
+  template <typename T>
+  using rebound_vector_type = std::vector<T, rebound_allocator_type<T>>;
+
+ public:
   MVPolynomial(const MVPolynomial& other)            = default;
   MVPolynomial& operator=(const MVPolynomial& other) = default;
   MVPolynomial(MVPolynomial&& other)                 = default;
@@ -187,10 +196,11 @@ class MVPolynomial final {
       default:
         auto max_pow2_under_exp = std::bit_floor(static_cast<unsigned int>(exp));
         auto max_bit_width      = std::bit_width(max_pow2_under_exp);
-        auto cache              = std::vector<MVPolynomial>(max_bit_width - 1);
-        cache.at(0)             = (*this) * (*this);
+        auto cache = rebound_vector_type<MVPolynomial>(rebound_allocator_type<MVPolynomial>(get_allocator()));
+        cache.reserve(max_bit_width - 1);
+        cache.push_back((*this) * (*this));
         for (int i = 2; i < max_bit_width; ++i) {
-          cache.at(i - 1) = cache.at(i - 2) * cache.at(i - 2);
+          cache.push_back(cache.at(i - 2) * cache.at(i - 2));
         }
 
         assert(!cache.empty());
@@ -210,7 +220,7 @@ class MVPolynomial final {
   }
 
   void DeleteZeroCoeffTerm() {
-    auto removed_term_indexes = std::vector<size_t>{};
+    auto removed_term_indexes = rebound_vector_type<size_type>(rebound_allocator_type<size_type>(get_allocator()));
     for (size_t i = 0; const auto& index_and_coeff : index2value_) {
       if (std::abs(index_and_coeff.second) < abs_tolerance) {
         removed_term_indexes.push_back(i);
@@ -226,7 +236,7 @@ class MVPolynomial final {
   }
 
   R operator()(const coord_type& x) const {
-    auto partial_sums = std::vector<R>{};
+    auto partial_sums = rebound_vector_type<R>(rebound_allocator_type<R>(get_allocator()));
     auto partial_sum  = rbegin()->second;
     for (auto it = rbegin(); it != std::prev(rend()); ++it) {
       const auto& index = it->first;
@@ -254,7 +264,7 @@ class MVPolynomial final {
   MVPolynomial operator()(const MVPolynomial& mvp, int axis) const {
     details::CheckAxis(dim, axis);
 
-    auto partial_sums = std::vector<MVPolynomial>{};
+    auto partial_sums = rebound_vector_type<MVPolynomial>(rebound_allocator_type<MVPolynomial>(get_allocator()));
     auto partial_sum  = MVPolynomial{rbegin()->second, get_allocator()};
     for (auto it = rbegin(); it != std::prev(rend()); ++it) {
       const auto& index = it->first;
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index afc5670..db8395b 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -1,7 +1,9 @@
 #include "mvPolynomial/mvPolynomial.hpp"
 
 #include <cstddef>
+#include <memory>
 #include <stdexcept>
+#include <utility>
 #include <vector>
 
 #include <catch2/catch_test_macros.hpp>
@@ -9,6 +11,41 @@
 using MP2 = mvPolynomial::MVPolynomial<int, double, 2>;
 using MP3 = mvPolynomial::MVPolynomial<int, double, 3>;
 
+template <typename T>
+class CountingAllocator {
+ public:
+  using value_type = T;
+
+  CountingAllocator() : allocations_(std::make_shared<std::size_t>(0)) {}
+
+  explicit CountingAllocator(std::shared_ptr<std::size_t> allocations) : allocations_(std::move(allocations)) {}
+
+  template <typename U>
+  CountingAllocator(const CountingAllocator<U>& other) : allocations_(other.allocations_) {}
+
+  [[nodiscard]] T* allocate(std::size_t n) {
+    *allocations_ += n;
+    return std::allocator<T>{}.allocate(n);
+  }
+
+  void deallocate(T* p, std::size_t n) noexcept { std::allocator<T>{}.deallocate(p, n); }
+
+  template <typename U>
+  bool operator==(const CountingAllocator<U>& other) const noexcept {
+    return allocations_ == other.allocations_;
+  }
+
+  template <typename U>
+  bool operator!=(const CountingAllocator<U>& other) const noexcept {
+    return !(*this == other);
+  }
+
+  std::shared_ptr<std::size_t> allocations_;
+};
+
+using CountingPairAllocator = CountingAllocator<std::pair<const Eigen::Array2i, double>>;
+using CountingMP2           = mvPolynomial::MVPolynomial<int, double, 2, CountingPairAllocator>;
+
 TEST_CASE("constructor", "[mvPolynomial]") {
   SECTION("default") {
     auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
@@ -186,6 +223,61 @@ TEST_CASE("pow", "[mvPolynomial]") {
   SECTION("15") { REQUIRE(m.pow(15) == m3 * m3 * m3 * m3 * m3); }
 }
 
+TEST_CASE("allocator propagation", "[mvPolynomial]") {
+  auto allocations = std::make_shared<std::size_t>(0);
+  auto alloc       = CountingPairAllocator(allocations);
+  auto p           = CountingMP2(
+      {
+          {{0, 0}, 1},
+          {{1, 0}, 2},
+          {{0, 1}, 3},
+  },
+      alloc
+  );
+
+  SECTION("pow uses rebound allocator") {
+    const auto before = *allocations;
+    auto       p7     = p.pow(7);
+
+    REQUIRE(*allocations > before);
+    REQUIRE(p7 == p.pow(7));
+  }
+
+  SECTION("point evaluation uses rebound allocator") {
+    const auto before = *allocations;
+    auto       value  = p(Eigen::Vector2d({2, 3}));
+
+    REQUIRE(*allocations > before);
+    REQUIRE(value == 14);
+  }
+
+  SECTION("polynomial composition uses rebound allocator") {
+    auto x = CountingMP2(
+        {
+            {{0, 0}, 1},
+            {{1, 0}, 2},
+            {{0, 1}, 3},
+    },
+        alloc
+    );
+
+    const auto before   = *allocations;
+    auto       composed = p(x, 1);
+
+    REQUIRE(*allocations > before);
+    REQUIRE(
+        composed == CountingMP2(
+                        {
+                            {{0, 0}, 4},
+                            {{1, 0}, 8},
+                            {{0, 1}, 9},
+    },
+                        alloc
+                    )
+    );
+  }
+}
+
 TEST_CASE("operator()", "[mvPolynomial]") {
   SECTION("point") {
     auto m = MP2({

From 2d5cf1d7421297df781cef14e4d3e9536920aec0 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 21:50:24 +0900
Subject: [PATCH 31/41] Update README.

---
 README.md | 63 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)

diff --git a/README.md b/README.md
index 61f4b02..f4e71cc 100644
--- a/README.md
+++ b/README.md
@@ -33,3 +33,66 @@ The examples exist in "test" directory.
 
 ## Multi-variable polynomial
 A class `MVPolynomial` implements multi-variable polynomials.
+
+```cpp
+#include "mvPolynomial/mvPolynomial.hpp"
+
+using MP2 = mvPolynomial::MVPolynomial<int, double, 2>;
+```
+
+You can construct a polynomial from an initializer list of `(index, coefficient)` pairs.
+
+```cpp
+auto p = MP2({
+    {{0, 0}, 1.0},
+    {{1, 0}, 2.0},
+    {{0, 1}, 3.0},
+});
+```
+
+The constant term is the zero index.
+
+```cpp
+double c = p.get(Eigen::Array2i::Zero());
+p.set(Eigen::Array2i({2, 0}), 5.0);
+```
+
+`set(index, 0.0)` removes the term, and zero-coefficient terms are normalized away automatically.
+
+The class provides basic polynomial arithmetic.
+
+```cpp
+auto q   = MP2({
+    {{0, 0}, 4.0},
+    {{1, 0}, 1.0},
+});
+auto sum = p + q;
+auto sub = p - q;
+auto mul = p * q;
+auto pw  = p.pow(3);
+```
+
+You can evaluate a polynomial at a point.
+
+```cpp
+double value = p(Eigen::Vector2d({2.0, 3.0}));
+```
+
+You can also substitute a polynomial for one axis.
+
+```cpp
+auto x = MP2({
+    {{0, 0}, 1.0},
+    {{1, 0}, 2.0},
+});
+auto composed = p(x, 1);
+```
+
+Differentiation and integration are free functions.
+
+```cpp
+auto dx = mvPolynomial::D(p, 0);
+auto iy = mvPolynomial::Integrate(p, 1);
+```
+
+Negative indices are not supported and cause `std::invalid_argument`.

From 5249a1f1e63067cad4547d00a9b8a380adbc19ff Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 22:18:47 +0900
Subject: [PATCH 32/41] Fix bugs and add random tests for operator().

---
 include/mvPolynomial/mvPolynomial.hpp |  12 +-
 test/mvPolynomial_test.cpp            | 156 ++++++++++++++++++++++++++
 2 files changed, 162 insertions(+), 6 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index d1e19bc..5c32421 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -274,16 +274,16 @@ class MVPolynomial final {
       const auto& next_index = next_it->first;
 
       index_type index_diff              = index - next_index;
-      index_type index_diff_without_axis = index_diff;
+      auto       index_diff_without_axis = index_diff;
       index_diff_without_axis[axis]      = 0;
 
       if ((index_diff > 0).all()) {
         auto tmp_mvp = MVPolynomial{{{index_diff_without_axis, 1}}, get_allocator()};
-        partial_sum  = next_coeff + partial_sum * (mvp.pow(index[axis])) * tmp_mvp;
+        partial_sum  = next_coeff + partial_sum * (mvp.pow(index_diff[axis])) * tmp_mvp;
       } else {
-        index_type index_without_axis = index;
-        index_without_axis[axis]      = 0;
-        auto tmp_mvp                  = MVPolynomial{{{index_without_axis, 1}}, get_allocator()};
+        auto index_without_axis  = index;
+        index_without_axis[axis] = 0;
+        auto tmp_mvp             = MVPolynomial{{{index_without_axis, 1}}, get_allocator()};
 
         partial_sum *= mvp.pow(index[axis]) * tmp_mvp;
         partial_sums.push_back(std::move(partial_sum));
@@ -292,7 +292,7 @@ class MVPolynomial final {
     }
     {
       const auto& first_index              = begin()->first;
-      index_type  first_index_without_axis = first_index;
+      auto        first_index_without_axis = first_index;
       first_index_without_axis[axis]       = 0;
       auto tmp_mvp                         = MVPolynomial{{{first_index_without_axis, 1}}, get_allocator()};
       partial_sum *= mvp.pow(first_index[axis]) * tmp_mvp;
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index db8395b..2e263a4 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -1,7 +1,10 @@
 #include "mvPolynomial/mvPolynomial.hpp"
 
+#include <cmath>
 #include <cstddef>
 #include <memory>
+#include <random>
+#include <sstream>
 #include <stdexcept>
 #include <utility>
 #include <vector>
@@ -11,6 +14,119 @@
 using MP2 = mvPolynomial::MVPolynomial<int, double, 2>;
 using MP3 = mvPolynomial::MVPolynomial<int, double, 3>;
 
+template <typename Poly>
+typename Poly::mapped_type NaivePointEval(const Poly& p, const typename Poly::coord_type& x) {
+  auto sum = typename Poly::mapped_type{0};
+  for (const auto& [index, coeff] : p) {
+    auto term = coeff;
+    for (int axis = 0; axis < Poly::dim; ++axis) {
+      term *= std::pow(x[axis], index[axis]);
+    }
+    sum += term;
+  }
+  return sum;
+}
+
+template <typename Poly>
+Poly NaiveCompose(const Poly& p, const Poly& mvp, int axis) {
+  auto composed = Poly(p.get_allocator());
+  for (const auto& [index, coeff] : p) {
+    auto index_without_axis  = index;
+    index_without_axis[axis] = 0;
+    composed += Poly(
+                    {
+                        {index_without_axis, coeff}
+    },
+                    p.get_allocator()
+                ) *
+                mvp.pow(index[axis]);
+  }
+  return composed;
+}
+
+template <typename Poly, typename URNG>
+Poly MakeRandomPolynomial(URNG& rng, int term_count, int max_degree, int max_abs_coeff) {
+  auto poly = Poly();
+
+  auto degree_dist = std::uniform_int_distribution<int>(0, max_degree);
+  auto coeff_dist  = std::uniform_int_distribution<int>(-max_abs_coeff, max_abs_coeff);
+  for (int term = 0; term < term_count; ++term) {
+    typename Poly::index_type index;
+    index.setZero();
+    for (int axis = 0; axis < Poly::dim; ++axis) {
+      index(axis) = degree_dist(rng);
+    }
+
+    auto coeff = coeff_dist(rng);
+    if (coeff == 0) {
+      continue;
+    }
+
+    if (poly.contains(index)) {
+      poly.set(index, poly.get(index) + coeff);
+    } else {
+      poly.set(index, coeff);
+    }
+  }
+
+  if (poly == Poly()) {
+    typename Poly::index_type index;
+    index.setZero();
+    index(0) = 1;
+    poly.set(index, 1);
+  }
+
+  return poly;
+}
+
+template <typename Poly, typename URNG>
+typename Poly::coord_type MakeRandomPoint(URNG& rng, int min_coord, int max_coord) {
+  typename Poly::coord_type point;
+  point.setZero();
+  auto coord_dist = std::uniform_int_distribution<int>(min_coord, max_coord);
+  for (int axis = 0; axis < Poly::dim; ++axis) {
+    point(axis) = static_cast<typename Poly::mapped_type>(coord_dist(rng));
+  }
+  return point;
+}
+
+template <typename Poly>
+std::string DescribePolynomial(const Poly& p) {
+  auto oss   = std::ostringstream{};
+  auto first = true;
+  oss << "{";
+  for (const auto& [index, coeff] : p) {
+    if (!first) {
+      oss << ", ";
+    }
+    first = false;
+    oss << "[";
+    for (int axis = 0; axis < Poly::dim; ++axis) {
+      if (axis > 0) {
+        oss << ",";
+      }
+      oss << index[axis];
+    }
+    oss << "]=" << coeff;
+  }
+  oss << "}";
+  return oss.str();
+}
+
+template <typename Poly>
+std::string DescribePoint(const typename Poly::coord_type& point) {
+  auto oss = std::ostringstream{};
+  oss << "[";
+  for (int axis = 0; axis < Poly::dim; ++axis) {
+    if (axis > 0) {
+      oss << ",";
+    }
+    oss << point[axis];
+  }
+  oss << "]";
+  return oss.str();
+}
+
 template <typename T>
 class CountingAllocator {
  public:
@@ -397,6 +513,46 @@ TEST_CASE("operator()", "[mvPolynomial]") {
     })
     );
   }
+
+  SECTION("point_randomized_matches_naive_sum") {
+    auto rng             = std::mt19937(20260517);
+    auto term_count_dist = std::uniform_int_distribution<int>(1, 12);
+
+    for (int trial = 0; trial < 200; ++trial) {
+      auto m        = MakeRandomPolynomial<MP3>(rng, term_count_dist(rng), 4, 5);
+      auto x        = MakeRandomPoint<MP3>(rng, -3, 3);
+      auto actual   = m(x);
+      auto expected = NaivePointEval(m, x);
+
+      CAPTURE(trial, DescribePolynomial(m), DescribePoint<MP3>(x));
+      REQUIRE(actual == expected);
+    }
+  }
+
+  SECTION("composition_randomized_matches_naive_sum") {
+    auto rng                   = std::mt19937(20260518);
+    auto outer_term_count_dist = std::uniform_int_distribution<int>(1, 10);
+    auto inner_term_count_dist = std::uniform_int_distribution<int>(1, 6);
+    auto axis_dist             = std::uniform_int_distribution<int>(0, MP3::dim - 1);
+
+    for (int trial = 0; trial < 120; ++trial) {
+      auto outer    = MakeRandomPolynomial<MP3>(rng, outer_term_count_dist(rng), 3, 4);
+      auto inner    = MakeRandomPolynomial<MP3>(rng, inner_term_count_dist(rng), 2, 3);
+      auto axis     = axis_dist(rng);
+      auto actual   = outer(inner, axis);
+      auto expected = NaiveCompose(outer, inner, axis);
+
+      CAPTURE(
+          trial,
+          axis,
+          DescribePolynomial(outer),
+          DescribePolynomial(inner),
+          DescribePolynomial(actual),
+          DescribePolynomial(expected)
+      );
+      REQUIRE(actual == expected);
+    }
+  }
 }
 
 TEST_CASE("D", "[mvPolynomial]") {

From 8db2044bbfe139db438f0f2a0cd41ae4d24cea4f Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 22:54:31 +0900
Subject: [PATCH 33/41] Unify responsibility for normalization.

---
 include/mvPolynomial/mvPolynomial.hpp | 226 ++++++++++++++------------
 test/mvPolynomial_test.cpp            |  95 ++++++++++-
 2 files changed, 210 insertions(+), 111 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 5c32421..8c60ebe 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -13,6 +13,7 @@
 #include <limits>
 #include <memory>
 #include <numeric>
+#include <optional>
 #include <ranges>
 #include <stdexcept>
 #include <type_traits>
@@ -33,23 +34,23 @@ inline void CheckAxis(int dim, int axis) {
 template <
     std::signed_integral IntType,
     std::floating_point  R,
-    int                  D,
-    class Allocator = std::allocator<std::pair<const IndexType<IntType, D>, R>>>
+    int                  dim_,
+    class Allocator = std::allocator<std::pair<const IndexType<IntType, dim_>, R>>>
 class MVPolynomial final {
  public:
-  static_assert(D > 0, "MVPolynomial: the dimension must be greater than 0.");
+  static_assert(dim_ > 0, "MVPolynomial: the dimension must be greater than 0.");
 
-  static constexpr int dim = D;
+  static constexpr int dim = dim_;
 
   // This setting is too strict, so I expect users to set tolerance.
   inline static R rel_tolerance = std::numeric_limits<R>::epsilon();
   inline static R abs_tolerance = std::numeric_limits<R>::min();
 
-  using index_type = IndexType<IntType, D>;
+  using index_type = IndexType<IntType, dim>;
   using coord_type = CoordType<R, dim>;
 
  private:
-  using IndexContainer = platanus::btree_map<index_type, R, IndexComparer<IntType, D>, Allocator>;
+  using IndexContainer = platanus::btree_map<index_type, R, IndexComparer<IntType, dim>, Allocator>;
 
  public:
   using key_type    = IndexContainer::key_type;
@@ -90,27 +91,27 @@ class MVPolynomial final {
   MVPolynomial& operator=(MVPolynomial&& other)      = default;
   ~MVPolynomial()                                    = default;
 
-  MVPolynomial() : index2value_(allocator_type{}) { AssignCoeffWithoutNormalization(index_type::Zero(), R{0}); }
+  MVPolynomial() : index2value_(allocator_type{}) { AssignCoeffRaw(index_type::Zero(), R{0}); }
 
   explicit MVPolynomial(const allocator_type& allocator) : index2value_(allocator) {
-    AssignCoeffWithoutNormalization(index_type::Zero(), R{0});
+    AssignCoeffRaw(index_type::Zero(), R{0});
   }
 
   template <typename InputIterator>
   MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e, allocator_type{}) {
     CheckSelfIndexes();
-    DeleteZeroCoeffTerm();
+    Normalize();
   }
 
   template <typename InputIterator>
   MVPolynomial(InputIterator s, InputIterator e, const allocator_type& allocator) : index2value_(s, e, allocator) {
     CheckSelfIndexes();
-    DeleteZeroCoeffTerm();
+    Normalize();
   }
 
   MVPolynomial(std::initializer_list<value_type> l, const allocator_type& a = allocator_type{}) : index2value_(l, a) {
     CheckSelfIndexes();
-    DeleteZeroCoeffTerm();
+    Normalize();
   }
 
   MVPolynomial(const MVPolynomial& m, const allocator_type& a) : index2value_(m.index2value_, a) {}
@@ -121,13 +122,13 @@ class MVPolynomial final {
     index2value_.clear();
     index2value_.insert(l.begin(), l.end());
     CheckSelfIndexes();
-    DeleteZeroCoeffTerm();
+    Normalize();
     return *this;
   }
 
   MVPolynomial(mapped_type r, const allocator_type& a = allocator_type{}) : index2value_(a) {
-    AssignCoeffWithoutNormalization(index_type::Zero(), r);
-    DeleteZeroCoeffTerm();
+    AssignCoeffRaw(index_type::Zero(), r);
+    Normalize();
   }
 
   allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); }
@@ -151,11 +152,16 @@ class MVPolynomial final {
 
   const mapped_type& get(const key_type& index) const { return index2value_.at(index); }
 
-  void set(const key_type& index, R coeff) {
-    AssignCoeffWithoutNormalization(index, coeff);
-    DeleteZeroCoeffTerm();
+  std::optional<const mapped_type&> try_get(const key_type& index) const {
+    if (contains(index)) {
+      return index2value_[index];
+    } else {
+      return std::nullopt;
+    }
   }
 
+  void set(const key_type& index, R coeff) { AssignCoeff(index, coeff); }
+
   void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); }
 
   const_iterator find(const key_type& i) const { return index2value_.find(i); }
@@ -219,22 +225,6 @@ class MVPolynomial final {
     }
   }
 
-  void DeleteZeroCoeffTerm() {
-    auto removed_term_indexes = rebound_vector_type<size_type>(rebound_allocator_type<size_type>(get_allocator()));
-    for (size_t i = 0; const auto& index_and_coeff : index2value_) {
-      if (std::abs(index_and_coeff.second) < abs_tolerance) {
-        removed_term_indexes.push_back(i);
-      }
-      ++i;
-    }
-    for (auto removed_index : removed_term_indexes | std::views::reverse) {
-      index2value_.erase(std::next(index2value_.begin(), removed_index));
-    }
-    if (index2value_.empty()) {
-      AssignCoeffWithoutNormalization(index_type::Zero(), R{0});
-    }
-  }
-
   R operator()(const coord_type& x) const {
     auto partial_sums = rebound_vector_type<R>(rebound_allocator_type<R>(get_allocator()));
     auto partial_sum  = rbegin()->second;
@@ -300,6 +290,52 @@ class MVPolynomial final {
     return std::reduce(partial_sums.cbegin(), partial_sums.cend()) + partial_sum;
   }
 
+  MVPolynomial D(int axis) const {
+    details::CheckAxis(dim, axis);
+
+    auto dp = MVPolynomial{get_allocator()};
+    dp.index2value_.clear();
+    auto it = begin();
+    while (it != end()) {
+      auto       value = it->second;
+      index_type index = it->first;
+      if (index[axis] == 0) {
+        auto d_end_it = end();
+        for (int ith_axis = 0; ith_axis <= axis; ++ith_axis) {
+          d_end_it = std::partition_point(it, d_end_it, [ith_axis, &index](const value_type& v) {
+            return v.first[ith_axis] == index[ith_axis];
+          });
+        }
+        // Skip indexes which axis-th element is zero.
+        it = d_end_it;
+      } else {
+        value *= index[axis]--;
+        dp.AssignCoeffRaw(index, value);
+        ++it;
+      }
+    }
+    dp.Normalize();
+
+    return dp;
+  }
+
+  MVPolynomial Integrate(int axis) const {
+    details::CheckAxis(dim, axis);
+
+    auto result = MVPolynomial(get_allocator());
+    result.index2value_.clear();
+    for (const auto& [index, value] : *this) {
+      auto new_index = index;
+      ++new_index[axis];
+      const auto new_value = value / new_index[axis];
+      result.AddCoeffRaw(new_index, new_value);
+    }
+
+    result.Normalize();
+
+    return result;
+  }
+
   MVPolynomial operator+() const { return *this; }
 
   MVPolynomial operator-() const& {
@@ -320,30 +356,30 @@ class MVPolynomial final {
   }
 
   MVPolynomial& operator+=(mapped_type r) {
-    AddCoeffWithoutNormalization(index_type::Zero(), r);
-    DeleteZeroCoeffTerm();
+    AddCoeff(index_type::Zero(), r);
     return *this;
   }
 
   MVPolynomial& operator+=(const MVPolynomial& r) {
-    for (const auto& [idx, coeff] : r) {
-      AddCoeffWithoutNormalization(idx, coeff);
+    auto        copied_r = (&r == this) ? std::optional<MVPolynomial>(MVPolynomial(r, get_allocator())) : std::nullopt;
+    const auto& rhs      = copied_r ? *copied_r : r;
+    for (const auto& [idx, coeff] : rhs) {
+      AddCoeff(idx, coeff);
     }
-    DeleteZeroCoeffTerm();
     return *this;
   }
 
   MVPolynomial& operator-=(mapped_type r) {
-    AddCoeffWithoutNormalization(index_type::Zero(), -r);
-    DeleteZeroCoeffTerm();
+    AddCoeff(index_type::Zero(), -r);
     return *this;
   }
 
   MVPolynomial& operator-=(const MVPolynomial& r) {
-    for (const auto& [idx, coeff] : r) {
-      AddCoeffWithoutNormalization(idx, -coeff);
+    auto        copied_r = (&r == this) ? std::optional<MVPolynomial>(MVPolynomial(r, get_allocator())) : std::nullopt;
+    const auto& rhs      = copied_r ? *copied_r : r;
+    for (const auto& [idx, coeff] : rhs) {
+      AddCoeff(idx, -coeff);
     }
-    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -352,7 +388,7 @@ class MVPolynomial final {
       auto& coeff = index_and_coeff.second;
       coeff *= r;
     }
-    DeleteZeroCoeffTerm();
+    Normalize();
     return *this;
   }
 
@@ -364,13 +400,13 @@ class MVPolynomial final {
       for (const auto& [index, coeff] : *this) {
         const auto new_index = index + r_index;
         const auto new_coeff = coeff * r_coeff;
-        result.AddCoeffWithoutNormalization(new_index, new_coeff);
+        result.AddCoeffRaw(new_index, new_coeff);
       }
+      result.Normalize();
       swap(result);
     } else {
       *this = *this * r;
     }
-    DeleteZeroCoeffTerm();
     return *this;
   }
 
@@ -443,22 +479,52 @@ class MVPolynomial final {
         const auto& [r_idx, r_v] = r_p;
         const auto idx           = l_idx + r_idx;
         const auto v             = l_v * r_v;
-        mul.AddCoeffWithoutNormalization(idx, v);
+        mul.AddCoeffRaw(idx, v);
       }
     }
 
-    mul.DeleteZeroCoeffTerm();
+    mul.Normalize();
 
     return mul;
   }
 
  private:
-  void AssignCoeffWithoutNormalization(const key_type& index, mapped_type coeff) {
+  void Normalize() {
+    auto removed_term_indexes = rebound_vector_type<size_type>(rebound_allocator_type<size_type>(get_allocator()));
+    for (size_t i = 0; const auto& index_and_coeff : index2value_) {
+      if (std::abs(index_and_coeff.second) < abs_tolerance) {
+        removed_term_indexes.push_back(i);
+      }
+      ++i;
+    }
+    for (auto removed_index : removed_term_indexes | std::views::reverse) {
+      index2value_.erase(std::next(index2value_.begin(), removed_index));
+    }
+    if (index2value_.empty()) {
+      AssignCoeffRaw(index_type::Zero(), R{0});
+    }
+  }
+
+  void AssignCoeff(const key_type& index, mapped_type coeff) {
+    AssignCoeffRaw(index, coeff);
+    Normalize();
+  }
+
+  void AddCoeff(const key_type& index, mapped_type delta) {
+    AddCoeffRaw(index, delta);
+    Normalize();
+  }
+
+  void AssignCoeffRaw(const key_type& index, mapped_type coeff) {
     CheckIndexIncludingNegative(index);
-    index2value_[index] = coeff;
+    if (auto it = index2value_.find(index); it != index2value_.end()) {
+      it->second = coeff;
+    } else {
+      index2value_[index] = coeff;
+    }
   }
 
-  void AddCoeffWithoutNormalization(const key_type& index, mapped_type delta) {
+  void AddCoeffRaw(const key_type& index, mapped_type delta) {
     CheckIndexIncludingNegative(index);
     if (auto it = index2value_.find(index); it != index2value_.end()) {
       it->second += delta;
@@ -482,60 +548,14 @@ class MVPolynomial final {
   IndexContainer index2value_;
 };
 
-template <std::signed_integral IntType, std::floating_point R, int Dim, class Allocator>
-auto D(const MVPolynomial<IntType, R, Dim, Allocator>& p, int axis) {
-  using MP            = MVPolynomial<IntType, R, Dim, Allocator>;
-  using Index         = typename MP::index_type;
-  using IndexAndCoeff = typename MP::value_type;
-
-  details::CheckAxis(MP::dim, axis);
-
-  auto dp   = MP{p.get_allocator()};
-  auto p_it = p.begin();
-  while (p_it != p.end()) {
-    auto value = p_it->second;
-    auto index = Index{p_it->first};
-    if (index[axis] == 0) {
-      auto d_end_it = p.end();
-      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]--;
-      dp.set(index, value);
-      ++p_it;
-    }
-  }
-  dp.DeleteZeroCoeffTerm();
-
-  return dp;
+template <std::signed_integral IntType, std::floating_point R, int dim, class Allocator>
+auto D(const MVPolynomial<IntType, R, dim, Allocator>& p, int axis) {
+  return p.D(axis);
 }
 
-template <std::signed_integral IntType, std::floating_point R, int D, class Allocator>
-auto Integrate(MVPolynomial<IntType, R, D, Allocator> p, int axis) {
-  using MP = MVPolynomial<IntType, R, D, Allocator>;
-
-  details::CheckAxis(D, axis);
-
-  auto result = MP(p.get_allocator());
-  for (const auto& [index, value] : p) {
-    auto new_index = index;
-    ++new_index[axis];
-    const auto new_value = value / new_index[axis];
-    if (result.contains(new_index)) {
-      result.set(new_index, result.get(new_index) + new_value);
-    } else {
-      result.set(new_index, new_value);
-    }
-  }
-
-  result.DeleteZeroCoeffTerm();
-
-  return result;
+template <std::signed_integral IntType, std::floating_point R, int dim, class Allocator>
+auto Integrate(MVPolynomial<IntType, R, dim, Allocator> p, int axis) {
+  return p.Integrate(axis);
 }
 
 }  // namespace mvPolynomial
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 2e263a4..e9a814e 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -264,6 +264,13 @@ TEST_CASE("constructor", "[mvPolynomial]") {
 }
 
 TEST_CASE("invariants", "[mvPolynomial]") {
+  auto require_canonical_zero = [](const MP2& p) {
+    REQUIRE(p.size() == 1);
+    REQUIRE(p == MP2());
+    REQUIRE(p.contains(Eigen::Array2i::Zero()));
+    REQUIRE(p.get(Eigen::Array2i::Zero()) == 0);
+  };
+
   SECTION("subtracting self stays canonical zero polynomial") {
     auto p = MP2({
         {{0, 0}, 1},
@@ -273,10 +280,7 @@ TEST_CASE("invariants", "[mvPolynomial]") {
 
     p -= p;
 
-    REQUIRE(p.size() == 1);
-    REQUIRE(p == MP2());
-    REQUIRE(p.contains(Eigen::Array2i::Zero()));
-    REQUIRE(p.get(Eigen::Array2i::Zero()) == 0);
+    require_canonical_zero(p);
   }
 
   SECTION("multiplying by zero stays canonical zero polynomial") {
@@ -288,10 +292,7 @@ TEST_CASE("invariants", "[mvPolynomial]") {
 
     p *= 0.0;
 
-    REQUIRE(p.size() == 1);
-    REQUIRE(p == MP2());
-    REQUIRE(p.contains(Eigen::Array2i::Zero()));
-    REQUIRE(p.get(Eigen::Array2i::Zero()) == 0);
+    require_canonical_zero(p);
   }
 
   SECTION("set rejects negative index for lvalue and rvalue keys") {
@@ -315,6 +316,46 @@ TEST_CASE("invariants", "[mvPolynomial]") {
     REQUIRE_FALSE(p.contains(Eigen::Array2i({1, 0})));
     REQUIRE(p.get(Eigen::Array2i::Zero()) == 1);
   }
+
+  SECTION("scalar addition and subtraction preserve canonical zero") {
+    auto p = MP2(2.0);
+
+    p += -2.0;
+    require_canonical_zero(p);
+
+    p -= 0.0;
+    require_canonical_zero(p);
+  }
+
+  SECTION("polynomial addition preserves canonical zero after cancellation") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+    });
+    auto q = MP2({
+        {{0, 0}, -1},
+        {{1, 0}, -2},
+    });
+
+    p += q;
+
+    require_canonical_zero(p);
+  }
+
+  SECTION("polynomial subtraction preserves canonical zero after cancellation") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{0, 1}, 3},
+    });
+    auto q = MP2({
+        {{0, 0}, 1},
+        {{0, 1}, 3},
+    });
+
+    p -= q;
+
+    require_canonical_zero(p);
+  }
 }
 
 TEST_CASE("pow", "[mvPolynomial]") {
@@ -592,6 +633,16 @@ TEST_CASE("D", "[mvPolynomial]") {
       REQUIRE(dm1.get(ans[i].first) == ans[i].second);
     }
   }
+
+  SECTION("zero result stays canonical zero polynomial") {
+    auto constant = MP2(7.0);
+    auto dx       = mvPolynomial::D(constant, 0);
+
+    REQUIRE(dx.size() == 1);
+    REQUIRE(dx == MP2());
+    REQUIRE(dx.contains(Eigen::Array2i::Zero()));
+    REQUIRE(dx.get(Eigen::Array2i::Zero()) == 0);
+  }
 }
 
 TEST_CASE("integral", "[mvPolynomial]") {
@@ -665,6 +716,16 @@ TEST_CASE("integral", "[mvPolynomial]") {
     REQUIRE((shi.lower_bound(Eigen::Array2i({4, 2}))->first == Eigen::Array2i({4, 2})).all());
     REQUIRE(shi.get(Eigen::Array2i({4, 2})) == 2);
   }
+
+  SECTION("zero polynomial stays canonical zero after integration") {
+    auto zero = MP2();
+    auto sz   = mvPolynomial::Integrate(zero, 0);
+
+    REQUIRE(sz.size() == 1);
+    REQUIRE(sz == MP2());
+    REQUIRE(sz.contains(Eigen::Array2i::Zero()));
+    REQUIRE(sz.get(Eigen::Array2i::Zero()) == 0);
+  }
 }
 
 TEST_CASE("multiply", "[mvPolynomial]") {
@@ -722,6 +783,24 @@ TEST_CASE("multiply", "[mvPolynomial]") {
     REQUIRE(inplace.lower_bound(Eigen::Array2i({2, 2})) != inplace.end());
     REQUIRE((inplace.lower_bound(Eigen::Array2i({2, 2}))->first == Eigen::Array2i({2, 2})).all());
   }
+
+  SECTION("multiplication normalizes zero coefficient collisions") {
+    auto cancel_l = MP2({
+        {{0, 0},  1},
+        {{1, 0}, -1},
+    });
+    auto cancel_r = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 1},
+    });
+
+    auto cancel_prod = cancel_l * cancel_r;
+
+    REQUIRE(cancel_prod.size() == 2);
+    REQUIRE(cancel_prod.get(Eigen::Array2i({0, 0})) == 1);
+    REQUIRE_FALSE(cancel_prod.contains(Eigen::Array2i({1, 0})));
+    REQUIRE(cancel_prod.get(Eigen::Array2i({2, 0})) == -1);
+  }
 }
 
 TEST_CASE("sum", "[mvPolynomial]") {

From ba7b5d9d1e30d906f26070f3ec44679cb58835e0 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 23:10:21 +0900
Subject: [PATCH 34/41] Add a func normalizing single term.

---
 include/mvPolynomial/mvPolynomial.hpp | 13 +++++-
 test/mvPolynomial_test.cpp            | 61 +++++++++++++++++++++++----
 2 files changed, 63 insertions(+), 11 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 8c60ebe..95fc8ca 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -505,14 +505,23 @@ class MVPolynomial final {
     }
   }
 
+  void NormalizeSingle(const key_type& index) {
+    if (std::abs(get(index)) < abs_tolerance) {
+      index2value_.erase(index);
+    }
+    if (index2value_.empty()) {
+      AssignCoeffRaw(index_type::Zero(), R{0});
+    }
+  }
+
   void AssignCoeff(const key_type& index, mapped_type coeff) {
     AssignCoeffRaw(index, coeff);
-    Normalize();
+    NormalizeSingle(index);
   }
 
   void AddCoeff(const key_type& index, mapped_type delta) {
     AddCoeffRaw(index, delta);
-    Normalize();
+    NormalizeSingle(index);
   }
 
   void AssignCoeffRaw(const key_type& index, mapped_type coeff) {
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index e9a814e..8d5e8d7 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -14,6 +14,12 @@
 using MP2 = mvPolynomial::MVPolynomial<int, double, 2>;
 using MP3 = mvPolynomial::MVPolynomial<int, double, 3>;
 
+template <typename Poly>
+bool IsCanonicalZeroPolynomial(const Poly& p) {
+  auto zero_index = Poly::index_type::Zero();
+  return p.size() == 1 && p.contains(zero_index) && p.get(zero_index) == typename Poly::mapped_type{0};
+}
+
 template <typename Poly>
 typename Poly::mapped_type NaivePointEval(const Poly& p, const typename Poly::coord_type& x) {
   auto sum = typename Poly::mapped_type{0};
@@ -46,7 +52,7 @@ Poly NaiveCompose(const Poly& p, const Poly& mvp, int axis) {
 
 template <typename Poly, typename URNG>
 Poly MakeRandomPolynomial(URNG& rng, int term_count, int max_degree, int max_abs_coeff) {
-  auto poly = Poly();
+  auto terms = std::vector<std::pair<typename Poly::index_type, typename Poly::mapped_type>>{};
 
   auto degree_dist = std::uniform_int_distribution<int>(0, max_degree);
   auto coeff_dist  = std::uniform_int_distribution<int>(-max_abs_coeff, max_abs_coeff);
@@ -62,21 +68,44 @@ Poly MakeRandomPolynomial(URNG& rng, int term_count, int max_degree, int max_abs
       continue;
     }
 
-    if (poly.contains(index)) {
-      poly.set(index, poly.get(index) + coeff);
-    } else {
-      poly.set(index, coeff);
+    auto found = false;
+    for (auto& [existing_index, existing_coeff] : terms) {
+      if ((existing_index == index).all()) {
+        existing_coeff += static_cast<typename Poly::mapped_type>(coeff);
+        found = true;
+        break;
+      }
+    }
+
+    if (!found) {
+      terms.emplace_back(index, static_cast<typename Poly::mapped_type>(coeff));
     }
   }
 
-  if (poly == Poly()) {
+  auto filtered_terms = std::vector<std::pair<typename Poly::index_type, typename Poly::mapped_type>>{};
+  filtered_terms.reserve(terms.size());
+  for (const auto& [index, coeff] : terms) {
+    if (coeff != typename Poly::mapped_type{0}) {
+      filtered_terms.emplace_back(index, coeff);
+    }
+  }
+
+  if (filtered_terms.empty()) {
     typename Poly::index_type index;
     index.setZero();
-    index(0) = 1;
-    poly.set(index, 1);
+    for (int axis = 0; axis < Poly::dim; ++axis) {
+      index(axis) = degree_dist(rng);
+    }
+
+    auto coeff = 0;
+    do {
+      coeff = coeff_dist(rng);
+    } while (coeff == 0);
+
+    filtered_terms.emplace_back(index, static_cast<typename Poly::mapped_type>(coeff));
   }
 
-  return poly;
+  return Poly(filtered_terms.begin(), filtered_terms.end());
 }
 
 template <typename Poly, typename URNG>
@@ -380,6 +409,20 @@ TEST_CASE("pow", "[mvPolynomial]") {
   SECTION("15") { REQUIRE(m.pow(15) == m3 * m3 * m3 * m3 * m3); }
 }
 
+TEST_CASE("random polynomial helper", "[mvPolynomial]") {
+  SECTION("never returns canonical zero polynomial") {
+    auto rng             = std::mt19937(20260519);
+    auto term_count_dist = std::uniform_int_distribution<int>(1, 12);
+
+    for (int trial = 0; trial < 200; ++trial) {
+      auto p = MakeRandomPolynomial<MP3>(rng, term_count_dist(rng), 4, 5);
+
+      CAPTURE(trial, DescribePolynomial(p));
+      REQUIRE_FALSE(IsCanonicalZeroPolynomial(p));
+    }
+  }
+}
+
 TEST_CASE("allocator propagation", "[mvPolynomial]") {
   auto allocations = std::make_shared<std::size_t>(0);
   auto alloc       = CountingPairAllocator(allocations);

From e6d0f1a1c2d62f7ed20e14add9444c974abbb2ab Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 23:21:59 +0900
Subject: [PATCH 35/41] Fix try_get and optimize.

---
 include/mvPolynomial/mvPolynomial.hpp | 37 ++++++++++++---------------
 test/mvPolynomial_test.cpp            | 27 +++++++++++++++++++
 2 files changed, 44 insertions(+), 20 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 95fc8ca..9d3ea7b 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -14,7 +14,6 @@
 #include <memory>
 #include <numeric>
 #include <optional>
-#include <ranges>
 #include <stdexcept>
 #include <type_traits>
 #include <vector>
@@ -150,17 +149,16 @@ class MVPolynomial final {
   size_type max_size() const noexcept { return index2value_.max_size(); }
   size_type capacity() const noexcept { return index2value_.capacity(); }
 
-  const mapped_type& get(const key_type& index) const { return index2value_.at(index); }
+  mapped_type get(const key_type& index) const { return index2value_.at(index); }
 
-  std::optional<const mapped_type&> try_get(const key_type& index) const {
-    if (contains(index)) {
-      return index2value_[index];
-    } else {
-      return std::nullopt;
+  std::optional<mapped_type> try_get(const key_type& index) const {
+    if (auto it = index2value_.find(index); it != index2value_.end()) {
+      return it->second;
     }
+    return std::nullopt;
   }
 
-  void set(const key_type& index, R coeff) { AssignCoeff(index, coeff); }
+  void set(const key_type& index, mapped_type coeff) { AssignCoeff(index, coeff); }
 
   void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); }
 
@@ -364,8 +362,9 @@ class MVPolynomial final {
     auto        copied_r = (&r == this) ? std::optional<MVPolynomial>(MVPolynomial(r, get_allocator())) : std::nullopt;
     const auto& rhs      = copied_r ? *copied_r : r;
     for (const auto& [idx, coeff] : rhs) {
-      AddCoeff(idx, coeff);
+      AddCoeffRaw(idx, coeff);
     }
+    Normalize();
     return *this;
   }
 
@@ -378,8 +377,9 @@ class MVPolynomial final {
     auto        copied_r = (&r == this) ? std::optional<MVPolynomial>(MVPolynomial(r, get_allocator())) : std::nullopt;
     const auto& rhs      = copied_r ? *copied_r : r;
     for (const auto& [idx, coeff] : rhs) {
-      AddCoeff(idx, -coeff);
+      AddCoeffRaw(idx, -coeff);
     }
+    Normalize();
     return *this;
   }
 
@@ -490,19 +490,16 @@ class MVPolynomial final {
 
  private:
   void Normalize() {
-    auto removed_term_indexes = rebound_vector_type<size_type>(rebound_allocator_type<size_type>(get_allocator()));
-    for (size_t i = 0; const auto& index_and_coeff : index2value_) {
-      if (std::abs(index_and_coeff.second) < abs_tolerance) {
-        removed_term_indexes.push_back(i);
+    auto normalized = IndexContainer(get_allocator());
+    for (const auto& [index, coeff] : index2value_) {
+      if (std::abs(coeff) >= abs_tolerance) {
+        normalized[index] = coeff;
       }
-      ++i;
     }
-    for (auto removed_index : removed_term_indexes | std::views::reverse) {
-      index2value_.erase(std::next(index2value_.begin(), removed_index));
-    }
-    if (index2value_.empty()) {
-      AssignCoeffRaw(index_type::Zero(), R{0});
+    if (normalized.empty()) {
+      normalized[index_type::Zero()] = R{0};
     }
+    index2value_.swap(normalized);
   }
 
   void NormalizeSingle(const key_type& index) {
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 8d5e8d7..067d687 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -300,6 +300,33 @@ TEST_CASE("invariants", "[mvPolynomial]") {
     REQUIRE(p.get(Eigen::Array2i::Zero()) == 0);
   };
 
+  SECTION("try_get returns value for existing term and nullopt otherwise") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+
+    auto constant = p.try_get(Eigen::Array2i::Zero());
+    auto x_term   = p.try_get(Eigen::Array2i({1, 0}));
+    auto missing  = p.try_get(Eigen::Array2i({2, 2}));
+
+    REQUIRE(constant.has_value());
+    REQUIRE(*constant == 1);
+    REQUIRE(x_term.has_value());
+    REQUIRE(*x_term == 2);
+    REQUIRE_FALSE(missing.has_value());
+  }
+
+  SECTION("try_get sees canonical zero term on zero polynomial") {
+    auto zero = MP2();
+    auto term = zero.try_get(Eigen::Array2i::Zero());
+
+    REQUIRE(term.has_value());
+    REQUIRE(*term == 0);
+    REQUIRE_FALSE(zero.try_get(Eigen::Array2i({1, 0})).has_value());
+  }
+
   SECTION("subtracting self stays canonical zero polynomial") {
     auto p = MP2({
         {{0, 0}, 1},

From 386aa033692d7876609e6e7f39ec11c3347470b4 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 23:25:47 +0900
Subject: [PATCH 36/41] Remove copy.

---
 include/mvPolynomial/mvPolynomial.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 9d3ea7b..fa338ed 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -560,7 +560,7 @@ auto D(const MVPolynomial<IntType, R, dim, Allocator>& p, int axis) {
 }
 
 template <std::signed_integral IntType, std::floating_point R, int dim, class Allocator>
-auto Integrate(MVPolynomial<IntType, R, dim, Allocator> p, int axis) {
+auto Integrate(const MVPolynomial<IntType, R, dim, Allocator>& p, int axis) {
   return p.Integrate(axis);
 }
 

From 6cb9b666e187fb79e6a8f78987bd6acb6c9a80e0 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 23:27:09 +0900
Subject: [PATCH 37/41] Update README.

---
 README.md | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index f4e71cc..f9e9ed8 100644
--- a/README.md
+++ b/README.md
@@ -57,6 +57,15 @@ double c = p.get(Eigen::Array2i::Zero());
 p.set(Eigen::Array2i({2, 0}), 5.0);
 ```
 
+If you want to read a term that may be absent, use `try_get`.
+
+```cpp
+auto maybe_xy = p.try_get(Eigen::Array2i({1, 1}));
+if (maybe_xy.has_value()) {
+  double coeff = *maybe_xy;
+}
+```
+
 `set(index, 0.0)` removes the term, and zero-coefficient terms are normalized away automatically.
 
 The class provides basic polynomial arithmetic.
@@ -88,9 +97,12 @@ auto x = MP2({
 auto composed = p(x, 1);
 ```
 
-Differentiation and integration are free functions.
+Differentiation and integration are available both as member functions and as free functions.
 
 ```cpp
+auto dx_member = p.D(0);
+auto iy_member = p.Integrate(1);
+
 auto dx = mvPolynomial::D(p, 0);
 auto iy = mvPolynomial::Integrate(p, 1);
 ```

From 7aac87980cbb21e7b29234ab3d044a5c0ad49230 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 23:32:04 +0900
Subject: [PATCH 38/41] Add binary op with scalar.

---
 include/mvPolynomial/mvPolynomial.hpp | 52 ++++++++++++++++++
 test/mvPolynomial_test.cpp            | 76 +++++++++++++++++++++++++++
 2 files changed, 128 insertions(+)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index fa338ed..9502766 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -392,6 +392,16 @@ class MVPolynomial final {
     return *this;
   }
 
+  MVPolynomial& operator/=(mapped_type r) {
+    CheckScalarDivisor(r);
+    for (auto& index_and_coeff : index2value_) {
+      auto& coeff = index_and_coeff.second;
+      coeff /= r;
+    }
+    Normalize();
+    return *this;
+  }
+
   MVPolynomial& operator*=(const MVPolynomial& r) {
     if (r.size() == 1) {
       const auto& [r_index, r_coeff] = *(r.begin());
@@ -452,6 +462,16 @@ class MVPolynomial final {
     return std::move(l);
   }
 
+  friend MVPolynomial operator+(MVPolynomial l, mapped_type r) {
+    l += r;
+    return l;
+  }
+
+  friend MVPolynomial operator+(mapped_type l, MVPolynomial r) {
+    r += l;
+    return r;
+  }
+
   friend MVPolynomial operator+(const MVPolynomial& l, MVPolynomial&& r) { return std::move(r) + l; }
 
   friend MVPolynomial operator+(MVPolynomial&& l, MVPolynomial&& r) { return std::move(l) + r; }
@@ -465,6 +485,17 @@ class MVPolynomial final {
     return std::move(l);
   }
 
+  friend MVPolynomial operator-(MVPolynomial l, mapped_type r) {
+    l -= r;
+    return l;
+  }
+
+  friend MVPolynomial operator-(mapped_type l, const MVPolynomial& r) {
+    auto result = MVPolynomial(l, r.get_allocator());
+    result -= r;
+    return result;
+  }
+
   friend MVPolynomial operator-(const MVPolynomial& l, MVPolynomial&& r) { return -std::move(r) + l; }
 
   friend MVPolynomial operator-(MVPolynomial&& l, MVPolynomial&& r) { return std::move(l) - r; }
@@ -488,6 +519,21 @@ class MVPolynomial final {
     return mul;
   }
 
+  friend MVPolynomial operator*(MVPolynomial l, mapped_type r) {
+    l *= r;
+    return l;
+  }
+
+  friend MVPolynomial operator*(mapped_type l, MVPolynomial r) {
+    r *= l;
+    return r;
+  }
+
+  friend MVPolynomial operator/(MVPolynomial l, mapped_type r) {
+    l /= r;
+    return l;
+  }
+
  private:
   void Normalize() {
     auto normalized = IndexContainer(get_allocator());
@@ -551,6 +597,12 @@ class MVPolynomial final {
     }
   }
 
+  void CheckScalarDivisor(mapped_type divisor) const {
+    if (std::abs(divisor) < abs_tolerance) {
+      throw std::invalid_argument("Division by zero scalar not supported!");
+    }
+  }
+
   IndexContainer index2value_;
 };
 
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 067d687..3896297 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -871,6 +871,55 @@ TEST_CASE("multiply", "[mvPolynomial]") {
     REQUIRE_FALSE(cancel_prod.contains(Eigen::Array2i({1, 0})));
     REQUIRE(cancel_prod.get(Eigen::Array2i({2, 0})) == -1);
   }
+
+  SECTION("binary scalar multiplication works on both sides") {
+    auto p        = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+    auto expected = MP2({
+        {{0, 0}, 3},
+        {{1, 0}, 6},
+        {{0, 1}, 9},
+    });
+
+    REQUIRE(p * 3.0 == expected);
+    REQUIRE(3.0 * p == expected);
+  }
+
+  SECTION("scalar division works and preserves canonical zero") {
+    auto p        = MP2({
+        {{0, 0}, 2},
+        {{1, 0}, 4},
+        {{0, 1}, 6},
+    });
+    auto expected = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+
+    auto inplace = p;
+    inplace /= 2.0;
+
+    REQUIRE(inplace == expected);
+    REQUIRE(p / 2.0 == expected);
+
+    auto zero = MP2();
+    zero /= 2.0;
+    REQUIRE(zero == MP2());
+  }
+
+  SECTION("scalar division by zero throws") {
+    auto p = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+    });
+
+    REQUIRE_THROWS_AS(p / 0.0, std::invalid_argument);
+    REQUIRE_THROWS_AS(p /= 0.0, std::invalid_argument);
+  }
 }
 
 TEST_CASE("sum", "[mvPolynomial]") {
@@ -900,6 +949,17 @@ TEST_CASE("sum", "[mvPolynomial]") {
   for (size_t i = 0; i < ans.size(); ++i) {
     REQUIRE(sum.get(ans[i].first) == ans[i].second);
   }
+
+  SECTION("binary scalar addition works on both sides") {
+    auto expected = MP2({
+        {{0, 0}, 4},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+    });
+
+    REQUIRE(l + 3.0 == expected);
+    REQUIRE(3.0 + l == expected);
+  }
 }
 
 TEST_CASE("sub", "[mvPolynomial]") {
@@ -927,4 +987,20 @@ TEST_CASE("sub", "[mvPolynomial]") {
   for (size_t i = 0; i < ans.size(); ++i) {
     REQUIRE(sub.get(ans[i].first) == ans[i].second);
   }
+
+  SECTION("binary scalar subtraction works on both sides") {
+    auto poly_minus_scalar = MP2({
+        {{0, 0}, -2},
+        {{1, 0},  2},
+        {{0, 1},  3},
+    });
+    auto scalar_minus_poly = MP2({
+        {{0, 0},  2},
+        {{1, 0}, -2},
+        {{0, 1}, -3},
+    });
+
+    REQUIRE(l - 3.0 == poly_minus_scalar);
+    REQUIRE(3.0 - l == scalar_minus_poly);
+  }
 }

From 5aa910bd97609902a84a30beee0fb64cdd926dd2 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Sun, 17 May 2026 23:36:02 +0900
Subject: [PATCH 39/41] Update README.

---
 README.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/README.md b/README.md
index f9e9ed8..7abc681 100644
--- a/README.md
+++ b/README.md
@@ -78,9 +78,12 @@ auto q   = MP2({
 auto sum = p + q;
 auto sub = p - q;
 auto mul = p * q;
+auto div = p / 2.0;
 auto pw  = p.pow(3);
 ```
 
+`p /= scalar` and `p / scalar` are supported. Division by zero scalar throws `std::invalid_argument`.
+
 You can evaluate a polynomial at a point.
 
 ```cpp

From 70cb874d48b0eac1d42067b0d2a705c489d89d88 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Mon, 18 May 2026 01:03:12 +0900
Subject: [PATCH 40/41] Fix.

---
 include/mvPolynomial/mvPolynomial.hpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index 9502766..799f19f 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -15,6 +15,7 @@
 #include <numeric>
 #include <optional>
 #include <stdexcept>
+#include <string>
 #include <type_traits>
 #include <vector>
 
@@ -188,7 +189,7 @@ class MVPolynomial final {
 
   MVPolynomial pow(int exp) const {
     if (exp < 0) {
-      throw std::invalid_argument("Given exp must be positive.");
+      throw std::invalid_argument(std::string("Given exp must be positive."));
     }
     switch (exp) {
       case 0:
@@ -587,7 +588,7 @@ class MVPolynomial final {
 
   void CheckIndexIncludingNegative(const index_type& index) const {
     if ((index < 0).any()) {
-      throw std::invalid_argument("Negative index not supported!");
+      throw std::invalid_argument(std::string("Negative index not supported!"));
     }
   }
 
@@ -599,7 +600,7 @@ class MVPolynomial final {
 
   void CheckScalarDivisor(mapped_type divisor) const {
     if (std::abs(divisor) < abs_tolerance) {
-      throw std::invalid_argument("Division by zero scalar not supported!");
+      throw std::invalid_argument(std::string("Division by zero scalar not supported!"));
     }
   }
 

From 43c1ba2eabe6d4991db88837359f57b7b0802df2 Mon Sep 17 00:00:00 2001
From: Yuya Asano <64895419+sukeya@users.noreply.github.com>
Date: Mon, 18 May 2026 07:50:06 +0900
Subject: [PATCH 41/41] Set ASAN_OPTIONS=alloc_dealloc_mismatch=0.

---
 .github/workflows/test.yml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 65e22ca..90d8b76 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -33,7 +33,9 @@ jobs:
           cd ./build/debug/clang && ninja
       - name: Run tests
         run: |
-          ./build/debug/clang/test/mvPolynomial_test
+          # libc++ exception internals can trip ASan's alloc/dealloc mismatch check
+          # on Ubuntu runners even when the library code itself is fine.
+          ASAN_OPTIONS=alloc_dealloc_mismatch=0 ./build/debug/clang/test/mvPolynomial_test
 
   clang_memory:
     # libstdc++ and libc++ have uninitialized variables,