Skip to main content

Algebra

Language revisions

C++20
Concepts allow the algebraic hierarchy to be stated as a refinement chain rather than as documentation.
C++23
Extended floating-point types add carriers with distinct rounding behaviour to the same concepts.P1467R9
C++26
Contracts permit some algebraic preconditions to be checked at runtime where they cannot be checked statically.P2900

Definition​

An algebraic structure is a set with operations satisfying laws. The structures used in this reference form a chain, each adding operations or laws to the previous.

Carrier, operations, laws​

Three things are being stated separately, and conflating them is the source of most confusion when the structure is encoded in a type system.

The carrier is not the mathematical set​

The carrier is the set of values actually available. For double it is a finite subset of R∪{±∞,NaN}\mathbb{R} \cup \{\pm\infty, \mathrm{NaN}\} with about 2642^{64} elements, not R\mathbb{R}. For int it is an interval of Z\mathbb{Z}. For a normalised rational type it is the set of pairs whose components fit in the underlying integer type, which is not closed under addition.

Every statement in this entry about a machine type is therefore a statement about a finite carrier, and the interesting question in each case is which laws survive the restriction.

Operations must be total​

The operations are total functions from the carrier to itself. Totality is a real requirement: division on a field is not an operation of the structure, because it is undefined at zero, which is why fields are defined with a partial inverse on S∖{0}S \setminus \{0\} rather than with a total division.

The same distinction appears in interface design. A Field concept requiring operator/ promises a total function that the mathematics does not provide, so the precondition y≠0y \neq 0 travels with every use and belongs in the documentation of the concept rather than in the documentation of each algorithm.

Closure and partiality​

Closure is part of the definition of an operation: x,y∈Sx, y \in S implies x⊕y∈Sx \oplus y \in S. Machine arithmetic satisfies closure only by adding elements to the carrier that the mathematics does not have, or by leaving the carrier entirely.

TypeOperationClosedStructure obtained
unsigned++, ×\timesyesThe ring Z/2nZ\mathbb{Z}/2^{n}\mathbb{Z}, exactly
int++, ×\timesnoUndefined behaviour on overflow; no structure
double++, ×\timesyes, with ∞\infty and NaNCommutative magma with identity
rational++, ×\timesnoField until the components overflow
modular<P>++, ×\timesyesRing, and a field when PP is prime

Unsigned integers are the instructive row: they are frequently described as a defect of the language, and they are the one arithmetic type in C++ that models a textbook structure exactly.

Equality and the substitution property​

A law is an equation, so a structure presupposes an equality relation. The requirement is that equality be an equivalence relation and that every operation respect it:

a=b  ⟹  f(a)=f(b)for every operation f.a = b \;\Longrightarrow\; f(a) = f(b) \quad\text{for every operation } f.

This is the substitution property, and it is exactly the regularity requirement of Types as Abstractions. It fails for floating point at one point only, and instructively: NaN != NaN, so equality is not reflexive on the full carrier, which is why every ordering-based algorithm has undefined behaviour in the presence of NaN.

The chain of structures​

Semigroups and monoids​

A semigroup is the weakest structure supporting a fold. A monoid is the weakest supporting a fold over a possibly empty sequence, and the weakest for which a parallel decomposition needs no special case.

The identity is not a convenience. In a decomposition into pp chunks, some chunk may be empty whenever n<pn < p or the partition is uneven, and without an identity that case has to be excluded by the interface, usually by returning an optional, which then propagates to every caller.

Groups and abelian groups​

A group adds inverses, and inverses are what make incremental computation possible: if SjS_j denotes a prefix aggregate, then the aggregate over [i,j)[i, j) is Sj⊖SiS_j \ominus S_i. Without them, range queries need a different data structure entirely, which is developed in the complexity table below.

Commutativity is a separate axis from invertibility and is often confused with it. Permutations under composition form a non-abelian group; the non-negative reals under addition form a commutative monoid that is not a group.

Rings​

A ring has two operations linked by distributivity, and that link is what supports the whole of polynomial and matrix algebra. Two conventions matter in practice. A ring here has a multiplicative identity, since every ring used computationally does. Multiplication is not assumed commutative, since the matrix ring is the standard example and assuming otherwise silently transposes results.

The distinction has computational content. In Z/6Z\mathbb{Z}/6\mathbb{Z} one has 2⋅3=02 \cdot 3 = 0, so a Gaussian elimination step may divide by a quantity that is non-zero and not invertible. Modular algorithms therefore use prime moduli, which makes the ring a field and every non-zero pivot invertible.

Fields​

A field is what almost every numerical algorithm is written against, and the fields that occur are few: Q\mathbb{Q}, R\mathbb{R}, C\mathbb{C}, and Z/pZ\mathbb{Z}/p\mathbb{Z}. The machine approximations to the first three are the subject of Numerical Analysis; the fourth is exact and is the reason modular methods are attractive.

Standard models​

StructureModels used in this reference
Semigroup(R,max⁡)(\mathbb{R}, \max), string concatenation, function composition
Monoid(N,+,0)(\mathbb{N}, +, 0), (R,×,1)(\mathbb{R}, \times, 1), (Fn×n,⋅,I)(F^{n\times n}, \cdot, I), permutations under composition
Group(Z,+,0)(\mathbb{Z}, +, 0), invertible matrices under multiplication, permutations
Abelian group(Z,+)(\mathbb{Z}, +), (Fn,+)(F^n, +), the additive group of any ring
RingZ\mathbb{Z}, Z/mZ\mathbb{Z}/m\mathbb{Z}, F[x]F[x], Fn×nF^{n \times n}
Integral domainZ\mathbb{Z}, F[x]F[x], Z/pZ\mathbb{Z}/p\mathbb{Z}
FieldQ\mathbb{Q}, R\mathbb{R}, C\mathbb{C}, Z/pZ\mathbb{Z}/p\mathbb{Z} for prime pp
Semiring(N,+,×)(\mathbb{N}, +, \times), (R∪{∞},min⁡,+)(\mathbb{R}\cup\{\infty\}, \min, +), Boolean ({0,1},∨,∧)(\{0,1\}, \lor, \land)

The matrix ring Fn×nF^{n\times n} is the standard example of a non-commutative ring, and the polynomial ring F[x]F[x] of a commutative ring that is not a field. Both matter here: the first because matrix algorithms must not assume commutativity, the second because polynomial arithmetic is the subject of Computational Mathematics.

Ordered structures​

Many algorithms need order as well as arithmetic. An ordered field is a field with a total order compatible with the operations:

x≤y⇒x+z≤y+z,0≤x∧0≤y⇒0≤x⋅y.x \le y \Rightarrow x + z \le y + z, \qquad 0 \le x \land 0 \le y \Rightarrow 0 \le x \cdot y.

R\mathbb{R} and Q\mathbb{Q} are ordered fields; C\mathbb{C} is not, and no order on C\mathbb{C} is compatible with its arithmetic. This is why pivoting strategies, which select by magnitude, use ∣z∣|z| rather than zz in the complex case, and why std::totally_ordered appears as a separate requirement from Field in the concepts below.

Semirings​

A semiring is a ring without additive inverses, and dropping inverses is precisely what admits the non-arithmetic examples: the tropical semiring (R∪{∞},min⁡,+)(\mathbb{R}\cup\{\infty\}, \min, +), the Boolean semiring ({0,1},∨,∧)(\{0,1\}, \lor, \land), and the max-plus semiring. Their computational importance is developed in the algebraic path problem below.

Structure-preserving maps​

Homomorphisms​

Homomorphisms are the reason algebraic reasoning transfers between representations. Three of them appear later in this reference.

Reduction modulo mm, Z→Z/mZ\mathbb{Z} \to \mathbb{Z}/m\mathbb{Z}, underlies modular arithmetic and the Chinese remainder technique for controlling coefficient growth.

Evaluation of a polynomial at a point, F[x]→FF[x] \to F, is a ring homomorphism, and it is exactly what Horner's rule computes.

The discrete Fourier transform is an isomorphism from the convolution algebra to the algebra of pointwise products, which is what makes fast convolution possible.

Kernels, ideals and quotients​

For a ring homomorphism ϕ:R→S\phi : R \to S, the kernel ker⁡ϕ={r:ϕ(r)=0}\ker\phi = \{r : \phi(r) = 0\} is an ideal, and R/ker⁡ϕ≅im⁡ϕR/\ker\phi \cong \operatorname{im}\phi. The construction is what makes modular arithmetic a ring rather than a convention: Z/mZ\mathbb{Z}/m\mathbb{Z} is the quotient of Z\mathbb{Z} by the ideal mZm\mathbb{Z}, and its arithmetic is inherited rather than defined afresh.

The computational reading is that a representative-based implementation is correct exactly when the operations are well defined on classes. Reducing after every operation, as the modular type below does, is the mechanical way to guarantee that.

Rounding is not a homomorphism​

