生成编译时数组结构 (c++17)

Generating compile-time array structure (c++17)

我有以下 C++17 代码,用于生成数组的编译时元组,其中零数组只是为了示例,在我的实现中它们将是完整的(使用 -std= 编译c++1z -fconcepts).

#include <array>
#include <tuple>
#include <cmath>

template <std::size_t nrA, std::size_t ncA, std::size_t nrB, std::size_t ncB,
          typename number=double>
constexpr auto operator *(const std::array<std::array<number,ncA>,nrA> & A,
                          const std::array<std::array<number,ncB>,nrB> & B)
{
  std::array<std::array<number,ncB>,nrA> res{} ;
  for (auto k=0u;k<ncB;++k)
    for (auto i=0u;i<nrA;++i)
      for (auto j=0u;j<nrB;++j)
        res[i][k] += A[i][j]*B[j][k];
  return res ;
}

constexpr auto logg2(const auto N)
{
  auto res = 0;
  auto n = N;
  while (n != 0) 
  {
    n /= 2;
    ++res;
  }
  return res;
}

template <std::size_t N,typename number=double>
constexpr auto create_R()
{
  return std::array<std::array<double,2*N>,N>{};
}

template <std::size_t N,typename number=double>
constexpr auto create_RT()
{
  return std::array<std::array<double,N>,2*N>{};
}

template <std::size_t N,std::size_t ...Is>
constexpr auto make_impl(const std::index_sequence<Is...>)
{
  return std::make_tuple(std::make_tuple(create_R<(N >> Is)>()...),
                         std::make_tuple(create_RT<(N >> Is)>()...));
}

template <std::size_t N,typename number=double>
constexpr auto make()
{
  return make_impl<N/2>(std::make_index_sequence<logg2(N/2) - 1>());
}

int main(int argc, char *argv[])
{
  const auto n = 4u;
  const auto A = std::array<std::array<double,2*n>,2*n>{};
  const auto [R,RT] = make<2*n>();
}

我想将 make<>() 修改为 make<>(A),并且 return 一个结构化绑定 [R,RT,As] 其中 As 是一个包含以下数组的元组它

                              A,
               std::get<0>(R)*A*std::get<0>(RT),
std::get<1>(R)*std::get<0>(R)*A*std::get<0>(RT)*std::get<1>(RT)
                             ...

我已经尝试了一段时间,没有找到解决方案。

有什么想法吗?

编辑 1

根据@MaxLanghof 的要求,以下打印矩阵:

template <std::size_t nr, std::size_t nc, typename number=double>
constexpr auto print(const std::array<std::array<number,nc>,nr> & A)
{
  for (auto i=0u;i<nr;++i)
    {
      for (auto j=0u;j<nc;++j)
        std::cout << std::right << std::setw(12) << A[i][j];
      std::cout << std::endl ;
    }
  std::cout << std::endl ;
}

并将以下行添加到 main()

print(A);
print(std::get<0>(R)*A*std::get<0>(RT));
print(std::get<1>(R)*std::get<0>(R)*A*std::get<0>(RT)*std::get<1>(RT));

得到如下输出

       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0
       0           0           0           0           0           0           0           0

       0           0           0           0
       0           0           0           0
       0           0           0           0
       0           0           0           0

       0           0
       0           0

解决方案涉及为 As 元组的每个条目创建另一个索引序列,然后使用它来折叠表示乘法。为了可读性,我冒昧地将 std::array<std::array<T, Cols>, Rows> 包装在一个类型中(这最终也是必要的,见下文)。每次调用 makeMul 都会产生一个 As 元组元素(原来的 A 是单独添加的)。

template <std::size_t Rows, std::size_t Cols, typename T = double>
struct Mat2D : std::array<std::array<T, Cols>, Rows> {};

template <class T_Rs, std::size_t... Is>
constexpr auto mult_lhs(T_Rs Rs, std::index_sequence<Is...>) {
  return (std::get<sizeof...(Is) - Is - 1>(Rs) * ...);
}

template <class T_RTs, std::size_t... Is>
constexpr auto mult_rhs(T_RTs RTs, std::index_sequence<Is...>) {
  return (std::get<Is>(RTs) * ...);
}

template <class T_A, class T_Rs, class T_RTs, std::size_t... Is>
constexpr auto makeMul_impl(T_A A, T_Rs Rs, T_RTs RTs,
                            std::index_sequence<Is...> is) {
  return mult_lhs(Rs, is) * A * mult_rhs(RTs, is);
}

template <std::size_t Index, class T_A, class T_Rs, class T_RTs>
constexpr auto makeMul(T_A A, T_Rs Rs, T_RTs RTs) {
  return makeMul_impl(A, Rs, RTs, std::make_index_sequence<Index + 1>());
}

template <std::size_t N, std::size_t... Is, typename T = double>
constexpr auto make_impl(const Mat2D<2 * N, 2 * N, T>& A,
                         std::index_sequence<Is...>) {
  auto Rs = std::make_tuple(create_R<(N >> Is)>()...);
  auto RTs = std::make_tuple(create_RT<(N >> Is)>()...);
  auto As = std::make_tuple(A, makeMul<Is>(A, Rs, RTs)...);
  return std::make_tuple(Rs, RTs, As);
}

template <std::size_t N, typename T = double>
constexpr auto make(const Mat2D<N, N, T>& A) {
  return make_impl<N / 2>(A, std::make_index_sequence<logg2(N / 2) - 1>());
}

int main(int argc, char* argv[]) {
  const auto n = 4u;
  const auto A = Mat2D<2 * n, 2 * n, double>{};
  const auto [Rs, RTs, As] = make(A);
}

Demo

重要的是要注意 std 类型的重载运算符在这里是一个问题,至少对于 clang(更严格地遵循标准):一旦您尝试使用重载的 operator* 在模板中,它不会被发现,因为 ADL 在 namespace std 中查找它(我最初将 Mat2D 作为别名而不是 struct 继承的东西) - 和您不能向 namespace std 添加内容(一些特定的自定义点除外)。至少我是这样理解的 this error.

std 类型继承总体来说非常糟糕,但我假设您的矩阵实际上是用户定义的类型,因此 none 这很重要。

最后,我强烈建议您为元组类型命名。当你有元组的元组(数组的数组)时,每个 reader 都必须投入大量时间甚至掌握代码。如果你例如,它已经有帮助了。将每个 RRT 分组在一个结构中:

template<std::size_t N, typename T = double>
struct R_RT {
  Mat2D<N, 2 * N, T> R;
  Mat2D<2 * N, N, T> RT;
};

and/or 例如一个

template<class TupleOfR, class TupleOfRT, class TupleOfAs>
struct Rs_RTs_As {
  TupleOfR Rs;
  TupleOfRT RTs;
  TupleOfAs As;
};

即使这在技术上允许任何三种类型,它仍然记录了您应该期望在那里找到的内容,并且通过概念,您实际上可以适当地限制所有内容(包括诸如“As 将具有的乐趣比 RsRTs" 多一个元素",这可能需要大多数 reader 一段时间才能从纯元组代码中实现。