How can I create the cartesian product of a vector of vectors?

Question

I've a vector of vectors say vector<vector<int>> of different sizes as follows:

vector<vector<int>> items = {
    { 1, 2, 3 },
    { 4, 5 },
    { 6, 7, 8 }
};

I want to create combinations in terms of Cartesian product of these vectors like

1,4,6
1,4,7
1,4,8
1,5,6
// ...
3,5,7
3,5,8

How can I do that?

Matron

Related question, but requires a fixed number of vectors. – Amund

Answer 1

First, I'll show you a recursive version.

typedef std::vector<int> Vi;
typedef std::vector<Vi> Vvi;

// recursive algorithm to to produce cart. prod.
// At any given moment, "me" points to some Vi in the middle of the
// input data set. 
//   for int i in *me:
//      add i to current result
//      recurse on next "me"
void cart_product(
    Vvi& rvvi,  // final result
    Vi&  rvi,   // current result 
    Vvi::const_iterator me, // current input
    Vvi::const_iterator end) // final input
{
    if(me == end) {
        // terminal condition of the recursion. We no longer have
        // any input vectors to manipulate. Add the current result (rvi)
        // to the total set of results (rvvvi).
        rvvi.push_back(rvi);
        return;
    }

    // need an easy name for my vector-of-ints
    const Vi& mevi = *me;
    for(Vi::const_iterator it = mevi.begin(); it != mevi.end(); it++) {
        // final rvi will look like "a, b, c, ME, d, e, f"
        // At the moment, rvi already has "a, b, c"
        rvi.push_back(*it);  // add ME
        cart_product(rvvi, rvi, me+1, end); add "d, e, f"
        rvi.pop_back(); // clean ME off for next round
    }
}

See full code at Compiler Explorer.

Now, I'll show you the recursive iterative version that I shamelessly stole from @John :

// Seems like you'd want a vector of iterators
// which iterate over your individual vector<int>s.
struct Digits {
    Vi::const_iterator begin;
    Vi::const_iterator end;
    Vi::const_iterator me;
};
typedef std::vector<Digits> Vd;
void cart_product(
    Vvi& out,  // final result
    Vvi& in)  // final result
 
{
    Vd vd;
    
    // Start all of the iterators at the beginning.
    for(Vvi::const_iterator it = in.begin();
        it != in.end();
        ++it) {
        Digits d = {(*it).begin(), (*it).end(), (*it).begin()};
        vd.push_back(d);
    }
    
    while(1) {
        // Construct your first product vector by pulling 
        // out the element of each vector via the iterator.
        Vi result;
        for(Vd::const_iterator it = vd.begin();
            it != vd.end();
            it++) {
            result.push_back(*(it->me));
        }
        out.push_back(result);

        // Increment the rightmost one, and repeat.

        // When you reach the end, reset that one to the beginning and
        // increment the next-to-last one. You can get the "next-to-last"
        // iterator by pulling it out of the neighboring element in your
        // vector of iterators.
        for(Vd::iterator it = vd.begin(); ; ) {
            // okay, I started at the left instead. sue me
            ++(it->me);
            if(it->me == it->end) {
                if(it+1 == vd.end()) {
                    // I'm the last digit, and I'm about to roll
                    return;
                } else {
                    // cascade
                    it->me = it->begin;
                    ++it;
                }
            } else {
                // normal
                break;
            }
        }
    }
}

Plumbaginaceous

Take a close look at main. All of the results are in the output variable. It happens to be a std::vector<std::vector<int> >, but you could easily modify it to be a std::set<std::vector<int> >. You'd need to change rvvi.push_back() to rvvi.insert(). – Boson

++ Thanks for expressing my algorithm in code. And no worries, I won't sue you. ;) – Undry

What you call struct Digits is really just a Digit (singular). I know it seems minor but it could already be hard to follow what is a vector and what isn't (and what is a vector of vectors). – Monopolize

