FFT学习笔记

前置知识你可以装作不会

复数相乘:

\((a+bi)*(c+di)=(ac-bd)+(bc+ad)i​\)

欧拉定理:
\(e^{ix}=cosx+isinx​\)

单位根:

\(\omega_n^k=e^{\frac{2\pi ki}{n}}=cos\frac{2k\pi}{n}+sin\frac{2k\pi}{n}i\)

\(\omega_n^k表示x^n=1在复数域上的第k个解,其中\omega_n^0=\omega_n^n=1\)

神奇性质:

\(1.\omega_n^k=\omega_{2n}^{2k}\)

\(2.\sum_{i=0}^{n-1}(\omega_n^k)^i=[n|k]​\)

\(3.\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)

那么现在我们有一个\(f(x)=\sum_{i=0}^{n-1}a_ix^i\)的多项式,于是我们可以搞它(

------------------------来了-----------------------

Cooley-Tukey算法:

DFT(系数转点值):

\[A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\\ =\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_n^{2ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_n^{2ki}\\ =A'(\omega_n^{2ki})+\omega_n^kA''(\omega_n^{2ki})\\ 因为\omega_n^{\frac{n}{2}}=-1\\ 所以A(\omega_n^{k+\frac{n}{2}})=A'(\omega_n^{2ki})-\omega_n^kA''(\omega_n^{2ki})\\ 于是可以二分\\ 观察到二分时是按奇偶分类\\ 因此按照蝴蝶定理\\ 奇偶分类一直到最下层就是按二进制反序排列\\ 因此可以预处理\\ \]

IDFT(点值转系数):

\[把\omega_n^k改成\omega_n^{-k},然后把dft得到的点值A(\omega_n^i)作为a_i,再做一遍dft得到n次函数B(x)\\ 使C(x)=\frac{1}{n}B(x)\\ 则C(i)为A(x)的i次项系数 \]

即得易见平凡,仿照上例显然.留作习题答案略,读者自证不难.
反之亦然同理,推论自然成立,略去过程QED,由上可知证毕.

洛谷P1919 【模板】A*B Problem升级版（FFT快速傅里叶）

#include<bits/stdc++.h>
#define N 200005
#define pi 3.1415926535
using namespace std;
struct Complex{double x,y;};
Complex operator+(Complex a,Complex b){return(Complex){a.x+b.x,a.y+b.y};}
Complex operator-(Complex a,Complex b){return(Complex){a.x-b.x,a.y-b.y};}
Complex operator*(Complex a,Complex b){return(Complex){a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};}
Complex a[N],b[N],c[N];
int rota[N],cnt;
void pre(int n,int m){
    cnt=1;
    int high=0;
    while(cnt<=n+m)cnt<<=1,high++;
    for(int i=0;i<cnt;i++)
        rota[i]=(rota[i>>1]>>1)|((i&1)<<(high-1));
}
void fft(int lim,Complex *buf,int dft){
    for(int i=0;i<lim;i++)if(i<rota[i])swap(buf[i],buf[rota[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        Complex wn=(Complex){cos(dft*pi/mid),sin(dft*pi/mid)};
        for(int s=0;s<lim;s+=(mid<<1)){
            Complex w=(Complex){1,0};
            for(int k=s;k<s+mid;k++,w=w*wn){
                Complex x,y;
                x=buf[k],y=w*buf[k+mid];
                buf[k]=x+y;
                buf[k+mid]=x-y;
            }
        }
    }
    if(dft==-1)for(int i=0;i<lim;i++)buf[i].x/=lim;
}
int n,ans[N];
char ch[N];
signed main(){
    scanf("%d",&n);
    pre(n,n);
    scanf("%s",ch);
    for(int i=0;i<n;i++)a[i].x=ch[n-1-i]-'0';
    scanf("%s",ch);
    for(int i=0;i<n;i++)b[i].x=ch[n-1-i]-'0';
    fft(cnt,a,1);
    fft(cnt,b,1);
    for(int i=0;i<=cnt;i++)c[i]=a[i]*b[i];
    fft(cnt,c,-1);
    for(int i=0;i<=cnt;i++){
        ans[i]+=(int)(c[i].x+0.1) ;
        if(ans[i]>=10)
            ans[i+1]+=ans[i]/10,ans[i]%=10,cnt+=(i==cnt);
    }
    while(!ans[cnt]&&cnt>=1)cnt--;
    for(int i=cnt;i>=0;i--)printf("%d",ans[i]);
    return 0;
}
/*
2
25
25

*/