Library

since C++98

Intermediate

Numeric Algorithms

"C++ numeric algorithms: accumulate, reduce, transform_reduce, prefix scans, inner_product, iota, gcd, lcm, midpoint, and lerp with parallel execution notes."

<numeric> Algorithmssince C++98

The <numeric> header provides fold, prefix-scan, inner-product, and sequence-fill algorithms over ranges; C++17 added parallel-safe reduce, transform_reduce, and scan variants; C++20 added midpoint (in <numeric>) and lerp (in <cmath>).

Overview

<numeric> algorithms fall into four groups:

Folds: accumulate (C++98, ordered), reduce (C++17, reorderable/parallel)
Map-reduce: inner_product (C++98), transform_reduce (C++17)
Prefix scans: partial_sum (C++98), inclusive_scan/exclusive_scan/transform_*_scan (C++17)
Utilities: adjacent_difference (C++98), iota (C++11), gcd/lcm (C++17), midpoint (C++20), lerp in <cmath> (C++20)

Parallel overloads accepting std::execution policies require <execution> and linking against an appropriate parallel backend (e.g., TBB on Linux, built-in on MSVC).

std::accumulate vs std::reduce

accumulate (C++98) is a strict left fold — it evaluates elements left-to-right in order, making it correct for non-associative and non-commutative operations but preventing parallelism.

cpp

#include <numeric>
#include <functional>

std::vector<double> v{1.5, 2.5, 3.0};

// Init type determines accumulator type — use 0.0 for floating-point ranges
double sum  = std::accumulate(v.begin(), v.end(), 0.0);           // 7.0
double prod = std::accumulate(v.begin(), v.end(), 1.0, std::multiplies<>{});

// String concat: order matters — accumulate is the correct tool
std::vector<std::string> words{"foo", "bar", "baz"};
std::string s = std::accumulate(words.begin(), words.end(), std::string{},
    [](std::string acc, const std::string& x) { return std::move(acc) + x; });
// "foobarbaz"

reduce (C++17) permits reordering of operations — correct when the operation is associative and commutative. It exposes execution policy overloads:

cpp

#include <numeric>
#include <execution>

std::vector<double> data(1'000'000, 1.0);

double s1 = std::reduce(data.begin(), data.end(), 0.0);                    // C++17
double s2 = std::reduce(std::execution::par_unseq,                         // C++17
                        data.begin(), data.end(), 0.0);

// WRONG: string concat order is unspecified with reduce
// auto bad = std::reduce(words.begin(), words.end());  // UB / wrong result

std::transform_reduce

One-pass map-reduce — avoids materialising an intermediate transformed range. The two-range overload computes a generalised inner product:

cpp

#include <numeric>

std::vector<double> a{1, 2, 3};
std::vector<double> b{4, 5, 6};

// Dot product: multiply elements pairwise, then sum                        // C++17
double dot = std::transform_reduce(a.begin(), a.end(), b.begin(), 0.0);
// 1*4 + 2*5 + 3*6 = 32.0

// L2 norm squared — unary transform overload
double norm_sq = std::transform_reduce(
    a.begin(), a.end(), 0.0,
    std::plus{},
    [](double x) { return x * x; });
// 1 + 4 + 9 = 14.0

// Parallel weighted sum — safe: + and * are commutative+associative on reals
std::vector<double> weights{0.5, 0.3, 0.2};
double wsum = std::transform_reduce(
    std::execution::par_unseq,
    a.begin(), a.end(), weights.begin(), 0.0);

Prefer transform_reduce over inner_product in C++17+ code — it composes with execution policies.

Prefix Scans

C++17 splits partial_sum (C++98) into explicit inclusive/exclusive variants that accept execution policies and a custom binary operation:

cpp

#include <numeric>

std::vector<int> v{1, 2, 3, 4, 5};
std::vector<int> out(v.size());

