Generalized Cartesian Product

Much to my chagrin, neither STL nor Boost has a cartesian product. Namely, given as arguments one or more iterables, create an iterator producing std::tuples of every combination of elements (with one drawn from each iterator in the order of the arguments).

My goal here was to make cartesian_product to behave exactly like python's itertools.product (with the exception of the optional repeat kwarg). Ideally, it would be a zero cost abstraction (when compared to using nested for loops).

As an example:

#include <vector>

std::vector<int> as = {1, 2};
std::vector<char> bs = {'a', 'b'};
std::vector<float> cs = {1.5, 2.5};

for (auto [a, b, c] : cartesian_product(as, bs, cs)) {
std::cout << "(a = " << a << ", b = " << b << ", c = " << c << ")" << std::endl;
}


Should produce:

(a = 1, b = a, c = 1.5)
(a = 1, b = a, c = 2.5)
(a = 1, b = b, c = 1.5)
(a = 1, b = b, c = 2.5)
(a = 2, b = a, c = 1.5)
(a = 2, b = a, c = 2.5)
(a = 2, b = b, c = 1.5)
(a = 2, b = b, c = 2.5)


To make it general and support an arbitrary number of args, I had to use parameter packs and do some template pattern matching.

Here's what I arrived at:

#pragma once

#include <tuple>

template<typename... Ts>
class product_iterator;

template<typename... Ts>
class product;

template<typename T>
class product<T> {
public:
explicit product(const T &x) : m_x(x) {}

product_iterator<T> begin() const;
product_iterator<T> end() const;

protected:
const T &m_x;
};

template<typename T, typename... Ts>
class product<T, Ts...> {
public:
product(const T &x, const Ts&... xs) : m_x(x), m_xs(product<Ts...>(xs...)) {}

product_iterator<T, Ts...> begin() const;
product_iterator<T, Ts...> end() const;

protected:
const T &m_x;
product<Ts...> m_xs;
};

template<typename T>
class product_iterator<T> {
friend class product<T>;

public:
std::tuple<typename T::value_type> operator*() const {
return std::make_tuple(*m_it);
}

const product_iterator<T> &operator++() {
m_it++;
return *this;
}

bool operator==(const product_iterator &other) const {
return m_it == other.m_it;
}

bool operator!=(const product_iterator &other) const {
return !(*this == other);
}

protected:
typedef typename T::const_iterator t_iterator;

product_iterator(t_iterator it, t_iterator end) : m_it(it), m_end(end) {}

t_iterator m_it;
t_iterator m_end;
};

template<typename T, typename... Ts>
class product_iterator<T, Ts...> {
friend class product<T, Ts...>;

public:
decltype(auto) operator*() const {
return std::tuple_cat(std::make_tuple(*m_x), *m_xs);
}

const product_iterator<T, Ts...> &operator++() {
if (++m_xs == m_xs_end && ++m_x != m_x_end) {
m_xs = m_xs_begin;
}

return *this;
}

bool operator==(const product_iterator &other) const {
return m_x == other.m_x && m_xs == other.m_xs;
}

bool operator!=(const product_iterator &other) const {
return !(*this == other);
}

protected:
typedef typename T::const_iterator t_iterator;
typedef product_iterator<Ts...> ts_iterator;

product_iterator(t_iterator x, t_iterator x_end, ts_iterator xs,
ts_iterator xs_begin, ts_iterator xs_end)
: m_x(x), m_x_end(x_end), m_xs(xs), m_xs_begin(xs_begin), m_xs_end(xs_end) {}

t_iterator m_x;
t_iterator m_x_end;
ts_iterator m_xs;
ts_iterator m_xs_begin;
ts_iterator m_xs_end;
};

template<typename T>
product_iterator<T> product<T>::begin() const {
return product_iterator<T>(m_x.begin(), m_x.end());
}

template<typename T>
product_iterator<T> product<T>::end() const {
return product_iterator<T>(m_x.end(), m_x.end());
}

template<typename T, typename... Ts>
product_iterator<T, Ts...> product<T, Ts...>::begin() const {
return product_iterator<T, Ts...>(m_x.begin(), m_x.end(), m_xs.begin(),
m_xs.begin(), m_xs.end());
}

template<typename T, typename... Ts>
product_iterator<T, Ts...> product<T, Ts...>::end() const {
return product_iterator<T, Ts...>(m_x.end(), m_x.end(), m_xs.end(), m_xs.begin(),
m_xs.end());
}

template<typename... Ts>
product<Ts...> cartesian_product(Ts&... xs) {
return product<Ts...>(xs...);
}


I've tested for correctness and also for speed. Both clang 6 and gcc 8 struggle to optimize this to be equivalent to nested for loops. For empirical results, see this gist with a reproducible benchmark. On my machine, I consistently get around the following:

$g++-8 --version g++-8 (Homebrew GCC 8.1.0) 8.1.0 Copyright (C) 2018 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.$ clang++-6 --version
clang version 6.0.1 (tags/RELEASE_601/final)
Target: x86_64-apple-darwin17.7.0
InstalledDir: /usr/local/opt/llvm/bin

$make clean && CXX=g++-8 make benchmark # ... snip ... time ./run_cartesian 32.53 real 32.42 user 0.03 sys time ./run_loop 32.39 real 32.29 user 0.03 sys$ make clean && CXX=clang++-6 make benchmark
# ... snip ...
time ./run_cartesian
30.31 real        30.24 user         0.02 sys
time ./run_loop
27.30 real        27.24 user         0.02 sys


