けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

TCO 2018 Round 3B Med TestProctoring

楽しい!!!!!!!!!すごく勉強になったん!!!!!
大きく 2 つのやり方があって、高速ゼータ変換を用いた O(3n) から O(n2n) への高速化や、期待値に関する重要な考察をするなど色々やれるん。

問題概要

 n 人がテストを行う。 毎秒ごとに  i 番目の人がテストを終える確率は  p_i/q_i で与えられる。全員がテストを終えるまでに要する時間の期待値を求めよ。

制約

  •  1 \le n \le 20

解法 1: O(3n) な期待値 DP を O(n2n) に高速化する

まず、1 人の場合は有名問題で、毎ステップごとに確率  p で終わる事象が終わるまでのステップ数の期待値は  \frac{1}{p} になる。これは期待値の方程式を立てればわかる:

 E = (1-p)(E+1) + p

これを解いて、 E = \frac{1}{p} が導出される。
それでは  m 人の場合 ( S で表される集合の人たち) の場合の式を書いてみる。ここで簡単のために

  •  P_{S} = \prod_{i \in S} p_{i}
  •  Q_{S} = \prod_{i \in S} (1 - p_{i})

とおいておく。

 E_{S} = Q_{S}(E_{S} + 1) + \sum_{T \subsetneq S} ( Q_{T}P_{S-T}(E_{T} + 1))

