Polynomial expressions
This page describes how libcommute represent polynomial expressions built out of noncommuting quantummechanical operators and what features such expressions support.
Class definitions

template<typename ScalarType, typename ...IndexTypes>
class expression Defined in <libcommute/expression/expression.hpp>
A polynomial expression, the main building block of libcommute’s DSL.
Expressions are finite ordered sums of monomials \(M^{(n)}_{i_1 i_2 \ldots i_n}\) accompanied by their respective coefficients \(C_{i_1 i_2 \ldots i_n}\),
\[E = C M^{(0)} + \sum_i C_i M^{(1)}_i + \sum_{ij} C_{ij} M^{(2)}_{ij} + \ldots\]The monomials are in turn canonically ordered products of algebra generators (operators \(c^\dagger\)/\(c\), \(a^\dagger\)/\(a\), etc).
ScalarType
is the type of coefficients \(C_{i_1 i_2 \ldots i_n}\). The most common choices ofScalarType
are double for expressions with real coefficients and std::complex<double> for the complex expressions. There are two convenience type aliases in the libcommute::static_indices namespace for these particular scalar types,static_indices::expr_real
andstatic_indices::expr_complex
. Use of custom scalar types – subject to some requirements – is also allowed. See Section Custom scalar types for more details.Indices \(i\), \(j\), etc in the definition above are compound indices with a fixed number of components. Types of the components are given by the template parameter pack
IndexTypes
. All index types must be lesscomparable and form strictly ordered sets.// Polynomial expressions with real coefficients. Monomials in 'expr1' // and 'expr2' are products of operators with one integer index. libcommute::expression<double, int> expr1; libcommute::static_indices::expr_real<int> expr2; // Polynomial expressions with complex coefficients. Monomials in 'expr3' // and 'expr4' are products of operators with one integer // and one string index. libcommute::expression<std::complex<double>, int, std::string> expr3; libcommute::static_indices::expr_complex<int, std::string> expr4;
Member type aliases

using scalar_type = ScalarType

using index_types = std::tuple<IndexTypes...>

using monomial_t = monomial<IndexTypes...>

using monomials_map_t = std::map<monomial_t, ScalarType>
Constructors

expression() = default
Construct a trivial expression, i.e. an expression containing no monomials.

template<typename S>
expression(expression<S, IndexTypes...> const &x)
Construct a copy of an expression x with a different scalar type S by converting its coefficients to
ScalarType
.Construct a constant expression \(E = x M^{(0)}\).

template<typename S>
explicit expression(S const &x, monomial_t const &monomial)
Construct an expression made of exactly one monomial, \(E = x M^{(n)}_{i_1 i_2 \ldots i_n}\).
The constructors listed here have limited functionality. One is supposed to use the factory functions to build expressions most of the time.
Copy/moveconstructors and assignments

expression(expression const&) = default

expression(expression&&) noexcept = default

expression &operator=(expression const&) = default

expression &operator=(expression&&) noexcept = default

template<typename S>
expression &operator=(expression<S, IndexTypes...> const &x)
Arithmetic operations
Two expression objects can be added, subtracted and multiplied as long as their
IndexTypes
agree and the corresponding arithmetic operation is defined for their scalar types. For example, it is possible to mix real and complex expressions as operands of+
,
and*
. The unary minus is defined for an expression type iff it is defined for expression’s scalar type. Another possibility is to mix expression objects with objects of other arbitrary types in binary operations. Object of the nonexpression types are treated as constants, i.e. contributions to \(C M^{(0)}\).The following table shows how result types of arithmetic operations are calculated.
Arithmetic operation
Result type
expression<S1, IT...>{} + expression<S2, IT...>{}
expression<decltype(S1{} + S2{}), IT...>
expression<S1, IT...>{}  expression<S2, IT...>{}
expression<decltype(S1{}  S2{}), IT...>
expression<S1, IT...>{} * expression<S2, IT...>{}
expression<decltype(S1{} * S2{}), IT...>
expression<S, IT...>{}
expression<decltype(S{}), IT...>
expression<S1, IT...>{} + S2{}
expression<decltype(S1{} + S2{}), IT...>
expression<S1, IT...>{}  S2{}
expression<decltype(S1{}  S2{}), IT...>
expression<S1, IT...>{} * S2{}
expression<decltype(S1{} * S2{}), IT...>
S1{} + expression<S2, IT...>{}
expression<decltype(S1{} + S2{}), IT...>
S1{}  expression<S2, IT...>{}
expression<decltype(S1{}  S2{}), IT...>
S1{} * expression<S2, IT...>{}
expression<decltype(S1{} * S2{}), IT...>
Compound assignments
+=
,=
,*=
are available under the same scalar type compatibility conditions between LHS and RHS. If the RHS is of a nonexpression type S, libcommute will attempt to select the optimized compound operator ScalarType::operator+=(S const &x) first (similarly for=
,*=
). If it fails, ScalarType::operator+(S const &x) and the regular assignment will be called instead.Iteration interface and transformations

