naoya_t@hatenablog

いわゆるチラシノウラであります

TCO 18 Round3B Medium (500): TestProctoring

去年の TCO 18 Round3Bの2問目(Medium:500点問題)。
https://community.topcoder.com/stat?c=problem_statement&pm=14766&rd=17205

問題はこんなやつ

  • あなたは試験監督。
  • N人の生徒(最大20人)が受験中である。
  • 毎秒(毎分だったか)、生徒ip_i の確率で回答を終えて教室を出ていく。(入力としてはint/intの分数で与えられる)
  • 全員出ていけばあなたの仕事は終わり。
  • あなたの待ち時間の期待値は?

そんなに難しそうには見えない?サンプルケースの最初の2件はこんな感じ*1

  • Example 0: { 2/4 } → 2.0
  • Example 1: { 1/2, 2/4 } → 2.6666666666666665

1人1人の生徒の退室までの時間の期待値を独立に求めて済む問題ではない、という話。


例えば3人いてそれぞれの(出ていく)確率をp,q,rとしたら、
E_{111} = (1-p)(1-q)(1-r)\cdot(1+E_{111}) + (1-p)(1-q)r\cdot(1+E_{110}) + (1-p)q(1-r)\cdot(1+E_{101}) + (p(1-q)(1-r)\cdot(1+E_{011})+(1-p)qr\cdot(1+E_{100}) + p(1-q)r\cdot(1+E_{010}) + pq(1-r)\cdot(1+E_{001}) + pqr\cdot(1+E_{000})
のような式になるのだけれど、本番でO(3^N)解を落とされてその後ちゃんと通してなかったので、どうすればO(N2^N)の形に持っていけるか考えた。

cf.
naoyat.hatenablog.jp

スタバで式をこねこねしても綺麗な(高速ゼータ変換に持っていけるような)形が見えてこなくて、コーヒーとシュガードーナツ1つ(と時間)を溶かした。E_{s} ではなく F_{s}=(something)\times(1+E_{s}) のような項で累積していくのではというところまでは考えついたのだけれどあと一息だった。(1-p)(1-q)r のような何かを掛けたり割ったりするのだろうと思ったけれどうまく導出できず。

けんちょん先生のヒント(というか模範解答)を読む。

P_s,Q_s,R_s を集合sのビットの立っている項のそれぞれ p_i, 1-p_i, \frac{1-p_i}{p_i} の積とすると \displaystyle E_s=\frac{P_s}{1-Q_s}(R_s+\sum_{t⊆s}R_t(1+E_t) のように Σ の中身が部分集合 t のみに依存した形になり、ここを \displaystyle F_s=\sum_{t⊆s}R_t(1+E_t) とすることで高速ゼータ変換が使えるようになる。

なるほど…[tex:S,S-T,T}と集合を分けて綺麗にくくり出している。(自分みたいに当てずっぽうじゃなく、けんちょん先生みたいに筋道立てて考えられると嬉しいのん)

とはいえこの場合

(2) rep(i,N) rep(j,1<<N) if (j&(1<<i)) f[j]+=f[j^(1<<i)];

の形では都合が悪い。同じ配列を再利用しているのでメモリ効率が良いのだけれど、すべての部分集合を足し込んだ後でないと f[j] が定まらないので出来れば j, i の順でループして、jが来た時点で部分集合の和が求まったところで補正しながら進みたい。

これは f[i][j] のような2次元配列で普通にDPと思って書けば出来る。i で下から何ビット目まで足し込んだかを示し、j は集合を表している。iを回した後に f[N][j] に部分集合の総和が求まっているようにすればよい。

あとは補正。F_{s} を求める過程では E_s=0 を仮定しているので、E_s が確定し次第 R_s(1+E_s)F_s に加算していく。(これがやりたいためのループ順序の変更である)。

計算量としてはO(N2^N)に収まった。2次元配列で O(N2^N) のメモリを食うがここが問題で、TopCoderの256MB制限だと long double でここでMLEになる*2ため double しか使えない。(問題で求められる精度は double で十分なので大丈夫)

→AC

#include <bits/stdc++.h>
using namespace std;

#define rep(var,n)  for(int var=0;var<(n);++var)

class TestProctoring {
public:
    double expectedTime(vector<int> numer, vector<int> denom) {
        vector<double> p;
        rep(i, numer.size()) {
            if (numer[i] == denom[i]) continue;
            double pi = (double)numer[i] / denom[i];  // 0.0 < pi < 1.0
            p.push_back(pi);
        }
        if (p.empty()) return 1.0;

        int N = p.size(), P = 1 << N;

        vector<double> e(P, 0);
        vector<vector<double>> f(N+1, vector<double>(P, 0));

        rep(pat,P) {
            double P = 1.0, Q = 1.0, R = 1.0;
            rep(i,N) {
                if (pat & (1<<i)) {
                    f[i+1][pat] = f[i][pat] + f[i][pat ^ (1<<i)];
                    P *= p[i]; Q *= (1.0 - p[i]); R *= (1.0 - p[i])/p[i];
                } else {
                    f[i+1][pat] = f[i][pat];
                }
            }
            if (pat > 0) {
                e[pat] = P/(1.0 - Q)*(R + f[N][pat]);
            }
            rep(i,N+1) f[i][pat] += R * (1.0 + e[pat]);
        }

        return e[P-1];
    }
};

だがしかし

この問題、こういう形じゃなくても解けるんだよね
kmjp.hatenablog.jp

f(i) := i秒目の段階で未合格の人がいる確率
とすると、求める解はf(1)+f(2)+...となる。

ここで、f(i)は包除原理の要領で、2^N通りの項の加減算で求めることができる。
この時、f(1), f(2), ...における各項の値は等比数列を成すので、2^N通りの項それぞれでその和を求めればよい。

TopCoderOpen 2018 Round3B Medium TestProctoring - kmjp's blog

あ、はい…

f(i) = (i秒目に1人でも終わっていないやつがいる確率) = 1 - (i秒目に全員終わっている確率) で、
また (i秒目に全員終わっている確率) = Π(ある人がi秒目に終わっている) である。

ここで
(ある人がi秒目に終わっている) = 1 - (ある人がi秒居残っている) = 1 - (ある人の居残り率)^i
但し (居残り率) = 1 - (離脱率)

f(i) = 1 - Π(1 - (ある人の居残り率)^i)

これを展開して各項が見えると、求めたい f(1)+f(2)+... で居残り率の等比級数の和が出てきて、これは簡単な分数で表されるのであとは包除原理の要領で… というのは大体(けんちょん先生の記事のお陰で)理解できた。

あとでこちらの解法でも通しておこう。

*1:実際には分子と分母は別々のvectorで来る。

*2:sizeof(long double)=16である。10バイトしか使ってないにも関わらず。