• 微信公众号:美女很有趣。 工作之余,放松一下,关注即送10G+美女照片!

[gym102268E]Expected Value

开发技术 开发技术 2周前 (04-30) 5次浏览

令$X$为移动次数,答案即$sum_{i=0}^{infty}P(X>i)$,后者记作$S_{i}$

关于$S_{i}$,令$f_{i,j}$表示走了$i$步后位于$j$且未到达过$k$的概率,即有$S_{i}=sum_{jin V,jne t}f_{i,j}$

初始状态即$f_{0,s}=1$,转移即$f_{i,j}=sum_{(k,j)in E}w_{(k,j)}f_{i-1,k}$(特别的,强制$w_{(t,i)}=0$)

令矩阵$M$为临接矩阵,即$M_{i,j}=w_{(i,j)}$,显然有$f_{i}=f_{i-1}M$

根据矩阵乘法的结合律,有$f_{i,j}=M^{i}_{s,j}$,那么$S_{i}$即$M^{i}$中第$s$行除第$t$列以外所有元素之和

定义${a_{0},a_{1},…,a_{n}}$的最短递推数列:满足$mge 1$,$r_{0}ne 0$且$forall mle jle n,sum_{i=0}^{m}r_{i}a_{j-i}=0$的数列中,$m$最小的数列${r_{0},r_{1},…,r_{m}}$

考虑求出${S_{i}}$的最短递推数列${r_{0},r_{1},…,r_{m}}$,即使得$forall jge m,sum_{i=0}^{m}r_{i}S_{j-i}=0$

对于这个$sum_{i=0}^{m}r_{i}S_{j+i}=0$,也即$M^{j}sum_{i=0}^{m}r_{i}M^{i}$中第$s$行除第$t$列以外的元素求和

因此,$sum_{i=0}^{m}r_{i}M^{i}=0$(这个$0$指矩阵中每一个元素都为0)是$forall jge m,sum_{i=0}^{m}r_{i}S_{j-i}=0$的充分条件

称多项式$f(x)=sum_{i=0}^{n}a_{i}x^{i}$为矩阵$M$的零化多项式,当且仅当$f(M)=0$

结论1(Cayley-Hamilton定理):令$I$为单位矩阵,$det(xI-M)$即一个关于$x$的$|V|$次多项式,此即为矩阵$M$的零化多项式(之一)

但如果直接根据此定理求出矩阵$M$的零化多项式作为${r_{i}}$,通过拉格朗日插值法选择若干个$x$代入来求,每一次$o(|V|^{3})$的高斯消元,时间复杂度为$o(|V|^{4})$

(另外,我们要求的是最短递推数列,其也不满足此性质,虽然这并不重要)

虽然不能直接计算,但从这个结论中,可以得到$mle |V|$

下面,考虑Berlekam-Masey算法——

Berlekamp-Massey算法

记数列${r^{(j)}_{0},r^{(j)}_{1},…,r^{(j)}_{m^{(j)}}}$为${S_{0},S_{1},…,S_{j}}$的最短递推数列

考虑依次计算,先令${r^{(-1)}}={1}$,接下来考虑从${r^{(i-1)}}$推出${r^{(i)}}$

若$sum_{j=0}^{m^{(i-1)}}r^{(i-1)}_{j}S_{i-j}=0$,令${r^{(i)}}={r^{(i-1)}}$即可,否则有以下结论——

结论2:当$sum_{j=0}^{m^{(i-1)}}r^{(i-1)}_{j}S_{i-j}ne 0$,有$m^{(i)}ge max(m^{(i-1)},i+1-m^{(i-1)})$

反证法,即证明当$m^{(i)}<max(m^{(i-1)},i+1-m^{(i-1)})$时,必然有$sum_{j=0}^{m^{(i-1)}}r^{(i-1)}_{j}S_{i-j}=0$

当$m^{(i)}<m^{(i-1)}$时令${r^{(i-1)}}={r^{(i)}}$即矛盾,因此可得$m^{(i)}<i-m^{(i-1)}$

此时,根据${r^{(i)}}$的定义,有$sum_{j=0}^{m^{(i)}}r^{(i)}_{j}S_{i-j}=0$

