FFT
蒟蒻学了 FFT 赶紧做个总结,要不然就忘了...
多项式的表示
首先学习两种多项式的表达法
- 系数表达法
- 点值表达法
系数表达法
比如一个 \(n\) 次多项式 \(A(x)\) 他有 \(n + 1\) 项 于是 设每一项的系数为 \(a_i\) 则有 \(A(x) = \sum _{i = 0}^ n a_i x^ i\)
利用这种方法计算多项式乘法复杂度为 \(O(n^2)\) 即利用(第一个多项式中每个系数都需要与第二个多项式的每个系数相乘)
用这种方法,计算 \(A(x_0)\) 的值是非常好用的,我们使用霍纳法则即可。
点值表达法
我们知道如果给出 \(n + 1\) 个互不相同的点,他们确定一个 \(n\) 次多项式,这个结论是显然的
我们在下文中对于次数的规定稍微放开,规定 A 的次数界 (即最大次数 < 次数界)为 \(K_a\) 而同理,有 \(B_k\) 。通常的,把 \(A(x)*B(x)\) 的次数界规定为 \(K_a + K_b\)
对于许多多项式的操作,点值表达式十分方便,例如有一个点值表达式 A
还有一个点值表达式 B
那么 \(A + B\) 的表达就是
表达法END
复数浅讲
上过高中的请跳过,我是初中生,会讲锅的。
定义一个复数 \(u = a + bi\) 其中 \(i = \sqrt{-1}\) 。你其实可以按向量那么理解。
然后我们复数乘法显得巧妙
比如说一个
复数 \(b\) * 复数 \(c\) 可以理解为就是旋转,比如这个图
看到了吗,旋转,模长相乘。
我们还有一类特殊的向量——\(n\) 次本源单位根 \(\omega\)
-.-。我谔谔,表示什么,就是方程 \(x^n = i\) 的解
这就很显然了
有几个关于本源单位根的定理
- 消去引理: \(\omega ^{dk}_{dn} = \omega ^k_n\)
- 对于任意偶数 \(n > 0\) 有 \(\omega ^{\frac {n}{2}}_n = \omega^1_2 = -1\)
- 折半引理:如果 \(n\) 为偶数,则有 \((\omega^k_n)^2 = (\omega ^{k + n/2}_n)^2\) 从复数相乘的角度上看这是显然的
- 求和引理:对于任意整数 \(n\ge1\) 且 不能被 \(n\) 整除的 \(k\),有 \(\sum _{j=0}^{n-1} (w^k_n )^j = 0\) 为啥,由于我们可以用等比数列求和公式来计算为 \(\sum ^{n-1}_0 (\omega_n^k)^j = \frac {(w^k_n)^n - 1}{w_n^k -1} = 0\)
DFT
离散小波变换 离散傅里叶变换
计算次数界为 \(n\) 的 的多项式 \(A(x) = \sum _0^{n-1} a_j x^j\) 我们不妨带入每个单位根得到 \(y_k = A(\omega_n^k) = \sum _{j = 0} ^ {n-1}a_j \omega _n ^{kj}\)
FFT
快速傅里叶变换、利用单位复数根的性质,我们可以在 \(O(n\log n)\) 的时间内计算 出 \(DFT_n\)
不妨假设都为 \(n\) 以次数界 , 且 \(n\) 为 2 的幂次。
不妨计算 \(A(\omega_n^k)\) 然后按下标的奇偶分类为 \(A_1A_2\)
我们通过观察 有
所以,我们在计算上一个时,可以 \(O(1)\) 得到下一个 于是就很明显了
于是明天再写 逆FFT变换
#include <bits/stdc++.h>
template <typename T>
inline T read()
{
T x = 0;
char ch = getchar();
bool f = 0;
while(ch < '0' || ch > '9')
{
f = (ch == '-');
ch = getchar();
}
while(ch <= '9' && ch >= '0')
{
x = (x << 1) + (x << 3) + ch - '0';
ch = getchar();
}
return f ? -x: x;
}
template <typename T>
void put(T x)
{
if(x < 0) x = -x, putchar('-');
if(x < 10)
{
putchar(x + 48);
return ;
}
put(x / 10);
putchar(x % 10 + 48);
return;
}
#define rd read <int>
#define pt(x) put <int> (x), putchar(' ')
const int Maxn = 1e6 + 111;
int n, m, limit;
typedef double db;
namespace myspace
{
struct complex
{
db x, y;
complex(db xx, db yy){ x = xx, y = yy; }
complex(){ x = 0, y = 0; }
complex operator * (const complex &a) { return complex(x * a.x - y * a.y, x * a.y + y * a.x); }
complex operator + (const complex &a) { return complex(x + a.x , y + a.y); }
complex operator - (const complex &a) { return complex(x - a.x , y - a.y); }
db abs() { return sqrt(x * x + y * y); }
};
};
myspace::complex a[Maxn], b[Maxn];
const db Pi = 3.1415926535897;
void fft(int limit, myspace::complex *a, int type)
{
if(limit == 1) return;
myspace::complex a1[(limit >> 1)], a2[(limit >> 1)];
// 写的很巧妙,注意一定是 < limit
for(int i = 0; i < limit; i += 2)
{
a1[(i >> 1)] = a[i];
a2[(i >> 1)] = a[i + 1];
}
fft(limit >> 1, a1, type);
fft(limit >> 1, a2, type);
myspace::complex w(1, 0), wn(cos(2.0 * Pi / limit), type * sin(2.0 * Pi / limit));
// 记录本源单位根,同时计算a
for(int i = 0; i < limit >> 1; ++i, w = w * wn)
a[i] = a1[i] + w * a2[i], a[i + (limit >> 1)] = a1[i] - w * a2[i];
return;
}
int main()
{
#ifdef _DEBUG
freopen("in.txt", "r", stdin);
// freopen("out.txt", "w", stdout);
#endif
n = rd();
m = rd();
for(int i = 0; i <= n; ++i) a[i].x = rd();
for(int i = 0; i <= m; ++i) b[i].x = rd();
limit = 1;
while(limit <= n + m) limit <<= 1;
fft(limit, a, 1);
fft(limit, b, 1);
for(int i = 0; i <= limit; ++i) a[i] = a[i] * b[i];
fft(limit, a, -1);
for(int i = 0; i <= n + m; ++i) pt((int) (a[i].x / limit + 0.5));
putchar('\n');
return 0;
}