Thanks for these. Has anyone benchmarked them against each other? – Cruise

Looks like the recursive version is a lot faster in my (extremely) limited testing. – Cruise

Answer 2

Here is a solution in C++11.

The indexing of the variable-sized arrays can be done eloquently with modular arithmetic.

The total number of lines in the output is the product of the sizes of the input vectors. That is:

N = v[0].size() * v[1].size() * v[2].size()

Therefore the main loop has n as the iteration variable, from 0 to N-1. In principle, each value of n encodes enough information to extract each of the indices of v for that iteration. This is done in a subloop using repeated modular arithmetic:

#include <cstdlib>
#include <iostream>
#include <numeric>
#include <vector>
using namespace std;

void cartesian( vector<vector<int> >& v ) {
  auto product = []( long long a, vector<int>& b ) { return a*b.size(); };
  const long long N = accumulate( v.begin(), v.end(), 1LL, product );
  vector<int> u(v.size());
  for( long long n=0 ; n<N ; ++n ) {
    lldiv_t q { n, 0 };
    for( long long i=v.size()-1 ; 0<=i ; --i ) {
      q = div( q.quot, v[i].size() );
      u[i] = v[i][q.rem];
    }
    // Do what you want here with u.
    for( int x : u ) cout << x << ' ';
    cout << '\n';
  }
}

int main() {
  vector<vector<int> > v { { 1, 2, 3 },
                           { 4, 5 },
                           { 6, 7, 8 } };
  cartesian(v);
  return 0;
}

Output:

1 4 6 
1 4 7 
1 4 8 
...
3 5 8

Galloway

why #include <cstdlib>? – Fuge

it works now in g++ without explicitly #includeing <cstdlib> header! – Fuge

That's because header files will include other header files. If you try to compile int main() { return div(12,3).rem; } without including any header file, then you'll get a message like error: ‘div’ was not declared in this scope – Galloway

This is amazing! Loved it so much, made a variadic templated version of it: https://mcmap.net/q/282684/-how-can-i-create-the-cartesian-product-of-a-vector-of-vectors – Deliver

Answer 3

Shorter code:

vector<vector<int>> cart_product (const vector<vector<int>>& in) {
    // note: this creates a vector containing one empty vector
    vector<vector<int>> results = {{}};
    for (const auto& new_values : in) {
        // helper vector so that we don't modify results while iterating
        vector<vector<int>> next_results;
        for (const auto& result : results) {
            for (const auto value : new_values) {
                next_results.push_back(result);
                next_results.back().push_back(value);
            }
        }
        results = move(next_results);
    }
    return result;
}

Dip

Answer 4

Seems like you'd want a vector of iterators which iterate over your individual vector<int>s.

Start all of the iterators at the beginning. Construct your first product vector by pulling out the element of each vector via the iterator.

Increment the rightmost one, and repeat.

When you reach the end, reset that one to the beginning and increment the next-to-last one. You can get the "next-to-last" iterator by pulling it out of the neighboring element in your vector of iterators.

Continue cycling through until both the last and next-to-last iterators are at the end. Then, reset them both, increment the third-from-last iterator. In general, this can be cascaded.

It's like an odometer, but with each different digit being in a different base.

Undry

Could you provide an example for the loop ? – Matron

I can explain the principles but it would take me a bit to code it up, as I'm not an STL ninja yet. :) – Undry

Answer 5

Here's my solution. Also iterative, but a little shorter than the above...

void xp(const vector<vector<int>*>& vecs, vector<vector<int>*> *result) {
  vector<vector<int>*>* rslts;
  for (int ii = 0; ii < vecs.size(); ++ii) {
    const vector<int>& vec = *vecs[ii];
    if (ii == 0) {
      // vecs=[[1,2],...] ==> rslts=[[1],[2]]
      rslts = new vector<vector<int>*>;
      for (int jj = 0; jj < vec.size(); ++jj) {
        vector<int>* v = new vector<int>;
        v->push_back(vec[jj]);
        rslts->push_back(v);
      }
    } else {
      // vecs=[[1,2],[3,4],...] ==> rslts=[[1,3],[1,4],[2,3],[2,4]]
      vector<vector<int>*>* tmp = new vector<vector<int>*>;
      for (int jj = 0; jj < vec.size(); ++jj) {  // vec[jj]=3 (first iter jj=0)
        for (vector<vector<int>*>::const_iterator it = rslts->begin();
             it != rslts->end(); ++it) {
          vector<int>* v = new vector<int>(**it);       // v=[1]
          v->push_back(vec[jj]);                        // v=[1,3]
          tmp->push_back(v);                            // tmp=[[1,3]]
        }
      }
      for (int kk = 0; kk < rslts->size(); ++kk) {
        delete (*rslts)[kk];
      }
      delete rslts;
      rslts = tmp;
    }
  }
  result->insert(result->end(), rslts->begin(), rslts->end());
  delete rslts;
}

I derived it with some pain from a haskell version I wrote:

xp :: [[a]] -> [[a]]
xp [] = []
xp [l] = map (:[]) l
xp (h:t) = foldr (\x acc -> foldr (\l acc -> (x:l):acc) acc (xp t)) [] h

Whereon

Thanks for taking the effort. I appreciate you help ! :-) – Matron

In haskell, I would have wrote xp = sequence – Chockablock

Answer 6

Since I needed the same functionality, I implemented an iterator which computes the Cartesian product on the fly, as needed, and iterates over it.

It can be used as follows.

#include <forward_list>
#include <iostream>
#include <vector>
#include "cartesian.hpp"

int main()
{
    // Works with a vector of vectors
    std::vector<std::vector<int>> test{{1,2,3}, {4,5,6}, {8,9}};
    CartesianProduct<decltype(test)> cp(test);
    for(auto const& val: cp) {
        std::cout << val.at(0) << ", " << val.at(1) << ", " << val.at(2) << "\n";
    }

    // Also works with something much less, like a forward_list of forward_lists
    std::forward_list<std::forward_list<std::string>> foo{{"boo", "far", "zab"}, {"zoo", "moo"}, {"yohoo", "bohoo", "whoot", "noo"}};
    CartesianProduct<decltype(foo)> bar(foo);
    for(auto const& val: bar) {
        std::cout << val.at(0) << ", " << val.at(1) << ", " << val.at(2) << "\n";
    }
}

The file cartesian.hpp looks like this.

#include <cassert>

#include <limits>
#include <stdexcept>
#include <vector>

#include <boost/iterator/iterator_facade.hpp>

//! Class iterating over the Cartesian product of a forward iterable container of forward iterable containers
template<typename T>
class CartesianProductIterator : public boost::iterator_facade<CartesianProductIterator<T>, std::vector<typename T::value_type::value_type> const, boost::forward_traversal_tag>
{
    public:
        //! Delete default constructor
        CartesianProductIterator() = delete;

        //! Constructor setting the underlying iterator and position
        /*!
         * \param[in] structure The underlying structure
         * \param[in] pos The position the iterator should be initialized to.  std::numeric_limits<std::size_t>::max()stands for the end, the position after the last element.
         */
        explicit CartesianProductIterator(T const& structure, std::size_t pos);

    private:
        //! Give types more descriptive names
        // \{
        typedef T OuterContainer;
        typedef typename T::value_type Container;
        typedef typename T::value_type::value_type Content;
        // \}

        //! Grant access to boost::iterator_facade
        friend class boost::iterator_core_access;

        //! Increment iterator
        void increment();

        //! Check for equality
        bool equal(CartesianProductIterator<T> const& other) const;

        //! Dereference iterator
        std::vector<Content> const& dereference() const;

        //! The part we are iterating over
        OuterContainer const& structure_;

        //! The position in the Cartesian product
        /*!
         * For each element of structure_, give the position in it.
         * The empty vector represents the end position.
         * Note that this vector has a size equal to structure->size(), or is empty.
         */
        std::vector<typename Container::const_iterator> position_;