再根据${r^{(i-1)}}$的定义以及$i-m^{(i)}>m^{(i-1)}$,有$sum_{j=0}^{m^{(i)}}r^{(i)}_{j}sum_{k=0}^{m^{(i-1)}}r^{(i-1)}_{k}S_{i-j-k}=0$

将其交换枚举顺序,有$sum_{k=0}^{m^{(i-1)}}r_{k}^{(i-1)}sum_{j=0}^{m^{(i)}}r_{j}^{(i)}S_{i-j-k}=0$

进而后者根据定义,即$S_{i-k}$,有$sum_{j=0}^{m^{(i-1)}}r_{j}^{(i-1)}S_{i-j}=0$,即所求证

下面,考虑构造${r^{(i)}}$以取到下限$m^{(i)}=max(m^{(i-1)},i+1-m^{(i-1)})$——

若此时是第一次出现此类情况,即有$m^{(i-1)}=0$,此时$m^{(i)}ge i+1$,取$r^{(i)}={1,0,0,…,0}$(其中$0$的个数为$i+1$)

若此时不是第一次,任取$p<i$使得时$sum_{j=0}^{m^{(p-1)}}r^{(p-1)}_{j}S_{p-j}ne 0$,则按如下方式构造${r^{(i)}}$
$$
r^{(i)}_{j}=begin{cases}r^{(i-1)}_{j}&(j<i-p)\r^{(i-1)}_{j}-wr^{(p-1)}_{j-(i-p)}&(i-ple jle max(m^{(i-1)},i-p+m^{(p-1)}))end{cases}
$$
(其中$w=frac{sum_{j=0}^{m^{(i-1)}}r_{j}^{(i-1)}S_{i-j}}{sum_{j=0}^{m^{(p-1)}}r^{(p-1)}_{j}S_{p-j}}$)

关于这一做法的正确性,即要求$forall max(m^{(i-1)},i-p+m^{(p-1)})le kle i$
$$
sum_{j=0}^{i-p-1}r_{j}^{(i-1)}S_{k-j}+sum_{j=i-p}^{i-p+m^{(p-1)}}(r^{(i-1)}_{j}-wr^{(p-1)}_{j-(i-p)})S_{k-j}=0
$$
将其按照${r^{(i-1)}}$和${r^{(p-1)}}$分为两类,分类讨论:

1前者即$sum_{j=0}^{m^{(i-1)}}r_{j}S_{k-j}$,根据定义在$k<i$时为0,在$k=i$时为$sum_{j=0}^{m^{(i-1)}}r_{j}S_{k-j}$

2.后者即$wsum_{j=i-p}^{i-p+m^{(p-1)}}r^{(p-1)}_{j-(i-p)}S_{k-j}=wsum_{j=0}^{m^{(p-1)}}r^{(p-1)}_{j}S_{k-(i-p)-j}$,根据定义在$k-(i-p)<p$(也即$k<i$)时为0,在$k=i$时为$wsum_{j=0}^{m^{(p-1)}}r^{(p-1)}_{j}S_{p-j}=sum_{j=0}^{m^{(i-1)}}r_{j}S_{k-j}$

由此,两者相减在$max(m^{(i-1)},i-p+m^{(p-1)})le kle i$时恒为0,即得证

同时,取$p$为上一次出现此类情况且$m^{(p)}=p+1-m^{(p-1)}$(即尽量大),即有$m^{(p)}=m^{(i-1)}$

进一步的,归纳$m^{(p)}=p+1-m^{(p-1)}$,则$m^{(i)}=i-p+m^{(p-1)}=i+1-m^{(i-1)}$,取到下限

综上,通过Berlekam-Masey算法,对于长度为$n$的数列${S_{0},S_{1},…,S_{n}}$,可以在$o(n^{2})$的时间内求出其最短递推数列

剩下的问题

但${S_{i}}$是无限数列,Berlekamp-Massey仍然不行,考虑通过结论2,还有以下结论——

结论3:对于无限数列${S_{i}}$,若其最短递推数列${r_{0},r_{1},…,r_{m}}$满足$mle n$,则${S_{i}}$的最短递推数列即为${S_{0},S_{1},…,S_{2n-1}}$的最短递推数列

