頂点数と次数を指定して正則グラフを生成する

概要

AtCoder Heuristic Contest 016 に参加したときに、頂点数  n と次数  d とを任意に指定して正則グラフを生成するアルゴリズムを調べて実装しました。せっかくなので、それについて書きます。生成されるグラフのランダム性などは考えず、とりあえず所望の正則グラフがひとつ得られたら満足するものとします。この記事は文献 https://ipsj.ixsq.nii.ac.jp/ej/action=repository_uri&item_id=32553&file_id=1&file_no=1 で紹介されているアルゴリズムの、目的に必要な部分だけを紹介するものですが、正当性の証明なども補足したいと思います。正則グラフの生成アルゴリズムの研究に関してしっかり調べたわけではないので、もっとよいアルゴリズムがあるかもしれません。

アルゴリズム

 n 頂点の  d -正則グラフを生成するアルゴリズムです。ただし、 d\lt n かつ  nd は偶数とします
(各辺は2回ずつ次数の総和  nd に寄与するので、偶数でなければならないことがわかります)。逆に、この条件を満たせば正則グラフが構成できることが以下のアルゴリズムにより保証されます。

手順1

完全グラフ  K_{d+1} を、頂点数が超過しない限り並べる。最大で  \lfloor{\frac{n}{d+1}}\rfloor 個並べることができる。こうしてできるグラフは、 d -正則グラフである。

手順2

 d -正則グラフであることを保ったまま頂点を増やしていき、目的の頂点数  n d -正則グラフを得る。 d の偶奇によって操作が異なる。

手順2-a ( d が偶数のとき)

頂点を  1 つずつ増やしていく。
 d/2 本の互いに隣接しない辺を選ぶ。辺の端点すべてに対して、新たに加える頂点  v との間に辺を張る。そして、選んだ辺をすべて取り除く。

手順2-b ( d が奇数のとき)

頂点を  2 つずつ増やしていく(頂点数が偶数でなければならないので、 1 つだけ増やすことは不可能)。
長さ  d の単純パス  (x_0,x_1,\dots,x_d) をとる。新たに加える頂点を  v_1,v_2 として、 v_1 x_0,\dots,x_{d-1} に、 v_2 x_1,\dots,x_d に接続する。そして、パス上の辺をすべて取り除く。


手順2に関しては、文献にわかりやすい図が載っているので参照してください。

アルゴリズムの正当性

 d/2 本の互いに隣接しない辺を選ぶには( d は偶数)

まだ辺の端点として選ばれていない頂点の集合が誘導する部分グラフから、任意に1つの辺を選ぶことを繰り返せばよいです。これが必ず可能なことを示します。
(証明)まだ辺の端点として選ばれていない頂点の集合を  S とする。 v\in S を任意にとる。このとき、ある  u\in S が存在して 辺  (v,u) が存在することを示す。このような  u が存在しない状況を考える。 v が隣接する頂点がすべて既に辺の端点として選ばれていることになるが、 v の次数は  d なので、既に  d/2 本の互いに隣接しない辺が選ばれていることになる。 □

長さ  d の単純パスをとるには( d は奇数)

はじめに任意に頂点を選び、1頂点からなるパスとします。そのあと、現時点でのパスの終点と、パスに含まれない頂点との間に辺を張ることを繰り返せばよいです。これが必ず可能なことを示します。
(証明)各時点において、パスに含まれない頂点であって、パスの終点との間に辺が張られたものが存在することを示す。このような頂点が存在しない状況を考える。パスの終点を  v とする。 v に隣接する頂点がすべてパスに含まれることになるが、 v の次数は  d なので、既に  d+1 頂点からなる単純パスが得られていることになる。 □

実装例

1000頂点くらいまでなら高速に動作します。

#include <vector>