class const_iterator
Constant bidirectional iterator over monomialcoefficient pairs \((M,C)\) in a polynomial expression. Given an iterator it, it>monomial returns a constant reference to the
monomial
object \(M\), and it>coeff is a constant reference to the respective coefficient \(C\).

const_iterator begin() const noexcept

const_iterator cbegin() const noexcept
Constant iterator to the first monomialcoefficient pair.

const_iterator end() const noexcept

const_iterator cend() const noexcept
Constant pasttheend iterator.

template<typename F, typename NewScalarType>
friend expression<NewScalarType, IndexTypes...> transform(expression const &expr, F &&f) Apply functor
f
to each monomialcoefficient pair inexpr
. Return a new expression obtained by replacing coefficients inexpr
with respective values returned byf
. The expected signature off
isNewScalarType f(monomial_t const&, ScalarType const&)
The transformed scalar type
NewScalarType
is automatically deduced fromf
’s return type. Whenf
returns a zero for a certain monomial, that monomial is excluded from the resulting expression.
Other methods and friend functions

std::map<monomial_t, ScalarType> const &get_monomials() const
Direct readonly access to the list of monomials.

std::size_t size() const
Number of monomials in this polynomial expression.

void clear()
Set expression to zero by removing all monomials.

friend bool operator==(expression const &e1, expression const &e2)
Check if
e1
ande2
contain identical lists of monomials.

friend bool operator!=(expression const &e1, expression const &e2)
Check if
e1
ande2
contain different lists of monomials.

friend expression conj(expression const &expr)
Return Hermitian conjugate of
expr
.

friend std::ostream &operator<<(std::ostream &os, expression const &expr)
Output stream insertion operator.

using scalar_type = ScalarType

template<typename ...IndexTypes>
using static_indices::expr_real = expression<double, IndexTypes...>; Declared in <libcommute/expression/expression.hpp>
Shorthand type for expressions with real coefficients and statically typed indices.

