Computational Mathematics

Language revisions

C++20

Class types as non-type template parameters allow exact rational and modular constants to be template arguments.

C++23

A larger constexpr subset permits exact symbolic manipulation during translation.P1938R3

C++26

Reflection permits algebraic expressions to be inspected and rewritten at translation time.P2996

Definition

Computational mathematics is the study of algorithms for mathematical objects, such as polynomials, integers, matrices and algebraic structures, with cost measured in operations of the underlying algebra rather than in machine cycles.

Cost models

The arithmetic model

DefinitionArithmetic complexity

Fix a ring $R$ and count only ring operations, each at unit cost. The arithmetic complexity of a computation is the minimum number of such operations over all algorithms computing it. Bounds in this model are independent of the representation of elements, and translate to a machine cost by multiplying by the cost of an operation in $R$.

The model is useful because it separates two questions that are otherwise entangled: how many operations an algorithm performs, and how expensive each operation is. For fixed-precision arithmetic the second factor is constant and the model predicts running time directly. For multiprecision arithmetic it is not.

The bit model

The complementary model counts bit operations, and it is the one that matters when operands grow.

Operation on $n$-bit integers	Bit complexity
Addition	$\Theta(n)$
Schoolbook multiplication	$\Theta(n^2)$
Karatsuba	$\Theta(n^{1.585})$
Toom-3	$\Theta(n^{1.465})$
Schönhage-Strassen	$O(n \log n \log\log n)$
Harvey and van der Hoeven	$O(n \log n)$
Division, gcd	$O(M(n)\log n)$

An algorithm with a fixed arithmetic complexity can therefore have wildly different bit complexities depending on how its intermediate operands grow, which is the subject of coefficient growth below. The convention $M(n)$ for the cost of multiplication is used because most other operations reduce to it with a logarithmic factor.

Which model applies where

Setting	Predictive model
Fixed-precision numerical work	Arithmetic, with a memory term
Modular arithmetic in machine words	Arithmetic
Rational or multiprecision integer work	Bit
Polynomial work over a small field	Arithmetic
Polynomial work over $\mathbb{Q}$	Bit, with coefficient growth

A bound quoted without its model is not usable, and the mismatch between the two models is the most common source of surprise in exact computation.

Exactness

DefinitionExact computation

A computation is exact when its result is the mathematically correct value, as opposed to an approximation with a bounded error. Exactness requires a representation closed under the operations used: $\mathbb{Z}$ with unbounded integers, $\mathbb{Q}$ with normalised pairs, $\mathbb{Z}/m\mathbb{Z}$ with machine words, or $F[x]$ with coefficient sequences over an exact $F$.

Exactness is not a stronger form of accuracy; it is a different property with a different cost structure. An exact computation has no error bound to state, and its cost grows with the size of the values rather than being fixed by the format.

The three regimes

Regime	Guarantee	Cost
Floating point	Relative error $O(\kappa u)$	Fixed per operation
Interval arithmetic	Rigorous enclosure	2 to 4 times floating point, widening
Exact	Correct value	Grows with operand size

Interval arithmetic is the middle term and is underused: it gives a rigorous bound rather than an estimate, at a constant factor, and its weakness is that naive interval evaluation widens the enclosure at every step, so the bound degrades unless the expression is rearranged to reduce dependency.

Motivation

Two distinct questions are addressed by the same body of theory.

When approximation is not acceptable

Some computations must be exact: determining whether a matrix is singular, deciding the sign of a geometric predicate, verifying an algebraic identity, computing with polynomials whose coefficients must not drift.

The geometric case

A floating-point orientation test on three nearly collinear points can return inconsistent answers for different orderings of the same points. A convex hull algorithm fed inconsistent predicates does not produce a slightly wrong hull: it loops, crashes, or produces a self-intersecting object.

The failure mode is qualitatively different from numerical inaccuracy, because the algorithm's control flow depends on the predicate. This is the general reason that discrete decisions taken from floating-point comparisons are dangerous in a way that floating-point values are not.

Cryptographic and combinatorial computation

Modular exponentiation, polynomial factorisation over finite fields and lattice reduction are exact by nature: an approximate answer is not a worse answer but a meaningless one. Their cost models are the bit model throughout.

When structure beats brute force

Fast algorithms, namely Karatsuba, Strassen and the fast Fourier transform, are algebraic identities that reduce the operation count below the obvious bound. Whether they are worth using is a question about their constants, their stability, and the size at which the asymptotic advantage overcomes them, and the answers differ for each.

Formal formulation

Reductions between problems

The fundamental reductions are between problems, and the reductions are cheap.

Convolution and the Fourier transform

Polynomial multiplication of degree $n$ reduces to the discrete Fourier transform of length $2n$, since the transform diagonalises convolution:

$$\widehat{(f \cdot g)} = \hat{f} \odot \hat{g}, \qquad f \cdot g = \mathrm{DFT}^{-1}\big( \mathrm{DFT}(f) \odot \mathrm{DFT}(g) \big),$$

giving $\Theta(n \log n)$ operations against $\Theta(n^2)$ for the direct method. The statement is an isomorphism of algebras, and it is why the same transform serves signal processing, integer multiplication and polynomial arithmetic.

The number-theoretic transform

Over $\mathbb{Z}/p\mathbb{Z}$ with $p$ chosen so that a primitive $2^k$-th root of unity exists, the same transform is exact: there is no rounding, so the convolution is correct rather than approximately correct. This is the construction used for exact polynomial multiplication and in cryptographic implementations, and it is the reason primes of the form $c\cdot 2^k + 1$ are tabulated in such libraries.

Karatsuba and Toom-Cook

Karatsuba's identity splits a product into three half-sized products rather than four:

$$(a_1 x + a_0)(b_1 x + b_0) = a_1 b_1 x^2 + \big((a_1 + a_0)(b_1 + b_0) - a_1b_1 - a_0b_0\big)x + a_0 b_0,$$

giving $T(n) = 3T(n/2) + \Theta(n)$ and $T(n) = \Theta(n^{\log_2 3})$. Toom-Cook generalises the split to $k$ parts with $2k-1$ products, giving $\Theta(n^{\log_k (2k-1)})$; the exponents decrease towards one as $k$ grows while the constants grow, which is the standard shape of an asymptotic hierarchy with a crossover per level.

Strassen and the exponent of matrix multiplication

Strassen's identity uses seven block products instead of eight:

$$T(n) = 7T(n/2) + \Theta(n^2) \;\Rightarrow\; T(n) = \Theta(n^{\log_2 7}) = \Theta(n^{2.807}).$$

The exponent $\omega = \inf\{c : \text{multiplication in } O(n^c)\}$ satisfies $2 \le \omega < 2.372$, with the upper bounds below Strassen's obtained by methods whose constants make them impractical. Whether $\omega = 2$ is open, and the practical exponent remains $2.807$ or $3$.

Modular methods

Modular arithmetic supplies the standard technique for controlling coefficient growth.

TheoremChinese remainder

For pairwise coprime $m_1, \dots, m_k$ with $M = \prod m_i$, the map

$$\mathbb{Z}/M\mathbb{Z} \to \prod_i \mathbb{Z}/m_i\mathbb{Z}, \qquad x \mapsto (x \bmod m_1, \dots, x \bmod m_k)$$

is a ring isomorphism. A computation whose result is bounded in absolute value by $M/2$ may therefore be performed independently modulo each $m_i$ and reconstructed.

The three practical gains

Each modular computation uses machine-word arithmetic rather than multiprecision. The $k$ computations are independent, hence embarrassingly parallel. And the number of primes needed is determined by an a priori bound on the result, such as Hadamard's bound $|\det A| \le \prod_i \lVert a_i \rVert_2$ on a determinant, so the work is known before the computation begins.

Coefficient growth

Exact arithmetic over $\mathbb{Q}$ suffers explosive intermediate growth. Gaussian elimination over $\mathbb{Q}$ produces entries whose bit lengths grow linearly in the elimination step, so the cost is $\Theta(n^5)$ bit operations rather than the $\Theta(n^3)$ arithmetic operations the model counts.

Two remedies exist: fraction-free elimination, which keeps entries integral and bounded by minors (the Bareiss algorithm), and modular computation followed by reconstruction. Both replace an uncontrolled growth by a controlled one, which is the general shape of every technique in this area.

Exact predicates and filters

For a predicate that is a polynomial in the inputs, such as the sign of a determinant, the standard construction evaluates it in floating point together with a running error bound. If the computed value exceeds the bound in magnitude, its sign is certified and the exact path is not taken; otherwise the exact computation runs.