Not sure what's going on with gcc's timings here (hopefully this benchmark isn't borked!), but there's a noticeable difference with clang. Furthermore, looking at godbolt for both gcc-8 and clang-6 shows both appear to be unable to optimize away some of the abstraction from product_iterator.

Furthermore, if you replace dummy with an accumulator and do an actual dot product like:

uint32_t dot(const std::vector<uint32_t> &as, const std::vector<uint32_t> &bs) {
uint32_t acc = 0;

for (auto [a, b] : cartesian_product(as, bs)) {
acc += a * b;
}

return acc;
}


And:

uint32_t dot(const std::vector<uint32_t> &as, const std::vector<uint32_t> &bs) {
uint32_t acc = 0;

for (auto a : as) {
for (auto b : bs) {
acc += a * b;
}
}

return acc;
}


The different becomes incredibly noticeable. The loop version runtime drops to about 1 second and the cartesian product remains at 30s on my machine. Looking at godbolt for this you can see very clearly that both gcc and clang are able to vectorize the nested loop version but not the cartesian version.

I'm curious about the following things:

• How idomatic is this code?
• How pluggable is this code? (ie. is it compatible with all of the contexts that it should be valid in?)
• Have I implemented everything that I should for this iterator?
• Is there a better way to express this such that gcc and clang can better optimize loops using cartesian_product? Ideally, cartesian_product should be a zero cost abstraction (compared to nested for loops).

Naming

Using m_x for the container (in product) as well as the iterator (in product_iterator) gets confusing. Similar for m_xs.

Performance / Optimization

The compiler relies on the "nestedness" of loops to optimize. The product_iterator removes this, so the compiler is unable to reason about it as perfectly.

For comparison: Clang and GCC are able to generate nearly identical (mostly some different XMM register allocation) assembly for these 2 alternatives (dot product example) as it does for the nested loops one.

Snippet 1

uint32_t dot(const std::vector<uint32_t>& as, const std::vector<uint32_t>& bs) {
uint32_t acc = 0;

acc = std::accumulate(
std::begin(as),
std::end(as),
acc,
[&bs](auto&& acc, auto&& a) {
return std::accumulate(
std::begin(bs),
std::end(bs),
acc,
[&a](auto&& acc, auto&& b) { return acc + a*b; });
});

return acc;
}


Snippet 2

template<typename Callable, typename Cont>
auto cartesian_product(Callable&& op, const Cont& cont) {
for(auto&& x : cont) {
op(x);
}
}

template<typename Callable, typename Cont, typename... Conts>
void cartesian_product(Callable&& op, const Cont& cont, const Conts&... conts) {
for(auto&& x : cont) {
cartesian_product([&x, &op](auto&&... args) { op(x, args...); }, conts...);
}
}

uint32_t dot(const std::vector<uint32_t>& as, const std::vector<uint32_t>& bs) {
uint32_t acc = 0;

cartesian_product([&acc](auto&& a, auto&& b){ acc += a*b; }, as, bs);

return acc;
}


Implementation

• Prefer uniform initialization (using {}) over ().

• product_iterator<T>::m_end isn't needed.

• product_iterator<T, Ts...>::m_x_end isn't needed. (The comparison in operator++ can be removed.)

• Members should be private if possible.

• Is it necessary for product_iterators constructors to be protected? Making them public removes the necessity for those friend class declarations.

• Prefer non-member std::begin/std::end over the corresponding member functions. This allows usage of types where those member functions don't exist (e.g. C style arrays, like int arr[4]).

Questions

• It generally seems OK to me.

• There are some things missing on the iterators for general usage (e.g. with standard library):

• Type aliases iterator_category, value_type, reference, difference_type and pointer (used to deduce iterator traits)

• pointer dereference (operator->())

• post increment (operator++(int))

• Usually a default constructor is provided for iterators (= end iterator).

• Also, you might want to think about keeping a const product<Ts...>* in the general product_iterator instead of m_xs_begin and m_xs_end. After all, some other code (e.g. the one using the iterators) might modify the underlying container(s), which for some containers is well defined behavior.

• See Snippet 2 above.

• Thanks for your answer! I see how the naming could be confusing. I was going for the typical functional list destructuring into (x:xs) and carried it through to indicate how the pieces worked together, but perhaps it was too far removed. I'll take a deeper look into your snippets in a bit, but duly noted. I was hoping that clang might be able to deduce the nested loop given the similar output, but perhaps vectorization occurs earlier. Thanks for the list, I will go about fixing those! Jul 14, 2018 at 7:52
• @BaileyParker: The problem is that the compiler doesn't see the "nestedness" in product_iterator. productor_iterator produces a linear result set, so that's all the compiler can reason about. And since this result set isn't contiguous, vectorization isn't possible. (Maybe this could be improved using coroutines in the future, basically keeping the nested structure but co_yielding individual combinations. Jul 14, 2018 at 7:58
• What do you think would be less confusing for x/xs. I had a hard time coming up with good concise names (current and rest?). To the points you just added: uniform initialization for both product and product_iterator? any particular reason why? You're right about m_end. m_x_end as well, it seems. The constructors could be private as well. I figured they shouldn't be public as only a product should be able to construct a product_iterator. Good tip about std::begin and std::end, I'll use them instead! Jul 14, 2018 at 7:59
• @BaileyParker: My confusion with m_x/m_xs was mostly that I couldn't tell whether it refered to the container or an element inside the container on a glance (in different contexts). Just renaming all instances of m_x tocurrent won't fix this. (Maybe append _iter to your iterator name of choice?) Jul 14, 2018 at 8:03
• Also since std::tuple has no members or member functions (well it has operator= and swap but neither of those make sense, right?), what is the purpose of operator->? Jul 14, 2018 at 8:04