LOADING

加载过慢请开启缓存 浏览器默认开启

FFT 学习笔记

定理:

把多项式 \(A(x)\) 的离散傅里叶变换结果作为另一个多项式 \(B(x)\) 的系数
取单位根的倒数即 \(\omega^{0}_{n},\omega^{-1}_{n},\omega^{-2}_{n},\dots ,\omega^{-(n-1)}_{n}\) 作为x代入 \(B(x)\)
得到的每个数再除以 \(n\) ,得到的是 \(A(x)\) 的各项系数。
(我并不想证明)(挺显然的其实(挨打))
> 单位根:将一个单位圆等分为几分后所得的几个点的值(是一个数,高端表达: \(\in\mathbb{C}\))
>> 单位圆:半径为1的圆
>> 复数:由实部和虚部组成
>>> 实部:实数部分
>>> 虚部:虚数部分
>>>> 虚数:理论中的数,由\(i\)表示
>>>> \(i=\sqrt{-1}\)

\(\omega _{n}^{k}(\cos \cfrac{k}{n} 2\pi ,\sin \cfrac{k}{n} 2\pi )\) (简单的三角函数())
复数为 $ 2 +i 2 $

typedef complex<double> fs;
const double pi=acos(-1);
omega[k]=fs(cos(k/n*2*pi),sin(k/n*2*pi));
inv[k]=conj(omega[k]);//omega[k]的共轭复数(即复数的虚部取反),当复数模长为1时,结果等于原复数的倒数

有:

\[ \large \begin{aligned} \omega ^{k+\frac{n}{2} }_{n}&=- \omega ^{k}_{n}\\ \omega^{k}_{n}&=\omega^{k}_{\frac{n}{2}}(结合图形思考)\\ f(x)&=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+\dots +a_{n-1}x^{n-1}\\ &=(a_{0}+a_{2}x^{2}+a_{4}x^{4}+\dots +a_{n-2}x^{n-2})+(a_{1}x+a_{3}x^{3}+a_{5}x^{5}+\dots +a_{n-1}x^{n-1})\\ \omega^{n}_{n}&=1\\ \end{aligned} \]

令:

\[ \large \begin{aligned} f1(x)&=(a_{0}+a_{2}x+a_{4}x^{2}+\dots +a_{n-2}x^{\frac{n}{2}})\\ f2(x)&=(a_{1}+a_{3}x+a_{5}x^{2}+\dots +a_{n-1}x^{\frac{n}{2}})\\ \end{aligned} \]

则有:

\[ \large \begin{aligned} f(x)&=f1(x^{2})+x\cdot f2(x^{2})\\ f(\omega ^{k}_{n})&=f1(\omega ^{2k}_{n})+\omega _{n}^{k+\frac{n}{2}}\cdot f2(\omega ^{2k}_{n})\\ &=f1(\omega ^{k}_{\frac{n}{2} })+\omega _{n}^{k}\cdot f2(\omega ^{k}_{\frac{n}{2} })\\ f(\omega _{n}^{k+\frac{n}{2} })&=f1(\omega _{n}^{2k+n})+\omega _{n}^{k+\frac{n}{2}}\cdot f2(\omega _{n}^{2k+n})\\ &=f1(\omega^{2k}_{n}\cdot \omega^{n}_{n})+\omega _{n}^{k+\frac{n}{2}}\cdot f2(\omega^{2k}_{n}\cdot \omega^{n}_{n})\\ &=f1(\omega _{\frac{n}{2} }^{k}\cdot \omega _{n}^{n})+\omega _{n}^{k+\frac{n}{2}}\cdot f2(\omega _{\frac{n}{2} }^{k}\cdot \omega _{n}^{n})\\ &=f1(\omega _{\frac{n}{2} }^{k}\cdot \omega _{n}^{n})-\omega _{n}^{k}\cdot f2(\omega _{\frac{n}{2} }^{k}\cdot \omega _{n}^{n})\\ &=f1(\omega _{\frac{n}{2} }^{k}\cdot 1)-\omega _{n}^{k}\cdot f2(\omega _{\frac{n}{2} }^{k}\cdot 1)\\ &=f1(\omega _{\frac{n}{2} }^{k})-\omega _{n}^{k}\cdot f2(\omega _{\frac{n}{2} }^{k})\\ \end{aligned} \]

综上:

