[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)$,可以通过
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 125 4 #define mod 1000000009 5 #define w 115381398 6 #define w2 884618610 7 #define ll long long 8 int n,m,x,y,C[N][N],h[N][N],g[N][N],f[N][N],F0[N][N],F[N][N]; 9 ll k; 10 int qpow(int n,ll m){ 11 int s=n,ans=1; 12 while (m){ 13 if (m&1)ans=(ll)ans*s%mod; 14 s=(ll)s*s%mod,m>>=1; 15 } 16 return ans; 17 } 18 void turn_a(){ 19 memset(h,0,sizeof(h)); 20 for(int i=0;i<=n;i++) 21 for(int j=0;j<=n-i;j++){ 22 h[i][j]=g[i][j]; 23 if (i)h[i][j]=(h[i][j]+(ll)w*g[i-1][j]-h[i-1][j]+mod)%mod; 24 if (j)h[i][j]=(h[i][j]+(ll)w2*g[i][j-1]-h[i][j-1]+mod)%mod; 25 } 26 memcpy(g,h,sizeof(g)); 27 } 28 void turn_b(){ 29 memset(h,0,sizeof(h)); 30 for(int i=0;i<=n;i++) 31 for(int j=0;j<=n-i;j++){ 32 h[i][j]=f[i][j]; 33 if (i)h[i][j]=(h[i][j]+(ll)w2*f[i-1][j]-h[i-1][j]+mod)%mod; 34 if (j)h[i][j]=(h[i][j]+(ll)w*f[i][j-1]-h[i][j-1]+mod)%mod; 35 } 36 memcpy(f,h,sizeof(f)); 37 } 38 void FWT(int p=0){ 39 memset(g,0,sizeof(g)); 40 memset(F0,0,sizeof(F0)); 41 for(int i=0;i<=n;i++) 42 for(int j=0;j<=n-i;j++){ 43 g[i][j]=(ll)C[n][i]*C[n-i][j]%mod; 44 F[i][j]=(ll)F[i][j]*g[i][j]%mod; 45 } 46 for(int i=0;i<=n;i++){ 47 memcpy(f,g,sizeof(f)); 48 for(int j=0;j<=n-i;j++){ 49 for(int x=0;x<=n;x++) 50 for(int y=0;y<=n-x;y++)F0[x][y]=(F0[x][y]+(ll)F[i][j]*f[x][y])%mod; 51 turn_b(); 52 } 53 turn_a(); 54 } 55 for(int i=0;i<=n;i++) 56 for(int j=0;j<=n-i;j++){ 57 int s=(ll)C[n][i]*C[n-i][j]%mod; 58 F[i][j]=(ll)F0[i][j]*qpow(s,mod-2)%mod; 59 } 60 if (p){ 61 int s=qpow(qpow(3,n),mod-2); 62 for(int i=0;i<=n;i++) 63 for(int j=0;(j<i)&&(j<=n-i);j++)swap(F[i][j],F[j][i]); 64 for(int i=0;i<=n;i++) 65 for(int j=0;j<=n-i;j++)F[i][j]=(ll)F[i][j]*s%mod; 66 } 67 } 68 int main(){ 69 for(int i=0;i<N;i++){ 70 C[i][0]=C[i][i]=1; 71 for(int j=1;j<i;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod; 72 } 73 scanf("%d%d%lld",&n,&m,&k); 74 for(int i=1;i<=m;i++){ 75 scanf("%d%d",&x,&y); 76 F[x][y]++; 77 } 78 FWT(); 79 for(int i=0;i<=n;i++) 80 for(int j=0;j<=n-i;j++)F[i][j]=qpow(F[i][j],k); 81 FWT(1); 82 for(int i=0;i<=n;i++) 83 for(int j=0;j<=n-i;j++)F[i][j]=(ll)F[i][j]*C[n][i]%mod*C[n-i][j]%mod; 84 for(int i=0;i<=n;i++){ 85 for(int j=0;j<=n-i;j++)printf("%d ",F[i][j]); 86 printf("\n"); 87 } 88 return 0; 89 }