题意：有一张顶点数为 (L+1)×n(L+1)\times n 的有向图。这张图的每个顶点由一个二元组 (u,v)(u,v) 表示 (0≤u≤L,1≤v≤n)(0\le u\le L,1\le v\le n)。这张图不是简单图，对于任意两个顶点 (u1,v1)(u2,v2)(u_1,v_1)(u_2,v_2)，如果 u1<u2u_1< u_2，则从 (u1,v1)(u_1,v_1) 到 (u2,v2)(u_2,v_2) 一共有 wv1,v2w _{v_1,v_2} 条不同的边，如果 u1≥u2u_1\ge u_2 则没有边。白兔将在这张图上上演一支舞曲。白兔初始时位于该有向图的顶点 (0,x)(0,x)。白兔将会跳若干步。每一步，白兔会从当前顶点沿任意一条出边跳到下一个顶点。白兔可以在任意时候停止跳舞（也可以没有跳就直接结束）。当到达第一维为 LL 的顶点就不得不停止，因为该顶点没有出边。假设白兔停止时，跳了 mm 步，白兔会把这只舞曲给记录下来成为一个序列。序列的第 ii 个元素为它第 ii 步经过的边。问题来了：给定正整数 kk 和 yy（1≤y≤n1\le y\le n），对于每个 tt（0≤t<k0\le t< k），求有多少种舞曲（假设其长度为 mm）满足 m mod k=tm \bmod k=t，且白兔最后停在了坐标第二维为 yy 的顶点？两支舞曲不同定义为它们的长度（mm）不同或者存在某一步它们所走的边不同。输出的结果对 pp 取模。

数据范围与约定：

对于全部数据，pp 为一个质数，108<p<23010^{8}<p<2^{30}，1≤n≤31\le n\le 3，1≤x≤n1\le x\le n，1≤y≤n1\le y\le n，0≤wi,j<p0\le w_{i,j}<p，1≤k≤655361\le k\le 65536，kk 为 p−1p-1 的约数，1≤L≤1081\le L\le 10^8。

对于每组测试点，特殊限制如下：

  • 测试点 1,21,2：L≤105L\le 10^5；
  • 测试点 33：n=1,w(1,1)=1n=1,w(1,1)=1，kk 的最大质因子为 22；
  • 测试点 44：n=1n=1，kk 的最大质因子为 22；
  • 测试点 55：n=1,w(1,1)=1n=1,w(1,1)=1；
  • 测试点 66：n=1n=1；
  • 测试点 7,87,8：kk 的最大质因子为 22。

解析：一道很难的多项式题。

我们先考虑 n=1n=1 的情况，这时候从任意一个点到另外一个点都有 a=w1,1a=w _{1,1} 条边，因此我们可以得到答案为：

anst=∑i mod k=tLai(Li)ans _{t}=\sum _{i\bmod k=t} ^{L}a ^{i}\dbinom{L}{i}

其中 ii 枚举的是白兔跳了多少步，每跳一步有 aa 种方案，根据乘法原理方案为 aia ^{i}，由于跳了 ii 步，故除起点外有 ii 个经过的点，这 ii 个点不确定，选择的方案有 (Li)\dbinom{L}{i} 种，故总方案数为 ai(Li)a ^{i}\dbinom{L}{i}。

我们转化一下，可得：

anst=∑i=0L[k∣(i−t)]ai(Li)\begin{aligned} ans _{t}=\sum _{i=0} ^{L}\left[k|(i-t)\right]a^{i}\dbinom{L}{i} \end{aligned}

其中出现了 [k∣(i−t)][k|(i-t)] 这个式子，接下来就是关键的一步。

单位根反演

一个比较冷门的算法，它的主要形式为：

[n∣k]=1n∑i=0n−1ωnki[n|k]=\dfrac{1}{n}\sum _{i=0} ^{n-1}\omega _{n} ^{ki}

证明分为两种情况：

  • 若 n∣kn\mid k，我们设 k=npk=np，则：

  • [n∣k]=1n∑i=0n−1ωnnpi=1n∑i=0n−1ωn0=1n×n=1[n\mid k]=\dfrac{1}{n}\sum _{i=0} ^{n-1}\omega _{n} ^{npi}=\dfrac{1}{n}\sum _{i=0} ^{n-1}\omega _{n} ^{0}=\dfrac{1}{n}\times n=1

  • 若 n∤kn\nmid k，根据等比数列求和公式，则：

  • [n∤k]=1n∑i=0n−1(ωnk)i=1n(ωn0(1−ωnkn)1−ωnk)=0[n\nmid k]=\dfrac{1}{n}\sum _{i=0} ^{n-1}(\omega _{n} ^{k}) ^{i}=\dfrac{1}{n}\left(\dfrac{\omega _{n} ^{0}(1-\omega _{n} ^{kn})}{1-\omega _{n} ^{k}}\right)=0

