[HDOJ 4936] Rainbow Island [动态规划+高斯消元]

一共有20个岛，你现在在1号岛，每次到一个岛，都有pi的概率在两个岛之间连上一条无向边(可能连重边)，接着你会等概率的传送到集合Si中的一个岛上。

问把所有的岛连成一个联通块需要传送的次数的期望。

dp[s][i]为在s的状态下，现在在第i个岛，还需要的期望。其中s描述了现有每个联通块的规模。显然，dp[20][i]=0，然后我们要求的是dp[1,1,...,1][i]。

注意到dp[s][i]只会转移到联通块个数更少的状态，或者和自己联通块规模相同的状态。于是可以分层做高斯消元，把联通块个数更少的状态作为常数带入下一层的方程。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
#include <map>

using namespace std;

const double eps=1e-20;
int n;

struct Matrix {
double a[20][21];
double ans[20];
Matrix() {
for (int i=0;i<n;i++)
for (int j=0;j<=n;j++)
a[i][j]=0;
for (int i=0;i<n;i++) {
a[i][i]=-1;
a[i]
=1;
ans[i]=0;
}
}
void solve() {
int i,j,k;
for (j=0;j<n;j++) {
int maxi=j;
for (i=j+1;i<n;i++) if (fabs(a[i][j])>fabs(a[maxi][j])) maxi=i;
for (k=j;k<=n;k++) swap(a[j][k],a[maxi][k]);
for (i=j+1;i<n;i++) {
if (fabs(a[i][j])>eps) {
double tmp=a[i][j]/a[j][j];
for (k=j+1;k<=n;k++)
a[i][k]-=tmp*a[j][k];
}
}
}
for (i=n-1;i>=0;i--) {
double tmp=-a[i]
;
for (j=i+1;j<n;j++) tmp-=ans[j]*a[i][j];
ans[i]=tmp/a[i][i];
}
}
};

map<vector<int>,Matrix> dp;
vector<int> blueTo[20];
double rainbowP[20];

void cal(const vector<int> &s) {
int m=s.size(),i,j,k,l;
Matrix &dps=dp[s];
for (i=0;i<m;i++) {
double ranP=(double)(s[i]*(s[i]-1))/(n*(n-1));
for (k=0;k<n;k++) {
double p=ranP*rainbowP[k]*(1/(double)blueTo[k].size());
for (l=0;l<blueTo[k].size();l++) {
dps.a[k][blueTo[k][l]]+=p;
}
}
for (j=i+1;j<m;j++) {
vector<int> tmp(s);
ranP=(double)(s[i]*s[j])/(n*(n-1)/2);
tmp[i]+=tmp[j];
for (k=j+1;k<tmp.size();k++) tmp[k-1]=tmp[k];
tmp.resize(m-1);
sort(tmp.begin(),tmp.end());
if (dp.find(tmp)==dp.end()) cal(tmp);
Matrix &dpt=dp[tmp];
for (k=0;k<n;k++) {
double p=ranP*rainbowP[k]*(1/(double)blueTo[k].size());
for (l=0;l<blueTo[k].size();l++) {
dps.a[k]
+=p*dpt.ans[blueTo[k][l]];
}
}
}
}
for (k=0;k<n;k++) {
double p=(1-rainbowP[k])*(1/(double)blueTo[k].size());
for (l=0;l<blueTo[k].size();l++) {
dps.a[k][blueTo[k][l]]+=p;
}
}
//printf("Calculate:");
//for (i=0;i<m;i++) printf(" %d",s[i]);
//printf("\n");
//printf(" Matrix:\n");
//for (k=0;k<n;k++) {
//for (l=0;l<=n;l++)
//printf(" %lf",dps.a[k][l]);
//printf("\n");
//}
dps.solve();
//printf(" Answer:");
//for (i=0;i<n;i++) printf(" %lf",dps.ans[i]);
//printf("\n");
//printf("\n");
}

int main() {
int t,tt,i,m,j;
scanf("%d",&t);
for (tt=1;tt<=t;tt++) {
scanf("%d",&n);
dp.clear();
for (i=0;i<n;i++) scanf("%lf",&rainbowP[i]);
for (i=0;i<n;i++) {
blueTo[i].clear();
scanf("%d",&m);
for (j=0;j<m;j++) {
int x;
scanf("%d",&x);
blueTo[i].push_back(x-1);
}
}
dp[vector<int>(1,n)]=Matrix();
if (n!=1) cal(vector<int>(n,1));
printf("Case #%d: %.6lf\n",tt,dp[vector<int>(n,1)].ans[0]);
}
return 0;
}