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: '^' - 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 ... 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/.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, 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++] diff --git a/CMakeLists.txt b/CMakeLists.txt index fe3e4be..45e274f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -63,7 +63,7 @@ set(gcc "$") set(clang "$") set(msvc "$") -set(gcc_like $) if (MVPOLYNOMIAL_BUILD_TESTS) enable_testing() diff --git a/README.md b/README.md index 1aef24c..7abc681 100644 --- a/README.md +++ b/README.md @@ -33,17 +33,81 @@ 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. +```cpp +#include "mvPolynomial/mvPolynomial.hpp" -## 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. +using MP2 = mvPolynomial::MVPolynomial; +``` + +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); +``` + +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. + +```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 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 +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 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); +``` -## 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. +Negative indices are not supported and cause `std::invalid_argument`. 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. diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp index c574bba..799f19f 100644 --- a/include/mvPolynomial/mvPolynomial.hpp +++ b/include/mvPolynomial/mvPolynomial.hpp @@ -6,10 +6,18 @@ #include #include +#include #include +#include #include +#include +#include +#include +#include #include +#include #include +#include #include "fmt/core.h" #include "platanus/btree_map.hpp" @@ -18,48 +26,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) - ); - } -} - -template -auto OfImpl(Iterator begin, Iterator end, int dim, int axis, const Coord& x) { - using Index = std::remove_cvref_t; - - CheckAxis(dim, axis); - - if (axis == dim - 1) { - auto last_coeff = begin->second; - auto last_index = Index{begin->first}; - for (auto it = std::next(begin); it != end; ++it) { - const auto& [next_index, next_coeff] = *it; - last_coeff *= std::pow(x[axis], last_index[axis] - next_index[axis]); - last_coeff += next_coeff; - last_index = next_index; - } - last_coeff *= x.pow(last_index.template cast()).prod(); - return last_coeff; - } else { - auto sum = typename Iterator::value_type::second_type(0); - while (true) { - const auto& [first_index, first_coeff] = *begin; - auto partition_point = std::partition_point( - begin, - end, - [axis, &first_index](const typename Iterator::value_type& pair) { - return pair.first[axis] == first_index[axis]; - } - ); - sum += OfImpl(begin, partition_point, dim, axis + 1, x); - if (partition_point == end) { - // The calculation ends. - break; - } - begin = partition_point; - } - return sum; + throw std::runtime_error(fmt::format("CheckAxis: Given axis {} must be in [0, {}).", axis, dim)); } } } // namespace details @@ -67,23 +34,23 @@ auto OfImpl(Iterator begin, Iterator end, int dim, int axis, const Coord& x) { template < std::signed_integral IntType, std::floating_point R, - int D, - class Allocator = std::allocator, R>>> + int dim_, + class Allocator = std::allocator, 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. - static R tolerance; + inline static R rel_tolerance = std::numeric_limits::epsilon(); + inline static R abs_tolerance = std::numeric_limits::min(); - using index_type = IndexType; + using index_type = IndexType; using coord_type = CoordType; private: - using Comparer = IndexComparer; - using IndexContainer = platanus::btree_map; + using IndexContainer = platanus::btree_map, Allocator>; public: using key_type = IndexContainer::key_type; @@ -110,69 +77,67 @@ class MVPolynomial final { using reverse_iterator = IndexContainer::reverse_iterator; using const_reverse_iterator = IndexContainer::const_reverse_iterator; + private: + template + using rebound_allocator_type = typename std::allocator_traits::template rebind_alloc; + + template + using rebound_vector_type = std::vector>; + + public: MVPolynomial(const MVPolynomial& other) = default; MVPolynomial& operator=(const MVPolynomial& other) = default; MVPolynomial(MVPolynomial&& other) = default; MVPolynomial& operator=(MVPolynomial&& other) = default; ~MVPolynomial() = default; - MVPolynomial() - : index2value_({ - {index_type::Zero(), 0} - }) {} + MVPolynomial() : index2value_(allocator_type{}) { AssignCoeffRaw(index_type::Zero(), R{0}); } - explicit MVPolynomial(const allocator_type& allocator) - : index2value_( - { - {index_type::Zero(), 0} - }, - allocator - ) {} + explicit MVPolynomial(const allocator_type& allocator) : index2value_(allocator) { + AssignCoeffRaw(index_type::Zero(), R{0}); + } template - MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e) { + MVPolynomial(InputIterator s, InputIterator e) : index2value_(s, e, allocator_type{}) { CheckSelfIndexes(); + Normalize(); } template - 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(); + Normalize(); } - MVPolynomial(std::initializer_list l, const allocator_type& a = allocator_type{}) - : index2value_(l, a) { + MVPolynomial(std::initializer_list l, const allocator_type& a = allocator_type{}) : index2value_(l, a) { CheckSelfIndexes(); + Normalize(); } 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 l) { index2value_.clear(); index2value_.insert(l.begin(), l.end()); CheckSelfIndexes(); + Normalize(); 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) { + AssignCoeffRaw(index_type::Zero(), r); + Normalize(); + } allocator_type get_allocator() const noexcept { return index2value_.get_allocator(); } + key_compare key_comp() const noexcept { return index2value_.key_comp(); } - Comparer 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(); } @@ -181,48 +146,25 @@ 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(); } - mapped_type& operator[](const key_type& index) { return index2value_[index]; } - mapped_type& operator[](key_type&& index) { return index2value_[index]; } + mapped_type get(const key_type& index) const { return index2value_.at(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); } - - void swap(MVPolynomial& m) { index2value_.swap(m.index2value_); } - - void clear() { index2value_.clear(); } - - // I didn't want to add insert and erase, but for efficiency, I did. - std::pair insert(const value_type& x) { return index2value_.insert(x); } - std::pair insert(value_type&& x) { return index2value_.insert(std::move(x)); } - - iterator insert(iterator hint, const value_type& x) { return index2value_.insert(hint, x); } - iterator insert(iterator hint, value_type&& x) { return index2value_.insert(hint, std::move(x)); } - - template - void insert(InputIterator b, InputIterator e) { - index2value_.insert(b, e); + std::optional try_get(const key_type& index) const { + if (auto it = index2value_.find(index); it != index2value_.end()) { + return it->second; + } + return std::nullopt; } - void insert(std::initializer_list list) { - index2value_.insert(list.begin(), list.end()); - } + void set(const key_type& index, mapped_type coeff) { AssignCoeff(index, coeff); } - 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); } + 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 - iterator find(const K& i) { - return index2value_.find(i); - } - template const_iterator find(const K& i) const { return index2value_.find(i); } @@ -233,31 +175,21 @@ 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 - iterator lower_bound(const K& i) { - return index2value_.lower_bound(i); - } - template 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 - iterator upper_bound(const K& i) { - return index2value_.upper_bound(i); - } - template const_iterator upper_bound(const K& i) const { return index2value_.upper_bound(i); } 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: @@ -269,10 +201,11 @@ class MVPolynomial final { default: auto max_pow2_under_exp = std::bit_floor(static_cast(exp)); auto max_bit_width = std::bit_width(max_pow2_under_exp); - auto cache = std::vector(max_bit_width - 1, get_allocator()); - cache.at(0) = (*this) * (*this); + auto cache = rebound_vector_type(rebound_allocator_type(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()); @@ -291,22 +224,122 @@ 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 partial_sums = rebound_vector_type(rebound_allocator_type(get_allocator())); + auto partial_sum = rbegin()->second; + 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())).prod(); + } else { + partial_sum *= (x.array().pow(index.template cast())).prod(); + partial_sums.push_back(partial_sum); + partial_sum = next_coeff; + } + } + { + const auto& first_index = begin()->first; + partial_sum *= (x.array().pow(first_index.template cast())).prod(); + } + return std::reduce(partial_sums.cbegin(), partial_sums.cend()) + partial_sum; + } + + MVPolynomial operator()(const MVPolynomial& mvp, int axis) const { + details::CheckAxis(dim, axis); + + auto partial_sums = rebound_vector_type(rebound_allocator_type(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; + 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_diff[axis])) * tmp_mvp; + } else { + 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)); + partial_sum = MVPolynomial{next_coeff, get_allocator()}; + } + } + { + const auto& first_index = begin()->first; + 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; + } + return std::reduce(partial_sums.cbegin(), partial_sums.cend()) + partial_sum; + } + + MVPolynomial D(int axis) const { + details::CheckAxis(dim, axis); - MVPolynomial operator()(const MVPolynomial& x, int axis) const { + 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 composed_mvp = MVPolynomial{get_allocator()}; - composed_mvp.clear(); - OfImpl(composed_mvp, rbegin(), rend(), 0, axis, x); - return composed_mvp; + 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& { 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; } @@ -314,7 +347,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; } @@ -322,65 +355,66 @@ class MVPolynomial final { } MVPolynomial& operator+=(mapped_type r) { - auto idx = index_type::Zero(); - if (contains(idx)) { - (*this)[idx] += r; - } else { - (*this)[idx] = r; - } + AddCoeff(index_type::Zero(), r); return *this; } MVPolynomial& operator+=(const MVPolynomial& r) { - for (const auto& [idx, coeff] : r) { - if (contains(idx)) { - (*this)[idx] += coeff; - } else { - (*this)[idx] = coeff; - } + auto copied_r = (&r == this) ? std::optional(MVPolynomial(r, get_allocator())) : std::nullopt; + const auto& rhs = copied_r ? *copied_r : r; + for (const auto& [idx, coeff] : rhs) { + AddCoeffRaw(idx, coeff); } + Normalize(); return *this; } MVPolynomial& operator-=(mapped_type r) { - auto idx = index_type::Zero(); - if (contains(idx)) { - (*this)[idx] -= r; - } else { - (*this)[idx] = r; - } + AddCoeff(index_type::Zero(), -r); return *this; } MVPolynomial& operator-=(const MVPolynomial& r) { - for (const auto& [idx, coeff] : r) { - if (contains(idx)) { - (*this)[idx] -= coeff; - } else { - (*this)[idx] = -coeff; - } + auto copied_r = (&r == this) ? std::optional(MVPolynomial(r, get_allocator())) : std::nullopt; + const auto& rhs = copied_r ? *copied_r : r; + for (const auto& [idx, coeff] : rhs) { + AddCoeffRaw(idx, -coeff); } + Normalize(); return *this; } MVPolynomial& operator*=(mapped_type r) { - auto idx = index_type::Zero(); - for (auto& index_and_coeff : *this) { + for (auto& index_and_coeff : index2value_) { auto& coeff = index_and_coeff.second; - coeff *= 2; + coeff *= r; } + Normalize(); + 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()); - for (auto& index_and_coeff : *this) { - auto& index = const_cast(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; + result.AddCoeffRaw(new_index, new_coeff); } + result.Normalize(); + swap(result); } else { *this = *this * r; } @@ -402,8 +436,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; @@ -422,10 +463,18 @@ class MVPolynomial final { return std::move(l); } - friend MVPolynomial operator+(const MVPolynomial& l, MVPolynomial&& r) { - return std::move(r) + 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; } friend MVPolynomial operator-(const MVPolynomial& l, const MVPolynomial& r) { @@ -437,17 +486,24 @@ class MVPolynomial final { return std::move(l); } - friend MVPolynomial operator-(const MVPolynomial& l, MVPolynomial&& r) { - return -std::move(r) + 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; } friend MVPolynomial operator*(const MVPolynomial& l, const MVPolynomial& r) { - 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; @@ -455,163 +511,110 @@ 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[idx] += v; - } else { - mul[idx] = v; - } + mul.AddCoeffRaw(idx, v); } } + mul.Normalize(); + 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 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 Normalize() { + auto normalized = IndexContainer(get_allocator()); + for (const auto& [index, coeff] : index2value_) { + if (std::abs(coeff) >= abs_tolerance) { + normalized[index] = coeff; + } + } + if (normalized.empty()) { + normalized[index_type::Zero()] = R{0}; } + index2value_.swap(normalized); } - 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; - - 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]); + 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}); + } + } - // The calculation ends. - if (partition_point == end) { - break; - } - begin = partition_point; - mvp.clear(); - } + void AssignCoeff(const key_type& index, mapped_type coeff) { + AssignCoeffRaw(index, coeff); + NormalizeSingle(index); + } + + void AddCoeff(const key_type& index, mapped_type delta) { + AddCoeffRaw(index, delta); + NormalizeSingle(index); + } + + void AssignCoeffRaw(const key_type& index, mapped_type coeff) { + CheckIndexIncludingNegative(index); + if (auto it = index2value_.find(index); it != index2value_.end()) { + it->second = coeff; } 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; - } + index2value_[index] = coeff; } } - IndexContainer index2value_; -}; - -template -R MVPolynomial::tolerance = - std::ldexp(std::numeric_limits::epsilon(), std::numeric_limits::min_exponent); - -template -auto D(const MVPolynomial& p, int axis) { - using MP = MVPolynomial; - 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; + void AddCoeffRaw(const key_type& index, mapped_type delta) { + CheckIndexIncludingNegative(index); + if (auto it = index2value_.find(index); it != index2value_.end()) { + it->second += delta; } else { - value *= index[axis]--; - dp[index] = value; - ++p_it; + index2value_[index] = delta; } } - return dp; -} -template -auto Integrate(MVPolynomial&& p, int axis) { - using MP = MVPolynomial; + void CheckIndexIncludingNegative(const index_type& index) const { + if ((index < 0).any()) { + throw std::invalid_argument(std::string("Negative index not supported!")); + } + } - details::CheckAxis(D, axis); + void CheckSelfIndexes() const { + for (const auto& index_and_coeff : index2value_) { + CheckIndexIncludingNegative(index_and_coeff.first); + } + } - for (auto& index_and_value : p) { - auto& value = index_and_value.second; - auto& index = const_cast(index_and_value.first); - value /= ++index[axis]; + void CheckScalarDivisor(mapped_type divisor) const { + if (std::abs(divisor) < abs_tolerance) { + throw std::invalid_argument(std::string("Division by zero scalar not supported!")); + } } - return std::move(p); + + IndexContainer index2value_; +}; + +template +auto D(const MVPolynomial& p, int axis) { + return p.D(axis); } -template -auto Integrate(const MVPolynomial& p, int axis) { - return Integrate(MVPolynomial(p, p.get_allocator()), axis); +template +auto Integrate(const MVPolynomial& p, int axis) { + return p.Integrate(axis); } } // namespace mvPolynomial 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>" diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp index 4db4ad0..3896297 100644 --- a/test/mvPolynomial_test.cpp +++ b/test/mvPolynomial_test.cpp @@ -1,6 +1,12 @@ #include "mvPolynomial/mvPolynomial.hpp" +#include #include +#include +#include +#include +#include +#include #include #include @@ -8,6 +14,183 @@ using MP2 = mvPolynomial::MVPolynomial; using MP3 = mvPolynomial::MVPolynomial; +template +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::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 +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 +Poly MakeRandomPolynomial(URNG& rng, int term_count, int max_degree, int max_abs_coeff) { + auto terms = std::vector>{}; + + auto degree_dist = std::uniform_int_distribution(0, max_degree); + auto coeff_dist = std::uniform_int_distribution(-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; + } + + auto found = false; + for (auto& [existing_index, existing_coeff] : terms) { + if ((existing_index == index).all()) { + existing_coeff += static_cast(coeff); + found = true; + break; + } + } + + if (!found) { + terms.emplace_back(index, static_cast(coeff)); + } + } + + auto filtered_terms = std::vector>{}; + 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(); + 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(coeff)); + } + + return Poly(filtered_terms.begin(), filtered_terms.end()); +} + +template +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(min_coord, max_coord); + for (int axis = 0; axis < Poly::dim; ++axis) { + point(axis) = static_cast(coord_dist(rng)); + } + return point; +} + +template +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 +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 +class CountingAllocator { + public: + using value_type = T; + + CountingAllocator() : allocations_(std::make_shared(0)) {} + + explicit CountingAllocator(std::shared_ptr allocations) : allocations_(std::move(allocations)) {} + + template + CountingAllocator(const CountingAllocator& other) : allocations_(other.allocations_) {} + + [[nodiscard]] T* allocate(std::size_t n) { + *allocations_ += n; + return std::allocator{}.allocate(n); + } + + void deallocate(T* p, std::size_t n) noexcept { std::allocator{}.deallocate(p, n); } + + template + bool operator==(const CountingAllocator& other) const noexcept { + return allocations_ == other.allocations_; + } + + template + bool operator!=(const CountingAllocator& other) const noexcept { + return !(*this == other); + } + + std::shared_ptr allocations_; +}; + +using CountingPairAllocator = CountingAllocator>; +using CountingMP2 = mvPolynomial::MVPolynomial; + TEST_CASE("constructor", "[mvPolynomial]") { SECTION("default") { auto ans = std::vector>({ @@ -16,8 +199,8 @@ TEST_CASE("constructor", "[mvPolynomial]") { auto m = MP2(); REQUIRE(m.size() == ans.size()); - for (auto i = 0; i < ans.size(); ++i) { - REQUIRE(m[ans[i].first] == ans[i].second); + for (size_t i = 0; i < ans.size(); ++i) { + REQUIRE(m.get(ans[i].first) == ans[i].second); } } @@ -35,10 +218,21 @@ 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); } } + SECTION("empty range is canonical zero polynomial") { + auto empty = std::vector>{}; + 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>({ {{0, 0}, 1}, @@ -60,9 +254,164 @@ 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); } } + + 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({ + {{-1, 0}, 1}, + }), + std::invalid_argument + ); + } +} + +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("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}, + {{1, 0}, 2}, + {{0, 1}, 3}, + }); + + p -= p; + + require_canonical_zero(p); + } + + SECTION("multiplying by zero stays canonical zero polynomial") { + auto p = MP2({ + {{0, 0}, 1}, + {{1, 0}, 2}, + {{0, 1}, 3}, + }); + + p *= 0.0; + + require_canonical_zero(p); + } + + SECTION("set rejects negative index for lvalue and rvalue keys") { + auto p = MP2(); + auto bad_key = Eigen::Array2i({-1, 0}); + + 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); + } + + 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]") { @@ -87,6 +436,75 @@ 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(1, 12); + + for (int trial = 0; trial < 200; ++trial) { + auto p = MakeRandomPolynomial(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(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({ @@ -102,6 +520,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}, @@ -115,11 +546,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 +558,30 @@ 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}, + }) + ); + } + + 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}, }) ); } @@ -150,14 +599,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,17 +614,56 @@ 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}, }) ); } + + SECTION("point_randomized_matches_naive_sum") { + auto rng = std::mt19937(20260517); + auto term_count_dist = std::uniform_int_distribution(1, 12); + + for (int trial = 0; trial < 200; ++trial) { + auto m = MakeRandomPolynomial(rng, term_count_dist(rng), 4, 5); + auto x = MakeRandomPoint(rng, -3, 3); + auto actual = m(x); + auto expected = NaivePointEval(m, x); + + CAPTURE(trial, DescribePolynomial(m), DescribePoint(x)); + REQUIRE(actual == expected); + } + } + + SECTION("composition_randomized_matches_naive_sum") { + auto rng = std::mt19937(20260518); + auto outer_term_count_dist = std::uniform_int_distribution(1, 10); + auto inner_term_count_dist = std::uniform_int_distribution(1, 6); + auto axis_dist = std::uniform_int_distribution(0, MP3::dim - 1); + + for (int trial = 0; trial < 120; ++trial) { + auto outer = MakeRandomPolynomial(rng, outer_term_count_dist(rng), 3, 4); + auto inner = MakeRandomPolynomial(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]") { @@ -199,7 +686,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); } } @@ -213,9 +700,19 @@ 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); } } + + 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]") { @@ -241,7 +738,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); } } @@ -258,9 +755,47 @@ 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); } } + + 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.get(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.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]") { @@ -292,7 +827,98 @@ 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") { + 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()); + } + + 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); + } + + 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); } } @@ -321,7 +947,18 @@ 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); + } + + 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); } } @@ -347,7 +984,23 @@ TEST_CASE("sub", "[mvPolynomial]") { auto sub = l - r; REQUIRE(sub.size() == ans.size()); - for (auto i = 0; i < ans.size(); ++i) { - REQUIRE(sub[ans[i].first] == ans[i].second); + 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); } }