ZOJ3822 ACM-ICPC 2014 亞洲區域賽牡丹江賽區現場賽D題Domination 概率DP
Edward is the headmaster of Marjar University. He is enthusiastic about chess and often plays chess with his friends. What's more, he bought a large decorative chessboard with N rows and M columns.
Every day after work, Edward will place a chess piece on a random empty cell. A few days later, he found the chessboard was dominated
"That's interesting!" Edward said. He wants to know the expectation number of days to make an empty chessboard of N × M dominated. Please write a program to help him.
Input
There are multiple test cases. The first line of input contains an integer T indicating the number of test cases. For each test case:
There are only two integers N and M (1 <= N, M <= 50).
Output
For each test case, output the expectation number of days.
Any solution with a relative or absolute error of at most 10-8
Sample Input
2 1 3 2 2
Sample Output
3.000000000000 2.666666666667
這道題算是一個概率DP的裸題吧,很裸的演算法,題目意思是告訴你有一個n*m的區域,佔領一個格子可以控制這個行和列,求控制所有行列所需要佔領格子的期望,直接推匯出狀態轉移方程即可,對於算概率,我設dp[i][j][k]表示的意思為當我用k步控制了i行j列的概率,那麼我們可以知道,這時候下一步則有4不狀態轉移,分別為繼續取已控制的i行j列的格子,取未控制的i行,已控制的j列的格子,取已控制的i行和未控制的j列的格子,取未控制的i行j列的格子四種不同的情況,對於每種情況,狀態轉移方程如下:
1、繼續取已控制的i行j列的格子
dp[i][j][k+1]+=dp[i][j][k]*(i*j-k)/(n*m-k);
2、取未控制的i行,已控制的j列的格子dp[i+1][j][k+1]+=dp[i][j][k]*(n-i)*j/(n*m-k);
3、取已控制的i行和未控制的j列的格子dp[i][j+1][k+1]+=dp[i][j][k]*(m-j)*i/(n*m-k);
4、取未控制的i行j列的格子dp[i+1][j+1][k+1]+=dp[i][j][k]*(n-i)*(m-j)/(n*m-k);
從而最終求出控制n行m列需要步數的概率,由期望的定義式EX=n1p1+n2p2+...即可求出最終期望,具體程式如下:
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#include<map>
#include<vector>
#include<queue>
using namespace std;
double dp[55][55][55*55],a[55][55];
const double eps=1e-8;
int maxx(int a,int b)
{
if(a>b)return a;
return b;
}
int main()
{
// freopen("in.txt","r",stdin);
int t;
cin>>t;
while(t--)
{
int n,m;
scanf("%d%d",&n,&m);
memset(dp,0,sizeof(dp));
memset(a,0,sizeof(a));
dp[1][1][1]=1;
a[1][1]=1;
for(int i=1;i<=n;i++)
{
for(int j=1;j<=m;j++)
{
for(int k=1;k<=a[i][j];k++)
{
if(i==n&&j==m)
break;
dp[i][j][k+1]+=dp[i][j][k]*(i*j-k)/(n*m-k);
if(i*j>k)
a[i][j]=maxx(a[i][j],k+1);
}
for(int k=1;k<=a[i][j];k++)
{
dp[i+1][j][k+1]+=dp[i][j][k]*(n-i)*j/(n*m-k);
a[i+1][j]=maxx(a[i+1][j],k+1);
dp[i][j+1][k+1]+=dp[i][j][k]*(m-j)*i/(n*m-k);
a[i][j+1]=maxx(a[i][j+1],k+1);
dp[i+1][j+1][k+1]+=dp[i][j][k]*(n-i)*(m-j)/(n*m-k);
a[i+1][j+1]=maxx(a[i+1][j+1],k+1);
}
}
}
double sum=0;
for(int i=1;i<=a[n][m];i++)
sum+=dp[n][m][i]*i;
printf("%.12f\n",sum);
}
return 0;
}