目的
多项式 F ( x ) = ∑ i = 0 n f i x i F(x)=\sum_{i=0}^n f_ix^i F(x)=∑i=0nfixi , G ( x ) = ∑ i = 0 m g i x i ,G(x)=\sum_{i=0}^m g_ix^i ,G(x)=∑i=0mgixi若要实现乘法,朴素情况下显然可以发现是 O ( n 2 ) O(n^2) O(n2)级别的。
然而我们知道,一个次数为n的多项式,又可以被n个点值构成的集合表示。
通过巧妙选取点值(FFT),将多项式乘法转化为点值乘法,就可以将复杂度压缩至 O ( n l o g n ) O(nlogn) O(nlogn)级别。
简明理论背景
ω n 1 \omega ^1 _n ωn1称为n次单位根,将圆n等分,第k个点就记为 ω n k \omega ^k _n ωnk
ω n 1 = c o s 1 n 2 π + i s i n 1 n 2 π \omega ^1 _n=cos\frac{1}{n}2\pi+isin\frac{1}{n}2\pi ωn1=cosn12π+isinn12π
ω n k = c o s k n 2 π + i s i n k n 2 π = ( ω n 1 ) k \omega ^k _n=cos\frac{k}{n}2\pi+isin\frac{k}{n}2\pi = (\omega^1_n)^k ωnk=cosnk2π+isinnk2π=(ωn1)k
还要用到三个前置性质
ω n k = ω 2 n 2 k \omega^k_n=\omega^{2k}_{2n}\\ ωnk=ω2n2k
ω n k + n 2 = − ω n k ω n k + n = ω n k \omega^{k+\frac{n}{2}}_n=-\omega^k_n\\ \omega^{k+n}_n=\omega^k_n ωnk+2n=−ωnkωnk+n=ωnk
代入上式即可证,第二个式子可以理解为表示的向量等大反向(n/2即圆n等分的一半,即向量转了半圈),第三个就是转了一圈
FFT分治实现
F ( x ) = ∑ i = 0 n − 1 a i x i = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 = ( 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 ) F(x)=\sum_{i=0}^{n-1} a_ix^i=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} \\=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2}) F(x)=i=0∑n−1aixi=a0+a1x+a2x2+...+an−1xn−1=(a0+a2x2+...+an−2xn−2)+x(a1+a3x2+...+an−1xn−2)
令 F L ( x ) = ∑ i = 0 n 2 − 1 a 2 i x i , F R ( x ) = ∑ i = 0 n 2 − 1 a 2 i + 1 x i 则 F ( x ) = F L ( x 2 ) + x F R ( x 2 ) 令\displaystyle F_L(x)=\sum_{i=0}^{{n\over 2}-1}a_{2i}x^i, F_R(x)=\sum_{i=0}^{{n\over 2}-1}a_{2i+1}x^i\\则F(x)=F_L(x^2)+xF_R(x^2) 令FL(x)=i=0∑2n−1a2ixi,FR(x)=i=0∑2n−1a2i+1xi则F(x)=FL(x2)+xFR(x2)
利用单位根的性质进行代入
F ( ω n i ) = F L ( ω n 2 i ) + ω n i F R ( ω n 2 i ) = F L ( ω n 2 i ) + ω n i F R ( ω n 2 i ) , 0 ≤ i ≤ n 2 − 1 , i ∈ Z F ( ω n i + n 2 ) = F L ( ω n 2 i ) − ω n i F R ( ω n 2 i ) , n 2 ≤ i ≤ n − 1 , i ∈ Z \displaystyle F(\omega_n^i)=F_L(\omega_n^{2i})+\omega_n^iF_R(\omega_n^{2i})=F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i),0≤i≤ \frac{n}{2} −1,i∈Z \\ \displaystyle F(\omega_n^{i+{n\over 2}})=F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i),\frac{n}{2} ≤i≤n-1 ,i∈Z F(ωni)=FL(ωn2i)+ωniFR(ωn2i)=FL(ω2ni)+ωniFR(ω2ni),0≤i≤2n−1,i∈ZF(ωni+2n)=FL(ω2ni)−ωniFR(ω2ni),2n≤i≤n−1,i∈Z
发现 F L ( ω n 2 i ) , F R ( ω n 2 i ) F_L(\omega_{n\over 2}^i), F_R(\omega_{n\over 2}^i) FL(ω2ni),FR(ω2ni)只有前面的系数有所差异。即我们只要得到前n/2的点就可以推后n/2个点(知道前半圆就可以知道一整个圆)
这又显然符合分治结构。而求解 F L ( ω n 2 i ) , F R ( ω n 2 i ) F_L(\omega_{n\over 2}^i), F_R(\omega_{n\over 2}^i) FL(ω2ni),FR(ω2ni)后的合并是 O ( n ) O(n) O(n)的。则有 T ( n ) = 2 T ( n 2 ) + Θ ( n ) = Θ ( n log n ) \displaystyle T(n)=2T({n\over 2})+\Theta(n)=\Theta(n\log n) T(n)=2T(2n)+Θ(n)=Θ(nlogn)
递归实现
参考了http://www.longluo.me/blog/2022/04/02/FFT
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x, y = _y;
}
Complex operator -(const Complex &a) const {
return Complex(x - a.x, y - a.y);
}
Complex operator +(const Complex &a) const {
return Complex(x + a.x, y + a.y);
}
Complex operator *(const Complex &a) const {
return Complex(x * a.x - y * a.y, x * a.y + y * a.x);
}
}; // 封装了复数
const double PI = acos(-1.0);
vector<Complex> FFT(vector<Complex> &a) {
// 实部为系数
int n = a.size();
if (n == 1) return a;
vector<Complex> FL(n / 2), FR(n / 2);
for (int i = 0; 2 * i < n; i++) { // 拆分奇偶
FL[i] = a[2 * i];
FR[i] = a[2 * i + 1];
}
vector<Complex> GL = FFT(FL), GR = FFT(FR); // 分治
vector<Complex> G(n);
double angle = 2 * PI / n; // omega^1_n
Complex omega = Complex(cos(angle), sin(angle));
Complex curRoot = Complex(1, 0);
for (int i = 0; 2 * i < n; i++) { // 上述递推式
G[i] = GL[i] + curRoot * GR[i];
G[i + n / 2] = GL[i] - curRoot * GR[i];
curRoot = curRoot * omega;
}
return G;
}
实例应用
大数乘法https://leetcode.cn/problems/multiply-strings/
#include <iostream>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x, y = _y;
}
Complex operator -(const Complex &a) const {
return Complex(x - a.x, y - a.y);
}
Complex operator +(const Complex &a) const {
return Complex(x + a.x, y + a.y);
}
Complex operator *(const Complex &a) const {
return Complex(x * a.x - y * a.y, x * a.y + y * a.x);
}
};
const double PI = acos(-1.0);
vector<Complex> FFT(vector<Complex> &a, bool flag) {
// 实部为系数
int n = a.size();
if (n == 1) return a;
vector<Complex> FL(n / 2), FR(n / 2);
for (int i = 0; 2 * i < n; i++) {
FL[i] = a[2 * i];
FR[i] = a[2 * i + 1];
}
vector<Complex> GL = FFT(FL, flag), GR = FFT(FR, flag);
vector<Complex> G(n);
double angle = 2 * PI / n * (flag ? -1 : 1);
Complex omega(cos(angle), sin(angle));
Complex curRoot(1, 0);
for (int i = 0; 2 * i < n; i++) {
G[i] = GL[i] + curRoot * GR[i];
G[i + n / 2] = GL[i] - curRoot * GR[i];
curRoot = curRoot * omega;
}
return G;
}
string multiply(string num1, string num2) {
if (num1 == "0" || num2 == "0") return "0";
int len1 = num1.size(), len2 = num2.size();
int n = 1;
while (n < len1 + len2) n <<= 1;
vector<Complex> a(n), b(n);
for (int i = len1 - 1; i >= 0; i--) {
a[i] = Complex(num1[len1 - 1 - i] - '0', 0);
}
for (int i = len2 - 1; i >= 0; i--) {
b[i] = Complex(num2[len2 - 1 - i] - '0', 0);
}
a = FFT(a, 0), b = FFT(b, 0);
for (int i = 0; i < n; i++)
a[i] = a[i] * b[i];
a = FFT(a, 1);
string ans;
int carry = 0;
for (int i = 0; i < n; i++) {
int sum = round(round(a[i].x) / n) + carry;
carry = sum / 10;
ans += sum % 10 + '0';
}
if (carry > 0) {
ans += carry % 10 + '0';
}
int ind = ans.size() - 1;
while (ans[ind] == '0' && ind > 0) {
ind --;
}
ans = ans.substr(0, ind + 1);
reverse(ans.begin(), ans.end());
return ans;
}
int main() {
string a, b;
cin >> a >> b;
cout << multiply(a, b) << endl;
return 0;
}
参考
- http://www.longluo.me/blog/2022/04/02/FFT
- https://www.cnblogs.com/JustinRochester/p/17036805.html
- https://blog.csdn.net/enjoy_pascal/article/details/81478582/