它的主要应用是对于一个数列 aia _{i}，求：

∑i=0n[k∣i]ai\sum _{i=0} ^{n}[k|i]a _{i}

我们首先构造 aia _{i} 的生成函数 f(x)=∑i=0naixi\displaystyle{f(x)=\sum _{i=0} ^{n}a _{i}x ^{i}}，那么：

∑i=0n[k∣i]ai=1k∑i=0k−1f(ωki)\sum _{i=0} ^{n}[k|i]a _{i}=\dfrac{1}{k}\sum _{i=0} ^{k-1}f(\omega _{k}^{i})

证明如下：

1k∑i=0k−1f(ωki)=1k∑i=0k−1∑j=0naj(ωki)j=∑j=0naj(1k∑i=0k−1ωkji)=∑j=0naj[k∣j]=∑i=0n[k∣i]ai\begin{aligned} \dfrac{1}{k}\sum _{i=0} ^{k-1}f(\omega _{k} ^{i}) & =\dfrac{1}{k}\sum _{i=0} ^{k-1}\sum _{j=0} ^{n}a _{j}(\omega _{k} ^{i}) ^{j}\\ & =\sum _{j=0} ^{n}a _{j}\left(\dfrac{1}{k}\sum _{i=0} ^{k-1}\omega _{k} ^{ji}\right)\\ & =\sum _{j=0} ^{n}a _{j}[k|j]\\ & =\sum _{i=0} ^{n}[k|i]a _{i} \end{aligned}

同理，若我们要求：

∑i=0n[i mod k=p]ai\sum _{i=0} ^{n}[i\bmod k=p]a _{i}

我们把 f(x)f(x) 平移 pp 即可，即：

f(x)=apx0+ap+1x1+⋯f(x)=a _{p}x ^{0}+a _{p+1}x ^{1}+\cdots

再用上面的式子即可求得所有 i mod k=pi\bmod k=p 的 aia _{i} 和。

了解单位根反演后，代入原式，得：

anst=∑i=0L1k∑j=0k−1(ωk(i−t)j)ai(Li)=1k∑j=0k−1ωk−tj∑i=0Lωkijai(Li)\begin{aligned} ans _{t} & =\sum _{i=0} ^{L}\dfrac{1}{k}\sum _{j=0} ^{k-1}(\omega _{k} ^{(i-t)j})a ^{i}\dbinom{L}{i}\\ & =\dfrac{1}{k}\sum _{j=0} ^{k-1}\omega _{k} ^{-tj}\sum _{i=0} ^{L}\omega _{k} ^{ij}a ^{i}\dbinom{L}{i}\\ \end{aligned}

后面这个和式 ∑i=0Lωkijai(Li)\displaystyle{\sum _{i=0} ^{L}\omega _{k} ^{ij}a ^{i}\dbinom{L}{i}} 可以继续用二项式定理化简，即：

∑i=0Lωkijai(Li)=∑i=0L(Li)(ωkja)i1L−i=(ωkja+1)L\begin{aligned} \sum _{i=0} ^{L}\omega _{k} ^{ij}a ^{i}\dbinom{L}{i} & =\sum _{i=0} ^{L}\dbinom{L}{i}(\omega _{k} ^{j}a) ^{i}1 ^{L-i}\\ & =(\omega _{k} ^{j}a+1) ^{L} \end{aligned}

发现 (ωkja+1)L(\omega _{k} ^{j}a+1) ^{L} 只与 jj 有关，设为 pjp _{j}，代回原式：

anst=1k∑j=0k−1ωk−tjpj\begin{aligned} ans _{t}=\dfrac{1}{k}\sum _{j=0} ^{k-1}\omega _{k} ^{-tj}p _{j} \end{aligned}

接下来是一个神仙化简操作，利用如下恒等式：

−tj=(t2)+(j2)−(t+j2)-tj=\dbinom{t}{2}+\dbinom{j}{2}-\dbinom{t+j}{2}

代回原式：

anst=1k∑j=0k−1ωk(t2)+(j2)−(t+j2)pj=1kωk(t2)∑j=0k−1ωk(j2)pjωk−(t+j2)\begin{aligned} ans _{t} & =\dfrac{1}{k}\sum_{j=0} ^{k-1}\omega _{k} ^{\binom{t}{2}+\binom{j}{2}-\binom{t+j}{2}}p _{j}\\ & =\dfrac{1}{k}\omega _{k} ^{\binom{t}{2}}\sum _{j=0} ^{k-1}\omega _{k} ^{\binom{j}{2}}p _{j}\omega _{k} ^{-\binom{t+j}{2}} \end{aligned}