Rounding, fl:R→F\mathrm{fl} : \mathbb{R} \to \mathbb{F}, satisfies neither fl(x+y)=fl(x)+fl(y)\mathrm{fl}(x + y) = \mathrm{fl}(x) + \mathrm{fl}(y) nor any weakened version of it that would make the diagram commute. Its failure is the entire subject of Numerical Analysis, and the quantitative substitute for the homomorphism property is the standard model,

fl(x∘y)=(x∘y)(1+δ),∣δ∣≤u.\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta), \qquad |\delta| \le u.

Reading the standard model as a bounded failure of homomorphy is the cleanest way to see why algebraic identities transfer approximately: each application of a law introduces one factor of (1+δ)(1 + \delta), so a rewrite of depth dd costs γd\gamma_d.

Motivation​

Each law licenses a specific algorithmic transformation, and the correspondence is exact enough to be used as a checklist.

What each law licenses​

Associativity: reordering the parenthesisation​

Associativity licenses reassociation, and therefore every divide-and-conquer reduction, every tree-structured parallel sum, and every vectorized accumulation with multiple partial sums. Without it, the only correct evaluation of x1+x2+⋯+xnx_1 + x_2 + \cdots + x_n is the sequential one, and the algorithm is a strictly sequential fold with a dependence chain of length nn.

Commutativity: reordering the operands​

Commutativity licenses reordering, and therefore work-stealing schedulers, unordered reductions over hash containers, and atomic accumulation of contributions in arrival order. An operation that is associative but not commutative, such as matrix multiplication or string concatenation, still admits tree reduction, but the tree must respect the original order of its leaves.

Identity: admitting the empty case​

An identity licenses the empty case. A reduction over an empty range returns the identity, and a parallel decomposition may assign an empty subrange to a worker without special handling. A semigroup without an identity forces every reduction interface to distinguish the empty input, and that distinction propagates through every caller.

Inverses: cancellation and incremental update​

Invertibility licenses cancellation, and therefore sliding-window updates, incremental statistics, and prefix-sum-based range queries: the sum over [i,j)[i, j) is Sj−SiS_j - S_i only because the additive group has inverses. For a monoid without inverses, such as (R,max⁡)(\mathbb{R}, \max), the same query requires a sparse table or a segment tree, at Θ(nlog⁡n)\Theta(n \log n) preprocessing instead of Θ(n)\Theta(n).

Distributivity: factoring and strength reduction​

Distributivity licenses factoring, and therefore Horner's rule, common-subexpression elimination across a product, and Karatsuba-style splitting. The rewrite

a0+a1x+a2x2+⋯+anxn=a0+x(a1+x(a2+⋯ ))a_0 + a_1x + a_2x^2 + \cdots + a_nx^n = a_0 + x(a_1 + x(a_2 + \cdots))

reduces 2n2n multiplications to nn and is valid in any commutative ring.

Idempotence: overlapping decompositions​

An operation is idempotent when a⊕a=aa \oplus a = a, as max⁡\max and min⁡\min are. Idempotence licenses overlapping decompositions: a range query may be answered by combining two overlapping precomputed blocks, which is what makes the sparse table an O(1)O(1)-query structure. For a non-idempotent monoid the blocks must partition the range, and the query costs Θ(log⁡n)\Theta(\log n).

Choosing the weakest sufficient requirements​

The dropping test​

An algorithm is not written for a type. It is written for the weakest structure whose laws its correctness argument uses. The practical test is mechanical: for each law, ask what breaks if it is dropped. If nothing breaks, the law does not belong in the requirements.

Applied to a summation routine: dropping associativity breaks the multiple-accumulator form, dropping the identity breaks the empty case, dropping commutativity breaks nothing, since the implementation combines adjacent partial sums. The requirement is therefore a monoid, not a commutative monoid, and the interface says so.

The cost of over-constraining​

Requiring more than is needed excludes models for no benefit. A reduction constrained to Field rejects (R,max⁡)(\mathbb{R}, \max), string concatenation, matrix products and every semiring application, none of which has a division and all of which reduce correctly. The cost of the mistake is invisible until someone needs one of the excluded models, at which point the interface must change, which is a breaking change for every existing caller.

Where the structures come from in practice​

The structures in this entry are not chosen for elegance; each is the requirement set of a family of computations that recurs in this reference.

StructureWhere it appears
MonoidEvery reduction, scan and parallel aggregation
GroupSliding windows, incremental statistics, prefix-difference queries
RingPolynomial and matrix arithmetic, modular computation
FieldGaussian elimination, least squares, every division-based method
Ordered fieldPivoting, sorting, interval methods, optimisation
SemiringShortest paths, transitive closure, dynamic programming, parsing