考虑对${S_{0},S_{1},…,S_{2n-1}}$求出其最短递推数列${r’_{0},r’_{1},…,r’_{m’}}$,必然有$m’le m$

同时,若其不是最短递推数列,即存在$ige m’$使得$sum_{j=0}^{m’}r_{j}S_{i-j}ne 0$

取其中最小的$i$,显然有$ige 2n$,而$m^{(i-1)}=m’$,即$m^{(i)}ge i+1-m’>n$,即与$mle n$矛盾

由此,再根据前面$mle |V|$的条件,即只需要求${S_{0},S_{1},…,S_{2|V|-1}}$的最短递推数列即可

由于$M$中非0的元素个数是$o(|V|+|E|)$的,暴力$o(|V|^{2}+|V||E|)$的求出$S_{0},S_{1},…,S_{2|V|-1}$,并通过Berlekamp-Massey算法在$o(|V|^{2})$求出最短递推数列

当求出${S_{i}}$的最短递推数列${r_{0},r_{1},…,r_{m}}$后,即可快速求出$sum_{i=0}^{infty}S_{i}$——

令$F(x)=sum_{i=0}^{infty}S_{i}x^{i}$,$G(x)=sum_{i=0}^{m}r_{i}x^{i}$,即有$(Ftimes G)(x)=sum_{i=0}^{m-1}(sum_{j=0}^{i}r_{j}S_{i-j})x^{i}$(根据前面的性质,$m$次以后都为0)

问题即求$F(1)$,也即$frac{sum_{i=0}^{m-1}sum_{j=0}^{i}r_{j}S_{i-j}}{sum_{i=0}^{m}r_{i}}$,计算复杂度为$o(|V|^{2})$

(特别的,可以认为$sum_{i=0}^{m}r_{i}$在模$P$意义为0概率即$frac{1}{P}$,此概率可以忽略)

综上,我们就得到了一个$o(|V|^{2}+|V||E|)$的做法

另一个做法

考虑直接高斯消元的做法,事实上,方程可以看作一个$|V|times |V|$的矩阵$M$和$1times |V|$的矩阵$b$,并构造出$1times |V|$的矩阵$f$使得$fM=b$

(注意$M$的每一列作为一个方程,而不是每一行,这样可以使得所有矩阵都是$1times |V|$的)

同时右乘$M^{-1}$,即$f=bM^{-1}$,问题即求$M^{-1}$

考虑求出$M$的零化多项式$f(x)=sum_{i=0}^{n}a_{i}x^{i}$(定义见前面),则有$M^{-1}=frac{-sum_{i=1}^{n}a_{i}M^{i-1}}{a_{0}}$

但矩阵的零化多项式如果使用高斯消元+拉格朗日插值法,复杂度为$o(|V|^{4})$,需要优化

事实上,这个问题也即求数列$S_{i}=M^{i}$的最短递推数列

随机一个$1times |V|$的矩阵$A$和$|V|times 1$的矩阵$B$,此时$AM^{i}B$即为一个整数

结论4(Schwartz-Zippel引理):在模$P$意义下,对于$AM^{i}B$的最短递推数列${r_{0},r_{1},…,r_{m}}$,其有$frac{2|V|}{P}$的概率是${S_{i}}$的最短递推数列

由此,考虑随机$A$和$B$并计算$AM^{i}B$,同样根据结论3只需要计算$0le i<2|V|$,且计算完成后再通过Berlekamp-Massey算法即可做到$o(|V|^{2})$的复杂度

关于如何计算$AM^{i}B$,先求出$S’_{i}=AM^{i}=S’_{i-1}M$,由于$M$中非0的元素仅有$o(|E|)$个,因此计算$S’_{i}$的复杂度为$o(|V||E|)$,进而再计算$S’_{i}B$(复杂度为$o(|V|)$),即总复杂度为$o(|V|^{2}+|V||E|)$

