分治 FFT学习笔记

先给一道luogu板子题:P4721 【模板】分治 FFT


今天模拟有道题的部分分做法是分治fft,于是就学了一下。感觉不是很难,国赛上如果推出式子的话应该能写出来。


分治fft用来解决这么一个式子$$f(i) = \sum _ {j = 1} ^ {i} f(j) * g(i - j)$$
如果暴力fft的话,复杂度\(O(n ^ 2logn)\)还没有暴力优秀。


我们可以用cdq分治的思想,对于区间\([L, R]\),假设\([L, mid]\)已经求出,现在要算\([mid + 1, R]\)。
那么我们考虑对于\([mid + 1, R]\)中的某一项\(x\),前面\([L, mid]\)对他的贡献:$$f(x) = \sum _{i = L} ^ {mid} f(i) * g(x - i)$$
为了方便,我们把循环到\(mid\)改为\(x\),反正\(mid\)~\(x\)这些项还没有被计算,暂且为0。于是有

\[\begin{align*} f(x) &= \sum _ {i = L} ^ {x} f(i)* g(x - i) \\ &= \sum _ {i = 0} ^ {x - L} f(i + L) * g(x - L - i) \end{align*}\]

我们令\(A(i) = f(i + L), B(i) = g(i)\),于是式子就变成了$$f(x) = C(x - L) = \sum _ {i = 0} ^{x - L} A(i) *B(x - L - i)$$
这个时候能看出后面是一个卷积的形式,直接FFT即可。然后把\(C(x - L)\)加到\(f(x)\)上,左边对右边的贡献就算完了。再递归右边。


递归的每一层是\(O(len log len)\)的(区间长度),有\(logn\)层,因此总时间复杂度为\(O(nlog ^ 2n)\)。

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<queue>
#include<assert.h>
#include<ctime>
using namespace std;
#define enter puts("") 
#define space putchar(' ')
#define Mem(a, x) memset(a, x, sizeof(a))
#define In inline
#define forE(i, x, y) for(int i = head[x], y; ~i && (y = e[i].to); i = e[i].nxt)
typedef long long ll;
typedef double db;
const int INF = 0x3f3f3f3f;
const db eps = 1e-8;
const int maxn = 3e5 + 5;
const ll mod = 998244353;
const ll G = 3;
inline ll read()
{
	ll ans = 0;
	char ch = getchar(), last = ' ';
	while(!isdigit(ch)) last = ch, ch = getchar();
	while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar();
	if(last == '-') ans = -ans;
	return ans;
}
inline void write(ll x)
{
	if(x < 0) x = -x, putchar('-');
	if(x >= 10) write(x / 10);
	putchar(x % 10 + '0');
}
In void MYFILE()
{
#ifndef mrclr
	freopen("ha.in", "r", stdin);
	freopen("ha.out", "w", stdout);
#endif
}

int n;
ll f[maxn], g[maxn];

In ll inc(ll a, ll b) {return a + b < mod ? a + b : a + b - mod;}
In ll quickpow(ll a, ll b)
{
	ll ret = 1;
	for(; b; b >>= 1, a = a * a % mod)
		if(b & 1) ret = ret * a % mod;
	return ret;
}

int rev[maxn];
ll A[maxn], B[maxn];
In void ntt(ll* a, int len, bool flg)
{
	for(int i = 0; i < len; ++i) if(i > rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1; i < len; i <<= 1)
	{
		ll gn = quickpow(G, (mod - 1) / (i << 1));
		for(int j = 0; j < len; j += (i << 1))
		{
			ll g = 1;
			for(int k = 0; k < i; ++k, g = g * gn % mod)
			{
				ll tp1 = a[j + k], tp2 = a[j + k + i] * g % mod;
				a[j + k] = inc(tp1, tp2), a[j + k + i] = inc(tp1, mod - tp2);
			}
		}
	}
	if(flg) return;
	reverse(a + 1, a + len); ll inv = quickpow(len, mod - 2);
	for(int i = 0; i < len; ++i) a[i] = a[i] * inv % mod;
}

In void fftSolve(int L, int R)
{
	if(L == R) return;
	int mid = (L + R) >> 1;
	fftSolve(L, mid);
	int n = R - L + 1, len = 1, lim = 0;
	while(len <= n + n) len <<= 1, ++lim;
	for(int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1));
	fill(A, A + len, 0), fill(B, B + len, 0);
	for(int i = L; i <= mid; ++i) A[i - L] = f[i];
	for(int i = 0; i < n; ++i) B[i] = g[i];
	ntt(A, len, 1), ntt(B, len, 1);
	for(int i = 0; i < len; ++i) A[i] = A[i] * B[i] % mod;
	ntt(A, len, 0);
	for(int i = mid + 1; i <= R; ++i) f[i] = inc(f[i], A[i - L]);
	fftSolve(mid + 1, R);
}

int main()
{
//	MYFILE();
	n = read();
	for(int i = 1; i < n; ++i) g[i] = read();
	f[0] = 1, fftSolve(0, n - 1);
	for(int i = 0; i < n; ++i) write(f[i]), space; enter;
	return 0;	
}
上一篇:分治fft


下一篇:2021杭电多校第7场1011(任意模数,分治FFT)