From 98bb9384ee6199eb647347981b7990ef722a87f8 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Wed, 12 Aug 2026 16:48:31 +0200 Subject: [PATCH 01/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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 From 92354e0723a052bfd9ea9a4fc11fe977e7dd01c9 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 16:46:58 +0200 Subject: [PATCH 22/53] Replace virtual clone with virtual unique_copy The new function returns a unique_ptr instead of an ExprPtr. Reason being that this gives much more flexibility such as moving the ownership of the object out of the smart pointer or simply using as a unique_ptr. Since unique_ptr is implicitly convertible to a shared_ptr (via move ctor), conversion to ExprPtr is trivially possible. To retain compatibility with existing interface, Expr now implements a clone() function by means of the new unique_copy(). --- SeQuant/core/expressions/constant.cpp | 8 ++++++-- SeQuant/core/expressions/constant.hpp | 6 ++++-- SeQuant/core/expressions/expr.cpp | 2 ++ SeQuant/core/expressions/expr.hpp | 4 +++- SeQuant/core/expressions/power.cpp | 14 ++++++++------ SeQuant/core/expressions/power.hpp | 8 ++++++-- SeQuant/core/expressions/product.cpp | 9 ++++----- SeQuant/core/expressions/product.hpp | 9 ++++----- SeQuant/core/expressions/sum.cpp | 17 ++++++++++------- SeQuant/core/expressions/sum.hpp | 6 ++++-- SeQuant/core/expressions/tensor.cpp | 4 ++++ SeQuant/core/expressions/tensor.hpp | 5 +++-- SeQuant/core/expressions/variable.cpp | 5 ++++- SeQuant/core/expressions/variable.hpp | 6 ++++-- SeQuant/core/op.hpp | 21 ++++++++++++++------- SeQuant/domain/mbpt/op.hpp | 7 ++++--- SeQuant/domain/mbpt/op.ipp | 6 ++++-- tests/unit/test_expr.cpp | 12 +++++++----- 18 files changed, 95 insertions(+), 54 deletions(-) diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp index 04bb759c54..d44ac2ddff 100644 --- a/SeQuant/core/expressions/constant.cpp +++ b/SeQuant/core/expressions/constant.cpp @@ -4,6 +4,8 @@ #include #include +#include + namespace sequant { std::wstring Constant::to_latex() const { @@ -14,8 +16,6 @@ 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(); @@ -61,6 +61,10 @@ bool Constant::is_zero(scalar_type v) { return v.is_zero(); } bool Constant::is_zero() const { return is_zero(this->value()); } +std::unique_ptr Constant::unique_copy() const { + return std::make_unique(this->value()); +} + Expr::hash_type Constant::memoizing_hash() const { if (!hash_value_) { hash_value_ = hash::value(value_); diff --git a/SeQuant/core/expressions/constant.hpp b/SeQuant/core/expressions/constant.hpp index 9bade3b542..991ece7b28 100644 --- a/SeQuant/core/expressions/constant.hpp +++ b/SeQuant/core/expressions/constant.hpp @@ -8,6 +8,7 @@ #include +#include #include namespace sequant { @@ -73,8 +74,6 @@ class Constant : public Expr { bool is_scalar() const override; - ExprPtr clone() const override; - /// @brief adjoint of a Constant is its complex conjugate virtual void adjoint() override; @@ -91,6 +90,9 @@ class Constant : public Expr { /// @return `Constant::is_zero(this->value())` bool is_zero() const final; + protected: + std::unique_ptr unique_copy() const override; + private: scalar_type value_; diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index cb9655e12b..55117846af 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -60,6 +60,8 @@ std::wstring Expr::to_latex() const { throw Exception("to_latex not implemented for " + type_name()); } +ExprPtr Expr::clone() const { return unique_copy(); } + bool proportional_to::operator()(const ExprPtr &expr1, const ExprPtr &expr2) const { if (expr1->type_id() != diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index e1f7ba1870..03fe15726f 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -75,7 +75,7 @@ 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 - virtual ExprPtr clone() const = 0; + ExprPtr clone() const; /// like Expr::shared_from_this, but returns ExprPtr /// @return a shared_ptr to this object wrapped into ExprPtr, if this object @@ -418,6 +418,8 @@ class Expr : public std::enable_shared_from_this { return true; } + virtual std::unique_ptr unique_copy() const = 0; + private: /// @return returns next type id in the grand class list static type_id_type get_next_type_id() { diff --git a/SeQuant/core/expressions/power.cpp b/SeQuant/core/expressions/power.cpp index f3920dac24..9e7cfc7450 100644 --- a/SeQuant/core/expressions/power.cpp +++ b/SeQuant/core/expressions/power.cpp @@ -5,6 +5,8 @@ #include #include +#include + namespace sequant { Power::Power(ExprPtr base, exponent_type exponent) @@ -118,12 +120,6 @@ 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) { @@ -158,6 +154,12 @@ Power& Power::operator*=(const Expr& that) { throw Exception("Power::operator*=(that): not valid for that"); } +std::unique_ptr Power::unique_copy() const { + auto cloned = std::make_unique(base_, exponent_); + if (conjugated_) cloned->as().conjugate(); + return cloned; +} + Expr::hash_type Power::memoizing_hash() const { auto compute_hash = [this]() { if (exponent_ == 1 && !conjugated_) return hash::value(*base_); diff --git a/SeQuant/core/expressions/power.hpp b/SeQuant/core/expressions/power.hpp index 29fb69e505..b9f594cf5b 100644 --- a/SeQuant/core/expressions/power.hpp +++ b/SeQuant/core/expressions/power.hpp @@ -7,6 +7,9 @@ #include #include +#include +#include + namespace sequant { /// @brief Represents base^exponent where base is a scalar (Constant or @@ -80,8 +83,6 @@ class Power : public Expr { bool is_scalar() const override; - ExprPtr clone() const override; - /// @brief adjoint of Power: flips the conjugation flag. void adjoint() override; @@ -95,6 +96,9 @@ class Power : public Expr { /// @throw Exception if @p that is not combinable. Power& operator*=(const Expr& that); + protected: + std::unique_ptr unique_copy() const override; + private: ExprPtr base_; exponent_type exponent_; diff --git a/SeQuant/core/expressions/product.cpp b/SeQuant/core/expressions/product.cpp index 48ac8a2684..c312fe05c0 100644 --- a/SeQuant/core/expressions/product.cpp +++ b/SeQuant/core/expressions/product.cpp @@ -221,6 +221,10 @@ ExprPtr Product::rapid_canonicalize(CanonicalizeOptions opt) { bool Product::static_commutativity() const { return false; } +std::unique_ptr Product::unique_copy() const { + return std::make_unique(deep_copy()); +} + std::wstring Product::to_latex() const { return to_latex(false); } std::wstring Product::to_latex(bool negate) const { @@ -251,11 +255,6 @@ 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) { diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 1c12ed2728..53de097186 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -11,6 +11,7 @@ #include +#include #include #include @@ -306,11 +307,6 @@ class Product : public Expr { 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; - Product deep_copy() const; Product &operator*=(const Expr &that); @@ -338,6 +334,9 @@ class Product : public Expr { CanonicalizeOptions::default_options().copy_and_set( CanonicalizationMethod::Rapid)) override; + protected: + std::unique_ptr unique_copy() const override; + private: scalar_type scalar_ = {1, 0}; factors_type factors_{}; diff --git a/SeQuant/core/expressions/sum.cpp b/SeQuant/core/expressions/sum.cpp index 880c47d2dc..13133f4102 100644 --- a/SeQuant/core/expressions/sum.cpp +++ b/SeQuant/core/expressions/sum.cpp @@ -4,6 +4,8 @@ #include #include +#include + namespace sequant { Sum::Sum(ExprPtrList summands) { @@ -144,13 +146,6 @@ std::wstring Sum::to_latex() const { 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) { @@ -254,6 +249,14 @@ ConstExprIterator Sum::end_subexpr() const { return ConstExprIterator{summands_.data() + summands_.size()}; } +std::unique_ptr Sum::unique_copy() const { + auto cloned_summands = + summands() | + ranges::views::transform([](const ExprPtr &ptr) { return ptr->clone(); }); + return std::make_unique(ranges::begin(cloned_summands), + ranges::end(cloned_summands)); +} + Expr::hash_type Sum::memoizing_hash() const { auto compute_hash = [this]() { if (summands_.size() == 1) diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index a0f3186c45..0f5e9f85ba 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -126,8 +127,6 @@ class Sum : public Expr { 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; @@ -143,6 +142,9 @@ class Sum : public Expr { ConstExprIterator end_subexpr() const override; + protected: + std::unique_ptr unique_copy() const override; + private: summands_type summands_{}; std::optional diff --git a/SeQuant/core/expressions/tensor.cpp b/SeQuant/core/expressions/tensor.cpp index 0495e2dd22..ea546c8c43 100644 --- a/SeQuant/core/expressions/tensor.cpp +++ b/SeQuant/core/expressions/tensor.cpp @@ -45,4 +45,8 @@ ExprPtr Tensor::canonicalize(CanonicalizeOptions) { return canonicalizer_ptr ? canonicalizer_ptr->apply(*this) : ExprPtr{}; } +std::unique_ptr Tensor::unique_copy() const { + return std::make_unique(*this); +} + } // namespace sequant diff --git a/SeQuant/core/expressions/tensor.hpp b/SeQuant/core/expressions/tensor.hpp index 07b3b7d5c2..9023bc0af8 100644 --- a/SeQuant/core/expressions/tensor.hpp +++ b/SeQuant/core/expressions/tensor.hpp @@ -754,8 +754,6 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { type_id_type type_id() const override { return get_type_id(); }; - ExprPtr clone() const override { return ex(*this); } - void reset_tags() const { ranges::for_each(slots(), [](const auto &idx) { idx.reset_tag(); }); } @@ -775,6 +773,9 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { return false; // TODO do we compare typeid? labels? probably the latter } + protected: + std::unique_ptr unique_copy() const override; + private: std::wstring label_{}; sequant::bra bra_{}; diff --git a/SeQuant/core/expressions/variable.cpp b/SeQuant/core/expressions/variable.cpp index 58b15d00e9..15d614f478 100644 --- a/SeQuant/core/expressions/variable.cpp +++ b/SeQuant/core/expressions/variable.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -60,7 +61,9 @@ std::wstring Variable::to_latex() const { return result; } -ExprPtr Variable::clone() const { return ex(*this); } +std::unique_ptr Variable::unique_copy() const { + return std::make_unique(*this); +} void Variable::adjoint() { conjugate(); } diff --git a/SeQuant/core/expressions/variable.hpp b/SeQuant/core/expressions/variable.hpp index 4943845ac4..de13956ae3 100644 --- a/SeQuant/core/expressions/variable.hpp +++ b/SeQuant/core/expressions/variable.hpp @@ -5,6 +5,7 @@ #include #include +#include #include #include @@ -50,11 +51,12 @@ class Variable : public Expr, public MutatableLabeled { bool is_scalar() const override; - ExprPtr clone() const override; - /// @brief adjoint of a Variable is its complex conjugate virtual void adjoint() override; + protected: + std::unique_ptr unique_copy() const override; + private: std::wstring label_; bool conjugated_ = false; diff --git a/SeQuant/core/op.hpp b/SeQuant/core/op.hpp index 4cc47e993d..b304f39853 100644 --- a/SeQuant/core/op.hpp +++ b/SeQuant/core/op.hpp @@ -389,7 +389,10 @@ class Operator : public container::svector>, public Expr { type_id_type type_id() const override { return get_type_id(); }; - ExprPtr clone() const override { return std::make_shared(*this); } + protected: + std::unique_ptr unique_copy() const override { + return std::make_unique(*this); + } private: base_type make_ops(Action action, IndexList indices) { @@ -711,10 +714,6 @@ class NormalOperator : public Operator, return Expr::get_type_id(); }; - ExprPtr clone() const override { - return std::make_shared(*this); - } - virtual void adjoint() override { // same as base adjoint(), but updates extra state Operator::adjoint(); @@ -739,6 +738,11 @@ class NormalOperator : public Operator, return mutated; } + protected: + std::unique_ptr unique_copy() const override { + return std::make_unique(*this); + } + private: Vacuum vacuum_; std::size_t ncreators_ = 0; @@ -1068,8 +1072,6 @@ 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() && @@ -1077,6 +1079,11 @@ class NormalOperatorSequence : public container::svector>, static_cast(nopseq2); } + protected: + std::unique_ptr unique_copy() const override { + return std::make_unique(*this); + } + private: Vacuum vacuum_ = Vacuum::Physical; /// ensures that all operators use same vacuum, and sets vacuum_ diff --git a/SeQuant/domain/mbpt/op.hpp b/SeQuant/domain/mbpt/op.hpp index 658e1bacba..45e909f841 100644 --- a/SeQuant/domain/mbpt/op.hpp +++ b/SeQuant/domain/mbpt/op.hpp @@ -43,7 +43,7 @@ #include #include #include -#include +#include #include #include #include @@ -950,6 +950,9 @@ class Operator : public Operator { /// @brief returns the perturbation order of this operator [[nodiscard]] size_t order() const { return order_; } + protected: + std::unique_ptr unique_copy() const override; + private: std::function qn_action_; @@ -963,8 +966,6 @@ class Operator : public Operator { Expr::type_id_type type_id() const override; - ExprPtr clone() const override; - std::wstring to_latex() const override; Expr::hash_type memoizing_hash() const override; diff --git a/SeQuant/domain/mbpt/op.ipp b/SeQuant/domain/mbpt/op.ipp index 8f28062b51..2030799552 100644 --- a/SeQuant/domain/mbpt/op.ipp +++ b/SeQuant/domain/mbpt/op.ipp @@ -12,6 +12,8 @@ #include #include +#include + namespace sequant { namespace mbpt { @@ -192,8 +194,8 @@ Expr::type_id_type Operator::type_id() const { }; template -ExprPtr Operator::clone() const { - return ex(*this); +std::unique_ptr Operator::unique_copy() const { + return std::make_unique(*this); } // Expresses general operators in human interpretable form. for example: diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index e9a1e0924d..bd48e26445 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -37,7 +37,9 @@ struct Dummy : public sequant::Expr { virtual ~Dummy() = default; 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(); } + std::unique_ptr unique_copy() const override { + return std::make_unique(); + } void adjoint() override {} bool static_equal(const sequant::Expr &) const override { return true; } }; @@ -110,8 +112,8 @@ struct VecExpr : public std::vector, public sequant::Expr { static_cast(static_cast(that)); } - sequant::ExprPtr clone() const override { - return sequant::ex(this->begin(), this->end()); + std::unique_ptr unique_copy() const override { + return std::make_unique(this->begin(), this->end()); } }; @@ -123,8 +125,8 @@ struct Adjointable : public sequant::Expr { return L"{\\text{Adjointable}{" + std::to_wstring(v) + L"}}"; } type_id_type type_id() const override { return get_type_id(); }; - sequant::ExprPtr clone() const override { - return sequant::ex(v); + std::unique_ptr unique_copy() const override { + return std::make_unique(v); } bool static_equal(const sequant::Expr &that) const override { return v == that.as().v; From 4e150e46e5054967ee048513b8fb4c25565bb2bd Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 16:57:00 +0200 Subject: [PATCH 23/53] Add more arithmetic operators for Constants --- SeQuant/core/expressions/constant.cpp | 60 +++++++++++++++++++++------ SeQuant/core/expressions/constant.hpp | 7 ++++ 2 files changed, 55 insertions(+), 12 deletions(-) diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp index d44ac2ddff..7bf3854572 100644 --- a/SeQuant/core/expressions/constant.cpp +++ b/SeQuant/core/expressions/constant.cpp @@ -21,40 +21,52 @@ void Constant::adjoint() { reset_hash_value(); } +Constant &Constant::operator*=(const Constant &that) { + value_ *= that.value(); + + reset_hash_value(); + + return *this; +} + Constant &Constant::operator*=(const Expr &that) { - if (that.is()) { - value_ *= that.as().value(); - } else { + if (!that.is()) { throw Exception("Constant::operator*=(that): not valid for that"); } + return *this *= that.as(); +} + +Constant &Constant::operator+=(const Constant &that) { + value_ += that.value(); + reset_hash_value(); return *this; } Constant &Constant::operator+=(const Expr &that) { - if (that.is()) { - value_ += that.as().value(); - } else { + if (!that.is()) { throw Exception("Constant::operator+=(that): not valid for that"); } + return *this += that.as(); +} + +Constant &Constant::operator-=(const Constant &that) { + value_ -= that.value(); + reset_hash_value(); return *this; } Constant &Constant::operator-=(const Expr &that) { - if (that.is()) { - value_ -= that.as().value(); - } else { + if (!that.is()) { throw Exception("Constant::operator-=(that): not valid for that"); } - reset_hash_value(); - - return *this; + return *this -= that.as(); } bool Constant::is_zero(scalar_type v) { return v.is_zero(); } @@ -78,4 +90,28 @@ bool Constant::static_equal(const Expr &that) const { return value() == static_cast(that).value(); } +Constant operator*(const Constant &lhs, const Constant &rhs) { + Constant result(lhs); + + result *= rhs; + + return result; +} + +Constant operator+(const Constant &lhs, const Constant &rhs) { + Constant result(lhs); + + result += rhs; + + return result; +} + +Constant operator-(const Constant &lhs, const Constant &rhs) { + Constant result(lhs); + + result -= rhs; + + return result; +} + } // namespace sequant diff --git a/SeQuant/core/expressions/constant.hpp b/SeQuant/core/expressions/constant.hpp index 991ece7b28..9d433f32d2 100644 --- a/SeQuant/core/expressions/constant.hpp +++ b/SeQuant/core/expressions/constant.hpp @@ -77,10 +77,13 @@ class Constant : public Expr { /// @brief adjoint of a Constant is its complex conjugate virtual void adjoint() override; + Constant &operator*=(const Constant &that); Constant &operator*=(const Expr &that); + Constant &operator+=(const Constant &that); Constant &operator+=(const Expr &that); + Constant &operator-=(const Constant &that); Constant &operator-=(const Expr &that); /// @param[in] v a scalar @@ -102,6 +105,10 @@ class Constant : public Expr { }; // class Constant +Constant operator*(const Constant &lhs, const Constant &rhs); +Constant operator+(const Constant &lhs, const Constant &rhs); +Constant operator-(const Constant &lhs, const Constant &rhs); + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_CONSTANT_HPP From 0f5a767e2071239eb2fa17a717342eeba528db94 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 17:18:10 +0200 Subject: [PATCH 24/53] Fix unused variable warning --- SeQuant/domain/mbpt/op.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SeQuant/domain/mbpt/op.cpp b/SeQuant/domain/mbpt/op.cpp index 50f7b295d4..171981c375 100644 --- a/SeQuant/domain/mbpt/op.cpp +++ b/SeQuant/domain/mbpt/op.cpp @@ -580,14 +580,14 @@ ExprPtr OpMaker::operator()( if (!dep && csv) { if (opclass == OpClass::Ex) { if constexpr (assert_enabled()) { - for (auto&& s : cre_spaces_) { + for ([[maybe_unused]] const auto& s : cre_spaces_) { SEQUANT_ASSERT(isr->contains_unoccupied(s)); } } dep = UseDepIdx::Bra; } else if (opclass == OpClass::Deex) { if constexpr (assert_enabled()) { - for (auto&& s : ann_spaces_) { + for ([[maybe_unused]] const auto& s : ann_spaces_) { SEQUANT_ASSERT(isr->contains_unoccupied(s)); } } From 09cace5d5dd54cc8560716639b80eeff17d0cbfd Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 19:17:05 +0200 Subject: [PATCH 25/53] Add and use Expr & overloads --- SeQuant/core/expressions/expr_algorithms.cpp | 10 +++++++--- SeQuant/core/expressions/expr_algorithms.hpp | 2 ++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/SeQuant/core/expressions/expr_algorithms.cpp b/SeQuant/core/expressions/expr_algorithms.cpp index 982f2f09dd..df21424cbc 100644 --- a/SeQuant/core/expressions/expr_algorithms.cpp +++ b/SeQuant/core/expressions/expr_algorithms.cpp @@ -19,10 +19,10 @@ namespace sequant { -std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align, +std::wstring to_latex_align(const Expr& expr, size_t max_lines_per_align, size_t max_terms_per_line) { - std::wstring result = io::latex::to_string(exprptr); - if (exprptr->is()) { + std::wstring result = io::latex::to_string(expr); + if (expr.is()) { result.erase(0, 7); // remove leading "{ \bigl" result.replace(result.size() - 8, 8, L")"); // replace trailing "\bigr) }" with ")" @@ -78,6 +78,10 @@ std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align, result += L"\n\\end{align}"; return result; } +std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align, + size_t max_terms_per_line) { + return to_latex_align(*exprptr, max_lines_per_align, max_terms_per_line); +} std::size_t size(const Expr& expr) { return ranges::size(expr); } diff --git a/SeQuant/core/expressions/expr_algorithms.hpp b/SeQuant/core/expressions/expr_algorithms.hpp index 1bd0454e54..c21a9deade 100644 --- a/SeQuant/core/expressions/expr_algorithms.hpp +++ b/SeQuant/core/expressions/expr_algorithms.hpp @@ -22,6 +22,8 @@ namespace sequant { /// @param max_lines_per_align the maximum number of lines in the align before /// starting new align block (if zero, will produce single align block) /// @param max_terms_per_line the maximum number of terms per line +std::wstring to_latex_align(const Expr& expr, size_t max_lines_per_align = 0, + size_t max_terms_per_line = 1); std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align = 0, size_t max_terms_per_line = 1); From 93dda4a4f838c8f6b24590d7d627c84c2e0e7fef Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 19:18:41 +0200 Subject: [PATCH 26/53] Deprecate {shared,weak}_from_this This is in the way of having expression objects that are not managed by shared_ptr --- SeQuant/core/expressions/expr.cpp | 18 +++++++++ SeQuant/core/expressions/expr.hpp | 59 ++++++++++++++++++++++------ SeQuant/core/expressions/product.cpp | 2 +- SeQuant/core/expressions/product.hpp | 24 +++++------ SeQuant/core/expressions/sum.cpp | 15 ++++--- SeQuant/core/expressions/sum.hpp | 2 +- 6 files changed, 87 insertions(+), 33 deletions(-) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index 55117846af..34b0c4ed54 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -8,6 +8,8 @@ #include #include +#include + namespace sequant { ExprIterator Expr::begin() { return begin_subexpr(); } @@ -62,6 +64,22 @@ std::wstring Expr::to_latex() const { ExprPtr Expr::clone() const { return unique_copy(); } +std::shared_ptr Expr::shared_from_this() { + return std::enable_shared_from_this::shared_from_this(); +} + +std::shared_ptr Expr::shared_from_this() const { + return std::enable_shared_from_this::shared_from_this(); +} + +std::weak_ptr Expr::weak_from_this() { + return std::enable_shared_from_this::weak_from_this(); +} + +std::weak_ptr Expr::weak_from_this() const { + return std::enable_shared_from_this::weak_from_this(); +} + bool proportional_to::operator()(const ExprPtr &expr1, const ExprPtr &expr2) const { if (expr1->type_id() != diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 03fe15726f..b95dab52a8 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -77,11 +77,25 @@ class Expr : public std::enable_shared_from_this { /// @return a clone of this object, i.e. an object that is equal to @c this ExprPtr clone() const; + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + shared_ptr + shared_from_this(); + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + shared_ptr + shared_from_this() const; + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + weak_ptr + weak_from_this(); + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + weak_ptr + weak_from_this() const; + /// like Expr::shared_from_this, but returns ExprPtr /// @return a shared_ptr to this object wrapped into ExprPtr, if this object /// is already managed by a shared_ptr, else returns a shared_ptr to a clone /// of this object wrapped into ExprPtr - ExprPtr exprptr_from_this() { + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] ExprPtr + exprptr_from_this() { if (weak_from_this().use_count() == 0) return this->clone(); else @@ -92,7 +106,8 @@ class Expr : public std::enable_shared_from_this { /// @return a shared_ptr to this object wrapped into ExprPtr, if this object /// is already managed by a shared_ptr, else returns a shared_ptr to a clone /// of this object wrapped into ExprPtr - ExprPtr exprptr_from_this() const { + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] ExprPtr + exprptr_from_this() const { if (weak_from_this().use_count() == 0) return this->clone(); else @@ -237,9 +252,16 @@ class Expr : public std::enable_shared_from_this { /// Expr::memoizing_hash /// @return the hash value for this Expr hash_type hash_value( - std::function &)> hasher = {}) - const { - return hasher ? hasher(shared_from_this()) : memoizing_hash(); + std::function hasher = {}) const { + return hasher ? hasher(*this) : memoizing_hash(); + } + + [[deprecated( + "Use a hashing function that takes a const Expr & instead of " + "shared_ptr")]] hash_type + hash_value(std::function &)> + hasher) const { + return hasher ? hasher(this->clone()) : memoizing_hash(); } /// Computes and returns the derived type identifier @@ -351,9 +373,10 @@ class Expr : public std::enable_shared_from_this { typename E, typename Visitor, typename = std::enable_if_t, Expr>>> static bool visit_impl(E &&expr, Visitor &&visitor, const bool atoms_only) { - if (expr.weak_from_this().use_count() == 0) - throw Exception( - "Expr::visit: cannot visit expressions not managed by shared_ptr"); + constexpr bool visitor_uses_exprptr = + std::is_invocable_r_v, + ExprPtr &>; + for (auto &subexpr_ptr : expr.expr()) { const auto subexpr_is_an_atom = subexpr_ptr->is_atom(); const auto need_to_visit_subexpr = !atoms_only || subexpr_is_an_atom; @@ -362,18 +385,32 @@ class Expr : public std::enable_shared_from_this { visited = visit_impl(*subexpr_ptr, std::forward(visitor), atoms_only); // call on the subexpression itself, if not yet done so - if (need_to_visit_subexpr && !visited) visitor(subexpr_ptr); + if (need_to_visit_subexpr && !visited) { + if constexpr (visitor_uses_exprptr) { + visitor(subexpr_ptr); + } else { + visitor(*subexpr_ptr); + } + } } + // N.B. can only visit itself if visitor is nonmutating! bool this_visited = false; if constexpr (std::is_invocable_r_v, const ExprPtr &>) { if (!atoms_only || expr.is_atom()) { - const ExprPtr this_exprptr = expr.exprptr_from_this(); - visitor(this_exprptr); + visitor(expr.clone()); + this_visited = true; + } + } else if constexpr (std::is_invocable_r_v, + const Expr &>) { + if (!atoms_only || expr.is_atom()) { + visitor(std::as_const(expr)); this_visited = true; } } + return this_visited; } diff --git a/SeQuant/core/expressions/product.cpp b/SeQuant/core/expressions/product.cpp index c312fe05c0..e904eba3c0 100644 --- a/SeQuant/core/expressions/product.cpp +++ b/SeQuant/core/expressions/product.cpp @@ -269,7 +269,7 @@ Product Product::deep_copy() const { Product &Product::operator*=(const Expr &that) { if (!that.is()) { - this->append(1, const_cast(that).shared_from_this()); + this->append(1, that.clone()); } else { scalar_ *= that.as().value(); } diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 53de097186..d40c957a61 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -61,7 +61,7 @@ class Product : public Expr { if constexpr (rng_is_expr || rng_is_exprptr) { ExprPtr rng_as_exprptr; if constexpr (rng_is_expr) { - rng_as_exprptr = rng.exprptr_from_this(); + rng_as_exprptr = rng.clone(); } else { rng_as_exprptr = rng; } @@ -180,10 +180,10 @@ class Product : public Expr { typename = std::enable_if_t>> Product &append(T scalar, Factor &&factor, Flatten flatten_tag = Flatten::Yes) { - return this->append(scalar, - std::static_pointer_cast( - std::forward(factor).shared_from_this()), - flatten_tag); + return this->append( + scalar, + std::static_pointer_cast(std::forward(factor).clone()), + flatten_tag); } /// (post-)multiplies the product by@c factor @@ -199,9 +199,9 @@ class Product : public Expr { /// @warning if @p factor is a Product, it is flattened recursively template >> Product &append(Factor &&factor, Flatten flatten_tag = Flatten::Yes) { - return this->append(std::static_pointer_cast( - std::forward(factor).shared_from_this()), - flatten_tag); + return this->append( + std::static_pointer_cast(std::forward(factor).clone()), + flatten_tag); } /// (pre-)multiplies the product by @c scalar times @c factor @@ -252,10 +252,10 @@ class Product : public Expr { typename = std::enable_if_t>> Product &prepend(T scalar, Factor &&factor, Flatten flatten_tag = Flatten::Yes) { - return this->prepend(scalar, - std::static_pointer_cast( - std::forward(factor).shared_from_this()), - flatten_tag); + return this->prepend( + scalar, + std::static_pointer_cast(std::forward(factor).clone()), + flatten_tag); } const scalar_type &scalar() const; diff --git a/SeQuant/core/expressions/sum.cpp b/SeQuant/core/expressions/sum.cpp index 13133f4102..51fbaa7c81 100644 --- a/SeQuant/core/expressions/sum.cpp +++ b/SeQuant/core/expressions/sum.cpp @@ -156,8 +156,8 @@ void Sum::adjoint() { 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; + std::wcout << "Sum::canonicalize_impl: input = " << to_latex_align(*this) + << std::endl; const auto npasses = multipass ? 2 : 1; for (auto pass = 0; pass != npasses; ++pass) { @@ -191,7 +191,7 @@ ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { if (Logger::instance().canonicalize) std::wcout << "Sum::canonicalize_impl (pass=" << pass << "): after canonicalizing summands = " - << to_latex_align(shared_from_this()) << std::endl; + << to_latex_align(*this) << std::endl; HashingAccumulator acc; for (auto &summand : summands_) { @@ -208,15 +208,15 @@ ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { if (Logger::instance().canonicalize) std::wcout << "Sum::canonicalize_impl (pass=" << pass - << "): after reducing summands = " - << to_latex_align(shared_from_this()) << std::endl; + << "): after reducing summands = " << to_latex_align(*this) + << std::endl; } return {}; // side effects are absorbed into summands } Sum &Sum::operator+=(const Expr &that) { - this->append(const_cast(that).shared_from_this()); + this->append(that.clone()); return *this; } @@ -224,8 +224,7 @@ 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()})); + this->append(ex(-1, ExprPtrList{that.clone()})); return *this; } diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index 0f5e9f85ba..d050a933d8 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -65,7 +65,7 @@ class Sum : public Expr { if constexpr (rng_is_expr || rng_is_exprptr) { ExprPtr rng_as_exprptr; if constexpr (rng_is_expr) { - rng_as_exprptr = rng.exprptr_from_this(); + rng_as_exprptr = rng.clone(); } else { rng_as_exprptr = rng; } From d644b7a635d0efc45af0f3a09079480ac82123c4 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Mon, 17 Aug 2026 19:20:31 +0200 Subject: [PATCH 27/53] Implement ExprContainer --- CMakeLists.txt | 2 + SeQuant/core/expr.hpp | 1 + SeQuant/core/expressions/expr.hpp | 2 + SeQuant/core/expressions/expr_container.cpp | 133 ++++++++++++++++++++ SeQuant/core/expressions/expr_container.hpp | 58 +++++++++ tests/unit/test_expr.cpp | 120 ++++++++++++++++++ 6 files changed, 316 insertions(+) create mode 100644 SeQuant/core/expressions/expr_container.cpp create mode 100644 SeQuant/core/expressions/expr_container.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f119f10dbc..1083c18082 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -312,6 +312,8 @@ set(SeQuant_symb_src SeQuant/core/expressions/expr.hpp SeQuant/core/expressions/expr_algorithms.cpp SeQuant/core/expressions/expr_algorithms.hpp + SeQuant/core/expressions/expr_container.cpp + SeQuant/core/expressions/expr_container.hpp SeQuant/core/expressions/expr_operators.hpp SeQuant/core/expressions/expr_ptr.cpp SeQuant/core/expressions/expr_ptr.hpp diff --git a/SeQuant/core/expr.hpp b/SeQuant/core/expr.hpp index 9c8e4ed763..31a09b779a 100644 --- a/SeQuant/core/expr.hpp +++ b/SeQuant/core/expr.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index b95dab52a8..8b1d98dc44 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -59,6 +59,8 @@ static const wchar_t adjoint_label = L'\u207A'; /// @endcode class Expr : public std::enable_shared_from_this { public: + friend class ExprContainer; + using hash_type = std::size_t; using type_id_type = int; // to speed up comparisons diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp new file mode 100644 index 0000000000..c69019fea2 --- /dev/null +++ b/SeQuant/core/expressions/expr_container.cpp @@ -0,0 +1,133 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace sequant { + +ExprPtr to_expr_ptr(ExprContainer &&container) { + return std::shared_ptr(std::move(container.expr_)); +} + +ExprContainer::ExprContainer(const ExprContainer &container) + : ExprContainer(container->unique_copy()) {} + +ExprContainer::ExprContainer(std::unique_ptr expr) + : expr_(std::move(expr)) {} + +ExprContainer::ExprContainer(const Expr &expr) + : ExprContainer(expr.unique_copy()) {} + +ExprContainer::ExprContainer(Expr &&expr) + : ExprContainer(std::move(expr).unique_copy()) {} + +ExprContainer &ExprContainer::operator=(const ExprContainer &container) { + // Copy-and-swap + ExprContainer copy(container); + + swap(*this, copy); + + return *this; +} + +ExprContainer &ExprContainer::operator=(const Expr &expr) { + // Copy-and-swap + ExprContainer copy(expr); + + swap(*this, copy); + + return *this; +} + +ExprContainer &ExprContainer::operator=(Expr &&expr) { + // Copy-and-swap + ExprContainer copy(std::move(expr)); + + swap(*this, copy); + + return *this; +} + +ExprContainer ExprContainer::copy() const { return expr_->unique_copy(); } + +ExprContainer::operator const Expr &() const { return *expr_; } + +ExprContainer::operator Expr &() & { return *expr_; } +ExprContainer::operator Expr &&() && { return std::move(*expr_); } + +const Expr &ExprContainer::operator*() const { return *expr_; } + +Expr &ExprContainer::operator*() { return *expr_; } + +const Expr *ExprContainer::operator->() const { return expr_.get(); } + +Expr *ExprContainer::operator->() { return expr_.get(); } + +ExprContainer &ExprContainer::operator+=(const Expr &expr) { + if (expr_->is()) { + expr_->as() += expr; + } else if (expr_->is() && expr.is()) { + *this = expr_->as() + expr.as(); + } else { + *this = Sum(ExprPtrList{to_expr_ptr(std::move(*this)), expr.clone()}); + } + + return *this; +} + +ExprContainer &ExprContainer::operator-=(const Expr &expr) { + if (expr_->is()) { + expr_->as() -= expr; + } else if (expr_->is() && expr.is()) { + *this = expr_->as() - expr.as(); + } else { + *this = Sum(ExprPtrList{to_expr_ptr(std::move(*this)), + ex(-1, ExprPtrList{expr.clone()})}); + } + + return *this; +} + +ExprContainer &ExprContainer::operator*=(const Expr &expr) { + if (expr_->is()) { + expr_->as() *= expr; + } else if (expr_->is() && expr.is()) { + *this = expr_->as() * expr.as(); + } else { + *this = Product(ExprPtrList{to_expr_ptr(std::move(*this)), expr.clone()}); + } + return *this; +} + +void swap(ExprContainer &lhs, ExprContainer &rhs) { + std::swap(lhs.expr_, rhs.expr_); +} + +ExprContainer operator+(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont += rhs; + + return cont; +} + +ExprContainer operator-(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont -= rhs; + + return cont; +} + +ExprContainer operator*(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont *= rhs; + + return cont; +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp new file mode 100644 index 0000000000..0a74a80772 --- /dev/null +++ b/SeQuant/core/expressions/expr_container.hpp @@ -0,0 +1,58 @@ +#ifndef SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP +#define SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP + +#include + +#include + +namespace sequant { + +class ExprContainer { + public: + explicit ExprContainer(const ExprContainer &container); + ExprContainer(ExprContainer &&container) = default; + + explicit ExprContainer(const Expr &expr); + ExprContainer(Expr &&expr); + + ExprContainer &operator=(const ExprContainer &container); + ExprContainer &operator=(ExprContainer &&container) = default; + + ExprContainer &operator=(const Expr &expr); + ExprContainer &operator=(Expr &&expr); + + ~ExprContainer() = default; + + ExprContainer copy() const; + + operator const Expr &() const; + operator Expr &() &; + operator Expr &&() &&; + + const Expr &operator*() const; + Expr &operator*(); + + const Expr *operator->() const; + Expr *operator->(); + + ExprContainer &operator+=(const Expr &expr); + ExprContainer &operator-=(const Expr &expr); + ExprContainer &operator*=(const Expr &expr); + + friend void swap(ExprContainer &, ExprContainer &); + + private: + std::unique_ptr expr_; + + ExprContainer(std::unique_ptr expr); + + friend ExprPtr to_expr_ptr(ExprContainer &&container); +}; + +ExprContainer operator+(const Expr &lhs, const Expr &rhs); +ExprContainer operator-(const Expr &lhs, const Expr &rhs); +ExprContainer operator*(const Expr &lhs, const Expr &rhs); + +} // namespace sequant + +#endif // SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index bd48e26445..73d4e00189 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -1119,6 +1119,126 @@ TEST_CASE("expr", "[elements]") { } } + SECTION("ExprContainer") { + SECTION("Constructors") { + SECTION("from conrete") { + ExprContainer cont1(Constant(1)); + REQUIRE(cont1->is()); + REQUIRE(cont1->as() == Constant(1)); + + ExprContainer cont2(Variable("bla")); + REQUIRE(cont2->is()); + REQUIRE(cont2->as() == Variable("bla")); + + Product prod(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont3(prod); + REQUIRE(cont3->is()); + REQUIRE(cont3->as() == prod); + + Sum sum(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont4(sum); + REQUIRE(cont4->is()); + REQUIRE(cont4->as() == sum); + } + SECTION("from base") { + ExprPtr expr = ex(42); + + ExprContainer cont(*expr); + REQUIRE(cont->is()); + REQUIRE(cont->as().value() == 42); + } + SECTION("conversion via 'assignment'") { + ExprContainer cont1 = Constant(1); + REQUIRE(cont1->is()); + REQUIRE(cont1->as() == Constant(1)); + + ExprContainer cont2 = Variable("bla"); + REQUIRE(cont2->is()); + REQUIRE(cont2->as() == Variable("bla")); + + // Note: copy-ctor is explicit so in order for this "assignment" to + // work, we need to assign rvalues + Product prod(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont3 = Product(prod); + REQUIRE(cont3->is()); + REQUIRE(cont3->as() == prod); + + Sum sum(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont4 = Sum(sum); + REQUIRE(cont4->is()); + REQUIRE(cont4->as() == sum); + } + } + SECTION("Assignment") { + ExprContainer cont = Constant(1); + REQUIRE(cont->is()); + REQUIRE(cont->as() == Constant(1)); + + cont = Variable("bla"); + REQUIRE(cont->is()); + REQUIRE(cont->as() == Variable("bla")); + + Product prod(ExprPtrList{ex("bla"), ex(2)}); + cont = Product(prod); + REQUIRE(cont->is()); + REQUIRE(cont->as() == prod); + + Sum sum(ExprPtrList{ex("bla"), ex(2)}); + cont = Sum(sum); + REQUIRE(cont->is()); + REQUIRE(cont->as() == sum); + } + SECTION("value semantics") { + ExprContainer cont = Constant(1); + ExprContainer copy(cont); + copy = Variable("test"); + + REQUIRE(copy->is()); + REQUIRE(cont->is()); + } + SECTION("conversion to Expr &") { + bool passed1 = false; + + auto func1 = [&passed1](Expr &) { passed1 = true; }; + + ExprContainer expr = Constant(5); + func1(expr); + + REQUIRE(passed1); + + bool passed2 = false; + auto func2 = [&passed2](const Expr &) { passed2 = true; }; + + func2(std::as_const(expr)); + + REQUIRE(passed2); + + bool passed3 = false; + auto func3 = [&passed3](Expr &&) { passed3 = true; }; + + func3(std::move(expr)); + + REQUIRE(passed3); + } + SECTION("freestanding Expr arithmetic") { + // This allows to use arithmetic directly on Expr & instances (instead of + // requiring ExprPtr or ExprContainer wrappers) + ExprContainer res = Constant(1) + Variable("One"); + REQUIRE_THAT(res, EquivalentTo("1 + One")); + + res = Variable("A") * Tensor("T", bra({"a1"}), ket()) - Constant(42); + REQUIRE_THAT(res, EquivalentTo("A * T{a1} - 42")); + } + SECTION("In-place ExprContainer arithmetic") { + ExprContainer res = Constant(1); + res += Variable("A"); + res -= Variable("B"); + res *= Constant(3); + + REQUIRE_THAT(res, EquivalentTo("(1 + A - B) * 3")); + } + } + SECTION("ResultExpr") { SECTION("accessors") { SECTION("as_variable") { From 08b1d514cf5b1de0347785538a154932411f91a9 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 09:55:45 +0200 Subject: [PATCH 28/53] Suppress deprecated warnings --- SeQuant/core/expressions/expr.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 8b1d98dc44..7cda1df709 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -98,10 +98,12 @@ class Expr : public std::enable_shared_from_this { /// of this object wrapped into ExprPtr [[deprecated("Expr objects may no longer be managed by shared_ptr")]] ExprPtr exprptr_from_this() { + SEQUANT_PRAGMA_IGNORE_DEPRECATED_BEGIN if (weak_from_this().use_count() == 0) return this->clone(); else return static_cast(this->shared_from_this()); + SEQUANT_PRAGMA_IGNORE_DEPRECATED_END } /// like Expr::shared_from_this, but returns ExprPtr @@ -110,11 +112,13 @@ class Expr : public std::enable_shared_from_this { /// of this object wrapped into ExprPtr [[deprecated("Expr objects may no longer be managed by shared_ptr")]] ExprPtr exprptr_from_this() const { + SEQUANT_PRAGMA_IGNORE_DEPRECATED_BEGIN if (weak_from_this().use_count() == 0) return this->clone(); else return static_cast( std::const_pointer_cast(this->shared_from_this())); + SEQUANT_PRAGMA_IGNORE_DEPRECATED_END } /// Canonicalizes @c this and returns the byproduct of canonicalization (e.g. From c5fcbbd932e7d25f0ea1efadfa7b6f86c736e0af Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 10:32:59 +0200 Subject: [PATCH 29/53] Implement support for non-commuting product for ExprContainer --- SeQuant/core/expressions/expr_container.cpp | 25 +++++++++++++++++++++ SeQuant/core/expressions/expr_container.hpp | 2 ++ 2 files changed, 27 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index c69019fea2..2047ff5605 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -105,6 +105,24 @@ ExprContainer &ExprContainer::operator*=(const Expr &expr) { return *this; } +ExprContainer &ExprContainer::operator^=(const Expr &other) { + auto this_is_product = expr_->is(); + auto other_is_product = other.is(); + if (!this_is_product && !other_is_product) { + *this = + NCProduct(ExprPtrList{to_expr_ptr(std::move(*this)), other.clone()}); + } else if (this_is_product) { + *this = NCProduct(std::move(expr_->as())); + expr_->as().append(1, other.clone()); + } else { // other_is_product + NCProduct result(other.clone().as()); + result.prepend(1, to_expr_ptr(std::move(*this))); + *this = std::move(result); + } + + return *this; +} + void swap(ExprContainer &lhs, ExprContainer &rhs) { std::swap(lhs.expr_, rhs.expr_); } @@ -130,4 +148,11 @@ ExprContainer operator*(const Expr &lhs, const Expr &rhs) { return cont; } +ExprContainer operator^(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont ^= rhs; + + return cont; +} + } // namespace sequant diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 0a74a80772..7be1c90223 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -38,6 +38,7 @@ class ExprContainer { ExprContainer &operator+=(const Expr &expr); ExprContainer &operator-=(const Expr &expr); ExprContainer &operator*=(const Expr &expr); + ExprContainer &operator^=(const Expr &expr); friend void swap(ExprContainer &, ExprContainer &); @@ -52,6 +53,7 @@ class ExprContainer { ExprContainer operator+(const Expr &lhs, const Expr &rhs); ExprContainer operator-(const Expr &lhs, const Expr &rhs); ExprContainer operator*(const Expr &lhs, const Expr &rhs); +ExprContainer operator^(const Expr &lhs, const Expr &rhs); } // namespace sequant From 146f4d6901daac8e08c76328dfce657112de9534 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 10:34:49 +0200 Subject: [PATCH 30/53] Allow unique_ptr to be moved out of ExprContainer --- SeQuant/core/expressions/expr_container.cpp | 4 +++- SeQuant/core/expressions/expr_container.hpp | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index 2047ff5605..ac4ee33f9b 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -12,7 +12,7 @@ namespace sequant { ExprPtr to_expr_ptr(ExprContainer &&container) { - return std::shared_ptr(std::move(container.expr_)); + return std::move(container).take_expr(); } ExprContainer::ExprContainer(const ExprContainer &container) @@ -56,6 +56,8 @@ ExprContainer &ExprContainer::operator=(Expr &&expr) { ExprContainer ExprContainer::copy() const { return expr_->unique_copy(); } +std::unique_ptr ExprContainer::take_expr() && { return std::move(expr_); } + ExprContainer::operator const Expr &() const { return *expr_; } ExprContainer::operator Expr &() & { return *expr_; } diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 7be1c90223..0c2536b24f 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -25,6 +25,8 @@ class ExprContainer { ExprContainer copy() const; + std::unique_ptr take_expr() &&; + operator const Expr &() const; operator Expr &() &; operator Expr &&() &&; @@ -46,8 +48,6 @@ class ExprContainer { std::unique_ptr expr_; ExprContainer(std::unique_ptr expr); - - friend ExprPtr to_expr_ptr(ExprContainer &&container); }; ExprContainer operator+(const Expr &lhs, const Expr &rhs); From 5496e2456ff6252de4b632d51e499db77ae862a3 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 10:35:37 +0200 Subject: [PATCH 31/53] Add ExprContainer::operator* dereferencing to Expr && --- SeQuant/core/expressions/expr_container.cpp | 6 ++++-- SeQuant/core/expressions/expr_container.hpp | 5 +++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index ac4ee33f9b..831abb847b 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -63,9 +63,11 @@ ExprContainer::operator const Expr &() const { return *expr_; } ExprContainer::operator Expr &() & { return *expr_; } ExprContainer::operator Expr &&() && { return std::move(*expr_); } -const Expr &ExprContainer::operator*() const { return *expr_; } +const Expr &ExprContainer::operator*() const & { return *expr_; } -Expr &ExprContainer::operator*() { return *expr_; } +Expr &ExprContainer::operator*() & { return *expr_; } + +Expr &&ExprContainer::operator*() && { return std::move(*expr_); } const Expr *ExprContainer::operator->() const { return expr_.get(); } diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 0c2536b24f..47b4c3d4b8 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -31,8 +31,9 @@ class ExprContainer { operator Expr &() &; operator Expr &&() &&; - const Expr &operator*() const; - Expr &operator*(); + const Expr &operator*() const &; + Expr &operator*() &; + Expr &&operator*() &&; const Expr *operator->() const; Expr *operator->(); From 4a651e76ec8aa80d9c0c07cdc073f14d8c0c44bd Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 11:16:36 +0200 Subject: [PATCH 32/53] Make ExprContainer iterable --- SeQuant/core/expressions/expr_container.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 47b4c3d4b8..662954cef6 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -27,6 +27,13 @@ class ExprContainer { std::unique_ptr take_expr() &&; + auto begin() { return expr_->begin(); } + auto end() { return expr_->begin(); } + auto begin() const { return expr_->begin(); } + auto end() const { return expr_->begin(); } + auto cbegin() const { return expr_->begin(); } + auto cend() const { return expr_->begin(); } + operator const Expr &() const; operator Expr &() &; operator Expr &&() &&; From a1d0df0da2f1c436c1cbc47bcdeb2362acb911a2 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 11:17:45 +0200 Subject: [PATCH 33/53] Ensure specialized arithmetic operators only apply to ExprPtr ignoring any conversion --- SeQuant/core/expressions/expr_operators.hpp | 29 +++++++++++---------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/SeQuant/core/expressions/expr_operators.hpp b/SeQuant/core/expressions/expr_operators.hpp index 2f83df592a..4db7740876 100644 --- a/SeQuant/core/expressions/expr_operators.hpp +++ b/SeQuant/core/expressions/expr_operators.hpp @@ -19,83 +19,84 @@ namespace sequant { template requires(std::constructible_from) -ExprPtr operator+(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator+(const std::same_as auto &lhs, T &&rhs) { return lhs + ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator+(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator+(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) + rhs; } template requires(std::constructible_from) -ExprPtr operator-(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator-(const std::same_as auto &lhs, T &&rhs) { return lhs - ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator-(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator-(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) - rhs; } template requires(std::constructible_from) -ExprPtr operator*(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator*(const std::same_as auto &lhs, T &&rhs) { return lhs * ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator*(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator*(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) * rhs; } template requires(std::is_arithmetic_v) -ExprPtr operator/(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator/(const std::same_as auto &lhs, T &&rhs) { return lhs * ex(rational(1, std::forward(rhs))); } -inline ExprPtr operator/(const ExprPtr &lhs, const Constant &rhs) { +inline ExprPtr operator/(const std::same_as auto &lhs, + const Constant &rhs) { return lhs * ex(1.0 / rhs.value()); } template requires(std::constructible_from) -ExprPtr operator+(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator+(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) + rhs; } template requires(std::constructible_from) -ExprPtr operator+(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator+(const std::same_as auto &lhs, T &&rhs) { return lhs + ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator-(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator-(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) - rhs; } template requires(std::constructible_from) -ExprPtr operator-(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator-(const std::same_as auto &lhs, T &&rhs) { return lhs - ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator*(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator*(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) * rhs; } template requires(std::constructible_from) -ExprPtr operator*(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator*(const std::same_as auto &lhs, T &&rhs) { return lhs * ex(std::forward(rhs)); } From 30fed002d019f7f67027cb1a4f23694a3e5b9ac6 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 11:18:32 +0200 Subject: [PATCH 34/53] Disambiguate arithmetic ops involving ExprContainer --- SeQuant/core/expressions/expr_container.cpp | 32 +++++++++++++++++++++ SeQuant/core/expressions/expr_container.hpp | 11 +++++++ 2 files changed, 43 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index 831abb847b..3393314b31 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -138,6 +138,14 @@ ExprContainer operator+(const Expr &lhs, const Expr &rhs) { return cont; } +ExprContainer operator+(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) + rhs; +} + +ExprContainer operator+(const Expr &lhs, const ExprContainer &rhs) { + return lhs + static_cast(rhs); +} + ExprContainer operator-(const Expr &lhs, const Expr &rhs) { ExprContainer cont(lhs); cont -= rhs; @@ -145,6 +153,14 @@ ExprContainer operator-(const Expr &lhs, const Expr &rhs) { return cont; } +ExprContainer operator-(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) - rhs; +} + +ExprContainer operator-(const Expr &lhs, const ExprContainer &rhs) { + return lhs - static_cast(rhs); +} + ExprContainer operator*(const Expr &lhs, const Expr &rhs) { ExprContainer cont(lhs); cont *= rhs; @@ -152,6 +168,14 @@ ExprContainer operator*(const Expr &lhs, const Expr &rhs) { return cont; } +ExprContainer operator*(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) * rhs; +} + +ExprContainer operator*(const Expr &lhs, const ExprContainer &rhs) { + return lhs * static_cast(rhs); +} + ExprContainer operator^(const Expr &lhs, const Expr &rhs) { ExprContainer cont(lhs); cont ^= rhs; @@ -159,4 +183,12 @@ ExprContainer operator^(const Expr &lhs, const Expr &rhs) { return cont; } +ExprContainer operator^(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) ^ rhs; +} + +ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs) { + return lhs ^ static_cast(rhs); +} + } // namespace sequant diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 662954cef6..a49c026934 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -59,9 +59,20 @@ class ExprContainer { }; ExprContainer operator+(const Expr &lhs, const Expr &rhs); +ExprContainer operator+(const ExprContainer &lhs, const Expr &rhs); +ExprContainer operator+(const Expr &lhs, const ExprContainer &rhs); + ExprContainer operator-(const Expr &lhs, const Expr &rhs); +ExprContainer operator-(const ExprContainer &lhs, const Expr &rhs); +ExprContainer operator-(const Expr &lhs, const ExprContainer &rhs); + ExprContainer operator*(const Expr &lhs, const Expr &rhs); +ExprContainer operator*(const ExprContainer &lhs, const ExprContainer &rhs); +ExprContainer operator*(const Expr &lhs, const Expr &rhs); + ExprContainer operator^(const Expr &lhs, const Expr &rhs); +ExprContainer operator^(const ExprContainer &lhs, const Expr &rhs); +ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs); } // namespace sequant From 22feea537281c82749c2fdeffb9ac7bd80267255 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 11:19:25 +0200 Subject: [PATCH 35/53] Enable conversions from ExprContainer to ExprPtr --- SeQuant/core/expressions/expr_ptr.cpp | 6 ++++++ SeQuant/core/expressions/expr_ptr.hpp | 4 ++++ tests/unit/test_expr.cpp | 18 ++++++++++++++++++ 3 files changed, 28 insertions(+) diff --git a/SeQuant/core/expressions/expr_ptr.cpp b/SeQuant/core/expressions/expr_ptr.cpp index b1e21a0bb5..9276b92db7 100644 --- a/SeQuant/core/expressions/expr_ptr.cpp +++ b/SeQuant/core/expressions/expr_ptr.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -11,6 +12,11 @@ namespace sequant { +ExprPtr::ExprPtr(const ExprContainer &container) : ExprPtr(container.copy()) {} + +ExprPtr::ExprPtr(ExprContainer &&container) + : ExprPtr(std::move(container).take_expr()) {} + ExprPtr ExprPtr::clone() const & { if (!*this) return {}; return ExprPtr(as_shared_ptr()->clone()); diff --git a/SeQuant/core/expressions/expr_ptr.hpp b/SeQuant/core/expressions/expr_ptr.hpp index 39b9ebbe40..b453c16f6a 100644 --- a/SeQuant/core/expressions/expr_ptr.hpp +++ b/SeQuant/core/expressions/expr_ptr.hpp @@ -10,6 +10,8 @@ namespace sequant { +class ExprContainer; + /// @brief ExprPtr is a multiple-owner smart pointer to Expr /// It can be used mostly interchangeably with `std::shared_ptr`, but @@ -23,6 +25,8 @@ class ExprPtr : public std::shared_ptr { ExprPtr() = default; ExprPtr(const ExprPtr &) = default; ExprPtr(ExprPtr &&) = default; + explicit ExprPtr(const ExprContainer &container); + ExprPtr(ExprContainer &&container); template , Expr> || std::is_base_of_v>>> diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index 73d4e00189..4eef88c766 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -1237,6 +1237,24 @@ TEST_CASE("expr", "[elements]") { REQUIRE_THAT(res, EquivalentTo("(1 + A - B) * 3")); } + SECTION("Conversion to ExprPtr") { + ExprContainer cont = Constant(3); + + ExprPtr ptr = std::move(cont); + REQUIRE(ptr->is()); + REQUIRE(ptr->as().value() == 3); + + cont = Variable("A"); + // Copy-conversion-ctor is explicit + ptr = ExprPtr(cont); + REQUIRE(ptr->is()); + REQUIRE(ptr->as().label() == L"A"); + + // Ensure that ptr actually points to a copy + ptr->as().set_label(L"B"); + REQUIRE(cont->as().label() == L"A"); + REQUIRE(ptr->as().label() == L"B"); + } } SECTION("ResultExpr") { From 6c2a52418bea2aaf348de78e143c191c5a3cbd32 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 11:24:46 +0200 Subject: [PATCH 36/53] Fix deprecation warning in test_expr --- tests/unit/test_expr.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index 4eef88c766..b9f068ae87 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -911,10 +911,14 @@ TEST_CASE("expr", "[elements]") { REQUIRE(hash_value(ex(1)) == hash_value(ex(1))); - auto hasher = [](const std::shared_ptr &) -> unsigned int { + auto hasher1 = [](const std::shared_ptr &) -> unsigned int { return 0; }; - REQUIRE_NOTHROW(ex(1)->hash_value(hasher) == 0); + auto hasher2 = [](const Expr &) -> unsigned int { return 2; }; + SEQUANT_PRAGMA_IGNORE_DEPRECATED_BEGIN + REQUIRE_NOTHROW(ex(1)->hash_value(hasher1) == 0); + SEQUANT_PRAGMA_IGNORE_DEPRECATED_END + REQUIRE_NOTHROW(ex(1)->hash_value(hasher2) == 2); { // Power const auto c2 = ex(rational{1, 2}); From a0c1701d06a19a523db059d01970ceb9cb231a60 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 12:02:15 +0200 Subject: [PATCH 37/53] Add ExprContainer to fwd declarations --- SeQuant/core/expr_fwd.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SeQuant/core/expr_fwd.hpp b/SeQuant/core/expr_fwd.hpp index 1dbb36493c..05575986fe 100644 --- a/SeQuant/core/expr_fwd.hpp +++ b/SeQuant/core/expr_fwd.hpp @@ -12,6 +12,7 @@ namespace sequant { class Expr; class ResultExpr; class ExprPtr; +class ExprContainer; class Labeled; class Constant; From a65e1ba93293b2906a568ca9c10e55e34cc6016d Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 12:02:53 +0200 Subject: [PATCH 38/53] Enable ExprPtr -> ExprContainer conversion --- SeQuant/core/expressions/expr_container.cpp | 3 +++ SeQuant/core/expressions/expr_container.hpp | 4 ++++ tests/unit/test_expr.cpp | 13 +++++++++++++ 3 files changed, 20 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index 3393314b31..d32c508b3a 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -27,6 +27,9 @@ ExprContainer::ExprContainer(const Expr &expr) ExprContainer::ExprContainer(Expr &&expr) : ExprContainer(std::move(expr).unique_copy()) {} +ExprContainer::ExprContainer(const ExprPtr &ptr) + : ExprContainer(ptr->unique_copy()) {} + ExprContainer &ExprContainer::operator=(const ExprContainer &container) { // Copy-and-swap ExprContainer copy(container); diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index a49c026934..1fccbd998f 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -7,6 +7,8 @@ namespace sequant { +class ExprPtr; + class ExprContainer { public: explicit ExprContainer(const ExprContainer &container); @@ -15,6 +17,8 @@ class ExprContainer { explicit ExprContainer(const Expr &expr); ExprContainer(Expr &&expr); + explicit ExprContainer(const ExprPtr &ptr); + ExprContainer &operator=(const ExprContainer &container); ExprContainer &operator=(ExprContainer &&container) = default; diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index b9f068ae87..ecd573d024 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -1259,6 +1259,19 @@ TEST_CASE("expr", "[elements]") { REQUIRE(cont->as().label() == L"A"); REQUIRE(ptr->as().label() == L"B"); } + SECTION("Conversion from ExprPtr") { + // This conversion is always explicit as it always has to perform a copy. + // Even a moved-from ExpPtr might point to an object that is co-owned by + // another ExprPtr and thus "resource stealing" is not possible. + ExprPtr ptr = ex(2); + ExprContainer cont(ptr); + REQUIRE(cont->is()); + REQUIRE(cont->as().value() == 2); + + ptr->as() = Constant(3); + REQUIRE(ptr->as().value() == 3); + REQUIRE(cont->as().value() == 2); + } } SECTION("ResultExpr") { From 0b58236202ffaae986e2109286457649713fe9f0 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 13:17:12 +0200 Subject: [PATCH 39/53] Include supporting header --- SeQuant/core/expressions/expr_container.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index d32c508b3a..c1c6a36c0a 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include From 1abcf20766ea8340eae065b51ccedfac03773b36 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 14:52:35 +0200 Subject: [PATCH 40/53] Consistently use operator-> --- SeQuant/core/export/export.hpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/SeQuant/core/export/export.hpp b/SeQuant/core/export/export.hpp index 4b3627218a..86f627fecd 100644 --- a/SeQuant/core/export/export.hpp +++ b/SeQuant/core/export/export.hpp @@ -539,10 +539,10 @@ void track_usage(const EvalNode &node, PreprocessResult &result) { handle_variable(expr.as()); } else if (expr.is()) { const Power &power = expr.as(); - if (power.base().is()) { - handle_tensor(power.base().as()); - } else if (power.base().is()) { - handle_variable(power.base().as()); + if (power.base()->is()) { + handle_tensor(power.base()->as()); + } else if (power.base()->is()) { + handle_variable(power.base()->as()); } } } @@ -762,10 +762,10 @@ class PreprocessVisitor { handle_variable(node.left()->as_variable()); } else if (node.left()->is_power()) { const Power &power = node.left()->as_power(); - if (power.base().is()) { - handle_tensor(power.base().as()); - } else if (power.base().is()) { - handle_variable(power.base().as()); + if (power.base()->is()) { + handle_tensor(power.base()->as()); + } else if (power.base()->is()) { + handle_variable(power.base()->as()); } } @@ -775,10 +775,10 @@ class PreprocessVisitor { handle_variable(node.right()->as_variable()); } else if (node.right()->is_power()) { const Power &power = node.right()->as_power(); - if (power.base().is()) { - handle_tensor(power.base().as()); - } else if (power.base().is()) { - handle_variable(power.base().as()); + if (power.base()->is()) { + handle_tensor(power.base()->as()); + } else if (power.base()->is()) { + handle_variable(power.base()->as()); } } } From 9a3bd74bb62fa44cc269dd1e52d70f9396eb9723 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 14:53:08 +0200 Subject: [PATCH 41/53] Add expr_holder concept --- SeQuant/core/expressions/traits.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/SeQuant/core/expressions/traits.hpp b/SeQuant/core/expressions/traits.hpp index efeb818d82..76c793bf81 100644 --- a/SeQuant/core/expressions/traits.hpp +++ b/SeQuant/core/expressions/traits.hpp @@ -4,6 +4,9 @@ #include #include +#include +#include + namespace sequant { template @@ -36,6 +39,10 @@ constexpr bool is_a_power_v = meta::is_base_of_v; template constexpr bool is_power_v = meta::is_same_v; +template +concept expr_holder = std::same_as, ExprPtr> || + std::same_as, ExprContainer>; + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_TRAITS_HPP From 4bbeede5a84b883191f0010920047f88d3eb0e75 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 14:54:01 +0200 Subject: [PATCH 42/53] Make ExprContainer comparable with ExprPtr --- SeQuant/core/expressions/expr_container.cpp | 8 ++++++++ SeQuant/core/expressions/expr_container.hpp | 3 +++ 2 files changed, 11 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index c1c6a36c0a..aaf7f85097 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -195,4 +195,12 @@ ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs) { return lhs ^ static_cast(rhs); } +bool operator==(const ExprContainer &lhs, const ExprPtr &rhs) { + return *lhs == *rhs; +} + +bool operator==(const ExprPtr &lhs, const ExprContainer &rhs) { + return *lhs == *rhs; +} + } // namespace sequant diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 1fccbd998f..5a224b4fd0 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -78,6 +78,9 @@ ExprContainer operator^(const Expr &lhs, const Expr &rhs); ExprContainer operator^(const ExprContainer &lhs, const Expr &rhs); ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs); +bool operator==(const ExprContainer &lhs, const ExprPtr &rhs); +bool operator==(const ExprPtr &lhs, const ExprContainer &rhs); + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP From de2ef0abc246383dc4fc74660dbb53650756c08f Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 14:54:17 +0200 Subject: [PATCH 43/53] Add not-null assertion --- SeQuant/core/expressions/expr_container.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index aaf7f85097..f5ae2eb0b0 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -20,7 +20,9 @@ ExprContainer::ExprContainer(const ExprContainer &container) : ExprContainer(container->unique_copy()) {} ExprContainer::ExprContainer(std::unique_ptr expr) - : expr_(std::move(expr)) {} + : expr_(std::move(expr)) { + SEQUANT_ASSERT(expr_ != nullptr); +} ExprContainer::ExprContainer(const Expr &expr) : ExprContainer(expr.unique_copy()) {} From 6f96a42124b180b741f2a4b67781d891f73c6b5d Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 14:55:15 +0200 Subject: [PATCH 44/53] Switch Power to using ExprContainer --- SeQuant/core/export/itf.hpp | 2 +- .../core/export/julia_tensor_operations.hpp | 2 +- SeQuant/core/export/python_einsum.hpp | 2 +- SeQuant/core/export/text_generator.hpp | 2 +- SeQuant/core/export/utils.cpp | 5 +- SeQuant/core/export/utils.hpp | 6 +- SeQuant/core/expressions/power.cpp | 74 ++++++++++++------- SeQuant/core/expressions/power.hpp | 16 ++-- 8 files changed, 64 insertions(+), 45 deletions(-) diff --git a/SeQuant/core/export/itf.hpp b/SeQuant/core/export/itf.hpp index 3c1d0803fb..a7db26cb07 100644 --- a/SeQuant/core/export/itf.hpp +++ b/SeQuant/core/export/itf.hpp @@ -215,7 +215,7 @@ class ItfGenerator : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); // ITF can only express powers of Constants if (!base->is()) { throw Exception( diff --git a/SeQuant/core/export/julia_tensor_operations.hpp b/SeQuant/core/export/julia_tensor_operations.hpp index 70342fa2dc..8140ee2c66 100644 --- a/SeQuant/core/export/julia_tensor_operations.hpp +++ b/SeQuant/core/export/julia_tensor_operations.hpp @@ -130,7 +130,7 @@ class JuliaTensorOperationsGenerator : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); std::string base_str = to_julia_expr(*base, ctx); if (base->is() && base->as().conjugated()) { base_str = wrap_conj(std::move(base_str)); diff --git a/SeQuant/core/export/python_einsum.hpp b/SeQuant/core/export/python_einsum.hpp index ca1a6e8501..49cf3f84f3 100644 --- a/SeQuant/core/export/python_einsum.hpp +++ b/SeQuant/core/export/python_einsum.hpp @@ -198,7 +198,7 @@ class PythonEinsumGeneratorBase : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); std::string base_str = stringify_scalar(*base, ctx); if (base->is() && base->as().conjugated()) { base_str = wrap_conj(std::move(base_str)); diff --git a/SeQuant/core/export/text_generator.hpp b/SeQuant/core/export/text_generator.hpp index 0a478fc5dd..c99d9f3d64 100644 --- a/SeQuant/core/export/text_generator.hpp +++ b/SeQuant/core/export/text_generator.hpp @@ -121,7 +121,7 @@ class TextGenerator : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); std::string base_str = stringify(*base, ctx); if (base->is() && base->as().conjugated()) { base_str = wrap_conj(std::move(base_str)); diff --git a/SeQuant/core/export/utils.cpp b/SeQuant/core/export/utils.cpp index 182cc1977b..fd0841e6cc 100644 --- a/SeQuant/core/export/utils.cpp +++ b/SeQuant/core/export/utils.cpp @@ -4,8 +4,7 @@ #include -#include -#include +#include #include #include @@ -30,7 +29,7 @@ std::string format_power_exponent(const Power::exponent_type &exponent, return ss.str(); } -std::string format_power_base(const ExprPtr &base, std::string base_str) { +std::string format_power_base(const ExprContainer &base, std::string base_str) { if (base->is()) { const auto &v = base->as().value(); if (v.imag() == 0 && diff --git a/SeQuant/core/export/utils.hpp b/SeQuant/core/export/utils.hpp index 14798f42b5..103f217ac0 100644 --- a/SeQuant/core/export/utils.hpp +++ b/SeQuant/core/export/utils.hpp @@ -5,9 +5,7 @@ #ifndef SEQUANT_CORE_EXPORT_UTILS_HPP #define SEQUANT_CORE_EXPORT_UTILS_HPP -#include -#include -#include +#include #include @@ -27,7 +25,7 @@ std::string format_power_exponent(const Power::exponent_type &exponent, /// @param base_str @p base already rendered to a string by the caller /// @return @p base_str, wrapped in parens iff @p base is a Constant whose /// value is a non-integer or negative real -std::string format_power_base(const ExprPtr &base, std::string base_str); +std::string format_power_base(const ExprContainer &base, std::string base_str); } // namespace sequant::detail diff --git a/SeQuant/core/expressions/power.cpp b/SeQuant/core/expressions/power.cpp index 9e7cfc7450..ea8eeb259d 100644 --- a/SeQuant/core/expressions/power.cpp +++ b/SeQuant/core/expressions/power.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -9,19 +10,18 @@ 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(); +Power::Power(const ExprPtr& base, exponent_type exponent) + : Power(ExprContainer(base), std::move(exponent)) {} + +Power::Power(ExprContainer base, exponent_type exponent) + : base_{std::move(base)}, exponent_{std::move(exponent)} { + SEQUANT_ASSERT(base_->is() || base_->is()); // 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 ExprContainer& Power::base() const { return base_; } const Power::exponent_type& Power::exponent() const { return exponent_; } @@ -37,30 +37,38 @@ bool Power::is_zero() const { base_->as().is_zero(); } -void Power::flatten(ExprPtr& expr) { - if (!expr || !expr->is()) return; - const auto& pw = expr->as(); +template +void flatt_impl(E& expr) { + auto create_constant = [](const auto& val) { + if constexpr (std::same_as, ExprContainer>) { + return Constant(val); + } else { + return ex(val); + } + }; + + const auto& pw = expr->template as(); // b^1 = b and conjugate if needed - if (pw.exponent_ == 1) { - auto lifted = pw.base_->clone(); - if (pw.conjugated_) lifted->adjoint(); + if (pw.exponent() == 1) { + auto lifted = pw.base().copy(); + 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}); + if (pw.exponent() == 0) { + expr = create_constant(Constant::scalar_type{1}); return; } - if (!pw.base_->is()) return; + if (!pw.base()->template is()) return; using scalar_type = Constant::scalar_type; - const auto& base_val = pw.base_->as().value(); + const auto& base_val = pw.base()->template as().value(); // 1^k = 1 for any rational k. if (base_val == scalar_type{1}) { - expr = ex(scalar_type{1}); + expr = create_constant(scalar_type{1}); return; } @@ -78,11 +86,11 @@ void Power::flatten(ExprPtr& expr) { // initialize the base scalar_type base{0}; - auto exp_nr = numerator(pw.exponent_); // numerator of exponent + auto exp_nr = numerator(pw.exponent()); - if (denominator(pw.exponent_) == 1) { + if (denominator(pw.exponent()) == 1) { base = base_val; - } else if (denominator(pw.exponent_) == 2 && base_val.imag() == 0 && + } 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, @@ -112,8 +120,20 @@ void Power::flatten(ExprPtr& expr) { } if (negate) value = scalar_type{1} / value; - if (pw.conjugated_) value = conj(value); - expr = ex(std::move(value)); + if (pw.conjugated()) value = conj(value); + expr = create_constant(std::move(value)); +} + +void Power::flatten(ExprPtr& expr) { + if (!expr || !expr->is()) return; + + flatt_impl(expr); +} + +void Power::flatten(ExprContainer& expr) { + if (!expr->is()) return; + + flatt_impl(expr); } Expr::type_id_type Power::type_id() const { return get_type_id(); } @@ -155,9 +175,9 @@ Power& Power::operator*=(const Expr& that) { } std::unique_ptr Power::unique_copy() const { - auto cloned = std::make_unique(base_, exponent_); - if (conjugated_) cloned->as().conjugate(); - return cloned; + auto copy = std::make_unique(base_.copy(), exponent_); + if (conjugated_) copy->as().conjugate(); + return copy; } Expr::hash_type Power::memoizing_hash() const { diff --git a/SeQuant/core/expressions/power.hpp b/SeQuant/core/expressions/power.hpp index b9f594cf5b..86b9840156 100644 --- a/SeQuant/core/expressions/power.hpp +++ b/SeQuant/core/expressions/power.hpp @@ -3,7 +3,9 @@ #include #include +#include #include +#include #include #include @@ -27,25 +29,24 @@ 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); + Power(const ExprPtr& base, exponent_type exponent); + Power(ExprContainer base, exponent_type exponent); /// @overload constructs a `Variable` base from @p label template - requires std::constructible_from && - (!std::convertible_to) + requires(std::constructible_from && !expr_holder) Power(L&& label, exponent_type exponent) : Power(ex(std::forward(label)), std::move(exponent)) {} /// @overload constructs a `Constant` base from scalar @p value template - requires(!std::constructible_from && - !std::convertible_to && + requires(!std::constructible_from && !expr_holder && std::constructible_from) Power(V&& value, exponent_type exponent) : Power(ex(std::forward(value)), std::move(exponent)) {} /// @return the base expression - const ExprPtr& base() const; + const ExprContainer& base() const; /// @return the rational exponent const exponent_type& exponent() const; @@ -78,6 +79,7 @@ class Power : public Expr { /// 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); + static void flatten(ExprContainer& expr); type_id_type type_id() const override; @@ -100,7 +102,7 @@ class Power : public Expr { std::unique_ptr unique_copy() const override; private: - ExprPtr base_; + ExprContainer base_; exponent_type exponent_; bool conjugated_ = false; From 2d2a994de4b702649fb4b905eece23e608e5575a Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 15:02:06 +0200 Subject: [PATCH 45/53] Get rid of unnecessary heap allocations --- tests/unit/test_expr.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index ecd573d024..dd39cf1d80 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -294,9 +294,9 @@ TEST_CASE("expr", "[elements]") { // convenience ctors: REQUIRE(Power(L"x", 2) == Power(vx, rational{2})); REQUIRE(Power(L"x", rational{1, 2}) == Power(vx, rational{1, 2})); - REQUIRE(Power(2, 3) == Power(ex(2), rational{3})); + REQUIRE(Power(2, 3) == Power(Constant(2), rational{3})); REQUIRE(Power(rational{2, 3}, 2) == - Power(ex(rational{2, 3}), rational{2})); + Power(Constant(rational{2, 3}), rational{2})); if constexpr (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { // base must be a Constant or Variable; Power-of-Power is not allowed @@ -304,20 +304,20 @@ TEST_CASE("expr", "[elements]") { REQUIRE_THROWS(Power(inner, rational{2, 3})); // 0^n is defined only for n >= 0 - REQUIRE_THROWS(Power(ex(0), rational{-1})); + REQUIRE_THROWS(Power(Constant(0), rational{-1})); } } { // accessors Power p(c2, rational{1, 2}); - REQUIRE(p.base() == ex(rational{1, 2})); + REQUIRE(p.base() == Constant(rational{1, 2})); REQUIRE(p.exponent() == rational{1, 2}); // is_zero: base == 0, exponent > 0 - Power pz(ex(0), rational{2}); + Power pz(Constant(0), rational{2}); REQUIRE(pz.is_zero()); // 0^0 is not zero by our convention - Power pz2(ex(0), rational{0}); + Power pz2(Constant(0), rational{0}); REQUIRE(!pz2.is_zero()); REQUIRE(!p.is_zero()); } @@ -371,8 +371,8 @@ TEST_CASE("expr", "[elements]") { REQUIRE(pc_conj2.exponent() == rational{3, 2}); // 2^{1/2} * 2^{1/2} = 2 - Power pe(ex(2), rational{1, 2}); - Power pf(ex(2), rational{1, 2}); + Power pe(Constant(2), rational{1, 2}); + Power pf(Constant(2), rational{1, 2}); pe *= pf; REQUIRE(pe.exponent() == rational{1}); REQUIRE(to_latex(pe) == Constant(2).to_latex()); From 7955bbd81b54a081464f18909b77c28971436ce2 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 15:22:28 +0200 Subject: [PATCH 46/53] Use C++ templating facilities --- SeQuant/core/expressions/sum.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index d050a933d8..7e36c924b1 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -16,8 +16,10 @@ #include #include +#include #include #include +#include #include namespace sequant { @@ -53,15 +55,14 @@ class Sum : public Expr { /// construct a Sum out of a range of summands /// @param rng a range - template - requires(meta::is_range_v> && - !meta::is_same_v, ExprPtrList>) + template + requires(!std::same_as, ExprPtrList>) explicit Sum(Range &&rng) { // N.B. use append to flatten out Sum summands - constexpr auto rng_is_expr = - meta::is_base_of_v>; + constexpr auto rng_is_expr = is_an_expr_v>; constexpr auto rng_is_exprptr = - meta::is_same_v>; + std::same_as>; + if constexpr (rng_is_expr || rng_is_exprptr) { ExprPtr rng_as_exprptr; if constexpr (rng_is_expr) { From 726ab37ee7546f74bca060dcf468c079d7ed8765 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 15:22:50 +0200 Subject: [PATCH 47/53] Change Sum::filter to use predicate on cosnt Expr & --- SeQuant/core/expressions/sum.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index 7e36c924b1..b70fe70686 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -109,13 +110,13 @@ class Sum : public Expr { /// offset ExprPtr take_n(size_t offset, size_t count) const; - /// @tparam Filter a boolean predicate type, such `Filter(const ExprPtr&)` - /// evaluates to true - /// @param f an object of Filter type - /// Selects elements {`e`} for which `f(e)` is true - template + /// @param f Boolean predicate + /// @returns A sum containing only the summands for which f was true. + template Filter> ExprPtr filter(Filter &&f) const { - return ex(summands_ | ranges::views::filter(f)); + return ex(summands_ | + ranges::views::transform([](const auto &e) { return *e; }) | + ranges::views::filter(f)); } /// @return true if the number of factors is zero From 96d4c94e3f9a117377432253e1153f0a4974aa74 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Tue, 18 Aug 2026 15:27:25 +0200 Subject: [PATCH 48/53] Use Catch2 sections instead of comments --- tests/unit/test_expr.cpp | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index dd39cf1d80..9fbf6424e5 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -287,7 +287,7 @@ TEST_CASE("expr", "[elements]") { const auto c2 = ex(rational{1, 2}); const auto vx = ex(L"x"); - { // constructors + SECTION("constructors") { REQUIRE_NOTHROW(Power(c2, rational{1, 2})); REQUIRE_NOTHROW(Power(vx, rational{3, 1})); @@ -308,7 +308,7 @@ TEST_CASE("expr", "[elements]") { } } - { // accessors + SECTION("accessors") { Power p(c2, rational{1, 2}); REQUIRE(p.base() == Constant(rational{1, 2})); REQUIRE(p.exponent() == rational{1, 2}); @@ -322,7 +322,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(!p.is_zero()); } - { // comparison + SECTION("comparison") { Power p1(c2, rational{1, 2}); Power p2(c2, rational{1, 2}); REQUIRE(p1 == p2); @@ -341,7 +341,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(!(plt_a < plt_c)); } - { // operator*= + SECTION("operator*=") { // b^e1 *= b^e2 -> b^(e1+e2) Power pa(vx, rational{1, 2}); Power pb(vx, rational{1, 3}); @@ -378,7 +378,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(to_latex(pe) == Constant(2).to_latex()); } - { // Power should NOT be absorbed into Product::scalar_ + SECTION("Don't absorb into Product::scalar") { auto p = ex(vx, rational{1, 2}); auto prod = ex(Product{}); prod->as().append(1, p, Product::Flatten::Yes); @@ -423,12 +423,12 @@ TEST_CASE("expr", "[elements]") { REQUIRE_NOTHROW(e->adjoint()); REQUIRE_NOTHROW(adjoint(e)); // check free-function adjoint } - { // Constant + SECTION("Constant") { const auto e = std::make_shared(Constant::scalar_type{1, 2}); REQUIRE_NOTHROW(e->adjoint()); REQUIRE(e->value() == Constant::scalar_type{1, -2}); } - { // Variable + SECTION("Variabkle") { const auto e = std::make_shared(L"q"); REQUIRE(e->conjugated() == false); REQUIRE_NOTHROW(e->adjoint()); @@ -437,7 +437,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE_NOTHROW(e->adjoint()); REQUIRE(e->conjugated() == false); } - { // Product + SECTION("Product") { // Product const auto e = std::make_shared(); e->append(Constant::scalar_type{2, -1}, ex()); e->append(1, ex(-2)); @@ -446,7 +446,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->factors()[0]->as().v == 2); REQUIRE(e->factors()[1]->as().v == -1); } - { // CProduct + SECTION("CProduct") { const auto e = std::make_shared(); e->append(Constant::scalar_type{2, -1}, ex()); e->append(1, ex(-2)); @@ -455,7 +455,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->factors()[0]->as().v == -1); REQUIRE(e->factors()[1]->as().v == 2); } - { // NCProduct + SECTION("NCProduct") { const auto e = std::make_shared(); e->append(Constant::scalar_type{2, -1}, ex()); e->append(1, ex(-2)); @@ -464,7 +464,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->factors()[0]->as().v == 2); REQUIRE(e->factors()[1]->as().v == -1); } - { // Sum + SECTION("Sum") { const auto e = std::make_shared(); e->append(ex()); e->append(ex(-2)); @@ -472,7 +472,8 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->summands()[0]->as().v == -1); REQUIRE(e->summands()[1]->as().v == 2); } - { // Power: adjoint flips the conjugation flag; base/exponent unchanged + SECTION("Power") { + // adjoint flips the conjugation flag; base/exponent unchanged Power pv(ex(L"z"), rational{1, 2}); REQUIRE(!pv.conjugated()); pv.adjoint(); From aec0a07ced9b263bfc87d9b5a5fd51cbcac76c68 Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Thu, 20 Aug 2026 12:46:35 +0200 Subject: [PATCH 49/53] Add Expr::proportional_to(const Expr &, const Expr &) overload --- SeQuant/core/expressions/expr.cpp | 20 ++++++++++++-------- SeQuant/core/expressions/expr.hpp | 1 + 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index 34b0c4ed54..59fdd2c2d3 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -80,17 +80,16 @@ std::weak_ptr Expr::weak_from_this() const { return std::enable_shared_from_this::weak_from_this(); } -bool proportional_to::operator()(const ExprPtr &expr1, - const ExprPtr &expr2) const { - if (expr1->type_id() != - expr2->type_id()) { // if expr1 is a Product with single factor == expr2, - // or vice versa +bool proportional_to::operator()(const Expr &expr1, const Expr &expr2) const { + if (expr1.type_id() != expr2.type_id()) { + // if expr1 is a Product with single factor == expr2, + // or vice versa if (expr1.is()) { return expr1.as().factors().size() == 1 && - expr1.as().factors().front() == expr2; + *expr1.as().factors().front() == expr2; } else if (expr2.is()) { return expr2.as().factors().size() == 1 && - expr2.as().factors().front() == expr1; + *expr2.as().factors().front() == expr1; } else return false; } @@ -101,10 +100,15 @@ bool proportional_to::operator()(const ExprPtr &expr1, return true; } if (expr1.is()) { - return expr1->hash_value() == expr2->hash_value() && + return expr1.hash_value() == expr2.hash_value() && expr1.as().factors() == expr2.as().factors(); } return expr1 == expr2; } +bool proportional_to::operator()(const ExprPtr &expr1, + const ExprPtr &expr2) const { + return (*this)(*expr1, *expr2); +} + } // namespace sequant diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 7cda1df709..ebbc21072a 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -502,6 +502,7 @@ struct proportional_to { /// @param[in] expr1 /// @param[in] expr2 /// @return true if @p expr1 is proportional to @p expr2 + bool operator()(const Expr &expr1, const Expr &expr2) const; bool operator()(const ExprPtr &expr1, const ExprPtr &expr2) const; }; From 07db221c6c664bb396ac7679e4a641279346a51e Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Thu, 20 Aug 2026 12:52:48 +0200 Subject: [PATCH 50/53] Move begin/end impl to cpp file (and fix it) --- SeQuant/core/expressions/expr_container.cpp | 25 +++++++++++++++++++++ SeQuant/core/expressions/expr_container.hpp | 13 ++++++----- 2 files changed, 32 insertions(+), 6 deletions(-) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index f5ae2eb0b0..28ca6d6916 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -64,6 +65,30 @@ ExprContainer ExprContainer::copy() const { return expr_->unique_copy(); } std::unique_ptr ExprContainer::take_expr() && { return std::move(expr_); } +ExprIterator ExprContainer::begin() { + SEQUANT_ASSERT(expr_); + return expr_->begin(); +} + +ExprIterator ExprContainer::end() { + SEQUANT_ASSERT(expr_); + return expr_->end(); +} + +ConstExprIterator ExprContainer::begin() const { + SEQUANT_ASSERT(expr_); + return std::as_const(*expr_).begin(); +} + +ConstExprIterator ExprContainer::end() const { + SEQUANT_ASSERT(expr_); + return std::as_const(*expr_).end(); +} + +ConstExprIterator ExprContainer::cbegin() const { return begin(); } + +ConstExprIterator ExprContainer::cend() const { return end(); } + ExprContainer::operator const Expr &() const { return *expr_; } ExprContainer::operator Expr &() & { return *expr_; } diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 5a224b4fd0..83a8618ccc 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -2,6 +2,7 @@ #define SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP #include +#include #include @@ -31,12 +32,12 @@ class ExprContainer { std::unique_ptr take_expr() &&; - auto begin() { return expr_->begin(); } - auto end() { return expr_->begin(); } - auto begin() const { return expr_->begin(); } - auto end() const { return expr_->begin(); } - auto cbegin() const { return expr_->begin(); } - auto cend() const { return expr_->begin(); } + ExprIterator begin(); + ExprIterator end(); + ConstExprIterator begin() const; + ConstExprIterator end() const; + ConstExprIterator cbegin() const; + ConstExprIterator cend() const; operator const Expr &() const; operator Expr &() &; From 52aede3e509c2a133ccc58d95a019ac9df32f3ef Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Thu, 20 Aug 2026 14:14:16 +0200 Subject: [PATCH 51/53] Use fwd decl of Expr --- SeQuant/core/expressions/expr_container.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 83a8618ccc..53bce22689 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -1,13 +1,13 @@ #ifndef SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP #define SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP -#include #include #include namespace sequant { +class Expr; class ExprPtr; class ExprContainer { From d5ecd7fe20b5519daf8b0d99f0c0b3fd0ec13c0a Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Thu, 20 Aug 2026 14:15:13 +0200 Subject: [PATCH 52/53] Define ExprContainerList --- SeQuant/core/expressions/expr_container.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 53bce22689..82808ec5fd 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -3,6 +3,7 @@ #include +#include #include namespace sequant { @@ -63,6 +64,8 @@ class ExprContainer { ExprContainer(std::unique_ptr expr); }; +using ExprContainerList = std::initializer_list; + ExprContainer operator+(const Expr &lhs, const Expr &rhs); ExprContainer operator+(const ExprContainer &lhs, const Expr &rhs); ExprContainer operator+(const Expr &lhs, const ExprContainer &rhs); From 09828aaf1c933ac95a2c79c19d311f614768cadc Mon Sep 17 00:00:00 2001 From: Robert Adam Date: Thu, 20 Aug 2026 14:15:54 +0200 Subject: [PATCH 53/53] Add free-standing adjoint() overload for ExprContainer --- SeQuant/core/expressions/expr_container.cpp | 6 ++++++ SeQuant/core/expressions/expr_container.hpp | 2 ++ 2 files changed, 8 insertions(+) diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp index 28ca6d6916..7aaefc3e8c 100644 --- a/SeQuant/core/expressions/expr_container.cpp +++ b/SeQuant/core/expressions/expr_container.cpp @@ -230,4 +230,10 @@ bool operator==(const ExprPtr &lhs, const ExprContainer &rhs) { return *lhs == *rhs; } +ExprContainer adjoint(const ExprContainer &cont) { + ExprContainer copy = cont.copy(); + copy->adjoint(); + return copy; +} + } // namespace sequant diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp index 82808ec5fd..23d97c8424 100644 --- a/SeQuant/core/expressions/expr_container.hpp +++ b/SeQuant/core/expressions/expr_container.hpp @@ -85,6 +85,8 @@ ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs); bool operator==(const ExprContainer &lhs, const ExprPtr &rhs); bool operator==(const ExprPtr &lhs, const ExprContainer &rhs); +ExprContainer adjoint(const ExprContainer &expr); + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP