Cartesian product for multiple sets at compile time

Viewed 408

I am struggling with an implementation of the Cartesian product for multiple indices with a given range 0,...,n-1.

The basic idea is to have a function:

cartesian_product<std::size_t range, std::size_t sets>()

with an output array that contains tuples that hold the different products

[(0,..,0), (0,...,1), (0,...,n-1),...., (n-1, ..., n-1)]

An simple example would be the following:

auto result = cartesian_product<3, 2>();

with the output type std::array<std::tuple<int, int>, (3^2)>:

[(0,0), (0,1), (0,2), (1,0), (1,1), (1,2), (2,0), (2,1), (2,2)]

My main problem is that my version of the Cartesian product is slow and creates a stack overflow if you choose to have more than 5 sets. I believe that my code has too many recursions and temporary variables.

My implementation (C++17) can be found here: cartesian_product

#include <stdio.h>
#include <iostream>

#include <tuple>


template<typename T, std::size_t ...is>
constexpr auto flatten_tuple_i(T tuple, std::index_sequence<is...>) {

    return std::tuple_cat(std::get<is>(tuple)...);

}

template<typename T>
constexpr auto flatten_tuple(T tuple) {
    return flatten_tuple_i(tuple, std::make_index_sequence<std::tuple_size<T>::value>{});
}

template<std::size_t depth, typename T>
constexpr auto recursive_flatten_tuple(T tuple){

    if constexpr(depth <= 1){
        return tuple;
    }else{
        return recursive_flatten_tuple<depth-1>(flatten_tuple(tuple));
    }
}

template<std::size_t depth, typename T, std::size_t ...is>
constexpr auto wdh(T&& tuple, std::index_sequence<is...>){

    if constexpr (depth == 0) {
        return tuple;
    }else{
        //return (wdh<depth-1>(std::tuple_cat(tuple, std::make_tuple(is)),std::make_index_sequence<sizeof...(is)>{})...);
        return std::make_tuple(wdh<depth-1>(std::tuple_cat(tuple, std::make_tuple(is)), std::make_index_sequence<sizeof...(is)>{})...);
    }
}

template<std::size_t sets, typename T, std::size_t ...is>
constexpr auto to_array(T tuple, std::index_sequence<is...>){

    if constexpr (sets == 0){
        auto t = (std::make_tuple(std::get<is>(tuple)),...);
        std::array<decltype(t), sizeof...(is)> arr = {std::make_tuple(std::get<is>(tuple))...};
        //decltype(arr)::foo = 1;
        return arr;
    }else{
        auto t = ((std::get<is>(tuple)),...);
        std::array<decltype(t), sizeof...(is)> arr = {std::get<is>(tuple)...};
        return arr;
    }
}

template<std::size_t sets, std::size_t ...is>
constexpr auto ct_i(std::index_sequence<is...>){

    if constexpr (sets == 0){

        auto u = std::tuple_cat(wdh<sets>(std::make_tuple(is), std::make_index_sequence<sizeof...(is)>{})...);
        auto arr = to_array<sets>(u, std::make_index_sequence<std::tuple_size<decltype(u)>::value>{});

        return arr;

    }else {

        auto u = std::tuple_cat(wdh<sets>(std::make_tuple(is), std::make_index_sequence<sizeof...(is)>{})...);

        auto r = recursive_flatten_tuple<sets>(u);

        auto d = to_array<sets>(r, std::make_index_sequence<std::tuple_size<decltype(r)>::value>{});


        return d;
    }

}

template<std::size_t range, std::size_t sets>
constexpr auto cartesian_product(){

    static_assert( (range > 0), "lowest input must be cartesian<1,1>" );
    static_assert( (sets > 0), "lowest input must be cartesian<1,1>" );
    return ct_i<sets-1>(std::make_index_sequence<range>{});
}


int main()
{
    constexpr auto crt = cartesian_product<3, 2>();

    for(auto&& ele : crt){

        std::cout << std::get<0>(ele) << " " << std::get<1>(ele) << std::endl;

    }

    return 0;
}
5 Answers