// partial_sum (C++98) — inclusive, strictly ordered
std::partial_sum(v.begin(), v.end(), out.begin());
// out: {1, 3, 6, 10, 15}

// inclusive_scan (C++17) — parallel-capable, same result for commutative ops
std::inclusive_scan(std::execution::par, v.begin(), v.end(), out.begin());
// out: {1, 3, 6, 10, 15}

// exclusive_scan (C++17) — identity first, last element excluded
std::exclusive_scan(v.begin(), v.end(), out.begin(), 0);
// out: {0, 1, 3, 6, 10}

// Running product with init value
std::inclusive_scan(v.begin(), v.end(), out.begin(), std::multiplies{}, 1);
// out: {1, 2, 6, 24, 120}

// transform_exclusive_scan (C++17): squares then exclusive prefix sum
std::transform_exclusive_scan(v.begin(), v.end(), out.begin(),
    0, std::plus{}, [](int x) { return x * x; });
// transforms: {1, 4, 9, 16, 25}  →  out: {0, 1, 5, 14, 30}

Exclusive scan is the standard building block for work-offset tables in parallel scatter/gather: each element gets the sum of all preceding elements, useful for packing results into a flat array.

std::inner_product and std::adjacent_difference

Two C++98 algorithms with no C++17 parallel equivalents — use transform_reduce and manual loops instead when parallelism matters.

cpp

#include <numeric>

std::vector<int> a{1, 2, 3};
std::vector<int> b{4, 5, 6};

// inner_product (C++98): generalised dot product
int dot  = std::inner_product(a.begin(), a.end(), b.begin(), 0);
// 1*4 + 2*5 + 3*6 = 32

// Manhattan distance with custom ops
int dist = std::inner_product(a.begin(), a.end(), b.begin(), 0,
    std::plus{},
    [](int x, int y) { return std::abs(x - y); });
// |1-4| + |2-5| + |3-6| = 9

// adjacent_difference (C++98): discrete derivative
// out[0] = in[0]; out[i] = op(in[i], in[i-1])
std::vector<int> prices{100, 103, 97, 110};
std::vector<int> deltas(prices.size());
std::adjacent_difference(prices.begin(), prices.end(), deltas.begin());
// deltas: {100, 3, -6, 13}  (first element copied, rest are differences)

// Custom op: per-mille change
std::adjacent_difference(prices.begin(), prices.end(), deltas.begin(),
    [](int cur, int prev) { return (cur - prev) * 1000 / prev; });
// deltas: {100, 30, -58, 134}

In-place (d_first == first) is valid for adjacent_difference; a shifted-by-one overlap is not.

std::iota, std::gcd, std::lcm

cpp

#include <numeric>

// iota (C++11): fill with consecutive values — classic argsort pattern
std::vector<size_t> idx(data.size());
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(),
    [&](size_t i, size_t j) { return data[i] < data[j]; });
// idx is now a permutation that stably sorts data

// gcd / lcm (C++17): work on integers, take absolute values internally
std::gcd(12, 8)    // 4
std::gcd(-12, 8)   // 4
std::gcd(0, 5)     // 5
std::lcm(4, 6)     // 12

auto reduce_fraction = [](int n, int d) -> std::pair<int, int> {
    int g = std::gcd(n, d);
    return {n / g, d / g};
};
// reduce_fraction(6, 9) → {2, 3}

// GCD of an entire range via accumulate
std::vector<int> vals{12, 18, 30};
int g = std::accumulate(vals.begin(), vals.end(), 0,
    [](int a, int b) { return std::gcd(a, b); });
// g = 6

std::midpoint and std::lerp (C++20)

midpoint lives in <numeric>; lerp lives in <cmath>. Both are C++20.

cpp

#include <numeric>  // midpoint
#include <cmath>    // lerp

// midpoint (C++20): overflow-safe; rounds toward first argument for integers
std::midpoint(0, 7)                    // 3 (toward 0)
std::midpoint(7, 0)                    // 4 (toward 7)
std::midpoint(INT_MAX - 1, INT_MAX)    // INT_MAX - 1  (no overflow)
std::midpoint(0.0, 1.0)               // 0.5

