【清华集训2016】石家庄的工人阶级队伍比较坚强

【清华集训】石家庄的工人阶级队伍比较坚强

给定 \(m\)，然后令 \(n=3^m\)，给定数组 \(f_{0,0},f_{0,1}...f_{0,3^m-1}\)

给定数组 \(b\)，然后定义卷积：

\[f_{k,x}=\sum_{y}f_{k-1,y}b_{u,v} \]

其中 \(u\) 定义为 \(W(x,y)\)，\(v\) 定义为 \(L(x,y)\)

其中 \(W(x,y)\) 被定义为将三进制的每一位进行比较，出现一组 10/21/02 就累加一次贡献，\(L(x,y)\) 定义为 \(W(y,x)\)

你需要计算 \(f_{t}\)，\(t\) 给定。

答案对 \(p\) 取模，保证不存在 \(x,y\) 满足 \(\frac{1}{x}+\frac{1}{y}=\frac{3}{p}\)

\(m\le 12,t\le 10^9,p\le 10^9\)

Solution

又来水题解了。

考虑三进制不退位减法：

\[\begin{matrix} 0\oplus 0=0&0\oplus 1=2&0\oplus 2=1 \\1\oplus 1=0&1\oplus 2=2&1\oplus 0=1 \\2\oplus 2=0&2\oplus 0=2&2\oplus 1=1 \end{matrix}\]

我们发现结果为 \(1\) 等价于 win，结果为 \(2\) 等价于 lose

设 \(bit_k(x)\) 表示 \(x\) 中的 \(k\) 的数量。

那么我们有转移形如：

\[f_{k,x}=\sum_y f_{k-1,y}b_{bit_1(x\oplus y),bit_2(x\oplus y)} \]

不难发现贡献实际上只和 \(x\oplus y\) 相关。

设 \(\odot\) 表示三进制不进位加法，那么不难发现：

\[f_{k,x}=\sum_{y\odot z=x} f_{k-1,y}b_{z} \]

于是相当于计算 3 进制下异或卷积的 \(k\) 次幂。

我们搞出 FWT 的点值即可处理了。

注意到 \(k\) 进制 FWT 相当于在做循环卷积，这样我们只需要将给 \(k\) 进制单位根代入进去即可。

具体式子：

\[FWT(i)=\sum_{j=0} c(i,j)A_j \]

