Generalized Cartesian Product

Question

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
Thread model: posix
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).

Answer 1

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.