bzoj1919: [Ctsc2010]性能优化

bzoj1919: [Ctsc2010]性能优化

权限题贴个题面。
luogu上也有。

分析

循环卷积裸题。
由于$w_n^i=w_n^{i\mathrm{m}\mathrm{o}\mathrm{d}n}$
所以说$w_n^i{w}_n^j=w_n^{(i+j)\mathrm{m}\mathrm{o}\mathrm{d}n}$
发现这刚好是循环卷积的形式。
所以说做长度为$n$的$DFT$之后再$IDFT$就好了。
重点在于这题的$n$的长度不是$2^n$，模数也非常诡异。
所以要介绍Cooley–Tukey FFT algorithm
注意到模数比较小而且是质数，所以直接映射单位复根到模意义下的原根。
回忆快速数论变换的过程。
$e^{\frac{2\pi i}{n}}\equiv g^{\frac{P-1}{n}}$
这个时候我们不这么搞，直接
$e^{\frac{2\pi i}{n}}\equiv g$
带点值的时候直接把$g^k$带进去就好了。
然后最关键的问题来了，$n$不是$2^n$
但是题目给了一个条件$n=2^{k_1}{3}^{k_2}{5}^{k_3}{7}^{k_4}$
$FFT$原本是即将序列分成两半，然后合并。
现在就吧序列分成$k=2,3,5,7$块，然后合并。
比较麻烦的事情就来了——式子得重推。
对于一次分治，假设当前是分成$p$块，当前总长度为$L$，那么分治之后每块的大小就是$d=\frac{L}{p}$。
$A(g^k)=a_0+a_1{g}^k+a_2{g}^{2k}\cdots a_{L-1}{g}^{k(L-1)}$
为了方便，接下来的式子中令$X=g^k$
按照$\mathrm{m}\mathrm{o}\mathrm{d}p$的剩余系分类
$=(a_0+a_p{X}^p+\cdots a_{p(d-1)}{X}^{p(d-1))})+(a_1{X}^1+a_{p+1}{X}^{p+1}\cdots a_{p(d-1)+1}{X}^{p(d-1)+1})+\cdots +(a_{p-1}{X}^{p-1}+a_{p-1+d}{X}^{p-1+d}\cdots a_{pd-1}{X}^{pd-1})$
$$=\sum _{0\le b<p}\sum _{0\le a<d}{a}_{ap+b}{X}^{ap+b}=\sum _{0\le b<p}{X}^b\sum _{0\le a<d}{a}_{ap+b}{X}^{ap}$$
考虑分治往下一层，这个时候稍微注意一下下一层的点值用$w=g^p$的次幂去带。
如果我们已经得到了分治之后的结果，也就是
$A_0(w^m)=a_0+a_p{w}^m+\cdots a_{p(d-1)}{w}^{m(d-1)}={\sum }_{0\le a<d}{a}_{ap}{w}^{ma}$
$A_b(w^m)={\sum }_{0\le a<d}{a}_{ap+b}{w}^{ma}$
$m\in [0,d)$
那么
$$A(g^k)=\sum _{0\le b<p}{g}^{kb}\sum _{0\le a<d}{a}_{ap+b}{g}^{akp}=\sum _{0\le b<p}{g}^{kb}{A}_b(g^{kp})=\sum _{0\le b<p}{g}^{kb}{A}_b(w^k)$$
一波推导，我们得到了一个式子。
$$A(g^k)=\sum _{0\le b<p}{g}^{kb}{A}_b(w^k)$$
美中不足的是$k\in [0,L)$
这个时候令$k=ud+v$
那么
$$A(g^{ud+v})=\sum _{0\le b<p}{g}^{(ud+v)b}{A}_b(w^{(ud+v)})$$
注意到$w^d=g^{pd}=g^L=1$
于是就有
$$A(g^{ud+v})=\sum _{0\le b<p}{g}^{(ud+v)b}{A}_b(w^v)$$
搞定！
果然还是太菜了，网上的大佬似乎都是直接秒出这个式子，只有蒟蒻用markdown推了半天。。。

代码

#include<bits/stdc++.h>
const int N = 5e5 + 10;
int ri() {
	char c = getchar(); int x = 0, f = 1; for(;c < '0' || c > '9'; c = getchar()) if(c == '-') f = -1;
	for(;c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) - '0' + c; return x * f;
}
int a[N], b[N], w[N], R[N], pr[101], tot, P, G, n;
int add(int &a, int b) {return a += b, a >= P ? a - P : a;}
int mul(int a, int b) {return 1LL * a * b % P;}
int Pow(int x, int k) {
	int r = 1;
	for(;k; x = mul(x, x), k >>= 1)
		if(k & 1)
			r = mul(r, x);
	return r;
}
bool Ck(int x, int t) {
	for(int i = 1;i <= tot; ++i)
		if(Pow(x, t / pr[i]) == 1)
			return true;
	return false;
}
void Findori() {
	int x = P - 1;
	for(int i = 2;i <= 7; ++i)
		for(;!(x % i); pr[++tot] = i, x /= i) ;
	G = 2;
	for(;Ck(G, n);) ++G;
}
int Pos(int x, int s, int de, int len) {
	if(de == tot + 1) return s;
	int d = len / pr[de], r = x % pr[de];
	return Pos(x / pr[de], s + d * r, de + 1, d);
}
void DFT(int *A) {
	static int B[N];
	for(int i = 0;i < n; ++i)
		B[R[i]] = A[i];
	int *cur = B, *nxt = A;
	for(int d = 1, de = tot; de; d *= pr[de--], std::swap(cur, nxt)) {
		int pd = d * pr[de], p = n / pd; //w(n,m)^p=w(n/p,m)
		for(int st = 0; st < n; st += pd)
			for(int a = 0; a < pd; a += d) // a * d
				for(int b = 0;b < d; ++b) {
					int &ps = nxt[st + a + b], o = (a + b) * p; ps = 0;
					for(int r = 0; r < pr[de]; ++r)
						ps = add(ps, mul(w[1LL * o * r % n], cur[st + r * d + b]));
				}
	}
	if(cur != A)
		memcpy(A, cur, sizeof(int) * n);
}
int main() {
	n = ri(); int C = ri() % n; P = n + 1;
	Findori();
	for(int i = 0;i < n; ++i)
		a[i] = ri() % P;
	for(int i = 0;i < n; ++i)
		b[i] = ri() % P;
	w[0] = 1;
	for(int i = 1;i < n; ++i)
		w[i] = mul(w[i - 1], G);
	for(int i = 1;i < n; ++i)
		R[i] = Pos(i, 0, 1, n);
	DFT(a); DFT(b);
	for(int i = 0;i < n; ++i)
		a[i] = mul(a[i], Pow(b[i], C));
	DFT(a);
	std::reverse(a + 1, a + n);
	for(int i = 0;i < n; ++i)
		printf("%d\n", mul(a[i], n)); // n * n = 1 (mod n + 1)
	return 0;
}