\[FWT(i)=\sum_{j=0}^{k-1} c(i_{\mathbf{M}},j)\sum c(i',j')A_{j'+j\times k^{\mathbf{M}}} \]

\[FWT(i)=\sum_{j=0}^{k-1} \omega_{k}^{i_{\mathbf{M}}j}\sum c(i',j')A_{j'+j\times k^{\mathbf{M}}} \]

这样递归即可，由于只需要乘以单项式，所以复杂度是 \(\mathcal O(k^{m+3})\) 的。

由于要次幂一下，会掉精度，所以只能扩域，我们找一个代数单位来表示 \(\omega_3\)

根据单位根的性质，我们注意到 \(\omega_3^2+\omega_3+1=0\)，这样的话就有 \(\omega_3^2=-\omega_3-1\)

这样我们先相当于拿着 \(ax^2+bx+c\) 类型的多项式在模(?)上做运算，最后我们再把 \(a\) 表示成 \(b,c\) 相关消一消。

• 坑：这个模数保证了 \(3\) 存在逆元，但是不能直接套质数情况的求法...

然后发现这个题卡这个做法，然后你会发现没必要 FWT 的时候取模，所以可以最后取模。

然后手动循环展开一下，就可以过了。最后放一下比较丑的代码

\(Code:\)

#include<bits/stdc++.h>
using namespace std ;
#define Next( i, x ) for( register int i = head[x]; i; i = e[i].next )
#define rep( i, s, t ) for( register int i = (s); i <= (t); ++ i )
#define drep( i, s, t ) for( register int i = (t); i >= (s); -- i )
#define re register
#define int long long
int gi() {
	char cc = getchar() ; int cn = 0, flus = 1 ;
	while( cc < '0' || cc > '9' ) {  if( cc == '-' ) flus = - flus ; cc = getchar() ; }
	while( cc >= '0' && cc <= '9' )  cn = cn * 10 + cc - '0', cc = getchar() ;
	return cn * flus ;
}
const int N = 1e6 + 5 ; 
const int M = 20 + 5 ; 
int lim, m, t, P, b[M][M] ; 
struct gf {
	int a[3] ; 
	void init() { a[0] = a[1] = a[2] = 0 ; }
} f[N], g[N], Nx[N] ;
int fpow(int x, int k) {
	int ans = 1, base = x ;
	while(k) {
		if(k & 1) ans = 1ll * ans * base % P ;
		base = 1ll * base * base % P, k >>= 1 ;
	} return ans ;
}
void inc(int &x, int y) { (x += y) ; }
int bit1(int x) { int ans = 0 ; while( x ) ans += (x % 3 == 1), x /= 3 ; return ans ; }
int bit2(int x) { int ans = 0 ; while( x ) ans += (x % 3 == 2), x /= 3 ; return ans ; }
gf Move(gf x, int l) {
	gf ans ; ans.init() ; 
	ans.a[(0 + l) % 3] = x.a[0], ans.a[(1 + l) % 3] = x.a[1], ans.a[(2 + l) % 3] = x.a[2] ;
	return ans ; 
}
gf operator * (gf x, gf y) {
	gf ans ; ans.init() ; 
	ans.a[0] = ( x.a[0] * y.a[0] + x.a[1] * y.a[2] + x.a[2] * y.a[1] ) % P ; 
	ans.a[1] = ( x.a[0] * y.a[1] + x.a[1] * y.a[0] + x.a[2] * y.a[2] ) % P ; 
	ans.a[2] = ( x.a[0] * y.a[2] + x.a[1] * y.a[1] + x.a[2] * y.a[0] ) % P ; 
	return ans ; 
}
gf operator * (gf x, int y) { gf ans ; rep(i, 0, 2) ans.a[i] = x.a[i] * y % P ; return ans ; }
gf operator + (gf x, gf y) { gf ans ; rep(i, 0, 2) ans.a[i] = (x.a[i] + y.a[i]) ; return ans ; }
gf fpow(gf x, int k) {
	gf ans, base = x ; ans.init(), ans.a[0] = 1 ; 
	while(k) {
		if(k & 1) ans = ans * base ;
		base = base * base, k >>= 1 ; 
	} return ans ; 
}
gf operator % (gf x, int p) { rep( i, 0, 2 ) x.a[i] %= p ; return x ; }
void FWT(gf *a) {
	for(re int k = 1; k < lim; k *= 3 ) for(re int i = 0; i < lim; i += (k * 3) )
	for(re int j = i; j < i + k; ++ j) {
		rep( l, 0, 2 ) Nx[l].init() ; 
		rep( l, 0, 2 ) rep( o, 0, 2 ) 
		Nx[l] = (Nx[l] + Move(a[j + o * k], o * l % 3)) ; 
		rep( l, 0, 2 ) a[j + l * k] = Nx[l] ; 
	}
}
void IFWT(gf *a) {
	for(re int k = 1; k < lim; k *= 3 ) for(re int i = 0; i < lim; i += (k * 3) )
	for(re int j = i; j < i + k; ++ j) {
		rep( l, 0, 2 ) Nx[l].init() ; 
		rep( l, 0, 2 ) rep( o, 0, 2 ) 
		Nx[l] = (Nx[l] + Move(a[j + o * k], 3 - o * l % 3)) ; 
		rep( l, 0, 2 ) a[j + l * k] = Nx[l] ; 
	}
}
int phi(int x) {
	int ans = x ; 
	for(re int i = 2; i <= sqrt(x); ++ i) {
		if(x % i) continue ;
		ans = (ans - ans / i) ;
		while(x % i == 0) x /= i ; 
	} if( x != 1 ) ans = ans - ans / x ;
	return ans ;  
}
signed main()
{
	m = gi(), t = gi(), P = gi(), lim = 1 ; 
	rep( i, 1, m ) lim = lim * 3 ; 
	rep( i, 0, lim - 1 ) f[i].a[0] = gi() ; 
	rep( i, 0, m ) rep( j, i, m ) b[i][j - i] = gi() ; 
	rep( i, 0, lim - 1 ) g[i].a[0] = b[bit1(i)][bit2(i)] ; 
	FWT(f), FWT(g) ;
	rep( i, 0, lim - 1 ) f[i] = f[i] % P, g[i] = g[i] % P ; 
	rep( i, 0, lim - 1 ) f[i] = f[i] * fpow(g[i], t) ;
	int iv = fpow(lim, phi(P) - 1) ; IFWT(f) ;
	rep( i, 0, lim - 1 ) f[i] = f[i] % P ; 
	rep( i, 0, lim - 1 ) f[i] = f[i] * iv ;
	rep( i, 0, lim - 1 ) {
		int d = (f[i].a[0] - f[i].a[2] + P) % P ; 
		printf("%lld\n", d ) ; 
	}
	return 0 ;
}