#3058. 「HNOI2019」白兔之舞(单位根反演)
题目描述
单位根反演
因为ω太难写了所以用w代替
有\([n|k]=\frac{1}{n}\sum_{i=0}^{n-1} w_n^{ik}\)
证明:
当n|k时显然是1,否则\(\frac{w_n^{nk}-1}{w^n-1}=0\)
题解
一开始想矩乘存多项式然后快速幂循环卷积,然后多乘了一个log直接炸飞
先暴力找p-1的原根,然后即可求出k次单位根w,A是读入的矩阵
\(ans_t=\sum_{i=0}^L [k|(i-t)]A_{x,y}^i \binom{L}{i}\)
\(=\frac{1}{k}\sum_{i=0}^L \sum_{j=0}^{k-1}w^{j(i-t)} A_{x,y}^i \binom{L}{i}\)
\(=\frac{1}{k}\sum_{j=0}^{k-1}w^{-jt} \sum_{i=0}^L w^{ij} A_{x,y}^i \binom{L}{i}\)
\(=\frac{1}{k}\sum_{j=0}^{k-1}w^{-jt} (w^jA+1)^L_{x,y}\)
\(=\frac{1}{k}\sum_{j=0}^{k-1}w^{\binom{j}{2}}w^{\binom{t}{2}}w^{-\binom{j+t}{2}} (w^jA+1)^L_{x,y}\)
把j变成k-1-j即可卷积,要写mtt
注意mtt求a-bi的点值时的关系是共轭而不是相等
code
#include <bits/stdc++.h>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define add(a,b) a=((a)+(b))%mod
#define ll long long
#define D 10000
//#define file
using namespace std;
int N,len,n,K,L,x,y,mod,Mod,i,j,k,l;
ll a[262144],b[262144],g,w,s;
struct type{double x,y;} AA[262144];
type operator + (type a,type b) {return {a.x+b.x,a.y+b.y};}
type operator - (type a,type b) {return {a.x-b.x,a.y-b.y};}
type operator * (type a,type b) {return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
struct mat{ll a[3][3]; void clear() {memset(a,0,sizeof(a));}} I,A,B;
mat operator * (mat a,mat b)
{
static mat c;
int i,j,k;
fo(i,0,n-1)
{
fo(j,0,n-1)
{
c.a[i][j]=0;
fo(k,0,n-1)
add(c.a[i][j],a.a[i][k]*b.a[k][j]);
}
}
return c;
}
mat Qpower(mat a,ll b) {mat ans=I; while (b) {if (b&1) ans=ans*a;a=a*a;b>>=1;} return ans;}
ll qpower(ll a,ll b) {ll ans=1; while (b) {if (b&1) ans=ans*a%mod;a=a*a%mod;b>>=1;} return ans;}
void dft(type *a,int tp)
{
int i,j,k,l,S=N,s1=2,s2=1;
type w,W,u,v;
fo(i,0,N-1)
{
j=i,k=0;
fo(l,1,len)
k=k*2+(j&1),j>>=1;
AA[i]=a[k];
}
memcpy(a,AA,N*sizeof(type));
fo(i,1,len)
{
S>>=1;
fo(j,0,S-1)
{
fo(k,0,s2-1)
{
W={cos(2*M_PI*k/s1),sin(2*M_PI*k/s1)*tp};
u=a[j*s1+k],v=a[j*s1+k+s2]*W;
a[j*s1+k]=u+v;
a[j*s1+k+s2]=u-v;
}
}
s1<<=1,s2<<=1;
}
}
void mul(ll *a,ll *b)
{
static type A[262144],B[262144],P[262144],Q[262144];
static ll ans[262144];
int i,j,k,l;
fo(i,0,N-1) A[i]={a[i]/D,a[i]%D},B[i]={b[i]/D,b[i]%D};
dft(A,1),dft(B,1);
fo(i,0,N-1) P[i]=A[i]*B[i],Q[i]=type{A[(N-i)%N].x,-A[(N-i)%N].y}*B[i];
dft(P,-1),dft(Q,-1);
fo(i,0,N-1) P[i].x/=N,Q[i].x/=N,P[i].y/=N,Q[i].y/=N;
fo(i,0,N-1)
ans[i]=(((ll)floor((P[i].x+Q[i].x)/2+0.5)%mod*D%mod*D)+(ll)floor(P[i].y+0.5)*D+(ll)floor((Q[i].x-P[i].x)/2+0.5))%mod;
memcpy(a,ans,N*8);
}
bool pd(ll g)
{
int i,j,k,l,s=floor(sqrt(mod-1));
fo(i,1,s)
if (!((mod-1)%i))
{
if (qpower(g,i)==1) return 0;
if (i>1 && qpower(g,(mod-1)/i)==1) return 0;
}
return 1;
}
int main()
{
#ifdef file
freopen("loj3058.in","r",stdin);
#endif
scanf("%d%d%d%d%d%d",&n,&K,&L,&x,&y,&mod),--x,--y,len=ceil(log2(K*3)),N=pow(2,len),Mod=mod-2;
I.a[0][0]=I.a[1][1]=I.a[2][2]=1;
fo(i,0,n-1)
{
fo(j,0,n-1)
scanf("%lld",&A.a[i][j]);
}
for (g=1; !pd(g); ++g);
w=qpower(g,(mod-1)/K);
fo(i,0,K*2-1) a[i]=qpower(qpower(w,1ll*i*(i-1)/2),Mod);
fo(i,0,K-1)
{
B=A,s=qpower(w,i);
fo(j,0,n-1)
{
fo(k,0,n-1)
B.a[j][k]=B.a[j][k]*s%mod+(j==k);
}
B=Qpower(B,L);
b[(K-1)-i]=qpower(w,1ll*i*(i-1)/2)*B.a[x][y]%mod;
}
mul(a,b);
fo(i,0,K-1)
printf("%lld\n",(a[i+(K-1)]*qpower(K,Mod)%mod*qpower(w,1ll*i*(i-1)/2)%mod+mod)%mod);
fclose(stdin);
fclose(stdout);
return 0;
}