Static and dynamic filters

A static filter uses a bound computed in advance from the range of the inputs, costs one comparison, and fails whenever the inputs exceed the assumed range. A dynamic filter computes the bound alongside the value, costs a few extra operations, and adapts. Production geometry kernels use a cascade: static filter, then dynamic, then exact.

Adaptive precision

Shewchuk's construction goes further and computes the exact result incrementally, stopping as soon as the sign is determined. The cost is then proportional to the difficulty of the instance rather than to the worst case, which for the near-degenerate inputs that filters fail on is the difference between a usable and an unusable implementation.

C++ representation

Exact rationals

#include <cstdint>
#include <numeric>
#include <stdexcept>

namespace hpc {

// Invariant: denominator > 0 and gcd(|numerator|, denominator) == 1.
// The invariant is established in the constructor and preserved by every
// operation, so equality is representation equality.
class rational {
public:
    constexpr rational(std::int64_t p = 0, std::int64_t q = 1) : num_{p}, den_{q} {
        if (q == 0) { throw std::domain_error{"zero denominator"}; }
        if (q < 0) { num_ = -num_; den_ = -den_; }
        const std::int64_t g = std::gcd(num_ < 0 ? -num_ : num_, den_);
        if (g > 1) { num_ /= g; den_ /= g; }
    }

    constexpr std::int64_t numerator() const noexcept { return num_; }
    constexpr std::int64_t denominator() const noexcept { return den_; }

    constexpr rational operator+(const rational& r) const {
        return rational{num_ * r.den_ + r.num_ * den_, den_ * r.den_};
    }

    constexpr rational operator*(const rational& r) const {
        return rational{num_ * r.num_, den_ * r.den_};
    }

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

private:
    std::int64_t num_;
    std::int64_t den_;
};

}  // namespace hpc

The partiality of the type

The type satisfies Field from Algebra on the subset of values where no intermediate overflows, and that qualification is the entire practical difficulty: exact rational arithmetic over unbounded integers is a field, and over int64_t it is a partial one whose domain of validity depends on the computation rather than on the type.

Eager against lazy normalisation

Normalising in the constructor is a deliberate cost. It makes every operation more expensive by a gcd, which is $O(\log \min(p,q))$, and makes equality a bit comparison. Normalising lazily makes equality expensive and permits unbounded growth in the interim. For a type whose values are compared often, the eager choice is right; for one used in long chains of arithmetic before a single comparison, the lazy choice can be faster.

Modular arithmetic with a compile-time modulus

#include <cstdint>

namespace hpc {

// With P a template parameter, the compiler replaces the division by a
// multiply-shift sequence (Barrett or Montgomery reduction), roughly an order
// of magnitude faster than a hardware divide.
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{static_cast<std::uint64_t>(
            (static_cast<__uint128_t>(v_) * r.v_) % P)};
    }

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

private:
    std::uint64_t v_;
};

}  // namespace hpc

Why the modulus belongs in the type

A runtime modulus forces a hardware division per multiplication, at 20 to 40 cycles. A compile-time modulus lets the compiler emit a multiply and a shift, at 3 to 5 cycles, and the resulting loop vectorises. This is one of the clearest instances of the general principle in Type-Level Programming: moving a value into the type system pays when it unlocks a strength reduction, and here the payoff is nearly an order of magnitude.

Cutoffs as named, measured constants

#include <cstddef>
#include <span>
#include <vector>

namespace hpc {

// Cutoffs determined by measurement on the target, not by the asymptotics.
// The values are recorded with the measurement that produced them.
inline constexpr std::size_t karatsuba_cutoff = 32;
inline constexpr std::size_t toom3_cutoff = 128;
inline constexpr std::size_t fft_cutoff = 1024;

// Dispatches by size across the hierarchy above.
std::vector<double> multiply_polynomials(std::span<const double> a,
                                         std::span<const double> b);

}  // namespace hpc

A literal inside a recursion is a measurement that has lost its provenance. Naming the constants and recording the measurement in Benchmarks keeps the tuning auditable and makes retuning on a new machine a matter of changing three values.

Implementation

Choosing a representation

The representation is chosen from the operations required, in the same way a storage format is chosen from the access pattern in Sparse Computation:

