FFT(快速傅里叶变换)
学习了FFT用来求多项式的乘法，看了算导上的介绍，上面讲的非常明白，概括一下FFT的原理就是，我们在计算多项式的乘法时，如果暴力模拟的话是n^2 复杂度的，就像小学学的竖式乘法一样，比如一个n位数乘上一个n位数，我们需要用竖式乘法计算要列n行，每一行有n个元素，然后相邻两行错开一位（很显然，竖式乘法谁都会），如果我们换一种形式呢？有一种表达是叫做点值表达，就好像我们上了初中学二次函数，如果已知函数图像上的三个不同的点坐标，我们可以写出函数的表达式，那么就是说利用函数图象上的点我们也可以确定一个函数（多项式），这样的好处在于点值表达式更加容易进行乘法运算，举个例子：
我们现在有两个函数： y1=2*x*x+5*x+1 y2=x*x+2*x+1 我们要求y1*y2的表达式。如果利用系数表达，我们利用竖式乘法可以手算出答案，复杂度n^2.
由于两个函数是二次的，乘起来是四次的，所以我们需要知道5个点的坐标。
| x | 1 | 2 | 3 | 4 | 5 |
| y1 | 8 | 19 | 34 | 53 | 76 |
| y2 | 4 | 9 | 16 | 25 | 36 |
| y | 32 | 171 | 544 | 1325 | 2736 |
我们现在已知了y的五个坐标，那么肯定可以求出y的表达式，那么多项式乘法也就做完了。
可以看出对于点值表达，它的乘法操作是O(N)的，只需要把对应的y值乘起来即可。
那么FFT算法就诞生了，他就是利用点值表达的优点，把系数表达转换成点值表达，然后再转化回去。
从系数表达到点值表达的过程交做求值，利用秦九韶公式我们可以O(N)的算出一个点对应的y值。
从点值表达到系数表达的过程叫做差值，著名的拉格朗日插值法复杂度是O(N^2)的，这对于我们多项式乘法没有太大帮助。
因为差值对于我们取哪些点并没有要求，所以我们选取一些特殊点，利用一些特殊性质可以优化插值的过程。
利用了复数根的知识，大概原理就是，我们可以分治的求一个系数向量，每一次递归求一半，然后可以O(N)的合并，在算导的第三十章都有讲到，推荐大家去看算导，里面讲的非常详细易懂。由于我太菜了，讲不清其中的奥秘。
所以只能在这里给出程序代码了，我没有进行递归的求解，进行了迭代操作，看一下算导上30-4的图就明白r数组的含义了，其实就是计算的顺序。 —— by VANE
- #include<bits/stdc++.h>
- #define N 3000005
- #define pi acos(-1)
- using namespace std;
- typedef complex<double> E;
- int n,m,l,r[N];
- E a[N],b[N];
- void FFT(E *a,int f){
- for(int i=0;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
- for(int i=1;i<n;i<<=1)
- {
- E wn(cos(pi/i),f*sin(pi/i));
- for(int p=i<<1,j=0;j<n;j+=p)
- {
- E w(1,0);
- for(int k=0;k<i;k++,w*=wn)
- {
- E x=a[j+k],y=w*a[j+k+i];
- a[j+k]=x+y;a[j+k+i]=x-y;
- }
- }
- }
- }
- inline int read()
- {
- char ch=getchar();int f=1,x=0;
- while(ch<\'0\'||ch>\'9\') {if(ch==\'-\') f=-f;ch=getchar();}
- while(ch>=\'0\'&&ch<=\'9\') x=x*10+ch-\'0\',ch=getchar();
- return x*f;
- }
- 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();
- m+=n;for(n=1;n<=m;n<<=1)l++;
- for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
- FFT(a,1);FFT(b,1);
- for(int i=0;i<=n;i++)a[i]=a[i]*b[i];
- FFT(a,-1);
- for(int i=0;i<=m;i++)printf("%d ",(int)(a[i].real()/n+0.5));
- }
附上一个高精度乘法，其实就是FFT的简单变形。
- // luogu-judger-enable-o2
- #include<bits/stdc++.h>
- #define N 300005
- #define pi acos(-1)
- using namespace std;
- typedef complex<double> E;
- int n,m,l,r[N],ans[N];
- E a[N],b[N];
- char s[N];
- void FFT(E *a,int f){
- for(int i=0;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
- for(int i=1;i<n;i<<=1)
- {
- E wn(cos(pi/i),f*sin(pi/i));
- for(int p=i<<1,j=0;j<n;j+=p)
- {
- E w(1,0);
- for(int k=0;k<i;k++,w*=wn)
- {
- E x=a[j+k],y=w*a[j+k+i];
- a[j+k]=x+y;a[j+k+i]=x-y;
- }
- }
- }
- }
- inline int read()
- {
- char ch=getchar();int f=1,x=0;
- while(ch<\'0\'||ch>\'9\') {if(ch==\'-\') f=-f;ch=getchar();}
- while(ch>=\'0\'&&ch<=\'9\') x=x*10+ch-\'0\',ch=getchar();
- return x*f;
- }
- int main(){
- n=read();n--;m=n;
- scanf("%s",s);
- for(int i=0;i<=n;i++)a[i]=s[n-i]-\'0\';
- scanf("%s",s);
- for(int i=0;i<=m;i++)b[i]=s[m-i]-\'0\';
- m+=n;for(n=1;n<=m;n<<=1)l++;
- for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
- FFT(a,1);FFT(b,1);
- for(int i=0;i<=n;i++)a[i]=a[i]*b[i];
- FFT(a,-1);
- for(int i=0;i<=m;i++)ans[i]=(int)(a[i].real()/n+0.5);
- for(int i=0;i<=m;++i)
- ans[i+1]+=ans[i]/10,ans[i]%=10;
- if(ans[m+1]) m++;
- while(ans[m]==0&&m) m--;
- for(int i=m;i>=0;--i) printf("%d",ans[i]);
- }