Nested parameter pack expansion

Question

Please see code snippets (implementation of matrix multiplication) below.
Is it possible to simplify them using nested pack expansion to have something like
{{((a[r][k] * b[k][c]) + ...)...}...}?

#include <array>
#include <utility>

template<typename T, size_t R, size_t C>
using Matrix = std::array<std::array<T, C>, R>;

template<typename A, typename B>
using mul_el_t = decltype(std::declval<A>()[0][0] * std::declval<B>()[0][0]);

Helper to compute single element.

template<size_t R1, size_t C2, size_t... C1_R2, typename A, typename B>
auto _mat_mul_element(const A &a, const B &b, std::index_sequence<C1_R2...>)
{
    return ((a[R1][C1_R2] * b[C1_R2][C2]) + ...);
}

Helper to compute particular row.

template<size_t R1, size_t... C2, typename C1_R2, typename A, typename B>
auto _mat_mul_row(const A &a, const B &b, std::index_sequence<C2...>, C1_R2 c1_r2)
-> std::array<mul_el_t<A, B>, sizeof...(C2)>
{
    return {_mat_mul_element<R1, C2>(a, b, c1_r2)...};
}

This computes whole matrix using parameters packs.

template<size_t... R1, typename C2, typename C1_R2, typename A, typename B>
auto _mat_mul(const A &a, const B &b, std::index_sequence<R1...>, C2 c2, C1_R2 c1_r2)
-> Matrix<mul_el_t<A, B>, sizeof...(R1), C2::size()>
{
    return {_mat_mul_row<R1>(a, b, c2, c1_r2)...};
}

And actual interface.

template<typename T, size_t R1, size_t C1_R2, size_t C2>
Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1_R2> &a, const Matrix<T, C1_R2, C2> &b)
{
    return _mat_mul(
            a, b,
            std::make_index_sequence<R1>{},
            std::make_index_sequence<C2>{},
            std::make_index_sequence<C1_R2>{}
    );
};

UPDATE (looks like I was not clear about the actual problem I have)

When I am trying to replace _mat_mul with:

template<size_t... R1, size_t...  C2, size_t...  C1_R2, typename A, typename B>
auto _mat_mul(const A &a, const B &b,
                std::index_sequence<R1...>,
                std::index_sequence<C2...>,
                std::index_sequence<C1_R2...>)
-> Matrix<mul_el_t<A, B>, sizeof...(R1), sizeof...(C2)>
{
    return {{((a[R1][C1_R2] * b[C1_R2][C2]) + ...)...}...};
}

using Apple LLVM version 9.1.0 (clang-902.0.39.1) compilation fails with:

[ 50%] Building CXX object CMakeFiles/main.cpp.o
main.cpp:38:51: error: pack expansion does not contain any unexpanded parameter packs
    return {{((a[R1][C1_R2] * b[C1_R2][C2]) + ...)...}...};
             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^

I think the failure is expected since compiler doesn't know which pack to expand (R1, C2 or C1_R2) in each expansion block. How can I hint the compiler in this situation (note, I can use any compiler)?

Answer 1

According to documentation nested pack expansion can be seen as an iterative process which starts with innermost pack expansion [3 dots]. Each pack expansion expands all parameter packs within subexpression which is contained by that pack expansion.

Thus {{((a[R1][C1_R2] * b[C1_R2][C2]) + ...)...}...} after first step becomes
{{(a[0][0] * b[0][0] + a[1][1] * b[1][1])...}...} (R1/C2/C1_R2 are index_sequence<2>). So next two pack expansions just have nothing to expand.

Naive solution

It's possible to move each parameter pack to required subexpression at the same time leaving actual value it carries in the desired place. One can use analogue of FP's let ... in expression:

auto let = [](auto a, auto f) { return f(a); };

Thus original expression becomes:

{let(R1, [&](auto r1) {
    return std::array<T, sizeof...(C2)>{let(C2, [&](auto c2) {
        return ((a[r1][C1_R2] * b[C1_R2][c2]) + ...);
    })...};
})...};

That is probably good enough, but may take some time to decipher what's going on there. Also one would still need to introduce those parameter packs in the scope.

Better solution

One could try to improve readability by abstracting over parameter pack introduction/expansion. Using following utility function.

template<typename H, typename F, typename T, T... I>
decltype(auto) repack(std::integer_sequence<T, I...>, H h, F f)
{
    return h(f(std::integral_constant<T, I>{})...);
}

This function takes value which carries some pack (one could make an overload for something other than std::integer_sequence), function f which is applied to each element of the pack, and function h which is used to convert final pack to some value.

Thus full multiply routine becomes

template<typename T, size_t R1, size_t C1_R2, size_t C2>
Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1_R2> &a, const Matrix<T, C1_R2, C2> &b)
{
    std::make_index_sequence<R1> r1{};
    std::make_index_sequence<C2> c2{};
    std::make_index_sequence<C1_R2> c1_r2{};

    return repack(r1, ctor<Matrix<T, R1, C2>>(), [&](auto r1) {
        return repack(c2, ctor<std::array<T, C2>>(), [&](auto c2) {
            return repack(c1_r2, sum, [&](auto c1_r2) {
                return a[r1][c1_r2] * b[c1_r2][c2];
            });
        });
    });
}

where ctor is

template<typename H>
auto ctor()
{
    return [](auto... xs) { return H{xs...}; };
}

and sum = [](auto... xs) { return (xs +...); };.

Crazy solution

One might have spot pattern in expression with 3 nested repack's, thus multiplication routine may become:

template<typename T, size_t R1, size_t C1_R2, size_t C2>
Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1_R2> &a, const Matrix<T, C1_R2, C2> &b)
{
    auto item = [&](auto r1, auto c2, auto c1_r2) { return a[r1][c1_r2] * b[c1_r2][c2]; };
    auto curried_repack = curry(POLY(repack));
    auto m = curried_repack(std::make_index_sequence<R1>{}, ctor<Matrix<T, R1, C2>>());
    auto r = curried_repack(std::make_index_sequence<C2>{}, ctor<std::array<T, C2>>());
    auto e = curried_repack(std::make_index_sequence<C1_R2>{}, sum);

    auto op = [](auto w, auto f) {
        return compose(w, curry(f));
    };
    return foldr(op, m, r, e, item)();
}

With utilities:

template<typename F>
auto curry(F f)
{
    return [=](auto... a) {
        return [=](auto... b) { return f(a..., b...); };
    };
};

template<typename F, typename G>
auto compose(F f, G g)
{
    return [=](auto... xs) {
        return f(g(xs...));
    };
};

And macro to convert template function into value

#define POLY(f) ([](auto... a){ return f(a...); })

And foldr which is left as a homework.

Compilation note

All solutions are equivalent in sense they produce the same binary.