题目描述
略
Sol
考场上暴力 (O(L)) 50分真良心。
简单的推一下式子，对于一个 t 来说，答案就是：
[sum_{i=0}^{L} [k|(i-t)] {Lchoose i}F(i)
]
就是对于所有 mod k 的结果是 t 的 i 的后面那一坨东西的和。
(F(i)) 表示走了 (i) 步从纵坐标 x 到纵坐标为 y 的方案数，这个东西显然非常好递推。
所以 (O(L)) 的暴力做法就直接递推组合数就行了。
然后是正解。
([k|(i-t)])这玩意可以套用单位根的性质: ([k|n]=frac{1}{k}sum_{i=0}^kw_k^{in})
正确性结合单位根性质和等比数列求和就可以得到。
所以式子就变成了:
[frac{1}{k}sum_{i=0}^{L} {Lchoose i}F(i)sum_{j=0}^k w_k^{j(i-t)}
]
显而易见的把里面拆了然后交换求和顺序:
[frac{1}{k}sum_{j=0}^k w_k^{-jt}sum_{i=0}^{L}w_k^{ji} {Lchoose i}F(i)
]
后面的东西是个套路，具体见 bzoj3328 PYXFIB
我们把 (F(i)) 写乘一个矩阵幂次的形式,因为是可以矩阵乘法递推的。
假设转移矩阵是 (T)
[frac{1}{k}sum_{j=0}^k w_k^{-jt}sum_{i=0}^{L}w_k^{ji} {Lchoose i}T^i
]
真正的 (F(i)) 就是 (T^i[x][y])
后面就是一个组合数+幂次项，是一个的二项式定理的形式，合并一下就是：
[frac{1}{k}sum_{j=0}^k w_k^{-jt}(Tw_k^j+I)^L
]
(I) 是单位矩阵。
后面的东西对于一个 (j) 而言显然是固定的， 我们设 (G(w_k^j)=(Tw_k^j+I)^L)
所以要求的是:
[frac{1}{k}sum_{j=0}^k (w_k^{-t})^j G(w_k^j)
]
这玩意看上去好眼熟啊。不就是个 (IDFT) 吗?
我们 (IDFT) 就是对于求出点值后的多项式代入单位根的倒数,最后结果乘上 (frac{1}{n})。
所以就是这么回事了。就是让我们对求完后的 (G(x)) 给 (IDFT) 一下回去就是所有 (t) 的答案了
所以问题变成求解任意长度的 (DFT) ，这里我们就没有办法用 性能优化 那题的方法 (DFT) 了
然后就可以使用 (Bluestein’s; Algorithm) 了。(具体可参见我的性能优化的题解)
注意矩阵里面的运算要卡常，先用 long long 存下答案再取模，不然常数大不开 O2 过不了。(出这种毒瘤题也就算了怎么还卡常啊QAQ)
code:
#include<bits/stdc++.h>
using namespace std;
#define Set(a,b) memset(a,b,sizeof(a))
template<class T>inline void init(T&x){
x=0;char ch=getchar();bool t=0;
for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
if(t) x=-x;return;
}typedef long long ll;
int mod;const int MAXN=65536;
template<class T>inline void Inc(T&x,int y){x+=y;if(x>=mod) x-=mod;}
template<class T>inline void Dec(T&x,int y){x-=y;if(x < 0 ) x+=mod;}
template<class T>inline int fpow(int x,T k){int ret=1;for(;k;k>>=1,x=(ll)x*x%mod)if(k&1) ret=(ll)ret*x%mod;return ret;}
inline int Sum(int x,int y){x+=y;if(x>=mod) x-=mod;return x;}
inline int Dif(int x,int y){x-=y;if(x < 0 ) x+=mod;return x;}
int n,k,L,x,y,w[3][3];
int F[MAXN],g,wn[MAXN],iwn[MAXN];
inline void Getroot(){
int phi=mod-1;int x=phi;
static int pri[50],cnt=0;
for(int i=2;i*i<=x;++i) {
if(x%i==0) {pri[++cnt]=i,x/=i;while(x%i==0) x/=i;}
}
for(g=2;;++g){bool fl=1;
for(int i=1;i<=cnt;++i) if(fpow(g,phi/pri[i])==1) {fl=0;break;}
if(fl)return;
}
}
namespace Work1{
struct matrix{
ll a[3][3];matrix(){Set(a,0);}
inline ll* operator [](int x){return a[x];}
inline matrix operator +(matrix b){matrix c;for(int i=0;i<n;++i) for(int j=0;j<n;++j) c[i][j]=Sum(a[i][j],b[i][j]);return c;}
inline matrix operator *(matrix b){
matrix c;
for(int i=0;i<n;++i)
for(int j=0;j<n;++j){
for(int k=0;k<n;++k) c[i][j]+=a[i][k]*b[k][j];
c[i][j]%=mod;
}
return c;
}
}T,I;
inline matrix matrixfpow(matrix x,int k){matrix ret=I;for(;k;k>>=1,x=x*x)if(k&1)ret=ret*x;return ret;}
inline void work(){
for(int i=0;i<n;++i) I[i][i]=1;
for(int i=0;i<k;++i) {T=matrix();
for(int u=0;u<n;++u)for(int v=0;v<n;++v) T[u][v]=(ll)w[u][v]*wn[i]%mod;
T=T+I;T=matrixfpow(T,L);F[i]=T[x][y];
}
return;
}
}
namespace Work2{
const int MAXM=MAXN<<2;
typedef double db;int rader[MAXM];
inline int Init(int n){
int len=1,up=-1;while(len<=n) len<<=1,++up;
for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
return len;
}
namespace MTT{
const db PI=acos(-1);
struct Complex{
db x,y;Complex(db _x=0.0,db _y=0.0){x=_x,y=_y;}
inline Complex operator +(const Complex B){return Complex(x+B.x,y+B.y);}
inline Complex operator -(const Complex B){return Complex(x-B.x,y-B.y);}
inline Complex operator *(const Complex B){return Complex(x*B.x-y*B.y,x*B.y+y*B.x);}
}w[MAXM];
inline void Calc(int n){for(int i=1;i<n;i<<=1) for(int j=0;j<i;++j) w[n/i*j]=Complex(cos(PI/i*j),sin(PI/i*j));return;}
inline void FFT(Complex*A,int n,int f){
for(int i=0;i<n;++i) if(rader[i]>i) swap(A[rader[i]],A[i]);
for(int i=1;i<n;i<<=1)
for(int j=0,p=i<<1;j<n;j+=p)
for(int k=0;k<i;++k){
Complex X=A[j|k],Y=A[j|k|i]* ((~f)? w[n/i*k]:Complex(w[n/i*k].x,-w[n/i*k].y));
A[j|k]=X+Y,A[j|k|i]=X-Y;
}
if(!~f) for(int i=0;i<n;++i) A[i].x/=(db)n;return;
}
inline void Mul(int*A,int*B,int*C,int len){
static Complex A1[MAXM],A2[MAXM],B1[MAXM],B2[MAXM];
int MO=sqrt(mod);
for(int i=0;i<len;++i) A1[i]=Complex(A[i]/MO,0.0),B1[i]=Complex(A[i]%MO,0.0),A2[i]=Complex(B[i]/MO,0.0),B2[i]=Complex(B[i]%MO,0.0);
FFT(A1,len,1),FFT(A2,len,1),FFT(B1,len,1),FFT(B2,len,1);
for(int i=0;i<len;++i) {Complex X;
X=A1[i]*A2[i],A2[i]=A2[i]*B1[i];
B1[i]=B1[i]*B2[i];B2[i]=B2[i]*A1[i];
A1[i]=X,A2[i]=A2[i]+B2[i];
}int MOD=MO*MO%mod;
FFT(A1,len,-1),FFT(B1,len,-1),FFT(A2,len,-1);
for(int i=0;i<len;++i) {
int X=(ll)(A1[i].x+0.5)%mod,Y=(ll)(B1[i].x+0.5)%mod,Z=(ll)(A2[i].x+0.5)%mod;
int ans=(ll)MOD*X%mod;Inc(ans,(ll)MO*Z%mod);Inc(ans,Y);
C[i]=ans;
}return;
}
}using MTT::Calc;
inline int Co(int x){return (ll)x*(x-1)/2%k;}
inline void DFT(int*A,int n,int len){
int m=2*n-1;static int F[MAXM],G[MAXM];
for(int i=0;i<n;++i) F[i]=(ll)A[i]*wn[Co(i)]%mod;for(int i=n;i<len;++i) F[i]=0;
for(int i=0;i<m;++i) G[i]=iwn[Co(i)];for(int i=m;i<len;++i) G[i]=0;
reverse(F,F+n);MTT::Mul(F,G,F,len);
for(int k=0,i=n-1;i<m;++i,++k) A[k]=(ll)F[i]*wn[Co(k)]%mod;
for(int i=0,inv=fpow(n,mod-2);i<n;++i) A[i]=(ll)A[i]*inv%mod;
return;
}
void work(){
int len=Init(3*k-3);Calc(len);DFT(F,k,len);
for(int i=0;i<k;++i) printf("%d
",F[i]);
return;
}
}
int main()
{
freopen("dance.in","r",stdin);
freopen("dance.out","w",stdout);
init(n),init(k),init(L),init(x),init(y),init(mod);--x,--y;Getroot();
wn[0]=iwn[0]=1,wn[1]=fpow(g,(mod-1)/k),iwn[1]=fpow(wn[1],mod-2);
for(int i=2;i<k;++i) wn[i]=(ll)wn[i-1]*wn[1]%mod,iwn[i]=(ll)iwn[i-1]*iwn[1]%mod;
for(int i=0;i<n;++i) for(int j=0;j<n;++j) init(w[i][j]);
Work1::work();Work2::work();
return 0;
}