【HNOI2019】白兔之舞（单位根反演）（矩阵快速幂）（MTT）（循环卷积优化常数）

最新推荐文章于 2022-06-25 21:26:36 发布

zxyoi_dreamer

最新推荐文章于 2022-06-25 21:26:36 发布

阅读量209

点赞数

分类专栏：单位根反演矩阵快速幂

本文链接：https://blog.csdn.net/zxyoi_dreamer/article/details/100528003

版权

矩阵快速幂同时被 2 个专栏收录

18 篇文章 0 订阅

订阅专栏

单位根反演

5 篇文章 0 订阅

订阅专栏

传送门

题解：

~~说个笑话，这道题很休闲的~~

~~真正的好题，会让你算t等于所有情况的答案~~

考虑走 $i$ 步的答案，显然就是 $A^i[x][y]{L\choose i}$ ，其中 $A^i$ 表示原矩阵的 $i$ 次幂，其中 $[x] [y]$ 表示取矩阵中间的某个元素， ${L\choose i}$ 是因为需要在中间选择 $i$ 个地方落脚。

那么考虑模 $k$ 为某个定值的方案，显然是单位根反演，为了方便书写，以下矩阵中某个元素求和部分省略，以矩阵求和代替。

$\begin{aligned} ans_t=&\sum_{i\%k=t}A^i{L\choose i}\\ =&\sum_{i=0}^L[k\mid i-t]A^i{L\choose i}\\ =&\sum_{i=0}^L\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{j(i-t)}A^i{L\choose i}\\ =&\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-jt}\sum_{i=0}^L(\omega_k^jA)^i{L\choose i}\\ =&\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-jt}(\omega_k^jA+I)^L \end{aligned}$

显然后面的显然后面的矩阵可以直接快速幂，并且我们只会用到一个位置上的值，设 $a_j=(\omega_k^jA+I)^L[x][y]$

继续化简：

$\begin{aligned} ans_t=&\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{-jt}a_j \end{aligned}$

这个已经可以 $O(k^2)$ 愉快搞掉了。多点求值是被卡了的。

现在需要考虑对于所有 $t$ 算这个东西。

显然我们需要把 $t$ 从指数上面和 $j$ 搞分开。

考虑BlueStein里面构造单位根的方式 $ij=\frac{i^2+j^2-(i-j)^2}{2}$ ，然而我们并不敢保证 $i^2/2$ 一定是一个整数。

似乎陷入了僵局。

然鹅，我们考虑这个式子 $-jt={j\choose 2}+{t\choose 2}-{j+t\choose 2}$

好的，组合数永远是整数， ~~喜闻乐见，我们赢了~~。

所以上面这个东西：

$ans_t=\frac{t\choose 2}{k}\sum_{i=0}^{k-1}\omega^{j\choose 2}a_j\cdot \omega_k^{-{j+t\choose 2}}$

好的，翻转，构造卷积，~~我们赢了~~

考虑一个不重要的常数优化。

由于是一个 $k$ 次多项式和一个 $2 * k - 2$ 次多项式做卷积。

只会用到结果的 $k$ 到 $2 k - 1$ 次，可以把卷积数组开到 $2 k$ ，循环卷积到前面的值可以直接丢掉。

~~然后写一个悠闲的MTT，一发入魂~~

细节不是很多，但是我矩阵乘法写挂了，~~叫你循环展开~~。。。

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;

int n,k,L,x,y,mod;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}

cs int SIZE=1<<17|1;
int ome[SIZE],iome[SIZE],C[SIZE],gr;


inline void get_gr(){
	int phi=mod-1;std::vector<int> p;
	for(int re i=2;i*i<=phi;++i)if(phi%i==0){
		p.push_back(i);
		while(phi%i==0)phi/=i;
	}
	if(phi>1)p.push_back(phi);phi=mod-1;
	for(gr=2;;++gr){
		bool flag=true;
		for(int t:p)if(power(gr,phi/t)==1){
			flag=false;break;
		}
		if(flag)break;
	}
}

struct matrix{
	int data[3][3];
	matrix(){}
	cs int *operator[](int of)cs{return data[of];}
	int *operator[](int of){return data[of];}
};
#define p mod
#define int64 ll
inline matrix operator+(cs matrix &A,cs matrix &B){
	matrix C;
	switch(n){
		case 3:
			C[0][0]=add(A[0][0],B[0][0]);C[0][1]=add(A[0][1],B[0][1]);C[0][2]=add(A[0][2],B[0][2]);
			C[1][0]=add(A[1][0],B[1][0]);C[1][1]=add(A[1][1],B[1][1]);C[1][2]=add(A[1][2],B[1][2]);
			C[2][0]=add(A[2][0],B[2][0]);C[2][1]=add(A[2][1],B[2][1]);C[2][2]=add(A[2][2],B[2][2]);
			break;
		case 2:
			C[0][0]=add(A[0][0],B[0][0]);C[0][1]=add(A[0][1],B[0][1]);
			C[1][0]=add(A[1][0],B[1][0]);C[1][1]=add(A[1][1],B[1][1]);
			break;
		case 1:
			C[0][0]=add(A[0][0],B[0][0]);
	}
	return C;
}