// n頂点のd-正則グラフを生成し、その隣接行列を返す
// ただし、d*nは偶数でd<nとする
std::vector<std::vector<bool>> generate_regular_graph(int n, int d) {
    std::vector<std::vector<bool>> graph(n, std::vector<bool>(n, false));

    // d+1頂点の完全グラフを可能な限り作る
    const int num{n / (d + 1)};
    for (int i = 0; i < num; ++i) {
        const int begin{i * (d + 1)};
        const int end{(i + 1) * (d + 1)};
        for (int a = begin; a < end; ++a) {
            for (int b = a + 1; b < end; ++b) {
                graph[a][b] = graph[b][a] = true;
            }
        }
    }

    // d-正則グラフのまま頂点数を増やしていく
    int now_n{num * (d + 1)};
    if (d % 2 == 0) {  // dが偶数のときは1頂点ずつ増やす
        std::vector<int> from(d / 2), to(d / 2);
        std::vector<bool> used(n, false);
        while (now_n < n) {
            std::fill(std::begin(used), std::end(used), false);

            // d/2本の互いに隣接しない辺を選ぶ
            int a{0};
            for (int i = 0; i < d / 2; ++i) {
                while (used[a]) ++a;
                for (int b = a + 1; b < now_n; ++b) {
                    if (used[b]) continue;
                    if (graph[a][b]) {
                        from[i] = a, to[i] = b;
                        used[a] = used[b] = true;
                        break;
                    }
                }
            }

            // 頂点を増やしてグラフを組み替える
            for (int i = 0; i < d / 2; ++i) {
                const int a{from[i]}, b{to[i]};
                graph[a][b] = graph[b][a] = false;
                graph[a][now_n] = graph[now_n][a] = true;
                graph[b][now_n] = graph[now_n][b] = true;
            }
            ++now_n;
        }
    } else {  // dが奇数のときは2頂点ずつ増やす
        std::vector<int> path;
        std::vector<bool> used(n, false);
        while (now_n + 1 < n) {
            std::fill(std::begin(used), std::end(used), false);
            path = {0};  // 頂点0から始める
            used[0] = true;

            // 長さdの単純パスを選ぶ
            for (int i = 0; i < d; ++i) {
                const int from{path.back()};
                for (int to = 0; to < now_n; ++to) {
                    if (used[to] || !graph[from][to]) continue;
                    path.emplace_back(to);
                    used[to] = true;
                    break;
                }
            }

            // 頂点を増やしてグラフを組み替える
            for (int i = 0; i < d; ++i) {
                const int a{path[i]}, b{path[i + 1]};
                graph[a][b] = graph[b][a] = false;
            }
            for (int i = 0; i < d; ++i) {
                graph[path[i]][now_n] = graph[now_n][path[i]] = true;
            }
            for (int i = 1; i < d + 1; ++i) {
                graph[path[i]][now_n + 1] = graph[now_n + 1][path[i]] = true;
            }
            now_n += 2;
        }
    }

    return graph;
}

感想

思ってたよりも簡単なアルゴリズムで正則グラフが生成できて驚き。

XORベクトル空間の基底の数え上げ

 \mathbb{F}_2 = \{0,1\} 上のベクトル空間  \mathbb{F}_2^n について。
線形独立な集合  S\subset\mathbb{F}_2^n とする。このとき、任意の  x\not\in\mathrm{span}\;S に対して  S\cup\{x\} もまた線形独立になる。ただし、 \mathrm{span}\;\emptyset = \{0\} である。また、線形独立な集合の部分集合もまた線形独立である。線形代数で基本的なことだけど、これらにより次がわかる。

  • 線形独立で大きさ  k\in\mathbb{N} の集合の個数は、 \prod_{i=0}^{k-1}{(2^n-2^i)}\;/k!
  • 特に、基底(順序を考えない)の個数は、 \prod_{i=0}^{n-1}{(2^n-2^i)}\;/n!

関連

A053601 - OEIS
ARC146 C - Even XOR (この問題をきっかけに記事を書いた)

感想

ちょっと前にも少し考えていたのだけど、そのときは思いつかなかった。サイズを  1 ずつ大きくすることを考えたらすぐなのだけど。(そしてこのブログ、やけにこの手の話題が好きだな...)

ARC145 C - Split and Maximize 解法の正当性

ARC145 C - Split and Maximize公式解説があっさりしていて説明不足だと思ったので、行間を埋めていたら結構たいへんだった。

スコアが最大となる構造の分析

 A\lt B\lt C\lt D のとき、 AD+BC\lt AC+BD\lt AB+CD が成り立つのは公式解説の通りです。
 1 2N を使った  N ペアの作り方が任意に与えられたとします。ペア  \{1,X\} \{2,Y\} があるとき、 \{1, 2\} \{X, Y\} にペアを組み替えます。するとスコアを減少させずに  \{1,2\} のペアが作れます。続けて  3 4 5 6 \cdots 2N-1 2N に注目して同じことを繰り返すと、スコアを減少させずに  \{1, 2\}, \{3, 4\}, \dots, \{2N-1, 2N\} というペアに組み替えることができます。任意のペアの作り方についてこれが言えるので、 \{1, 2\}, \{3, 4\}, \dots, \{2N-1, 2N\} というペアを作ったときが最大スコアとなります。