Requirement	Representation
Exact integer arithmetic, moderate size	Fixed-width with overflow checking
Exact integer arithmetic, unbounded	Multiprecision limbs
Exact field arithmetic, bounded result	$\mathbb{Z}/p\mathbb{Z}$ with machine words
Exact rational arithmetic	Normalised pairs, or a single final division
Certified sign only	Floating-point filter with exact fallback
Rigorous enclosure	Interval arithmetic with directed rounding
Polynomial arithmetic, dense	Coefficient array
Polynomial arithmetic, sparse	Term list sorted by exponent

Filters in practice

#include <cmath>

namespace hpc {

// Sign of the 2x2 determinant, with a static error bound. Returns the sign
// when it is certified and calls the exact path otherwise. On generic input
// the exact path is taken rarely, so the amortised cost is that of the
// floating-point evaluation.
inline int orientation(double ax, double ay, double bx, double by,
                       double cx, double cy)
{
    const double detleft  = (ax - cx) * (by - cy);
    const double detright = (ay - cy) * (bx - cx);
    const double det = detleft - detright;

    const double sum = std::abs(detleft) + std::abs(detright);
    constexpr double relative_bound = 3.0 * 1.1102230246251565e-16;  // 3u
    if (std::abs(det) > relative_bound * sum) {
        return det > 0.0 ? 1 : -1;
    }
    return exact_orientation(ax, ay, bx, by, cx, cy);
}

}  // namespace hpc

The bound is derived from the standard model of Numerical Analysis, and deriving it correctly is the whole content of the technique: a bound that is too small returns wrong signs, and one that is too large makes the exact path common and the filter pointless.

Modular determinant

#include <cstdint>
#include <span>
#include <vector>

namespace hpc {

// Exact determinant of an integer matrix by modular computation:
//   1. Bound |det| by Hadamard's inequality.
//   2. Choose primes p_1 ... p_k with prod p_i > 2 |det|_max.
//   3. Compute det mod p_i by Gaussian elimination over Z/p_iZ (independent).
//   4. Reconstruct by the Chinese remainder theorem and symmetric lift.
//
// Cost: Theta(k n^3) word operations, fully parallel over the primes, against
// Theta(n^5) bit operations for elimination over the rationals.
std::int64_t exact_determinant(std::span<const std::int64_t> A, std::size_t n);

}  // namespace hpc

The structure is characteristic of the whole area: an a priori bound on the answer determines how much work is needed, the work decomposes into independent word-sized pieces, and reconstruction is cheap.

Interval arithmetic

#include <cfenv>

namespace hpc {

// A rigorous enclosure. Correctness requires directed rounding, so the
// rounding mode is set for the duration and restored, and the compiler must be
// prevented from reordering arithmetic across the mode changes.
struct interval {
    double lo;
    double hi;

    friend interval operator+(interval a, interval b) noexcept {
        std::fesetround(FE_DOWNWARD);
        const double lo = a.lo + b.lo;
        std::fesetround(FE_UPWARD);
        const double hi = a.hi + b.hi;
        return {lo, hi};
    }
};

}  // namespace hpc

The dependency problem limits the technique: evaluating $x - x$ over an interval gives an enclosure of width twice that of $x$ rather than zero, because the two occurrences are treated as independent. Rearranging expressions to reduce repeated occurrences is the main skill in using interval methods, and it is why naive interval evaluation of a long computation produces useless bounds.

Complexity

Complexity

Operation	Time	Space	Remarks
Schoolbook polynomial multiplication	Θ(n²)	Θ(n)	Best below the cutoff
Karatsuba	Θ(n^{1.585})	Θ(n log n)	Crossover typically n = 32 to 64
Toom-3	Θ(n^{1.465})	Θ(n)	Crossover above Karatsuba
FFT-based convolution	Θ(n log n)	Θ(n)	Crossover about n = 10³; error O(u log n)
Number-theoretic transform	Θ(n log n)	Θ(n)	Exact; requires a suitable prime
Schoolbook matrix multiplication	2n³ flops	Θ(1)	Near-peak throughput when blocked
Strassen	Θ(n^{2.807})	Θ(n²) extra	Crossover n = 500 to 2000; weaker stability
Modular determinant, k primes	Θ(k n³) word ops	Θ(n²)	Independent and parallel; exact
Rational Gaussian elimination	Θ(n⁵) bit ops	Growing	Coefficient explosion
Fraction-free (Bareiss) elimination	Θ(n³) ops on bounded entries	Θ(n²)	Entries are minors
Extended gcd of n-bit integers	O(M(n) log n)	Θ(n)	Half-gcd; Θ(n²) by Euclid
Modular exponentiation, n-bit	Θ(n) multiplications	Θ(n)	Square and multiply
Filtered predicate, generic input	Θ(1) floating point	Θ(1)	Exact path taken rarely
Filtered predicate, degenerate input	Exact path cost	Θ(1)	Adversarial input defeats the filter
Interval arithmetic	2 to 4× floating point	2× storage	Rigorous; widens with dependency

