【LOJ3058】【HNOI2019】白兔之舞

Description

https://loj.ac/problem/3058

Solution

首先答案长这样子：anst=∑i=0L[k∣(i−t)]Ai(Li)ans_t=\sum\limits_{i=0}^L[k|(i-t)]A^i\binom{L}{i}anst=i=0∑L[k∣(i−t)]Ai(iL)，AAA是读入的矩阵，最后取(x,y)(x,y)(x,y)的值
套上单位根反演就是：∑i=0LAx,yi(Li)1k∑j=0k−1ωk(i−t)j\sum\limits_{i=0}^LA^i_{x,y}\binom{L}{i}\dfrac{1}{k}\sum\limits_{j=0}^{k-1}\omega_k^{(i-t)j}i=0∑LAx,yi(iL)k1j=0∑k−1ωk(i−t)j，其中ωk=gp−1k\omega_k=g^{\frac{p-1}{k}}ωk=gkp−1，ggg是原根。
交换主体再化简：1k∑j=0k−1ωk−tj∑i=0L(Li)ωkijAx,yi=1k∑j=0k−1ωk−tj(ωkjA+1)L\dfrac{1}{k}\sum\limits_{j=0}^{k-1}\omega_k^{-tj}\sum\limits_{i=0}^L\binom{L}{i}\omega_k^{ij}A_{x,y}^i=\dfrac{1}{k}\sum\limits_{j=0}^{k-1}\omega_k^{-tj}(\omega_{k}^jA+1)^Lk1j=0∑k−1ωk−tji=0∑L(iL)ωkijAx,yi=k1j=0∑k−1ωk−tj(ωkjA+1)L
然后把−tj-tj−tj拆开：−tj=(t2)+(j2)−(t+j2)-tj=\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}−tj=(2t)+(2j)−(2t+j)。
最后式子就是：anst=1kωk(t2)∑j=0k−1ωk(j2)(ωkjA+1)Lωk(t+j2)ans_t=\dfrac{1}{k}\omega_k^{\binom{t}{2}}\sum\limits_{j=0}^{k-1}\omega_k^{\binom{j}{2}}(\omega_k^jA+1)^L\omega_k^{\binom{t+j}{2}}anst=k1ωk(2t)j=0∑k−1ωk(2j)(ωkjA+1)Lωk(2t+j)
把ωk(t+j2)\omega_k^{\binom{t+j}{2}}ωk(2t+j)和anstans_tanst翻转一下就是一个卷积形式了。

