「HNOI 2019」白兔之舞
作者：互联网
一道清真的数论题
题解
考虑$ n=1$的时候怎么做
设$ s$为转移的方案数
设答案多项式为$\sum\limits_{i=0}^L (sx)^i\binom{L}{i}=(sx+1)^L$
答案相当于这个多项式模$ k$的各项系数的和
发现这和LJJ学二项式定理几乎一模一样
然而直接搞是$ k^2$的，无法直接通过本题
以下都用$ w$表示$ k$次单位根
设$ F_i$为次数模$ k$为$ i$的项的系数和
单位根反演一下得到$F(i)=\sum\limits_{j=0}^{k-1}w^{-ij}(sw^j+1)^L$
组合数有个非常优美的性质
$$ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}$$
推波式子
$$
\begin{aligned}
F(i)&=\sum_{j=0}^{k-1}w^{-ij}(sw^j+1)^L\\
&=\sum_{j=0}^{k-1}w^{\binom{i}{2}+\binom{j}{2}-\binom{i+j}{2}}(sw^j+1)^L\\
&=w^\binom{i}{2}\sum_{j=0}^{k-1}w^{-\binom{i+j}{2}}w^\binom{j}{2}(sw^j+1)^L\\
\end{aligned}
$$
发现这是一个差卷积的形式
扔一个$ MTT$上去就能拿那$ 40$分了
考虑$ n>1$的情况
相当于把$ s$从一个数变成了矩阵，把$ 1$变成单位矩阵
$ (sw^j+1)^L$这个矩阵我们只需要关注一个位置上的值
因此可以乘出来然后取[x][y]这个位置上的数即可
这样就可以通过此题
复杂度：大常数单 $\log$
代码
#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
ll x=0;char zf=1;char ch=getchar();
while(ch!='-'&&!isdigit(ch))ch=getchar();
if(ch=='-')zf=-1,ch=getchar();
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();return x*zf;
}
void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
void writeln(const ll y){write(y);putchar('\n');}
int k,m,n,x,y,cnt,p,L;
int w[65539],val[65539];
int f[65539],g[65539],F[65539];
int ksm(int x,int y=p-2){
int ans=1;
for(;y;y>>=1,x=1ll*x*x%p)if(y&1)ans=1ll*ans*x%p;
return ans;
}
struct mat{
ll a[3][3];
void print(){
for(rt i=0;i<n;i++)for(rt j=0;j<n;j++)cout<<a[i][j]<<" \n"[j==n-1];
}
inline mat operator*(const mat s)const{
mat ret={0};
for(rt i=0;i<n;i++)
for(rt k=0;k<n;k++)
for(rt j=0;j<n;j++)
(ret.a[i][j]+=a[i][k]*s.a[k][j]);
for(rt i=0;i<n;i++)for(rt j=0;j<n;j++)ret.a[i][j]%=p;
return ret;
}
mat operator*(const int s)const{
mat ret;
for(rt i=0;i<n;i++)
for(rt j=0;j<n;j++)
ret.a[i][j]=a[i][j]*s%p;
return ret;
}
mat operator+(const mat s)const{
mat ret;memset(ret.a,0,sizeof(ret.a));
for(rt i=0;i<n;i++)
for(rt j=0;j<n;j++)
ret.a[i][j]=(a[i][j]+s.a[i][j])%p;
return ret;
}
}I,zy;
mat ksm(mat x,int y){
mat ans=I;
for(;y;y>>=1,x=x*x)if(y&1)ans=ans*x;
return ans;
}
int V(int x){
return 1ll*x*(x-1)/2%k;
}
const double PI=acos(-1.0);
struct cp{
double x,y;
cp operator +(const cp &s)const{return {x+s.x,y+s.y};}
cp operator -(const cp &s)const{return {x-s.x,y-s.y};}
cp operator *(const cp &s)const{return {x*s.x-y*s.y,x*s.y+y*s.x};}
}W[18][1<<17],A[270010],B[270010],C[270010],D[270010];
int R[270010];
void FFT(const int n,cp *A){
for(rt i=0;i<n;i++)if(i>R[i])swap(A[i],A[R[i]]);
for(rt j=0;j<n;j+=2){
const cp x=A[j],y=A[j+1];
A[j]=x+y,A[j+1]=x-y;
}
for(rt i=2,s=1;i<n;i<<=1,s++)
for(rt j=0;j<n;j+=i<<1)
for(rt k=0;k<i;k+=2){
cp x=A[j+k],y=W[s][k]*A[i+j+k];
A[j+k]=x+y;A[i+j+k]=x-y;
x=A[j+k+1],y=W[s][k+1]*A[i+j+k+1];
A[j+k+1]=x+y;A[i+j+k+1]=x-y;
}
}
int yg(int n){
for(rt i=1;i<=n;i++){
for(rt j=2;j*j<=n;j++)if((n-1)%j==0)
if(ksm(i,(n-1)/j)==1)goto end;
return i;end:;
}
}
void subtask(){
w[0]=1;w[1]=ksm(yg(p),(p-1)/k);
for(rt i=2;i<k;i++)w[i]=1ll*w[i-1]*w[1]%p;mat s=zy;
for(rt i=0;i<k;i++){
g[i]=ksm(s*w[i]+I,L).a[x-1][y-1]*w[V(i)]%p;
}
for(rt i=0;i<k+k;i++)f[i]=w[(k-V(i))%k];
for(rt i=0;i<k/2;i++)swap(g[i],g[k-i]);
n=k+k-1;m=k;
int lim=1;while(lim<=n+m)lim<<=1;
for(rt i=0;(1<<i)<lim;i++){
W[i][0]={1,0};
W[i][1]={cos(PI/(1<<i)),sin(PI/(1<<i))};
for(rt j=2;j<1<<i;j++)
if(j%32==0)W[i][j]={cos(PI*j/(1<<i)),sin(PI*j/(1<<i))};
else W[i][j]=W[i][j-1]*W[i][1];
}
for(rt i=0;i<=n;i++){
A[i].x=f[i]>>15;
A[i].y=f[i]&32767;
}
for(rt i=0;i<=m;i++){
B[i].x=g[i]>>15;
B[i].y=g[i]&32767;
}
for(rt i=1;i<lim;i++)R[i]=(R[i>>1]>>1)|(i&1)*(lim>>1);
FFT(lim,A);FFT(lim,B);
for(rt i=0;i<lim;i++){
const int pl=(lim-1)&(lim-i);
const cp ca={A[pl].x,-A[pl].y},cb={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};
}
FFT(lim,C);FFT(lim,D);
for(rt i=k;i<k+k;i++){
const int vv=1ll*w[V(i-k)]*ksm(k,p-2)%p;
ll a=C[i].x/lim+0.5,b=C[i].y/lim+0.5,c=D[i].x/lim+0.5,d=D[i].y/lim+0.5;
a=((a%p)<<30)+(((b+c)%p)<<15)+d;a=(a%p+p)%p;
writeln(1ll*a*vv%p);
}
exit(0);
}
int jc[100010],njc[100010],inv[100010],ff[100010][3],ans[100010];
int c(int x,int y){
return 1ll*jc[x]*njc[y]%p*njc[x-y]%p;
}
void bl(){
for(rt i=0;i<2;i++)jc[i]=njc[i]=inv[i]=1;
for(rt i=2;i<=L;i++){
jc[i]=1ll*jc[i-1]*i%p;
inv[i]=1ll*inv[p%i]*(p-p/i)%p;
njc[i]=1ll*njc[i-1]*inv[i]%p;
}
ff[0][x-1]=1;if(x==y)ans[0]=1;
for(rt i=1;i<=L;i++)
for(rt j=0;j<n;j++){
for(rt k=0;k<n;k++)
(ff[i][j]+=1ll*zy.a[k][j]*ff[i-1][k]%p)%=p;
if(j==y-1)(ans[i%k]+=1ll*ff[i][j]*c(L,i)%p)%=p;
}
for(rt i=0;i<k;i++)writeln(ans[i]);exit(0);
}
int main(){
cin>>n>>k>>L>>x>>y>>p;
for(rt i=0;i<n;i++)
for(rt j=0;j<n;j++)
cin>>zy.a[i][j];
if(L<=100000)bl();
for(rt i=0;i<n;i++)I.a[i][i]=1;
subtask();
return 0;
}
标签：ch,const,int,HNOI,2019,return,binom,include,之舞 来源： https://www.cnblogs.com/DreamlessDreams/p/10672461.html