inline matrix operator*(cs matrix &A,int b){
	matrix C;
	switch(n){
		case 3:
			C[0][0]=mul(A[0][0],b);C[0][1]=mul(A[0][1],b);C[0][2]=mul(A[0][2],b);
			C[1][0]=mul(A[1][0],b);C[1][1]=mul(A[1][1],b);C[1][2]=mul(A[1][2],b);
			C[2][0]=mul(A[2][0],b);C[2][1]=mul(A[2][1],b);C[2][2]=mul(A[2][2],b);
			break;
		case 2:
			C[0][0]=mul(A[0][0],b);C[0][1]=mul(A[0][1],b);
			C[1][0]=mul(A[1][0],b);C[1][1]=mul(A[1][1],b);
			break;
		case 1:
			C[0][0]=mul(A[0][0],b);
	}
	return C;
}

inline matrix operator*(cs matrix &A,cs matrix &B){
	matrix C;
	switch(n){
		case 3:
			C[0][0]=((ll)A[0][0]*B[0][0]+(ll)A[0][1]*B[1][0]+(ll)A[0][2]*B[2][0])%mod;
			C[0][1]=((ll)A[0][0]*B[0][1]+(ll)A[0][1]*B[1][1]+(ll)A[0][2]*B[2][1])%mod;
			C[0][2]=((ll)A[0][0]*B[0][2]+(ll)A[0][1]*B[1][2]+(ll)A[0][2]*B[2][2])%mod;
			C[1][0]=((ll)A[1][0]*B[0][0]+(ll)A[1][1]*B[1][0]+(ll)A[1][2]*B[2][0])%mod;
			C[1][1]=((ll)A[1][0]*B[0][1]+(ll)A[1][1]*B[1][1]+(ll)A[1][2]*B[2][1])%mod;
			C[1][2]=((ll)A[1][0]*B[0][2]+(ll)A[1][1]*B[1][2]+(ll)A[1][2]*B[2][2])%mod;
			C[2][0]=((ll)A[2][0]*B[0][0]+(ll)A[2][1]*B[1][0]+(ll)A[2][2]*B[2][0])%mod;
			C[2][1]=((ll)A[2][0]*B[0][1]+(ll)A[2][1]*B[1][1]+(ll)A[2][2]*B[2][1])%mod;
			C[2][2]=((ll)A[2][0]*B[0][2]+(ll)A[2][1]*B[1][2]+(ll)A[2][2]*B[2][2])%mod;
			break;
		case 2:
			C[0][0]=((ll)A[0][0]*B[0][0]+(ll)A[0][1]*B[1][0])%mod;
			C[0][1]=((ll)A[0][0]*B[0][1]+(ll)A[0][1]*B[1][1])%mod;
			C[1][0]=((ll)A[1][0]*B[0][0]+(ll)A[1][1]*B[1][0])%mod;
			C[1][1]=((ll)A[1][0]*B[0][1]+(ll)A[1][1]*B[1][1])%mod;
			break;
		case 1:
			C[0][0]=(ll)A[0][0]*B[0][0]%mod;
	}
	return C;
}

inline int power(matrix A,int b,int r,int c){
	int res[3]={};res[r]=1;
	while(b){
		if(b&1){
			int tp[3]={};
			switch(n){
				case 3:
					tp[0]=((ll)res[0]*A[0][0]+(ll)res[1]*A[1][0]+(ll)res[2]*A[2][0])%mod;
					tp[1]=((ll)res[0]*A[0][1]+(ll)res[1]*A[1][1]+(ll)res[2]*A[2][1])%mod;
					tp[2]=((ll)res[0]*A[0][2]+(ll)res[1]*A[1][2]+(ll)res[2]*A[2][2])%mod;
					res[0]=tp[0];res[1]=tp[1];res[2]=tp[2];
					break;
				case 2:
					tp[0]=((ll)res[0]*A[0][0]+(ll)res[1]*A[1][0])%mod;
					tp[1]=((ll)res[0]*A[0][1]+(ll)res[1]*A[1][1])%mod;
					res[0]=tp[0];res[1]=tp[1];
					break;
				case 1:
					res[0]=(ll)res[0]*A[0][0]%mod;
			}
		}
		A=A*A;b>>=1;
	}
	return res[c];
}

