也就对应的是圆上的几个点。
所以我们可以得到一个结论。多项式\(A(x)\)的离散傅里叶变换的另一个多项式\(B(x)\)的系数,取单位根的倒数\(w_n^0,w_n^{-1},...,x_n^{-(n-1)}\)作为\(x\)代入\(B(x)\),得到的每个数再除以\(n\),得到的是\(A(x)\)的各项系数。这就实现了傅里叶逆变换。
至此\(FFT\)的理论基础已结束。
根据上述分析可得,一个序列需要划分成两部分后分治递归即可。
但是可以再度优化。
可以发现原序列和后序列分组其实是按照原序列下标的二进制翻转。
因此按照下标进行奇偶性分类是没有必要的,于是我们可以免去递归的过程。
对于二进制翻转要怎么做呢?
这是一个\(trick:\)蝴蝶定理。
假设即将反转的数字为\(i\),在\(i\)之前的数字都已经翻转好了。
那么对于\(i/2\),就相当于\(i\)右移了一位后翻转完成的对应的那个数右移一位,如果被移走的那位是\(1\)的话,在最高位补上\(1\)。(想不通就看代码后手摸几个例子)
之后直接枚举子区间后向上合并即可。
枚举分割的中点的时间复杂度为\(O(logn)\),合并复杂度为\(O(n)\),总时间复杂度为\(O(nlogn)\)。
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e7 + 10;
const double PI = acos(-1.0);
int n, m, limit, l, r[maxn];
inline int read() {
char c = getchar(); int x = 0, f = 1;
while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();}
while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();}
return x * f;
}
//手写复数类
struct Complex
{
double x, y;
Complex (double xx=0, double yy=0){
x = xx; y = yy;
}
Complex operator + (const Complex b) const{
return Complex(x+b.x, y+b.y);
}
Complex operator - (const Complex b) const{
return Complex(x-b.x, y-b.y);
}
Complex operator * (const Complex b) const{
return Complex(x*b.x-y*b.y, x*b.y+y*b.x);
}
}a[maxn], b[maxn];
void fft(Complex c[], int type)
{
for(int i = 0; i < limit; i++)
if(i < r[i]) swap(c[i], c[r[i]]);
//枚举待合并区间的中点的长度
for(int mid = 1; mid < limit; mid <<= 1)
{
//设立单位根
Complex wn(cos(PI/mid), type*sin(PI/mid));
//R是区间的长度,j表示当前已经到哪个位置了,而且是左端点
for(int R = mid<<1, j = 0; j < limit; j += R)
{
Complex w(1, 0); //初始幂次
for(int k = 0; k < mid; k++, w = w*wn) //枚举左半部分
{
Complex x = c[j+k], y = w*c[j+mid+k];
c[j+k] = x+y;
c[j+mid+k] = x-y; //右半部分减去即可
}
}
}
}
int main()
{
//1e6的读入,需要写快读
n = read(), m = read();
for(int i = 0; i <= n; i++) a[i].x = read();
for(int i = 0; i <= m; i++) b[i].x = read();
limit = 1;
//要求limit一定是2的整次幂
while(limit <= n+m) limit <<= 1, l++;
for(int i = 0; i < limit; i++)
r[i] = (r[i>>1]>>1)|((i&1)<<(l-1));
//对a序列和b序列分别处理
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 ", (int)(a[i].x/limit+0.5));
return 0;
}
原文:https://www.cnblogs.com/zxytxdy/p/12078627.html