括弧列との対応付け

ここが説明不足だと感じました。
スコア  M を達成する順列  P をとります。そして、そのスコア  M を達成する  A,B の分割をひとつとります。 P の先頭から  A に振り分けられたものを ( と書き、 B に振り分けられたものを ) と書いて並べることで、ひとつの括弧列を得ます。ここで、この括弧列はバランスがとれているとは限らないことに注意してください。
この括弧列を先頭から見ていって、初めてバランスが壊れる)の場所を見つけ、それ以降の()を反転させます。これは  A に振り分けものと  B に振り分けるものを交換することを意味し、操作後も分割のスコアは  M のままです。この操作を括弧列のバランスがとれるまで繰り返すことで、スコア  M を達成する  A,B の分割であって、対応する括弧列のバランスがとれているものが得られます。
つまり、スコア  M の順列があったとき、スコア  Mを達成し、対応する括弧列がバランスのとれているような分割がただひとつ存在します。よって、バランスのとれた括弧列すべてに対して、対応する分割のスコアが  M になるような順列を数え上げればよく、あとは公式解説の通りです。

感想

解法を考えるより、解法の正当性を示すほうが楽しい(ときもある)。

ダブリングの抽象化

競プロでダブリングと呼ばれるアルゴリズムの抽象化を考えてみました。ダブリングのアルゴリズム自体の解説はしません。

扱う問題

有限の状態集合  T と、半群  (S, \cdot) を考えます。2つの写像  t:T\to T v: T\to S とが与えられたとき、任意の  s\in T k\in\mathbb{N} に対して、 v(s)\cdot v(t(s))\cdot v(t^2(s))\cdot\cdots\cdot v(t^{k}(s)) を求めるという問題を扱います。
簡潔にいうと、状態遷移しながら状態がもつ値を拾っていってその積を求めるということです。半群というのは集合の上に二項演算が定義されていて、その演算が結合律を満たすものをいいます。これは、結合律を満たさないとダブリングができないためです。
また、この問題において  |T|=1 の場合が繰り返し二乗法であると見做すこともできます。

実装

template <typename S, int POW_MAX>
struct Doubling {
   private:
    int n;
    vector<int> to;
    vector<S> value;
    function<S(S, S)> op;
    vector<vector<int>> next;
    vector<vector<S>> dp;

    void init() {
        next.resize(POW_MAX, vector<int>(n));
        dp.resize(POW_MAX, vector<S>(n));
        for (int start = 0; start < n; ++start) {
            next[0][start] = to[start];
            dp[0][start] = value[to[start]];
        }
        for (int pow = 1; pow < POW_MAX; ++pow) {
            for (int start = 0; start < n; ++start) {
                next[pow][start] = next[pow - 1][next[pow - 1][start]];
                dp[pow][start] =
                    op(dp[pow - 1][start], dp[pow - 1][next[pow - 1][start]]);
            }
        }
    }

   public:
    Doubling(int n, const vector<int>& to, const vector<S>& value,
             function<S(S, S)> op)
        : n(n), to(to), value(value), op(op) {
        init();
    }
    S prod(int start, long long step) const {
        S prod{value[start]};
        int now{start};
        for (int pow = 0; pow < POW_MAX; ++pow) {
            if (step & (1LL << pow)) {
                prod = op(prod, dp[pow][now]);
                now = next[pow][now];
            }
        }
        return prod;
    }
};

POW_MAXというのはprodに投げるstepの最大値を2冪で指定します。例えば  10^{18} まで取りうるならば60とすればよいです。ダブリングの配列の添字を逆にしたい気持ちになったのですが、添字の順番を逆にすると激遅になります。キャッシュを意識するのは大事ですね。


抽象化パワーのお手並み拝見といきましょう

ABC167 D - Teleporter

ダブリングが使える基本問題ですね。 今回の問題設定に帰着させるには、

  • 状態集合は  [N]=\{1,2,\dots,N\}
  • 半群の台集合は  [N]、演算は  a\cdot b = b
  •  t:[N]\to[N] t(i) = A_i
  •  v:[N]\to[N] v(i) = i

とすればよいです。この演算が結合律 (a\cdot b)\cdot c = a\cdot (b\cdot c) を満たすのは簡単に確認できます。この問題では状態の遷移先だけ考えればいいのですが、冗長になってしまいました。

