[ACM]CCF CSP [201612-5]E题 卡牌游戏【75分的程序】

思路：马尔可夫过程，转为求解方程组：每一种状态有一个概率，把概率转移方程列出来，用高斯消元法求解。

由于一共有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;}