// Works with pointers into the same array
int arr[10];
int* mid = std::midpoint(arr, arr + 10);  // arr + 5

// lerp (C++20): a + t*(b-a), but with exact boundary guarantees
// lerp(a, b, 0) == a exactly; lerp(a, b, 1) == b exactly
std::lerp(0.0, 10.0, 0.25)  // 2.5
std::lerp(0.0, 10.0, 1.0)   // exactly 10.0

// Monotone for t in [0,1] when a <= b — safe for animation blending
struct Vec3 { float x, y, z; };
auto lerp3 = [](Vec3 a, Vec3 b, float t) -> Vec3 {
    return {std::lerp(a.x, b.x, t),
            std::lerp(a.y, b.y, t),
            std::lerp(a.z, b.z, t)};
};
// t outside [0,1] extrapolates — valid but check your domain

Best Practices

Match init type to result type: accumulate(v.begin(), v.end(), 0) yields int even if v is vector<double>. Always pass 0.0, 0LL, or the correct typed value.
Use reduce for arithmetic, accumulate when order matters (string building, running state, non-commutative ops).
Prefer transform_reduce over inner_product for dot products in new C++17 code — it participates in execution policies.
Use exclusive_scan for offset tables — its output[i] holds the sum of all preceding elements, the right shape for parallel compaction indices.
Use iota + sort with a lambda for argsort rather than maintaining parallel index arrays manually.

Common Pitfalls

Integer init silently truncates:

cpp

std::vector<double> v{1.5, 2.5};
auto bad  = std::accumulate(v.begin(), v.end(), 0);    // int: 3
auto good = std::accumulate(v.begin(), v.end(), 0.0);  // double: 4.0

reduce with non-commutative ops gives unspecified order:

cpp

std::vector<std::string> w{"a", "b", "c"};
// UB/unspecified: result could be "bac", "abc", "cab"
auto s = std::reduce(w.begin(), w.end());
// Fix: std::accumulate

Parallel reduce is not bit-identical to sequential accumulate: floating-point addition is not associative, so parallel reductions may produce slightly different rounding. This is expected, correct, and usually acceptable for scientific workloads — but do not use parallel reduce to reproduce a specific numerical sequence.

adjacent_difference first element is always copied: out[0] == in[0] regardless of your custom op. Account for this in downstream consumers.

Quick Reference

Algorithm	Header	Since	Purpose
`accumulate`	`<numeric>`	C++98	Ordered left fold
`inner_product`	`<numeric>`	C++98	Generalised dot product
`partial_sum`	`<numeric>`	C++98	Inclusive prefix sum (ordered)
`adjacent_difference`	`<numeric>`	C++98	Discrete derivative
`iota`	`<numeric>`	C++11	Fill with consecutive values
`reduce`	`<numeric>`	C++17	Reorderable fold, parallel-safe
`transform_reduce`	`<numeric>`	C++17	Map then reduce
`inclusive_scan`	`<numeric>`	C++17	Inclusive scan + parallel
`exclusive_scan`	`<numeric>`	C++17	Exclusive scan + parallel
`transform_inclusive_scan`	`<numeric>`	C++17	Map then inclusive scan
`transform_exclusive_scan`	`<numeric>`	C++17	Map then exclusive scan
`gcd`	`<numeric>`	C++17	Greatest common divisor
`lcm`	`<numeric>`	C++17	Least common multiple
`midpoint`	`<numeric>`	C++20	Overflow-safe midpoint
`lerp`	`<cmath>`	C++20	Linear interpolation

std::gcd and std::lcm (C++17)

std::gcd and std::lcm replace hand-rolled Euclidean loops. Both accept any pair of integers (including signed/unsigned mixes that share a common type), take absolute values internally, and are constexpr.

