在O(nlogn)
时间内将系数表示多项式转化为点值表示,
n
n
n个点确定一条
n
−
1
n-1
n−1次的多项式曲线,后者可以O(n)
时间内完成高精度乘法。
而这 n n n个点的取法是有讲究的,均匀取自复平面的单位圆,对,将 x x x的取值扩展到复数域,这些点也称为单位根,定义为 ω n k = c o s ( k n 2 π ) + i s i n ( k n 2 π ) \omega_n^k=cos(\frac{k}{n}2\pi)+isin(\frac{k}{n}2\pi) ωnk=cos(nk2π)+isin(nk2π),默认 n n n为2的整数次幂。
直接将这些特殊的点作为x
代入计算多项式的时间复杂度是
O
(
n
2
)
O(n^2)
O(n2)的,而本文主角fft
通过分治的技巧将复杂度降为
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)。
对于多项式
A
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
.
.
.
a
n
−
1
x
n
−
1
A(x)=a_0+a_1x+a_2x^2+...a_{n-1}x^{n-1}
A(x)=a0+a1x+a2x2+...an−1xn−1,
将项根据下标的奇偶性分别为两部分,
A
(
x
)
=
(
a
0
+
a
2
x
2
+
.
.
.
+
a
n
−
2
x
n
−
2
)
+
x
(
a
1
+
a
3
x
2
+
.
.
.
+
a
n
−
1
x
n
−
2
)
A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})
A(x)=(a0+a2x2+...+an−2xn−2)+x(a1+a3x2+...+an−1xn−2),
设
A
1
(
x
)
=
a
0
+
a
2
x
+
a
4
x
2
+
.
.
.
+
a
n
−
2
x
n
2
−
1
A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1}
A1(x)=a0+a2x+a4x2+...+an−2x2n−1,
A
2
(
x
)
=
a
1
+
a
3
x
+
a
5
x
2
+
.
.
.
+
a
n
−
1
x
n
2
−
1
A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1}
A2(x)=a1+a3x+a5x2+...+an−1x2n−1,
则可以写成
A
(
x
)
=
A
1
(
x
2
)
+
x
A
2
(
x
2
)
A(x)=A_1(x^2)+xA_2(x^2)
A(x)=A1(x2)+xA2(x2),
设
k
<
n
2
k<\frac{n}{2}
k<2n,将
ω
n
k
\omega_n^k
ωnk作为
x
x
x代入
A
(
x
)
A(x)
A(x),利用
单位根的性质,得
A
(
ω
n
k
)
A(\omega_n^k)
A(ωnk)
=
A
1
(
(
ω
n
k
)
2
)
+
ω
n
k
A
2
(
(
ω
n
k
)
2
)
=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2)
=A1((ωnk)2)+ωnkA2((ωnk)2)
=
A
1
(
ω
n
2
k
)
+
ω
n
k
A
2
(
ω
n
2
k
)
=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})
=A1(ωn2k)+ωnkA2(ωn2k)
=
A
1
(
ω
n
2
k
)
+
ω
n
k
A
2
(
ω
n
2
k
)
=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_{\frac{n}{2}}^{k})
=A1(ω2nk)+ωnkA2(ω2nk),
同理代入
ω
n
k
+
n
2
\omega_n^{k+\frac{n}{2}}
ωnk+2n得到
A
(
ω
n
k
+
n
2
)
=
A
1
(
ω
n
2
k
)
−
ω
n
k
A
2
(
ω
n
2
k
)
A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_2(\omega_{\frac{n}{2}}^{k})
A(ωnk+2n)=A1(ω2nk)−ωnkA2(ω2nk),
以上两个多项式后一项仅有符号不同,因此,已知 A 1 ( ω n 2 k ) A_1(\omega_{\frac{n}{2}}^k) A1(ω2nk)和 A 2 ( ω n 2 k ) A_2(\omega_{\frac{n}{2}}^k) A2(ω2nk)可以推出 A ( ω n k ) A(\omega_n^k) A(ωnk)和 A ( ω n k + n 2 ) A(\omega_n^{k+\frac{n}{2}}) A(ωnk+2n)。进一步,可以通过递归分治实现fft,,将序列根据奇偶性分为两部分分别处理,回溯时通过上述两个关系求出序列。
代码模板
struct Complex
{
double x, y;
Complex(double x = 0, double y = 0) : x(x), y(y) {}
} a[N], b[N];
Complex operator + (Complex a, Complex b) {
return Complex(a.x + b.x, a.y + b.y);
}
Complex operator - (Complex a, Complex b) {
return Complex(a.x - b.x, a.y - b.y);
}
Complex operator * (Complex a, Complex b) {
return Complex(a.x * b.x - a.y * b.y, a.x * b.y + b.x * a.y);
}
void fft(Complex *a, int n, int inv) {
if (n == 1) return;
int mid = n / 2;
static Complex tmp[N];
for (int i = 0; i < mid; i++) {
tmp[i] = a[i * 2];
tmp[i + mid] = a[i * 2 + 1];
}
for (int i = 0; i < n; i++) a[i] = tmp[i];
fft(a, mid, inv), fft(a + mid, mid, inv);
// 默认n为2的次幂
for(int i = 0; i < mid; i++) {
Complex x(cos(2 * pi * i / n), inv * sin(2 * pi * i / n));
tmp[i] = a[i] + x * a[i + mid];
tmp[i + mid] = a[i] - x * a[i + mid];
}
for (int i = 0; i < n; i++) a[i] = tmp[i];
}
void solve() {
int n, m; // 这里是次数
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; i++) { scanf("%lf", &a[i].x); }
for (int i = 0; i <= m; i++) { scanf("%lf", &b[i].x); }
int k = 1;
while (k <= 2 * max(n, m)) { k *= 2; }
// printf("k=%d\n", k);
// 将n调整为2的整数次幂(用k表达)
fft(a, k, 1), fft(b, k, 1);
for (int i = 0; i < k; i++) a[i] = a[i] * b[i];
fft(a, k, -1);
for (int i = 0; i <= (n + m); i++) { // 最终是(n+m+1)项
printf("%d ", (int)(a[i].x / k + 0.5));
} printf("\n");
}
坑,循环优化写法。
reference
十分简明易懂的FFT(快速傅里叶变换)
FFT模板-附详细注释