template<typename ...IndexTypes>
using static_indices::expr_complex = expression<std::complex<double>, IndexTypes...>; Declared in <libcommute/expression/expression.hpp>
Shorthand type for expressions with complex coefficients and statically typed indices.
Custom scalar types
Choosing double or std::complex<double> as the scalar type of
expressions covers the vast majority of practically important cases.
Nonetheless, sometimes it may be desirable to go beyond and pass a
userdefined type as ScalarType
.
For instance, one may want to use types from an arbitraryprecision arithmetic
library to represent expansion coefficients of a quantummechanical operator.
Another potential use  making coefficients depend on a parameter, such as time
or an external field. The polynomial
class from the Boost Math Toolkit or
a similar type would allow to represent the functional dependence on
the parameter while also implementing addition, subtraction and multiplication
operations.
Mathematically speaking, instances of
ScalarType
are assumed to form
a ring without a multiplicative identity
(a.k.a. rng). More
specifically, the set of ScalarType
instances must be
An abelian group under addition (with binary
operator+
andoperator
). In particular, there must be a welldefined zero element.A semigroup under multiplication (
operator*
).Multiplication must be distributive with respect to addition.
Let us say we have a type S with the required algebraic properties. Before using it as a scalar type, we must define a specialization of structure scalar_traits in the namespace libcommute to teach libcommute how to deal with the new type.
namespace libcommute {
template<> struct scalar_traits<S> {
// Test whether 's' is the zero element
static bool is_zero(S const& s) { ... }
// Make a constant of type 'S' from a double value 'x'
static S make_const(double x) { ... }
// OPTIONAL: Complex conjugate of 's'
static S conj(S const& s) { ... }
};
}
The static member scalar_traits<S>::conj() is optional and will only be
called by the Hermitian conjugation function
conj()
.
Note
For the builtin floatingpoint types, the zerovalue test method is implemented as
static bool is_zero(S const& s) {
return std::abs(s) < 100 * std::numeric_limits<S>::epsilon();
}
One can adjust the test and change the constant 100 to something else by defining a special macro LIBCOMMUTE_FLOATING_POINT_TOL_EPS.
[C++17] Dynamically typed index sequences
On the most basic level, index sequences of all generators found in a single
expression must agree in types with the IndexTypes
template parameter pack. In many
situations, however, it is more natural to have generators with different
numbers/types of indices mixed in one expression. This is where the dynamically
typed index sequences step in. They are instantiations of the
dyn_indices_generic class template defined in a special nested namespace
libcommute::dynamic_indices.

template<typename ...IndexTypes>
class dynamic_indices::dyn_indices_generic A wrapper around std::vector<std::variant<IndexTypes...>>.

using indices_t = std::vector<std::variant<IndexTypes...>>
Underlying dynamically typed sequence of indices.

dyn_indices_generic() = default
Construct an empty index sequence.

dyn_indices_generic(indices_t &&indices)

dyn_indices_generic(indices_t const &indices)
Construct an index sequence from a vector of indices.

template<typename ...Args>
dyn_indices_generic(Args&&... args) Construct an index sequence from an argument pack of indices.

dyn_indices_generic(dyn_indices_generic const&) = default

dyn_indices_generic(dyn_indices_generic&&) noexcept = default

dyn_indices_generic &operator=(dyn_indices_generic const&) = default

dyn_indices_generic &operator=(dyn_indices_generic&&) noexcept = default
Copy/moveconstructors and assignments

std::size_t size() const
Number of indices in the sequence.

friend bool operator==(dyn_indices_generic const &ind1, dyn_indices_generic const &ind2)

friend bool operator!=(dyn_indices_generic const &ind1, dyn_indices_generic const &ind2)

friend bool operator<(dyn_indices_generic const &ind1, dyn_indices_generic const &ind2)

friend bool operator>(dyn_indices_generic const &ind1, dyn_indices_generic const &ind2)
Compare two dynamic index sequences ind1 and ind2. These operators compare sequences’ lengths first, and in the case of equal lengths call the corresponding methods of std::vector.

friend std::ostream &operator<<(std::ostream &os, dyn_indices_generic const &ind)
Output stream insertion operator.

using indices_t = std::vector<std::variant<IndexTypes...>>
#include <libcommute/expression/dyn_indices.hpp>
// A userdefined index type
enum spin {up, down};
// Our own dynamically typed index sequence (mind the namespace!)
using my_dyn_indices =
libcommute::dynamic_indices::dyn_indices_generic<int, std::string, spin>;
// An expression with dynamically typed indices
libcommute::expression<double, my_dyn_indices> dyn_expr;
dyn_expr is an expression with dynamically typed indices. It can contain generators that effectively carry a variable number of indices, each of type int, std::string or of the userdefined enumeration type spin.
There is also a special set of factory functions defined in namespaces nested under libcommute::dynamic_indices. Those return commonly used QM operators with the dynamically typed indices. Some related type aliases are declared in the same namespace for the sake of convenience.

using dynamic_indices::dyn_indices = dyn_indices_generic<int, std::string>
Declared in <libcommute/expression/dyn_indices.hpp>
Dynamic mixture of integer and string indices.

