Cartesian for-each at runtime

Useful when the number of arguments to a carteisan product is determined at runtime.

Cartesian product

In mathematics, a cartesian product, also called the product set, set direct product, or cross product is defined for two sets, \(A\) and \(B\), as a set: \[A\times B := \{(x,y): \forall x \in A\ \forall y \in B\}\] Let’s say we have two vectors in C++, \(\{0,1,2\}\) and \(\{3,4\}\), then their cartesian product is a vector of two-tuples : \(\{(0,3),(0,4),(1,3),(1,4),(2,3),(2,4)\}\). The cross product of \(n\) vectors is relatively easy to generate either by using nested for loops or using template meta-programming as long as the value of \(n\) is known at compile time.

Cartesian product at runtime

Sometimes one might need to generate the cartesian product of \(n\) vectors where \(n\) is determined at runtime. If we know the max value \(n\) is ever going to take at runtime, we could still resort to template meta-programming to generate all possible function definitions for the value of \(n\) from \(2\) through the max vlaue—of course this technique has its pros and cons.

Here we focus on how to dynamically generate a cartesian product of \(n\) vectors. To keep the post simple we speak in terms of vectors, however, the scheme should apply to other ‘vector-like’ data types.

Breadth-first algorithm (pseudo-code)

Breadth-first algorithm is straight-forward to think of.

// procedure: extend_cp
// input: cps, es
// output: cps extended

result <- {}
for p in cps
  for e in es
    p_ <- p
    p_.append(e)
    result.append(p_)
return result

// procedure: cartProdBf
// input: vov vector<vector<T>> vector of vectors (for eg.)
// output: vov_ vector<vector<T>> cart prod result
if vov.empty
  return {}

//
// <| front |>
//  |   .   | . | . | . . . | . | (vov)
//          |<      tail       >|
//
init <- {{e} for e in vov.front} // cart prod seed

// see https://en.cppreference.com/w/cpp/algorithm/accumulate
return accumulate(vov.tail, init, extend_cp) 

Better approach: depth first algorithm

The problem with the breadth-first algorithm is that all the members of the cartesian product are generated explosively at once. It would be better if we could generate one member of the product at a time and do something with it (for example call a function on it). I have discovered an algorithm that involves neat tricks using while loops and iterator arithmetics. I did not come up with this myself—borrowed from my colleague Andrey Asadchev’s work in TiledArray.

//
// filename: cartesian_foreach.hpp
//
// required headers:
//   - vector
//   - utility
//   - functional

template <typename R, typename F>
void cartesian_foreach(std::vector<R> const& rs, F&& call_back) {
  using It = decltype(std::begin(rs[0]));
  using T = typename R::value_type;

  if (rs.empty()) return;

  std::vector<It> its, ends;
  its.reserve(rs.size());
  ends.reserve(rs.size());

  for (auto&& r : rs) {
    its.push_back(std::begin(r));
    ends.push_back(std::end(r));
  }

  while (its.front() != ends.front()) {
    std::vector<T> p;
    p.reserve(its.size());

    for (auto&& it: its)
      p.push_back(*it);

    // do something with the cartesian product
    std::invoke(std::forward<F>(call_back), p);

    size_t i = its.size();
    while (i > 0) {
      --i;
      ++its[i];
      if (i == 0) break;
      if (its[i] != ends[i]) break;
      its[i] = std::begin(rs[i]);
    }
  }
}

Example use

//
// file name: main.cpp
//
// required headers:
//   - cartesian_foreach.hpp
//   - vector
//   - iterator
//   - iostream

// prints a vector of printables
template <typename T>
std::ostream& operator<<(std::ostream& os, std::vector<T> const& vec) {
  os << "{";
  if (!vec.empty()) {
    os << vec.front();
    for (auto it = std::next(vec.begin(),1); it != vec.end(); ++it)
      os << ", " << *it;
  }
  os << "}";
  return os;
}

//
// main function
//
int main(int argc, char* argv[]) {
  using std::cout;
  using std::endl;

  auto vov = std::vector<std::vector<int>>{{1, 2}, {3, 4}, {5, 6}};
  
  auto print_vec = [](auto const& v) {
    cout << v << endl;
  };

  cartesian_foreach(vov, print_vec);

  cartesian_foreach(decltype(vov){}, print_vec); // empty vov, no output

  return 0;
}

Output:

{1, 3, 5}
{1, 3, 6}
{1, 4, 5}
{1, 4, 6}
{2, 3, 5}
{2, 3, 6}
{2, 4, 5}
{2, 4, 6}