Performance considerations

Asymptotically fast is not fast

Asymptotically fast algorithms lose to direct ones at small sizes, and the crossover is a property of the machine. Strassen's advantage at $n = 4096$ is real, roughly 20 percent against a tuned matrix multiplication, and it requires $\Theta(n^2)$ temporary storage and gives up the elementwise backward stability of the direct method. Whether the trade is acceptable depends on the application; in this reference it is not taken without stating the change.

Operand growth, not operation count, dominates exact arithmetic

This is why the arithmetic model mispredicts so badly here: an algorithm performing $\Theta(n^3)$ ring operations on operands whose size grows linearly performs $\Theta(n^5)$ bit operations. Modular methods and fraction-free elimination both address the growth, the former by keeping every intermediate in a machine word.

Compile-time moduli

Templating on the modulus converts a division per operation into a multiply and a shift, and the resulting kernels vectorise. The payoff is a factor of about ten on the inner loop, which is larger than most memory-hierarchy optimisations and requires only that the modulus be known at translation time.

Allocation, not bandwidth

Memory behaviour differs from the numerical case. Exact algorithms with growing operands are allocation-bound rather than bandwidth-bound, and the effective optimisation is arena allocation with reuse rather than layout or blocking. Profiling such code by cache miss rate is looking at the wrong counter; the right one is allocations per second.

Parallelism is free in the modular setting

The $k$ modular computations are entirely independent, with no communication until reconstruction. This is the rare case where a parallel implementation is trivial and scales to as many workers as there are primes, and it is a strong argument for the modular formulation independently of its sequential cost.

Limitations

Exact arithmetic is unavailable for transcendental quantities: $\pi$, $e^x$ and the values of special functions have no exact finite representation, and computations involving them are approximations regardless of how the surrounding algebra is carried out.

Fast algorithms are numerically weaker. FFT-based convolution has error $O(u \log n)$ against $O(un)$ for the direct method in the worst case but loses the componentwise bound. Strassen is norm-wise stable but not componentwise stable, so small entries of a product may be computed with no correct digits.

Arithmetic complexity ignores the cost of an operation, which for multiprecision operands grows with operand size. Bounds in the model are correct and must be converted before they predict running time.

Filtered predicates are defeated by adversarial or highly degenerate input, where the exact path is taken every time and the amortised argument fails. Applications that must bound worst-case latency cannot rely on the filter.

Interval arithmetic produces rigorous but often useless bounds, because the dependency problem widens the enclosure at every step. Making it useful requires rewriting expressions, which is manual work that does not transfer between problems.

Related entries

Algebra for the structures; Numerical Analysis for the stability of the fast methods and for the filter bounds; Computational Complexity for the machine-level cost models; Type-Level Programming for compile-time parameters such as the modulus; Matrix Multiplication for the kernel Strassen competes against.

References

• J. von zur Gathen and J. Gerhard, Modern Computer Algebra, 3rd ed., Cambridge University Press, 2013.
• D. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., Addison-Wesley, 1997.
• V. Strassen, Gaussian Elimination is not Optimal, Numerische Mathematik 13, 1969.
• J. Alman and V. Vassilevska Williams, A Refined Laser Method and Faster Matrix Multiplication, SODA 2021.
• D. Harvey and J. van der Hoeven, Integer multiplication in time O(n log n), Annals of Mathematics 193(2), 2021.
• N. J. Higham, Exploiting Fast Matrix Multiplication Within the Level 3 BLAS, ACM TOMS 16(4), 1990.
• J. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete and Computational Geometry 18, 1997.
• E. Bareiss, Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination, Mathematics of Computation 22, 1968.