其实现在就依稀可以看出后面那个和式是一个卷积形式，由于模数不一定要用 MTT 处理。为了更清楚地判断后面那个式子是不是卷积，我们继续化简。设 fj=ωk(j2)pj,gj=ωk−(j2)f _{j}=\omega _{k} ^{\binom{j}{2}}p _{j},g _{j}=\omega _{k} ^{-\binom{j}{2}}，则：

anst=1kωk(t2)∑j=0k−1fjgt+j\begin{aligned} ans _{t} & =\dfrac{1}{k}\omega _{k} ^{\binom{t}{2}}\sum _{j=0} ^{k-1}f _{j}g _{t+j} \end{aligned}

我们把枚举 jj 转为枚举 t+jt+j，则：

anst=1kωk(t2)∑j=tk−1+tfj−tgjans _{t}=\dfrac{1}{k}\omega _{k} ^{\binom{t}{2}}\sum _{j=t} ^{k-1+t}f _{j-t}g _{j}

然后经典套路，设 f(k−1)−i′=fif' _{(k-1)-i}=f _{i}，则：

anst=1kωk(t2)∑j=tk−1+tf(k−1+t)−j′gjans _{t}=\dfrac{1}{k}\omega _{k} ^{\binom{t}{2}}\sum _{j=t} ^{k-1+t}f' _{(k-1+t)-j}g _{j}

这时候能很轻易地看出后面的和式即为卷积形式，由于 tt 是变化的，且 t∈[0,k−1]t\in [0,k-1]，所以 f′f' 的下标范围恒为 [0,k−1][0,k-1]，gg 的下标范围为 [0,2k−2][0,2k-2]。故乘积后长度为 3k−33k-3。直接 MTT 处理多项式乘法即可。注意单位根 ωki\omega _{k} ^{i} 需要预处理，因此需要先求 pp 的原根。

n=1n=1 的形式处理完后，n=3n=3 的形式就不难了，实际上就是把 n=1n=1 的 aa 换成输入的矩阵，同时改一下 pjp _{j} 的求法，先算出 st×(ωkja+e)Lst\times (\omega _{k} ^{j}a+e) ^{L}，其中 stst 是初始状态矩阵，根据题意只有 (1,x)(1,x) 处值为 11，ee 则为单位矩阵，这里就要用到矩阵乘法。根据题意可得结束状态为 (1,y)(1,y)，故 pjp _{j} 为上面这个式子算出的矩阵 (1,y)(1,y) 处的值。

代码中将矩阵的横纵坐标都减去了 11，即从 00 开始。注意 MTT 要记得开 long double，包括常量 π\pi。

