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

md为了做这题恶补了之前一知半解的很多多项式相关姿势，才总算是写掉了这题

首先考虑一个暴力的矩乘做法，我们考虑设一个矩阵$f$，其中$f_{i,j}$即为$i$轮后状态为$i$的人的分数，很显然我们可以找出一个矩阵$B$来进行转移，那么最后要求的就是$B^t\times f_0$

然后我们进行一个转化，令$\oplus$为三进制下不进位加法，$\ominus$为三进制下不退位减法，它们两个互为逆运算

然后我们根据剪刀石头布的定义可以知道当$k\in [0,3^m)$时，有$B_{i,j}=B_{i\oplus k,j\oplus k}=B_{i\ominus k,j\ominus k}$

因此我们有$B_{i,j}=B_{i\ominus j,0}$，然后它也可以用归纳法扩展到$B^n$

所以$f_t=B^t\times f_0$的第$i$项就是：

$$f_{n,i}=\sum_{k} B_{k,i}^t \times f_{0,k}=\sum_{k} B_{0,i\ominus k}^t \times f_{0,k}=\sum_{x\oplus y=i} B_{0,y}^t\times f_{0,x}$$

然后我们发现这个式子就是一个卷积的形式，而且只需要$B$的第一行即可

那么我们再化一下式子：

$$B_{0,i}^t =\sum_{k} B_{0,k}^{t-1} \times B_{k,i}=\sum_{k} B_{0,k}^{t-1}\times B_{0,i\ominus k}=\sum_{x\oplus y=i} B_{0,x}^{t-1}\times B_{0,y}$$

那么我们发现其实只要把$B$的第一行和自己做$t$次卷积即可，然后再卷上$f_0$就是答案

那么现在我们要找到一个三进制不进位加法的卷积变换，那么它其实是一个循环卷积，它的一般形式是：

$$A\times B=\sum_{i}\sum_{(j+k) \mod n=i} A_j\times B_k$$

对于循环卷积的话其实就是对于它的卷积的矩阵$T$的每一行$x$都有：

$$x_k\times x_t=x_{k\oplus t}$$

因此$T$的每一行都是方程组$x_k\times x_t=x_{k\oplus t}$的一组解

但是我们都知道如果卷积没有逆变换（就叫鸡卷好了）就GG了，因此我们的$T$还必须有逆矩阵

然后这也就说明$T$的行列式不为$0$，因此方程组$x_k\times x_t=x_{k\oplus t}$至少要有$n$组不同的解

所以这个矩阵$T$应该怎么构造呢？然后如果你认真学了FFT的姿势就会知道这个矩阵就是范德蒙德矩阵

因此矩阵长这样：

$$T=\left[\begin{matrix} \omega_n^0&\omega_n^0&\cdots&\omega_n^0\\ \omega_n^0&\omega_n^1&\cdots&\omega_n^{n-1}\\ \omega_n^0&\omega_n^2&\cdots&\omega_n^{2(n-1)}\\ \vdots &\vdots&\ddots&\vdots\\ \omega_n^0&\omega_n^{n-1}&\cdots&\omega_n^{(n-1)(n-1)}\\ \end{matrix}\right]$$

即$T_{i,j}=\omega_n^{ij}$

那么它的逆矩阵为：

$$T^{-1}=\frac{1}{n}\left[\begin{matrix} \omega_n^0&\omega_n^0&\cdots&\omega_n^0\\ \omega_n^0&\omega_n^{-1}&\cdots&\omega_n^{-(n-1)}\\ \omega_n^0&\omega_n^{-2}&\cdots&\omega_n^{-2(n-1)}\\ \vdots &\vdots&\ddots&\vdots\\ \omega_n^0&\omega_n^{-(n-1)}&\cdots&\omega_n^{-(n-1)(n-1)}\\ \end{matrix}\right]$$

然后因为这题是三进制不进位加法，因此我们选三次单位根$\omega$，然后我们发现所有的复数都可以被表示成$a+b\omega$的形式

又由单位根的一个性质$\omega^2+\omega+1=0$，因此有$\omega^2=-\omega-1$，因此我们可以将乘出来的结果降幂，因此我们发现现在复数的乘法就是：

$$(a+b\omega)\times (c+d\omega)=ac+(ad+bc)\omega+bd(-\omega-1)=(ac-bd)+(ad+bc-bd)\omega$$

然后我们由上面的可知：

$$T=\left[\begin{matrix} 1&1&1\\ 1&\omega&\omega^2\\ 1&\omega^2&\omega \end{matrix}\right]$$

同时我们知道：

$$\left[\begin{matrix} 1&1&1\\ 1&\omega&\omega^2\\ 1&\omega^2&\omega \end{matrix}\right] \times \left[\begin{matrix} 1&1&1\\ 1&\omega^2&\omega\\ 1&\omega&\omega^2 \end{matrix}\right]=3I$$