求出零化多项式$f(x)=sum_{i=0}^{n}a_{i}x^{i}$后,暴力计算$M^{-1}$复杂度仍为$o(|V|^{3})$,考虑将$b$左乘提到后面,即变为$bM^{-1}=frac{-sum_{i=1}^{n}a_{n-i}bM^{i-1}}{a_{n}}$

(注意要将最短递推数列翻转,即若其为${r_{0},r_{1},…,r_{m}}$,则$f(x)=sum_{i=0}^{m}r_{m-i}x^{i}$)

计算$bM^{i-1}$,类似前面$S’_{i}$也可以做到$o(|V||E|)$的复杂度,再求和复杂度为$o(|V|^{2})$

综上,我们又得到了一个$o(|V|^{2}+|V||E|)$的做法

(事实上,这个做法可以通用于求稀疏矩阵的方程组)

然而,这个做法似乎常数更大,其无法通过

[gym102268E]Expected Value[gym102268E]Expected Value

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 #define N 3005
 4 #define mod 998244353
 5 struct Edge{
 6     int nex,to;
 7 }edge[N*6];
 8 int E,n,m,x,y,ans,head[N],deg[N],S[N<<1],f[N<<1][N],l[N<<1],r[N<<1][N<<1];
 9 int pow(int n,int m){
10     int s=n,ans=1;
11     while (m){
12         if (m&1)ans=1LL*ans*s%mod;
13         s=1LL*s*s%mod;
14         m>>=1;
15     }
16     return ans;
17 }
18 void add(int x,int y){
19     edge[E].nex=head[x];
20     edge[E].to=y;
21     head[x]=E++;
22 }
23 int main(){
24     scanf("%d",&n);
25     for(int i=1;i<=n;i++)scanf("%*d%*d");
26     memset(head,-1,sizeof(head));
27     scanf("%d",&m);
28     for(int i=1;i<=m;i++){
29         scanf("%d%d",&x,&y);
30         if (x!=n){
31             deg[x]++;
32             add(x,y);
33         }
34         if (y!=n){
35             deg[y]++;
36             add(y,x);
37         }
38     }
39     for(int i=1;i<=n;i++)deg[i]=pow(deg[i],mod-2);
40     f[0][1]=S[0]=1;
41     for(int i=1;i<2*n;i++){
42         for(int j=1;j<=n;j++)
43             for(int k=head[j];k!=-1;k=edge[k].nex)
44                 f[i][edge[k].to]=(f[i][edge[k].to]+1LL*deg[j]*f[i-1][j])%mod;
45         for(int j=1;j<n;j++)S[i]=(S[i]+f[i][j])%mod;
46     }
47     int lst=-1;
48     l[0]=0,r[0][0]=1;
49     for(int i=0;i<2*n;i++){
50         int s=0;
51         for(int j=0;j<=l[i];j++)s=(s+1LL*r[i][j]*S[i-j])%mod;
52         if (!s){
53             l[i+1]=l[i];
54             memcpy(r[i+1],r[i],sizeof(r[i]));
55         }
56         else{
57             if (lst<0){
58                 l[i+1]=i+1;
59                 r[i+1][0]=1;
60             }
61             else{
62                 l[i+1]=max(l[i],i+1-l[i]);
63                 memcpy(r[i+1],r[i],sizeof(r[i]));
64                 int ss=0;
65                 for(int j=0;j<=l[lst];j++)ss=(ss+1LL*r[lst][j]*S[lst-j])%mod;
66                 s=1LL*s*pow(ss,mod-2)%mod;
67                 for(int j=i-lst;j<=i-lst+l[lst];j++)r[i+1][j]=(r[i+1][j]-1LL*s*r[lst][j-(i-lst)]%mod+mod)%mod;
68             }
69             if (l[i+1]==i+1-l[i])lst=i;
70         }
71     }
72     for(int i=0;i<l[2*n];i++)
73         for(int j=0;j<=i;j++)ans=(ans+1LL*r[2*n][j]*S[i-j])%mod;
74     int s=0;
75     for(int i=0;i<=l[2*n];i++)s=(s+r[2*n][i])%mod;
76     ans=1LL*ans*pow(s,mod-2)%mod;
77     printf("%d",ans);
78 }

View Code

