快速傅里叶变换求多项式乘法
一、例题
题目背景
这是一道多项式乘法模板题。
题目描述
给定一个 \(n\) 次多项式 \(F(x)\)，和一个 \(m\) 次多项式 \(G(x)\)。
请求出 \(F(x)\) 和 \(G(x)\) 的卷积。
输入格式
第一行两个整数 \(n,m\)。
接下来一行 \(n+1\) 个数字，从低到高表示 \(F(x)\) 的系数。
接下来一行 \(m+1\) 个数字，从低到高表示 \(G(x)\) 的系数。
输出格式
一行 \(n+m+1\) 个数字，从低到高表示 \(F(x) \cdot G(x)\) 的系数。
输入数据
1 2
1 2
1 2 1
输出数据
1 4 5 2
提示
保证输入中的系数大于等于 \(0\) 且小于等于 \(9\)。
对于 \(100\%\) 的数据：\(1 \le n, m \leq {10}^6\)。
二、公式推导
1、思路
可分为 \(3\) 步求解：
（1）分别对两个多项式插入 \(n+m+1\) 个点，求出各自多项式的点值（FFT）
（2）将对应的点值相乘，得到相乘多项式的 \(n+m+1\) 个点值
（3）通过点值还原该多项式（逆FFT）
通过 \(FFT\) 可以在 \(O(nlog\;n)\) 的时间内完成 \(n\) 个点的多项式插值计算，当然这 \(n\) 个点是 \(1\) 在复平面上的 \(n\) 次方根，满足幂的 \(n\) 次周期循环关系，利用该关系，通过二分来实现算法。
2、向量计算
给定向量
-
向量和
\(a+b=(a_0+b_0,a_1+b_1,\cdots,a_{n-1}+b_{n-1})\)
-
内积
\(a\cdot b=a_0b_0+a_1b_1+\cdots+a_{n-1}b_{n-1}\)
-
卷积
\(a\ast b=(c_0,c_1,\cdots,c_{2n-2})\)
\(c_k=\sum_{i+j=k(i,j\lt n)}a_ib_j,k=0,1,\cdots,2n-2\)
3、FFT 变换——分解
任给一个多项式 \(A(x)=a_0+a_1x+a_2x^2\dots+a_{n-1}x^{n-1}\)，令
则有
可得 \(A(x)=A_1(x^2)+xA_2(x^2)\)
在复数域上对 \(1\) 开 \(2n\) 次方，得到 \(2n\) 个复数根，可表示为：
考虑将 \(\omega_{2n}^k\) 作为 \(x\) 代入 \(A(x)\) ，
因此，可以通过分治策略，通过二分归并求解。
4、FFT 逆变换——构造
令 \(y_k=A(\omega_n^k)\)，可以构造
由等比数列求和公式分析可得，
因此，可得 \(c_k=n\cdot a_j\)，通过构造 \(B(x)=y_0+y_1x+y_2x^2+\cdots y_{n-1}x^{n-1}\)，可知 \(c_k=B(\omega_n^{-k})\)，从而求出 \(a_k\) 即可。
三、代码
#include <iostream>
#include <cmath>
#include <algorithm>
const double PI = acos(-1);
const int SZ = 3e6 + 7;
template <typename T>
struct complex {
T r, m;
complex(T _r=0, T _m=0): r(_r), m(_m) {}
complex operator + (const complex &b) {
return complex(r+b.r, m+b.m);
}
complex operator - (const complex &b) {
return complex(r-b.r, m-b.m);
}
//(a+ib)*(c+id)=ac-bd + i(ad+bc)
complex operator * (const complex &b) {
return complex(r*b.r-m*b.m, r*b.m+m*b.r);
}
};
complex<double> a[SZ], b[SZ], c[SZ<<1];
int pos[SZ<<1];
//op为1是正变换，-1是逆变换
void FFT(complex<double> a[], int n, int op) {
for (int i = 1; i < n; ++i)
if (i < pos[i]) std::swap(a[i], a[pos[i]]);
for (int i = 2; i <= n; i <<= 1) {
complex<double> wn(cos(2*PI*op/i), sin(2*PI*op/i));
for (int j = 0; j < n; j += i) {
complex<double> w(1, 0);
for (int k = j, end = j+i/2; k < end; ++k) {
complex<double> t1 = a[k], t2 = w * a[k+i/2];
a[k] = t1 + t2, a[k+i/2] = t1 - t2;
w = w * wn;
}
}
}
if (op == -1) for (int i = 0; i < n; ++i) a[i].r /= n;
}
int main() {
int n, m, p = 0;
//输入多项式的最高幂，以及系数
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; ++i) scanf("%lf", &a[i]);
for (int i = 0; i <= m; ++i) scanf("%lf", &b[i]);
m += n, n = 1;
//寻找合适的2的幂，赋值为n
while (n <= m) n <<= 1, ++p; --p;
for (int i = 1; i < n; ++i) pos[i] = (pos[i>>1]>>1) | ((i&1)<<p);
//正变换
FFT (a, n, 1), FFT (b, n, 1);
//0~n-1,n为2的整数幂
for (int i = 0; i < n; ++i) c[i] = a[i] * b[i];
//逆变换
FFT (c, n, -1);
for (int i = 0; i <= m; ++i) printf("%d ", (int)(c[i].r+0.5));
return 0;
}
🎯 目标：给定
计算 \(A(x)\ast B(x)\) 的乘积。
直接计算：
期望结果系数：\([2, 5, 2]\)
第一步：确定长度 \(N\)
\(deg(A) = 1，deg(B) = 1 \Rightarrow deg(C) = 2\)
需要至少 \(3\) 个系数 \(\Rightarrow\) 取最小 2 的幂 \(\ge 3\Rightarrow N=4\)
补零后：
\(a = [1, 2, 0, 0]\) （对应 \(1 + 2x\)）
\(b = [2, 1, 0, 0]\) （对应 \(2 + x\)）
第二步：手动执行 \(FFT(a) = FFT([1, 2, 0, 0])\)
我们将显式写出每一层递归。
🔁 \(FFT(a):\) 输入 \([1, 2, 0, 0]\)，长度 \(n = 4\)
▶ Level 0（\(n=4\)）
偶下标（0,2）：\(a_{\text{even}} = [1, 0]\)
奇下标（1,3）：\(a_{\text{odd}} = [2, 0]\)
→ 递归计算 \(FFT([1,0])\) 和 \(FFT([2,0])\)
▶ Level 1（\(n=2\)）
子问题 1：\(FFT([1, 0])\)
长度 2，\(\omega_2 = -1\)
输出：\([1+0, 1-0] = [1, 1]\)
子问题 2：\(FFT([2, 0])\)
输出：\([2+0, 2-0] = [2, 2]\)
所以：
\(A_e = [1, 1]\)
\(A_o = [2, 2]\)
▶ 合并（\(n=4，\omega_4 = i\)）
需要 \(\omega_4^0 = 1，\omega_4^1 = i\)
对 \(k = 0, 1：\)
\(k = 0:\)
\(X[0] = A_e[0] + \omega_4^0 \cdot A_o[0] = 1 + 1·2 = 3\)
\(X[2] = A_e[0] - \omega_4^0 \cdot A_o[0] = 1 - 2 = -1\)
\(k = 1:\)
\(X[1] = A_e[1] + \omega_4^1 \cdot A_o[1] = 1 + i·2 = 1 + 2i\)
\(X[3] = A_e[1] - \omega_4^1 \cdot A_o[1] = 1 - 2i\)
✅ 所以：
\(\text{FFT}(a) = [3, 1+2i, -1, 1-2i]\)
✨ 注意：偶部是 \([1,0]\)，奇部是 \([2,0]\) —— 完全不同！清晰体现分治。
第三步：手动执行 \(FFT(b) = FFT([2, 1, 0, 0])\)
▶ Level 0（\(n=4\)）
偶部（0,2）：\(b_e = [2, 0]\)
奇部（1,3）：\(b_o = [1, 0]\)
→ 递归计算 \(FFT([2,0])\) 和 \(FFT([1,0])\)
▶ Level 1（\(n=2\)）
\(FFT([2,0]) = [2, 2]\)
\(FFT([1,0]) = [1, 1]\)
所以：
\(B_e = [2, 2]\)
\(B_o = [1, 1]\)
▶ 合并（\(n=4，\omega_4 = i\)）
k = 0:
\(Y[0] = 2 + 1·1 = 3\)
\(Y[2] = 2 - 1 = 1\)
k = 1:
\(Y[1] = 2 + i·1 = 2 + i\)
\(Y[3] = 2 - i\)
✅ 所以：
\(\text{FFT}(b) = [3, 2+i, 1, 2-i]\)
第四步：频域逐点相乘
\(\hat{c}[k] = \text{FFT}(a)[k] \times \text{FFT}(b)[k]\)
逐项计算：
\(k=0: 3 × 3 = 9\)
\(k=1: (1+2i)(2+i) = 1·2 + 1·i + 2i·2 + 2i·i = 2 + i + 4i + 2i^2 = 2 + 5i - 2 = 5i\)
\(k=2: (-1) × 1 = -1\)
\(k=3: (1-2i)(2-i) = 2 - i -4i + 2i^2 = 2 -5i -2 = -5i\)
所以：
\(\hat{c} = [9, 5i, -1, -5i]\)
第五步：\(IFFT(\hat{c})\) → 得到系数
我们用 FFT 实现 IFFT：
\(IFFT(x) = (1/N) · conjugate( FFT( conjugate(x) ) )\)
先取共轭（因为有虚数）：
\(\overline{\hat{c}} = [9, -5i, -1, 5i]\)
现在计算 \(FFT([9, -5i, -1, 5i])\)，再共轭，再除以 \(4\)。
执行 \(FFT([9, -5i, -1, 5i])\)
Level 0（\(n=4\)）
偶部：\([9, -1]\)
奇部：\([-5i, 5i]\)
Level 1（\(n=2\)）
\(FFT([9, -1]) = [9 + (-1), 9 - (-1)] = [8, 10]\)
\(FFT([-5i, 5i]) = [-5i + 5i, -5i - 5i] = [0, -10i]\)
合并（\(\omega_4 = i\)）
\(k=0:\)
\(Z[0] = 8 + 1·0 = 8\)
\(Z[2] = 8 - 0 = 8\)
\(k=1:\)
\(Z[1] = 10 + i·(-10i) = 10 + (-10i^2) = 10 + 10 = 20\)
\(Z[3] = 10 - i·(-10i) = 10 - 10 = 0\)
所以：
\(\text{FFT}(\overline{\hat{c}}) = [8, 20, 8, 0]\)
取共轭（实数，不变）→ \([8, 20, 8, 0]\)
除以 N=4：
\(c = [2, 5, 2, 0]\)
✅ 最终系数：\([2, 5, 2, 0]\) → 对应多项式 \(2 + 5x + 2x^2\)
完美匹配直接乘法结果！
📌 总结流程（带分治细节）
补零至 \(N=4\)
\(FFT(a):\)
拆：\([1,2,0,0]\) → 偶=\([1,0]\), 奇=\([2,0]\)
递归：\(FFT([1,0])=[1,1], FFT([2,0])=[2,2]\)
合并：用 \(\omega_4^k\) 得 \([3, 1+2i, -1, 1−2i]\)
\(FFT(b):\) 类似，得 \([3, 2+i, 1, 2−i]\)
点乘：得 \([9, 5i, −1, −5i]\)
\(IFFT：\)通过共轭+FFT+缩放，得 \([2,5,2,0]\)
浙公网安备 33010602011771号