int main() {
    int N;
    long long K;
    cin >> N >> K;
    vector<int> to(N), value(N);
    for (auto& A : to) {
        cin >> A;
        --A;
    }
    iota(begin(value), end(value), 1);
    auto op = [](int a, int b) { return b; };
    Doubling<int, 60> doubling(N, to, value, op);
    cout << doubling.prod(0, K) << '\n';
}

ABC179 E - Sequence Sum

これは今回の問題設定そのまんまという感じです。 [M]_0=\{0,1,2,\dots,M-1\} と書くことにします。

  • 状態集合は  [M]_0
  • 半群の台集合は  \mathbb{N}、演算は  a\cdot b = a + b \mathbb{N} における標準的な和)
  •  t:[M]_0\to[M]_0 t(i) = i^2 \mathrm{mod} M
  •  v:[M]_0\to \mathbb{N} v(i) = i

とすればよいです。

int main() {
    long long n;
    int x, m;
    cin >> n >> x >> m;
    vector<int> to(m);
    for (long long from = 0; from < m; ++from) {
        to[from] = from * from % m;
    }
    vector<long long> value(m);
    iota(begin(value), end(value), 0LL);
    auto op = [](long long a, long long b) { return a + b; };
    Doubling<long long, 34> doubling(m, to, value, op);
    cout << doubling.prod(x, n - 1) << '\n';
}

ABC175 D - Moving Piece

実装でバグり散らかしそうな問題です。今回の問題設定に帰着させるには、

  • 状態集合は  [N]=\{1,2,\dots,N\}
  • 半群の台集合は  \mathbb{Z}^2、演算は  (a_1,a_2)\cdot(b_1,b_2) = (a_1+b_1,\max(a_2,a_1+b_2))
  •  t:[N]\to[N] t(i)=P_i
  •  v:[N]\to\mathbb{Z}^2 v(i) = (C_i,C_i)

とすればよいです。 (a_1,a_2)\in\mathbb{Z}^2の気持ちとしては、 a_1 がスコアの累積和、 a_2がスコアの累積和の累積最大値を表しています。演算が結合律を満たすことを確かめましょう。

 \begin{align}
( (a_1,a_2)\cdot(b_1,b_2))\cdot(c_1,c_2) &= (a_1+b_1,\max(a_2,a_1+b_2))\cdot(c_1,c_2) \\
&= (a_1+b_1+c_1,\max(\max(a_2,a_1+b_2),a_1+b_1+c_2)) \\
&= (a_1+b_1+c_1,\max(a_2,a_1+b_2,a_1+b_1+c_2))
\end{align}
 \begin{align}