cpp

#include <numeric>  // C++17

// Basic usage
static_assert(std::gcd(12, 8)   == 4);
static_assert(std::gcd(-12, 8)  == 4);   // sign-agnostic
static_assert(std::gcd(0, 7)    == 7);   // gcd(0,n) == n
static_assert(std::lcm(4, 6)    == 12);
static_assert(std::lcm(0, 6)    == 0);   // lcm(0,n) == 0

// Fraction reduction at compile time
constexpr auto reduce(int n, int d) {
    int g = std::gcd(n, d);
    return std::pair{n / g, d / g};
}
static_assert(reduce(6, 9) == std::pair{2, 3});

// GCD of an entire range
std::vector<int> vals{12, 18, 30, 42};
int g = std::accumulate(vals.begin(), vals.end(), 0,
    [](int a, int b) { return std::gcd(a, b); });   // C++17
// g == 6

// LCM of a range — watch for overflow on large inputs
std::vector<int> small{2, 3, 5};
int l = std::accumulate(small.begin(), small.end(), 1,
    [](int a, int b) { return std::lcm(a, b); });
// l == 30

std::lcm may overflow for large inputs even when the mathematical result fits in the type; prefer std::int64_t or __int128 when inputs are large.

std::midpoint and std::lerp (C++20)

std::midpoint and std::lerp close two long-standing correctness gaps in numeric code. midpoint is overflow-safe for integers and pointers; lerp guarantees exact results at the boundaries t=0 and t=1.

cpp

#include <numeric>  // midpoint  — C++20
#include <cmath>    // lerp      — C++20

// midpoint: rounds toward the first argument for integers
static_assert(std::midpoint(0, 7) == 3);   // toward 0
static_assert(std::midpoint(7, 0) == 4);   // toward 7

// No overflow — contrast with the classic (a+b)/2 bug:
int hi = INT_MAX, lo = INT_MAX - 4;
int naive = (hi + lo) / 2;             // UB: signed overflow
int safe  = std::midpoint(lo, hi);     // INT_MAX - 2, no overflow

// Works with pointers into the same array
int arr[10]{};
int* mid = std::midpoint(arr, arr + 10);   // arr + 5

// Binary search without hand-coded midpoint arithmetic
auto bsearch = [](auto* first, auto* last, auto target) {
    while (first < last) {
        auto* m = std::midpoint(first, last);
        if (*m < target) first = m + 1;
        else if (*m > target) last = m;
        else return m;
    }
    return last;
};

// lerp: linear interpolation with boundary exactness
static_assert(std::lerp(0.0, 10.0, 0.0) == 0.0);   // exact
static_assert(std::lerp(0.0, 10.0, 1.0) == 10.0);  // exact
// std::lerp(0.0, 10.0, 0.25) == 2.5

// Per-component 3D vector lerp for animation blending
struct Vec3 { float x, y, z; };
auto lerp3 = [](Vec3 a, Vec3 b, float t) -> Vec3 {
    return {std::lerp(a.x, b.x, t),
            std::lerp(a.y, b.y, t),
            std::lerp(a.z, b.z, t)};
};

lerp with t outside [0, 1] extrapolates — valid but ensure the domain permits it. Neither function is constrained to [0, 1].

std::iota and range fills

std::iota is the standard way to fill a range with a consecutive arithmetic sequence. It writes value, value+1, value+2, ... using operator++ on the init value, so it works with any incrementable type.

cpp

#include <numeric>   // iota     — C++11
#include <vector>
#include <algorithm>

// Simple integer sequence
std::vector<int> v(6);
std::iota(v.begin(), v.end(), 0);   // {0, 1, 2, 3, 4, 5}
std::iota(v.begin(), v.end(), 1);   // {1, 2, 3, 4, 5, 6}