Code

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,j,k) for(int i=j;i<=k;++i)
#define fd(i,j,k) for(int i=j;i>=k;--i)
using namespace std;
typedef long long ll;
const int N=1e5+10,M=3e5+10,qmo=31622;
int n,P;
int w[N];
int qpow(int x,int y){int s=1;for(;y;y>>=1,x=(ll)x*x%P) if(y&1) s=(ll)s*x%P;return s;
}
void inc(int &x,int y){x=x+y>=P?x+y-P:x+y;
}
struct mat{int a[4][4];void clear(){memset(a,0,sizeof(a));fo(i,1,n) a[i][i]=1;}friend mat operator *(mat x,mat y){mat z;memset(z.a,0,sizeof(z.a));//z.a[1] z.a[2] z.a[3]fo(i,1,n)fo(j,1,n){ll t=0;fo(k,1,n) t+=(ll)x.a[i][k]*y.a[k][j];z.a[i][j]=t%P;}return z;}void mul(int x){fo(i,1,n)fo(j,1,n) a[i][j]=(ll)a[i][j]*x%P;}
}A,S,z;
void matpow(int y){S.clear();for(;y;y>>=1,A=A*A)if(y&1) S=S*A;
}
int pr[50];
int getrt(){int x=P-1;for(int i=2;i*i<=x;++i) if(x%i==0){pr[++pr[0]]=i;for(;x%i==0;x/=i);}fo(i,2,P-1){bool flag=1;fo(j,1,pr[0]) if(qpow(i,(P-1)/pr[j])==1){flag=0;break;}if(flag) return i;}
}
struct node{double x,y;node(double _x=0,double _y=0) {x=_x,y=_y;}
}a[M],b[M],c1[M],c2[M],c3[M],W[M];
const double pi=acos(-1);
node operator+(node x,node y) {return node(x.x+y.x,x.y+y.y);}
node operator-(node x,node y) {return node(x.x-y.x,x.y-y.y);}
node operator*(node x,node y) {return node(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
node conj(node a) {return node(a.x,-a.y);}
int rev[M],g[M];
void DFT(node *a,int fn,int sig) {fo(i,0,fn-1) if(i<rev[i]) swap(a[i],a[rev[i]]);for(int m=2;m<=fn;m<<=1){int half=m/2,t=fn/m;fo(i,0,half-1){node w=sig>0?W[t*i]:W[fn-t*i];for(int j=i;j<fn;j+=m){node u=a[j],v=w*a[j+half];a[j]=u+v,a[j+half]=u-v;}}}if(sig==-1) fo(i,0,fn-1) a[i].x/=fn;
}
void mul(int *A,int *B,int nl,int nr){int fn,cnt=0;for(fn=1;fn<=nl+nr;fn<<=1) ++cnt;fo(i,0,fn) W[i]=node(cos(2*i*pi/fn),sin(2*i*pi/fn));fo(i,1,fn-1) rev[i]=rev[i>>1]>>1|(i&1)<<(cnt-1);fo(i,0,nl) a[i]=node(A[i]/qmo,A[i]%qmo);fo(i,nl+1,fn-1) a[i]=node(0,0);fo(i,0,nr) b[i]=node(B[i]/qmo,B[i]%qmo);fo(i,nr+1,fn-1) b[i]=node(0,0);DFT(a,fn,1),DFT(b,fn,1);fo(i,0,fn-1){int j=(fn-i)&(fn-1);node p1,p2,q1,q2;p1=(a[i]+conj(a[j]))*node(0.5,0);p2=(a[i]-conj(a[j]))*node(0,-0.5);q1=(b[i]+conj(b[j]))*node(0.5,0);q2=(b[i]-conj(b[j]))*node(0,-0.5);c1[i]=p1*q1,c2[i]=p1*q2+p2*q1,c3[i]=p2*q2;}DFT(c1,fn,-1),DFT(c2,fn,-1),DFT(c3,fn,-1);fo(i,0,fn-1){g[i]=((ll)(c1[i].x+0.5)*qmo%P*qmo%P+(ll)(c2[i].x+0.5)*qmo%P+(ll)(c3[i].x+0.5))%P;}fo(i,0,nl+nr) A[i]=g[i];fo(i,nl+nr+1,fn-1) A[i]=0;
}
int k;
int C2(int x){return (ll)x*(x-1)/2%k;
}
int c[M],d[M];
int main()
{freopen("dance.in","r",stdin);freopen("dance.out","w",stdout);int L,x,y;scanf("%d %d %d %d %d %d",&n,&k,&L,&x,&y,&P);int nk=qpow(k,P-2);w[0]=1,w[1]=qpow(getrt(),(P-1)/k);fo(i,2,k-1) w[i]=(ll)w[i-1]*w[1]%P;fo(i,1,n)fo(j,1,n) scanf("%d",&z.a[i][j]);fo(i,0,k-1){A=z;A.mul(w[i]);fo(j,1,n) inc(A.a[j][j],1);matpow(L);c[i]=(ll)S.a[x][y]*w[C2(i)]%P;}fo(i,0,2*k-2) d[i]=qpow(w[C2(i)],P-2);/*fo(i,0,k-1){int tmp=0;fo(j,0,k-1) inc(tmp,(ll)c[j]*d[i+j]%P);tmp=(ll)tmp*nk%P*w[C2(i)]%P;printf("%d\n",tmp);}*/reverse(d,d+2*k-1);mul(c,d,k-1,2*k-2);reverse(c,c+2*k-1);fo(i,0,k-1){c[i]=(ll)c[i]*nk%P*w[C2(i)]%P;printf("%d\n",c[i]);}
}