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

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

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

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

然后我们根据剪刀石头布的定义可以知道当(kin [0,3^m))时，有(B_{i,j}=B_{ioplus k,joplus k}=B_{iominus k,jominus k})

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

所以(f_t=B^t imes f_0)的第(i)项就是：

[f_{n,i}=sum_{k} B_{k,i}^t imes f_{0,k}=sum_{k} B_{0,iominus k}^t imes f_{0,k}=sum_{xoplus y=i} B_{0,y}^t imes f_{0,x} ]

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

那么我们再化一下式子：

[B_{0,i}^t =sum_{k} B_{0,k}^{t-1} imes B_{k,i}=sum_{k} B_{0,k}^{t-1} imes B_{0,iominus k}=sum_{xoplus y=i} B_{0,x}^{t-1} imes B_{0,y} ]

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

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

[A imes B=sum_{i}sum_{(j+k) mod n=i} A_j imes B_k ]

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

[x_k imes x_t=x_{koplus t} ]

因此(T)的每一行都是方程组(x_k imes x_t=x_{koplus t})的一组解

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

然后这也就说明(T)的行列式不为(0)，因此方程组(x_k imes x_t=x_{koplus t})至少要有(n)组不同的解

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

因此矩阵长这样：

[T=left[egin{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} ight] ]

即(T_{i,j}=omega_n^{ij})

那么它的逆矩阵为：

[T^{-1}=frac{1}{n}left[egin{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} ight] ]

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

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

[(a+bomega) imes (c+domega)=ac+(ad+bc)omega+bd(-omega-1)=(ac-bd)+(ad+bc-bd)omega ]

然后我们由上面的可知：

[T=left[egin{matrix} 1&1&1\ 1&omega&omega^2\ 1&omega^2&omega end{matrix} ight] ]

同时我们知道：

[left[egin{matrix} 1&1&1\ 1&omega&omega^2\ 1&omega^2&omega end{matrix} ight] imes left[egin{matrix} 1&1&1\ 1&omega^2&omega\ 1&omega&omega^2 end{matrix} ight]=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 ot | p)

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

[frac{1}{k+1}+frac{1}{k(k+1)}=frac{1}{k}=frac{p}{3} ]

与假设矛盾！因此得证(3 ot | 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
",f[i].x);
	return 0;
}