// Argsort pattern: stable indirect sort via index vector
std::vector<double> data{3.1, 1.4, 1.5, 9.2, 2.6};
std::vector<size_t> idx(data.size());
std::iota(idx.begin(), idx.end(), size_t{0});       // C++11
std::stable_sort(idx.begin(), idx.end(),
    [&](size_t i, size_t j) { return data[i] < data[j]; });
// idx: {1, 2, 4, 0, 3}  — permutation that sorts data ascending

// Works with any type supporting operator++, e.g. enum class
enum class Dir { North, East, South, West };
std::vector<Dir> dirs(4);
std::iota(dirs.begin(), dirs.end(), Dir::North);

// C++23 ranges::iota_view / std::views::iota (lazy, no storage)
// auto seq = std::views::iota(0, 100);  // [0, 100) on demand

For lazy sequences without storage, prefer std::views::iota (C++20) — it composes with range pipelines and avoids the up-front allocation.

Exclusive and inclusive prefix sums

Prefix scans convert a sequence of values into a running-total sequence, and the choice between inclusive and exclusive shifts determines where the identity element lands. Inclusive scan includes element i in output i; exclusive scan excludes it, placing the identity at position 0.

cpp

#include <numeric>    // C++17
#include <execution>

std::vector<int> v{1, 2, 3, 4, 5};
std::vector<int> out(v.size());

// inclusive_scan: out[i] = v[0] + ... + v[i]
std::inclusive_scan(v.begin(), v.end(), out.begin());
// out: {1, 3, 6, 10, 15}

// exclusive_scan: out[i] = v[0] + ... + v[i-1], out[0] = init
std::exclusive_scan(v.begin(), v.end(), out.begin(), 0);
// out: {0, 1, 3, 6, 10}

// Custom binary op — running product
std::inclusive_scan(v.begin(), v.end(), out.begin(),
    std::multiplies<int>{}, 1);   // C++17, init overload
// out: {1, 2, 6, 24, 120}

// Parallel inclusive scan — requires commutative+associative op
std::inclusive_scan(std::execution::par, v.begin(), v.end(), out.begin());

// transform_exclusive_scan: map then exclusive scan in one pass
std::transform_exclusive_scan(v.begin(), v.end(), out.begin(),
    0, std::plus<>{}, [](int x) { return x * x; });   // C++17
// squares: {1,4,9,16,25}  → exclusive sum → out: {0,1,5,14,30}

// Parallel scatter/gather: exclusive scan gives write offsets
// Each thread i writes its chunk starting at out[i]

Exclusive scan is the canonical building block for parallel compaction. Given a boolean mask, exclusive-scan it to get the flat-array write offset for each element that passed the filter.

std::inner_product

std::inner_product is the C++98 generalised dot product — it pairs elements from two ranges with a binary transform, then reduces with another binary op. For new code in C++17+, prefer transform_reduce which supports execution policies.

cpp

#include <numeric>   // C++98

std::vector<int> a{1, 2, 3};
std::vector<int> b{4, 5, 6};

// Default: sum of element-wise products (dot product)
int dot = std::inner_product(a.begin(), a.end(), b.begin(), 0);
// 1*4 + 2*5 + 3*6 = 32

// Custom ops: (reduce_op, transform_op)
// Manhattan distance: sum(|a[i] - b[i]|)
int manhattan = std::inner_product(
    a.begin(), a.end(), b.begin(), 0,
    std::plus<>{},
    [](int x, int y) { return std::abs(x - y); });
// |1-4| + |2-5| + |3-6| = 9

// Count mismatches (Hamming distance on equal-length sequences)
std::string s1 = "AGTCG";
std::string s2 = "AATCC";
int mismatches = std::inner_product(
    s1.begin(), s1.end(), s2.begin(), 0,
    std::plus<>{},
    std::not_equal_to<>{});
// 2 (positions 1 and 4 differ)

// Prefer transform_reduce for parallel-capable equivalent:   // C++17
double dot17 = std::transform_reduce(a.begin(), a.end(), b.begin(), 0.0);

