loj#3058 「HNOI2019」白兔之舞 MTT+单位根反演+矩阵快速幂

最新推荐文章于 2019-05-22 22:00:00 发布

olahiuj

最新推荐文章于 2019-05-22 22:00:00 发布

阅读量209

点赞数

分类专栏： c++ MTT 矩阵快速幂单位根反演文章标签： loj

本文链接：https://blog.csdn.net/jpwang8/article/details/89219514

版权

c++ 同时被 3 个专栏收录

1073 篇文章 0 订阅

订阅专栏

单位根反演

2 篇文章 0 订阅

订阅专栏

MTT

1 篇文章 0 订阅

订阅专栏

Description

自己看吧

Solution

考场上只想到了n=1的O(L)做法，死于姿势不足不会单位根反演.jpg
正解是真的难写，套一个MTT是最恶心的。。

考虑20分怎么做，一个比较naive的dp就是f[i,j]表示走了i步到了第j列的方案数，nL²转移就可以成功签到了

考虑n=1的情况怎么做。当w=1时，可以发现走恰好d步的方案就是一个组合数，那么答案实际上就是 $f\left(d\right)=\sum_{i=0}^{\lfloor\frac{L}{k}\rfloor}{\binom{L}{ki+d}}$
这样就可以获得0分的好成绩

可以发现上柿等价于 $\sum_{i=0}^{L}{{\left[k|\left(i-d\right)\right]}{\binom{L}{i}}w^{i}}$
套一个单位根反演就有
$\frac{1}{k}\sum_{j=0}^{k-1}{\omega_k^{-dj}\sum_{i=0}^{L}{W^i{\omega_k}^{ij}\binom{L}{i}}}$
可以发现后面那一坨可以二项式定理搞成矩阵快速幂，那么就只和j有关了。记 $C_i={\left({\omega_k}^{j}W+I\right)}^{L}$
很显然这里I是单位矩阵，W是给出的邻接矩阵。实际操作我们只需要取出W[x][y]就可以了

现在考虑怎么算 $\sum_{i=0}^{k-1}{{w_k}^{-di}C_i}$
有个骚操作就是构造 $ab=\binom{a+b}{2}-\binom{a}{2}-\binom{b}{2}$ 证明的话可以考虑拆开来直接算

那么上柿子就可以变成
${\omega_k}^{\binom{d}{2}}\sum_{i=0}^{k-1}{\omega_k}^{\binom{i}{2}}C_i\cdot {\omega_k}^{-\binom{d+i}{2}}$

可以发现这样就能构造卷积了，把其中一个翻转然后结果的第k+i位对应原本的第i位

Code

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define fill(x,t) memset(x,t,sizeof(x))

typedef long long LL;
typedef long double ld;
const ld pi=acos(-1);
const int N=524588;
const int M=32767;

struct Mat {
	LL r[4][4];
	Mat() {fill(r,0);}
	LL* operator [](int x) {
		return r[x];
	}
} I,W;

struct com {
	ld r,i;
} ;

int rv[N],v[N],w[N],MOD,n,k;
int A[N],B[N],C[N],D[N];

void upd(LL &x,LL v) {
	x+=v; (x>=MOD)?(x-=MOD):0;
}

com operator +(const com &a,const com &b) {return (com) {a.r+b.r,a.i+b.i};}
com operator -(const com &a,const com &b) {return (com) {a.r-b.r,a.i-b.i};}
com operator *(const com &a,const com &b) {return (com) {a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r};}
com conj(com a) {return (com) {a.r,-a.i};}

int read() {
	int x=0,v=1; char ch=getchar();
	for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):v,ch=getchar());
	for (;ch<='9'&&ch>='0';x=x*10+ch-'0',ch=getchar());
	return x*v;
}

Mat operator *(Mat A,Mat B) {
	Mat C; rep(i,1,3) rep(k,1,3) {
		rep(j,1,3) upd(C[i][j],A[i][k]*B[k][j]%MOD);
	}
	return C;
}

Mat operator *(Mat A,LL b) {
	Mat C;
	rep(i,1,3) rep(j,1,3) C[i][j]=A[i][j]*b%MOD;
	return C;
}

Mat operator +(Mat A,Mat B) {
	Mat C;
	rep(i,1,3) rep(j,1,3) C[i][j]=(A[i][j]+B[i][j])%MOD;
	return C;
}

Mat ksm(Mat A,int dep) {
	Mat C=I;
	for (;dep;dep>>=1,A=A*A) {
		if (dep&1) C=C*A;
	}
	return C;
}

int ksm(int x,int dep) {
	int res=1;
	for (;dep;dep>>=1,x=(LL)x*x%MOD) {
		(dep&1)?(res=(LL)res*x%MOD):0;
	}
	return res;
}

int P(int n) {
	return ((LL)n*(n-1)/2)%k;
}

