公式

H - Stroll 解説 by Nyaan


この問題は、数列の値を高速に求めるための強力な武器となる 分割統治 FFT と通称されるアルゴリズムを取り上げた問題です。
FFT / 畳み込みを使ったアルゴリズムの中では難しい部類に属しますが、理解さえしていれば 20 行前後で実装できるアルゴリズムなので、数え上げに関心がある人は是非とも勉強してみてください。

ナイーブなDP解

初めに愚直な DP 解法を考えてみましょう。 DP テーブル \(d_{s,t}\) を次のように定めます。

  • \(d_{s, t} :=\) \(t\) キロメートル歩いて地点 \(s\) にいる場合の数

求める答えは \(d_{1, T}\) です。
このとき \(d_{s,t}\) は全ての \(u \lt t\) に対して \(d_{\ast, u}\) がわかっていたら計算できるので、 \(t\) の昇順に DP テーブルを埋めていくことができます。
より厳密に説明すると、 \(d_{s,t}\) は問題文の変数を利用して以下の式で表すことができます。

\[d_{s,t} = \sum_{\substack{b_i = s \\ 1 \leq i \leq M }} \left( \sum_{0 \leq u \lt t} d_{a_i, u} \times p_{i, t - u} \right) + \sum_{\substack{a_i = s \\ 1 \leq i \leq M }} \left( \sum_{0 \leq u \lt t} d_{b_i, u} \times p_{i, t - u} \right) \]

上の DP の計算量は \(\mathrm{O}(M T^2)\) になりますが(詳しくは後述します)、与えられた制約では \(MT^2 = 1.6 \times 10^{10}\) になり定数倍もよくないため TLE してしまいます。この \(T^2\) の部分をいかに計算量を落とすかが、これから説明するアルゴリズムの肝となっています。

なぜ \(\mathrm{O}(T^2)\) かかるのか?

アルゴリズムを説明する前に、この問題を \(T\) に対して準線形で解くことがなぜ本質的に難しいのかを説明していきたいと思います。

以下の章では説明を簡単にするため \(N = 2\) の場合を説明します。\(N \geq 3\) の場合も辺の本数の \(\mathrm{O}(M)\) が計算量に掛かることを除けば本質的に同じです。

\(N = 2\) のときの DP テーブル \(d_{1,t}, d_{2,t}\) は、\(1 - 2\) 間の長さ \(u\) の辺の本数を \(p_u\) とおくと次の式で与えられます。

\[ d_{1,t} = \sum_{0 \leq u \lt t} d_{2,u} \times p_{t-u} \]

\[ d_{2,t} = \sum_{0 \leq u \lt t} d_{1,u} \times p_{t-u} \]

最近の ABC を解いている方の中には、もしかしたら「この式は畳み込みの式に似ている!」と気づいた方がいるかもしれません。もしも上の式の右辺の \(d_{1,u}, d_{2,u}\) が定数だった場合、 FFT を利用した畳み込みを用いて \(\mathrm{O}(T \log T)\) で高速に計算することができます。
一方、この問題で FFT をナイーブに適用しようとすると、 \(d_1\) を計算するには \(d_2\) の値があらかじめ必要で、 \(d_2\) を計算するには \(d_1\) があらかじめ必要で…という堂々巡りに陥ってしまいます。

この問題のように DP テーブルを埋めるのに DP テーブルが必要になる DP は オンライン DP と呼ばれています。
上式を見るとわかる通り、オンライン DP は他の DP テーブルの値に依存するため高速化が難しく、一般には DP テーブルを DAG 順に愚直に埋めていく方法でしか解くことができません。そして、上の式を愚直に計算すると DP テーブルを埋めるのに \(\mathrm{O}(T^2)\) かかってしまいます。

さて、ここまではこの問題の持つ厄介な性質を説明してきました。困ってしまったかのようですが、実はオンライン DP のうちいくつかの良い性質を持つケースに関しては、愚直な解法から計算量を改善できることが知られています。
そしてこの問題では寄与が畳み込みの形で表せる性質を利用すると、分割統治法と畳み込みを組み合わせることで「オンラインでも計算可能な畳み込み」を実現することができるのです!

分割統治 FFT

さて、本題に入りましょう。この章では 分割統治 FFT と呼ばれるアルゴリズムを説明していきます。( なお、名称は他にも cdq 分割統治 + FFT / オンライン FFT / オンライン畳み込み / relaxed multiplication など資料によって多様ですが、この文章では「分割統治 FFT 」で統一します。)

まずはじめに、区間 \([l, r)\) における DP テーブルの値を計算する DP の遷移を図で表したものを以下に示します。ここで一番左の丸は \(d_{\ast, l}\) の値に、一番右の丸は \(d_{\ast, r-1}\) の値に対応しています。

この形の DP を愚直に計算すると頂点数の 2 乗の計算量がかかってしまうことは上のグラフが無向完全グラフに向きをつけたものであることからもわかると思います。

分割統治 FFT では、上記の遷移の形をした DP を以下の 3 ステップを経て段階的に DP を計算するアルゴリズムになっています。

以下の説明では \(m = \left\lfloor \frac{l+r}{2} \right\rfloor\) とします。また不変条件として次の条件が常に成り立つことを保証します。

  • 区間 \([l, r)\) の値を計算する時点で次の2つの条件が満たされている。
    • \([0, l)\) の DP テーブルの値が確定している。
    • \([l ,r)\) の DP テーブル上には \([0, l)\) から \([l, r)\) への寄与の和が載っている。

1. 区間 \([l, m)\) の値を再帰的に計算する

