Simplify a longer expanded loop expression by metaprogramming constexpr or template

Question

It looks like I have a longer expression (expanded loop), such as the code below, several times in software that inflates it with a few thousand lines.

Since it polytakes template arguments for performance (the second argument corresponds to the values ​​of the z loops ), I wonder if it is possible to simplify the code below by metaphorically programming the template and build something like a loop through recursion. The syntax of the expression is apparently for each x=bx (a + b + c * by * bz) + ..

I think it would be simpler if it were polyn’t for the function of the template, but would fulfill the functional parameters.

void calc(float mat[3][3][3], float fS, float fT, float fU)
{
    const float bs20_u = poly<2, 0>(fU);
    const float bs21_u = poly<2, 1>(fU);
    const float bs22_u = poly<2, 2>(fU);

    const float bs20_s = poly<2, 0>(fS);
    const float bs21_s = poly<2, 1>(fS);
    const float bs22_s = poly<2, 2>(fS);

    const float bs20_t = poly<2, 0>(fT);
    const float bs21_t = poly<2, 1>(fT);
    const float bs22_t = poly<2, 2>(fT);

    float result = 
       ((mat[0][0][0] * bs20_u + mat[0][0][1] * bs21_u + mat[0][0][2] * bs22_u) * bs20_t
      + (mat[0][1][0] * bs20_u + mat[0][1][1] * bs21_u + mat[0][1][2] * bs22_u) * bs21_t
      + (mat[0][2][0] * bs20_u + mat[0][2][1] * bs21_u + mat[0][2][2] * bs22_u) * bs22_t)
      * bs20_s
      +
       ((mat[1][0][0] * bs20_u + mat[1][0][1] * bs21_u + mat[1][0][2] * bs22_u) * bs20_t
      + (mat[1][1][0] * bs20_u + mat[1][1][1] * bs21_u + mat[1][1][2] * bs22_u) * bs21_t
      + (mat[1][2][0] * bs20_u + mat[1][2][1] * bs21_u + mat[1][2][2] * bs22_u) * bs22_t)
      * bs21_s
      +
       ((mat[2][0][0] * bs20_u + mat[2][0][1] * bs21_u + mat[2][0][2] * bs22_u) * bs20_t
      + (mat[2][1][0] * bs20_u + mat[2][1][1] * bs21_u + mat[2][1][2] * bs22_u) * bs21_t
      + (mat[2][2][0] * bs20_u + mat[2][2][1] * bs21_u + mat[2][2][2] * bs22_u) * bs22_t)
      * bs22_s;
}

Answer 1

template<std::size_t...Is>
auto index_over(std::index_sequence<Is...>) {
  return [](auto&& f)->decltype(auto){
    return decltype(f)(f)( std::integral_constant< std::size_t, Is >{}... );
  };
}
template<std::size_t N>
auto index_upto(std::integral_constant<std::size_t, N> ={}) {
  return index_over( std::make_index_sequence<N>{} );
}
inline float sum() { return 0.0f; }
template<class...Args>
float sum( float a, Args... args ) {
  return a + sum(args...);
}

, .

auto Z_f = [&](auto X, auto Y)->float {
  return index_upto<2>()( [&](auto...Zs)->float{
    return sum((poly<2,Zs>(fU) * mat[X][Y][Zs])...);
  });
};
auto Y_f = [&](auto X)->float {
  return index_upto<2>()( [&](auto...Ys)->float{
    return sum( (poly<2,Ys>(fT) * Z_f(X, Ys))... );
  });
};
auto X_f = [&]()->float {
  return index_upto<2>()( [&](auto...Xs)->float{
    return sum( (poly<2,Xs>(fS) * Y_f(Xs))... );
  });
};

float val = X_f();

, - , , , X_f, Y_f Z_f .

clang .

++ 14 (index_sequence make_index_sequence), ++ 11.

auto lambdas, ++ 14. ++ 11 , , .

sum (0.f + ... + args) ++ 17.

, .

godbolt, , .

, , 3 Xs, Ys Zs, .

auto contribution = [&](auto X, auto Y, auto Z) {
  return mat[X][Y][Z] * poly<2,X>(fS) * poly2<2,Y>(fT) * poly2<2,Z>(fU);
};

, 3 , .

auto summer_1d = [](auto...Vals)->decltype(auto){
  return sum(Vals...);
};
template<std::size_t X_max, std::size_t Y_max, std::size_t Z_max, class Sum = decltype(summer_1d)>
auto sumup_3d(Sum sum = summer_1d) {
  return [](auto&& f)->decltype(auto) {
    auto Z_part = [&](auto X, auto Y)->decltype(auto) {
      return index_upto<Z_max>()([&](auto...Zs)->decltype(auto){
        return sum( f(X,Y,Zs)... );
      });
    };
    auto Y_part = [&](auto X)->decltype(auto) {
      return index_upto<Y_max>()([&](auto...Ys)->decltype(auto){
        return sum( Z_part(X, Ys)... );
      });
    };
    return index_upto<X_max>()([&](auto...Xs)->decltype(auto){
      return sum( Y_part(Xs)... );
    });
  };
};

auto val = sumup_3d<3,3,3>()(contribution);

somesuch.