FFT学习笔记
总结
单位根
$n$次单位根$w_n^i$满足$(w_{n}^i)^n=1$
$\text{DFT}$的加速版本
设:
$$ f(x)=\sum_{i=1}^{n-1}a_ix^i \\ f^{[0]}(x)=a_0 + a_2x+...+a_{n-2}x^{n/2-1} \\ f^{[1]}(x)=a_1 + a_3x+...+a_{n-1}x^{n/2-1} \\ \therefore f(x) = f^{[0]}(x^2)+xf^{[1]}(x^2) \\ \therefore f(w_n^k)=f^{[0]}(w_{n/2}^k)+w_n^kf^{[1]}(w_{n/2}^k) $$
同理:
$$ \therefore f(w_n^{k+n/2})=f^{[0]}(w_{n/2}^k)-w_n^kf^{[1]}(w_{n/2}^k) $$
这样可以递归计算。
$\text{NTT}$
$$ \because w_n \equiv g^{\frac{p-1}n} $$
于是将单位根换成原根即可
代码
$\text{FFT}$
#include<cstdio>
#include<cstring>
#include<cctype>
#include<cmath>
#include<algorithm>
#define RG register
#define file(x) freopen(#x".in", "r", stdin);freopen(#x".out", "w", stdout);
#define clear(x, y) memset(x, y, sizeof(x))
inline int read()
{
int data = 0, w = 1; char ch = getchar();
while(ch != '-' && (!isdigit(ch))) ch = getchar();
if(ch == '-') w = -1, ch = getchar();
while(isdigit(ch)) data = data * 10 + (ch ^ 48), ch = getchar();
return data * w;
}
const int maxn(3e6 + 10);
const double pi(acos(-1));
struct complex { double x, y; } a[maxn], b[maxn];
inline complex operator + (const complex &lhs, const complex &rhs)
{ return (complex) {lhs.x + rhs.x, lhs.y + rhs.y}; };
inline complex operator - (const complex &lhs, const complex &rhs)
{ return (complex) {lhs.x - rhs.x, lhs.y - rhs.y}; };
inline complex operator * (const complex &lhs, const complex &rhs)
{
return (complex) {lhs.x * rhs.x - lhs.y * rhs.y,
lhs.y * rhs.x + lhs.x * rhs.y};
}
int n, m, r[maxn], P;
template<int opt> void FFT(complex *p)
{
for(RG int i = 0; i < n; i++) if(i < r[i]) std::swap(p[i], p[r[i]]);
for(RG int i = 1; i < n; i <<= 1)
{
complex rot = (complex) {cos(pi / i), opt * sin(pi / i)};
for(RG int j = 0; j < n; j += (i << 1))
{
complex w = (complex) {1, 0};
for(RG int k = 0; k < i; ++k, w = w * rot)
{
complex x = p[j + k], y = w * p[i + j + k];
p[j + k] = x + y, p[i + j + k] = x - y;
}
}
}
}
int main()
{
#ifndef ONLINE_JUDGE
file(cpp);
#endif
n = read(), m = read();
for(RG int i = 0; i <= n; i++) a[i].x = read();
for(RG int i = 0; i <= m; i++) b[i].x = read();
for(m += n, n = 1; n <= m; n <<= 1, ++P);
for(RG int i = 0; i < n; i++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (P - 1));
FFT<1>(a), FFT<1>(b);
for(RG int i = 0; i < n; i++) a[i] = a[i] * b[i];
FFT<-1>(a);
for(RG int i = 0; i <= m; i++) printf("%d ", (int)(a[i].x / n + .5));
return 0;
}
$\text{NTT}$
#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#define RG register
#define file(x) freopen(#x".in", "r", stdin);freopen(#x".out", "w", stdout);
#define clear(x, y) memset(x, y, sizeof(x))
inline int read()
{
int data = 0, w = 1; char ch = getchar();
while(ch != '-' && (!isdigit(ch))) ch = getchar();
if(ch == '-') w = -1, ch = getchar();
while(isdigit(ch)) data = data * 10 + (ch ^ 48), ch = getchar();
return data * w;
}
const int maxn(3e6 + 10), g(3), Mod(998244353), phi(Mod - 1);
inline int fastpow(int x, int y)
{
int ans = 1;
while(y)
{
if(y & 1) ans = 1ll * ans * x % Mod;
x = 1ll * x * x % Mod, y >>= 1;
}
return ans;
}
int n, m, r[maxn], P, a[maxn], b[maxn];
template<int opt> void FFT(int *p)
{
for(RG int i = 0; i < n; i++) if(i < r[i]) std::swap(p[i], p[r[i]]);
for(RG int i = 1; i < n; i <<= 1)
{
int rot = fastpow(g, phi / (i << 1));
for(RG int j = 0; j < n; j += (i << 1))
{
int w = 1;
for(RG int k = 0; k < i; ++k, w = 1ll * w * rot % Mod)
{
int x = p[j + k], y = 1ll * w * p[i + j + k] % Mod;
p[j + k] = (x + y) % Mod, p[i + j + k] = (x - y + Mod) % Mod;
}
}
}
if(opt == -1) std::reverse(p + 1, p + n);
}
int main()
{
#ifndef ONLINE_JUDGE
file(cpp);
#endif
n = read(), m = read();
for(RG int i = 0; i <= n; i++) a[i] = read();
for(RG int i = 0; i <= m; i++) b[i] = read();
for(m += n, n = 1; n <= m; n <<= 1, ++P);
for(RG int i = 0; i < n; i++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (P - 1));
FFT<1>(a), FFT<1>(b);
for(RG int i = 0; i < n; i++) a[i] = 1ll * a[i] * b[i] % Mod;
FFT<-1>(a);
int inv = fastpow(n, Mod - 2);
for(RG int i = 0; i <= m; i++) printf("%lld ", 1ll * a[i] * inv % Mod);
return 0;
}