首先,考虑直接DP,显然我们需要求出两个数组:f[i][j]:走一条路,与出来是i,步数为j的方案数;g[i][j]:走若干条路,与出来是i,步数为j的方案数。
显然f可以直接通过矩阵乘法n^3m做出来,但是这样并不能通过,因此我们考虑使用特征多项式来优化,原理就是我们可以用[0,n-1]次的矩阵的线性组合来表示出n次的,然后对于n+t次的矩阵,我们只需要给等式的右边乘上这个矩阵就好了。对于这道题,如果我们枚举了一个与出来的值i,那么我们可以搞出来一个行(列)向量以及它的转置,然后给等式左边和右边乘上这个向量及转置,之后就有递推式了,形式化一点,设\(B_{n+i} = X A^{n+i} X^{T}\),那么有\(B_{n+i}=\sum _{j=0}^{n-1}h_{j}X A^{j+i} X^{T} = \sum _{j=0} ^{n-1} h_{j} B{j+i} \),之后,我们可以暴力O(n^4)算出前n个矩阵,然后O(n^2m)递推剩下的。
对于特征多项式的求解,由于这道题并不是什么常系数递推,我们只能暴力求解,具体就是用定义式\(char(A) = | A - \lambda E | \),带入不同的λ,求解行列式的值,然后高斯消元得到系数h,复杂度O(n^4+n^3)。
之后我们就能求出f了,考虑如何求g。我们可以先对第一维进行FWT,之后对第二维进行NTT,这样单次卷积的复杂度就是O(nm)了,考虑到最后的结果其实就是k次卷积之和,我们可以用等比数列求和公式转化到多项式求逆上,然后可以O(nmlogm)的复杂度得到结果。考试的时候想的是用倍增来做,不过这样需要每次都DFT回去进行取模,复杂度好像是两个log的。
然后我们只需要卡卡常数就可以过去了,对于32位的机子,只有在卡取模的常数的时候才要用long long,其他的还是全程int吧。多项式取模现算长度应该会快一点?
总的来说,考试的时候只会\(O(n^3 m + n m log^2 m logn)\)的做法,但是觉得前面的n^3m会让这个算法直接沦为最弱的暴力,于是就直接写了\(O(n^3 m + n^2 m^2)\)的,成功被subtask卡掉了。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 |
#pragma GCC optimize("-Ofast") #pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native") #include <bits/stdc++.h> #define ll long long #define mod 998244353 #define P 998244353 #define mxl 65536 using namespace std; inline int read(){ int s=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-')f=-1; ch=getchar();} while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar(); return s*f; } inline ll qpow(ll a,int b=mod-2){ ll ret=1; while(b){ if(b&1) ret=ret*a%mod; a=a*a%mod,b>>=1; } return ret; } int rev[70005],G[2][70005]; inline void ntt_init(){ register ll w=1,wn=qpow(3,(P-1)/mxl); for(int i=0;i<=mxl;++i) G[0][i]=w,w=w*wn%mod; } inline void ntt_prep(int n){ for(int i=0;i<n;++i) { if(i&1) rev[i]=(rev[i>>1]>>1)|(n>>1); else rev[i]=rev[i>>1]>>1; } } inline void ntt(int *s,int n,int o){ register int i,j,l,m; register int t,u; for(int i=0;i<n;++i) if(i>rev[i]) swap(s[i],s[rev[i]]); for(l=1,m=2;m<=n;l=m,m<<=1) { for(i=0;i<l;++i) G[1][i]=G[0][o==1?(mxl/m*i):(mxl-mxl/m*i)]; for(i=0;i<n;i+=m) { for(j=0;j<l;++j) { t=s[i+j],u=1ll*G[1][j]*s[i+j+l]%mod,s[i+j]=(t+u)%mod,s[i+j+l]=(t-u+mod)%mod; } } } if(o==-1) { register ll inv=qpow(n); for(i=0;i<n;++i) s[i]=1ll*s[i]*inv%mod; } } int tp[70005],pf[70005],pg[70005]; inline int glen(int x){ int ret=1; while(ret<x) ret<<=1; return ret; } inline void pl_inv(int *a,int *b,int M){ if(M==1) {b[0]=qpow(a[0]); return;} register int N=glen(M); pl_inv(a,b,(M+1)>>1); memcpy(tp,a,sizeof(tp[0])*N),memset(tp+N,0,sizeof(tp[0])*N); ntt_prep(N<<1); ntt(tp,N<<1,1),ntt(b,N<<1,1); for(int i=0;i<(N<<1);++i) tp[i]=1ll*b[i]*(2ll-(1ll*tp[i]*b[i])%mod)%mod,tp[i]<0?tp[i]+=mod:0; ntt(tp,N<<1,-1); memcpy(b,tp,sizeof(tp[0])*N),memset(b+N,0,sizeof(b[0])*N); } #define add(a,b) (((a)+=(b)),((a)>=mod?(a)-=mod:0)) #define del(a,b) (((a)-=(b)),((a)<0?(a)+=mod:0)) int n,m,bit[32]; struct matrix{ int s[65][65]; matrix(){memset(s,0,sizeof(s));} }A[67],S; ll tmp[65][65]; inline matrix operator * (const matrix &a,const matrix &b) { matrix c; for(int i=0;i<n;++i) for(int j=0;j<n;++j) tmp[i][j]=0; for(int i=0;i<n;++i) { for(int k=0;k<n;++k) if(a.s[i][k]) { for(int j=0;j<n;++j) if(b.s[k][j]) { tmp[i][j]+=1ll*a.s[i][k]*b.s[k][j]; tmp[i][j]>=8000000000000000000ll?tmp[i][j]%=mod:0; } } } for(int i=0;i<n;++i) for(int j=0;j<n;++j) c.s[i][j]=tmp[i][j]%mod; return c; } inline matrix operator + (const matrix &a,const matrix &b) { matrix c; for(int i=0;i<n;++i) for(int j=0;j<n;++j) c.s[i][j]=a.s[i][j]+b.s[i][j],c.s[i][j]>=mod?c.s[i][j]-=mod:0; return c; } int gs[75][75],h[75]; inline void gauss(){ register int i,j,k,r; for(i=0;i<=n;++i) { r=i; for(j=i+1;j<=n;++j) if(gs[j][i]>gs[r][i]) r=j; if(!gs[r][i]) continue; if(r!=i) for(j=i;j<=n+1;++j) swap(gs[r][j],gs[i][j]); register ll inv=qpow(gs[i][i]); for(j=i+1;j<=n;++j) { register ll tm=1ll*inv*gs[j][i]%mod; for(k=i;k<=n+1;++k) gs[j][k]=(gs[j][k]-1ll*gs[i][k]*tm%mod+mod)%mod; } } for(i=n;i>=0;--i) { for(j=n;j>i;--j) gs[i][n+1]+=mod-((1ll*gs[i][j]*h[j])%mod),gs[i][n+1]>=mod?gs[i][n+1]-=mod:0; gs[i][n+1]=h[i]=1ll*gs[i][n+1]*qpow(gs[i][i])%mod; } } int dt[65][65]; inline ll det(int N,int x){ int ret=1; for(int i=0;i<N;++i) for(int j=0;j<N;++j) { dt[i][j]=S.s[i][j]; if(i==j) dt[i][j]-=x; } for(int i=0;i<N;++i) { if(dt[i][i]<0) { ret=-ret; for(int j=i;j<N;++j) dt[i][j]=-dt[i][j]; } for(int j=i+1;j<N;++j) { while(dt[j][i]) { if(dt[j][i]<0) { ret=-ret; for(int k=i;k<n;++k) dt[j][k]=-dt[j][k]; } int t=dt[i][i]/dt[j][i]; for(int k=i;k<N;++k) dt[i][k]=(dt[i][k]-1ll*dt[j][k]*t)%mod,swap(dt[i][k],dt[j][k]); ret=-ret; } } ret=1ll*ret*dt[i][i]%mod; } return ret; } ll g[65][70005],f[65][70005]; int t[65][70005],d[65][70005]; int len=1; signed main(){ ntt_init(); #ifndef ONLINE_JUDGE freopen("4_10.in","r",stdin); #endif for(int i=0;i<=30;++i) bit[i]=1<<i; n=read(),m=read(); while(len<=(m+1)) len<<=1; for(int i=0;i<n;++i) A[0].s[i][i]=1; for(int i=0;i<n;++i) for(int j=0;j<n;++j) S.s[i][j]=read(); for(int i=1;i<=n;++i) A[i]=A[i-1]*S; for(int i=0;i<=n;++i) { ll pw=1; for(int j=0;j<=n;++j) gs[i][j]=pw,pw=pw*(i+1)%mod; gs[i][n+1]=(det(n,i+1)%mod+mod)%mod; } gauss(); for(int j=0;j<n;++j) h[j]=mod-h[j]; for(int i=0;i<=n;++i) for(int j=0;j<n;++j) for(int k=0;k<n;++k) (g[j&k][i]+=A[i].s[j][k])%=mod; for(int t=1;n+t<=m;++t) { for(int i=0;i<n;++i) { for(int j=0;j<n;++j) { g[i][n+t]+=1ll*h[j]*g[i][j+t]; g[i][n+t]>=1000000000000000000ll?g[i][n+t]%=mod:0; } g[i][n+t]%=mod; } } for(int i=0;i<n;++i) g[i][0]=0; memcpy(f,g,sizeof(f)); for(int L=1,M=2;M<=n;L=M,M<<=1) { for(int i=0;i<n;i+=M) { for(int j=0;j<L;++j) { for(int k=1;k<=m;++k) add(f[i+j][k],f[i+j+L][k]); } } } //cout<<(double)clock()/CLOCKS_PER_SEC<<endl; for(int i=0;i<n;++i) { memset(pf,0,sizeof(pf)),memset(pg,0,sizeof(pg)); pf[0]++; for(int j=1;j<=m;++j) pf[j]=mod-f[i][j]; pl_inv(pf,pg,m+1); for(int j=1;j<=m;++j) d[i][j]=pg[j]; } for(int L=1,M=2;M<=n;L=M,M<<=1) { for(int i=0;i<n;i+=M) { for(int j=0;j<L;++j) { for(int k=1;k<=m;++k) del(d[i+j][k],d[i+j+L][k]); } } } //cout<<(double)clock()/CLOCKS_PER_SEC<<endl; int ans=0; for(int i=0;i<n;++i) for(int j=1;j<=m;++j) ans^=d[i][j]; printf("%d",ans); } |
叨叨几句... NOTHING