LOJ#3058. 「HNOI2019」白兔之舞
LOJ#3058. 「HNOI2019」白兔之舞
分析
-
不妨设转移矩阵为\(f\)，且用\(f^i\)来表示走\(i\)步从\(x\)到\(y\)的方案数。
-
跳了\(i\)步的方案数为\(\sum\limits_{j=i}^L\binom{j-1}{i-1}=\binom{L}{i}\)
-
直接推式子\(ans_t=\sum\limits_{i=0}^L[i\bmod k=t]f^i\binom{L}{i}\)
-
发现这是个挺套路的单位根反演
-
\(ans_t\times k=\sum\limits_{i=0}^{L}\sum\limits_{j=0}^{k-1}w^{j(i-t)}f^i\binom{L}{i}\)
- $ =\sum\limits_{j=0}{k-1}w\sum\limits_{i=0}Lwf^i\binom{L}{i}$
- \(=\sum\limits_{j=0}^{k-1}w^{-jt}(w^jf+I)\)
-
把乘积化成和的形式，有\(ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}\)
-
\(ans_t\times k=\sum\limits_{j=0}^{k-1}w^{-\binom{j+t}{2}+\binom{j}{2}+\binom{t}{2}}g_j\)
-
\(ans_t\times k\times w^{-\binom{t}{2}}=\sum\limits_{j=0}^{k-1}w^{-\binom{j+t}{2}}w^{\binom{j}{2}}g_j\)
-
是一个卷积的形式，直接用\(mtt\)做就行了，用https://cmxrynp.github.io/2019/01/07/fft-optimization/中的\(4\)次或\(3.5\)次dft再加上一点点卡常即可通过此题。
代码
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;
typedef long double f2;
#define N 262190
#define db(x) cerr<<#x<<" = "<<x<<endl
namespace groot {
int pri[100],cnt;
ll qp(ll x,ll y,ll p) {
ll re=1;for(;y;y>>=1,x=x*x%p)if(y&1)re=re*x%p; return re;
}
bool check(int x,int phi) {
int i;
for(i=1;i<=cnt;i++) if(qp(x,phi/pri[i],phi+1)==1) return 0;
return 1;
}
int get_g(int p,int phi) {
int x=phi,i;
for(i=2;i*i<=x;i++) {
if(x%i==0) {
pri[++cnt]=i;
while(x%i==0) x/=i;
}
}
if(x!=1) pri[++cnt]=x;
int g=2;while(!check(g,phi)) g++;
return g;
}
}
int mod,n,K,L,X,Y,G;
ll g[N],w[N],f[N];
ll qp(ll x,ll y) {
ll re=1; for(;y;y>>=1,x=x*x%mod) if(y&1) re=re*x%mod; return re;
}
ll C2(ll x) {return x*(x-1)/2%K;}
namespace matrix {
struct Mat {
ll a[3][3];
Mat() {memset(a,0,sizeof(a));}
Mat operator * (const Mat &u) const {
int i,j,k; Mat re;
for(i=0;i<n;i++)for(k=0;k<n;k++)for(j=0;j<n;j++) {
re.a[i][j]+=a[i][k]*u.a[k][j];
}for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]%=mod;
return re;
}
Mat operator * (const ll &x) const {
Mat re=*this; int i,j;
for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]=re.a[i][j]*x%mod;
return re;
}
Mat operator + (const Mat &u) const {
Mat re=*this; int i,j;
for(i=0;i<n;i++)for(j=0;j<n;j++)re.a[i][j]=(re.a[i][j]+u.a[i][j])%mod;
return re;
}
}I,O;
Mat qp(Mat x,int y) {
Mat re=I;
for(;y;y>>=1,x=x*x)if(y&1)re=re*x; return re;
}
void solve() {
int i;
for(i=0;i<n;i++) I.a[i][i]=1;
for(i=0;i<K;i++) {
g[i]=qp(O*w[i]+I,L).a[X-1][Y-1]*w[C2(i)]%mod;
}
}
}
ll tmp[N];
namespace fft {
const f2 pi=acos(-1);
struct cp {
f2 x,y;
cp() {}
cp(f2 x_,f2 y_) {x=x_,y=y_;}
cp operator + (const cp &u) const {return cp(x+u.x,y+u.y);}
cp operator - (const cp &u) const {return cp(x-u.x,y-u.y);}
cp operator * (const cp &u) const {return cp(x*u.x-y*u.y,x*u.y+y*u.x);}
}A[N],B[N],C[N],D[N];
void dft(cp *a,int len) {
int i,j,k,t; cp tmp,w,wn;
for(i=k=0;i<len;i++) {
if(i>k) swap(a[i],a[k]);
for(j=len>>1;(k^=j)<j;j>>=1);
}
for(k=2;k<=len;k<<=1) {
wn=cp(cos(2*pi/k),sin(2*pi/k));
for(t=k>>1,i=0;i<len;i+=k) {
for(w=cp(1,0),j=i;j<i+t;j++) {
tmp=a[j+t]*w; a[j+t]=a[j]-tmp; a[j]=a[j]+tmp; w=w*wn;
}
}
}
}
int solve() {
int i;
w[K]=1;
for(i=0;i<K+K;i++) f[i]=w[K-C2(i)];
reverse(g,g+K+1);
int n=K+K-1,m=K;
int l=1;
while(l<=(n+m))l<<=1;
//for(i=0;i<l;i++) printf("%lld %lld\n",f[i],g[i]);
for(i=0;i<=n;i++) {
A[i].x=f[i]>>15;
A[i].y=f[i]&32767;
}
for(i=0;i<=m;i++) {
B[i].x=g[i]>>15;
B[i].y=g[i]&32767;
}
dft(A,l),dft(B,l);
for(i=0;i<l;i++) {
int pl=(l-1)&(l-i);
const cp ca=cp(A[pl].x,-A[pl].y),cb=cp(B[pl].x,-B[pl].y);
const cp a=(A[i]+ca)*cp(0.5,0),b=(A[i]-ca)*cp(0,-0.5),c=(B[i]+cb)*cp(0.5,0),d=(B[i]-cb)*cp(0,-0.5);
C[pl]=a*c+a*d*cp(0,1);
D[pl]=b*c+b*d*cp(0,1);
}
dft(C,l),dft(D,l);
ll invk=qp(K,mod-2);
for(i=0;i<l;i++) C[i].x/=l,C[i].y/=l,D[i].x/=l,D[i].y/=l;
//for(i=0;i<=n;i++) for(int j=0;j<=m;j++) {
//tmp[i+j]=(tmp[i+j]+f[i]*g[j])%mod;
//}
for(i=K;i<K+K;i++) {
ll a=C[i].x+0.5,b=C[i].y+0.5,c=D[i].x+0.5,d=D[i].y+0.5;
a=((a%mod)<<30)+(((b+c)%mod)<<15)+d;
a=(a%mod+mod)%mod;
//a=tmp[i];
//db(a);
ll t=w[C2(i-K)]*invk%mod;
printf("%lld\n",a*t%mod);
}
return 0;
}
}
int main() {
scanf("%d%d%d%d%d%d",&n,&K,&L,&X,&Y,&mod);
G=groot::get_g(mod,mod-1);
w[0]=1;w[1]=qp(G,(mod-1)/K);
int i,j;
for(i=2;i<K;i++) w[i]=w[i-1]*w[1]%mod;
//printf("%d\n",G);
for(i=0;i<n;i++)for(j=0;j<n;j++)scanf("%lld",&matrix::O.a[i][j]);
matrix::solve();
return fft::solve();
}