Since I was also working on a solution I thought I post it aswell (although very similar to Artyer's answer). Same premise, we replace the tuple with an array and just iterate over the elements, incrementing them one by one.

Note that the power function is broken, so if you need power values <0 or non-integer types you have to fix it.

template <typename It, typename T>
constexpr void setAll(It begin, It end, T value)
{
    for (; begin != end; ++begin)
        *begin = value;
}

template <typename T, std::size_t I>
constexpr void increment(std::array<T, I>& counter, T max)
{
    for (auto idx = I; idx > 0;)
    {
        --idx;
        if (++counter[idx] >= max)
        {
            setAll(counter.begin() + idx, counter.end(), 0); // because std::fill is not constexpr yet          
        }
        else
        {
            break;
        }
    }
}

// std::pow is not constexpr
constexpr auto power = [](auto base, auto power)
{
    auto result = base;
    while (--power)
        result *= base;
    return result;
};

template<std::size_t range, std::size_t sets>
constexpr auto cartesian_product()
{
    std::array<std::array<int, sets>, power(range, sets)> products{};
    std::array<int, sets> currentSet{};

    for (auto& product : products)
    {
        product = currentSet;
        increment(currentSet, static_cast<int>(range));
    }

    return products;
}

int main()
{
    constexpr auto crt = cartesian_product<5, 3>();

    for (auto&& ele : crt) 
    {
        for (auto val : ele)
            std::cout << val << " ";
        std::cout << "\n";
    }

    return 0;
}

Example

With Boost.Mp11, this is... alright, it's not a one-liner, but it's still not so bad:

template <typename... Lists>
using list_product = mp_product<mp_list, Lists...>;

template <typename... Ts>
constexpr auto list_to_tuple(mp_list<Ts...>) {
    return std::make_tuple(int(Ts::value)...);
}

template <typename... Ls>
constexpr auto list_to_array(mp_list<Ls...>) {
    return std::array{list_to_tuple(Ls{})...};
}

template <size_t R, size_t N>
constexpr auto cartesian_product()
{
    using L = mp_repeat_c<mp_list<mp_iota_c<R>>, N>; 
    return list_to_array(mp_apply<list_product, L>{});
}

With C++20, you can declare the two helper function templates as lambdas inside of cartesian_product, which makes this read nicer (top to bottom instead of bottom to top).

Explanation of what's going on, based on the OP example of cartesian_product<3, 2>:

  • mp_iota_c<R> gives us the list [0, 1, 2] (but as integral constant types)
  • mp_repeat_c<mp_list<mp_iota_c<R>>, N> gives us [[0, 1, 2], [0, 1, 2]]. We just repeat the list, but we want a list of lists (hence the extra mp_list in the middle).
  • mp_apply<list_product, L> does mp_product, which is a cartesian product of all the lists you pass in... sticking the result in an mp_list. This gives you [[0, 0], [0, 1], [0, 2], ..., [2, 2]], but as an mp_list of mp_list of integral constants.
  • At this point the hard part is over, we just have to convert the result back to an array of tuples. list_to_tuple takes an mp_list of integral constants and turns that into a tuple<int...> with the right values. And list_to_array takes an mp_list of mp_lists of integral constants and turns that into an std::array of tuples.

A slightly different approach using just the single helper function:

template <template <typename...> class L,
    typename... Ts, typename F>
constexpr auto unpack(L<Ts...>, F f) {
    return f(Ts{}...);
}

template <size_t R, size_t N>
constexpr auto cartesian_product()
{
    using P = mp_apply_q<
        mp_bind_front<mp_product_q, mp_quote<mp_list>>,
        mp_repeat_c<mp_list<mp_iota_c<R>>, N>>;

    return unpack(P{},
        [](auto... lists){
            return std::array{
                unpack(lists, [](auto... values){
                    return std::make_tuple(int(values)...);
                })...
            };
        });
}

This approach is harder to read though, but it's the same algorithm.

You can do this without recursion easily. Notice that each tuple is the digits of numbers from 0 to range ** sets in base range, so you could increment a counter (Or apply to a std::index_sequence) and calculate each value one after the other.

Here's an implementation (That returns a std::array of std::arrays, which works mostly the same as std::tuples as you can get<N>, tuple_size and tuple_element<N> on a std::array, though if you really wanted you can convert them to std::tuples):

#include <cstddef>
#include <array>

namespace detail {
    constexpr std::size_t ipow(std::size_t base, std::size_t exponent) noexcept {
        std::size_t p = 1;
        while (exponent) {
            if (exponent % 2 != 0) {
                p *= base;
            }
            exponent /= 2;
            base *= base;
        }
        return p;
    }
}

template<std::size_t range, std::size_t sets>
constexpr std::array<std::array<std::size_t, sets>, detail::ipow(range, sets)>
cartesian_product() noexcept {
    constexpr std::size_t size = detail::ipow(range, sets);
    std::array<std::array<std::size_t, sets>, size> result{};
    for (std::size_t i = 0; i < size; ++i) {
        std::size_t place = size;
        for (std::size_t j = 0; j < sets; ++j) {
            place /= range;
            result[i][j] = (i / place) % range;
        }
    }
    return result;
}

Here's a test link: https://www.onlinegdb.com/By_X9wbrI

Note that (empty_set)^0 is defined as a set containing an empty set here, but that can be changed by making ipow(0, 0) == 0 instead of 1

I was trying it out just for fun and I ended with pretty much the same idea as @Timo, just with a different format/style.

#include <iostream>
#include <array>
using namespace std;


template<size_t range, size_t sets>
constexpr auto cartesian_product() {
    // how many elements = range^sets
    constexpr auto size = []() {
        size_t x = range;
        size_t n = sets;
        while(--n != 0) x *= range;
        return x;
    }();

    auto products = array<array<size_t, sets>, size>();
    auto counter = array<size_t, sets>{}; // array of zeroes

    for (auto &product : products) {
        product = counter;

        // counter increment and wrapping/carry over
        counter.back()++;
        for (size_t i = counter.size()-1; i != 0; i--) {
            if (counter[i] == range) {
                counter[i] = 0;
                counter[i-1]++;
            }
            else break;
        }
    }
    return products;
}


int main() {
    auto prods = cartesian_product<3, 6>();
}

I basically have a counter array which I increment manually, like so:

// given cartesian_product<3, 4>
[0, 0, 0, 0]
[0, 0, 0, 1]
[0, 0, 0, 2]
[0, 0, 1, 0] // carry over
...
...
[2, 2, 2, 2]

Pretty much just how you would do it by hand.

Example

If you want it in compile-time, you should only employ compile-time evaluations over compile-time data structures. As @Barry pointed above, using Boost.Mp11 greatly facilitates it. Of course you can do reimplement the relevant fundamental functions in plain C++17 on your own:

#include <iostream>

template<class T> struct Box {
        using type = T;
};


template<class... Types> struct List {};


template<class Car, class Cdr> struct Cons;
template<class Car, class Cdr> using ConsT = typename Cons<Car, Cdr>::type;

template<class Car, class... Cdr> struct Cons<Car, List<Cdr...>>: Box<List<Car, Cdr...>> {};


using Nil = List<>;


template<std::size_t i, class L> struct Nth;
template<std::size_t i, class L> using NthT = typename Nth<i, L>::type;

template<std::size_t i, class... Ts> struct Nth<i, List<Ts...>>: std::tuple_element<i, std::tuple<Ts...>> {};


template<class L> struct Head;
template<class L> using HeadT = typename Head<L>::type;

template<class Car, class... Cdr> struct Head<List<Car, Cdr...>>: Box<Car> {};


template<class L> struct Tail;
template<class L> using TailT = typename Tail<L>::type;

template<class Car, class... Cdr> struct Tail<List<Car, Cdr...>>: Box<List<Cdr...>> {};


template<class... Lists> struct Concat;
template<class... Lists> using ConcatT = typename Concat<Lists...>::type;

template<class T, class... Rest> struct Concat<T, Rest...>: Cons<T, ConcatT<Rest...>> {};

template<class Head, class... Tail, class... Rest> struct Concat<List<Head, Tail...>, Rest...>: Cons<Head, ConcatT<List<Tail...>, Rest...>> {};

template<class... Rest> struct Concat<Nil, Rest...>: Concat<Rest...> {};

template<> struct Concat<>: Box<Nil> {};


template<class T, class Subspace> struct Prepend;
template<class T, class Subspace> using PrependT = typename Prepend<T, Subspace>::type;

template<class T, class... Points> struct Prepend<T, List<Points...>>: Box<List<ConsT<T, Points>...>> {};

template<class T> struct Prepend<T, Nil>: Box<List<List<T>>> {};


template<class Range, class Subspace> struct Product;
template<class Range, class Subspace> using ProductT = typename Product<Range, Subspace>::type;

template<class Range, class Subspace> struct Product: Concat<PrependT<HeadT<Range>, Subspace>, ProductT<TailT<Range>, Subspace>> {};

template<class Subspace> struct Product<Nil, Subspace>: Box<Nil> {};


template<std::size_t i> using IntValue = std::integral_constant<std::size_t, i>;


template<class Seq> struct IntegerSequence;
template<class Seq> using IntegerSequenceT = typename IntegerSequence<Seq>::type;

template<std::size_t... is> struct IntegerSequence<std::index_sequence<is...>>: Box<List<IntValue<is>...>> {};


template<std::size_t n> using Range = IntegerSequenceT<std::make_index_sequence<n>>;


template<std::size_t dimensions, std::size_t range> struct CartesianCube;
template<std::size_t dimensions, std::size_t range> using CartesianCubeT = typename CartesianCube<dimensions, range>::type;

template<std::size_t dimensions, std::size_t range> struct CartesianCube: Product<Range<range>, CartesianCubeT<dimensions - 1, range>> {};

template<std::size_t range> struct CartesianCube<0, range>: Box<Nil> {};


template<std::size_t i> std::ostream &operator<<(std::ostream &s, IntValue<i>) {
        return s << '<' << i << '>';
}

template<class... Ts> std::ostream &operator<<(std::ostream &s, List<Ts...>);

namespace detail_ {

template<class L, std::size_t... is> std::ostream &printList(std::ostream &s, L, std::index_sequence<is...>) {
        return ((s << (is == 0? "" : ", ") << NthT<is, L>{}), ...), s;
}

}

template<class... Ts> std::ostream &operator<<(std::ostream &s, List<Ts...>) {
        return detail_::printList(s << "List{", List<Ts...>{}, std::index_sequence_for<Ts...>{}) << '}';
}

int main() {
        std::cout << CartesianCubeT<2, 3>{} << '\n';
}

Note that CartesianCubeT here is actually a List of Lists of integral_constants. Once you have those, converting them into run-time values is trivial. Note that cartesian_product does not even have to be a function, since the whole data set is evaluated at compile-time it can be a templated value.

Related