(a_1,a_2)\cdot( (b_1,b_2)\cdot(c_1,c_2)) &= (a_1,a_2)\cdot(b_1+c_1,\max(b_2,b_1+c_2)) \\
&= (a_1+b_1+c_1,\max(a_2,a_1+\max(b_2,b_1+c_2)) \\
&= (a_1+b_1+c_1,\max(a_2,\max(a_1+b_2,a_1+b_1+c_2))) \\
&= (a_1+b_1+c_1,\max(a_2,a_1+b_2,a_1+b_1+c_2))
\end{align}

int main() {
    int N, K;
    cin >> N >> K;
    vector<int> to(N);
    for (auto& P : to) {
        cin >> P;
        --P;
    }
    using S = pair<long long, long long>;
    vector<S> value(N);
    for (auto& val : value) {
        long long C;
        cin >> C;
        val = {C, C};
    }
    auto op = [](S a, S b) -> S {
        return {a.first + b.first, max(a.second, a.first + b.second)};
    };
    Doubling<S, 30> doubling(N, to, value, op);
    long long ans{-(1LL << 62)};
    for (int i = 0; i < N; ++i) {
        ans = max(ans, doubling.prod(i, K - 1).second);
    }
    cout << ans << '\n';
}


ほかに適用できる問題を見つけたら随時追加していきます。


感想

抽象化したデータ構造やアルゴリズムを使って問題を解くと見通しはよくなりますが、持つべきデータと演算を定めるときに工夫が必要なときがあり、往々にしてパズルをするはめになるのがつらそうです。

ARC138 D - Differ by K bits の解法の正当性を示す

なんの記事?

N ビットグレイコードの、ハミング距離 K のバージョンは存在しますか?存在するならばそのひとつを構築してください、という問題がARC138 Dにて出題されました。

K=1 のときはグレイコードそのものなので、必ず構築できます。以下では、K\geq2 のときを考えることにします。所望の順列が存在するための自明な必要条件として、次の条件が考えられます。

  • K が奇数であり、かつ K\lt N

逆にこの条件を満たすならばいつでも構築可能なのですが、公式解説ではそれを示す過程で、

N bit の XOR の基底であって,popcount が
K になる値のみからなるものを取ることができます

という事実を使っています。これは(たぶん)非自明だと思うので証明を考えてみました、というのが本記事の内容です。

示したいこと

以下では単に「行列」と言ったら、体  \mathbb{F}_2=\{0,1\} 上の行列を考えるものとします。示すべきことは、「正整数 N と、K\lt N なる奇数K とを任意にとったとき、各行に 1 がちょうど K 個あるような N\times N 行列 M が存在して、M に対して行基本変形を施すことにより単位行列にできる」ということです。行基本変形は可逆なので、変形していく方向を逆にした、以下の命題が証明できればよいです。

命題

正整数 N と、K\lt N なる奇数 K とを任意にとる。このとき、N\times N 単位行列行基本変形を施すことにより、各行に 1 がちょうど K 個あるような行列にできる。

証明

行基本変形の手順を以下のように与えます。

  1. N\times N 単位行列からスタートする。K=1 ならば既に目標の行列になっている。
  2. 3,\dots,K+1 行目を 1 行目と 2 行目に足す。
  3. 1 行目と 2 行目を 3,\dots,N 行目に足す。
  4. ここで、右下の (N-2)\times(N-2) の領域に注目すると、(N',K')=(N-2,K-2) の部分問題が現れている(N',K' の条件は保たれる)。よって、再帰的に目標の行列を構成できる。

言葉だとわかりにくいので、(N,K)=(7,5) のときを例にとって、変形されていく様子を図解します。

matrix1

matrix2

matrix3

matrix4

matrix5

matrix6

matrix7

最後は、部分問題が (N,K)=(3,1) となって変形が終了しました。図で灰色になっている部分は、その時点で既に確定した部分です。既に確定したはずの部分が、後の変形で壊れたりすることがあるのではないかと思うかもしれませんが、K が奇数であることなどより、そのようなことは起こりえません(例を追ってみるとわかると思います)。だから安心して部分問題に帰着できるわけです。

おわりに

というわけで、証明できました。証明さえできれば、今回は N が小さいので、基底の構成は愚直に O(N2^N) で十分です。今回の証明での構成方法であれば、(基底の構成パートは)O(N^2) でできそうです。
もっと簡単な証明があったり、もっと綺麗な基底があったりしたら教えてください。

XORベクトル空間のアルゴリズムとその正当性

はじめに

自然数全体の集合  \mathbb{N}=\{0,1,2,\dots\} は、以下のように和とスカラー倍を定めることにより、体  \mathbb{F}_2=\{0,1\} 上のベクトル空間と見なせます。このベクトル空間を  \mathbb{N}_\mathrm{XOR} と書きます。

  • 和は、bitwiseのXOR演算(以下では単に  + と書く)
  • スカラー倍は、 0a=0,\;1a=a

本記事では、 \mathbb{N}_\mathrm{XOR} の部分空間の基底を管理するのに便利なアルゴリズムと、その正当性についてまとめています。

アルゴリズムの実装

struct XorVecSpace {
    // 部分空間の基底を蓄える配列
    vector<int> basis;

    // 部分空間の生成元を追加、basisのサイズ(部分空間の次元)が大きくなったらtrueを返す
    bool add(int v) {
        int mn = minimize(v);
        if (mn) basis.push_back(mn);
        return mn != 0;
    }

    // 部分空間の元であるかを判定
    bool is_element(int v) const { return minimize(v) == 0; }

    // 部分空間の元とvとの和の最小値を計算
    int minimize(int v) const {
        for (int eb : basis) {
            v = min(v, v ^ eb);
        }
        return v;
    }

    // 部分空間の元とvとの和の最大値を計算
    int maximize(int v) const {
        for (int eb : basis) {
            v = max(v, v ^ eb);
        }
        return v;
    }

    // 部分空間の次元を返す
    int dim() const { return basis.size(); }
};

以下では、自然数の集合  S\subset\mathbb{N} が生成する  \mathbb{N}_\mathrm{XOR} の部分空間を  \mathrm{span}(S) と書きます。

  • 配列 basis(以下では  b と書くことにします) に現在管理している部分空間の基底が蓄えられます。
  • add:部分空間の生成元として  v を追加し、次元が大きくなったか否かを返します。 このとき、渡した元そのものが  b に入るわけではないので、未加工の元からなる基底が欲しい場合は別途管理する必要があります。
  • is_element v\in\mathrm{span}(b) であるかを判定します。
  • minimize:集合  \{v+w\;|\;w\in\mathrm{span}(b)\} の最小元を返します。
  • maximize:集合  \{v+w\;|\;w\in\mathrm{span}(b)\} の最大元を返します。
  • dim:部分空間の次元を返します。

アルゴリズムの正当性

まず、add関数の正当性を示すことを考えます。add関数でやっていることは、集合  \{v+w\;|\;w\in\mathrm{span}(b)\} の最小元  m を求め、 m\neq 0 ならば  m\not\in\mathrm{span}(b) なので  b m を追加する、ということです。よって、minimize関数の正当性を示せばよいです。

 b=(b_1,b_2,\dots,b_n) とし、 b_i の最上位ビットを  d_i 桁目とします。minimize関数でやっていることは、 i=1,2,\dots,n の順に、以下で定義される操作  \mathrm{op}_i v に施すということです。

操作 \mathrm{op}_i v d_i 桁目が  1 ならば  v b_i を足す。

このアルゴリズムで集合  \{v+w\;|\;w\in\mathrm{span}(b)\} の最小元を求められるのは、add関数によって蓄えられる基底  b がもつ以下の性質によります。
性質

 b_j において  d_i\;(i\lt j) 桁目は  0 である。

(性質の証明)基底  b のサイズ  n に関する帰納法による。 b のサイズが  n-1 のとき性質が成り立つことを仮定する。add関数で  b b_n が追加されたとき、 b_n d_i\;(i=1,2,\dots,n-1) 桁目は  0 であることを示せばよい。minimize関数において、操作  \mathrm{op}_i の直後、 v d_i 桁目は  0 であるが、帰納法の仮定により、以降の操作で  d_i 桁目に干渉することはないため、 b_n においても  d_i 桁目は  0 になる。 □

さて、この性質を使うことでminimize関数の正当性を示すことができます。minimize関数でやりたいことは、 i=1,2,\dots,n の順に、 v b_i を足すor足さないの操作をすることで、 v の最終的な値を最小化するということです。性質より、 v d_i 桁目は  i 番目よりあとの操作で不変です。つまり、 i 番目の操作によって  v d_i 桁目の最終的な値を決めることができます。また、 i 番目の操作では  v d_i 桁目より上の桁には干渉しません。 これらを踏まえると、v の最終的な値を最小化するには、操作  \mathrm{op}_iのようにすればよいです。したがって、minimize関数の正当性がわかりました。

maximize関数の正当性もminimize関数の正当性と同様に理解することができます。

おまけ

bitsetバージョンの実装例です。bitsetのmsbを求める関数がないのはなぜでしょうか?std::bitset::_Find_firstでlsbは求まるので、次元を求めるだけならこれを使って実装してもよいと思います。ただし、最小化、最大化などは壊れてしまいます。

template <int N>
struct XorVecSpace {
    using bs = bitset<N>;

    // 部分空間の基底を蓄える配列
    vector<bs> basis;
    // 基底のmsbをもっておく
    vector<size_t> msb;

    // 部分空間の生成元を追加、basisのサイズ(部分空間の次元)が大きくなったらtrueを返す
    bool add(const bs& v) {
        bs mn{minimize(v)};
        if (mn.none()) return false;
        basis.emplace_back(mn);
        msb.emplace_back(N - 1);
        while (!basis.back().test(msb.back())) --msb.back();
        return true;
    }

    // 部分空間の元であるかを判定
    bool is_element(const bs& v) const { return minimize(v).none(); }

    // 部分空間の元とvとの和の最小値を計算
    bs minimize(bs v) const {
        for (int i = 0; i < (int)basis.size(); ++i) {
            if (v.test(msb[i])) v ^= basis[i];
        }
        return v;
    }

    // 部分空間の元とvとの和の最大値を計算
    bs maximize(bs v) const {
        for (int i = 0; i < (int)basis.size(); ++i) {
            if (!v.test(msb[i])) v ^= basis[i];
        }
        return v;
    }

    // 部分空間の次元を返す
    int dim() const { return basis.size(); }
};

おわりに

掃き出し法の文脈で出てくる「標準形」に対して、このアルゴリズムでは、行のswapをしない標準形、いわば「準標準形」で基底が蓄えられていくのが面白いです。このアイデア自体はXORの掃き出しをする場合に限ったものではないと思いますが、XORの場合は簡単に実装できるという話でした。