        //! The position just indexed by an integer
        std::size_t absolutePosition_ = 0;

        //! The begin iterators, saved for convenience and performance
        std::vector<typename Container::const_iterator> cbegins_;

        //! The end iterators, saved for convenience and performance
        std::vector<typename Container::const_iterator> cends_;

        //! Used for returning references
        /*!
         * We initialize with one empty element, so that we only need to add more elements in increment().
         */
        mutable std::vector<std::vector<Content>> result_{std::vector<Content>()};

        //! The size of the instance of OuterContainer
        std::size_t size_ = 0;
};

template<typename T>
CartesianProductIterator<T>::CartesianProductIterator(OuterContainer const& structure, std::size_t pos) : structure_(structure)
{
    for(auto & entry: structure_) {
        cbegins_.push_back(entry.cbegin());
        cends_.push_back(entry.cend());
        ++size_;
    }

    if(pos == std::numeric_limits<std::size_t>::max() || size_ == 0) {
        absolutePosition_ = std::numeric_limits<std::size_t>::max();
        return;
    }

    // Initialize with all cbegin() position
    position_.reserve(size_);
    for(std::size_t i = 0; i != size_; ++i) {
        position_.push_back(cbegins_[i]);
        if(cbegins_[i] == cends_[i]) {
            // Empty member, so Cartesian product is empty
            absolutePosition_ = std::numeric_limits<std::size_t>::max();
            return;
        }
    }

    // Increment to wanted position
    for(std::size_t i = 0; i < pos; ++i) {
        increment();
    }
}

template<typename T>
void CartesianProductIterator<T>::increment()
{
    if(absolutePosition_ == std::numeric_limits<std::size_t>::max()) {
        return;
    }

    std::size_t pos = size_ - 1;

    // Descend as far as necessary
    while(++(position_[pos]) == cends_[pos] && pos != 0) {
        --pos;
    }
    if(position_[pos] == cends_[pos]) {
        assert(pos == 0);
        absolutePosition_ = std::numeric_limits<std::size_t>::max();
        return;
    }
    // Set all to begin behind pos
    for(++pos; pos != size_; ++pos) {
        position_[pos] = cbegins_[pos];
    }
    ++absolutePosition_;
    result_.emplace_back();
}

template<typename T>
std::vector<typename T::value_type::value_type> const& CartesianProductIterator<T>::dereference() const
{
    if(absolutePosition_ == std::numeric_limits<std::size_t>::max()) {
        throw new std::out_of_range("Out of bound dereference in CartesianProductIterator\n");
    }
    auto & result = result_[absolutePosition_];
    if(result.empty()) {
        result.reserve(size_);
        for(auto & iterator: position_) {
            result.push_back(*iterator);
        }
    }

    return result;
}

template<typename T>
bool CartesianProductIterator<T>::equal(CartesianProductIterator<T> const& other) const
{
    return absolutePosition_ == other.absolutePosition_ && structure_ == other.structure_;
}

//! Class that turns a forward iterable container of forward iterable containers into a forward iterable container which iterates over the Cartesian product of the forward iterable containers
template<typename T>
class CartesianProduct
{
    public:
        //! Constructor from type T
        explicit CartesianProduct(T const& t) : t_(t) {}

        //! Iterator to beginning of Cartesian product
        CartesianProductIterator<T> begin() const { return CartesianProductIterator<T>(t_, 0); }

        //! Iterator behind the last element of the Cartesian product
        CartesianProductIterator<T> end() const { return CartesianProductIterator<T>(t_, std::numeric_limits<std::size_t>::max()); }

    private:
        T const& t_;
};

If someone has comments how to make it faster or better, I'd highly appreciate them.

Subjectify

I don't know why your answer was overlooked, but, at least to my eyes, it looks much more interesting, as it avoids the cost of storing the Cartesian product. Not tried your code yet, but that's what I need. – Kraft

