TieForMax [SRM747 Div1 Med]

T個の皿（tokenだけど）と、P個の置き場所がある。
T個の皿をP個の置き場に一様な確率で順番に置いていく。
このとき、最も多く置かれている置き場所が複数ある（1位タイになっている）確率を求めよ。

T, P≦50

前提知識

• 確率DP （自分の解法では微妙に違うかも）

解説

独立になるよう場合分けをして計算していく。
一位タイになっている置き場の個数tie, 一位になっている山の個数maを全探索する。
tie,maが異なる場合は独立になる。
　
この2つを固定するが、とりあえず最初の1～tie番目の置き場所を一位にしておく。
他の置き場所が一位であっても確率は全て対照で変わらないので、C(P, tie)をかければいい。
このとき、tie+1番目～の山はmaより少ない個数の山を作る必要がある。
この後半部分はDPで計算しよう。
　
まず、全ての置き場所に a1, a2, a3, ... , aP個置かれている確率を考える。
これは、どんな順番で置かれても確率は(1/P)^Pであり、置く順番についての組み合わせは同じものを含む順列なので、P!/(a1!a2!a3!...aP!)である。
つまり、P!/(a1!a2!a3!...aP!) * (1/P)^Pが確率になる。
(P! * (1/P)^P) * 1/(a1!a2!a3!...aP!)
のような感じで、1/(a1!a2!a3!...aP!)の部分だけ異なるので、ここをDPで計算出来ないか考える。
dp[i][use] := 一位じゃないi番目の置き場まで終わっていて、use個既に置いたときの1/(a1!a2!...ai)
これでi+1番目に例えば4個置いたとすれば、dp[i+1][use+4] += dp[i][use] * 1/4!と遷移していけばいい。
このDPで一位として使われていないt-tie*ma個をma以上を置かないように分配していけば答えが得られる。
　
アルゴリズムはこれで終わりなのだが、今回は誤差対策で少し工夫をしている
1. doubleの代わりにlong doubleを使っている
2. P!と(1/P)^Pの計算を最後に一気にやらず、適宜分配しながらやっている
2番目のせいで以下のコードがとても読みにくいのだが、50!を最後にやるのは誤差が怖くてやめた。
代わりにDPの遷移などで1/4!とかをするときに、既に使われている個数を見ながら、
1/4 * 1/3 * 1/2 * 1/1
において、ここで残り20個の皿が残っていたら、50!の50*49*...*21は処理済みなので、
20/4 * 19/3 * 18/2 * 17/1
と分配して、ついでに(1/P)^Pも分配して、
20/(4P) * 19/(3P) * 18/(2P) * 17/P
というふうにしている。
この辺の計算はpre[st][cn]で前計算している。
pre[st][cn] = P(st,cn)/(cn!*P^cn)
　
※ こんなふうにやったけど、もっといい考え方があるかも
DP遷移に使ったpre[st][cn]もC(st,cn)*(1/P)^cnで組み合わせ×確率で遷移してて、
普通の確率DPとして捉える方法ありそう（こんな面倒な式変形ではなくて）

typedef long double ld;
ld dp[51][51];
ld pre[51][51];
ld com[51][51];
struct TieForMax {
    double getProbability(int t, int p) {
        rep(st, 1, 51) rep(cn, 1, 51) {
            pre[st][cn] = 1;
            rep(i, 0, cn) pre[st][cn] *= 1.0 * (st - i) / (cn - i) / p;
        }

        rep(i, 0, 51) com[i][0] = com[i][i] = 1;
        rep(i, 1, 51) rep(j, 1, i) com[i][j] = com[i - 1][j - 1] + com[i - 1][j];

        ld ans = 0;
        rep(tie, 2, p + 1) rep(ma, 1, t + 1) {
            int other = (p - tie) * (ma - 1);
            int rest = t - tie * ma;
            if (other < rest or rest < 0) continue;
            
            int n = p - tie;
            rep(i, 0, n + 1) rep(j, 0, rest + 1) dp[i][j] = 0;
            dp[0][0] = 1;
            rep(tie2, 0, tie) dp[0][0] *= pre[t - ma * tie2][ma];

            rep(i, 0, n) rep(use, 0, rest + 1) {
                dp[i + 1][use] += dp[i][use];

                int cnt2 = rest - use;
                rep(nxt, 1, min(cnt2, ma - 1) + 1) {
                    dp[i + 1][use + nxt] += dp[i][use] * pre[cnt2][nxt];
                }
            }

            ans += com[p][tie] * dp[n][rest];
        }

        return ans;
    }
};