Formal formulation​

The refinement lattice​

The structures form a chain of concept refinements ordered by inclusion of their model classes:

Field⊂Ring⊂AbelianGroup⊂Group⊂Monoid⊂Semigroup.\mathcal{F}\text{ield} \subset \mathcal{R}\text{ing} \subset \mathcal{A}\text{belianGroup} \subset \mathcal{G}\text{roup} \subset \mathcal{M}\text{onoid} \subset \mathcal{S}\text{emigroup}.

Every field is a ring; not every ring is a field. Correspondingly, an algorithm requiring only Monoid accepts every model of Field, and an algorithm requiring Field rejects Z\mathbb{Z}.

Separating examples​

The inclusions are strict, and exhibiting a separating example for each is the standard check that a hierarchy is not redundant.

InclusionSeparating example
Semigroup, not monoid(R+,+)(\mathbb{R}^{+}, +) without zero
Monoid, not group(N,+,0)(\mathbb{N}, +, 0)
Group, not abelianPermutations of three elements
Abelian group, not ring(Zn,+)(\mathbb{Z}^n, +) with no product defined
Ring, not integral domainZ/6Z\mathbb{Z}/6\mathbb{Z}
Integral domain, not fieldZ\mathbb{Z}, and F[x]F[x]
Semiring, not ring(N,+,×)(\mathbb{N}, +, \times), and the tropical semiring

A hierarchy level with no separating example is a level that can be removed from the interface without loss, and removing it simplifies every algorithm that mentioned it.

Reduction and the monoid laws​

Generalised associativity​

The forward direction is proved by induction on nn. For n≤2n \le 2 there is one parenthesisation. For larger nn, any parenthesisation has the form L⊕RL \oplus R where LL covers x1,…,xkx_1, \dots, x_k and RR covers xk+1,…,xnx_{k+1}, \dots, x_n; by induction both equal their left-nested forms, and repeated application of the associativity law moves the split point to k=1k = 1. The converse is immediate from the definition of non-associativity.

accumulate against reduce​

The proposition explains a distinction in the standard library that is otherwise arbitrary. std::accumulate fixes the order of evaluation and requires nothing beyond a binary operation; std::reduce is permitted to reassociate and reorder, and its specification therefore requires associativity and commutativity. They are different algorithms on different structures, not two spellings of one, and substituting the second for the first changes the computed value for every floating-point input.

Powering in a monoid​

Linear recurrences as monoid powering​

A linear recurrence uk+1=∑jcjuk−ju_{k+1} = \sum_{j} c_j u_{k-j} of order dd is the action of a fixed d×dd \times d companion matrix, so its nn-th term is obtained by powering that matrix. The Θ(log⁡n)\Theta(\log n) Fibonacci algorithm is the case d=2d = 2 with the matrix whose rows are (1,1)(1,1) and (1,0)(1,0).

The same proposition, applied to a transition matrix over a probability semiring, gives Markov chain evolution; applied to the tropical semiring, it gives shortest paths with a bounded number of edges. One implementation, several unrelated applications, and one correctness argument.

Semirings and the algebraic path problem​

Matrix multiplication over a semiring is defined by the usual formula with (⊕,⊗)(\oplus, \otimes) in place of (+,×)(+, \times):

(A⊗B)ij=⨁k(aik⊗bkj).(A \otimes B)_{ij} = \bigoplus_{k} \left( a_{ik} \otimes b_{kj} \right).

The consequence is practical rather than decorative. A tuned matrix multiplication kernel, developed in Matrix Multiplication for the arithmetic semiring, applies unchanged to the tropical one: the blocking, the packing and the register tiling are all properties of the loop structure and not of the operations. Graph algorithms written this way inherit the memory-hierarchy behaviour of dense linear algebra, which is the entire premise of the GraphBLAS effort.

Where floating point sits​

The quantitative substitute​

The structure (F,⊕,0)(\mathbb{F}, \oplus, 0) is a commutative magma with an identity, and no more. What remains true is a bound: for any parenthesisation PP of x1,…,xnx_1, \dots, x_n,

∣P(x1,…,xn)−∑ixi∣≤γd(P)∑i∣xi∣,\left| P(x_1,\dots,x_n) - \sum_i x_i \right| \le \gamma_{d(P)} \sum_i |x_i|,

where d(P)d(P) is the depth of the parenthesisation tree and γk=ku/(1−ku)\gamma_k = ku/(1-ku). Associativity fails, but it fails by an amount controlled by the shape of the reassociation, which is why tree and multiple-accumulator reductions are not merely acceptable but more accurate than the sequential one.

Which laws survive​

LawHolds for double under ++Holds for unsigned under ++
Closureyes, with ∞\infty and NaNyes
Associativityno, error γd\gamma_dyes
Commutativityyesyes
Identityyes, with signed-zero caveatsyes
Inversesyes, except for ∞\infty and NaNyes
Equality is an equivalenceno, NaN is not reflexiveyes

The table is the reason this reference states obligations rather than assuming them: two types that satisfy the same concept can differ on five of six laws.

C++ representation​

The concept chain​

The chain is stated as a chain, so that each refinement conjoins the concept it refines and subsumption orders the overloads:

#include <concepts>

namespace hpc {

template <class T, class Op>
concept Magma =
std::regular<T> &&
std::invocable<Op, T, T> &&
std::same_as<std::invoke_result_t<Op, T, T>, T>;

// Semantic requirement: op is associative on T.
template <class T, class Op>
concept Semigroup = Magma<T, Op>;

// Semantic requirement: additionally, identity() is a two-sided identity.
template <class T, class Op, class Id>
concept Monoid =
Semigroup<T, Op> &&
std::invocable<Id> &&
std::same_as<std::invoke_result_t<Id>, T>;

template <class T>
concept AdditiveMonoid = requires(T x, T y) {
{ x + y } -> std::same_as<T>;
{ T{} } -> std::same_as<T>;
};

template <class T>
concept AdditiveGroup = AdditiveMonoid<T> && requires(T x, T y) {
{ -x } -> std::same_as<T>;
{ x - y } -> std::same_as<T>;
};

template <class T>
concept Ring = AdditiveGroup<T> && requires(T x, T y) {
{ x * y } -> std::same_as<T>;
{ T{1} } -> std::same_as<T>;
};

template <class T>
concept Field = Ring<T> && requires(T x, T y) {
{ x / y } -> std::same_as<T>;
};

template <class T>
concept OrderedField = Field<T> && std::totally_ordered<T>;

} // namespace hpc

Why Magma and Semigroup are syntactically identical​

The pair is identical syntactically and distinct semantically. The distinction exists to give the associativity requirement a name that algorithms can cite: a routine documented as requiring Semigroup including S1 is a routine that reassociates, and one documented as requiring Magma is one that does not. This is the situation described in Semantic Requirements, and it is why the obligations are numbered rather than described in prose.

The obligations, enumerated​

// Semantic requirements for Monoid<T, Op, Id>:
// S1. Associativity: op(op(a, b), c) == op(a, op(b, c)) for all a, b, c.
// S2. Identity: op(identity(), a) == a == op(a, identity()).
// S3. Purity: op is a function of its arguments alone.
// S4. Equality: == is an equivalence relation preserved by op.
//
// For Ring<T>, additionally:
// S5. Distributivity: a * (b + c) == a * b + a * c.
// S6. Additive inverses: a + (-a) == T{}.
//
// For Semiring<T, Plus, Times>, replace S6 by:
// S6'. Annihilation: zero() * a == a * zero() == zero().

Every algorithm in this reference that reassociates cites S1 explicitly, and floating-point instantiations are documented as violating it, with the bound stated.

Refinement discipline and subsumption​

Each concept above conjoins the one it refines rather than restating its requirements. That is not a stylistic preference: it is what makes the atomic constraints shared, and therefore what allows overload resolution to prefer the more constrained declaration. The mechanism, and the failure mode when requirements are restated instead, are developed in Constraints.

Operations as objects​

Passing the operation and the identity​

Passing the operation and the identity as objects, rather than fixing them to operator+ and T{}, is what lets one reduction serve (R,+,0)(\mathbb{R}, +, 0), (R,max⁡,−∞)(\mathbb{R}, \max, -\infty) and (Fn×n,⋅,I)(F^{n\times n}, \cdot, I):

#include <algorithm>
#include <concepts>
#include <limits>

namespace hpc {

template <class T>
struct plus_monoid {
constexpr T operator()(const T& a, const T& b) const { return a + b; }
static constexpr T identity() { return T{}; }
};

template <std::totally_ordered T>
struct max_monoid {
constexpr T operator()(const T& a, const T& b) const { return std::max(a, b); }
static constexpr T identity() { return std::numeric_limits<T>::lowest(); }
};

} // namespace hpc

max_monoid is a genuinely different model, which is the test that Monoid is a concept rather than a type in disguise.

The identity as a static member​

The identity is a static member rather than a data member because it is a property of the structure, not of an instance. That choice keeps the operation object empty, which matters: an empty class as a template parameter costs nothing, whereas a stateful one is copied into every algorithm frame and prevents the empty base optimisation in composed adaptors.

