高斯消元与概率相关

今天在集训时听大爷们讲题听到了sdoi2017r1硬币游戏的部分分做法，发现自己一点不会，赶紧学了一波；

首先一道板子题：bzoj3270 博物馆

这里其主要思想是根据边的转移方向列出一个方程，详细的方程解释这位大佬说的很清楚：题解

这里我主要再说一下烦了我很久的细节问题；

首先，直接通过转移得到的矩阵并不是最终要求的；

看一下转移方程的形式，以当前状态i为例；
一开始直接写成的是xi=a1x1+a2x2+a3x3….anxn;
但是常用的高斯消元等号的另一边必须是常数而不能是变量；
所以就要把xi挪到等号另一边；
这也就是为什么要把所有矩阵中满足(i,i)的位置都要减1的原因；

现在得到的矩阵仍然不能出解，因为常数项那一列全是0；
考虑在整个方程中哪里有常数？就是初始点；
初始点的概率一开始就为1；
设初始点为i(这里的初始点相当于一个二元组，就是题目中两人的起点）；
那么xi=a1x1+a2x2+a3x3….anxn+1;
然后把1和xi各移一个项，常数项那一列便不都是0了；
也就是把i对应的常数项那里设为-1；

这样高斯消元后就能出解了；

然后就是如何与ac自动机结合在一起；
其实也很简单，把ac自动机的每个点都假设有字符集那么多条出边，然后看每个点通过这些出边以及fail指针会转移到哪个点，像上边那个题一样做就可以；

代码：

bzoj3270
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define db double
bool mapp[21][21];
int du[21]; 
db p[21];
db jz[405][405];
db ans[401];
int n,m;
int a,b;
db eps=1e-10;
int lc(int a,int b){
    return (a-1)*n+b;
}
void init()
{
    for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)if(i!=j){
        jz[lc(i,j)][lc(i,j)]+=p[i]*p[j];
        for(int ki=1;ki<=n;ki++)if(mapp[i][ki]){
            jz[lc(ki,j)][lc(i,j)]+=(1.0-p[i])*p[j]/(db)du[i];
        }
        for(int kj=1;kj<=n;kj++)if(mapp[j][kj]){
            jz[lc(i,kj)][lc(i,j)]+=p[i]*(1.0-p[j])/(db)du[j];
        }
        for(int ki=1;ki<=n;ki++)for(int kj=1;kj<=n;kj++)if(mapp[j][kj]&&mapp[i][ki])
            jz[lc(ki,kj)][lc(i,j)]+=(1.0-p[i])*(1.0-p[j])/(db)(du[i]*du[j]);
    }
    for(int i=1;i<=n*n;i++)jz[i][i]-=1.0;
} 
void gauss(){
    for(int i=1;i<n*n;i++){
        int now=i;
        for(;(fabs(jz[now][i])<eps&&now<=n*n);now++);
        if(now>n*n)continue;
        else if(now!=i)
            std::swap(jz[i],jz[now]); 
        for(int j=i+1;j<=n*n;j++){
            if(fabs(jz[j][i])<eps)continue;
            db x=jz[j][i]/jz[i][i];
            for(int k=i;k<=n*n+1;k++)jz[j][k]=jz[j][k]-jz[i][k]*x;
        }
    }
    for(int i=n*n;i;i--){
        ans[i]=jz[i][n*n+1]/jz[i][i];
        for(int j=i-1;j;j--)jz[j][n*n+1]-=jz[j][i]*ans[i]; 
    }
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&a,&b);
    for(int i=1;i<=m;i++){
        int tp1,tp2;
        scanf("%d%d",&tp1,&tp2);
        mapp[tp1][tp2]=mapp[tp2][tp1]=true;
        du[tp1]++;
        du[tp2]++;
    }
    for(int i=1;i<=n;i++)scanf("%lf",&p[i]);
    init();
    jz[lc(a,b)][n*n+1]=-1.0;
    gauss();
    for(int i=1;i<=n;i++)printf("%.6lf ",fabs(ans[lc(i,i)])<eps ? 0.0:ans[lc(i,i)]);
    return 0;
}


bzoj1444 有趣的游戏

bzoj1444
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
const int maxn=103;
struct NODE{
    int next[26];
    int end;
    int fail;
}t[maxn];
db eps=1e-10;
int ntot=0;
void insert(char *s,int id){
    int len=strlen(s);
    int now=0;
    for(int i=0;i<len;i++){
        if(!t[now].next[s[i]-'A'])t[now].next[s[i]-'A']=++ntot;
        now=t[now].next[s[i]-'A'];
    }
    t[now].end=id;
}
int que[maxn],tail,head;
int n,l,m;
void bfs(){
    tail=head=0;
    que[tail++]=0;
    t[0].fail=-1;
    for(int now=que[head];head<tail;now=que[++head]){
        for(int i=0;i<m;i++)
        if(t[now].next[i]){
            que[tail++]=t[now].next[i];
            int k;
            for(k=t[now].fail;k!=-1;k=t[k].fail)if(t[k].next[i])break;
            t[t[now].next[i]].fail=k>=0 ? t[k].next[i] : 0;
        }
    }
}
char s[maxn];
db p[maxn];
db jz[maxn][maxn],jans[maxn],ans[maxn];
void getjz(){
    for(int i=0;i<=ntot;i++){
        if(t[i].end)continue;
        for(int k=0;k<m;k++){
            int now;
            for(now=i;now;now=t[now].fail)if(t[now].next[k])break;
            jz[t[now].next[k]][i]+=p[k];
        }
    }
    for(int i=0;i<=ntot;i++)jz[i][i]-=1.0;
    jz[0][ntot+1]=-1.0;
}
void gauss(){
    for(int i=0;i<=ntot;i++){
        int now=i;
        for(;fabs(jz[now][i])<eps&&now<=ntot;now++);
        if(now>ntot)continue;
        else if(now!=i)std::swap(jz[i],jz[now]);
        for(int j=i+1;j<=ntot;j++){
            db x=jz[j][i]/jz[i][i];
            for(int k=0;k<=ntot+1;k++)jz[j][k]-=x*jz[i][k]; 
        }
    }
    for(int i=ntot;i>=0;i--){
        jans[i]=jz[i][ntot+1]/jz[i][i];
        for(int j=i-1;j>=0;j--)
            jz[j][ntot+1]-=jz[j][i]*jans[i];
    }
}
int main()
{
    scanf("%d%d%d",&n,&l,&m);
    for(int i=0;i<m;i++){
        int tp1,tp2;
        scanf("%d%d",&tp1,&tp2);
        p[i]=(db)tp1/tp2;
    }
    for(int i=1;i<=n;i++){
        scanf("%s",s);
        insert(s,i);
    }
    bfs();
    getjz();
    gauss();
    for(int i=1;i<=ntot;i++)if(t[i].end)ans[t[i].end]=jans[i];
    for(int i=1;i<=n;i++)printf("%.2lf\n",fabs(ans[i])<eps? 0.00 :ans[i]);
    return 0;
}