void FFT(com *a,int n) {
	for (int i=0;i<n;++i) if (i<rv[i]) std:: swap(a[i],a[rv[i]]);
	for (int i=1;i<n;i<<=1) {
		com wn=(com) {cos(pi/i),sin(pi/i)};
		for (int j=0;j<n;j+=(i<<1)) {
			com w=(com) {1,0};
			for (int k=0;k<i;++k) {
				com u=a[j+k],v=a[j+k+i]*w;
				a[j+k]=u+v,a[j+k+i]=u-v;
				w=w*wn;
			}
		}
	}
}

void mul(int *A,int *B,int *C,int n) {
	static com a[N],b[N],Da[N],Db[N],Dc[N],Dd[N];
	for (int i=0;i<n;++i) A[i]=(A[i]+MOD)%MOD;
	for (int i=0;i<n;++i) B[i]=(B[i]+MOD)%MOD;
	for (int i=0;i<n;++i) a[i]=(com) {(ld)(A[i]&M),(ld)(A[i]>>15)};
	for (int i=0;i<n;++i) b[i]=(com) {(ld)(B[i]&M),(ld)(B[i]>>15)};
	FFT(a,n),FFT(b,n);
	for (int i=0;i<n;++i) {
		int j=(n-1)&(n-i);
		com da=(a[i]+conj(a[j]))*(com) {0.5,0};
		com db=(a[i]-conj(a[j]))*(com) {0,-0.5};
		com dc=(b[i]+conj(b[j]))*(com) {0.5,0};
		com dd=(b[i]-conj(b[j]))*(com) {0,-0.5};
		Da[j]=da*dc,Db[j]=da*dd,Dc[j]=db*dc,Dd[j]=db*dd;
	}
	for (int i=0;i<n;++i) a[i]=Da[i]+Db[i]*(com) {0,1};
	for (int i=0;i<n;++i) b[i]=Dc[i]+Dd[i]*(com) {0,1};
	FFT(a,n),FFT(b,n);
	for (int i=0;i<n;++i) {
		LL da=(LL)(a[i].r/n+0.5)%MOD;
		LL db=(LL)(a[i].i/n+0.5)%MOD;
		LL dc=(LL)(b[i].r/n+0.5)%MOD;
		LL dd=(LL)(b[i].i/n+0.5)%MOD;
		C[i]=((LL)da+((LL)((db+dc)%MOD)<<15)+((LL)(dd%MOD)<<30)%MOD)%MOD;
	}
}

int get_g() {
	int tmp=MOD-1; v[0]=0;
	for (int i=2;i*i<=tmp;++i) {
		if (tmp%i==0) {
			while (tmp%i==0) tmp/=i;
			v[++v[0]]=i;
		}
	}
	if (tmp!=1) v[++v[0]]=tmp;
	rep(g,2,MOD) {
		bool flag=false;
		rep(i,1,v[0]) {
			if (ksm(g,(MOD-1)/v[i])==1) {
				flag=true; break;
			}
		}
		if (!flag) return g;
	}
}

int main(void) {
	freopen("data.in","r",stdin);
	int n=read(); k=read(); int L=read();
	int x=read(),y=read(); MOD=read();
	rep(i,0,n) I[i][i]=1;
	rep(i,1,n) rep(j,1,n) W[i][j]=read();
	int g=get_g();
	w[0]=1,w[1]=ksm(g,(MOD-1)/k);
	rep(i,2,k-1) w[i]=(LL)w[i-1]*w[1]%MOD;
	rep(i,0,k-1) {
		C[i]=ksm(W*w[i]+I,L)[x][y];
		A[i]=(LL)w[P(i)]*C[i]%MOD;
	}
	rep(i,0,k*2-1) B[i]=w[(k-P(i))%k];
	rep(i,0,k/2) std:: swap(A[i],A[k-i]);
	int len=1,lg=0; for (;len<=k*3;) len<<=1,lg++;
	for (int i=0;i<len;++i) rv[i]=(rv[i>>1]>>1)|((i&1)<<(lg-1));
	mul(A,B,D,len);
	rep(i,0,k-1) {
		LL tmp=(LL)ksm(k,MOD-2)*w[P(i)]%MOD;
		LL ans=(LL)D[i+k]*tmp%MOD;
		ans=(ans%MOD+MOD)%MOD;
		printf("%lld\n", ans);
	}
	return 0;
}

olahiuj

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
loj#3058 「HNOI2019」白兔之舞 MTT+单位根反演+矩阵快速幂

Description自己看吧Solution考场上只想到了n=1的O(L)做法，死于姿势不足不会单位根反演.jpg正解是真的难写，套一个MTT是最恶心的。。考虑20分怎么做，一个比较naive的dp就是f[i,j]表示走了i步到了第j列的方案数，nL2转移就可以成功签到了考虑n=1的情况怎么做。当w=1时，可以发现走恰好d步的方案就是一个组合数，那么答案实际上就是f(d)=∑i=...
复制链接

扫一扫