From 98bb9384ee6199eb647347981b7990ef722a87f8 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 16:48:31 +0200 Subject: [PATCH 01/21] Move ExprPtr related things to expr_ptr.cpp --- CMakeLists.txt | 2 + SeQuant/core/expressions/expr.cpp | 90 ---------- SeQuant/core/expressions/expr_operators.hpp | 87 --------- SeQuant/core/expressions/expr_ptr.cpp | 189 ++++++++++++++++++++ SeQuant/core/expressions/expr_ptr.hpp | 9 + 5 files changed, 200 insertions(+), 177 deletions(-) create mode 100644 SeQuant/core/expressions/expr_ptr.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b04b376fc9..ce2b1195b4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -311,6 +311,8 @@ set(SeQuant_symb_src SeQuant/core/expressions/expr_algorithms.cpp SeQuant/core/expressions/expr_algorithms.hpp SeQuant/core/expressions/expr_operators.hpp + SeQuant/core/expressions/expr_ptr.cpp + SeQuant/core/expressions/expr_ptr.hpp SeQuant/core/expressions/expr_range.hpp SeQuant/core/expressions/result_expr.cpp SeQuant/core/expressions/result_expr.hpp diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index a8ce23862e..ded445b7cc 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -77,90 +77,6 @@ ExprPtr &Expr::back() { return at(size() - 1); } const ExprPtr &Expr::back() const { return at(size() - 1); } -ExprPtr ExprPtr::clone() const & { - if (!*this) return {}; - return ExprPtr(as_shared_ptr()->clone()); -} - -ExprPtr ExprPtr::clone() && noexcept { return std::move(*this); } - -ExprPtr::base_type &ExprPtr::as_shared_ptr() & { - return static_cast(*this); -} -const ExprPtr::base_type &ExprPtr::as_shared_ptr() const & { - return static_cast(*this); -} -ExprPtr::base_type &&ExprPtr::as_shared_ptr() && { - return static_cast(*this); -} - -Expr &ExprPtr::operator*() & { - SEQUANT_ASSERT(this->operator bool()); - return *(this->get()); -} - -const Expr &ExprPtr::operator*() const & { - SEQUANT_ASSERT(this->operator bool()); - return *(this->get()); -} - -Expr &&ExprPtr::operator*() && { - SEQUANT_ASSERT(this->operator bool()); - return std::move(*(this->get())); -} - -ExprPtr &ExprPtr::operator+=(const ExprPtr &other) { - if (!other) return *this; - - if (!*this) { - *this = other.clone(); - } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator+=(*other); - } else if (as_shared_ptr()->is() && other->is()) { - *this = ex(this->as().value() + - other->as().value()); - } else { - *this = ex(ExprPtrList{*this, other}); - } - return *this; -} - -ExprPtr &ExprPtr::operator-=(const ExprPtr &other) { - if (!other) return *this; - - if (!*this) { - *this = ex(-1) * other.clone(); - } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator-=(*other); - } else if (as_shared_ptr()->is() && other->is()) { - *this = ex(this->as().value() - - other->as().value()); - } else { - *this = ex(ExprPtrList{*this, ex(-1, ExprPtrList{other})}); - } - return *this; -} - -ExprPtr &ExprPtr::operator*=(const ExprPtr &other) { - if (!other) return *this; - - if (!*this) { - *this = other.clone(); - } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator*=(*other); - } else if (as_shared_ptr()->is() && other->is()) { - *this = ex(this->as().value() * - other->as().value()); - } else { - *this = ex(ExprPtrList{*this, other}); - } - return *this; -} - -std::size_t ExprPtr::size() const { return this->get()->size(); } - -std::wstring ExprPtr::to_latex() const { return as_shared_ptr()->to_latex(); } - Exception Expr::not_implemented(const char *fn) const { std::ostringstream oss; oss << "Expr::" << fn @@ -183,12 +99,6 @@ Expr &Expr::operator+=(const Expr &) { throw not_implemented("operator+="); } Expr &Expr::operator-=(const Expr &) { throw not_implemented("operator-="); } -ExprPtr adjoint(const ExprPtr &expr) { - auto result = expr->clone(); - result->adjoint(); - return result; -} - void Constant::adjoint() { value_ = conj(value_); reset_hash_value(); diff --git a/SeQuant/core/expressions/expr_operators.hpp b/SeQuant/core/expressions/expr_operators.hpp index e1e0967e3c..2f83df592a 100644 --- a/SeQuant/core/expressions/expr_operators.hpp +++ b/SeQuant/core/expressions/expr_operators.hpp @@ -14,96 +14,9 @@ #include #include -#include namespace sequant { -inline bool operator==(const ExprPtr &left, const ExprPtr &right) { - return *left == *right; -} - -inline ExprPtr operator*(const ExprPtr &left, const ExprPtr &right) { - if (left.is() && right.is()) { - auto c_ = left->clone(); - auto &c = c_.as(); - c *= right.as(); - return c_; - } - - auto left_is_product = left->is(); - auto right_is_product = right->is(); - if (!left_is_product && !right_is_product) { - return ex(ExprPtrList{left, right}); - } else if (left_is_product) { - auto result = std::static_pointer_cast(left->clone()); - result->append(1, right); - return result; - } else { // right_is_product - auto result = std::static_pointer_cast(right->clone()); - result->prepend(1, left); - return result; - } - - SEQUANT_UNREACHABLE; -} - -/// Unlike @code operator*(const ExprPtr&, const ExprPtr&) @endcode this -/// produces a non-commutative product (i.e. NCProduct) -inline ExprPtr operator^(const ExprPtr &left, const ExprPtr &right) { - auto left_is_product = left->is(); - auto right_is_product = right->is(); - if (!left_is_product && !right_is_product) { - return ex(ExprPtrList{left, right}); - } else if (left_is_product) { - auto result = std::make_shared(left->clone().as()); - result->append(1, right); - return result; - } else { // right_is_product - auto result = std::make_shared(right->clone().as()); - result->prepend(1, left); - return result; - } - - SEQUANT_UNREACHABLE; -} - -inline ExprPtr operator+(const ExprPtr &left, const ExprPtr &right) { - auto left_is_sum = left->is(); - auto right_is_sum = right->is(); - if (!left_is_sum && !right_is_sum) { - return ex(ExprPtrList{left, right}); - } else if (left_is_sum) { - auto result = std::static_pointer_cast(left->clone()); - result->append(right); - return result; - } else { // right_is_sum - auto result = std::static_pointer_cast(right->clone()); - result->prepend(left); - return result; - } - - SEQUANT_UNREACHABLE; -} - -inline ExprPtr operator-(const ExprPtr &left, const ExprPtr &right) { - auto left_is_sum = left->is(); - if (!left_is_sum) { - return ex(ExprPtrList{ - left, - (right->is() ? ex(-right->as().value()) - : ex(-1, ExprPtrList{right}))}); - } else if (left_is_sum) { - auto result = std::static_pointer_cast(left->clone()); - if (right->is()) - result->append(ex(-right->as().value())); - else - result->append(ex(-1, ExprPtrList{right})); - return result; - } - - SEQUANT_UNREACHABLE; -} - template requires(std::constructible_from) ExprPtr operator+(const ExprPtr &lhs, T &&rhs) { diff --git a/SeQuant/core/expressions/expr_ptr.cpp b/SeQuant/core/expressions/expr_ptr.cpp new file mode 100644 index 0000000000..d6f57fc710 --- /dev/null +++ b/SeQuant/core/expressions/expr_ptr.cpp @@ -0,0 +1,189 @@ +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace sequant { + +ExprPtr ExprPtr::clone() const & { + if (!*this) return {}; + return ExprPtr(as_shared_ptr()->clone()); +} + +ExprPtr ExprPtr::clone() && noexcept { return std::move(*this); } + +ExprPtr::base_type &ExprPtr::as_shared_ptr() & { + return static_cast(*this); +} +const ExprPtr::base_type &ExprPtr::as_shared_ptr() const & { + return static_cast(*this); +} +ExprPtr::base_type &&ExprPtr::as_shared_ptr() && { + return static_cast(*this); +} + +Expr &ExprPtr::operator*() & { + SEQUANT_ASSERT(this->operator bool()); + return *(this->get()); +} + +const Expr &ExprPtr::operator*() const & { + SEQUANT_ASSERT(this->operator bool()); + return *(this->get()); +} + +Expr &&ExprPtr::operator*() && { + SEQUANT_ASSERT(this->operator bool()); + return std::move(*(this->get())); +} + +ExprPtr &ExprPtr::operator+=(const ExprPtr &other) { + if (!other) return *this; + + if (!*this) { + *this = other.clone(); + } else if (as_shared_ptr()->is()) { + as_shared_ptr()->operator+=(*other); + } else if (as_shared_ptr()->is() && other->is()) { + *this = ex(this->as().value() + + other->as().value()); + } else { + *this = ex(ExprPtrList{*this, other}); + } + return *this; +} + +ExprPtr &ExprPtr::operator-=(const ExprPtr &other) { + if (!other) return *this; + + if (!*this) { + *this = ex(-1) * other.clone(); + } else if (as_shared_ptr()->is()) { + as_shared_ptr()->operator-=(*other); + } else if (as_shared_ptr()->is() && other->is()) { + *this = ex(this->as().value() - + other->as().value()); + } else { + *this = ex(ExprPtrList{*this, ex(-1, ExprPtrList{other})}); + } + return *this; +} + +ExprPtr &ExprPtr::operator*=(const ExprPtr &other) { + if (!other) return *this; + + if (!*this) { + *this = other.clone(); + } else if (as_shared_ptr()->is()) { + as_shared_ptr()->operator*=(*other); + } else if (as_shared_ptr()->is() && other->is()) { + *this = ex(this->as().value() * + other->as().value()); + } else { + *this = ex(ExprPtrList{*this, other}); + } + return *this; +} + +std::size_t ExprPtr::size() const { return this->get()->size(); } + +std::wstring ExprPtr::to_latex() const { return as_shared_ptr()->to_latex(); } + +ExprPtr adjoint(const ExprPtr &expr) { + auto result = expr->clone(); + result->adjoint(); + return result; +} + +bool operator==(const ExprPtr &left, const ExprPtr &right) { + return *left == *right; +} + +ExprPtr operator*(const ExprPtr &left, const ExprPtr &right) { + if (left.is() && right.is()) { + auto c_ = left->clone(); + auto &c = c_.as(); + c *= right.as(); + return c_; + } + + auto left_is_product = left->is(); + auto right_is_product = right->is(); + if (!left_is_product && !right_is_product) { + return ex(ExprPtrList{left, right}); + } else if (left_is_product) { + auto result = std::static_pointer_cast(left->clone()); + result->append(1, right); + return result; + } else { // right_is_product + auto result = std::static_pointer_cast(right->clone()); + result->prepend(1, left); + return result; + } + + SEQUANT_UNREACHABLE; +} + +/// Unlike @code operator*(const ExprPtr&, const ExprPtr&) @endcode this +/// produces a non-commutative product (i.e. NCProduct) +ExprPtr operator^(const ExprPtr &left, const ExprPtr &right) { + auto left_is_product = left->is(); + auto right_is_product = right->is(); + if (!left_is_product && !right_is_product) { + return ex(ExprPtrList{left, right}); + } else if (left_is_product) { + auto result = std::make_shared(left->clone().as()); + result->append(1, right); + return result; + } else { // right_is_product + auto result = std::make_shared(right->clone().as()); + result->prepend(1, left); + return result; + } + + SEQUANT_UNREACHABLE; +} + +ExprPtr operator+(const ExprPtr &left, const ExprPtr &right) { + auto left_is_sum = left->is(); + auto right_is_sum = right->is(); + if (!left_is_sum && !right_is_sum) { + return ex(ExprPtrList{left, right}); + } else if (left_is_sum) { + auto result = std::static_pointer_cast(left->clone()); + result->append(right); + return result; + } else { // right_is_sum + auto result = std::static_pointer_cast(right->clone()); + result->prepend(left); + return result; + } + + SEQUANT_UNREACHABLE; +} + +ExprPtr operator-(const ExprPtr &left, const ExprPtr &right) { + auto left_is_sum = left->is(); + if (!left_is_sum) { + return ex(ExprPtrList{ + left, + (right->is() ? ex(-right->as().value()) + : ex(-1, ExprPtrList{right}))}); + } else if (left_is_sum) { + auto result = std::static_pointer_cast(left->clone()); + if (right->is()) + result->append(ex(-right->as().value())); + else + result->append(ex(-1, ExprPtrList{right})); + return result; + } + + SEQUANT_UNREACHABLE; +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/expr_ptr.hpp b/SeQuant/core/expressions/expr_ptr.hpp index 5eb532d0d0..39b9ebbe40 100644 --- a/SeQuant/core/expressions/expr_ptr.hpp +++ b/SeQuant/core/expressions/expr_ptr.hpp @@ -157,6 +157,15 @@ using ExprPtrVector = container::svector; /// @return the adjoint of @p expr ExprPtr adjoint(const ExprPtr &expr); +ExprPtr operator*(const ExprPtr &left, const ExprPtr &right); + +/// Unlike @code operator*(const ExprPtr&, const ExprPtr&) @endcode this +/// produces a non-commutative product (i.e. NCProduct) +ExprPtr operator^(const ExprPtr &left, const ExprPtr &right); + +ExprPtr operator+(const ExprPtr &left, const ExprPtr &right); +ExprPtr operator-(const ExprPtr &left, const ExprPtr &right); + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_EXPR_PTR_HPP From 788c8ba2b4868496a814e9d999c99c455f7ff709 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 16:58:19 +0200 Subject: [PATCH 02/21] Also separate impl of Constant to dedicated file --- CMakeLists.txt | 2 + SeQuant/core/expressions/constant.cpp | 68 +++++++++++++++++++++++++++ SeQuant/core/expressions/constant.hpp | 59 ++++++----------------- SeQuant/core/expressions/expr.cpp | 5 -- 4 files changed, 84 insertions(+), 50 deletions(-) create mode 100644 SeQuant/core/expressions/constant.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ce2b1195b4..c762bedfae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -306,6 +306,8 @@ set(SeQuant_symb_src SeQuant/core/context.hpp SeQuant/core/expressions/abstract_tensor.cpp SeQuant/core/expressions/abstract_tensor.hpp + SeQuant/core/expressions/constant.cpp + SeQuant/core/expressions/constant.hpp SeQuant/core/expressions/expr.cpp SeQuant/core/expressions/expr.hpp SeQuant/core/expressions/expr_algorithms.cpp diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp new file mode 100644 index 0000000000..65d7d3b6a8 --- /dev/null +++ b/SeQuant/core/expressions/constant.cpp @@ -0,0 +1,68 @@ +#include +#include +#include +#include +#include + +namespace sequant { + +std::wstring Constant::to_latex() const { + return L"{" + io::latex::to_string(value()) + L"}"; +} + +Expr::type_id_type Constant::type_id() const { return get_type_id(); } + +bool Constant::is_scalar() const { return true; } + +ExprPtr Constant::clone() const { return ex(this->value()); } + +void Constant::adjoint() { + value_ = conj(value_); + reset_hash_value(); +} + +Expr &Constant::operator*=(const Expr &that) { + if (that.is()) { + value_ *= that.as().value(); + } else { + throw Exception("Constant::operator*=(that): not valid for that"); + } + return *this; +} + +Expr &Constant::operator+=(const Expr &that) { + if (that.is()) { + value_ += that.as().value(); + } else { + throw Exception("Constant::operator+=(that): not valid for that"); + } + return *this; +} + +Expr &Constant::operator-=(const Expr &that) { + if (that.is()) { + value_ -= that.as().value(); + } else { + throw Exception("Constant::operator-=(that): not valid for that"); + } + return *this; +} + +bool Constant::is_zero(scalar_type v) { return v.is_zero(); } + +bool Constant::is_zero() const { return is_zero(this->value()); } + +Expr::hash_type Constant::memoizing_hash() const { + if (!hash_value_) { + hash_value_ = hash::value(value_); + } else { + SEQUANT_ASSERT(*hash_value_ == hash::value(value_)); + } + return *hash_value_; +} + +bool Constant::static_equal(const Expr &that) const { + return value() == static_cast(that).value(); +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/constant.hpp b/SeQuant/core/expressions/constant.hpp index 9ae67a9cb2..88cc055ebc 100644 --- a/SeQuant/core/expressions/constant.hpp +++ b/SeQuant/core/expressions/constant.hpp @@ -3,8 +3,6 @@ #include #include -#include -#include #include #include @@ -14,6 +12,8 @@ namespace sequant { +class ExprPtr; + // implementation details of Constant; prefer sequant::detail over an unnamed // namespace in a header (see CppCoreGuidelines SF.21) namespace detail { @@ -67,68 +67,37 @@ class Constant : public Expr { throw Exception("Constant::value: cannot convert value to type T"); } - std::wstring to_latex() const override { - return L"{" + io::latex::to_string(value()) + L"}"; - } + std::wstring to_latex() const override; - type_id_type type_id() const override { return get_type_id(); } + type_id_type type_id() const override; - bool is_scalar() const override { return true; } + bool is_scalar() const override; - ExprPtr clone() const override { return ex(this->value()); } + ExprPtr clone() const override; /// @brief adjoint of a Constant is its complex conjugate virtual void adjoint() override; - virtual Expr &operator*=(const Expr &that) override { - if (that.is()) { - value_ *= that.as().value(); - } else { - throw Exception("Constant::operator*=(that): not valid for that"); - } - return *this; - } + Expr &operator*=(const Expr &that) override; - virtual Expr &operator+=(const Expr &that) override { - if (that.is()) { - value_ += that.as().value(); - } else { - throw Exception("Constant::operator+=(that): not valid for that"); - } - return *this; - } + virtual Expr &operator+=(const Expr &that) override; - virtual Expr &operator-=(const Expr &that) override { - if (that.is()) { - value_ -= that.as().value(); - } else { - throw Exception("Constant::operator-=(that): not valid for that"); - } - return *this; - } + virtual Expr &operator-=(const Expr &that) override; /// @param[in] v a scalar /// @return true if this is zero - static bool is_zero(scalar_type v) { return v.is_zero(); } + static bool is_zero(scalar_type v); /// @return `Constant::is_zero(this->value())` - bool is_zero() const final { return is_zero(this->value()); } + bool is_zero() const final; private: scalar_type value_; - hash_type memoizing_hash() const override { - if (!hash_value_) { - hash_value_ = hash::value(value_); - } else { - SEQUANT_ASSERT(*hash_value_ == hash::value(value_)); - } - return *hash_value_; - } + hash_type memoizing_hash() const override; + + bool static_equal(const Expr &that) const override; - bool static_equal(const Expr &that) const override { - return value() == static_cast(that).value(); - } }; // class Constant } // namespace sequant diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index ded445b7cc..d84875eace 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -99,11 +99,6 @@ Expr &Expr::operator+=(const Expr &) { throw not_implemented("operator+="); } Expr &Expr::operator-=(const Expr &) { throw not_implemented("operator-="); } -void Constant::adjoint() { - value_ = conj(value_); - reset_hash_value(); -} - std::wstring_view Variable::label() const { return label_; } void Variable::set_label(std::wstring label) { From 516504d9651990e13acde4f7c67908a2e50adc26 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 17:40:30 +0200 Subject: [PATCH 03/21] Remove in-place arithmetic operators from Expr interface They aren't universally supported by all expression subclasses. Therefore, them being defined in the general Expr API doesn't make too much sense. It's better to be notified of a missing operator via a compiler error than via a runtime exception. --- SeQuant/core/expressions/constant.cpp | 6 ++--- SeQuant/core/expressions/constant.hpp | 6 ++--- SeQuant/core/expressions/expr.cpp | 8 ------- SeQuant/core/expressions/expr.hpp | 32 --------------------------- SeQuant/core/expressions/expr_ptr.cpp | 6 ++--- SeQuant/core/expressions/power.hpp | 2 +- SeQuant/core/expressions/product.hpp | 2 +- SeQuant/core/expressions/sum.hpp | 11 ++++----- SeQuant/core/tensor_network/v1.cpp | 2 +- SeQuant/core/wick.impl.hpp | 2 +- 10 files changed, 19 insertions(+), 58 deletions(-) diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp index 65d7d3b6a8..29ec1a8201 100644 --- a/SeQuant/core/expressions/constant.cpp +++ b/SeQuant/core/expressions/constant.cpp @@ -21,7 +21,7 @@ void Constant::adjoint() { reset_hash_value(); } -Expr &Constant::operator*=(const Expr &that) { +Constant &Constant::operator*=(const Expr &that) { if (that.is()) { value_ *= that.as().value(); } else { @@ -30,7 +30,7 @@ Expr &Constant::operator*=(const Expr &that) { return *this; } -Expr &Constant::operator+=(const Expr &that) { +Constant &Constant::operator+=(const Expr &that) { if (that.is()) { value_ += that.as().value(); } else { @@ -39,7 +39,7 @@ Expr &Constant::operator+=(const Expr &that) { return *this; } -Expr &Constant::operator-=(const Expr &that) { +Constant &Constant::operator-=(const Expr &that) { if (that.is()) { value_ -= that.as().value(); } else { diff --git a/SeQuant/core/expressions/constant.hpp b/SeQuant/core/expressions/constant.hpp index 88cc055ebc..9bade3b542 100644 --- a/SeQuant/core/expressions/constant.hpp +++ b/SeQuant/core/expressions/constant.hpp @@ -78,11 +78,11 @@ class Constant : public Expr { /// @brief adjoint of a Constant is its complex conjugate virtual void adjoint() override; - Expr &operator*=(const Expr &that) override; + Constant &operator*=(const Expr &that); - virtual Expr &operator+=(const Expr &that) override; + Constant &operator+=(const Expr &that); - virtual Expr &operator-=(const Expr &that) override; + Constant &operator-=(const Expr &that); /// @param[in] v a scalar /// @return true if this is zero diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index d84875eace..704fae12d9 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -91,14 +91,6 @@ ExprPtr Expr::clone() const { throw not_implemented("clone"); } void Expr::adjoint() { throw not_implemented("adjoint"); } -Expr &Expr::operator*=(const Expr &) { throw not_implemented("operator*="); } - -Expr &Expr::operator^=(const Expr &) { throw not_implemented("operator^="); } - -Expr &Expr::operator+=(const Expr &) { throw not_implemented("operator+="); } - -Expr &Expr::operator-=(const Expr &) { throw not_implemented("operator-="); } - std::wstring_view Variable::label() const { return label_; } void Variable::set_label(std::wstring label) { diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 7bee7e8bc3..7501278e1f 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -330,38 +330,6 @@ class Expr : public std::enable_shared_from_this { return boost::core::demangle(typeid(*this).name()); } - /** @name in-place arithmetic operators - * Virtual in-place arithmetic operators to be overridden in expressions for - * which these make sense. - */ - ///@{ - - /// @brief in-place multiply @c *this by @c that - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator*=(const Expr &that); - - /// @brief in-place non-commutatively-multiply @c *this by @c that - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator^=(const Expr &that); - - /// @brief in-place add @c that to @c *this - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator+=(const Expr &that); - - /// @brief in-place subtract @c that from @c *this - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator-=(const Expr &that); - - ///@} - ExprIterator begin(); ExprIterator end(); ConstExprIterator begin() const; diff --git a/SeQuant/core/expressions/expr_ptr.cpp b/SeQuant/core/expressions/expr_ptr.cpp index d6f57fc710..062ddbfaa3 100644 --- a/SeQuant/core/expressions/expr_ptr.cpp +++ b/SeQuant/core/expressions/expr_ptr.cpp @@ -48,7 +48,7 @@ ExprPtr &ExprPtr::operator+=(const ExprPtr &other) { if (!*this) { *this = other.clone(); } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator+=(*other); + as() += *other; } else if (as_shared_ptr()->is() && other->is()) { *this = ex(this->as().value() + other->as().value()); @@ -64,7 +64,7 @@ ExprPtr &ExprPtr::operator-=(const ExprPtr &other) { if (!*this) { *this = ex(-1) * other.clone(); } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator-=(*other); + as() -= *other; } else if (as_shared_ptr()->is() && other->is()) { *this = ex(this->as().value() - other->as().value()); @@ -80,7 +80,7 @@ ExprPtr &ExprPtr::operator*=(const ExprPtr &other) { if (!*this) { *this = other.clone(); } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator*=(*other); + as() *= *other; } else if (as_shared_ptr()->is() && other->is()) { *this = ex(this->as().value() * other->as().value()); diff --git a/SeQuant/core/expressions/power.hpp b/SeQuant/core/expressions/power.hpp index f5731dabf3..1e94f17dfe 100644 --- a/SeQuant/core/expressions/power.hpp +++ b/SeQuant/core/expressions/power.hpp @@ -193,7 +193,7 @@ class Power : public Expr { /// and the effective conjugation parities align. For a Constant base /// only the fully unconjugated case combines. /// @throw Exception if @p that is not combinable. - Expr& operator*=(const Expr& that) override { + Power& operator*=(const Expr& that) { // b^e1 *= b^e2 -> b^(e1+e2) if (that.is()) { const auto& other = that.as(); diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 8d9e1a629a..eafd8cb79a 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -362,7 +362,7 @@ class Product : public Expr { return result; } - virtual Expr &operator*=(const Expr &that) override { + virtual Product &operator*=(const Expr &that) { if (!that.is()) { this->append(1, const_cast(that).shared_from_this()); } else { diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index 4cfcda2c8b..9b12855444 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -105,7 +105,7 @@ class Sum : public Expr { auto summand_constant = summand.as_shared_ptr(); if (constant_summand_idx_) { // add up to the existing constant ... SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); - *summands_[*constant_summand_idx_] += *summand_constant; + summands_[*constant_summand_idx_].as() += *summand_constant; do_erase = true; } else { // or memorize the position of the constant constant_summand_idx_ = pos; @@ -134,7 +134,7 @@ class Sum : public Expr { if (constant_summand_idx_) { SEQUANT_ASSERT( summands_.at(*constant_summand_idx_)->is()); - *(summands_[*constant_summand_idx_]) += *summand; + summands_[*constant_summand_idx_].as() += *summand; } else { summands_.push_back(summand->clone()); constant_summand_idx_ = summands_.size() - 1; @@ -162,7 +162,8 @@ class Sum : public Expr { if (constant_summand_idx_) { // add up to the existing constant ... SEQUANT_ASSERT( summands_.at(*constant_summand_idx_)->is()); - *summands_[*constant_summand_idx_] += *summand_constant; + summands_[*constant_summand_idx_].as() += + *summand_constant; } else { // or include the nonzero constant and update // constant_summand_idx_ summands_.insert(summands_.begin(), summand->clone()); @@ -261,12 +262,12 @@ class Sum : public Expr { /// @brief adjoint of a Sum is a sum of adjoints of its factors virtual void adjoint() override; - virtual Expr &operator+=(const Expr &that) override { + Sum &operator+=(const Expr &that) { this->append(const_cast(that).shared_from_this()); return *this; } - virtual Expr &operator-=(const Expr &that) override { + Sum &operator-=(const Expr &that) { if (that.is()) this->append(ex(-that.as().value())); else diff --git a/SeQuant/core/tensor_network/v1.cpp b/SeQuant/core/tensor_network/v1.cpp index cee3abce17..f351250420 100644 --- a/SeQuant/core/tensor_network/v1.cpp +++ b/SeQuant/core/tensor_network/v1.cpp @@ -458,7 +458,7 @@ ExprPtr TensorNetworkV1::canonicalize( nondefault_canonizer_ptr ? nondefault_canonizer_ptr.get() : &default_tensor_canonizer; auto bp = tensor_canonizer->apply(*tensor); - if (bp) *canon_byproduct *= *bp; + if (bp) canon_byproduct.as() *= *bp; } } edges_.clear(); diff --git a/SeQuant/core/wick.impl.hpp b/SeQuant/core/wick.impl.hpp index 89042c6f1f..f57579ee13 100644 --- a/SeQuant/core/wick.impl.hpp +++ b/SeQuant/core/wick.impl.hpp @@ -831,7 +831,7 @@ ExprPtr WickTheorem::compute(const bool count_only, nopseq->push_back(factor->template as>()); } else { SEQUANT_ASSERT(factor->is_cnumber()); - *prefactor *= *factor; + prefactor.as() *= *factor; } } init_input(nopseq); From 7c0211d42bc2f57bd85df05634fc1677b845f347 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 17:46:12 +0200 Subject: [PATCH 04/21] Make certain Expr API functions pure virtual This avoids situations in which important functions are not implemented for a given expression type (as was the case with NormalOperatorSequence). Thus, this gives a compiler-enforced guarantee that these functions will not just remain at the (useless) base implementations that just throw. --- SeQuant/core/expressions/expr.cpp | 14 ++------------ SeQuant/core/expressions/expr.hpp | 22 ++++------------------ SeQuant/core/op.hpp | 2 ++ tests/unit/test_expr.cpp | 7 +++---- 4 files changed, 11 insertions(+), 34 deletions(-) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index 704fae12d9..730003eb19 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -77,20 +77,10 @@ ExprPtr &Expr::back() { return at(size() - 1); } const ExprPtr &Expr::back() const { return at(size() - 1); } -Exception Expr::not_implemented(const char *fn) const { - std::ostringstream oss; - oss << "Expr::" << fn - << " not implemented in this derived class (type_name=" << type_name() - << ")"; - return Exception(oss.str()); +std::wstring Expr::to_latex() const { + throw Exception("to_latex not implemented for " + type_id()); } -std::wstring Expr::to_latex() const { throw not_implemented("to_latex"); } - -ExprPtr Expr::clone() const { throw not_implemented("clone"); } - -void Expr::adjoint() { throw not_implemented("adjoint"); } - std::wstring_view Variable::label() const { return label_; } void Variable::set_label(std::wstring label) { diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 7501278e1f..3d23259853 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -77,7 +77,7 @@ class Expr : public std::enable_shared_from_this { /// @return a clone of this object, i.e. an object that is equal to @c this /// @note - must be overridden in the derived class. /// - the default implementation throws an exception - virtual ExprPtr clone() const; + virtual ExprPtr clone() const = 0; /// like Expr::shared_from_this, but returns ExprPtr /// @return a shared_ptr to this object wrapped into ExprPtr, if this object @@ -230,7 +230,7 @@ class Expr : public std::enable_shared_from_this { /// @brief changes this to its adjoint /// @note base implementation throws, must be reimplemented in the derived /// class - virtual void adjoint(); + virtual void adjoint() = 0; /// Computes and returns the hash value. If default @p hasher is used then the /// value will be memoized, otherwise @p hasher will be used to compute the @@ -250,14 +250,7 @@ class Expr : public std::enable_shared_from_this { /// @note this function must be overridden in the derived class /// @sa Expr::get_type_id /// @return the hash value for this Expr - virtual type_id_type type_id() const -#if __GNUG__ - { - abort(); - } -#else - = 0; -#endif + virtual type_id_type type_id() const = 0; friend inline bool operator==(const Expr &a, const Expr &b); @@ -409,14 +402,7 @@ class Expr : public std::enable_shared_from_this { /// @note @c that is guaranteed to be of same type as @c *this, hence can be /// statically cast /// @return true if @c that is equivalent to *this - virtual bool static_equal([[maybe_unused]] const Expr &that) const -#if __GNUG__ - { - abort(); - } -#else - = 0; -#endif + virtual bool static_equal(const Expr &that) const = 0; /// @param that an Expr object /// @note @c that is guaranteed to be of same type as @c *this, hence can be diff --git a/SeQuant/core/op.hpp b/SeQuant/core/op.hpp index 5225429122..b07281ec2b 100644 --- a/SeQuant/core/op.hpp +++ b/SeQuant/core/op.hpp @@ -1061,6 +1061,8 @@ class NormalOperatorSequence : public container::svector>, return Expr::get_type_id(); }; + ExprPtr clone() const override { return ex(*this); } + friend bool operator==(const NormalOperatorSequence &nopseq1, const NormalOperatorSequence &nopseq2) { return nopseq1.vacuum() == nopseq2.vacuum() && diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index b9112c7ee6..e9a1e0924d 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -38,6 +38,7 @@ struct Dummy : public sequant::Expr { std::wstring to_latex() const override { return L"{\\text{Dummy}}"; } type_id_type type_id() const override { return get_type_id(); }; sequant::ExprPtr clone() const override { return sequant::ex(); } + void adjoint() override {} bool static_equal(const sequant::Expr &) const override { return true; } }; @@ -69,6 +70,8 @@ struct VecExpr : public std::vector, public sequant::Expr { type_id_type type_id() const override { return get_type_id>(); }; + void adjoint() override {} + sequant::ConstExprIterator begin_subexpr() const override { if constexpr (sequant::Expr::is_shared_ptr_of_expr::value) { return sequant::ConstExprIterator{base_type::data()}; @@ -413,10 +416,6 @@ TEST_CASE("expr", "[elements]") { } SECTION("adjoint") { - { // not implemented by default - const auto e = std::make_shared(); - REQUIRE_THROWS_AS(e->adjoint(), Exception); - } { // implemented in Adjointable const auto e = std::make_shared(); REQUIRE_NOTHROW(e->adjoint()); From 0f7724a2d05a4a391f3e698b4d8bc8dd96ac0503 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 17:55:54 +0200 Subject: [PATCH 05/21] Fix error message assembly --- SeQuant/core/expressions/expr.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index 730003eb19..aae04a52c9 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -78,7 +78,7 @@ ExprPtr &Expr::back() { return at(size() - 1); } const ExprPtr &Expr::back() const { return at(size() - 1); } std::wstring Expr::to_latex() const { - throw Exception("to_latex not implemented for " + type_id()); + throw Exception("to_latex not implemented for " + type_name()); } std::wstring_view Variable::label() const { return label_; } From 926eca9e3ea7af0d4babfc29a36a30874057dfe5 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 18:27:45 +0200 Subject: [PATCH 06/21] Clean up Variable impl --- CMakeLists.txt | 2 + SeQuant/core/expressions/expr.cpp | 22 --------- SeQuant/core/expressions/variable.cpp | 64 +++++++++++++++++++++++++++ SeQuant/core/expressions/variable.hpp | 38 +++++----------- 4 files changed, 76 insertions(+), 50 deletions(-) create mode 100644 SeQuant/core/expressions/variable.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index c762bedfae..bcb90b7bd7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -320,6 +320,8 @@ set(SeQuant_symb_src SeQuant/core/expressions/result_expr.hpp SeQuant/core/expressions/tensor.cpp SeQuant/core/expressions/tensor.hpp + SeQuant/core/expressions/variable.cpp + SeQuant/core/expressions/variable.hpp SeQuant/core/hash.cpp SeQuant/core/hash.hpp SeQuant/core/hugenholtz.hpp diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index aae04a52c9..b23a76ba04 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -26,7 +26,6 @@ #include #include -#include #include namespace sequant { @@ -81,27 +80,6 @@ std::wstring Expr::to_latex() const { throw Exception("to_latex not implemented for " + type_name()); } -std::wstring_view Variable::label() const { return label_; } - -void Variable::set_label(std::wstring label) { - label_ = std::move(label); - reset_hash_value(); -} - -void Variable::conjugate() { conjugated_ = !conjugated_; } - -bool Variable::conjugated() const { return conjugated_; } - -std::wstring Variable::to_latex() const { - std::wstring result = L"{" + io::latex::utf_to_string(label_) + L"}"; - if (conjugated_) result = L"{" + result + L"^*" + L"}"; - return result; -} - -ExprPtr Variable::clone() const { return ex(*this); } - -void Variable::adjoint() { conjugate(); } - bool Product::is_commutative() const { bool result = true; const auto nfactors = size(); diff --git a/SeQuant/core/expressions/variable.cpp b/SeQuant/core/expressions/variable.cpp new file mode 100644 index 0000000000..1f57611722 --- /dev/null +++ b/SeQuant/core/expressions/variable.cpp @@ -0,0 +1,64 @@ +#include +#include +#include +#include +#include + +#include +#include + +namespace sequant { + +Variable::Variable(std::wstring label) + : label_(std::move(label)), conjugated_(false) {} + +Variable::Variable(const std::string &label) + : label_(sequant::toUtf16(label)), conjugated_(false) {} + +Expr::type_id_type Variable::type_id() const { return get_type_id(); } + +bool Variable::is_scalar() const { return true; } + +Expr::hash_type Variable::memoizing_hash() const { + auto compute_hash = [this]() { + auto val = hash::value(label_); + hash::combine(val, conjugated_); + return val; + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + + return *hash_value_; +} + +bool Variable::static_equal(const Expr &that) const { + return label_ == static_cast(that).label_ && + conjugated_ == static_cast(that).conjugated_; +} + +std::wstring_view Variable::label() const { return label_; } + +void Variable::set_label(std::wstring label) { + label_ = std::move(label); + reset_hash_value(); +} + +void Variable::conjugate() { conjugated_ = !conjugated_; } + +bool Variable::conjugated() const { return conjugated_; } + +std::wstring Variable::to_latex() const { + std::wstring result = L"{" + io::latex::utf_to_string(label_) + L"}"; + if (conjugated_) result = L"{" + result + L"^*" + L"}"; + return result; +} + +ExprPtr Variable::clone() const { return ex(*this); } + +void Variable::adjoint() { conjugate(); } + +} // namespace sequant diff --git a/SeQuant/core/expressions/variable.hpp b/SeQuant/core/expressions/variable.hpp index 650ce64a39..4943845ac4 100644 --- a/SeQuant/core/expressions/variable.hpp +++ b/SeQuant/core/expressions/variable.hpp @@ -2,17 +2,16 @@ #define SEQUANT_EXPRESSIONS_VARIABLE_HPP #include -#include #include -#include #include -#include #include #include namespace sequant { +class ExprPtr; + /// This is represented as a "run-time" complex rational number class Variable : public Expr, public MutatableLabeled { public: @@ -29,10 +28,9 @@ class Variable : public Expr, public MutatableLabeled { std::constructible_from) explicit Variable(U &&label) : label_(std::forward(label)) {} - Variable(std::wstring label) : label_(std::move(label)), conjugated_(false) {} + Variable(std::wstring label); - Variable(const std::string &label) - : label_(sequant::toUtf16(label)), conjugated_(false) {} + Variable(const std::string &label); /// @return variable label /// @warning conjugation does not change it @@ -48,9 +46,9 @@ class Variable : public Expr, public MutatableLabeled { std::wstring to_latex() const override; - type_id_type type_id() const override { return get_type_id(); } + type_id_type type_id() const override; - bool is_scalar() const override { return true; } + bool is_scalar() const override; ExprPtr clone() const override; @@ -61,26 +59,10 @@ class Variable : public Expr, public MutatableLabeled { std::wstring label_; bool conjugated_ = false; - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - auto val = hash::value(label_); - hash::combine(val, conjugated_); - return val; - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - - return *hash_value_; - } - - bool static_equal(const Expr &that) const override { - return label_ == static_cast(that).label_ && - conjugated_ == static_cast(that).conjugated_; - } + hash_type memoizing_hash() const override; + + bool static_equal(const Expr &that) const override; + }; // class Variable } // namespace sequant From 5279d12f1ceba1d4f0e2992c710df58444d1dc48 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 18:27:58 +0200 Subject: [PATCH 07/21] Clean up Power impl --- CMakeLists.txt | 2 + SeQuant/core/expressions/power.cpp | 190 +++++++++++++++++++++++++++++ SeQuant/core/expressions/power.hpp | 184 +++------------------------- 3 files changed, 207 insertions(+), 169 deletions(-) create mode 100644 SeQuant/core/expressions/power.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index bcb90b7bd7..ac9120c7d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -316,6 +316,8 @@ set(SeQuant_symb_src SeQuant/core/expressions/expr_ptr.cpp SeQuant/core/expressions/expr_ptr.hpp SeQuant/core/expressions/expr_range.hpp + SeQuant/core/expressions/power.cpp + SeQuant/core/expressions/power.hpp SeQuant/core/expressions/result_expr.cpp SeQuant/core/expressions/result_expr.hpp SeQuant/core/expressions/tensor.cpp diff --git a/SeQuant/core/expressions/power.cpp b/SeQuant/core/expressions/power.cpp new file mode 100644 index 0000000000..f3920dac24 --- /dev/null +++ b/SeQuant/core/expressions/power.cpp @@ -0,0 +1,190 @@ +#include +#include +#include +#include +#include +#include + +namespace sequant { + +Power::Power(ExprPtr base, exponent_type exponent) + : base_{}, exponent_{std::move(exponent)} { + SEQUANT_ASSERT(base); + SEQUANT_ASSERT(base->is() || base->is()); + // clone on construction so that external + // mutations of the input cannot invalidate our memoized hash + base_ = base->clone(); + // 0^n is defined only for n >= 0 (0^0 = 1 by convention) + SEQUANT_ASSERT(!base_->is() || !base_->as().is_zero() || + exponent_ >= 0); +} + +const ExprPtr& Power::base() const { return base_; } + +const Power::exponent_type& Power::exponent() const { return exponent_; } + +bool Power::conjugated() const { return conjugated_; } + +void Power::conjugate() { + conjugated_ = !conjugated_; + reset_hash_value(); +} + +bool Power::is_zero() const { + return exponent_ > 0 && base_->is() && + base_->as().is_zero(); +} + +void Power::flatten(ExprPtr& expr) { + if (!expr || !expr->is()) return; + const auto& pw = expr->as(); + + // b^1 = b and conjugate if needed + if (pw.exponent_ == 1) { + auto lifted = pw.base_->clone(); + if (pw.conjugated_) lifted->adjoint(); + expr = std::move(lifted); + return; + } + // b^0 = 1 for any base (the ctor rejects 0^(negative) + if (pw.exponent_ == 0) { + expr = ex(Constant::scalar_type{1}); + return; + } + if (!pw.base_->is()) return; + + using scalar_type = Constant::scalar_type; + const auto& base_val = pw.base_->as().value(); + + // 1^k = 1 for any rational k. + if (base_val == scalar_type{1}) { + expr = ex(scalar_type{1}); + return; + } + + // Both remaining fold cases share one shape — `rational base raised to + // an integer exponent` — so we normalize to that shape and run a single + // exp-by-squaring loop. Anything else is left untouched. + // + // Case A: integer exponent (any Constant base, real or complex). + // `base` is just `base_val`. + // Case B: half-integer exponent on a non-negative real rational base + // `p/q` with both `p` and `q` perfect squares. Then + // (p/q)^(m/2) = (sqrt(p)/sqrt(q))^m, + // so we replace `base` with `sqrt(p)/sqrt(q)` (still a rational) and + // keep `exp_int = m`. + + // initialize the base + scalar_type base{0}; + auto exp_nr = numerator(pw.exponent_); // numerator of exponent + + if (denominator(pw.exponent_) == 1) { + base = base_val; + } else if (denominator(pw.exponent_) == 2 && base_val.imag() == 0 && + base_val.real() >= 0) { + intmax_t p = numerator(base_val.real()); + intmax_t q = denominator(base_val.real()); // > 0 by Boost's convention, + // sign is with the numerator + + // check for perfect squares + intmax_t p_rem{0}, q_rem{0}; + intmax_t p_root = boost::multiprecision::sqrt(p, p_rem); + intmax_t q_root = boost::multiprecision::sqrt(q, q_rem); + // fold if p and q are perfect squares, else return + if (p_rem != 0 || q_rem != 0) return; + base = scalar_type{rational(p_root) / rational(q_root)}; + } else { + return; + } + + // Standard exp-by-squaring; for negative exponents we power the + // magnitude and invert at the end. + const bool negate = exp_nr < 0; + if (negate) exp_nr = -exp_nr; + scalar_type value{1}; + scalar_type b = base; + while (exp_nr > 0) { + if (exp_nr % 2 != 0) value *= b; + exp_nr /= 2; + if (exp_nr > 0) b *= b; + } + if (negate) value = scalar_type{1} / value; + + if (pw.conjugated_) value = conj(value); + expr = ex(std::move(value)); +} + +Expr::type_id_type Power::type_id() const { return get_type_id(); } + +bool Power::is_scalar() const { return true; } + +ExprPtr Power::clone() const { + auto cloned = ex(base_, exponent_); + if (conjugated_) cloned->as().conjugate(); + return cloned; +} + +void Power::adjoint() { conjugate(); } + +Power& Power::operator*=(const Expr& that) { + // b^e1 *= b^e2 -> b^(e1+e2) + if (that.is()) { + const auto& other = that.as(); + if (conjugated_ == other.conjugated_ && *base_ == *other.base_) { + exponent_ += other.exponent_; + reset_hash_value(); + return *this; + } + } + // (b^e)* *= b* -> (b^(e+1))* + else if (base_->is() && that.is()) { + // check effective conjugation of Variable in this and that, if valid + // operation iff they match + const auto& base_var = base_->as(); + const auto& that_var = that.as(); + if (base_var.label() == that_var.label() && + (base_var.conjugated() ^ conjugated_) == that_var.conjugated()) { + exponent_ += rational{1}; + reset_hash_value(); + return *this; + } + } + // C^e *= C -> C^(e+1) + else if (!conjugated_ && *base_ == that) { + exponent_ += rational{1}; + reset_hash_value(); + return *this; + } + throw Exception("Power::operator*=(that): not valid for that"); +} + +Expr::hash_type Power::memoizing_hash() const { + auto compute_hash = [this]() { + if (exponent_ == 1 && !conjugated_) return hash::value(*base_); + auto val = hash::value(*base_); + hash::combine(val, hash::value(exponent_)); + hash::combine(val, conjugated_); + return val; + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + return *hash_value_; +} + +bool Power::static_equal(const Expr& that) const { + const auto& other = static_cast(that); + return exponent_ == other.exponent_ && conjugated_ == other.conjugated_ && + *base_ == *other.base_; +} + +bool Power::static_less_than(const Expr& that) const { + const auto& other = static_cast(that); + if (*base_ != *other.base_) return *base_ < *other.base_; + if (exponent_ != other.exponent_) return exponent_ < other.exponent_; + return conjugated_ < other.conjugated_; +} +} // namespace sequant diff --git a/SeQuant/core/expressions/power.hpp b/SeQuant/core/expressions/power.hpp index 1e94f17dfe..29fb69e505 100644 --- a/SeQuant/core/expressions/power.hpp +++ b/SeQuant/core/expressions/power.hpp @@ -5,10 +5,7 @@ #include #include #include -#include -#include #include -#include namespace sequant { @@ -27,17 +24,7 @@ class Power : public Expr { /// @param[in] base the base expression; must be a Constant or Variable. /// @param[in] exponent rational exponent - Power(ExprPtr base, exponent_type exponent) - : base_{}, exponent_{std::move(exponent)} { - SEQUANT_ASSERT(base); - SEQUANT_ASSERT(base->is() || base->is()); - // clone on construction so that external - // mutations of the input cannot invalidate our memoized hash - base_ = base->clone(); - // 0^n is defined only for n >= 0 (0^0 = 1 by convention) - SEQUANT_ASSERT(!base_->is() || !base_->as().is_zero() || - exponent_ >= 0); - } + Power(ExprPtr base, exponent_type exponent); /// @overload constructs a `Variable` base from @p label template @@ -55,29 +42,23 @@ class Power : public Expr { : Power(ex(std::forward(value)), std::move(exponent)) {} /// @return the base expression - const ExprPtr& base() const { return base_; } + const ExprPtr& base() const; /// @return the rational exponent - const exponent_type& exponent() const { return exponent_; } + const exponent_type& exponent() const; /// @return whether this Power has been complex-conjugated via adjoint() /// @note Conjugation is tracked as a flag because, in general, /// `conj(base^exponent) != conj(base)^exponent` - bool conjugated() const { return conjugated_; } + bool conjugated() const; /// @brief toggles the conjugation flag - void conjugate() { - conjugated_ = !conjugated_; - reset_hash_value(); - } + void conjugate(); /// @return true if the base is zero and the exponent is positive /// @note Construction rejects all undefined 0^n cases; 0^0 is legal and /// treated as 1. - bool is_zero() const override { - return exponent_ > 0 && base_->is() && - base_->as().is_zero(); - } + bool is_zero() const override; /// @brief Attempts to flatten a Power, mutating @p expr in place. Folds /// when @p expr holds a Power and any of: @@ -93,97 +74,16 @@ class Power : public Expr { /// @note Only square-root exponents are folded (that is the only /// case needed in practice right now). Extending to general n-th roots only /// requires replacing the integer-square-root step with an integer n-th-root. - static void flatten(ExprPtr& expr) { - if (!expr || !expr->is()) return; - const auto& pw = expr->as(); + static void flatten(ExprPtr& expr); - // b^1 = b and conjugate if needed - if (pw.exponent_ == 1) { - auto lifted = pw.base_->clone(); - if (pw.conjugated_) lifted->adjoint(); - expr = std::move(lifted); - return; - } - // b^0 = 1 for any base (the ctor rejects 0^(negative) - if (pw.exponent_ == 0) { - expr = ex(Constant::scalar_type{1}); - return; - } - if (!pw.base_->is()) return; + type_id_type type_id() const override; - using scalar_type = Constant::scalar_type; - const auto& base_val = pw.base_->as().value(); + bool is_scalar() const override; - // 1^k = 1 for any rational k. - if (base_val == scalar_type{1}) { - expr = ex(scalar_type{1}); - return; - } - - // Both remaining fold cases share one shape — `rational base raised to - // an integer exponent` — so we normalize to that shape and run a single - // exp-by-squaring loop. Anything else is left untouched. - // - // Case A: integer exponent (any Constant base, real or complex). - // `base` is just `base_val`. - // Case B: half-integer exponent on a non-negative real rational base - // `p/q` with both `p` and `q` perfect squares. Then - // (p/q)^(m/2) = (sqrt(p)/sqrt(q))^m, - // so we replace `base` with `sqrt(p)/sqrt(q)` (still a rational) and - // keep `exp_int = m`. - - // initialize the base - scalar_type base{0}; - auto exp_nr = numerator(pw.exponent_); // numerator of exponent - - if (denominator(pw.exponent_) == 1) { - base = base_val; - } else if (denominator(pw.exponent_) == 2 && base_val.imag() == 0 && - base_val.real() >= 0) { - intmax_t p = numerator(base_val.real()); - intmax_t q = denominator(base_val.real()); // > 0 by Boost's convention, - // sign is with the numerator - - // check for perfect squares - intmax_t p_rem{0}, q_rem{0}; - intmax_t p_root = boost::multiprecision::sqrt(p, p_rem); - intmax_t q_root = boost::multiprecision::sqrt(q, q_rem); - // fold if p and q are perfect squares, else return - if (p_rem != 0 || q_rem != 0) return; - base = scalar_type{rational(p_root) / rational(q_root)}; - } else { - return; - } - - // Standard exp-by-squaring; for negative exponents we power the - // magnitude and invert at the end. - const bool negate = exp_nr < 0; - if (negate) exp_nr = -exp_nr; - scalar_type value{1}; - scalar_type b = base; - while (exp_nr > 0) { - if (exp_nr % 2 != 0) value *= b; - exp_nr /= 2; - if (exp_nr > 0) b *= b; - } - if (negate) value = scalar_type{1} / value; - - if (pw.conjugated_) value = conj(value); - expr = ex(std::move(value)); - } - - type_id_type type_id() const override { return get_type_id(); } - - bool is_scalar() const override { return true; } - - ExprPtr clone() const override { - auto cloned = ex(base_, exponent_); - if (conjugated_) cloned->as().conjugate(); - return cloned; - } + ExprPtr clone() const override; /// @brief adjoint of Power: flips the conjugation flag. - void adjoint() override { conjugate(); } + void adjoint() override; /// @brief Combines exponents when effective bases match: /// - `b^e1 *= b^e2` → `b^(e1+e2)` when this and @p that share the same @@ -193,37 +93,7 @@ class Power : public Expr { /// and the effective conjugation parities align. For a Constant base /// only the fully unconjugated case combines. /// @throw Exception if @p that is not combinable. - Power& operator*=(const Expr& that) { - // b^e1 *= b^e2 -> b^(e1+e2) - if (that.is()) { - const auto& other = that.as(); - if (conjugated_ == other.conjugated_ && *base_ == *other.base_) { - exponent_ += other.exponent_; - reset_hash_value(); - return *this; - } - } - // (b^e)* *= b* -> (b^(e+1))* - else if (base_->is() && that.is()) { - // check effective conjugation of Variable in this and that, if valid - // operation iff they match - const auto& base_var = base_->as(); - const auto& that_var = that.as(); - if (base_var.label() == that_var.label() && - (base_var.conjugated() ^ conjugated_) == that_var.conjugated()) { - exponent_ += rational{1}; - reset_hash_value(); - return *this; - } - } - // C^e *= C -> C^(e+1) - else if (!conjugated_ && *base_ == that) { - exponent_ += rational{1}; - reset_hash_value(); - return *this; - } - throw Exception("Power::operator*=(that): not valid for that"); - } + Power& operator*=(const Expr& that); private: ExprPtr base_; @@ -232,35 +102,11 @@ class Power : public Expr { /// @return hash of this Power /// @note when exponent is 1 and not conjugated the hash matches the base's - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - if (exponent_ == 1 && !conjugated_) return hash::value(*base_); - auto val = hash::value(*base_); - hash::combine(val, hash::value(exponent_)); - hash::combine(val, conjugated_); - return val; - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - return *hash_value_; - } + hash_type memoizing_hash() const override; - bool static_equal(const Expr& that) const override { - const auto& other = static_cast(that); - return exponent_ == other.exponent_ && conjugated_ == other.conjugated_ && - *base_ == *other.base_; - } + bool static_equal(const Expr& that) const override; - bool static_less_than(const Expr& that) const override { - const auto& other = static_cast(that); - if (*base_ != *other.base_) return *base_ < *other.base_; - if (exponent_ != other.exponent_) return exponent_ < other.exponent_; - return conjugated_ < other.conjugated_; - } + bool static_less_than(const Expr& that) const override; }; } // namespace sequant From 8e476fc4dc2753f309f4c3faf55ba0414a7762fc Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 18:30:26 +0200 Subject: [PATCH 08/21] Add declaration of specialization of NormalOperator::labels() --- SeQuant/core/op.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/SeQuant/core/op.hpp b/SeQuant/core/op.hpp index b07281ec2b..4cc47e993d 100644 --- a/SeQuant/core/op.hpp +++ b/SeQuant/core/op.hpp @@ -955,6 +955,13 @@ class NormalOperator : public Operator, } }; +template <> +const container::svector & +NormalOperator::labels(); +template <> +const container::svector & +NormalOperator::labels(); + static_assert( is_tensor>, "The NormalOperator class does not fulfill the " From 873e0f1423e4f33973220dcb5f9ff548d04ad45e Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 09:29:33 +0200 Subject: [PATCH 09/21] Remove outdated comments --- SeQuant/core/expressions/expr.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 3d23259853..36b952501e 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -75,8 +75,6 @@ class Expr : public std::enable_shared_from_this { virtual std::wstring to_latex() const; /// @return a clone of this object, i.e. an object that is equal to @c this - /// @note - must be overridden in the derived class. - /// - the default implementation throws an exception virtual ExprPtr clone() const = 0; /// like Expr::shared_from_this, but returns ExprPtr @@ -228,8 +226,6 @@ class Expr : public std::enable_shared_from_this { } /// @brief changes this to its adjoint - /// @note base implementation throws, must be reimplemented in the derived - /// class virtual void adjoint() = 0; /// Computes and returns the hash value. If default @p hasher is used then the @@ -247,7 +243,6 @@ class Expr : public std::enable_shared_from_this { } /// Computes and returns the derived type identifier - /// @note this function must be overridden in the derived class /// @sa Expr::get_type_id /// @return the hash value for this Expr virtual type_id_type type_id() const = 0; From 0bb04beac2caf9461e4b32a51b7b5d6bbab0d06e Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 09:32:13 +0200 Subject: [PATCH 10/21] Make Product::operator* non-virtual --- SeQuant/core/expressions/product.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index eafd8cb79a..427d9462f6 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -362,7 +362,7 @@ class Product : public Expr { return result; } - virtual Product &operator*=(const Expr &that) { + Product &operator*=(const Expr &that) { if (!that.is()) { this->append(1, const_cast(that).shared_from_this()); } else { From 06c0b93e659027e390d263db2c8fd26643a4e415 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 09:36:06 +0200 Subject: [PATCH 11/21] Reset memoized hash after modification --- SeQuant/core/expressions/constant.cpp | 9 +++++++++ SeQuant/core/expressions/variable.cpp | 5 ++++- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp index 29ec1a8201..04bb759c54 100644 --- a/SeQuant/core/expressions/constant.cpp +++ b/SeQuant/core/expressions/constant.cpp @@ -27,6 +27,9 @@ Constant &Constant::operator*=(const Expr &that) { } else { throw Exception("Constant::operator*=(that): not valid for that"); } + + reset_hash_value(); + return *this; } @@ -36,6 +39,9 @@ Constant &Constant::operator+=(const Expr &that) { } else { throw Exception("Constant::operator+=(that): not valid for that"); } + + reset_hash_value(); + return *this; } @@ -45,6 +51,9 @@ Constant &Constant::operator-=(const Expr &that) { } else { throw Exception("Constant::operator-=(that): not valid for that"); } + + reset_hash_value(); + return *this; } diff --git a/SeQuant/core/expressions/variable.cpp b/SeQuant/core/expressions/variable.cpp index 1f57611722..58b15d00e9 100644 --- a/SeQuant/core/expressions/variable.cpp +++ b/SeQuant/core/expressions/variable.cpp @@ -47,7 +47,10 @@ void Variable::set_label(std::wstring label) { reset_hash_value(); } -void Variable::conjugate() { conjugated_ = !conjugated_; } +void Variable::conjugate() { + conjugated_ = !conjugated_; + reset_hash_value(); +} bool Variable::conjugated() const { return conjugated_; } From e37010c8675c077eb0b168de7258337bc940151b Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 09:36:48 +0200 Subject: [PATCH 12/21] Remove declaration of non-existent function --- SeQuant/core/expressions/expr.hpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 36b952501e..e1f7ba1870 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -432,12 +432,6 @@ class Expr : public std::enable_shared_from_this { static type_id_type type_id = get_next_type_id(); return type_id; } - - private: - /// @input[in] fn the name of function that is missing in this class - /// @return an Exception object containing a message describing that @p - /// fn is missing from this type - Exception not_implemented(const char *fn) const; }; // class Expr static_assert(std::ranges::sized_range); From 6ac7b1a581885c2319fa10565f7d4b673d6a2db9 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 09:59:42 +0200 Subject: [PATCH 13/21] Separate Product impl --- CMakeLists.txt | 2 + SeQuant/core/expressions/expr.cpp | 205 -------------- SeQuant/core/expressions/product.cpp | 384 +++++++++++++++++++++++++++ SeQuant/core/expressions/product.hpp | 157 ++--------- 4 files changed, 412 insertions(+), 336 deletions(-) create mode 100644 SeQuant/core/expressions/product.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ac9120c7d5..a84a9fd92b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -318,6 +318,8 @@ set(SeQuant_symb_src SeQuant/core/expressions/expr_range.hpp SeQuant/core/expressions/power.cpp SeQuant/core/expressions/power.hpp + SeQuant/core/expressions/product.cpp + SeQuant/core/expressions/product.hpp SeQuant/core/expressions/result_expr.cpp SeQuant/core/expressions/result_expr.hpp SeQuant/core/expressions/tensor.cpp diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index b23a76ba04..c4c8ac1b1b 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -26,8 +26,6 @@ #include #include -#include - namespace sequant { ExprIterator Expr::begin() { return begin_subexpr(); } @@ -80,209 +78,6 @@ std::wstring Expr::to_latex() const { throw Exception("to_latex not implemented for " + type_name()); } -bool Product::is_commutative() const { - bool result = true; - const auto nfactors = size(); - for (size_t f = 0; f != nfactors; ++f) { - for (size_t s = 1; result && s != nfactors; ++s) { - result &= factors_[f]->commutes_with(*factors_[s]); - } - } - return result; -} - -ExprPtr Product::canonicalize_impl(CanonicalizeOptions opts) { - // recursively canonicalize non-tensor subfactors (tensors will be - // canonicalized as part of the TN built of all tensor factors of this) ... - ranges::for_each(factors_, [this, opts](auto &factor) { - if (factor.template is()) { - return; - } - auto bp = factor->canonicalize(opts); - if (bp) { - SEQUANT_ASSERT(bp->template is()); - this->scalar_ *= std::static_pointer_cast(bp)->value(); - } - }); - - if (Logger::instance().canonicalize) { - std::wcout << "Product canonicalization(" << to_wstring(opts.method) - << ") input: " << to_latex() << std::endl; - } - - // pull out all scalar factors to the front - auto is_scalar = [](const auto &factor) { return factor->is_scalar(); }; - auto scalars = - factors_ | ranges::views::filter(is_scalar) | ranges::to_vector; - // scalars commute, so we can reorder them freely - ranges::sort(scalars, [](const auto &first, const auto &second) { - return *first < *second; - }); - - factors_ = factors_ | ranges::views::filter([&is_scalar](const auto &factor) { - return !is_scalar(factor); - }) | - ranges::to; - - // if there are no factors, insert scalars back and return - if (factors_.empty()) { - factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); - return {}; - } - - auto contains_nontensors = ranges::any_of(factors_, [](const auto &factor) { - return std::dynamic_pointer_cast(factor) == nullptr; - }); - if (!contains_nontensors) { // tensor network canonization is a special case - // that's done in - // TensorNetwork - auto make_canonical_tn = [this, &opts](auto *tn_null_ptr) { - using TN = std::decay_t>; - ExprPtr canon_factor; - TN tn(this->factors_); - if constexpr (TN::version() == 3) { - canon_factor = tn.canonicalize( - TensorCanonicalizer::cardinal_tensor_labels(), opts); - } else { - using NamedIndexSet = tensor_network::NamedIndexSet; - std::shared_ptr named_indices = - !opts.named_indices - ? nullptr - : std::make_shared(opts.named_indices->begin(), - opts.named_indices->end()); - canon_factor = tn.canonicalize( - TensorCanonicalizer::cardinal_tensor_labels(), - opts.method == CanonicalizationMethod::Rapid, named_indices.get()); - } - return std::pair{std::move(tn), canon_factor}; - }; - using TN = TensorNetwork; - auto [tn, canon_factor] = make_canonical_tn(static_cast(nullptr)); - - const auto &tensors = tn.tensors(); - using std::size; - SEQUANT_ASSERT(size(tensors) == size(factors_)); - using std::begin; - using std::end; - std::transform(begin(tensors), end(tensors), begin(factors_), - [](const auto &tptr) { - auto exprptr = std::dynamic_pointer_cast(tptr); - SEQUANT_ASSERT(exprptr); - return exprptr; - }); - if (canon_factor) scalar_ *= canon_factor->template as().value(); - this->reset_hash_value(); - } else { // if contains non-tensors, do commutation-checking resort - - // comparer that respects cardinal tensor labels - auto &cardinal_tensor_labels = - TensorCanonicalizer::cardinal_tensor_labels(); - auto local_compare = [&cardinal_tensor_labels](const ExprPtr &first, - const ExprPtr &second) { - if (first->is() && second->is()) { - const auto first_label = first->as().label(); - const auto second_label = second->as().label(); - if (first_label == second_label) return *first < *second; - const auto first_is_cardinal_it = ranges::find_if( - cardinal_tensor_labels, - [&first_label](const std::wstring &l) { return l == first_label; }); - const auto first_is_cardinal = - first_is_cardinal_it != ranges::end(cardinal_tensor_labels); - const auto second_is_cardinal_it = ranges::find_if( - cardinal_tensor_labels, [&second_label](const std::wstring &l) { - return l == second_label; - }); - const auto second_is_cardinal = - second_is_cardinal_it != ranges::end(cardinal_tensor_labels); - if (first_is_cardinal && second_is_cardinal) - return first_is_cardinal_it < second_is_cardinal_it; - else if (first_is_cardinal && !second_is_cardinal) - return true; - else if (!first_is_cardinal && second_is_cardinal) - return false; - else { - SEQUANT_ASSERT(!first_is_cardinal && !second_is_cardinal); - return *first < *second; - } - } else - return *first < *second; - }; - - // ... then resort, respecting commutativity - using std::begin; - using std::end; - if (static_commutativity()) { - if (is_commutative()) { - std::stable_sort(begin(factors_), end(factors_), local_compare); - } - } else { - // must do bubble sort if not commuting to avoid swapping elements across - // a noncommuting element - bubble_sort( - begin(factors_), end(factors_), - [&local_compare](const ExprPtr &first, const ExprPtr &second) { - bool result = (first->commutes_with(*second)) - ? local_compare(first, second) - : false; - return result; - }); - } - } - // reinsert scalar factors at the front - factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); - - // TODO evaluate product of Tensors (turn this into Products of Products) - - if (Logger::instance().canonicalize) - std::wcout << "Product canonicalization(" << to_wstring(opts.method) - << ") result: " << to_latex() << std::endl; - - return {}; // side effects are absorbed into the scalar_ -} - -void Product::adjoint() { - SEQUANT_ASSERT(static_commutativity() == false); // assert no slicing - auto adj_scalar = conj(scalar()); - using namespace ranges; - auto adj_factors = - factors() | views::reverse | - views::transform([](auto &expr) { return ::sequant::adjoint(expr); }); - using std::swap; - *this = - Product(adj_scalar, ranges::begin(adj_factors), ranges::end(adj_factors)); -} - -ExprPtr Product::canonicalize(CanonicalizeOptions opt) { - return this->canonicalize_impl(opt); -} - -ExprPtr Product::rapid_canonicalize(CanonicalizeOptions opt) { - SEQUANT_ASSERT(opt.method == CanonicalizationMethod::Rapid); - return this->canonicalize_impl(opt); -} - -void CProduct::adjoint() { - auto adj_scalar = conj(scalar()); - using namespace ranges; - // no need to reverse for commutative product - auto adj_factors = factors() | views::transform([](auto &&expr) { - return ::sequant::adjoint(expr); - }); - *this = CProduct(adj_scalar, ranges::begin(adj_factors), - ranges::end(adj_factors)); -} - -void NCProduct::adjoint() { - auto adj_scalar = conj(scalar()); - using namespace ranges; - // no need to reverse for commutative product - auto adj_factors = - factors() | views::reverse | - views::transform([](auto &&expr) { return ::sequant::adjoint(expr); }); - *this = NCProduct(adj_scalar, ranges::begin(adj_factors), - ranges::end(adj_factors)); -} - void Sum::adjoint() { using namespace ranges; auto adj_summands = summands() | views::transform([](auto &&expr) { diff --git a/SeQuant/core/expressions/product.cpp b/SeQuant/core/expressions/product.cpp new file mode 100644 index 0000000000..238a7e969f --- /dev/null +++ b/SeQuant/core/expressions/product.cpp @@ -0,0 +1,384 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace sequant { + +Product::Product(ExprPtrList factors, Flatten flatten_tag) { + using std::begin; + using std::end; + for (auto it = begin(factors); it != end(factors); ++it) + append(1, *it, flatten_tag); +} + +Product &Product::append(ExprPtr factor, Flatten flatten_tag) { + return this->append(1, factor, flatten_tag); +} + +const Product::scalar_type &Product::scalar() const { return scalar_; } + +bool Product::is_zero() const { return Constant::is_zero(this->scalar()); } + +const Product::factors_type &Product::factors() const { return factors_; } +Product::factors_type &Product::factors() { return factors_; } + +const ExprPtr &Product::factor(size_t i) const { return factors_.at(i); } + +bool Product::empty() const { return factors_.empty(); } + +bool Product::is_commutative() const { + bool result = true; + const auto nfactors = size(); + for (size_t f = 0; f != nfactors; ++f) { + for (size_t s = 1; result && s != nfactors; ++s) { + result &= factors_[f]->commutes_with(*factors_[s]); + } + } + return result; +} + +ExprPtr Product::canonicalize_impl(CanonicalizeOptions opts) { + // recursively canonicalize non-tensor subfactors (tensors will be + // canonicalized as part of the TN built of all tensor factors of this) ... + ranges::for_each(factors_, [this, opts](auto &factor) { + if (factor.template is()) { + return; + } + auto bp = factor->canonicalize(opts); + if (bp) { + SEQUANT_ASSERT(bp->template is()); + this->scalar_ *= std::static_pointer_cast(bp)->value(); + } + }); + + if (Logger::instance().canonicalize) { + std::wcout << "Product canonicalization(" << to_wstring(opts.method) + << ") input: " << to_latex() << std::endl; + } + + // pull out all scalar factors to the front + auto is_scalar = [](const auto &factor) { return factor->is_scalar(); }; + auto scalars = + factors_ | ranges::views::filter(is_scalar) | ranges::to_vector; + // scalars commute, so we can reorder them freely + ranges::sort(scalars, [](const auto &first, const auto &second) { + return *first < *second; + }); + + factors_ = factors_ | ranges::views::filter([&is_scalar](const auto &factor) { + return !is_scalar(factor); + }) | + ranges::to; + + // if there are no factors, insert scalars back and return + if (factors_.empty()) { + factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); + return {}; + } + + auto contains_nontensors = ranges::any_of(factors_, [](const auto &factor) { + return std::dynamic_pointer_cast(factor) == nullptr; + }); + if (!contains_nontensors) { // tensor network canonization is a special case + // that's done in + // TensorNetwork + auto make_canonical_tn = [this, &opts](auto *tn_null_ptr) { + using TN = std::decay_t>; + ExprPtr canon_factor; + TN tn(this->factors_); + if constexpr (TN::version() == 3) { + canon_factor = tn.canonicalize( + TensorCanonicalizer::cardinal_tensor_labels(), opts); + } else { + using NamedIndexSet = tensor_network::NamedIndexSet; + std::shared_ptr named_indices = + !opts.named_indices + ? nullptr + : std::make_shared(opts.named_indices->begin(), + opts.named_indices->end()); + canon_factor = tn.canonicalize( + TensorCanonicalizer::cardinal_tensor_labels(), + opts.method == CanonicalizationMethod::Rapid, named_indices.get()); + } + return std::pair{std::move(tn), canon_factor}; + }; + using TN = TensorNetwork; + auto [tn, canon_factor] = make_canonical_tn(static_cast(nullptr)); + + const auto &tensors = tn.tensors(); + using std::size; + SEQUANT_ASSERT(size(tensors) == size(factors_)); + using std::begin; + using std::end; + std::transform(begin(tensors), end(tensors), begin(factors_), + [](const auto &tptr) { + auto exprptr = std::dynamic_pointer_cast(tptr); + SEQUANT_ASSERT(exprptr); + return exprptr; + }); + if (canon_factor) scalar_ *= canon_factor->template as().value(); + this->reset_hash_value(); + } else { // if contains non-tensors, do commutation-checking resort + + // comparer that respects cardinal tensor labels + auto &cardinal_tensor_labels = + TensorCanonicalizer::cardinal_tensor_labels(); + auto local_compare = [&cardinal_tensor_labels](const ExprPtr &first, + const ExprPtr &second) { + if (first->is() && second->is()) { + const auto first_label = first->as().label(); + const auto second_label = second->as().label(); + if (first_label == second_label) return *first < *second; + const auto first_is_cardinal_it = ranges::find_if( + cardinal_tensor_labels, + [&first_label](const std::wstring &l) { return l == first_label; }); + const auto first_is_cardinal = + first_is_cardinal_it != ranges::end(cardinal_tensor_labels); + const auto second_is_cardinal_it = ranges::find_if( + cardinal_tensor_labels, [&second_label](const std::wstring &l) { + return l == second_label; + }); + const auto second_is_cardinal = + second_is_cardinal_it != ranges::end(cardinal_tensor_labels); + if (first_is_cardinal && second_is_cardinal) + return first_is_cardinal_it < second_is_cardinal_it; + else if (first_is_cardinal && !second_is_cardinal) + return true; + else if (!first_is_cardinal && second_is_cardinal) + return false; + else { + SEQUANT_ASSERT(!first_is_cardinal && !second_is_cardinal); + return *first < *second; + } + } else + return *first < *second; + }; + + // ... then resort, respecting commutativity + using std::begin; + using std::end; + if (static_commutativity()) { + if (is_commutative()) { + std::stable_sort(begin(factors_), end(factors_), local_compare); + } + } else { + // must do bubble sort if not commuting to avoid swapping elements across + // a noncommuting element + bubble_sort( + begin(factors_), end(factors_), + [&local_compare](const ExprPtr &first, const ExprPtr &second) { + bool result = (first->commutes_with(*second)) + ? local_compare(first, second) + : false; + return result; + }); + } + } + // reinsert scalar factors at the front + factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); + + // TODO evaluate product of Tensors (turn this into Products of Products) + + if (Logger::instance().canonicalize) + std::wcout << "Product canonicalization(" << to_wstring(opts.method) + << ") result: " << to_latex() << std::endl; + + return {}; // side effects are absorbed into the scalar_ +} + +void Product::adjoint() { + SEQUANT_ASSERT(static_commutativity() == false); // assert no slicing + auto adj_scalar = conj(scalar()); + using namespace ranges; + auto adj_factors = + factors() | views::reverse | + views::transform([](auto &expr) { return ::sequant::adjoint(expr); }); + using std::swap; + *this = + Product(adj_scalar, ranges::begin(adj_factors), ranges::end(adj_factors)); +} + +ExprPtr Product::canonicalize(CanonicalizeOptions opt) { + return this->canonicalize_impl(opt); +} + +ExprPtr Product::rapid_canonicalize(CanonicalizeOptions opt) { + SEQUANT_ASSERT(opt.method == CanonicalizationMethod::Rapid); + return this->canonicalize_impl(opt); +} + +void CProduct::adjoint() { + auto adj_scalar = conj(scalar()); + using namespace ranges; + // no need to reverse for commutative product + auto adj_factors = factors() | views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = CProduct(adj_scalar, ranges::begin(adj_factors), + ranges::end(adj_factors)); +} + +void NCProduct::adjoint() { + auto adj_scalar = conj(scalar()); + using namespace ranges; + // no need to reverse for commutative product + auto adj_factors = + factors() | std::views::reverse | std::views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = NCProduct(adj_scalar, ranges::begin(adj_factors), + ranges::end(adj_factors)); +} + +bool Product::static_commutativity() const { return false; } + +std::wstring Product::to_latex() const { return to_latex(false); } + +std::wstring Product::to_latex(bool negate) const { + std::wstring result; + result = L"{"; + if (!scalar().is_zero()) { + const auto scal = negate ? -scalar() : scalar(); + if (!scal.is_identity()) { + // replace -1 prefactor by - + if (!(negate ? scalar() : -scalar()).is_identity()) { + result += io::latex::to_string(scal); + } else { + result += L"{-}"; + } + } + for (const auto &i : factors()) { + if (i->is()) + result += L"\\bigl(" + i->to_latex() + L"\\bigr)"; + else + result += i->to_latex(); + } + } + result += L"}"; + return result; +} + +Product::type_id_type Product::type_id() const { + return get_type_id(); +}; + +/// @return an identical clone of this Product (a deep copy allocated on the +/// heap) +/// @note this does not flatten the product +ExprPtr Product::clone() const { return ex(this->deep_copy()); } + +Product Product::deep_copy() const { + auto cloned_factors = + factors() | ranges::views::transform([](const ExprPtr &ptr) { + return ptr ? ptr->clone() : nullptr; + }); + Product result(this->scalar(), ExprPtrList{}); + ranges::for_each(cloned_factors, [&](const auto &cloned_factor) { + result.append(1, std::move(cloned_factor), Flatten::No); + }); + return result; +} + +Product &Product::operator*=(const Expr &that) { + if (!that.is()) { + this->append(1, const_cast(that).shared_from_this()); + } else { + scalar_ *= that.as().value(); + } + return *this; +} + +void Product::add_identical(const Product &other) { + SEQUANT_ASSERT(ranges::equal(this->factors(), other.factors())); + scalar_ += other.scalar_; +} + +void Product::add_identical(const std::shared_ptr &other) { + SEQUANT_ASSERT(ranges::equal(this->factors(), other->factors())); + scalar_ += other->scalar_; +} + +void Product::add_identical(const ExprPtr &other) { + if (other.is()) return this->add_identical(other.as()); + + // only makes sense if this has a single factor + SEQUANT_ASSERT(this->factors_.size() == 1 && this->factors_[0] == other); + scalar_ += 1; +} + +ExprIterator Product::begin_subexpr() { + if (!factors_.empty()) { + reset_hash_value(); + } + + return ExprIterator{factors_.data()}; +} + +ExprIterator Product::end_subexpr() { + return ExprIterator{factors_.data() + factors_.size()}; +} + +ConstExprIterator Product::begin_subexpr() const { + return ConstExprIterator{factors_.data()}; +} + +ConstExprIterator Product::end_subexpr() const { + return ConstExprIterator{factors_.data() + factors_.size()}; +} + +/// @return the hash of this object, by hashing only the factors, +/// not the scalar to make possible rapid finding of Products that only +/// differ by a factor +/// @note this ensures that hash of a Product involving a single factor is +/// identical to the hash of the factor itself. +Expr::hash_type Product::memoizing_hash() const { + auto compute_hash = [this]() { + if (factors().size() == 1) + return factors_[0]->hash_value(); + else { + auto deref_factors = + factors() | + ranges::views::transform( + [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); + auto value = + hash::range(ranges::begin(deref_factors), ranges::end(deref_factors)); + return value; + } + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + + return *hash_value_; +} + +bool Product::static_equal(const Expr &that) const { + const auto &that_cast = static_cast(that); + if (scalar() == that_cast.scalar() && + factors().size() == that_cast.factors().size()) { + if (this->empty()) return true; + // compare hash values first + if (this->hash_value() == + that.hash_value()) // hash values agree -> do full comparison + return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), + expr_ptr_comparer); + else + return false; + } else + return false; +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 427d9462f6..9c3ba1d784 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -35,6 +35,7 @@ class Product : public Expr { enum class Flatten { Once, Recursively, Yes = Recursively, No }; using scalar_type = Constant::scalar_type; + using factors_type = container::svector; Product() = default; virtual ~Product() = default; @@ -46,12 +47,7 @@ class Product : public Expr { /// construct a Product out of zero or more factors (multiplied by 1) /// @param factors the factors /// @param flatten_tag if Flatten::Yes, flatten the factors - Product(ExprPtrList factors, Flatten flatten_tag = Flatten::Yes) { - using std::begin; - using std::end; - for (auto it = begin(factors); it != end(factors); ++it) - append(1, *it, flatten_tag); - } + Product(ExprPtrList factors, Flatten flatten_tag = Flatten::Yes); /// construct a Product out of zero or more factors (multiplied by 1) /// @param rng a range of factors; if rng is a Product, it will be flattened @@ -198,9 +194,7 @@ class Product : public Expr { /// @param factor a factor by which to multiply the product /// @param flatten_tag specifies whether (and how) to flatten the argument(s) /// @return @c *this - Product &append(ExprPtr factor, Flatten flatten_tag = Flatten::Yes) { - return this->append(1, factor, flatten_tag); - } + Product &append(ExprPtr factor, Flatten flatten_tag = Flatten::Yes); /// (post-)multiplies the product by @c factor /// @param factor a factor by which to multiply the product @@ -268,13 +262,13 @@ class Product : public Expr { flatten_tag); } - const auto &scalar() const { return scalar_; } + const scalar_type &scalar() const; /// @return `Constant::is_zero(this->scalar())` - bool is_zero() const override { return Constant::is_zero(this->scalar()); } + bool is_zero() const override; - const auto &factors() const { return factors_; } - auto &factors() { return factors_; } + const factors_type &factors() const; + factors_type &factors(); /// @brief View view of factors that are scalars (anything for which /// Expr::is_scalar() returns true). @@ -294,10 +288,10 @@ class Product : public Expr { /// Factor accessor /// @param i factor index /// @return ith factor - const ExprPtr &factor(size_t i) const { return factors_.at(i); } + const ExprPtr &factor(size_t i) const; /// @return true if the number of factors is zero - bool empty() const { return factors_.empty(); } + bool empty() const; /// @brief checks commutativity recursively /// @return true if definitely commutative, false definitely not commutative @@ -312,135 +306,50 @@ class Product : public Expr { private: /// @return true if commutativity is decidable statically /// @sa CProduct::static_commutativity() and NCProduct::static_commutativity() - virtual bool static_commutativity() const { return false; } + virtual bool static_commutativity() const; public: - std::wstring to_latex() const override { return to_latex(false); } + std::wstring to_latex() const override; /// just like Expr::to_latex() , but can negate before conversion /// @param[in] negate if true, scalar will be before conversion - std::wstring to_latex(bool negate) const { - std::wstring result; - result = L"{"; - if (!scalar().is_zero()) { - const auto scal = negate ? -scalar() : scalar(); - if (!scal.is_identity()) { - // replace -1 prefactor by - - if (!(negate ? scalar() : -scalar()).is_identity()) { - result += io::latex::to_string(scal); - } else { - result += L"{-}"; - } - } - for (const auto &i : factors()) { - if (i->is()) - result += L"\\bigl(" + i->to_latex() + L"\\bigr)"; - else - result += i->to_latex(); - } - } - result += L"}"; - return result; - } + std::wstring to_latex(bool negate) const; - type_id_type type_id() const override { return get_type_id(); }; + type_id_type type_id() const override; /// @return an identical clone of this Product (a deep copy allocated on the /// heap) /// @note this does not flatten the product - ExprPtr clone() const override { return ex(this->deep_copy()); } - - Product deep_copy() const { - auto cloned_factors = - factors() | ranges::views::transform([](const ExprPtr &ptr) { - return ptr ? ptr->clone() : nullptr; - }); - Product result(this->scalar(), ExprPtrList{}); - ranges::for_each(cloned_factors, [&](const auto &cloned_factor) { - result.append(1, std::move(cloned_factor), Flatten::No); - }); - return result; - } + ExprPtr clone() const override; - Product &operator*=(const Expr &that) { - if (!that.is()) { - this->append(1, const_cast(that).shared_from_this()); - } else { - scalar_ *= that.as().value(); - } - return *this; - } + Product deep_copy() const; - void add_identical(const Product &other) { - SEQUANT_ASSERT(ranges::equal(this->factors(), other.factors())); - scalar_ += other.scalar_; - } + Product &operator*=(const Expr &that); - void add_identical(const std::shared_ptr &other) { - SEQUANT_ASSERT(ranges::equal(this->factors(), other->factors())); - scalar_ += other->scalar_; - } - - void add_identical(const ExprPtr &other) { - if (other.is()) return this->add_identical(other.as()); + void add_identical(const Product &other); - // only makes sense if this has a single factor - SEQUANT_ASSERT(this->factors_.size() == 1 && this->factors_[0] == other); - scalar_ += 1; - } + void add_identical(const std::shared_ptr &other); - ExprIterator begin_subexpr() override { - if (!factors_.empty()) { - reset_hash_value(); - } + void add_identical(const ExprPtr &other); - return ExprIterator{factors_.data()}; - } + ExprIterator begin_subexpr() override; - ExprIterator end_subexpr() override { - return ExprIterator{factors_.data() + factors_.size()}; - } + ExprIterator end_subexpr() override; - ConstExprIterator begin_subexpr() const override { - return ConstExprIterator{factors_.data()}; - } + ConstExprIterator begin_subexpr() const override; - ConstExprIterator end_subexpr() const override { - return ConstExprIterator{factors_.data() + factors_.size()}; - } + ConstExprIterator end_subexpr() const override; private: scalar_type scalar_ = {1, 0}; - container::svector factors_{}; + factors_type factors_{}; /// @return the hash of this object, by hashing only the factors, /// not the scalar to make possible rapid finding of Products that only /// differ by a factor /// @note this ensures that hash of a Product involving a single factor is /// identical to the hash of the factor itself. - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - if (factors().size() == 1) - return factors_[0]->hash_value(); - else { - auto deref_factors = - factors() | - ranges::views::transform( - [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); - auto value = hash::range(ranges::begin(deref_factors), - ranges::end(deref_factors)); - return value; - } - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - - return *hash_value_; - } + hash_type memoizing_hash() const override; ExprPtr canonicalize_impl(CanonicalizeOptions); @@ -454,21 +363,7 @@ class Product : public Expr { CanonicalizationMethod::Rapid)) override; private: - bool static_equal(const Expr &that) const override { - const auto &that_cast = static_cast(that); - if (scalar() == that_cast.scalar() && - factors().size() == that_cast.factors().size()) { - if (this->empty()) return true; - // compare hash values first - if (this->hash_value() == - that.hash_value()) // hash values agree -> do full comparison - return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), - expr_ptr_comparer); - else - return false; - } else - return false; - } + bool static_equal(const Expr &that) const override; }; // class Product class CProduct : public Product { From 0fb64e554c71d2924471fc70140a698ca7fb0dad Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 10:01:41 +0200 Subject: [PATCH 14/21] Cleanly separate public and private interface --- SeQuant/core/expressions/product.hpp | 29 +++++++++++++--------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 9c3ba1d784..403a0973a9 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -303,12 +303,6 @@ class Product : public Expr { /// factors, with complex-conjugated scalar virtual void adjoint() override; - private: - /// @return true if commutativity is decidable statically - /// @sa CProduct::static_commutativity() and NCProduct::static_commutativity() - virtual bool static_commutativity() const; - - public: std::wstring to_latex() const override; /// just like Expr::to_latex() , but can negate before conversion @@ -340,10 +334,23 @@ class Product : public Expr { ConstExprIterator end_subexpr() const override; + virtual ExprPtr canonicalize( + CanonicalizeOptions opt = + CanonicalizeOptions::default_options()) override; + + virtual ExprPtr rapid_canonicalize( + CanonicalizeOptions opts = + CanonicalizeOptions::default_options().copy_and_set( + CanonicalizationMethod::Rapid)) override; + private: scalar_type scalar_ = {1, 0}; factors_type factors_{}; + /// @return true if commutativity is decidable statically + /// @sa CProduct::static_commutativity() and NCProduct::static_commutativity() + virtual bool static_commutativity() const; + /// @return the hash of this object, by hashing only the factors, /// not the scalar to make possible rapid finding of Products that only /// differ by a factor @@ -353,16 +360,6 @@ class Product : public Expr { ExprPtr canonicalize_impl(CanonicalizeOptions); - public: - virtual ExprPtr canonicalize( - CanonicalizeOptions opt = - CanonicalizeOptions::default_options()) override; - virtual ExprPtr rapid_canonicalize( - CanonicalizeOptions opts = - CanonicalizeOptions::default_options().copy_and_set( - CanonicalizationMethod::Rapid)) override; - - private: bool static_equal(const Expr &that) const override; }; // class Product From 6e6b6251f80ea39aeef78abc763fed887d1093f8 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 10:05:18 +0200 Subject: [PATCH 15/21] Cleanly separate (N)CProduct impl --- SeQuant/core/expressions/product.cpp | 60 +++++++++++++++++----------- SeQuant/core/expressions/product.hpp | 16 ++++---- 2 files changed, 45 insertions(+), 31 deletions(-) diff --git a/SeQuant/core/expressions/product.cpp b/SeQuant/core/expressions/product.cpp index 238a7e969f..12911a53f5 100644 --- a/SeQuant/core/expressions/product.cpp +++ b/SeQuant/core/expressions/product.cpp @@ -217,29 +217,6 @@ ExprPtr Product::rapid_canonicalize(CanonicalizeOptions opt) { return this->canonicalize_impl(opt); } -void CProduct::adjoint() { - auto adj_scalar = conj(scalar()); - using namespace ranges; - // no need to reverse for commutative product - auto adj_factors = factors() | views::transform([](auto &&expr) { - return ::sequant::adjoint(expr); - }); - *this = CProduct(adj_scalar, ranges::begin(adj_factors), - ranges::end(adj_factors)); -} - -void NCProduct::adjoint() { - auto adj_scalar = conj(scalar()); - using namespace ranges; - // no need to reverse for commutative product - auto adj_factors = - factors() | std::views::reverse | std::views::transform([](auto &&expr) { - return ::sequant::adjoint(expr); - }); - *this = NCProduct(adj_scalar, ranges::begin(adj_factors), - ranges::end(adj_factors)); -} - bool Product::static_commutativity() const { return false; } std::wstring Product::to_latex() const { return to_latex(false); } @@ -381,4 +358,41 @@ bool Product::static_equal(const Expr &that) const { return false; } +CProduct::CProduct(const Product &other) : Product(other) {} +CProduct::CProduct(Product &&other) : Product(std::move(other)) {} + +bool CProduct::is_commutative() const { return true; } + +void CProduct::adjoint() { + auto adj_scalar = conj(scalar()); + using namespace ranges; + // no need to reverse for commutative product + auto adj_factors = factors() | views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = CProduct(adj_scalar, ranges::begin(adj_factors), + ranges::end(adj_factors)); +} + +bool CProduct::static_commutativity() const { return true; } + +NCProduct::NCProduct(const Product &other) : Product(other) {} +NCProduct::NCProduct(Product &&other) : Product(std::move(other)) {} + +bool NCProduct::is_commutative() const { return false; } + +void NCProduct::adjoint() { + auto adj_scalar = conj(scalar()); + using namespace ranges; + // no need to reverse for commutative product + auto adj_factors = + factors() | std::views::reverse | std::views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = NCProduct(adj_scalar, ranges::begin(adj_factors), + ranges::end(adj_factors)); +} + +bool NCProduct::static_commutativity() const { return true; } + } // namespace sequant diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 403a0973a9..a11278b34d 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -366,10 +366,10 @@ class Product : public Expr { class CProduct : public Product { public: using Product::Product; - CProduct(const Product &other) : Product(other) {} - CProduct(Product &&other) : Product(other) {} + CProduct(const Product &other); + CProduct(Product &&other); - bool is_commutative() const override { return true; } + bool is_commutative() const override; /// @brief adjoint of a CProduct is a product of adjoints of its factors, with /// complex-conjugated scalar @@ -377,23 +377,23 @@ class CProduct : public Product { virtual void adjoint() override; private: - bool static_commutativity() const override { return true; } + bool static_commutativity() const override; }; // class CProduct class NCProduct : public Product { public: using Product::Product; - NCProduct(const Product &other) : Product(other) {} - NCProduct(Product &&other) : Product(other) {} + NCProduct(const Product &other); + NCProduct(Product &&other); - bool is_commutative() const override { return false; } + bool is_commutative() const override; /// @brief adjoint of a NCProduct is a reserved product of adjoints of its /// factors, with complex-conjugated scalar virtual void adjoint() override; private: - bool static_commutativity() const override { return true; } + bool static_commutativity() const override; }; // class NCProduct } // namespace sequant From 6120935b8572b338485f892c1ca2a583e6270dde Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 10:07:48 +0200 Subject: [PATCH 16/21] Move some headers to impl file --- SeQuant/core/expressions/product.cpp | 2 ++ SeQuant/core/expressions/product.hpp | 5 ----- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/SeQuant/core/expressions/product.cpp b/SeQuant/core/expressions/product.cpp index 12911a53f5..48ac8a2684 100644 --- a/SeQuant/core/expressions/product.cpp +++ b/SeQuant/core/expressions/product.cpp @@ -1,8 +1,10 @@ #include #include #include +#include #include #include +#include #include #include #include diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index a11278b34d..1c12ed2728 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -4,17 +4,12 @@ #include #include #include -#include #include #include -#include #include #include -#include -#include #include -#include #include #include From 94b967e1e4763979781fcb3ec1b25af5bc65e862 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 10:10:08 +0200 Subject: [PATCH 17/21] Remove redundant swap function Since C++11 std::swap will use move semantics so this custom swap impl doesn't get us any benefit. --- SeQuant/core/expressions/expr.cpp | 3 ++- SeQuant/core/expressions/sum.hpp | 6 ------ 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index c4c8ac1b1b..7eaf896fb4 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -135,7 +135,8 @@ ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { // since need to sort in both cases auto new_sum = (pass == npasses - 1) ? acc.make_canonicalized_sum() : acc.make_sum(); - this->swap(*new_sum); + using std::swap; + swap(*this, *new_sum); if (Logger::instance().canonicalize) std::wcout << "Sum::canonicalize_impl (pass=" << pass diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index 9b12855444..18bb09c873 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -35,12 +35,6 @@ class Sum : public Expr { Sum &operator=(const Sum &) = default; Sum &operator=(Sum &&) = default; - void swap(Sum &other) { - Sum tmp = std::move(other); - other = std::move(*this); - *this = std::move(tmp); - } - /// construct a Sum out of zero or more summands /// @param summands an initializer list of summands Sum(ExprPtrList summands) { From 1dd073482324c2f6a88f365cc18e44d5fd01cca8 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 11:13:14 +0200 Subject: [PATCH 18/21] Add missing header --- SeQuant/core/expressions/expr_ptr.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SeQuant/core/expressions/expr_ptr.cpp b/SeQuant/core/expressions/expr_ptr.cpp index 062ddbfaa3..b1e21a0bb5 100644 --- a/SeQuant/core/expressions/expr_ptr.cpp +++ b/SeQuant/core/expressions/expr_ptr.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include From 2d902884c96f7f5ea1ba2d7bb175ee2cfda7944f Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 11:16:21 +0200 Subject: [PATCH 19/21] Separate Sum impl --- CMakeLists.txt | 2 + SeQuant/core/expressions/expr.cpp | 149 ------------ SeQuant/core/expressions/sum.cpp | 386 ++++++++++++++++++++++++++++++ SeQuant/core/expressions/sum.hpp | 243 +++---------------- 4 files changed, 416 insertions(+), 364 deletions(-) create mode 100644 SeQuant/core/expressions/sum.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a84a9fd92b..f119f10dbc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -322,6 +322,8 @@ set(SeQuant_symb_src SeQuant/core/expressions/product.hpp SeQuant/core/expressions/result_expr.cpp SeQuant/core/expressions/result_expr.hpp + SeQuant/core/expressions/sum.cpp + SeQuant/core/expressions/sum.hpp SeQuant/core/expressions/tensor.cpp SeQuant/core/expressions/tensor.hpp SeQuant/core/expressions/variable.cpp diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index 7eaf896fb4..e147806d7f 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -78,155 +78,6 @@ std::wstring Expr::to_latex() const { throw Exception("to_latex not implemented for " + type_name()); } -void Sum::adjoint() { - using namespace ranges; - auto adj_summands = summands() | views::transform([](auto &&expr) { - return ::sequant::adjoint(expr); - }); - *this = Sum(ranges::begin(adj_summands), ranges::end(adj_summands)); -} - -ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { - if (Logger::instance().canonicalize) - std::wcout << "Sum::canonicalize_impl: input = " - << to_latex_align(shared_from_this()) << std::endl; - - const auto npasses = multipass ? 2 : 1; - for (auto pass = 0; pass != npasses; ++pass) { - const auto rapid = (pass % 2 == 0); - - // canonicalizing TNs in a sum requires treating named indices as - // meaningful/distinct - auto opts_copy = opts; - opts_copy.ignore_named_index_labels = - CanonicalizeOptions::IgnoreNamedIndexLabel::No; - if (rapid) { - opts_copy.method = CanonicalizationMethod::Lexicographic; - } else - opts_copy.method = opts.method | CanonicalizationMethod::Topological; - - // recursively canonicalize summands ... - // using for_each and direct access to summands - sequant::for_each(summands_, [pass, &opts_copy, &rapid](ExprPtr &summand) { - ExprPtr bp; - if (rapid) { - bp = summand->rapid_canonicalize(opts_copy); - } else { - bp = summand->canonicalize(opts_copy); - } - if (bp) { - SEQUANT_ASSERT(bp->template is()); - summand = ex(std::static_pointer_cast(bp)->value(), - ExprPtrList{summand}); - } - }); - if (Logger::instance().canonicalize) - std::wcout << "Sum::canonicalize_impl (pass=" << pass - << "): after canonicalizing summands = " - << to_latex_align(shared_from_this()) << std::endl; - - HashingAccumulator acc; - for (auto &summand : summands_) { - acc.append(summand); - } - - // last pass? sort by hash then by Expr::operator< - // N.B. no point in differentiating between canonicalization methods here - // since need to sort in both cases - auto new_sum = - (pass == npasses - 1) ? acc.make_canonicalized_sum() : acc.make_sum(); - using std::swap; - swap(*this, *new_sum); - - if (Logger::instance().canonicalize) - std::wcout << "Sum::canonicalize_impl (pass=" << pass - << "): after reducing summands = " - << to_latex_align(shared_from_this()) << std::endl; - } - - return {}; // side effects are absorbed into summands -} - -HashingAccumulator &HashingAccumulator::append(ExprPtr summand, bool flatten) { - // flatten, if needed - if (flatten && summand.is()) { - for (auto &subsummand : summand.as().summands()) { - this->append(subsummand, flatten); - } - return *this; - } - - // process summand as a whole - auto it = summands_.find(summand); - if (it == summands_.end()) { - summands_.emplace(summand); - } else { // found existing term with the same hash - auto existing_summand = *it; - if (summand.template is()) { - if (existing_summand.is()) { - // both are products - add them - existing_summand.as().add_identical( - summand.template as()); - } else { - // convert existing term to product and add - auto product_copy = std::make_shared(summand->clone()); - product_copy->add_identical(existing_summand); - summands_.erase(it); - summands_.emplace(std::move(product_copy)); - } - } else { - if (existing_summand.is()) { - existing_summand.as().add_identical(summand); - } else { - // neither is a product - create new product - auto product_form = std::make_shared(); - product_form->append(2, summand.template as()); - summands_.erase(it); - summands_.emplace(std::move(product_form)); - } - } - } - - return *this; -} - -SumPtr HashingAccumulator::make_sum_impl(bool canonicalize) { - Sum::summands_type summands; - summands.reserve(summands_.size()); - for (auto summand : summands_) { - if (!summand->is_zero()) { - summands.push_back(summand); - } - } - - if (canonicalize) { - ranges::sort(summands, [](const auto &e1, const auto &e2) { - if (e1->hash_value() == e2->hash_value()) { - return e1 < e2; - } else { - return e1->hash_value() < e2->hash_value(); - } - }); - } - - return std::make_shared(std::move(summands), Sum::move_only_tag{}); -} - -SumPtr HashingAccumulator::make_sum() { return make_sum_impl(false); } - -SumPtr HashingAccumulator::make_canonicalized_sum() { - return make_sum_impl(true); -} - -ExprPtr HashingAccumulator::make_expr(bool canonicalize) { - if (summands_.size() == 0) { - return ex(0); - } else if (summands_.size() == 1) - return *(summands_.begin()); - else - return make_sum_impl(canonicalize); -} - bool proportional_to::operator()(const ExprPtr &expr1, const ExprPtr &expr2) const { if (expr1->type_id() != diff --git a/SeQuant/core/expressions/sum.cpp b/SeQuant/core/expressions/sum.cpp new file mode 100644 index 0000000000..880c47d2dc --- /dev/null +++ b/SeQuant/core/expressions/sum.cpp @@ -0,0 +1,386 @@ +#include +#include +#include +#include +#include + +namespace sequant { + +Sum::Sum(ExprPtrList summands) { + // use append to flatten out Sum summands + for (auto &&summand : summands) { + append(std::forward(summand)); + } +} + +Sum::Sum(summands_type &&summands, move_only_tag) + : summands_(std::move(summands)) { + std::size_t pos = 0; + for (auto it = summands_.begin(); it != summands_.end(); ++it) { + auto &summand = *it; + bool do_erase = false; + if (summand->is_zero()) { + do_erase = true; + } else if (summand->is()) { + auto summand_constant = summand.as_shared_ptr(); + if (constant_summand_idx_) { // add up to the existing constant ... + SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); + summands_[*constant_summand_idx_].as() += *summand_constant; + do_erase = true; + } else { // or memorize the position of the constant + constant_summand_idx_ = pos; + } + } + + // erase if needed + if (do_erase) { + summands_.erase(it); + it = summands_.begin(); + std::advance(it, pos); + } else + ++pos; + } +} + +Sum &Sum::append(ExprPtr summand) { + SEQUANT_ASSERT(summand); + if (!summand->is()) { + if (!summand->is_zero()) { // exclude zeros + if (summand->is()) { // add up constants + // immediately, if possible + auto summand_constant = summand.as_shared_ptr(); + if (constant_summand_idx_) { + SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); + summands_[*constant_summand_idx_].as() += *summand; + } else { + summands_.push_back(summand->clone()); + constant_summand_idx_ = summands_.size() - 1; + } + } else { + summands_.push_back(summand->clone()); + } + reset_hash_value(); + } + } else { // this recursively flattens Sum summands + for (auto &subsummand : *summand) this->append(subsummand); + } + return *this; +} + +Sum &Sum::prepend(ExprPtr summand) { + SEQUANT_ASSERT(summand); + if (!summand->is()) { + if (!summand->is_zero()) { + // exclude zeros + if (summand->is()) { + auto summand_constant = summand.as_shared_ptr(); + if (constant_summand_idx_) { // add up to the existing constant ... + SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); + summands_[*constant_summand_idx_].as() += *summand_constant; + } else { // or include the nonzero constant and update + // constant_summand_idx_ + summands_.insert(summands_.begin(), summand->clone()); + constant_summand_idx_ = 0; + } + } else { + summands_.insert(summands_.begin(), summand->clone()); + if (constant_summand_idx_) // if have a constant, update its position + ++*constant_summand_idx_; + } + reset_hash_value(); + } + } else { // this recursively flattens Sum summands + for (auto &subsummand : *summand) this->prepend(subsummand); + } + return *this; +} + +const Sum::summands_type &Sum::summands() const { return summands_; } + +const ExprPtr &Sum::summand(size_t i) const { return summands_.at(i); } + +ExprPtr Sum::take_n(size_t count) const { + const auto e = (count >= summands_.size() ? summands_.end() + : (summands_.begin() + count)); + return ex(summands_.begin(), e); +} + +ExprPtr Sum::take_n(size_t offset, size_t count) const { + const auto offset_plus_count = offset + count; + const auto b = (offset >= summands_.size() ? summands_.end() + : (summands_.begin() + offset)); + const auto e = (offset_plus_count >= summands_.size() + ? summands_.end() + : (summands_.begin() + offset_plus_count)); + return ex(b, e); +} + +bool Sum::empty() const { return summands_.empty(); } + +std::size_t Sum::size() const { return summands_.size(); } + +std::wstring Sum::to_latex() const { + std::wstring result; + result = L"{ \\bigl("; + std::size_t counter = 0; + for (const auto &i : summands()) { + const auto i_is_product = i->is(); + if (!i_is_product) { + result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); + } else { // i_is_product + const auto i_prod = i->as(); + const auto scalar = i_prod.scalar(); + if (scalar.real() < 0 || (scalar.real() == 0 && scalar.imag() < 0)) { + result += L" - " + i_prod.to_latex(true); + } else { + result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); + } + } + ++counter; + } + result += L"\\bigr) }"; + return result; +} + +Expr::type_id_type Sum::type_id() const { return Expr::get_type_id(); }; + +ExprPtr Sum::clone() const { + auto cloned_summands = + summands() | + ranges::views::transform([](const ExprPtr &ptr) { return ptr->clone(); }); + return ex(ranges::begin(cloned_summands), ranges::end(cloned_summands)); +} + +void Sum::adjoint() { + using namespace ranges; + auto adj_summands = summands() | views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = Sum(ranges::begin(adj_summands), ranges::end(adj_summands)); +} + +ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { + if (Logger::instance().canonicalize) + std::wcout << "Sum::canonicalize_impl: input = " + << to_latex_align(shared_from_this()) << std::endl; + + const auto npasses = multipass ? 2 : 1; + for (auto pass = 0; pass != npasses; ++pass) { + const auto rapid = (pass % 2 == 0); + + // canonicalizing TNs in a sum requires treating named indices as + // meaningful/distinct + auto opts_copy = opts; + opts_copy.ignore_named_index_labels = + CanonicalizeOptions::IgnoreNamedIndexLabel::No; + if (rapid) { + opts_copy.method = CanonicalizationMethod::Lexicographic; + } else + opts_copy.method = opts.method | CanonicalizationMethod::Topological; + + // recursively canonicalize summands ... + // using for_each and direct access to summands + sequant::for_each(summands_, [&opts_copy, &rapid](ExprPtr &summand) { + ExprPtr bp; + if (rapid) { + bp = summand->rapid_canonicalize(opts_copy); + } else { + bp = summand->canonicalize(opts_copy); + } + if (bp) { + SEQUANT_ASSERT(bp->template is()); + summand = ex(std::static_pointer_cast(bp)->value(), + ExprPtrList{summand}); + } + }); + if (Logger::instance().canonicalize) + std::wcout << "Sum::canonicalize_impl (pass=" << pass + << "): after canonicalizing summands = " + << to_latex_align(shared_from_this()) << std::endl; + + HashingAccumulator acc; + for (auto &summand : summands_) { + acc.append(summand); + } + + // last pass? sort by hash then by Expr::operator< + // N.B. no point in differentiating between canonicalization methods here + // since need to sort in both cases + auto new_sum = + (pass == npasses - 1) ? acc.make_canonicalized_sum() : acc.make_sum(); + using std::swap; + swap(*this, *new_sum); + + if (Logger::instance().canonicalize) + std::wcout << "Sum::canonicalize_impl (pass=" << pass + << "): after reducing summands = " + << to_latex_align(shared_from_this()) << std::endl; + } + + return {}; // side effects are absorbed into summands +} + +Sum &Sum::operator+=(const Expr &that) { + this->append(const_cast(that).shared_from_this()); + return *this; +} + +Sum &Sum::operator-=(const Expr &that) { + if (that.is()) + this->append(ex(-that.as().value())); + else + this->append(ex( + -1, ExprPtrList{const_cast(that).shared_from_this()})); + return *this; +} + +ExprIterator Sum::begin_subexpr() { + if (!summands_.empty()) { + reset_hash_value(); + } + + return ExprIterator{summands_.data()}; +} + +ExprIterator Sum::end_subexpr() { + return ExprIterator{summands_.data() + summands_.size()}; +} + +ConstExprIterator Sum::begin_subexpr() const { + return ConstExprIterator{summands_.data()}; +} + +ConstExprIterator Sum::end_subexpr() const { + return ConstExprIterator{summands_.data() + summands_.size()}; +} + +Expr::hash_type Sum::memoizing_hash() const { + auto compute_hash = [this]() { + if (summands_.size() == 1) + return summands_[0]->hash_value(); + else { + auto deref_summands = + summands() | + ranges::views::transform( + [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); + auto value = hash::range(ranges::begin(deref_summands), + ranges::end(deref_summands)); + return value; + } + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + + return *hash_value_; +} + +ExprPtr Sum::canonicalize(CanonicalizeOptions opt) { + return canonicalize_impl(true, opt); +} +ExprPtr Sum::rapid_canonicalize(CanonicalizeOptions opts) { + SEQUANT_ASSERT(opts.method == CanonicalizationMethod::Rapid); + return canonicalize_impl(false, opts); +} + +bool Sum::static_equal(const Expr &that) const { + const auto &that_cast = static_cast(that); + if (summands().size() == that_cast.summands().size()) { + if (this->empty()) return true; + // compare hash values first + if (this->hash_value() == + that.hash_value()) // hash values agree -> do full comparison + return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), + expr_ptr_comparer); + else + return false; + } else + return false; +} + +HashingAccumulator &HashingAccumulator::append(ExprPtr summand, bool flatten) { + // flatten, if needed + if (flatten && summand.is()) { + for (auto &subsummand : summand.as().summands()) { + this->append(subsummand, flatten); + } + return *this; + } + + // process summand as a whole + auto it = summands_.find(summand); + if (it == summands_.end()) { + summands_.emplace(summand); + } else { // found existing term with the same hash + auto existing_summand = *it; + if (summand.template is()) { + if (existing_summand.is()) { + // both are products - add them + existing_summand.as().add_identical( + summand.template as()); + } else { + // convert existing term to product and add + auto product_copy = std::make_shared(summand->clone()); + product_copy->add_identical(existing_summand); + summands_.erase(it); + summands_.emplace(std::move(product_copy)); + } + } else { + if (existing_summand.is()) { + existing_summand.as().add_identical(summand); + } else { + // neither is a product - create new product + auto product_form = std::make_shared(); + product_form->append(2, summand.template as()); + summands_.erase(it); + summands_.emplace(std::move(product_form)); + } + } + } + + return *this; +} + +SumPtr HashingAccumulator::make_sum_impl(bool canonicalize) { + Sum::summands_type summands; + summands.reserve(summands_.size()); + for (auto summand : summands_) { + if (!summand->is_zero()) { + summands.push_back(summand); + } + } + + if (canonicalize) { + ranges::sort(summands, [](const auto &e1, const auto &e2) { + if (e1->hash_value() == e2->hash_value()) { + return e1 < e2; + } else { + return e1->hash_value() < e2->hash_value(); + } + }); + } + + return std::make_shared(std::move(summands), Sum::move_only_tag{}); +} + +SumPtr HashingAccumulator::make_sum() { return make_sum_impl(false); } + +SumPtr HashingAccumulator::make_canonicalized_sum() { + return make_sum_impl(true); +} + +ExprPtr HashingAccumulator::make_expr(bool canonicalize) { + if (summands_.size() == 0) { + return ex(0); + } else if (summands_.size() == 1) + return *(summands_.begin()); + else + return make_sum_impl(canonicalize); +} + +bool HashingAccumulator::empty() const { return summands_.empty(); } + +} // namespace sequant diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index 18bb09c873..a0f3186c45 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -37,12 +37,7 @@ class Sum : public Expr { /// construct a Sum out of zero or more summands /// @param summands an initializer list of summands - Sum(ExprPtrList summands) { - // use append to flatten out Sum summands - for (auto &&summand : summands) { - append(std::forward(summand)); - } - } + Sum(ExprPtrList summands); /// construct a Sum out of a range of summands /// @param begin the begin iterator @@ -87,121 +82,30 @@ class Sum : public Expr { /// construct a Sum by moving in the summands, no flattening is performed, /// but zeros will be omitted and constants added up /// @param summands the summands to move in - explicit Sum(summands_type &&summands, move_only_tag) - : summands_(std::move(summands)) { - std::size_t pos = 0; - for (auto it = summands_.begin(); it != summands_.end(); ++it) { - auto &summand = *it; - bool do_erase = false; - if (summand->is_zero()) { - do_erase = true; - } else if (summand->is()) { - auto summand_constant = summand.as_shared_ptr(); - if (constant_summand_idx_) { // add up to the existing constant ... - SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); - summands_[*constant_summand_idx_].as() += *summand_constant; - do_erase = true; - } else { // or memorize the position of the constant - constant_summand_idx_ = pos; - } - } - - // erase if needed - if (do_erase) { - summands_.erase(it); - it = summands_.begin(); - std::advance(it, pos); - } else - ++pos; - } - } + explicit Sum(summands_type &&summands, move_only_tag); /// append a summand to the sum /// @param summand the summand - Sum &append(ExprPtr summand) { - SEQUANT_ASSERT(summand); - if (!summand->is()) { - if (!summand->is_zero()) { // exclude zeros - if (summand->is()) { // add up constants - // immediately, if possible - auto summand_constant = summand.as_shared_ptr(); - if (constant_summand_idx_) { - SEQUANT_ASSERT( - summands_.at(*constant_summand_idx_)->is()); - summands_[*constant_summand_idx_].as() += *summand; - } else { - summands_.push_back(summand->clone()); - constant_summand_idx_ = summands_.size() - 1; - } - } else { - summands_.push_back(summand->clone()); - } - reset_hash_value(); - } - } else { // this recursively flattens Sum summands - for (auto &subsummand : *summand) this->append(subsummand); - } - return *this; - } + Sum &append(ExprPtr summand); /// prepend a summand to the sum /// @param summand the summand - Sum &prepend(ExprPtr summand) { - SEQUANT_ASSERT(summand); - if (!summand->is()) { - if (!summand->is_zero()) { - // exclude zeros - if (summand->is()) { - auto summand_constant = summand.as_shared_ptr(); - if (constant_summand_idx_) { // add up to the existing constant ... - SEQUANT_ASSERT( - summands_.at(*constant_summand_idx_)->is()); - summands_[*constant_summand_idx_].as() += - *summand_constant; - } else { // or include the nonzero constant and update - // constant_summand_idx_ - summands_.insert(summands_.begin(), summand->clone()); - constant_summand_idx_ = 0; - } - } else { - summands_.insert(summands_.begin(), summand->clone()); - if (constant_summand_idx_) // if have a constant, update its position - ++*constant_summand_idx_; - } - reset_hash_value(); - } - } else { // this recursively flattens Sum summands - for (auto &subsummand : *summand) this->prepend(subsummand); - } - return *this; - } + Sum &prepend(ExprPtr summand); /// Summands accessor - const auto &summands() const { return summands_; } + const summands_type &summands() const; /// Summand accessor /// @param i summand index /// @return ith summand - const ExprPtr &summand(size_t i) const { return summands_.at(i); } + const ExprPtr &summand(size_t i) const; /// Takes the first @c count elements of the sum - ExprPtr take_n(size_t count) const { - const auto e = (count >= summands_.size() ? summands_.end() - : (summands_.begin() + count)); - return ex(summands_.begin(), e); - } + ExprPtr take_n(size_t count) const; /// Takes the first @c count elements of the sum starting with element @c /// offset - ExprPtr take_n(size_t offset, size_t count) const { - const auto offset_plus_count = offset + count; - const auto b = (offset >= summands_.size() ? summands_.end() - : (summands_.begin() + offset)); - const auto e = (offset_plus_count >= summands_.size() - ? summands_.end() - : (summands_.begin() + offset_plus_count)); - return ex(b, e); - } + ExprPtr take_n(size_t offset, size_t count) const; /// @tparam Filter a boolean predicate type, such `Filter(const ExprPtr&)` /// evaluates to true @@ -213,82 +117,31 @@ class Sum : public Expr { } /// @return true if the number of factors is zero - bool empty() const { return summands_.empty(); } + bool empty() const; /// @return the number of summands in a Sum - std::size_t size() const { return summands_.size(); } - - std::wstring to_latex() const override { - std::wstring result; - result = L"{ \\bigl("; - std::size_t counter = 0; - for (const auto &i : summands()) { - const auto i_is_product = i->is(); - if (!i_is_product) { - result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); - } else { // i_is_product - const auto i_prod = i->as(); - const auto scalar = i_prod.scalar(); - if (scalar.real() < 0 || (scalar.real() == 0 && scalar.imag() < 0)) { - result += L" - " + i_prod.to_latex(true); - } else { - result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); - } - } - ++counter; - } - result += L"\\bigr) }"; - return result; - } + std::size_t size() const; - Expr::type_id_type type_id() const override { - return Expr::get_type_id(); - }; + std::wstring to_latex() const override; - ExprPtr clone() const override { - auto cloned_summands = - summands() | ranges::views::transform( - [](const ExprPtr &ptr) { return ptr->clone(); }); - return ex(ranges::begin(cloned_summands), - ranges::end(cloned_summands)); - } + Expr::type_id_type type_id() const override; + + ExprPtr clone() const override; /// @brief adjoint of a Sum is a sum of adjoints of its factors virtual void adjoint() override; - Sum &operator+=(const Expr &that) { - this->append(const_cast(that).shared_from_this()); - return *this; - } + Sum &operator+=(const Expr &that); - Sum &operator-=(const Expr &that) { - if (that.is()) - this->append(ex(-that.as().value())); - else - this->append(ex( - -1, ExprPtrList{const_cast(that).shared_from_this()})); - return *this; - } + Sum &operator-=(const Expr &that); - ExprIterator begin_subexpr() override { - if (!summands_.empty()) { - reset_hash_value(); - } + ExprIterator begin_subexpr() override; - return ExprIterator{summands_.data()}; - } + ExprIterator end_subexpr() override; - ExprIterator end_subexpr() override { - return ExprIterator{summands_.data() + summands_.size()}; - } + ConstExprIterator begin_subexpr() const override; - ConstExprIterator begin_subexpr() const override { - return ConstExprIterator{summands_.data()}; - } - - ConstExprIterator end_subexpr() const override { - return ConstExprIterator{summands_.data() + summands_.size()}; - } + ConstExprIterator end_subexpr() const override; private: summands_type summands_{}; @@ -299,61 +152,21 @@ class Sum : public Expr { /// @return the hash of this object /// @note this ensures that hash of a Sum of a single summand is /// identical to the hash of the summand itself. - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - if (summands_.size() == 1) - return summands_[0]->hash_value(); - else { - auto deref_summands = - summands() | - ranges::views::transform( - [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); - auto value = hash::range(ranges::begin(deref_summands), - ranges::end(deref_summands)); - return value; - } - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - - return *hash_value_; - } + hash_type memoizing_hash() const override; /// @param multipass if true, will do a multipass canonicalization, with extra /// cleanup pass after the deep canonization pass ExprPtr canonicalize_impl(bool multipass, CanonicalizeOptions opt); - virtual ExprPtr canonicalize( - CanonicalizeOptions opt = - CanonicalizeOptions::default_options()) override { - return canonicalize_impl(true, opt); - } - virtual ExprPtr rapid_canonicalize( + ExprPtr canonicalize(CanonicalizeOptions opt = + CanonicalizeOptions::default_options()) override; + + ExprPtr rapid_canonicalize( CanonicalizeOptions opts = CanonicalizeOptions::default_options().copy_and_set( - CanonicalizationMethod::Rapid)) override { - SEQUANT_ASSERT(opts.method == CanonicalizationMethod::Rapid); - return canonicalize_impl(false, opts); - } + CanonicalizationMethod::Rapid)) override; - bool static_equal(const Expr &that) const override { - const auto &that_cast = static_cast(that); - if (summands().size() == that_cast.summands().size()) { - if (this->empty()) return true; - // compare hash values first - if (this->hash_value() == - that.hash_value()) // hash values agree -> do full comparison - return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), - expr_ptr_comparer); - else - return false; - } else - return false; - } + bool static_equal(const Expr &that) const override; }; // class Sum /// @brief utility for eagerly accumulating summands in a hash table @@ -375,7 +188,7 @@ class HashingAccumulator { /// zero summands), or the lone summand itself ExprPtr make_expr(bool canonicalize = true); - bool empty() const { return summands_.empty(); } + bool empty() const; private: /// @brief Common implementation for make_sum and make_canonicalized_sum From 8f2f1451d0f5527d9e372482ba822bfcae9bf186 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 11:20:26 +0200 Subject: [PATCH 20/21] Remove redundant headers --- SeQuant/core/expressions/expr.cpp | 23 ++--------------------- 1 file changed, 2 insertions(+), 21 deletions(-) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index e147806d7f..9c8e588150 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -2,29 +2,10 @@ // Created by Eduard Valeyev on 2019-02-06. // -#include -#include #include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include namespace sequant { From 92c8eeb28df801dd6449d844b2477344bb1d7a67 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 12:06:52 +0200 Subject: [PATCH 21/21] Add back missing header --- SeQuant/core/expressions/expr.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index 9c8e588150..cb9655e12b 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -3,6 +3,7 @@ // #include +#include #include #include #include