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]);
}