It also exposes a subtlety in max_monoid. With lowest() as the identity the structure is a monoid on the representable range, but that identity is not the mathematical −∞-\infty unless the type has one, so a reduction over an empty range returns a finite value the caller may not expect.

Semiring operation pairs​

#include <limits>

namespace hpc {

// The tropical semiring: min for addition, plus for multiplication.
// Used by the shortest-path formulation of matrix powering.
struct tropical {
static constexpr double zero() { return std::numeric_limits<double>::infinity(); }
static constexpr double one() { return 0.0; }

static constexpr double add(double a, double b) { return a < b ? a : b; }
static constexpr double multiply(double a, double b) { return a + b; }
};

// The arithmetic semiring, so that one kernel serves both.
struct arithmetic {
static constexpr double zero() { return 0.0; }
static constexpr double one() { return 1.0; }

static constexpr double add(double a, double b) { return a + b; }
static constexpr double multiply(double a, double b) { return a * b; }
};

} // namespace hpc

A matrix multiplication templated on the semiring compiles to two kernels with identical loop structure and different inner operations, and both are blocked and vectorized by the same code.

Testing the concepts​

Two models and a near miss​

A concept is validated by exhibiting its model class: at least two structurally different models, and at least one near miss that is correctly excluded.

#include <complex>
#include <string>

namespace hpc {

static_assert(Field<double>);
static_assert(Field<std::complex<double>>);
static_assert(!Field<int>); // no multiplicative inverses
static_assert(!Ring<std::string>); // concatenation has no additive inverse

static_assert(AdditiveMonoid<std::string>); // but it is a monoid under +

} // namespace hpc

The int case is the one worth stating explicitly: int supports /, and integer division is not the field operation. A Field concept that admitted int would license an algorithm to divide and expect exactness, which is the error the concept exists to prevent.

Implementation​

Reduction with several accumulators​

#include <cstddef>
#include <span>

namespace hpc {

// Requires: Monoid<T, Op> including S1 (associativity) and S2 (identity).
// Four accumulators break the sequential dependence chain, which is
// legitimate only under S1.
template <class T, class Op>
constexpr T reduce(std::span<const T> xs, Op op, T identity) {
T a0 = identity, a1 = identity, a2 = identity, a3 = identity;
std::size_t i = 0;
for (; i + 4 <= xs.size(); i += 4) {
a0 = op(a0, xs[i + 0]);
a1 = op(a1, xs[i + 1]);
a2 = op(a2, xs[i + 2]);
a3 = op(a3, xs[i + 3]);
}
T acc = op(op(a0, a1), op(a2, a3));
for (; i < xs.size(); ++i) {
acc = op(acc, xs[i]);
}
return acc;
}

} // namespace hpc

Applied to double and operator+, this computes a different value from the sequential sum, because floating-point addition is not associative. The difference is bounded: for the sequential sum the standard bound is

∣s^n−∑i=1nxi∣≤(n−1)u∑i=1n∣xi∣+O(u2),\left| \hat{s}_n - \sum_{i=1}^n x_i \right| \le (n-1) u \sum_{i=1}^n |x_i| + O(u^2),

while for the kk-accumulator form the factor (n−1)(n-1) is replaced by approximately n/k+kn/k + k, which is smaller. The unrolled form is both faster and more accurate, and neither fact is visible from the concept.

Exact summation when the bound is not enough​

#include <cmath>
#include <span>

namespace hpc {

// Neumaier compensated summation: accurate to about 2u relative to the sum of
// magnitudes, independently of n, at roughly four times the arithmetic.
constexpr double compensated_sum(std::span<const double> xs) noexcept {
double s = 0.0;
double c = 0.0;
for (const double x : xs) {
const double t = s + x;
c += (std::abs(s) >= std::abs(x)) ? ((s - t) + x) : ((x - t) + s);
s = t;
}
return s + c;
}

} // namespace hpc

The compensation term recovers the rounding error of each addition exactly, which is possible because the error of a correctly rounded sum of two floating-point numbers is itself representable. This is an error-free transformation, developed in Numerical Analysis.

Powering, written once​

#include <cstddef>

namespace hpc {

// Requires: Monoid<T, Op> including S1. Θ(log n) applications of op.
template <class T, class Op>
constexpr T power(T x, std::size_t n, Op op, T identity) {
T result = identity;
while (n > 0) {
if (n % 2 == 1) { result = op(result, x); }
x = op(x, x);
n /= 2;
}
return result;
}

} // namespace hpc