[gym102268E]Expected Value[gym102268E]Expected Value

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 #define N 3005
 4 #define mod 998244353
 5 vector<int>v[N];
 6 int n,m,x,y,deg[N],M[N][N],A[N],B[N],AA[N],S[N<<1],l[N<<1],r[N<<1][N],ans[N];
 7 int pow(int n,int m){
 8     int s=n,ans=1;
 9     while (m){
10         if (m&1)ans=1LL*ans*s%mod;
11         s=1LL*s*s%mod;
12         m>>=1;
13     }
14     return ans;
15 }
16 int main(){
17     srand(time(0));
18     scanf("%d",&n);
19     for(int i=1;i<=n;i++)scanf("%*d%*d");
20     scanf("%d",&m);
21     for(int i=1;i<=m;i++){
22         scanf("%d%d",&x,&y);
23         if (x!=n){
24             deg[x]++;
25             M[y][x]++;
26         }
27         if (y!=n){
28             deg[y]++;
29             M[x][y]++;
30         }
31     }
32     for(int i=1;i<=n;i++)deg[i]=pow(deg[i],mod-2);
33     for(int i=1;i<=n;i++){
34         for(int j=1;j<=n;j++)M[i][j]=mod-1LL*M[i][j]*deg[j]%mod;
35         M[i][i]=1;
36     }
37     for(int i=1;i<=n;i++)
38         for(int j=1;j<=n;j++)
39             if (M[i][j])v[j].push_back(i);
40     for(int i=1;i<=n;i++){
41         A[i]=1LL*rand()*rand()%mod;
42         B[i]=1LL*rand()*rand()%mod;
43     }
44     for(int i=0;i<2*n;i++){
45         for(int j=1;j<=n;j++)S[i]=(S[i]+1LL*A[j]*B[j])%mod;
46         for(int j=1;j<=n;j++){
47             AA[j]=0;
48             for(int k=0;k<v[j].size();k++)AA[j]=(AA[j]+1LL*A[v[j][k]]*M[v[j][k]][j])%mod;
49         }
50         memcpy(A,AA,sizeof(A));
51     }
52     int lst=-1;
53     l[0]=0,r[0][0]=1;
54     for(int i=0;i<2*n;i++){
55         int s=0;
56         for(int j=0;j<=l[i];j++)s=(s+1LL*r[i][j]*S[i-j])%mod;
57         if (!s){
58             l[i+1]=l[i];
59             memcpy(r[i+1],r[i],sizeof(r[i]));
60         }
61         else{
62             if (lst<0){
63                 l[i+1]=i+1;
64                 r[i+1][0]=1;
65             }
66             else{
67                 l[i+1]=max(l[i],i+1-l[i]);
68                 memcpy(r[i+1],r[i],sizeof(r[i]));
69                 int ss=0;
70                 for(int j=0;j<=l[lst];j++)ss=(ss+1LL*r[lst][j]*S[lst-j])%mod;
71                 s=1LL*s*pow(ss,mod-2)%mod;
72                 for(int j=i-lst;j<=i-lst+l[lst];j++)r[i+1][j]=(r[i+1][j]-1LL*s*r[lst][j-(i-lst)]%mod+mod)%mod;
73             }
74             if (l[i+1]==i+1-l[i])lst=i;
75         }
76     }
77     for(int i=1;i<n;i++)A[i]=1;
78     A[n]=0;
79     for(int i=1;i<=l[2*n];i++){
80         for(int j=1;j<=n;j++)ans[j]=(ans[j]+1LL*r[2*n][l[2*n]-i]*A[j])%mod;
81         for(int j=1;j<=n;j++){
82             AA[j]=0;
83             for(int k=0;k<v[j].size();k++)AA[j]=(AA[j]+1LL*A[v[j][k]]*M[v[j][k]][j])%mod;
84         }
85         memcpy(A,AA,sizeof(A));
86     }
87     int x=mod-pow(r[2*n][l[2*n]],mod-2);
88     for(int i=1;i<=n;i++)ans[i]=1LL*ans[i]*x%mod;
89     printf("%d",ans[1]);
90 }

View Code

 


程序员灯塔
转载请注明原文链接:[gym102268E]Expected Value
喜欢 (0)