假设$x,yin {0,1,2}$，则$x$能赢$y$（根据题中定义）当且仅当$x-yequiv 1(mod 3)$
定义$ominus$为两数3进制下不退位的减法，$S_{x}$表示$x$在3进制下1的个数，则$u=W(x,y)=S_{xominus y}$，类似地，有$v=W(y,x)=S_{yominus x}$
记$B_{x,y}=b_{S_{xominus y},S_{yominus x}}$，令$k=xominus y$，根据$ominus$的性质，有$kominus 0=xominus y$以及$0ominus k=yominus x$，代入后就可以得到$B_{x,y}=b_{S_{kominus 0},S_{0ominus k}}=B_{k,0}$，再记$g_{k}=B_{k,0}$，可以预处理出来
原题也就是$f_{i,x}=sum_{y=0}^{3^{m}-1}g_{xominus y}f_{i-1,y}$，再令$oplus$为3进制下按位加法，有$(xominus y)oplus y=x$，也就是一个$oplus$的卷积，记这个卷积为$otimes$运算，有$f_{i}=gotimes f_{i-1}$
根据卷积的结合律，有$f_{t}=g^{t}otimes f_{0}$（其中$g^{t}$也是上面的$otimes$），也就是要快速计算$otimes$运算
我们直接来考虑原问题的一个拓展，即令$oplus$为$B$进制下不进位加法，$otimes$为$oplus$卷积（$B$是一个较小的数），如何对两个序列快速计算$otimes$这个运算
为了方便起见，不妨假设序列长度$n=B^{k}$，不足可以补0
对于两个序列$a_{i}$和$b_{i}$，考虑构造一个$n imes n$的矩阵$A_{i,j}$（将序列看作一个$1 imes n$的矩阵），令$a'=aA$，$b'=bA$，再令$c'=a'*b'$（这里的*指对应位置相乘），最后$c=c'A^{-1}$
$forall n=B^{k}(kge 1)$，若矩阵$A$满足以下两个条件，则最终有$c=aotimes b$：
1.其左上角$B imes B$的矩阵在$n=B$时满足$c=aotimes b$
2.$A_{i,j}=A_{lfloorfrac{i}{B} floor,lfloorfrac{j}{B} floor}A_{i mod B,j mod B}$，且$A^{-1}$具有同样性质
证明如下——
注意到当$n=B$时，左上角$B imes B$的矩阵即为$A$本身，那么由条件1可知其正确
归纳$n=B^{k}(kge 1)$时正确，考虑$n=B^{k+1}$时正确性：
对于每一个$0le i<B^{k+1}$，都可以被表示为$i_{0}B+i_{1}$（其中$0le i_{0}<B^{k}$，$0le i_{1}<B$）
对于$a_{i}$和$b_{i}$，代入式子即有
$$
c_{i}=sum_{j=0}^{B^{k+1}-1}c'_{j}A^{-1}_{j,i}=sum_{j_{1}=0}^{B-1}sum_{j_{0}=0}^{B^{k}-1}c'_{j_{0}B+j_{1}}A^{-1}_{j_{0}B+j_{1},i_{0}B+i_{1}}=sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}sum_{j_{0}=0}^{B^{k}-1}c'_{j_{0}B+j_{1}}A^{-1}_{j_{0},i_{0}}
$$
考虑$c'_{i}$，同样可以代入得到
$$
c'_{i}=a'_{i}b'_{i}=(sum_{j_{1}=0}^{B-1}A_{j_{1},i_{1}}sum_{j_{0}=0}^{B^{k}-1}a_{j_{0}B+j_{1}}A_{j_{0},i_{0}})(sum_{j_{1}=0}^{B-1}A_{j_{1},i_{1}}sum_{j_{0}=0}^{B^{k}-1}b_{j_{0}B+j_{1}}A_{j_{0},i_{0}})
$$
（关于$a'_{i}$和$b'_{i}$用$a_{i}$和$b_{i}$表达，与$c_{i}$用$c'_{i}$表达相同，将$A^{-1}$改为$A$即可）
对$a_{i}$和$b_{i}$分别构造出$B$个序列，即$d_{j_{0},i}=a_{iB+j_{0}}$以及$e_{j_{0},i}=b_{iB+j_{0}}$，考虑$d_{x}$和$e_{y}$分别作为$a_{i}$和$b_{i}$通过前面的方式计算，可以算出$c'$，记作$f'_{x,y}$，即有
$$
f'_{x,y,i}=(sum_{j=0}^{B^{k}-1}d_{x,j}A_{j,i})(sum_{j=0}^{B^{k}-1}e_{y,j}A_{j,i})=(sum_{j=0}^{B^{k}-1}a_{jB+x}A_{j,i})(sum_{j=0}^{B^{k}-1}b_{jB+y}A_{j,i})
$$
可以发现这个式子与上面的$c'_{i}$有关，具体来说即
$$
c'_{i}=sum_{x=0}^{B-1}A_{x,i_{1}}sum_{y=0}^{B-1}A_{y,i_{1}}f'_{x,y,i_{0}}
$$
整体将$c'_{i}$的式子代入$c_{i}$中，也就得到了
$$
c_{i}=sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}sum_{j_{0}=0}^{B^{k}-1}A^{-1}_{j_{0},i_{0}}sum_{x=0}^{B-1}A_{x,j_{1}}sum_{y=0}^{B-1}A_{y,j_{1}}f'_{x,y,j_{0}}=sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}sum_{x=0}^{B-1}A_{x,j_{1}}sum_{y=0}^{B-1}A_{y,j_{1}}sum_{j_{0}=0}^{B^{k}-1}A^{-1}_{j_{0},i_{0}}f'_{x,y,j_{0}}
$$
显然最后一个式子也就是$f'_{x,y}A^{-1}$，记作$f_{x,y}$，根据归纳有$f_{x,y}=d_{x}otimes e_{y}$，联系其定义即
$$
f_{x,y,i}=sum_{0le p_{0},q_{0}<B^{k},p_{0}oplus q_{0}=i}d_{x,p_{0}}e_{y,q_{0}}
$$
再次代入$c_{i}$中，即
$$
c_{i}=sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}sum_{x=0}^{B-1}A_{x,j_{1}}sum_{y=0}^{B-1}A_{y,j_{1}}sum_{0le p_{0},q_{0}<B^{k},p_{0}oplus q_{0}=i_{0}}d_{x,p_{0}}e_{y,q_{0}}=sum_{0le p_{0},q_{0}<B^{k},p_{0}oplus q_{0}=i_{0}}sum_{j_{1}=0}^{B-1}A^{-1}_{j_{1},i_{1}}sum_{x=0}^{B-1}A_{x,j_{1}}d_{x,p_{0}}sum_{y=0}^{B-1}A_{y,j_{1}}e_{y,q_{0}}
$$
对于$p_{0}$和$q_{0}$，构造$u_{i}=d_{i,p_{0}}$以及$v_{i}=e_{i,q_{0}}$，将两者作为$a_{i}$和$b_{i}$进行计算，关于$x$的枚举也就是算出$u'_{j_{1}}$，第二个关于$y$的枚举也就是$v'_{j_{1}}$，两者对应位置相乘，也就是$c'_{j_{1}}$，再乘上前面的$A_{j_{1},i_{1}}^{-1}$，也就是$c_{i_{1}}$
（这里的$c'_{j_{1}}$和$c_{i_{1}}$指代入$u_{i}$和$v_{i}$后用前面的方式计算得到，与之前的$c'_{i}$与$c_{i}$无关）
进一步的，根据归纳有$c_{i_{1}}=(uotimes v)_{i_{1}}=sum_{0le p_{1},q_{1}<B,p_{1}oplus q_{1}=i_{1}}u_{p_{1}}v_{q_{1}}$，代入即
$$
c_{i}=sum_{0le p_{0},q_{0}<B^{k},p_{0}oplus q_{0}=i_{0}}sum_{0le p_{1},q_{1}<B,p_{1}oplus q_{1}=i_{1}}a_{p_{0}B+p_{1}}b_{q_{0}B+q_{1}}=sum_{0le p,q<B^{k+1},poplus q=i}a_{p}b_{q}=(aotimes b)_{i}
$$
（最后就利用到了$oplus$是按位运算的）
确定矩阵$A$的条件后，我们还要具体的构造出一个矩阵$A$，通过矩阵左上角$B imes B$的部分就可以根据第二个条件唯一确定$A$了，也就是如何构造$n=B$时的矩阵
不难发现这就是一个幂次对$B$取模的加法卷积，联系FFT，记$omega=cosfrac{2pi}{B}+sinfrac{2pi}{B}i$，取$A_{i,j}=omega^{ij}$，逆矩阵将则令$A^{-1}_{i,j}=frac{omega^{-(ij)}}{B}$，下面来证明一下：
当我们确定了$A_{i,j}=omega^{i}$，相当于将$a_{i}$看作一个$B$次多项式，求出其在$omega^{i}$的值
根据余数定理，当我们知道一个多项式在$x=omega^{i}$上的值，事实上也就是该多项式模$x-omega^{i}$后的结果，而求出了所有$omega^{i}$，那么也就知道了原多项式对$prod_{i=0}^{B-1}(x-omega^{i})$取模的结果
再将点值相乘， 也就确定了最终多项式对$prod_{i=0}^{B-1}(x-omega^{i})$取模后的结果
注意到$omega$满足$omega^{B}=1$且$forall 1le i<B,omega^{i} e 1$，则$prod_{i=0}^{B-1}(x-omega^{i})=x^{B}-1$，而$x^{n}equiv x^{n mod B}(mod x^{B}-1)$，因此恰好也就求出了下标对$B$取模的卷积结果
关于$prod_{i=0}^{B-1}(x-omega^{i})=x^{B}-1$，由于其可以看作$B$对点值${x=omega^{i},y=0}$，且每一个$x$各不相同，由此恰好可以确定一个$B$次多项式，且$x^{B}-1$恰好满足此性质，因此即相等
然而，本题需要对$p$取模，如果用实数表示的误差无法解决，代入$B=3$，可得$omega=frac{sqrt{3}}{2}i-frac{1}{2}$
考虑题中关于$p$的限制：不存在正整数$x,y$满足$frac{1}{x}+frac{1}{y}=frac{3}{p}$
若$pequiv 0(mod 3)$，取$x=y=frac{2p}{3}$即满足条件，因此$frac{1}{3}$在模$p$意义下有意义，但并不能保证3在模$p$意义下有二次剩余，换言之$sqrt{3}$不一定有意义（$frac{1}{2}$有意义也可以证，但不关心）
为了使系数为整数，考虑将$omega$作为一个整体，显然每一个虚数与$a+bomega$一一对应，同时根据$omega^{3}=1$以及$omega^{2}=-omega-1$（这两个条件也就是根据$prod_{i=0}^{B-1}(x-omega^{i})=x^{B}-1$得到），当产生高次项也都可以降为一次来表示
最终答案必然是整数，因此$b=0$，即为$a$
另外，计算矩阵乘法可以利用条件2来做，以由$a$计算$a'$为例子，即
$$
a'_{i}=sum_{j=0}^{B^{k}-1}a_{j}A_{j,i}=sum_{j_{1}=0}^{B-1}A_{j_{1},i_{1}}sum_{j_{0}=0}^{B^{k-1}-1}a_{j_{0}B+j_{1}}A_{j_{0},i_{0}}
$$
对于每一个$j_{1}$，后面也就是一个$B^{k-1}$规模的子问题（矩阵乘法），递归后将其求出后，再代入该式即可，这样单次矩乘时间复杂度为$o(kB^{k+1})$
由于要快速幂，总复杂度为$o((m+log_{2}t)3^{m})$可以通过
补充一下与本题无关的一些事情：
$B=2$时这个问题就是FWT，且做法与其完全相同
当$B$足够大时即为FFT，但时间复杂度会退化为$o(n^{2})$
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 600005 4 #define M 13 5 int n,m,t,mod,S[N],base[N][M],b[M][M]; 6 struct Complex{ 7 int a,b; 8 Complex(){ 9 a=b=0; 10 } 11 Complex(int x){ 12 a=x,b=0; 13 } 14 Complex(int x,int y){ 15 a=x,b=y; 16 } 17 Complex operator + (const Complex &k)const{ 18 return Complex((a+k.a)%mod,(b+k.b)%mod); 19 } 20 Complex operator * (const Complex &k)const{ 21 int s=mod-1LL*b*k.b%mod; 22 return Complex((1LL*a*k.a+s)%mod,(1LL*a*k.b+1LL*b*k.a+s)%mod); 23 } 24 }inv3,A[3][3],invA[3][3],a[N],ans[N]; 25 Complex pow(Complex n,int m){ 26 Complex s=n,ans=Complex(1); 27 while (m){ 28 if (m&1)ans=ans*s; 29 s=s*s; 30 m>>=1; 31 } 32 return ans; 33 } 34 void DFT(Complex *a){ 35 Complex aa[3]; 36 for(int i=0,s=1;i<m;i++,s*=3) 37 for(int j=0;j<n;j++) 38 if (!base[j][i]){ 39 for(int k=0;k<3;k++)aa[k]=Complex(); 40 for(int k=0;k<3;k++) 41 for(int l=0;l<3;l++)aa[k]=aa[k]+a[j+l*s]*A[l][k]; 42 for(int k=0;k<3;k++)a[j+k*s]=aa[k]; 43 } 44 } 45 void IDFT(Complex *a){ 46 Complex aa[3]; 47 for(int i=0,s=1;i<m;i++,s*=3) 48 for(int j=0;j<n;j++) 49 if (!base[j][i]){ 50 for(int k=0;k<3;k++)aa[k]=Complex(); 51 for(int k=0;k<3;k++) 52 for(int l=0;l<3;l++)aa[k]=aa[k]+a[j+l*s]*invA[l][k]; 53 for(int k=0;k<3;k++)a[j+k*s]=aa[k]; 54 } 55 } 56 int main(){ 57 scanf("%d%d%d",&m,&t,&mod); 58 n=1; 59 for(int i=1;i<=m;i++)n*=3; 60 for(int i=0;i<n;i++){ 61 S[i]=S[i/3]+(i%3==1); 62 base[i][0]=i%3; 63 for(int j=1;j<m;j++)base[i][j]=base[i/3][j-1]; 64 } 65 int k=mod,phi=mod; 66 for(int i=2;i*i<=k;i++) 67 if (k%i==0){ 68 phi=phi/i*(i-1); 69 while (k%i==0)k/=i; 70 } 71 if (k>1)phi=phi/k*(k-1); 72 inv3=pow(Complex(3),phi-1); 73 for(int i=0;i<3;i++) 74 for(int j=0;j<3;j++){ 75 A[i][j]=pow(Complex(0,1),i*j); 76 invA[i][j]=pow(Complex(0,1),6-i*j)*inv3; 77 } 78 for(int i=0;i<n;i++){ 79 scanf("%d",&k); 80 a[i]=Complex(k); 81 } 82 for(int i=0;i<=m;i++) 83 for(int j=0;j<=m-i;j++)scanf("%d",&b[i][j]); 84 for(int i=0;i<n;i++){ 85 k=0; 86 for(int j=m-1;j>=0;j--) 87 if (!base[i][j])k=k*3; 88 else k=k*3+3-base[i][j]; 89 ans[i]=Complex(b[S[i]][S[k]]); 90 } 91 DFT(a); 92 DFT(ans); 93 for(int i=0;i<n;i++)ans[i]=a[i]*pow(ans[i],t); 94 IDFT(ans); 95 for(int i=0;i<n;i++)printf("%d ",ans[i].a); 96 }