FFT

缘起

高精度乘法，以前给出过O(n^2)的朴素算法和O(n^log3)的karatsuba优化. 但是还有更好更快的优化——FFT（快速傅里叶变化(fast Fourier transform)） 它是基于 DFT（discrete fourier transform 离散傅里叶变换）

fft是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使dft时所需要的乘法次数大幅减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。fft的复杂度是 O(nlogn)的.

fft是dft的快速算法. 它是根据dft的奇、偶、虚、实等特性，对dft算法进行改进获得的. 所以说啊, 算法的根本性的改进总来自于数学理论的突破啊~ 所以学好数学真的重要——不然只能CRUD一辈子~

fft用来加速多项式乘法（高精度乘法只是顺便的结果而已）

分析

本文的主干借鉴了【1】（写的很好），将一个整数的二进制进行翻转主要参考了【2】，而另一个IFFT 变换参考了【3】. 感觉FFT 还是比较好懂和好写的（迫不及待想找道板题来搞了）

首先引入学习FFT之前必须掌握的一些概念

多项式的两种表示方法

n-1次多项式是（n>=1）

1
f(x)=a_{n-1}*x^{n-1}+a_{n-2}*x^{n-2}+....+a1*x+a0

或者使用f的系数唯一刻画f

1
f=(a_{n-1},...,a1,a0)

这就是多项式的系数表示法. 但是有一点线性代数知识的朋友都知道, n-1次实多项式其实还可以使用它（一条曲线）上面的n个点来唯一刻画——因为一旦知道了它穿过的n个点的话, 则就可以通过解方程（n个n元1次方程构成的方程组）来唯一刻画f. 这就是多项式的点值表示法.

高精度乘法和FFT之间的关系

为什么要引入这两种多项式的表示法? 系数表示法是最符合人类理解常识的表示法. 这不消多说什么——因为知道了一个多项式的各个系数的话, 等于就知道了该多项式在所有自变量上取值的表现行为. 这是极好的. 那么点值表示法有什么优势吗? 当然有! 不然我们引进它干什么? 就拿高精度乘法而言. 我们通过前面的学习知道了其实是把一个大整数bigint表示为

1
B = (10^{BASE_LEN*(n-1)})*a_{n-1}+...+(10^{BASE_LEN*0})*a_{0}

BASE_LEN 什么的参见我写的【4】即可. 这不就是f(x)在x=10^{BASE_LEN}处的取值吗? 所以对于两个大整数B1和B2， 我只需要记

1
2
B1 = (10^{BASE_LEN*(n-1)})*a_{n-1}+...+(10^{BASE_LEN*0})*a_{0}
B2 = (10^{BASE_LEN*(m-1)})*b_{m-1}+...+(10^{BASE_LEN*0})*b_{0}

B1是一个n-1次多项式f在x=10^BASE_LEN处的取值，B2是一个m-1次的多项式g在x=10^BASE_LEN处的取值.

所以要计算B1B2这种高精度乘法，只需要求出 f(x)\g(x) 的各个系数即可. 而且只需求出系数即可. 因为我在【4】中的大整数的数据结构就是使用10^BASE_LEN作为一段的.

其中

1
2
f(x) = a_{n-1}*x^{n-1}+...+a1*x+a0
g(x) = b_{m-1}*x^{m-1}+...+b1*x+b0

那么求h(x)= f(x)*g(x)的话最常见的思路就是使用卷积了. 

1
h(x)=f(x)*g(x)=c_{n+m-2}*x^{n+m-2}+....+c1*x+c0

其中

1
ck=a_{k-i}*b_{i}, 0<=i<=k 求和

那么不消说, 我们知道这样求两个多项式乘积的复杂度是O(n^2)的. 也是我们承受不起的. 有什么更好的方法么? 有! 如果使用多项式的点值表示法的话, 一切都很简单了——取n+m-1个点(因为n-1次多项式和m-1次多项式乘积是n+m-2次多项式)x0,…,x_{n+m-2}, 计算出f和g在这n+m-1个点的取值, 然后h在这n+m-1个点的取值就是O（n）完成的——

1
h(xi) = f(xi) * g(xi), 0<=i<n+m-1

完了吗? 完了呀~ 因为前面说了, n个一般意义下的点（所谓一般意义下指的是构成的n元一次方程组是可逆矩阵，有唯一解）就可以完全唯一刻画一个n-1次多项式. 我们现在在O(n)时间内得到了h在n+m-1个点下的值，当然完成了求解f和g乘积的任务.

是不是瞬间有种可以O(n) 暴艹所有高精度题目的冲动??? 其实细想——还差得远.

1
2
3
4
5
6
7
8
首先——计算f（g）在n+m-1个点的取值的复杂度就是O(n^2)的，就是DFT算法，这里假设 n>=m.

其次，得到了h在n+m-1个点的取值之后, 这是不直观的，我们更希望得到h的系数表示法. 所以还要解方程组得到h的
各个系数, 这个过程也是O(n^2)的. 这就是IDFT算法

