Mister雑記

競プロやります。

ABC113-D 「Number of Amidakuji」 発展解法

問題リンク

この問題については先日解説記事を書いたのだが、実はこれより遥かに速い解法が存在する。しかもDPもbit全探索も要らない(代わりにより高度なアルゴリズムと考察を要するが)。

本記事ではその解法の方針と実装例について述べようと思う。

概要

略。

解説

前の記事と同様に、 dp_{i, j} = i段からなり、1が jに移動するようなあみだくじの総数」とする。

またこれ以降、一次元配列 dp_i線型代数における「ベクトル」とみなすことがしばしばある。

あみだくじと行列

まず、あみだくじの形状というのは「行列」という線形代数の概念で表現できる。

これを文で説明していると無限に時間がかかるので、以下の図で納得していただきたい。

f:id:misteer:20181105211026j:plain

こんな具合に、「基本行列の、横線が入る場所だけ列が入れ替わったもの」があみだくじを表現する行列となる。

そして考えうる全部の1段からなるあみだくじについてこの行列を求めて、それぞれベクトル dp_iに掛けて足し合わせれば dp_{i + 1}が求められる。 さらに行列演算の線型性から、これらの行列を全て足し合わせてから dp_iに掛けても同じ結果が得られる。

f:id:misteer:20181105211039j:plain

さらにさらに、 iの値に依らず dp_iに掛ける行列は常に同じであることから、この行列を Mとおくと

 \displaystyle dp_{H} = M^H dp_0

となる。なお dp_0は定義から第一成分が1で他は全て0のベクトルである。

求める答えは dp_{H, K}なので、これは行列 M^H (K, 1)成分となる。

行列の累乗

さて、これで行列 M^Hを計算することで答えが求まることが分かった。普通に行列の累乗を計算していると、

  • 行列同士の積の計算量は O(W^3)
  • これを H回繰り返す

より、計算量は O(W^3 H)となる。しかし、これを O(W^3 \log H)に落とし込む方法が存在する。それが繰り返し二乗法である。

繰り返し二乗法を大まかに説明すると、

  •  Hが奇数なら M^{H - 1}Mにバラす
  •  Hが偶数なら \left({M^2}\right)^{H/2}に圧縮する

といった具合に再帰的に累乗を求めるアルゴリズムのことを指す。これで指数が高々2回毎に半分になるため、計算回数が O(\log H)に落とせる。1

あみだくじの数とFibonacci数列

後は Mを求めるだけとなった。横線の全パターンをbit全探索するという方針では、計算量は O(2^W)となる。

ところが実は、 Mはbit全探索なんてしなくても計算で求めることができる。鍵となるのは、「縦線 n本のあみだくじにおける横線のパターン数は fib(n)となる」という法則である。ここで fib(n)Fibonacci数列の第 n項である。

これは数学的帰納法で示せる。以下の図のように、 n本目と n-1本目の間に線を引くか否かで場合分けすればよい。

f:id:misteer:20181105211101j:plain

これで Mを計算する準備は整った。定義から M_{i, j} = i jに行き着くようなあみだくじの数」となる。段数が1なので |i - j| \gt 1なら M_{i, j} = 0となる。そうでない場合は、下の図のようにして求められる。

f:id:misteer:20181105211113j:plain

ただしindexに注意されたい。上の図は1-indexedの場合の値であり、0-indexedの場合は引数がズレることになる。

実装例

#include <iostream>
#include <vector>

using namespace std;

using ll = long long;

const ll MOD = 1e9 + 7;

ll fib[1000];

// Fibonacci数列を埋める
void precalc() {
    fib[0] = 1;
    fib[1] = 1;
    for(int i = 2; i < 1000; ++i) {
        fib[i] = (fib[i - 1] + fib[i - 2]) % MOD;
    }
    return;
}

// 正方行列の積
vector<vector<ll>> matmul(vector<vector<ll>> a, vector<vector<ll>> b) {
    int l = a[0].size();

    vector<vector<ll>> ret(l, vector<ll>(l, 0));

    for (int i = 0; i < l; ++i) {
        for (int j = 0; j < l; ++j) {
            for (int k = 0; k < l; ++k) {
                ret[i][j] += a[i][k] * b[k][j];
                ret[i][j] %= MOD;
            }
        }
    }
    return ret;
}

// ダブリングによる行列累乗
vector<vector<ll>> matpow(vector<vector<ll>> m, int n) {
    int l = m[0].size();

    // 単位行列を返す
    if (n == 0) {
        for (int i = 0; i < l; ++i) {
            for (int j = 0; j < l; ++j) {
                m[i][j] = (i == j);
            }
        }
        return m;
    }

    if (n == 1) return m;

    if (n % 2) {
        return matmul(matpow(m, n - 1), m);
    } else {
        return matpow(matmul(m, m), n / 2);
    }
}


int main() {
    // Fibonacci数列を先に求める
    precalc();

    int H, W, K;
    cin >> H >> W >> K;
    --K;

    vector<vector<ll>> M(W, vector<ll>(W, 0));
    // 以降Mは0-indexedであることに注意

    // 左に移動するパターン
    for (int i = 2; i <= W; ++i) {
        M[i - 2][i - 1] = fib[i - 2] * fib[W - i] % MOD;
    }

    // 右に移動するパターン
    for (int i = 1; i <= W - 1; ++i) {
        M[i][i - 1] = fib[i - 1] * fib[W - i - 1] % MOD;
    }

    // 真下に移動するパターン
    for (int i = 1; i <= W; ++i) {
        M[i - 1][i - 1] = fib[i - 1] * fib[W - i] % MOD;
    }

    auto res = matpow(M, H);
    cout << res[K][0] << endl;
    return 0;
}

なおこの実装による全体の計算量は O(W^3 \log H)となる。 W=100, H = 10^{18}のようなケースでもTLEにはならない。

感想

  • 制約が上くらい厳しければ700点はあると思う。
    • 700点の部分点として今回の制約が400点、でも良いかもしれない。
  • 横線のパターン数がFibonacci数であることに気づけなかった。

  1. 注意点として、これが使えるのは行列の積が結合法則を満たすためである。結合法則を満たさないような演算に対して繰り返し二乗法は使えない。