inner_product evaluates strictly left-to-right (like accumulate) — safe for non-commutative transform ops, but it cannot be parallelised.

Adjacent difference

std::adjacent_difference computes the discrete first derivative of a sequence — each output element is the difference between consecutive input elements, making it the inverse of inclusive_scan for the default addition case.

cpp

#include <numeric>   // C++98

std::vector<int> prices{100, 103, 97, 110, 108};
std::vector<int> deltas(prices.size());

// Default op (subtraction): out[0]=in[0], out[i]=in[i]-in[i-1]
std::adjacent_difference(prices.begin(), prices.end(), deltas.begin());
// deltas: {100, 3, -6, 13, -2}

// Custom op: per-mille percentage change
std::adjacent_difference(prices.begin(), prices.end(), deltas.begin(),
    [](int cur, int prev) { return (cur - prev) * 1000 / prev; });
// deltas: {100, 30, -58, 134, -18}

// Fibonacci via self-application (in-place is defined)
std::vector<int> fib(10, 1);
std::adjacent_difference(fib.begin(), fib.end() - 1, fib.begin() + 1,
    std::plus<>{});
// fib: {1, 1, 2, 3, 5, 8, 13, 21, 34, 55}

// Roundtrip: adjacent_difference then inclusive_scan recovers original
std::inclusive_scan(deltas.begin() + 1, deltas.end(),
    deltas.begin() + 1, std::plus<>{}, deltas.front());

out[0] is always a copy of in[0], regardless of the custom binary op. In-place operation (d_first == first) is explicitly allowed; a one-element shift overlap is not.

Floating-point numeric limits

std::numeric_limits<T> from <limits> is the standard way to query type properties at compile time — range, epsilon, infinity, NaN availability, and precision characteristics all live here.

cpp

#include <limits>    // C++98 (constexpr members since C++11)
#include <cmath>

using F = double;

constexpr F eps  = std::numeric_limits<F>::epsilon();     // ~2.22e-16
constexpr F inf  = std::numeric_limits<F>::infinity();    // +∞
constexpr F qnan = std::numeric_limits<F>::quiet_NaN();
constexpr F max  = std::numeric_limits<F>::max();         // ~1.8e+308
constexpr F min  = std::numeric_limits<F>::min();         // smallest positive normalised: ~2.2e-308
constexpr F low  = std::numeric_limits<F>::lowest();      // -max  (C++11)

// Digits of precision
constexpr int dig10 = std::numeric_limits<F>::digits10;   // 15 for double
constexpr int dig2  = std::numeric_limits<F>::digits;     // 53 for double (mantissa bits)

// Approximate equality using epsilon
bool approx_eq(double a, double b, int ulps = 4) {
    return std::abs(a - b)
        <= ulps * std::numeric_limits<double>::epsilon()
               * std::max(std::abs(a), std::abs(b));
}

// Integer limits remain useful alongside <climits>
constexpr int  imax = std::numeric_limits<int>::max();    // 2147483647
constexpr auto umax = std::numeric_limits<unsigned>::max();

// is_iec559: true when the type follows IEEE 754
static_assert(std::numeric_limits<float>::is_iec559);
static_assert(std::numeric_limits<double>::is_iec559);

// has_infinity / has_quiet_NaN guard before using those constants
static_assert(std::numeric_limits<double>::has_infinity);
static_assert(std::numeric_limits<double>::has_quiet_NaN);

Use lowest(), not min(), for the most-negative floating-point value. min() returns the smallest positive normalised number; lowest() (C++11) returns the negative counterpart of max(). For integer types they are identical.

See Also

Execution Policies (<execution>) — seq, par, par_unseq, unseq semantics and thread-safety requirements
Ranges Algorithms — C++20 range-based views composable with numeric algorithms
Fold Expressions — compile-time equivalent for parameter pack reductions
std::valarray — C++98 numeric array with built-in element-wise operators and slice views