\[ \Large f(\omega ^{k}_{n})= \left\{\begin{matrix} f1(\omega ^{k}_{\frac{n}{2} })+\omega _{n}^{k}\cdot f2(\omega ^{k}_{\frac{n}{2} })(0\le k< \frac{n}{2} )\\ \\ f1(\omega _{\frac{n}{2} }^{k-\frac{n}{2}})-\omega _{n}^{k-\frac{n}{2}}\cdot f2(\omega _{\frac{n}{2} }^{k-\frac{n}{2}})(\frac{n}{2}\le k<n)\\ \end{matrix}\right. \]

可是若直接这样进行递归,会较慢
可以先预处理出所有 \(\omega_{n}^{k}\)\(\omega_{n}^{-k}\)
可是这样依然较慢(会快一点,但不多())
所以我就没有预处理(

我们发现递归时转移是这样的(以8个为例)

\[ \omega^{0}_{8}\omega^{1}_{8}\omega^{2}_{8}\omega^{3}_{8}\omega^{4}_{8}\omega^{5}_{8}\omega^{6}_{8}\omega^{7}_{8}\\ ↓\\ \omega^{0}_{8}\omega^{2}_{8}\omega^{4}_{8}\omega^{6}_{8}│\omega^{1}_{8}\omega^{3}_{8}\omega^{5}_{8}\omega^{7}_{8}\\ ↓\\ \omega^{0}_{8}\omega^{4}_{8}┇\omega^{2}_{8}\omega^{6}_{8}│\omega^{1}_{8}\omega^{5}_{8}┇\omega^{3}_{8}\omega^{7}_{8}\\ ↓\\ \omega^{0}_{8}┇\omega^{4}_{8}│\omega^{2}_{8}┇\omega^{6}_{8}║\omega^{1}_{8}┇\omega^{5}_{8}│\omega^{3}_{8}┇\omega^{7}_{8}\\ \]

可能看不出什么
换成二进制看看

\[ 000┆001┆010┆011┆100┆101┆110┆111\\ ↓\\ 000┆010┆100┆110│001┆011┆101┆111\\ ↓\\ 000┆100│010┆110║001┆101│011┆111\\ ↓\\ 000│100║010│110║║001│101║011│111 \]

这下可以看出
对于任意一个数 \(a\) ,它在转移时,最终会被转移到它的二进制数倒转后对应的十进制数的位置
这样说可能有点绕

eg:

\[ 6_{(10)}=110_{(2)}\\ 110_{(2)}\to 011_{(2)}\\ 011_{(2)}=3_{(10)} \]

故原来的 6 经过转移后到了原来 3 的位置

那么怎么进行二进制翻转呢?

一、朴素法

对于每一个数对应的二进制数
我们可以一位一位地进行翻转

优点:方便理解,也方便打
缺点:较慢 \(O(n\log{n})\)

二、递推法

我觉得很人类智慧

\(R(x)\)\(x\) 所对应的二进制数翻转后的结果
如:\(R(6)=3\) (不知道为啥的话,见上文eg处)
因为我们是从小到大推的
所以我们在算 \(R(x)\) 时,已经算过了 \(R(\left \lfloor \frac{x}{2} \right \rfloor )\)
我们可以由 \(R{(\left \lfloor \frac{x}{2} \right \rfloor )}\) 转移到 \(R(x)\)
我们在得出 \(R{(\left \lfloor \frac{x}{2} \right \rfloor )}\) 后,会 惊奇的 发现 (这怎么发现的啊qwq)
$ $的值居然就是原来 \(R(x)\) 只有最低位没有翻转后的结果,所以我们只需要再处理一下即可
故:$R(x)= +2^{k-1}*x $
其中 \(k\) 为二进制数的长度,即 \(\log{n}\)
eg:

\[ \begin{aligned} 已知R(3)&=6\\ R(6)&=\left \lfloor \cfrac{ R(\left \lfloor \frac{6}{2} \right \rfloor )}{2} \right \rfloor +2^{3-1}*6\bmod 2\\ &=\left \lfloor \frac{ R(3)}{2} \right \rfloor +2^{2}*0\\ &=\left \lfloor \frac{ R(011_{(2)})}{2} \right \rfloor +0\\ &=\left \lfloor \frac{110_{(2)}}{2} \right \rfloor \\ &=\left \lfloor \frac{6}{2} \right \rfloor \\ &=3\\ \end{aligned} \]

代码如下:

for(int i=0;i<n;i++){
    r[i]=r[i/2]/2+pow(2,log2(n)-1)*(i&1);
}

ps:这个序列(注意,只是这个序列!!!)可以作为预处理来处理,因为进行 FFT 时, \(n\) 是不发生改变的

优点:较快 \(O(n)\)
缺点:比较难理解

这样的话,我们就可以写出递推,而不是递归
效率就会高很多
好,那我们怎么进行递推来得出答案呢?
参见oi-wiki蝶形运算优化

根据这个方法的性质
我们知道它是从小区间转移到大的区间 \((2\to 4\to 8\to \dots \to n)\)

所以我们 考虑区间DP 利用区间DP的思想(可以这么说吧)
设当前小区间的长度为 \(k\)
先枚举每个小区间的起点
每个小区间内:
\(\cfrac{k}{2}\) 个数对应的 \(f\) 就是原来所说的 \(f1\)
\(\cfrac{k}{2}\) 个数对应的 \(f\) 就是原来所说的 \(f2\)
再根据先前的结论

\[ \Large f(\omega ^{k}_{n})= \left\{\begin{matrix} f1(\omega ^{k}_{\frac{n}{2} })+\omega _{n}^{k}\cdot f2(\omega ^{k}_{\frac{n}{2} })(0\le k< \cfrac{n}{2} )\\ \\ f1(\omega _{\frac{n}{2} }^{k-\frac{n}{2}})-\omega _{n}^{k-\frac{n}{2}}\cdot f2(\omega _{\frac{n}{2} }^{k-\frac{n}{2}})(\cfrac{n}{2}\le k<n)\\ \end{matrix}\right. \]

可以将这 \(k\) 个数转移到新的位置上
代码如下

for(int k=2;k<=n;k*=2){
    fs omega=fs(cos(2*pi/k),sin(z*2*pi/k));
    for(int i=0;i<n;i+=k){
        fs w=fs(1,0);
        for(int j=i;j<i+k/2;j++){
            fs zs=w*c[j+k/2];
            c[j+k/2]=c[j]-zs;
            c[j]+=zs;
            w*=omega;
        }
    }
}

ps:在 IDFT(快速傅里叶逆变换)的时候,要乘上 \(-1\) 即:

if(z==-1){//z是判断是否为IDFT的变量,当z==-1是则当前为IDFT
    for(int i=0;i<n;i++){
        c[i].real(c[i].real()/n);
    }
}

最后只需要

fft(a,1);//将a数进行FFT
fft(b,1);//将b数进行FFT
for(int i=0;i<n;i++) a[i]*=b[i];//相乘
fft(a,-1);//将所得到的a进行逆FFT(IDFT)

最后的最后,只需要把 \(a\) 数组的答案输出就行了

for(int i=0;i<n;i++){
    ans[i]+=a[i].real()+0.5;
    ans[i+1]+=ans[i]/10;
    ans[i]%=10;
}

注意 \(ans\) 数组是倒着来的,需要倒着输出
而且因为 \(n\) 的值必须为 2 的整数次幂
所以会有前导零,需要去掉(我的方法较笨,不太美观,这里是可以优化的)
代码如下:

int q;
for(int i=n-1;i>=0;i--){
    if(ans[i]!=0){
        q=i;
        break;
    }
}
for(int i=q;i>=0;i--) cout<<ans[i];

好,到这里就结束了
AC 代码如下:

#include<bits/stdc++.h>
using namespace std;

const int MAXN=4e6+10;
const double pi=acos(-1);
typedef complex<double> fs;

char sa[MAXN],sb[MAXN];
double n=1;
int lena,lenb,r[MAXN<<1],ans[MAXN<<1];
fs a[MAXN],b[MAXN];

void fft(fs c[],int z){
    for(int i=0;i<n;i++) if(i<r[i]) swap(c[i],c[r[i]]);
    for(int k=2;k<=n;k*=2){
        fs omega=fs(cos(2*pi/k),sin(z*2*pi/k));
        for(int i=0;i<n;i+=k){
            fs w=fs(1,0);
            for(int j=i;j<i+k/2;j++){
                fs zs=w*c[j+k/2];
                c[j+k/2]=c[j]-zs;
                c[j]+=zs;
                w*=omega;
            }
        }
    }
    if(z==-1){
        for(int i=0;i<n;i++){
            c[i].real(c[i].real()/n);
        }
    }
}

int main(){
    cin>>sa>>sb;
    lena=strlen(sa);
    lenb=strlen(sb);
    while(n<lena+lenb) n*=2;
    for(int i=0;i<n;i++) r[i]=r[i/2]/2+pow(2,log2(n)-1)*(i&1);
    for(int i=0;i<lena;i++) a[i].real(sa[lena-i-1]-'0');
    for(int i=0;i<lenb;i++) b[i].real(sb[lenb-i-1]-'0');
    fft(a,1);
    fft(b,1);
    for(int i=0;i<n;i++) a[i]*=b[i];
    fft(a,-1);
    for(int i=0;i<n;i++){
        ans[i]+=a[i].real()+0.5;
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    int q;
    for(int i=n-1;i>=0;i--){
        if(ans[i]!=0){
            q=i;
            break;
        }
    }
    for(int i=q;i>=0;i--) cout<<ans[i];
    return 0;
}