まずはじめに、区間 \([l, m)\) の値を再帰的に計算して、区間 \([l, m)\) の DP テーブルの値を確定させます。不変条件より \([l, r)\) 外部からの寄与が DP テーブルに載っているので、 \([l, m)\) の値を問題なく計算することができます。
実装としては f(l, r) の内部で f(l, m) を呼び出す形になります。

2. \([l, m)\) から \([m, r)\) への寄与を計算する

次に \([l, m)\) から \([m, r)\) への寄与を計算します。
\(N = 2\) のとき、説明に対応する寄与を式で表すと次のようになります。( \(\leftarrow\) は代入の意味で用いています。)

\[ d_{1,t} \leftarrow d_{1,t} + \sum_{l \leq u \lt m} d_{2,u} \times p_{t-u} \]

\[ d_{2,t} \leftarrow d_{2,t} + \sum_{l \leq u \lt m} d_{1,u} \times p_{t-u} \]

ここで、先ほどは \(d_{1,u}, d_{2,u}\) の値が定まっていないせいで上の遷移を FFT で高速化することができない、という話をしました。
しかし今回はあらかじめ分割統治で \([l, m)\) の値が計算済みになっているので、上の式の \(\Sigma\) の部分を畳み込みでまとめて計算する高速化が適用できるようになっているのです!
よってこの部分は FFT を利用すると\(r - l = B\) として \(\mathrm{O} (B \log B)\) で計算することができます。

3. 区間 \([m, r)\) の値を再帰的に計算する

ここまで来ればあとはウイニングランです。最後に区間 \([m, r)\) の値を再帰的に計算して、区間 \([m, r)\) の DP テーブルの値を確定させます。

計算量は?

このアルゴリズムの計算量解析を考えてみましょう。幅 \(B\) の区間に分割統治 FFT を適用するときの計算量を \(T(B)\) と表すと、

\[T(B) = 2 T\left(\frac{B}{2}\right) + \mathrm{O} (B \log B)\]

という式で表せます。この式のように再帰が絡む形の計算量評価は厳密にやろうとすると難しいので、ここでは簡便な方法を紹介します。

\[b = \log B, U(b) = \frac{T(2^b)}{2^b} = \frac{T(B)}{B}\]

とおきます。元の式の両辺を \(B\) で割ると

\[ U(b)= U(b - 1) + \mathrm{O}(b) \]

を得るのでここから \(U(b) = \mathrm{O}(b^2) = \mathrm{O}(\log ^2 B)\) を示すことができます。よって

\[ T(B) = B \times U(b) = \mathrm{O}(B \log^2 B)\]

を得ることができます。

実装例

以上に説明したアルゴリズムを区間 \([0, T+1)\) に適用することで、この問題を \(\mathrm{O}(M T \log^2 T)\) で解くことができました。
実装例を以下に載せます。( 実装例のように AtCoder Library を使用すると実装がかなり楽になります。)

#include <iostream>
#include <utility>
#include <vector>

#include "atcoder/convolution.hpp"
#include "atcoder/modint.hpp"
using namespace std;
using mint = atcoder::modint998244353;

vector<mint> d[11], p[11][11];
vector<pair<int, int>> es;

vector<mint> get_vec(const vector<mint>& v, int l, int r) {
  return {begin(v) + l, begin(v) + r};
}
void online_convolution(int l, int r) {
  if (l + 1 == r) return;
  int m = (l + r) / 2;
  online_convolution(l, m);
  for (auto& [a, b] : es) {
    auto da = get_vec(d[a], l, m);
    auto pab = get_vec(p[a][b], 0, r - l);
    auto ra = atcoder::convolution(da, pab);
    for (int i = m; i < r; i++) d[b][i] += ra[i - l];
    auto db = get_vec(d[b], l, m);
    auto rb = atcoder::convolution(db, pab);
    for (int i = m; i < r; i++) d[a][i] += rb[i - l];
  }
  online_convolution(m, r);
}

int main() {
  ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int N, M, T, a, b, x;
  cin >> N >> M >> T;
  for (int i = 0; i < M; i++) {
    cin >> a >> b;
    --a, --b;
    es.emplace_back(a, b);
    p[a][b].resize(T + 1);
    for (int j = 1; j <= T; j++) cin >> x, p[a][b][j] = x;
  }
  for (int i = 0; i < N; i++) d[i].resize(T + 1);
  d[0][0] = 1;
  online_convolution(0, T + 1);
  cout << d[0][T].val() << "\n";
}

類題

AtCoder で過去に出題された関連問題を以下に何問か挙げるので練習にお使い下さい。

  • 第二回全国統一プログラミング王決定戦本戦 F - Jumping over 2 (日本語問題文のみ)

    • コンテスト中の AC は 5 人と少ないですが、分割統治 FFT を理解していればスムーズに考察できる問題だと思います。
  • ABC209 F - Deforestation (非想定解)

    • この問題は包除原理を利用した考察を行うことで分割統治 FFT が適用できる形になり \(\mathrm{O}(N \log^2 N)\) で解くことができます。(考察のパートは少し難しいと思います。)
      実際に実装する場合は、法が \(1000000007\) と FFT に適していないので、畳み込みを複数の mod で行い Garner のアルゴリズムを利用する必要があることに注意してください。
  • UTPC2020 D - Idol Group Costume (日本語問題文のみ)

    • この問題は解説で取り上げた分割統治 FFT とは少し違う形の遷移をしていますが、分割統治と FFT を組み合わせると解ける問題です。この形の DP もよく出題されるので抑えておくとよいでしょう。

投稿日時:
最終更新: