https://atcoder.jp/contests/abc199/tasks/abc199_f

前提知識

行列累乗

解説

https://atcoder.jp/contests/abc199/submissions/22041784

この問題は数学パートと行列累乗高速化パートに分かれる。
行列累乗を使ったDP高速化を知らない場合は解けないのでまずは学習が必要になる。
「行列累乗競技プログラミング」あたりで検索すると記事が出てくるかと思う。

何から始めるか

行列累乗に慣れていれば比較的方針が見えやすいが、そうでない場合は方針を出すことも難しいと思う。
今回の問題でまず難しい部分が、K回の操作後のK回が10⁹である部分である。
10⁹回ということは普通にループで回していくのは無理なので、数学を頑張ってO(1)にするか、なにかしらの方法でO(logK)にする必要性が出てくる。
この時点でだいぶ解法が絞られてくる。

次に注目されるのが、計算後にすべての頂点に対して答えを求める必要がある部分である。
それぞれ計算するという方針も当然考えるが、今回の問題はとある頂点に対する操作が他の部分に影響する。
なので、一気に計算が進められるのではないかと考えられる。

最後に制約を見ると、行列累乗である最後のピースを埋めるようにN≦100となっているので、ほぼ行列累乗かなと結論付く。

数学パート

正直、行列累乗で作るならこれしかないという感じで考えたので、ACによる証明をしている。

K回の操作後に頂点に書かれている数の期待値を求めよというのが問題であるが、1回の操作後をまずは考えてみよう。
とある頂点vに書かれる数の候補としては、

その頂点vとつながっていない辺が選択され、その頂点に書いてある数A[v]そのまま
その頂点vにつながっている辺が選択され、その辺で繋がっている数との平均値になる

となる。その頂点につながっていない辺の個数をcntとすると、場合1についてはcnt/Mの確率で起こることになる。
なので、cnt/Mの確率で数A[v]となる。
場面2については、つながっているすべての頂点wについて、(A[v]+A[w])/2となり、確率がそれぞれ1/Mとなる。
これを総合すると期待値計算は、

A[v]'
= A[v] * cnt/M + (A[v]+A[w1])/2 * 1/M + (A[v]+A[w2])/2 * 1/M + ...
= A[v] * (cnt + 1/2 + 1/2) / M + A[w1] * 1/(2M) + A[w2] * 1/(2M) + ...

となる。この計算をすべての頂点について行えば1回分の期待値を計算できる。
そしてK回分の期待値を計算するには同様の計算をK回行えばいい（←ACによる証明）
かつ、ここで重要なのが「この計算は常に同じ係数で計算される」ということである。

行列累乗高速化パート

前の立式を見るとA[v]'はA[v], A[w1], A[w2], ...に「一定の係数をかけて足し合わせたもの」になっている。
他の頂点についても同様のことがいえ、このような計算は行列とベクタを使うことで実現できることが知られている。
例えば、

A[2]' = x * A[0] + y * A[1] + z * A[2]

であるなら、

|     |   |     |   |A[0]|  
|     | = |     | * |A[1]|  
|A[2]'|   |x y z|   |A[2]|

のような行列とベクタによる計算に帰着できる。（他の部分は省略してある）
全ての頂点について同様の係数計算を行う事で1回の操作を行う行列Maを得ることができる。
問題で与えられる配列Aを初期状態としてMa^K * Aを行うことで最終的に得られるベクタが答えとなる。
そして、このMa^KはMaの行列のサイズがN×Nであるため、行列累乗を使用することでO(N^3logK)で求めることができる。
よって、これまでの計算をK回行う操作は行列累乗で高速化でき、ACすることができる。

行列累乗に渡す配列と遷移をDPのようにとらえて、行列累乗を使ったDP高速化と捉えることもできる。
色々な見方ができるようになると、色々な場面で引き出しやすくなるだろう。

int N, M, K, A[101];
int X[10101], Y[10101];
int cnt[101];
bool E[101][101];
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> N >> M >> K;
    rep(i, 0, N) cin >> A[i];
    rep(i, 0, M) cin >> X[i] >> Y[i], X[i]--, Y[i]--;

    Mat Ma(N, Vec(N));
    rep(i, 0, M) {
        cnt[X[i]]--;
        cnt[Y[i]]--;
        E[X[i]][Y[i]] = E[Y[i]][X[i]] = true;
    }

    rep(to, 0, N) {
        int cnt = M;
        rep(from, 0, N) if (to != from) {
            if (E[to][from]) {
                Ma[to][to] = add(Ma[to][to], mul(modinv(M), modinv(2)));
                Ma[to][from] = add(Ma[to][from], mul(modinv(M), modinv(2)));
                cnt--;
            }
        }
        Ma[to][to] = add(Ma[to][to], mul(cnt, modinv(M)));
    }

    Vec v(N, 0);
    rep(i, 0, N) v[i] = A[i];
    
    Ma = fastpow(Ma, K);
    v = mulMatVec(Ma, v);
    rep(i, 0, N) cout << v[i] << endl;
}

はまやんはまやんはまやん

hamayanhamayan's blog

Graph Smoothing [AtCoder Beginner Contest 199（Sponsored by Panasonic） F]

前提知識

解説

何から始めるか

数学パート

行列累乗高速化パート