という感じの式になる。この式で  T は、 S のうち 1 秒後も終わらずに残ったものたちを表す。この式を  E_{S} について解くと、

 E_{S} = \frac{1}{1 - Q_{S}} (Q_{S} + \sum_{T \subsetneq S} ( Q_{T}P_{S-T}(E_{T} + 1))

となる。これはちょうど bit DP の枠組みで解けるように見える。状態量が  2^{N} な bit DP で、状態遷移が submask 全体に行き渡るタイプの bit DP なので計算量は  O(3^{N}) になる。このままでは間に合わないので頑張って高速化する。

 E_{S} の計算をするのに毎回  \sum_{T \subsetneq S} の形の足し算をするのが時間かかるので、ここをなんとかしたい。もしこれの中身が  S には依存せずに  T のみに依存する式だったら、

  •  t_{S} = \sum_{T \subsetneq S} ()

みたいな補助変数を用意して、これを更新しながら計算すればよさそうである。DP を累積和を用いて高速化する話によく似ている。これができるようにするために、式変形をする。具体的には  P_{S-T} の部分をくくり出す。ここで

  •  R_{S} = \frac{Q_{S}}{P_{S}} = \prod_{i \in S} (\frac{1}{p_{S}} - 1)

としておく。そうすると

  •  E_{S} = \frac{P_{S}}{1 - Q_{S}} (R_{S} + \sum_{T \subsetneq S} ( R_{T}(E_{T} + 1))

という風に式変形できる。これでめでたく sum のところの中身が  T のみに依存する式になったので、アイディアとしては

  •  t_{S} = \sum_{T \subsetneq S} ( R_{T}(E_{T} + 1))

とおくことで、これを毎回計算しながら頑張ればよい。 t_{S} の更新は高速ゼータ変換の亜種 (naoya_t さんのブログ の (2) バージョン) の出番である!

ここでいう高速ゼータ変換は


集合に対する関数  f(S) があって

 g(S) = \sum_{T \subseteq S} f(T)

を計算したい


というものを考える ( T の向きが逆で  T \supseteq S の方を高速ゼータ変換と呼ぶことの方が多いかも)。

これは以下のようなシンプルなコードが知られているが

for (int i = 0; i < n; ++i)
  for (int bit = 0; bit < (1<<n); ++bit)
    if (bit & (1<<i))
      f[bit] += f[bit ^ (1<<i)];

これは配列再利用をしているため少し見えづらくなっている。結局以下のような DP をしている:


t[ i ][ bit ] := bit のうち、下から i 桁については bit の subbit 全体を動かし、それ以外の桁については bit と同じ値に固定したもの全体の総和

とすると

  • 初期化: t[ 0 ][ bit ] = f(bit)
  • 更新:

t[ i + 1 ][ bit ] =
t[ i ][ bit ] ( i が bit に含まれないとき)
t[ i ][ bit ] + t[ i ][ bit ^ (1<<i) ] ( i が bit に含まれるとき)

  • 答え: t[ n ][ bit ]

なお、こうして t を二次元配列として定義してあげると、

for (int i = 0; i < n; ++i)
  for (int bit = 0; bit < (1<<n); ++bit)
    if (bit & (1<<i))
      f[bit] += f[bit ^ (1<<i)];

の i と bit のループの順番を変更しても大丈夫になる。これを利用して、dp[bit] を求めながら、t[n][bit] も求めていく。もう 1 つの注意点として、本来は t の初期化は

t[ 0 ][ bit ] = f(bit)

という感じで行うのだが、f(bit) の値が予めわからないので、後で

for (int i = 0; i <= n; ++i) t[i][bit] += R*(dp[bit] + 1);

という感じで足している。

#include <iostream>
#include <vector>
using namespace std;

class TestProctoring {
public:
    double expectedTime(vector <int> a, vector <int> b) {
        int n = (int)a.size();
        vector<double> p(n);
        for (int i = 0; i < n; ++i) p[i] = (double)(a[i]) / b[i];
        vector<double> dp(1<<n, 0.0);
        vector<vector<double> > t(n+1, vector<double>(1<<n, 0.0));
        for (int bit = 0; bit < (1<<n); ++bit) {
            double P = 1.0, Q = 1.0, R = 1.0;
            for (int i = 0; i < n; ++i) {
                t[i+1][bit] = t[i][bit];
                if (bit & (1<<i)) {
                    P *= p[i], Q *= (1.0 - p[i]), R *= (1.0/p[i] - 1);
                    t[i+1][bit] += t[i][bit ^ (1<<i)];
                }
            }
            if (bit) dp[bit] = P / (1 - Q) * (R + t[n][bit]);
            for (int i = 0; i <= n; ++i) t[i][bit] += R*(dp[bit] + 1);
        }
        return dp[(1<<n)-1];
    }
};

解法 2: 試行回数の期待値 = (k 回以上かかる確率) の和

てんぷらさんに教えてもらった解法です。kmjp さんの解法とも一致しています。今回はちょうど  k 回で終わる確率を  P_k として、期待値の定義から

  •  \sum_{k = 1}^{\infty} k P_{k}

を計算したいです。しかし「ちょうど  k 回になる確率」というのはとても計算しづらいので、代わりに

  •  Q_{k} = k 回以上かかる確率

というのを代わりに考えます。このとき  P_{k} = Q_{k} - Q_{k+1} になるので、

 \sum_{k = 1}^{\infty} k P_{k}

 = \sum_{k = 1}^{\infty} k (Q_{k} - Q_{k +1})

 = (Q_1 - Q_2) + 2(Q_2 - Q_3) + 3(Q_3 - Q_4) + \dots

 = Q_1 + Q_2 + Q_3 + \dots

 = \sum_{k = 1}^{\infty} Q_k

という感じになります (収束性についてはきちんと議論しないとですが...)。標語的に言えば


試行回数の期待値 = ( k 回以上かかる確率) の  k についての無限和


ということになります。これは「ヒストグラムを書いて縦に集計するものを横に集計し直したもの」というようなものの見方ができます。結構汎用的なものの見方な気がします。AGC 023 C - Painting Machines でも出て来たとか。

さて、各  k について  Q_k を計算します。 k 回以上かかるというのはすなわち「 k-1 秒たってもまだ誰かしらが終わってない確率」とういことになります。これは余事象をとると「 k-1 秒で全員が終わる確率」です。 r_{i} = 1 - \frac{a_{i}}{b_{i}} とおくと、 i さんが  k 秒たっても終わらない確率は  r_{i}^{k-1} となることから、 i さんが  k 秒で終わる確率は  1 - r_{i}^{k-1} になります。したがって、

 Q_{k} = 1 - (1 - r_{1}^{k-1})(1 - r_{2}^{k-1})\dots(1 - r_{n}^{k-1})

になります。これを展開すると、次のように  2^{n} - 1 個の項が出て来ます。 S を { 1, 2, \dots, n} の空でない部分集合を動くとして、

 Q_{k} = \sum_{S \in 2^{n}, S \neq \emptyset} (-1)^{|S|+1} \prod_{i \in S}r_{i}^{k-1}

という感じになります。この式は包除原理で導いたものと考えることもできます (kmjp さんのブログでは包除原理)。

この項それぞれについて、 k = 1, 2, \dots, についての無限級数の総和をとったものを求めればよいです。ここで

  •  R = \prod_{i \in S}r_{i}

とおくと、無限級数の和は

 1 + R + R^{2} + R^{3} + \dots = \frac{1}{1 - R}

になります。あとはこれを各  S \in 2^{n}, S \neq \emptyset について足し引きすればよいです。

#include <iostream>
#include <vector>
using namespace std;

class TestProctoring {
public:
    double expectedTime(vector <int> a, vector <int> b) {
        double res = 0.0;
        int n = (int)a.size();
        for (int bit = 1; bit < (1<<n); ++bit) {
            double r = 1.0;
            int bitcount = 0;
            for (int i = 0; i < n; ++i) if (bit & (1<<i)) r *= 1.0 - (double)(a[i])/b[i], ++bitcount;
            double tmp = 1.0 / (1.0 - r);
            if (bitcount & 1) res += tmp;
            else res -= tmp;
        }
        return res;
    }
};