机房的人都嘲笑我博客只有两句话,我今天就要打他们的脸!!!我今天写一下FFT这个神奇的算法.首先,我们需要知道,FFT好像就是用来处理两个多项式乘积的,没有多么高端.首先思考暴力,就是暴力把n个数带入求值,然后相乘,复杂度O(n^2).我们考虑分治,分治奇偶,把奇数提出来,再进行分治.先放一个递归的简单易懂的代码:
#include<iostream> #include<cstdio> #include<cmath> #include<queue> #include<algorithm> #include<vector> #include <complex> #include<cstring> using namespace std; #define duke(i,a,n) for(int i = a;i <= n;i++) #define lv(i,a,n) for(int i = a;i >= n;i--) #define clean(a) memset(a,0,sizeof(a)) #define mp make_pair #define cp complex<db> #define enter puts("") const int INF = 1 << 30; const double eps = 1e-8; typedef long long ll; typedef double db; template <class T> void read(T &x) { char c; bool op = 0; while(c = getchar(), c < '0' || c > '9') if(c == '-') op = 1; x = c - '0'; while(c = getchar(), c >= '0' && c <= '9') x = x * 10 + c - '0'; if(op) x = -x; } template <class T> void write(T x) { if(x < 0) putchar('-'), x = -x; if(x >= 10) write(x / 10); putchar('0' + x % 10); } const int N = 4000005; const db PI = acos(-1); cp a[N],b[N],omg[N],inv[N]; int len = 1,ans[N]; bool inversed = false; cp omega(int n,int k) { if(inversed == true) return cp(cos(2 * PI * k / n),sin(2 * PI * k / n)); else return conj(cp(cos(2 * PI * k / n),sin(2 * PI * k / n))); } void fft(cp *a,int n) { if(n == 1) return; static cp buf[N]; int m = n / 2; duke(i,0,m - 1) { buf[i] = a[2 * i]; buf[i + m] = a[2 * i + 1]; } for(int i = 0;i < n;i++) { a[i] = buf[i]; } fft(a,m); fft(a + m,m); duke(i,0,m - 1) { cp x = omega(n,i); buf[i] = a[i] + x * a[i + m]; buf[i + m] = a[i] - x * a[i + m]; } duke(i,0,n - 1) a[i] = buf[i]; } int n,m; int main() { read(n);read(m); while(len < n + m + 1) len <<= 1; // cout<<len<<endl; int t; duke(i,0,n) { read(t); a[i].real(t); } duke(i,0,m) { read(t); b[i].real(t); } fft(a,len);fft(b,len); duke(i,0,len) { a[i] *= b[i]; } inversed = true; fft(a,len); duke(i,0,len) { ans[i] = floor(a[i].real() / len + 0.5); } duke(i,0,n + m) { printf("%d ",ans[i]); } return 0; }
这个代码简单易懂,但是luogu会T2个点,怎么办呢,我们尝试2进制优化,每次直接处理出现在这位的数,得到非递归版本的FFT,上一下代码:
#include<iostream> #include<cstdio> #include<cmath> #include<queue> #include<algorithm> #include<vector> #include <complex> #include<cstring> using namespace std; #define duke(i,a,n) for(int i = a;i <= n;i++) #define lv(i,a,n) for(int i = a;i >= n;i--) #define clean(a) memset(a,0,sizeof(a)) #define mp make_pair #define cp complex<db> #define enter puts("") const int INF = 1 << 30; const double eps = 1e-8; typedef long long ll; typedef double db; template <class T> void read(T &x) { char c; bool op = 0; while(c = getchar(), c < '0' || c > '9') if(c == '-') op = 1; x = c - '0'; while(c = getchar(), c >= '0' && c <= '9') x = x * 10 + c - '0'; if(op) x = -x; } template <class T> void write(T x) { if(x < 0) putchar('-'), x = -x; if(x >= 10) write(x / 10); putchar('0' + x % 10); } const int N = 4000005; const db PI = acos(-1); cp a[N],b[N],omg[N],inv[N]; int len = 1,ans[N]; void init(int n) { duke(i,0,n - 1) { omg[i] = cp(cos(2 * PI * i / n),sin(2 * PI * i / n)); inv[i] = conj(omg[i]); } } void fft(cp *a,int n) { int lim = 0; while((1 << lim) < n) lim++; duke(i,0,n - 1) { int t = 0; duke(j,0,lim - 1) { if((i >> j) & 1) t |= (1 << (lim - j - 1)); } if(i < t) swap(a[i],a[t]); } static cp buf[N]; for(int l = 2;l <= n;l *= 2) { int m = l / 2; for(int j = 0;j < n;j += l) { for(int i = 0;i < m;i++) { // cp x = omega(n / l,i); buf[j + i] = a[j + i] + omg[n / l * i] * a[j + i + m]; buf[j + i + m] = a[j + i] - omg[n / l * i] * a[j + i + m]; } } duke(i,0,n - 1) a[i] = buf[i]; } } int n,m; int main() { read(n);read(m); while(len < n + m + 1) len <<= 1; // cout<<len<<endl; int t; duke(i,0,n) { read(t); a[i].real(t); } duke(i,0,m) { read(t); b[i].real(t); } init(len); fft(a,len);fft(b,len); duke(i,0,len) { a[i] *= b[i]; } duke(i,0,len) { omg[i] = inv[i]; } fft(a,len); duke(i,0,len) { ans[i] = floor(a[i].real() / len + 0.5); } duke(i,0,n + m) { printf("%d ",ans[i]); } return 0; }
但还是很慢,还是会T两个点,我们需要用一个叫蝴蝶优化的东西,就是把操作顺序改变一下就行了.具体其实我不太懂,先背板子以后慢慢理解吧:
// luogu-judger-enable-o2 #include<iostream> #include<cstdio> #include<cmath> #include<queue> #include<algorithm> #include<vector> #include <complex> #include<cstring> using namespace std; #define duke(i,a,n) for(int i = a;i <= n;i++) #define lv(i,a,n) for(int i = a;i >= n;i--) #define clean(a) memset(a,0,sizeof(a)) #define mp make_pair #define cp complex<db> #define enter puts("") const int INF = 1 << 30; const double eps = 1e-8; typedef long long ll; typedef double db; template <class T> void read(T &x) { char c; bool op = 0; while(c = getchar(), c < '0' || c > '9') if(c == '-') op = 1; x = c - '0'; while(c = getchar(), c >= '0' && c <= '9') x = x * 10 + c - '0'; if(op) x = -x; } template <class T> void write(T x) { if(x < 0) putchar('-'), x = -x; if(x >= 10) write(x / 10); putchar('0' + x % 10); } const int N = 4000005; const db PI = acos(-1); cp a[N],b[N],omg[N],inv[N]; int len = 1,ans[N]; void init(int n) { duke(i,0,n - 1) { omg[i] = cp(cos(2 * PI * i / n),sin(2 * PI * i / n)); inv[i] = conj(omg[i]); } } void fft(cp *a,int n) { int lim = 0; while((1 << lim) < n) lim++; duke(i,0,n - 1) { int t = 0; duke(j,0,lim - 1) { if((i >> j) & 1) t |= (1 << (lim - j - 1)); } if(i < t) swap(a[i],a[t]); } for(int l = 2;l <= n;l *= 2) { int m = l / 2; for(cp *p = a;p != a + n;p += l) { for(int i = 0;i < m;i++) { cp t = omg[n / l * i] * p[i + m]; p[i + m] = p[i] - t; p[i] += t; } } } } int n,m; int main() { read(n);read(m); while(len < n + m + 1) len <<= 1; // cout<<len<<endl; int t; duke(i,0,n) { read(t); a[i].real(t); } duke(i,0,m) { read(t); b[i].real(t); } init(len); fft(a,len);fft(b,len); duke(i,0,len) { a[i] *= b[i]; } duke(i,0,len) { omg[i] = inv[i]; } fft(a,len); duke(i,0,len) { ans[i] = floor(a[i].real() / len + 0.5); } duke(i,0,n + m) { printf("%d ",ans[i]); } return 0; }