using dynamic_indices::expr_real = expression<double, dyn_indices>

using dynamic_indices::expr_complex = expression<std::complex<double>, dyn_indices>
Declared in <libcommute/expression/expression.hpp>
Real/complex expression shorthand types.
Iteration interface and transformations
Both expressions and their constituent monomials can be easily iterated over
using STLcompatible iterators
or
rangebased for
loops.
This allows for writing complex expression analysis algorithms.
using namespace libcommute;
// We are going to analyse the structure of this expression
expression<double, int> E;
//
// Fill expression 'E' ...
//
// Iterate over all monomialcoefficient pairs in 'E'
for(auto const& mc : E) {
std::cout << "Coefficient: " << mc.coeff << "\n";
std::cout << "Monomial: ";
// Iterate over algebra generators in current monomial
for(auto const& g : mc.monomial) {
if(is_fermion(g)) { // Print information about fermionic operators
auto const& f = dynamic_cast<generator_fermion<int> const&>(g);
std::cout << (f.dagger() ? " c^+(" : " c(");
std::cout << std::get<0>(g.indices());
std::cout << ")";
}
if(is_fermion(g)) { // Print information about bosonic operators
auto const& a = dynamic_cast<generator_boson<int> const&>(g);
std::cout << (a.dagger() ? " a^+(" : " a(");
std::cout << std::get<0>(g.indices());
std::cout << ")";
}
if(is_spin(g)) { // Print information about spin operators
auto const& s = dynamic_cast<generator_spin<int> const&>(g);
switch(s.component()) {
case plus:
std::cout << "S_+("; break;
case minus:
std::cout << "S_("; break;
case z:
std::cout << "S_z("; break;
}
std::cout << std::get<0>(g.indices());
std::cout << ")";
}
}
std::cout << '\n';
}
Note
Only the constant iterators are implemented by expression
and
monomial
.
Another common task is building a new expression out of an existing one by
changing some of monomial coefficients. Removing some monomials (for instance,
those corresponding to particle interaction terms) is a special case that
amounts to setting coefficients of the unwanted monomials to zero. In the
following example we show how to use function
libcommute::transform()
to change
Hamiltonian of a finite atomic chain
into the SuSchriefferHeeger (SSH) model, where the hopping constant is taken to be different on odd and even chain links.
using namespace libcommute;
// Hamiltonian of the atomic chain
expression<double, int> H;
using static_indices::c_dag;
using static_indices::c;
// Length of the chain
int const N = 10;
// Hopping matrix element
double const v = 1.0;
// Add hopping terms
for(int a = 0; a < N  1; ++a) {
H += v * (c_dag(a) * c(a + 1) + c_dag(a + 1) * c(a));
}
std::cout << "H = " << H << '\n';
// Construct the SuSchriefferHeeger (SSH) model by changing hopping
// constants on all even links of the chain.
// Hopping matrix element on the even links
double const w = 0.5;
// Transformation operation
auto update_hopping_element = [v, w](decltype(H)::monomial_t const& m,
double coeff) > double {
// The monomial 'm' we are expecting here is a product c^\dagger(a1) c(a2)
// Site index of c^\dagger
int a1 = std::get<0>(m[0].indices());
// Site index of c
int a2 = std::get<0>(m[1].indices());
// Index of link connecting a1 and a2
int link = std::min(a1, a2);
// Return the updated matrix element
return (link % 2 == 1 ? w : v);
};
auto H_SSH = transform(H, update_hopping_element);
std::cout << "H_SSH = " << H_SSH << '\n';
\(\pm H.c.\) notation
It is common to use the \(\pm H.c.\) (“plus/minus Hermitian conjugate”)
notation to abbreviate analytical expressions when writing down Hamiltonians
and other Hermitian operators. The same trick works with libcommute’s
expressions, thanks to the hc
singleton object.
Here is a simple usage example that constructs two operators
#include <libcommute/expression/hc.hpp>
auto H_A = c_dag(1) * c(2) + hc;
std::complex<double> I(0,1);
auto H_B = I * c_dag(1) * c(2)  hc;