@akim: Unfortunately, it must store it when it's being computed. This is because it needs to return a reference. It wouldn't be hard to change this, but then one could no longer use a standard iterator as far as I see, since they require a reference to be returned. So if you have huge cartesian products, you probably want to go this way and not have nice-to-have things like range based loops. – Subjectify

yes, I agree, some less cute solution is needed. Actually, because I need something with std::tuple of std::vector, I now use something similar to for_imp from this proposal: #13813507, but using C++14-like index_sequences. – Kraft

Answer 7

I was just forced to implement this for a project I was working on and I came up with the code below. It can be stuck in a header and it's use is very simple but it returns all of the combinations you can get from a vector of vectors. The array that it returns only holds integers. This was a conscious decision because I just wanted the indices. In this way, I could index into each of the vector's vector and then perform the calculations I/anyone would need... best to avoid letting CartesianProduct hold "stuff" itself, it is a mathematical concept based around counting not a data structure. I'm fairly new to c++ but this was tested in a decryption algorithm pretty thoroughly. There is some light recursion but overall this is a simple implementation of a simple counting concept.

// Use of the CartesianProduct class is as follows. Give it the number
// of rows and the sizes of each of the rows. It will output all of the 
// permutations of these numbers in their respective rows.
// 1. call cp.permutation() // need to check all 0s.
// 2. while cp.HasNext() // it knows the exit condition form its inputs.
// 3.   cp.Increment() // Make the next permutation
// 4.   cp.permutation() // get the next permutation

class CartesianProduct{
  public:
  CartesianProduct(int num_rows, vector<int> sizes_of_rows){
    permutation_ = new int[num_rows];
    num_rows_ = num_rows;
    ZeroOutPermutation();
    sizes_of_rows_ = sizes_of_rows;
    num_max_permutations_ = 1;
    for (int i = 0; i < num_rows; ++i){
      num_max_permutations_ *= sizes_of_rows_[i]; 
    }
  }

  ~CartesianProduct(){
    delete permutation_;
  }

  bool HasNext(){
    if(num_permutations_processed_ != num_max_permutations_) {
      return true;
    } else {
      return false;
    }
  }

 void Increment(){
    int row_to_increment = 0;
    ++num_permutations_processed_;
    IncrementAndTest(row_to_increment);
  }

  int* permutation(){
    return permutation_;
  }

  int num_permutations_processed(){
    return num_permutations_processed_;
  }
  void PrintPermutation(){
    cout << "( ";
    for (int i = 0; i < num_rows_; ++i){
      cout << permutation_[i] << ", ";
    }
    cout << " )" << endl;
  }

private:
  int num_permutations_processed_;
  int *permutation_;
  int num_rows_;
  int num_max_permutations_;
  vector<int> sizes_of_rows_;

  // Because CartesianProduct is called first initially with it's values
  // of 0 and because those values are valid and important output
  // of the CartesianProduct we increment the number of permutations
  // processed here when  we populate the permutation_ array with 0s.
  void ZeroOutPermutation(){
    for (int i = 0; i < num_rows_; ++i){
      permutation_[i] = 0;
    }

    num_permutations_processed_ = 1;
  }

  void IncrementAndTest(int row_to_increment){
    permutation_[row_to_increment] += 1;
    int max_index_of_row = sizes_of_rows_[row_to_increment] - 1;
    if (permutation_[row_to_increment] > max_index_of_row){
      permutation_[row_to_increment] = 0;
      IncrementAndTest(row_to_increment + 1);
    }
  }
};

Sclerometer

Answer 8

#include <iostream>
#include <vector>

