std::mdspan×空間充填曲線

C++2b(C++23)標準ライブラリの多次元ビューstd::mdspanクラステンプレート（→id:yohhoy:20230303）は、第3テンプレートパラメータLayoutPolicyを介して任意のメモリレイアウトをサポートする。*1

C++2b標準ライブラリは、多次元配列表現としてメジャーな行優先(row major)／列優先(column major)レイアウトstd::layout_right／std::layout_left、これらを汎化したストライドレイアウトstd::layout_strideのみを提供する。一方でstd::mdspanに指定するレイアウトポリシーの表現能力としては、タイル化レイアウト(tiled layout)や空間充填曲線(space-filling curve)といった複雑なメモリレイアウトも表現できる。

空間充填曲線の一種であるヒルベルト曲線(Hilbert curve)の実装例：*2

// C++2b(C++23)
#include <bit>
#include <mdspan>
#include <utility>

struct layout_hilbert {
  template <class Extents>
  class mapping {
  public:
    using extents_type = Extents;
    using index_type = typename extents_type::index_type;
    using rank_type = typename extents_type::rank_type;
    using layout_type = layout_hilbert;

    mapping(const extents_type& e)
      : extents_{e}
    {
      static_assert(extents_type::rank() == 2);
      assert(e.extent(0) == e.extent(1));
      assert(std::has_single_bit(e.extent(0)));
    }
    constexpr const extents_type& extents() const noexcept
      { return extents_; }
    template <class I0, class I1>
    constexpr index_type operator()(I0 i, I1 j) const noexcept
    {
      const index_type N = std::bit_ceil(extents_.extent(0));
      return xy2d(N, j, i);  // x,y=j,i
    }
    constexpr index_type required_span_size() const noexcept
    {
      const index_type N = std::bit_ceil(extents_.extent(0));
      return N * N;
    }
    static constexpr bool is_always_unique() noexcept { return true; }
    static constexpr bool is_always_exhaustive() noexcept { return true; }
    static constexpr bool is_always_strided() noexcept { return false; }
    static constexpr bool is_unique() noexcept { return true; }
    static constexpr bool is_exhaustive() noexcept { return true; }
    static constexpr bool is_strided() noexcept { return false; }
  private:
    // https://en.wikipedia.org/wiki/Hilbert_curve より
    static constexpr index_type xy2d(
      index_type n, index_type x, index_type y)
    {
      index_type rx, ry, s, d = 0;
      for (s = n / 2; s > 0; s /= 2) {
        rx = (x & s) > 0;
        ry = (y & s) > 0;
        d += s * s * ((3 * rx) ^ ry);
        rot(n, &x, &y, rx, ry);
      }
      return d;
    }
    static constexpr void rot(
      index_type n, index_type *x, index_type *y,
      index_type rx, index_type ry)
    {
      if (ry == 0) {
        if (rx == 1) {
          *x = n-1 - *x;
          *y = n-1 - *y;
        }
        std::swap(*x, *y);
      }
    }
  private:
    extents_type extents_{};
  };
};

int arr[16] = /*...*/;

// 4x4=16要素がヒルベルト曲線順に配置された2次元ビュー
using MatExts = std::dextents<size_t, 2>;
std::mdspan mat{arr, layout_hilbert::mapping{MatExts{4, 4}}};
// mat[i,j]⇒arr要素順
//   j | 0  1  2  3
// ----+------------
// i=0 | 1  2 15 16
//   1 | 4  3 14 13
//   2 | 5  8  9 12
//   3 | 6  7 10 11

layout_hilbert::mappingクラステンプレートによって、多次元インデクス(i,j)からメモリ領域へのオフセット位置への変換(operator())を定義している。変数mapはCTAD（→id:yohhoy:20230308）によりstd::mdspan<int, std::dextents<size_t, 2>, layout_hilbert>型*3に推論される。

yohhoyの日記

技術的メモをしていきたい日記

std::mdspan×空間充填曲線