最后, 得到了h的各个系数之后, 其实是各个段的值，但是显然可能比 10^{BASE_LEN} （10亿）大. 还要考虑进
位问题, 这个倒也简单——从低位开始到高位, 进的位由除以 10^{BASE_LEN}给出, 而余数才是这个段的值. 然后
输出这个大整型的时候，还要从高位开始去掉无意义的前导0.

上述”首先”，可以优化到 O(nlogn), 其中 log 是以2为底的对数. 这就是 FFT.

上述”其次”, 也可以优化到O(nlogn), 其中 log 是以2为底的对数. 这就是 IFFT. 其实后面会看到 FFT和IFFT除了参与的单位根是共轭复数之外没有别的区别. 所以IFFT本质就是FFT.

上述”最后” 到没什么，可以在O(n)内解决.不影响最终复杂度.

好了，讲清楚了 FFT 和高精度乘法之间的关系, 我们就知道了 为什么 高精度乘法可以优化到O(nlogn)了.

下面全力讲清楚 FFT. 以及 IFFT 如何转换为 FFT的.

FFT

再说一遍， FFT 解决的问题是

1
2
对于n-1次多项式 f(x)=a0+a1*x+...+a_{n-1}*x^{n-1},  给定n个点 x0,...,x_{n-1}, 如何在O(nlogn)时间内求出f
在这n个点的取值? 其中 x0,...,x_{n-1} 到底取什么值可以自定义. 其中n 是 2的幂次（对于实际问题, 如果n不是2的幂次, 我们可以预处理使得n是2的幂次——前导补0嘛~）

首先，这n个点 x0,…,x_{n-1}不能随意取.

傅里叶（Fourier）大神告诉我们——取

1
w_{n}^i=e^{i*2*pi/n};  0<=i<n

这n个点代入.为什么? 因为f(wi)(0<=i<n)之间存在关系而不是孤立的. 而且这n个点构成的矩阵是可逆的，保证了解的存在唯一性. 那么它们之间有什么关系呢?

首先知道w_n=e^{2*pi/n}(即上述i取值为1)具备以下两个性质

1
2
w_{n}^k = w_{2n}^{2k}; for 0<=k<n
w_{n}^{k+n/2}=-w_{n}^{k}; for 0<=k<n

这都是用高中数学就能知道的事情（欧拉公式嘛~）

下面进入高（神）能（奇）阶段

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f(x)=a0+...+a_{n-1}*x^{n-1};
=>
f(x)=(a0+a2*x^2+...+a_{n-2}*x^{n-2})+x*(a1+a3*x^2+...+a_{n-1}*x^{n-2})
=>
令
f1(x)=a0+a2*x+...+a_{n-2}*x^{n/2-1}，f2(x)=a1+a3*x+...+a_{n-1}*x^{n/2-1}
则
f(x)=f1(x^2)+x*f2(x^2)
=>
令x=w_{n}^k, 其中 k<n/2 则
f(w_{n}^k)=f1(（w_{n}^k）^2)+w_{n}^k*f2((w_{n}^k)^2)
		 =f1(w_{n}^{2k})+w_{n}^k*f2(w_{n}^{2k})
		 =f1(w_{n/2}^{k})+w_n^{k}*f2(w_{n/2}^{k})
令x=w_{n}^{k+n/2}, 其中 k<n/2, 则
f(w_{n}^{k+n/2})=f1(w_{n}^{2k+n})+w_{n}^{k+n/2}*f2(w_{n}^{2k+n})
		 =f1(w_{n}^{2k}*w_{n}^n)-w_{n}^k*f2(w_{n}^{2k}*w_{n}^n)
		 =f1(w_{n}^{2k})-w_n^{k}*f2(w_{n}^{2k})
		 =f1(w_{n/2}^{k})-w_n^{k}*f2(w_{n/2}^{k})

发现了什么吗?

1
2
3
4
f(w_{n}^k) 和 f(w_{n}^{k+n/2}) ,for 0<=k<n/2, 由两部分构成,
一部分是f1(w_{n/2}^{k})
另一部分是w_n^{k}*f2(w_{n/2}^{k})
相差的只是一个符号而已（根据这两部分求解f(w_{n}^k) 和 f(w_{n}^{k+n/2})的过程叫做蝴蝶变换）.

所以我们知道了

1
2
如果知道了 f1(w_{n/2}^{k}) 和f2(w_{n/2}^{k}) , 0<=k<n/2, 则就能知道f(w_n^k)和f(w_n^{k+n/2}) 的值. 
其中 f=(a0,...,a_{n-1}), f1=(a0,...,a_{n-2}), f2=(a1,.a3,...,a_{n-1})

于是, 我们看到了分治的希望. 而且不难知道FFT的复杂度是 O(nlogn)

IFFT

这里的数学理论只是简单的线性代数. 参考了【3】. 其实从根本上讲，IDFT是解一个n元一次方程组的问题.

写成矩阵形式如下（其实下图左边到右边就是FFT）

有点线性代数知识的都知道只需要两边乘以左边矩阵的逆矩阵就行了. 而左边矩阵的逆矩阵刚好就是

D/n, 其中D为如下矩阵

所以两边乘以逆矩阵得到