void cartesian (std::vector<std::vector<int>> const& items) {
  auto n = items.size();
  auto next = [&](std::vector<int> & x) {
      for ( int i = 0; i < n; ++ i ) 
        if ( ++x[i] == items[i].size() ) x[i] = 0; 
        else return true;
      return false;
  };
  auto print = [&](std::vector<int> const& x) {
    for ( int i = 0; i < n; ++ i ) 
      std::cout << items[i][x[i]] << ",";
    std::cout << "\b \n";
  };
  std::vector<int> x(n);
  do print(x); while (next(x)); // Shazam!
}

int main () {
  std::vector<std::vector<int>> 
    items { { 1, 2, 3 }, { 4, 5 }, { 6, 7, 8 } };
  cartesian(items);
  return 0;
}

The idea behind this is as follows.

Let n := items.size().
Let m_i := items[i].size(), for all i in {0,1,...,n-1}.
Let M := {0,1,...,m_0-1} x {0,1,...,m_1-1} x ... x {0,1,...,m_{n-1}-1}.

We first solve the simpler problem of iterating through M. This is accomplished by the next lambda. The algorithm is simply the "carrying" routine grade schoolers use to add 1, albeit with a mixed radix number system.

We use this to solve the more general problem by transforming a tuple x in M to one of the desired tuples via the formula items[i][x[i]] for all i in {0,1,...,n-1}. We perform this transformation in the print lambda.

We then perform the iteration with do print(x); while (next(x));.

Now some comments on complexity, under the assumption that m_i > 1 for all i:

• This algorithm requires O(n) space. Note that explicit construction of the Cartesian product takes O(m_0 m_1 m_2 ... m_{n-1}) >= O(2^n) space. So this is exponentially better on space than any algorithm which requires all tuples to be stored simultaneously in memory.
• The next function takes amortized O(1) time (by a geometric series argument).
• The print function takes O(n) time.
• Hence, altogether, the algorithm has time complexity O(n|M|) and space complexity O(n) (not counting the cost of storing items).

An interesting thing to note is that if print is replaced with a function which inspects on average only O(1) coordinates per tuple rather than all of them, then time complexity falls to O(|M|), that is, it becomes linear time with respect to the size of the Cartesian product. In other words, avoiding the copy of the tuple each iterate can be meaningful in some situations.

Penman

Ah, an example why the last note is relevant: suppose we store a multidimensional array as a contiguous array in memory. Now suppose we want to iterate through a "slice" of this array, i.e. restrict each coordinate to subset of the values. Then we want to calculate the address in the original array for each of the entries in the slice in O(1) time rather than O(n) time as we iterate through. – Penman

Answer 9

This version supports no iterators or ranges, but it is a simple direct implementation that uses the multiplication operator to represent the Cartesian product, and a lambda to perform the action.

The interface is designed with the particular functionality I needed. I needed the flexibility to choose vectors over which to apply the Cartesian product in a way that did not obscure the code.

int main()
{
    vector< vector<long> > v{ { 1, 2, 3 }, { 4, 5 }, { 6, 7, 8 } };
    (Cartesian<long>(v[0]) * v[1] * v[2]).ForEach(
        [](long p_Depth, long *p_LongList)
        {
            std::cout << p_LongList[0] << " " << p_LongList[1] << " " << p_LongList[2] << std::endl;
        }
    );
}

The implementation uses recursion up the class structure to implement the embedded for loops over each vector. The algorithm works directly on the input vectors, requiring no large temporary arrays. It is simple to understand and debug.

The use of std::function p_Action instead of void p_Action(long p_Depth, T *p_ParamList) for the lambda parameter would allow me to capture local variables, if I wanted to. In the above call, I don't.

But you knew that, didn't you. "function" is a template class which takes the type parameter of a function and makes it callable.

#include <vector>
#include <iostream>
#include <functional>
#include <string>
using namespace std;

