思路:马尔可夫过程,转为求解方程组:每一种状态有一个概率,把概率转移方程列出来,用高斯消元法求解。
由于一共有2^n-2个方程,n最大为15,高斯消元法解方程是三次方的,无法通过极限数据。方程组是稀疏的,不知道要怎么优化。
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<algorithm>
#include<queue>
#include<math.h>
#include<vector>
#include<iostream>
using namespace std;
#define N 15
#define pw(x) (1<<x)
int n,m;
double a[1<<N][N+1]; //
double p[N+1][N+1];
vector<int> qry;
int read(){
scanf("%d%d",&n,&m);
for (int i=1;i<n;i++)
for (int j=i+1;j<=n;j++){
scanf("%lf",&p[i][j]);
p[j][i]=1.0-p[i][j];
}
for (int i=1;i<=m;i++){
int S=0;
for (int j=1;j<=n;j++){
int k;
scanf("%d",&k);
if(k)S+=pw(j-1);
}
qry.push_back(S);
}
return 0;
}
#define OPPSTATE(S) (((1<<n)-1)^S) //对手状态
#define CONTAIN(S,c) ((S&(1<<(c-1)))>0) //S状态是否包含卡牌c
#define FORSTATE(S,c) for (int c=1;c<=n;c++)if (CONTAIN(S,c)) //枚举S中的卡牌
double p1[N];
double p2[N];
double sum[N];
void calcP1P2(int S,double* p1){ //计算在状态S时,出每张牌的概率
int S2=OPPSTATE(S);
for (int i=1;i<=n;i++) p1[i]=sum[i]=0.0;
double ssum=0;
FORSTATE(S,c1){
FORSTATE(S2,c2){
sum[c1]+=p[c1][c2];
}
ssum+=sum[c1];
}
FORSTATE(S,c1) p1[c1]=sum[c1]/ssum;
}
void calcA(){ //计算两两状态转移概率
memset(a,0,sizeof(a));
for (int S=1;S<(1<<n)-1;S++){ //状态0和(1<<n)-1为最终状态,不需要转移
int S2=OPPSTATE(S);
calcP1P2(S,p1);
calcP1P2(S2,p2);
//S 和 S2博弈
FORSTATE(S,c1)
FORSTATE(S2,c2){ //S 出c1,概率p1[c1],S2出c2概率p2[c2]
//S赢,得到c2
double prob_c1_win=p[c1][c2];
a[S][c2]+=p1[c1]*p2[c2]*prob_c1_win;
//S输,失去c1
a[S][c1]+=p1[c1]*p2[c2]*(1-prob_c1_win);
}
for (int i=1;i<=n;i++){
//printf("%d->%d %.2lf\n",S,S^pw(i-1),a[S][i]);
}
}
}
double mat[5000][5000];
void gauss(){//2^n-2个方程
memset(mat,0,sizeof(mat));
int sz=(1<<n)-2;
for (int S=1;S<=sz;S++){
for (int i=1;i<=n;i++){
int S2=S^pw(i-1);
if (S2==sz+1) mat[S][S2]=-a[S][i];
else
mat[S][S2]=a[S][i];
}
mat[S][S]=-1.0;
}
//for (int i=1;i<=sz;i++,putchar('\n'))
//for (int j=1;j<=sz+1;j++)printf("%.2lf ",mat[i][j]);
//
for (int i=1;i<=sz;i++){
int r=i;
for (int j=i+1;j<=sz;j++){
if (fabs(mat[j][i]) > fabs(mat[r][i])) r=j;
}
if (r!=i) for (int j=1;j<=sz+1;j++) swap(mat[i][j],mat[r][j]);
for (int j=1;j<=sz;j++)
if (j!=i){
double f=mat[j][i]/mat[i][i];
for (int k=i;k<=sz+1;k++) mat[j][k]-=f*mat[i][k];
}
}
//printf("\nafter gausss\n");for (int i=1;i<=sz;i++,putchar('\n'))
//for (int j=1;j<=sz+1;j++)printf("%.2lf ",mat[i][j]);
for (int i=1;i<=sz;i++){
mat[i][sz+1]/=mat[i][i];
// printf("STATE[%d] %.2lf\n",i,mat[i][sz+1]);
}
mat[0][sz+1]=0.0;
mat[sz+1][sz+1]=1.0;
}
int main(){
read();
calcA();
gauss();
for (int i=0;i<qry.size();i++){
printf("%.5lf\n",mat[qry[i]][(1<<n)-1]);
}
return 0;
}