yukicoder No.174 カードゲーム（Hard）

問題
No.174 カードゲーム（Hard） - yukicoder
A君とB君がカードでゲームを行う．１から１０００までの数字が書かれたカードがあり，その中から $N$ 枚ずつ配られる．試合は $N$ 回あり，試合ごとにプレイヤーがカードを出していき，大きい数字を出したほうがその数字の合計分だけ得点が得られる．また，同じカードは１度しか使えない． $N$ 試合終わった後，A君が持っている得点の期待値は幾つか．
A君，B君は出すカードに癖があり，持っている中で一番小さい数字を出す確率がそれぞれ $P_A,P_B$ である．残りのカードを出す確率は等しい．

解法
bitDPを使い，A君，B君がそれぞれ使ったカードの集合が $s$ の時の確率は幾つか求める．
$k=N$ の時，
　 $dp[s|1{<<}i]=dp[s]$
$i$ 番目の数字が一番小さい時，
　 $dp[s|1{<<}i]=dp[s]*P_A$
その他の時，
　 $dp[s|1{<<}i]=dp[s]*(1-P_A)/(N-k)$
この値をうまく使いA君が $k$ 回目に $i$ 番目の数字を出す確率を求める．この確率にカードの数字をかけて，勝った時に得られる数字の期待値を求める．これらの累積和を取り， $i$ 番目の数字以上の数字を出した時の確率，期待値を出しておく．
B君が使ったカードの集合を列挙しその後に，B君があるカードを使うときのA君が勝つために必要な最低のカードを二分探索で探し，
A君がそのカード以上を出す確率にB君の確率とB君のカードの数字をかけた，B君のカードから得られる得点の期待値
A君がそのカード以上を出す時の得られる得点の期待値にB君の確率をかけた，A君のカードから得られる得点の期待値
を足し合わせる．これらのすべての合計が求める期待値になる．
計算量は $O(2^NN\log N)$

ソースコード

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

#define RREP(i,s,e) for (i = s; i >= e; i--)
#define rrep(i,n) RREP(i,(int)(n)-1,0)
#define REP(i,s,e) for (i = s; i <= e; i++)
#define rep(i,n) REP(i,0,(int)(n)-1)
#define INF 100000000

typedef long long ll;

int n;
double pb;
double dp[1<<20], dp2[1<<20], p[20][20], e[20][20];
int a[20], b[20];

int main() {
    int i, j;
    double pa;
    cin >> n >> pa >> pb;
    rep (i,n) cin >> a[i];
    rep (i,n) cin >> b[i];
    dp[0] = dp2[0] = 1;
    rep (i,1<<n) {
        bool first = true;
        int cnt = __builtin_popcount(i);
        rep (j,n) {
            if (!(1<<j&i)) {
                if (cnt == n-1) {
                    dp[i|1<<j] += dp[i];
                    dp2[i|1<<j] += dp2[i];
                    p[cnt][j] += dp[i];
                }
                else if (first) {
                    dp[i|1<<j] += dp[i] * pa;
                    dp2[i|1<<j] += dp2[i] * pb;
                    first = false;
                    p[cnt][j] += dp[i] * pa;
                }
                else {
                    dp[i|1<<j] += dp[i] * (1-pa) / (n-cnt-1);
                    dp2[i|1<<j] += dp2[i] * (1-pb) / (n-cnt-1);
                    p[cnt][j] += dp[i] * (1-pa) / (n-cnt-1);
                }
            }
        }
    }
    sort(a,a+n);
    sort(b,b+n);
    rep (i,n) rep (j,n) e[i][j] = p[i][j] * a[j];
    rep (i,n) rrep (j,n-1) {
        p[i][j] += p[i][j+1];
        e[i][j] += e[i][j+1];
    }
    double ans = 0;
    rep (i,1<<n) {
        int cnt = __builtin_popcount(i);
        bool first = true;
        rep (j,n) {
            if (!(i&1<<j)) {
                int index = lower_bound(a,a+n,b[j]) - a;
                if (cnt == n-1)
                    ans += index < n ? dp2[i]*(e[cnt][index]+p[cnt][index]*b[j]) : 0;
                else if (first) {
                    ans += index < n ? dp2[i]*(pb*e[cnt][index]+p[cnt][index]*pb*b[j]) : 0;
                    first = false;
                }
                else
                    ans += index < n ? dp2[i]*((1-pb)/(n-cnt-1)*e[cnt][index]+p[cnt][index]*(1-pb)/(n-cnt-1)*b[j]) : 0;
            }
        }
    }
    cout << setprecision(10) << ans << endl;
    return 0;
}