前言:写得比较简陋~以后有时间一定来完善。
快速傅里叶变换,又叫FFT。
可以在O(n)内处理多项式乘法的算法。
一个n - 1多项式g(x) 可以用系数表示为$sum_{i = 0}^{n - 1} ai * x^{i}$
对于我们g(x)多项式中的每一项,我们可以看成有一个函数f(x)。
将每一项的x,x^1 ...都可以得到多项式中的对应的每一项。
那么每一项就可以看成一个点,在f(x)函数的图像上。
那么这n个点就可以确定唯一的一个f(x)函数。
这就是转化为点值法。
假设两个多项式点值法表示为:
F(x) = ((x0,f(x0)),(x1,f(x1)...)
G(x) = ((x0,g(x0)),(x1,g(x1)...)
那么他们的卷积H[x] = (x0,f(x0) * g(x0),....)
所以如果我们知道两个多项式的点值表示,那么就可以O(n)求出他们的卷积的点值表示。
然后引入单位根wn和欧拉定理:这里不赘述,可以自习百度。
然后根据推演:
对于wn^k 的k < n / 2的情况 = A + wn^k *B
对于wn^k 的k > n / 2的情况 = A - wn^k * B
两者只有一个常数系数的不同。
所以我们在算上面一项的时候可以O(1)算得下面一项的值。
所以我们可以折半分治它,借由蝴蝶操作。
最后我们算得两个多项式卷积的点值表示后,再有推论:ak = ck / n。
可得我们的系数表示 = 点值表示 / n。
#include <bits/stdc++.h> using namespace std; typedef long long LL; typedef pair<LL, int> pii; const int N = 2e6 + 5; const int M = 2000 + 5; const LL Mod = 998244353; #define pi acos(-1) #define INF 1e9 #define dbg(ax) cout << "now this num is " << ax << endl; namespace FASTIO { inline LL read() { LL x = 0, f = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') f = -1; c = getchar(); } while (c >= '0' && c <= '9') { x = (x << 1) + (x << 3) + (c ^ 48); c = getchar(); } return x * f; } } using namespace FASTIO; complex<double> a[N],b[N]; int limit = 1,n,m; void FFT(int len,complex<double> *a,int id) { if(len == 1) return ; complex<double> a1[len >> 1],a2[len >> 1]; for(int i = 0;i < len;i += 2) {//根据奇偶分类 a1[i >> 1] = a[i],a2[i >> 1] = a[i + 1]; } FFT(len >> 1,a1,id); FFT(len >> 1,a2,id); complex<double> wn(cos(2.0 * pi / len),id * sin(2.0 * pi / len)),w(1,0); // for(int i = 0;i < (len >> 1);++i,w *= wn) { a[i] = a1[i] + w * a2[i]; a[i + (len >> 1)] = a1[i] - w * a2[i]; } } int main() { n = read(),m = read(); for(int i = 0;i <= n;++i) a[i] = read(); for(int i = 0;i <= m;++i) b[i] = read(); 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) printf("%d%c",(int)(a[i].real() / limit + 0.5),i == n + m ? ' ' : ' '); //system("pause"); return 0; }
递归法的缺点很明显,因为要进行奇偶分类,所以申请内存会浪费很多时间。
由此,我们有了迭代法:对分治层操作的序列进行分析后,我们wn^k在最终操作后的序列中的位置就是2进制的反转。
这样我们先处理出反转位置就可以快速优化FFT。
#include <bits/stdc++.h> using namespace std; typedef long long LL; typedef pair<LL, int> pii; const int N = 4e6 + 5; const int M = 2000 + 5; const LL Mod = 998244353; #define pi acos(-1) #define INF 1e9 #define dbg(ax) cout << "now this num is " << ax << endl; namespace FASTIO { inline LL read() { LL x = 0, f = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') f = -1; c = getchar(); } while (c >= '0' && c <= '9') { x = (x << 1) + (x << 3) + (c ^ 48); c = getchar(); } return x * f; } } using namespace FASTIO; int limit = 1,n,m,r[N]; complex<double> a[N],b[N]; void FFT(complex<double> *a,int id) { for(int i = 0;i < limit;++i) { if(i < r[i]) swap(a[i],a[r[i]]); } for(int i = 1;i < limit;i <<= 1) { complex<double> x(cos(pi / i),id * sin(pi / i));//这里变化为了pi / i for(int j = 0;j < limit;j += (i << 1)) { complex<double> y(1,0); for(int k = 0;k < i;++k,y *= x) { complex<double> p = a[j + k],q = y * a[j + k + i]; a[j + k] = p + q; a[j + k + i] = p - q; } } } } int main() { n = read(),m = read(); for(int i = 0;i <= n;++i) a[i] = read(); for(int i = 0;i <= m;++i) b[i] = read(); int L = 0; while(limit <= (n + m)) limit <<= 1,++L; for(int i = 0;i < limit;++i) { r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1)); } FFT(a,1); FFT(b,1); for(int i = 0;i <= limit;++i) a[i] = a[i] * b[i]; FFT(a,-1); for(int i = 0;i <= n + m;++i) printf("%d%c",(int)(a[i].real() / limit + 0.5),i == n + m ? ' ' : ' '); // system("pause"); return 0; }