/*
 * @Author: clorf 
 * @Date: 2021-01-31 09:58:26 
 * @Last Modified by:   clorf 
 * @Last Modified time: 2021-01-31 09:58:26 
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 2e9
using namespace std;
const int maxn=140000;
const long double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
//6.常量数字与long long连用时要加ll！！！
//7.考虑无解和多解的情况
//8.两个int连成爆long long的一定要加1ll！！！！
//9.先除后乘！！！！！！
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)){f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,k,L,x,y,p;
int maxx=1,num,omega[maxn>>1];
inline int inc(int a,int b){return a+b>=p?a+b-p:a+b;}
namespace poly
{
    int rev[maxn*3];
    struct Cnum
    {
        long double x,y;
        Cnum(long double p=0,long double q=0){x=p,y=q;}
        Cnum conj(){return Cnum(x,-y);}
        friend Cnum operator + (Cnum a,Cnum b){return Cnum(a.x+b.x,a.y+b.y);}
        friend Cnum operator - (Cnum a,Cnum b){return Cnum(a.x-b.x,a.y-b.y);}
        friend Cnum operator * (Cnum a,Cnum b){return Cnum(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
    }a[maxn*3],b[maxn*3],c[maxn*3],d[maxn*3];
    inline void ready(int len)
    {
        while(maxx<=len)
        {
            maxx<<=1;
            num++;
        }
        for(int i=0;i<maxx;i++)
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(num-1));
    }
    inline void FFT(Cnum *a,int maxx,int type)
    {
        for(int i=0;i<maxx;i++)
            if(i<rev[i])
                swap(a[i],a[rev[i]]);
        for(int mid=1;mid<maxx;mid<<=1)
        {
            Cnum w1(cos(Pi/mid),type*sin(Pi/mid));
            for(int r=mid<<1,j=0;j<maxx;j+=r)
            {
                Cnum w(1,0);
                for(int k=0;k<mid;k++,w=w*w1)
                {
                    Cnum x=a[j|k];
                    Cnum y=a[j|k|mid]*w;
                    a[j|k]=x+y;
                    a[j|k|mid]=x-y;
                }
            }
        }
        if(!(~type))
            for(int i=0;i<maxx;i++)
            {
                a[i].x/=maxx;
                a[i].y/=maxx;
            }
    }
    inline void MTT(int *f,int *g,int *h,int maxx,int p)
    {
        for(int i=0;i<maxx;i++)
        {
            a[i].x=f[i]>>15;
            a[i].y=f[i]&32767;
            c[i].x=g[i]>>15;
            c[i].y=g[i]&32767;
        }
        FFT(a,maxx,1);
        FFT(c,maxx,1);
        for(int i=1;i<maxx;i++)
        {
            b[maxx-i]=a[i].conj();
            d[maxx-i]=c[i].conj();
        }
        b[0]=a[0].conj();
        d[0]=c[0].conj();
        for(int i=0;i<maxx;i++)
        {
            Cnum A=(a[i]+b[i])*Cnum(0.5,0);
            Cnum B=(a[i]-b[i])*Cnum(0,-0.5);
            Cnum C=(c[i]+d[i])*Cnum(0.5,0);
            Cnum D=(c[i]-d[i])*Cnum(0,-0.5);
            a[i]=A*C+Cnum(0,1)*(A*D+B*C);
            b[i]=B*D;
        }
        FFT(a,maxx,-1);
        FFT(b,maxx,-1);
        for(int i=0;i<maxx;i++)
        {
            int A=(long long)(a[i].x+0.5)%p;
            int B=(long long)(a[i].y+0.5)%p;
            int C=(long long)(b[i].x+0.5)%p;
            h[i]=(1ll*A*(1<<30)%p+1ll*B*(1<<15)%p+C+p)%p;
        }
    }
}
int f[maxn*3],g[maxn*3],h[maxn*3];
struct matrix
{
    int s[3][3];
    matrix(){memset(s,0,sizeof(s));}
    friend matrix operator + (matrix a,matrix b)
    {
        matrix c;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                c.s[i][j]=inc(a.s[i][j],b.s[i][j]);
        return c;
    }
    friend matrix operator * (matrix a,int b)
    {
        matrix c;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                c.s[i][j]=1ll*a.s[i][j]*b%p;
        return c;
    }
    friend matrix operator * (matrix a,matrix b)
    {
        matrix c;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                for(int k=0;k<3;k++)
                    c.s[i][j]=inc(c.s[i][j],1ll*a.s[i][k]*b.s[k][j]%p);
        return c;
    }
}e,st,w;
inline int power(int a,int b,int p)
{
    int ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=1ll*ans*a%p;
        a=1ll*a*a%p;
    }
    return ans;
}
inline matrix power(matrix a,int b)
{
    matrix ans=e;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a;
        a=a*a;
    }
    return ans;
}
inline int getroot(int x)
{
    static int fac[50],cnt=0;
    int num=x-1;
    for(int i=2;i<=x-1;i++)
    {
        if(num==1)
            break;
        if(!(num%i))
        {
            fac[++cnt]=i;
            while(!(num%i))
                num/=i;
        }
    }
    for(int g=2;;g++)
    {
        bool flag=1;
        for(int j=1;j<=cnt;j++)
            if(power(g,(x-1)/fac[j],x)==1)
            {
                flag=0;
                break;
            }
        if(flag)
            return g;
    }
}
int main()
{
    scanf("%d%d%d%d%d%d",&n,&k,&L,&x,&y,&p);
    x--,y--;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            scanf("%d",&w.s[i][j]);
    e.s[0][0]=e.s[1][1]=e.s[2][2]=1;
    st.s[0][x]=1;
    omega[0]=1;
    int G=getroot(p);
    omega[1]=power(G,(p-1)/k,p);
    for(int i=2;i<k;i++)
        omega[i]=1ll*omega[i-1]*omega[1]%p;
    int len=(k<<1)-2;
    for(int i=0;i<=len;i++)
        g[i]=omega[(k-1ll*i*(i-1)/2%k)%k];
    for(int i=0;i<k;i++)
        f[k-1-i]=1ll*omega[1ll*i*(i-1)/2%k]*(st*power(w*omega[i]+e,L)).s[0][y]%p;
    len+=(k-1);
    poly::ready(len);
    poly::MTT(f,g,h,maxx,p);
    int invk=power(k,p-2,p);
    for(int i=0;i<k;i++)
        printf("%lld\n",1ll*h[k-1+i]*invk%p*omega[1ll*i*(i-1)/2%k]%p);
    return 0;
}