[loj3653]抽奖机

用$n$位三进制数描述状态，记$F_{S}=\sum_{i=1}^{m}[a_{i}=x,b_{i}=y]$（其中$x$和$y$分别为$S$在三进制下1和2的个数）

定义$x\oplus y$为$x$和$y$在三进制下执行不进位加法的结果$,\odot $为$\oplus$对应的卷积，问题即求$F\odot F...\odot F$（共$k$个）

关于$\odot$，每一位可以理解为对每一位做一个模3的加法卷积，以三次单位根$\omega$做NTT即可

具体的，构造变换$f(A)_{S}=\sum_{T}\omega^{\sum_{i=0}^{n-1}trans(S_{i},T_{i})}A_{T}$（其中$S_{i}$和$T_{i}$分别指两者三进制下第$i$位）

关于$trans$，即NTT的系数中$\omega$的幂次，具体形式为$\left(\begin{array}{ll}0&0&0\\0&1&2\\0&2&1\end{array}\right)$

通过变换$f$，有$f(A\odot B)_{S}=f(A)_{S}\times f(B)_{S}$和$f^{-1}(A)_{S}=\frac{f(A)_{S\oplus S}}{3^{n}}$，具体均可代入证明

注意到$F_{S}$仅取决于$S$在三进制下1和2的个数，记$F_{a,b}$表示含有$a$个1和$b$个2的$F_{S}$的值，以此即可表示$F$

将此时$F_{a,b}$对$f(F)_{p,q}$的贡献系数构建二元生成函数$H_{a,b}(x,y)$​，看作一个二维背包，即有
$$
H_{a,b}(x,y)=\frac{{n\choose a,b}}{n\choose x,y}(1+\omega x+\omega^{2}y)^{a}(1+\omega^{2}x+\omega y)^{b}(1+x+y)^{n-a-b}
$$
提出常数项后，记其为$H'_{a,b}$，考虑其递推式，即
$$
(1+x+y)H'_{a,b}(x,y)=(1+\omega^{2}x+\omega y)H'_{a,b-1}(x,y)
$$
（对于$b=0$的情况，可以类似地对$a$递推，注意处理边界）

由于内存限制，上述递推需要滚动实现，那么在递推过程中实现变换$f$即可

最终，总复杂度为$o(n^{4}+n^{2}\log k)$，可以通过

#include<bits/stdc++.h>
using namespace std;
#define N 125
#define mod 1000000009
#define w 115381398
#define w2 884618610
#define ll long long
int n,m,x,y,C[N][N],h[N][N],g[N][N],f[N][N],F0[N][N],F[N][N];
ll k;
int qpow(int n,ll m){
    int s=n,ans=1;
    while (m){
        if (m&1)ans=(ll)ans*s%mod;
        s=(ll)s*s%mod,m>>=1;
    }
    return ans;
}
void turn_a(){
    memset(h,0,sizeof(h));
    for(int i=0;i<=n;i++)
        for(int j=0;j<=n-i;j++){
            h[i][j]=g[i][j];
            if (i)h[i][j]=(h[i][j]+(ll)w*g[i-1][j]-h[i-1][j]+mod)%mod;
            if (j)h[i][j]=(h[i][j]+(ll)w2*g[i][j-1]-h[i][j-1]+mod)%mod;
        }
    memcpy(g,h,sizeof(g));
}
void turn_b(){
    memset(h,0,sizeof(h));
    for(int i=0;i<=n;i++)
        for(int j=0;j<=n-i;j++){
            h[i][j]=f[i][j];
            if (i)h[i][j]=(h[i][j]+(ll)w2*f[i-1][j]-h[i-1][j]+mod)%mod;
            if (j)h[i][j]=(h[i][j]+(ll)w*f[i][j-1]-h[i][j-1]+mod)%mod;
        }
    memcpy(f,h,sizeof(f));
}
void FWT(int p=0){
    memset(g,0,sizeof(g));
    memset(F0,0,sizeof(F0));
    for(int i=0;i<=n;i++)
        for(int j=0;j<=n-i;j++){
            g[i][j]=(ll)C[n][i]*C[n-i][j]%mod;
            F[i][j]=(ll)F[i][j]*g[i][j]%mod;
        }
    for(int i=0;i<=n;i++){
        memcpy(f,g,sizeof(f));
        for(int j=0;j<=n-i;j++){
            for(int x=0;x<=n;x++)
                for(int y=0;y<=n-x;y++)F0[x][y]=(F0[x][y]+(ll)F[i][j]*f[x][y])%mod;
            turn_b();
        }
        turn_a();
    }
    for(int i=0;i<=n;i++)
        for(int j=0;j<=n-i;j++){
            int s=(ll)C[n][i]*C[n-i][j]%mod;
            F[i][j]=(ll)F0[i][j]*qpow(s,mod-2)%mod;
        }
    if (p){
        int s=qpow(qpow(3,n),mod-2);
        for(int i=0;i<=n;i++)
            for(int j=0;(j<i)&&(j<=n-i);j++)swap(F[i][j],F[j][i]);
        for(int i=0;i<=n;i++)
            for(int j=0;j<=n-i;j++)F[i][j]=(ll)F[i][j]*s%mod;
    }
}
int main(){
    for(int i=0;i<N;i++){
        C[i][0]=C[i][i]=1;
        for(int j=1;j<i;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
    scanf("%d%d%lld",&n,&m,&k);
    for(int i=1;i<=m;i++){
        scanf("%d%d",&x,&y);
        F[x][y]++;
    }
    FWT();
    for(int i=0;i<=n;i++)
        for(int j=0;j<=n-i;j++)F[i][j]=qpow(F[i][j],k);
    FWT(1);
    for(int i=0;i<=n;i++)
        for(int j=0;j<=n-i;j++)F[i][j]=(ll)F[i][j]*C[n][i]%mod*C[n-i][j]%mod;
    for(int i=0;i<=n;i++){
        for(int j=0;j<=n-i;j++)printf("%d ",F[i][j]);
        printf("\n");
    }
    return 0;
}