标签：直接 ret 乘法 color mem i++ 随机 相减 namespace
令$X$为移动次数，答案即$\sum_{i=0}^{\infty}P(X>i)$，后者记作$S_{i}$
关于$S_{i}$，令$f_{i,j}$表示走了$i$步后位于$j$且未到达过$k$的概率，即有$S_{i}=\sum_{j\in V,j\ne 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}\}$的最短递推数列：满足$m\ge 1$，$r_{0}\ne 0$且$\forall m\le j\le 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 j\ge 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 j\ge 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})$
（另外，我们要求的是最短递推数列，其也不满足此性质，虽然这并不重要）
虽然不能直接计算，但从这个结论中，可以得到$m\le |V|$
下面，考虑Berlekam-Masey算法——
记数列$\{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-p\le j\le \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 k\le 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.后者即$w\sum_{j=i-p}^{i-p+m^{(p-1)}}r^{(p-1)}_{j-(i-p)}S_{k-j}=w\sum_{j=0}^{m^{(p-1)}}r^{(p-1)}_{j}S_{k-(i-p)-j}$，根据定义在$k-(i-p)<p$（也即$k<i$）时为0，在$k=i$时为$w\sum_{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 k\le 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}\}$满足$m\le 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$
同时，若其不是最短递推数列，即存在$i\ge m‘$使得$\sum_{j=0}^{m‘}r_{j}S_{i-j}\ne 0$
取其中最小的$i$，显然有$i\ge 2n$，而$m^{(i-1)}=m‘$，即$m^{(i)}\ge i+1-m‘>n$，即与$m\le n$矛盾
由此，再根据前面$m\le |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}$，即有$(F\times 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$和$1\times |V|$的矩阵$b$，并构造出$1\times |V|$的矩阵$f$使得$fM=b$
（注意$M$的每一列作为一个方程，而不是每一行，这样可以使得所有矩阵都是$1\times |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}$的最短递推数列
随机一个$1\times |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只需要计算$0\le 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|)$的做法
（事实上，这个做法可以通用于求稀疏矩阵的方程组）
然而，这个做法似乎常数更大，其无法通过
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 }
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 }
标签：直接 ret 乘法 color mem i++ 随机 相减 namespace
原文地址：https://www.cnblogs.com/PYWBKTDA/p/14723438.html
