【LOJ3058】【HNOI2019】白兔之舞_综合

Description

https://loj.ac/problem/3058

Solution

首先答案长这样子： $ans_t=\sum\limits_{i=0}^L[k|(i-t)]A^i\binom{L}{i}$ ， $A$ 是读入的矩阵，最后取 $(x, y)$ 的值
套上单位根反演就是： $\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}$ ，其中 $\omega_k=g^{\frac{p-1}{k}}$ ， $g$ 是原根。
交换主体再化简： $\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)^L$
然后把 $? t j$ 拆开： $-tj=\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}$ 。
最后式子就是： $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}}$
把 $\omega_k^{\binom{t+j}{2}}$ 和 $ans_t$ 翻转一下就是一个卷积形式了。

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]);}
}