The implementation applies op(result, x) rather than op(x, result): the order is chosen so that the routine is correct for non-commutative monoids, where the two differ. A version that ignores this works for numbers and silently transposes matrix products.

An exact ring: modular arithmetic​

#include <cstdint>

namespace hpc {

// Z/PZ with P a template parameter, so the compiler replaces the division by a
// multiply-shift sequence. An exact model of Ring, and of Field when P is prime.
template <std::uint64_t P>
class modular {
public:
constexpr explicit modular(std::uint64_t v = 0) noexcept : v_{v % P} {}

constexpr modular operator+(modular r) const noexcept {
return modular{v_ + r.v_};
}

constexpr modular operator-(modular r) const noexcept {
return modular{v_ + P - r.v_};
}

constexpr modular operator*(modular r) const noexcept {
return modular{static_cast<std::uint64_t>(
(static_cast<__uint128_t>(v_) * r.v_) % P)};
}

constexpr bool operator==(const modular&) const = default;

constexpr std::uint64_t value() const noexcept { return v_; }

private:
std::uint64_t v_;
};

} // namespace hpc

This type satisfies S1 to S6 exactly, which makes it the right vehicle for testing generic algebraic algorithms. A reduction producing different results for different accumulator counts over modular<P> has a defect, whereas the same observation over double is expected.

A non-commutative model​

#include <array>

namespace hpc {

// 2x2 matrices over a ring: associative, not commutative, with identity.
// Used to test that generic algorithms do not assume commutativity.
template <class T>
struct matrix2 {
std::array<T, 4> a{};

friend constexpr matrix2 operator*(const matrix2& x, const matrix2& y) noexcept {
return matrix2{{
x.a[0] * y.a[0] + x.a[1] * y.a[2],
x.a[0] * y.a[1] + x.a[1] * y.a[3],
x.a[2] * y.a[0] + x.a[3] * y.a[2],
x.a[2] * y.a[1] + x.a[3] * y.a[3],
}};
}

static constexpr matrix2 identity() noexcept {
return matrix2{{T{1}, T{}, T{}, T{1}}};
}

constexpr bool operator==(const matrix2&) const = default;
};

} // namespace hpc

Instantiating power with matrix2<modular<P>> gives an exact, non-commutative test case, and it is the smallest configuration in which an unnoticed commutativity assumption produces a wrong answer rather than a rounding difference.

Property tests for the laws​

#include <cstddef>

namespace hpc::testing {

// Discharges S1 probabilistically over a generator of representative values.
// For double under addition this fails, correctly: the type does not satisfy
// S1, and the algorithms relying on it carry a bound instead.
template <class T, class Op, class Generator>
bool check_associativity(Op op, Generator& gen, std::size_t trials) {
for (std::size_t i = 0; i < trials; ++i) {
const T a = gen(), b = gen(), c = gen();
if (!(op(op(a, b), c) == op(a, op(b, c)))) { return false; }
}
return true;
}

template <class T, class Op, class Generator>
bool check_identity(Op op, T identity, Generator& gen, std::size_t trials) {
for (std::size_t i = 0; i < trials; ++i) {
const T a = gen();
if (!(op(identity, a) == a && op(a, identity) == a)) { return false; }
}
return true;
}

} // namespace hpc::testing

The generator must produce the awkward values deliberately: zero, one, the extremes of the range, and for floating point the infinities, the denormals and NaN. Laws fail at those values and rarely elsewhere, so a naive uniform generator establishes almost nothing.

Complexity​

Complexity; n elements, one operation counted at unit cost
OperationTimeSpaceRemarks
Sequential reductionΘ(n)Θ(1)Latency-bound: n · latency(op)
k-accumulator reductionΘ(n)Θ(k)Throughput-bound for k ≥ latency/throughput
Tree reduction, p processorsΘ(n/p + log p)Θ(p)Requires S1
Compensated summationΘ(n), about 4× flopsΘ(1)Error O(u), independent of n
Exponentiation by squaringΘ(log n) applicationsΘ(1)Requires S1 only
Linear recurrence, order d, term nΘ(d³ log n)Θ(d²)Companion-matrix powering
Prefix sums, sequentialΘ(n)Θ(1)Dependence chain of length n
Prefix sums, work-efficient parallel2n ops, Θ(log n) depthΘ(n)Requires S1
Range query with inversesΘ(1) after Θ(n) prefixΘ(n)Needs a group
Range query, idempotent monoidΘ(1) after Θ(n log n)Θ(n log n)Sparse table; overlapping blocks
Range query, general monoidΘ(log n) after Θ(n)Θ(n)Segment tree; disjoint blocks
All-pairs shortest paths by poweringΘ(n³ log n)Θ(n²)Tropical semiring; Θ(n³) by Floyd–Warshall