这样我们就发现了 IFFT 其实就是把FFT. 只是把FFT用的w_{n}换成 w_{n}_{-1}而已. 但是注意, 最后的结果要除以n（注意要进行丢失精度处理）.

算法实现

既然说到了是分治，那么递归是最爽的、最不负责任的写法了. 不难给出如下递归的写法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include<complex>
#define cp complex<double>
#define fo(i,s,n) for(int (i)=(s); (i)<=(n); ((i)++))

void fft(cp *a,int n,int inv)//inv=1就是fft, inv=-1就是ifft，n是数组a的长度
{
    if (n==1)return;
    int mid=n/2;
    static cp b[MAXN];  // 用static是为了防止递归爆栈
    fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1]; // b是a的偶系数和奇系数分离开来的
    fo(i,0,n-1)a[i]=b[i]; // 借助b作为中转, a将自己奇系数和偶系数分离开来了, 即0,1,2,3变成0,2,1,3
    fft(a,mid,inv),fft(a+mid,mid,inv);//分治,求出f1和f2
    fo(i,0,mid-1) // 根据f1和f2求出f
    {
        cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv取决是否取共轭复数
        b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid]; //蝴蝶变换
    }
    fo(i,0,n-1)a[i]=b[i];// b又做了一回中转站
}

使用上述板子进行多项式乘法的方法只需要像下面这样, 计算 多项式 a*多项式b, 结果存入c中, MAXN大于a的次数+b的次数（n等于这个次数之和）

1
2
3
4
5
6
cp a[MAXN],b[MAXN];
int c[MAXN];
fft(a,n,1),fft(b,n,1);//系数转点值，即fft
fo(i,0,n-1)a[i]*=b[i];
fft(a,n,-1);//点值转系数，即ifft
fo(i,0,n-1)c[i]=(int)(a[i].real()/n+0.5);//注意防止精度丢失

但是这块板子虽然好看，但是不中用. 容易吃T. 因为递归写法的常数太大了——所以必须要使用dp来进行常数优化.

我们注意到如下事实

发现了点儿什么? 递归到最后, 其实就是把按二进制的比特位倒序了. 即 011 变成 110（即a[3]递归最后变动到了a[6]）. 怎么实现呢? 硬写的话很简单的. 总能写出来. 但是网上有菊苣有更好的方法.

1
2
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); // 即rev[i]伊始=i, 经过此式变更, rev[i]就是递归最后最终的索引, 
例如rev[3]=3, 经过此式变更, rev[3]变成了6.

这个式子怎么理解呢? 第一次看到有点丈二和尚摸不着头脑的赶脚. 但是阅读了【2】之后豁然开朗.

rev[i]伊始=i, 其实可以分为 最后一位比特位和其他比特位(记为P).而P就等于i>>1 而采用dp的眼光来看，其实i的比特位逆序不就是P逆序之后再把i的最后一位比特位放到最高比特位么? 而我们假设P已经逆序了（dp假设），所以它逆序后的值不就是 rev[i>>1]么? 但为什么还要 >> 1呢? 因为P的位数比i少一位（所以以0填充最高位）, 所以P逆序之后会多一个后缀0（即P逆序之后的rev[i>>1]的最低位会有一个0），所以要截掉这个0，就需要 rev[i>>1]再次>>1. 而((i&1)<<(bit-1))是什么呢? 当然就是i的最后那个比特位（就是 i&1）移动到最高位来啦（<<(bit-1)的作用）~

所以我们一开始就使用上面的rev将元素归位到应该放置的位置，再一步一步向上合并（这就是dp，自底向上，而前面的递归算法是自顶向下）. 于是我们可以写出如下虽然不像分治递归那样优美易懂但是效率喜人的板子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void fft(cp *a,int n,int inv)
{
    int bit=0;
    while ((1<<bit)<n)bit++;
    fo(i,0,n-1) // 因为rev是基于dp的思想求出来的, 所以要从0到n-1
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就等于是没交换）
    }
    for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一，即如上图中合并a0,a4, 则mid=1,要合并的长度是2
    {
    	cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根，pi的系数2已经约掉了(2*pi/(2*mid)=pi/mid)
        for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度，i是长度为mid的合并序列的起点
		{
            cp omega(1,0);
            for (int j=0;j<mid;j++,omega*=temp)//只扫2*mid长度的区间的左半部分，得到右半部分的答案
            {
                cp x=a[i+j],y=omega*a[i+j+mid]; // 因为i是起点,所以i+j(0<=j<mid)就是区间的左半部分,i+j+mid就是区间的右半部分
                a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换
            }
        }
    }
}

下一篇文章我会用高精度乘法来测一下这个板子.

参考

【1】https://blog.csdn.net/enjoy_pascal/article/details/81478582

【2】https://blog.csdn.net/GGN_2015/article/details/69518685

【3】http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

【4】https://yfsyfs.github.io/2019/08/08/poj-1001-Exponentiation-%E9%AB%98%E7%B2%BE%E5%BA%A6%E4%B9%98%E6%B3%95%E4%B9%8Bkaratsuba%E7%AE%97%E6%B3%95/