matrix A,I;

namespace MTT{
	cs double PI=acos(-1);
	struct cp{
		double x,y;
		cp(){}
		cp(double _x,double _y):x(_x),y(_y){}
		cp conj()cs{return cp(x,-y);}
		friend cp operator+(cs cp &a,cs cp &b){return cp(a.x+b.x,a.y+b.y);}
		friend cp operator-(cs cp &a,cs cp &b){return cp(a.x-b.x,a.y-b.y);}
		friend cp operator*(cs cp &a,cs cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
		friend cp operator/(cs cp &a,double b){return cp(a.x/b,a.y/b);}
	};
	cp w[SIZE];int r[SIZE];
	inline void init_rev(int l){
		cp wn(cos(2*PI/l),sin(2*PI/l));
		for(int re i=0;i<l/2;++i)
		w[l/2+i]=(i&31)?w[l/2+i-1]*wn:cp(cos(2*PI*i/l),sin(2*PI*i/l));
		for(int re i=l/2-1;i;--i)w[i]=w[i<<1];
		for(int re i=1;i<l;++i)r[i]=r[i>>1]>>1|((i&1)?l>>1:0);
	}
	inline void DFT(cp *A,int l){
		for(int re i=1;i<l;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
		for(int re i=1;i<l;i<<=1)
		for(int re j=0;j<l;j+=i<<1)
		for(int re k=0;k<i;++k){
			cp t=A[j+k+i]*w[i+k];
			A[j+k+i]=A[j+k]-t;
			A[j+k]=A[j+k]+t;
		}
	}
	void mul(int *a,int *b,int len){
		static cp A[SIZE],B[SIZE],C[SIZE],D[SIZE];
		for(int re i=0;i<len;++i){
			A[i]=cp(a[i]&0x7fff,a[i]>>15);
			B[i]=cp(b[i]&0x7fff,b[i]>>15);
		} 
		init_rev(len),DFT(A,len),DFT(B,len);
		for(int re i=0;i<len;++i){
			int u=(len-i)&(len-1);
			C[i]=(A[i].conj()+A[u])*cp(0.5,0)*B[u];
			D[i]=(A[i].conj()-A[u])*cp(0,0.5)*B[u];
		}
		DFT(C,len),DFT(D,len);
		for(int re i=0;i<len;++i){
			ll x=C[i].x/len+0.5,y=(C[i].y+D[i].x)/len+0.5,z=D[i].y/len+0.5;
			x%=mod,y%=mod,z%=mod;
			a[i]=(x+(y<<15)+(z<<30))%mod;
		}
	}
}

int a[SIZE],b[SIZE],c[SIZE];

signed main(){
#ifdef zxyoi
	freopen("dance.in","r",stdin);
#endif
	scanf("%d%d%d%d%d%d",&n,&k,&L,&x,&y,&mod);--x,--y;
	for(int re i=0;i<n;++i)
	for(int re j=0;j<n;++j)scanf("%d",&A[i][j]),I[i][j]=(i==j);
	get_gr();
	iome[0]=ome[0]=1;ome[1]=power(gr,(mod-1)/k);
	for(int re i=2;i<k;++i)ome[i]=mul(ome[i-1],ome[1]);
	for(int re i=1;i<k;++i)iome[i]=ome[k-i];
	for(int re i=1;i<=2*k-2;++i)C[i]=(ll)i*(i-1)/2%k;
	for(int re i=0;i<k;++i)c[i]=mul(power(A*ome[i]+I,L,x,y),ome[C[i]]);
	std::reverse(c,c+k+1);
	for(int re i=0;i<=2*k-2;++i)a[i]=iome[C[i]];
	int len=1;while(len<2*k)len<<=1;
	MTT::mul(a,c,len);
	int inv_k=power(k,mod-2);
	for(int re i=0;i<k;++i)cout<<mul(a[k+i],mul(inv_k,ome[C[i]]))<<"\n";
	return 0;
}

zxyoi_dreamer

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
【HNOI2019】白兔之舞（单位根反演）（矩阵快速幂）（MTT）（循环卷积优化常数）

传送门题解：说个笑话，这道题很休闲的真正的好题，会让你算t等于所有情况的答案考虑走iii步的答案，显然就是Ai[x][y](Li)A^i[x][y]{L\choose i}Ai[x][y](iL)，其中AiA^iAi表示原矩阵的iii次幂，其中[x][y][x][y][x][y]表示取矩阵中间的某个元素，(Li){L\choose i}(iL)是因为需要在中间选择iii个地方落脚。...
复制链接

扫一扫