The last three rows are the clearest case of laws having complexity consequences: the same query costs Θ(1)\Theta(1), Θ(1)\Theta(1) with more preprocessing, or Θ(log⁡n)\Theta(\log n), according to whether inverses, idempotence, or neither is available.

Performance considerations​

Latency against throughput​

The gap between a latency-bound and a throughput-bound reduction is the ratio of the operation's latency to its reciprocal throughput. On current hardware that is roughly four cycles against one half cycle for floating-point addition, a factor of eight. Recovering it requires associativity, and compilers will not assume associativity for floating-point types unless told to, since it changes results.

This is why -ffast-math produces large speed-ups on reduction-heavy code, and why using it silently changes the error analysis of every reduction in the program. The disciplined alternative is to write the multiple-accumulator form explicitly: the transformation is then local, visible, and analysable, and the compiler needs no permission to change semantics globally.

Interaction with vector width​

Vectorization multiplies the accumulator count by the vector width. With W=8W = 8 lanes and k=4k = 4 accumulators the reduction carries 32 independent partial sums, which saturates the pipelines and changes the error bound from roughly γn\gamma_{n} to roughly γn/32+32\gamma_{n/32 + 32}. Both effects follow from S1, and neither is available without it.

Non-commutative operations​

Matrix multiplication is associative and not commutative, so tree reduction is available and reordering is not. A parallel product of a sequence of matrices must combine adjacent partial products only, which constrains the scheduler: work stealing remains usable, but the combination step must respect the original index order. Algorithms that assume commutativity silently produce transposed or otherwise permuted results, and the defect does not appear on symmetric test data.

Choosing the identity​

For max_monoid the identity is the lowest representable value, so reductions over empty ranges return a finite sentinel. Where that is unacceptable, the interface should return std::optional<T> for possibly-empty inputs and take a non-empty range otherwise, which moves the question from a runtime surprise to a type-level distinction.

Semiring kernels reuse the dense machinery​

Because a semiring matrix product has the loop structure of the arithmetic one, the blocking, packing and register tiling developed in Matrix Multiplication apply unchanged. A tropical-semiring product therefore runs at a comparable fraction of the machine's throughput, and graph algorithms expressed this way inherit the memory-hierarchy behaviour of dense linear algebra rather than the pointer-chasing behaviour of a classical graph traversal.

The qualification is that min⁡\min has no fused multiply-add analogue, so the peak is the issue-limited rate rather than the fused rate, and the achievable fraction is correspondingly lower.

Limitations​

Floating-point types model no algebraic structure exactly. double under addition is a commutative magma with an identity and approximate inverses; it is not a monoid, and every associativity claim about it is an approximation with a stated bound.

Integer types with wrapping arithmetic do form rings, namely Z/2nZ\mathbb{Z}/2^{n}\mathbb{Z} for unsigned types, but signed overflow is undefined, so signed integers model a ring only on a subset of their values, and the subset is not expressible as a concept.

Concepts cannot check laws, so Ring<T> is satisfied by types that violate distributivity. Property-based testing over generated values is the practical substitute and is treated in Correctness.

Equality is presupposed by every law and is itself unreliable: NaN breaks reflexivity, and a type with a lazy or non-normalised representation may have equality that is expensive or inconsistent with the operations.

The hierarchy presented here is the one this reference needs, not the one algebra provides. Modules, algebras over a ring, lattices, categories and universal algebra all have computational uses that are not developed here, and the omission is a choice of scope rather than a claim of irrelevance.

Linear Algebra builds on fields; Concepts supplies the mechanism; Constraints the refinement discipline; Semantic Requirements the treatment of the laws; Numerical Analysis the error bounds quoted above; Computational Mathematics exact arithmetic in these structures.

References​

  • A. Stepanov and D. Rose, From Mathematics to Generic Programming, Addison-Wesley, 2014.
  • A. Stepanov and P. McJones, Elements of Programming, Addison-Wesley, 2009.
  • S. Mac Lane and G. Birkhoff, Algebra, 3rd ed., AMS Chelsea, 1999.
  • N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, chapter 4.
  • A. Neumaier, Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen, ZAMM 54, 1974.
  • J. Kepner and J. Gilbert, Graph Algorithms in the Language of Linear Algebra, SIAM, 2011.
  • M. Gondran and M. Minoux, Graphs, Dioids and Semirings, Springer, 2008.