其中$I$是单位矩阵，因此我们就把后面那个当作逆矩阵来算就好了

然后由于DFT和IDFT的本质都是分治做向量对矩阵的乘法，它每次IDFT的时候都多乘了一个$3$，而一共有$\log_3 n$层，因此总共多乘了$3^{\log_3 n}=n$次，那么最后只要把所有元素都除以一个$n$就好了！

PS：是不是明白了一般的FFT的转移系数是哪里来的，以及做一般卷积的时候要把$n$扩大一倍（这样做循环卷积不会被膜到前面去），以及IDFT之后为什么要除以$n$这些问题

然后提到除以$n$我们就会想到是不是有逆元的问题，因为我们要求$n=3^m$有逆元，那么也就是$3\not | \ p$

假设$3|p$，那么我们令$k=\frac{p}{3}$，那么有：

$$\frac{1}{k+1}+\frac{1}{k(k+1)}=\frac{1}{k}=\frac{p}{3}$$

与假设矛盾！因此得证$3\not |\ p$

所以这题终于是做完了……

#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=531441,M=20;
int m,t,mod,n,x,b[M+5][M+5];
inline int sum(CI x,CI y)
{
	int t=x+y; return t>=mod?t-mod:t;
}
inline int sub(CI x,CI y)
{
	int t=x-y; return t<0?t+mod:t;
}
inline void exgcd(int& x,int& y,CI m,CI n)
{
	if (!n) return (void)(x=1,y=0);
	exgcd(y,x,n,m%n); y-=m/n*x;
}
inline int inv(CI x)
{
	int m,n; exgcd(m,n,x,mod);
	return (m%mod+mod)%mod;
}
struct Complex
{
	int x,y; //x+y*w
	inline Complex(CI X=0,CI Y=0)
	{
		x=X; y=Y;
	}
	friend inline Complex operator + (const Complex& A,const Complex& B)
	{
		return Complex(sum(A.x,B.x),sum(A.y,B.y));
	}
	friend inline Complex operator - (const Complex& A,const Complex& B)
	{
		return Complex(sub(A.x,B.x),sub(A.y,B.y));
	}
	friend inline Complex operator * (const Complex& A,const Complex& B)
	{
		return Complex(sub(1LL*A.x*B.x%mod,1LL*A.y*B.y%mod),
		sub(sum(1LL*A.x*B.y%mod,1LL*A.y*B.x%mod),1LL*A.y*B.y%mod));
	}
}f[N],g[N];
inline Complex quick_pow(Complex x,int p,Complex mul=Complex(1,0))
{
	for (;p;p>>=1,x=x*x) if (p&1) mul=mul*x;
	return mul;
}
inline Complex calc1(const Complex& A) //calc A*w
{
	return Complex(mod-A.y,sub(A.x,A.y));
}
inline Complex calc2(const Complex& A) //calc A*w^2
{
	return Complex(sub(A.y,A.x),mod-A.x);
}
inline void DFT(Complex *f)
{
	RI mid,j,k; for (mid=1;mid<n;mid*=3)
	for (j=0;j<n;j+=mid*3) for (k=0;k<mid;++k)
	{
		Complex x=f[j+k],y=f[j+k+mid],z=f[j+k+(mid<<1)];
		f[j+k]=x+y+z; f[j+k+mid]=x+calc1(y)+calc2(z);
		f[j+k+(mid<<1)]=x+calc2(y)+calc1(z);
	}
}
inline void IDFT(Complex *f)
{
	RI mid,j,k; for (mid=1;mid<n;mid*=3)
	for (j=0;j<n;j+=mid*3) for (k=0;k<mid;++k)
	{
		Complex x=f[j+k],y=f[j+k+mid],z=f[j+k+(mid<<1)];
		f[j+k]=x+y+z; f[j+k+mid]=x+calc2(y)+calc1(z);
		f[j+k+(mid<<1)]=x+calc1(y)+calc2(z);
	}
	int invn=inv(n); for (RI i=0;i<n;++i) f[i].x=1LL*f[i].x*invn%mod;
}
int main()
{
	RI i,j; for (scanf("%d%d%d",&m,&t,&mod),n=i=1;i<=m;++i) n*=3;
	for (i=0;i<n;++i) scanf("%d",&x),f[i].x=x;
	for (i=0;i<=m;++i) for (j=0;j<=m-i;++j) scanf("%d",&b[i][j]);
	for (i=0;i<n;++i)
	{
		int x=i,cw=0,cl=0; while (x)
		{
			int nw=x%3; if (nw==1) ++cw;
			else if (nw==2) ++cl; x/=3;
		}
		g[i].x=b[cw][cl];
	}
	for (DFT(f),DFT(g),i=0;i<n;++i) f[i]=f[i]*quick_pow(g[i],t);
	for (IDFT(f),i=0;i<n;++i) printf("%d\n",f[i].x);
	return 0;
}