template <class T>
class Cartesian
{
private:
    vector<T> &m_Vector;
    Cartesian<T> *m_Cartesian;
public:
    Cartesian(vector<T> &p_Vector, Cartesian<T> *p_Cartesian=NULL)
        : m_Vector(p_Vector), m_Cartesian(p_Cartesian)
    {};
    virtual ~Cartesian() {};
    Cartesian<T> *Clone()
    {
        return new Cartesian<T>(m_Vector, m_Cartesian ? m_Cartesian->Clone() : NULL);
    };
    Cartesian<T> &operator *=(vector<T> &p_Vector)
    {
        if (m_Cartesian)
            (*m_Cartesian) *= p_Vector;
        else
            m_Cartesian = new Cartesian(p_Vector);
        return *this;
    };
    Cartesian<T> operator *(vector<T> &p_Vector)
    {
        return (*Clone()) *= p_Vector;
    };
    long Depth()
    {
        return m_Cartesian ? 1 + m_Cartesian->Depth() : 1;
    };
    void ForEach(function<void (long p_Depth, T *p_ParamList)> p_Action)
    {
        Loop(0, new T[Depth()], p_Action);
    };
private:
    void Loop(long p_Depth, T *p_ParamList, function<void (long p_Depth, T *p_ParamList)> p_Action)
    {
        for (T &element : m_Vector)
        {
            p_ParamList[p_Depth] = element;
            if (m_Cartesian)
                m_Cartesian->Loop(p_Depth + 1, p_ParamList, p_Action);
            else
                p_Action(Depth(), p_ParamList);
        }
    };
};

Dracaena

Answer 10

The conventional (and fast) solution to this problem is to use a recursive and iterative function:

template <typename Consumer>
void cartesian_n_impl(const std::vector<std::vector<double>> &in,
                      std::vector<double> &stack,
                      std::size_t i,
                      Consumer& consume)
{
    // base case: there is nothing more to loop over, so we do something
    //            with the current stack
    if (i >= in.size()) {
        // std::as_const is necessary because it breaks our algorithm if the user
        // modifies the vector in their consume
        consume(std::as_const(stack));
        return;
    }
    // regular case: we loop over in[i] and recurse, which translates into an
    //               n-dimensional for-loop
    for (double x : in[i]) {
        stack.push_back(x);
        cartesian_n_pick(in, out, stack, i + 1, consume);
        stack.pop_back();
    }
}

// wrapper function with a cleaner interface
template <typename Consumer>
void cartesian_n(const std::vector<std::vector<double>> &in, Consumer consume)
{
    std::vector<double> stack;
    stack.reserve(in.size());
    cartesian_n_impl(in, stack, 0, consume);
}

Of course, you can generalize this and make it more universal with the following changes:

• accept std::span or any range, where appropriate
• don't create a helper std::vector but accept a forward iterator to a range you can use as a stack
• properly constrain these function templates

However, this would be quite verbose and I will leave it as an exercise to the reader. You use this function like:

const std::vector<std::vector<double>> vectors = /* ... */;
std::vecotr<std::vector<double>> items;

cartesian_n(vectors, [](const std::vector<double>& v) {
    items.push_back(v); // or do something else with v, like printing it
});

Ultann

Answer 11

I loved Matt's answer and made it into a template.

The only major difference between this implementation and Matt's is that here, due to how the template pack expansion works, the first argument will change more often than the second argument etc., reshuffling the resulting sequence.

I think it can be fixed but I don't need to for my case :)

template<class F, class... V, size_t... S>
  requires std::invocable<F, V...>
void cartesian_apply(F func, std::array<V, S> const & ... array_args) {
  // based on https://mcmap.net/q/282684/-how-can-i-create-the-cartesian-product-of-a-vector-of-vectors

  // the iteration count is a product of all array sizes
  const long long N = (S * ...);

  for (long long n = 0; n < N; ++n) {
    std::lldiv_t q { n, 0 };

    // we use parameter pack expansion as part of the brace initializer
    // to uniquely pick elements of the operands based on the iteration
    auto apply_tuple = std::tuple<V const &...> { 
      (q = div(q.quot, array_args.size()), array_args.at(q.rem)) 
      ... 
    };

    std::apply(func, apply_tuple);
  }
}

Deliver