狗都学的会得的FFT
\(FFT\)
一、前置知识
1、多项式的系数表示法
形如 \(f(x)=\sum_{i=0}^{n}a_ix^i\)的式子。
2、多项式的点值表示法
每当\(x\)取一个值\(x_i\)，多项式都会有对应\(f(x_i)\)，把\((x_i,f(x_i))\)这样的看作一个点。一个\(n\)次多项式可以被\(n\)个这样的点表达。
3、卷积
简单来说就是多项式乘法（当然实际上没那么简单）。
如果有多项式\(f(x),g(x)\)，需要求\(F(x) = f(x)*g(x)\)，不难想到\(O(nm)\)的暴力算法。
但事实上我们可以利用\(FFT\)将系数表示法变成点值表示法从而在\(O(nlogn)\)的时间复杂的下得到答案，最后再逆变换成系数表示输出。
4、欧拉公式
欧拉公式\(e^{iθ}=cosθ+isinθ\)，其中 \(e\) 是自然对数的底数。
当且仅当，\(θ\)=\(2k\pi\)时，\(e^{θ}=1\)。
5、单位根
即满足\(w^n=1\)的根，在复数域当中这样的解有很多。
由欧拉公式可知,\(w_i = e^{\frac{2\pi i}{n}}\)，对于这样的根我们记作\(w_n^i = e^{\frac{2\pi i}{n}}\)
6、单位根的性质
（1）折半定理:\(w_{2n}^{2k}=w_{k}^{k}\)
（2）消去引理:\(w_{n}^{k+\frac{n}{2}}=-w_{n}^{k}\)
二、\(DFT\)和\(IDFT\)
出于方便，我们将两个多项式都设置为\(n\)
对于长度为\(n\)的两个多项式，如果用系数表示法求其卷积的时间复杂度为\(O(n^2)\)，而如果使用点值表示法，只需要将点值相乘就可以了这样的时间复杂度是\(O(n)\)的。
现在我们只需要将系数表式法变成点值表示法就好了，这一过程就是\(DFT\),即离散傅里叶变化，而将点值表示法变成系数表式法的过程就是\(IDFT\),即离散傅里叶逆变化。
而\(FFT\)就是优化这两个部分
现在我们得到了一个大概算法流程就是：
- 构造
两个多项式补0，得到两个长度为\(2n\)的多项式，分别叫做\(a,b\)。
（为什么是\(2n\)?两个多项式乘积最高项为\(2n\)） - 求值
计算\(f_a = DFT(a),f_b = DFT(b)\) - 相乘
将\(f_a.f_b\)相乘，得到新的点值表示\(f\) - 逆变换（插值）
用\(FFT\)计算\(IDFT(f)\)。
三、快速傅里叶变化
算法的基本思想是分治，他的分治想法就是将奇偶分离。但这时候出现了一个问题，当多项式项数不是偶数怎么办？
所以在第一步构造的时候，我们需要补充到的不仅仅是2\(n\)，而是大于2\(n\)的一个2的幂次。
对于 \(f(x) = a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}+a_nx^n\),
分离为
\(f(x) = (a_0x^0+a_2x^2+...+a_{n-1}x^{n-1})+(a_1x^1+a_3x^3+...+a_{n}x^{n}\))
奇数部分提出一个\(x\),
\(f(x) = (a_0x^0+a_2x^2+...+a_{n-1}x^{n-1})+x(a_1x^0+a_3x^2+...+a_{n}x^{n-1})\)
令\(f_0 (x)= a_0x^0+a_2x^2+...+a_{n-1}x^{n-1}\)，\(f_1(x) = a_1x^0+a_3x^2+...+a_{n}x^{n-1}\)
显然\(f(x) = f_0(x^2)+xf_1(x^2)\)
带入\(x = w^k_n\)并运用折半定理得，
\(f(w^k_n)= f_1((w^k_n)^2)+w^k_nf_2((w^k_n)^2)\)
\(~~~~~~~~~~~~= f_1(w^{2k}_n)+w^k_nf_2(w^{2k}_n)\)
\(~~~~~~~~~~~~= f_1(w^{k}_{n/2})+w^k_nf_2(w^{k}_{n/2})\)
同样运用消去引理得
\(f(w^{k+n/2}_n)\)\(=f_1(w^{k}_{n/2})-w^k_nf_2(w^{k}_{n/2})\)
不难发现，\(k\)是需要小于\(n/2\)的。所以我们使用前一个公式计算左半部分，后一个公式计算右半部分。
时间复杂度为\(O(nlogn)\)。
四、快速傅里叶逆变换
现在我们已经有一个\(n\)维的点值向量\((A(x_0),A(x_1),...,A(x_{n-1}))\)推出\(n\)维系数向量\((a_0,a_1,...,a_{n-1})\)。
假设\((a_0,a_1,...,a_{n-1})\)的\(DFT\)为\((d_0,d_1,...,d_{n-1})\)。
构造一个这样的的多项式\(F(x) = d_0+d_1x+d_2x^2...,d_{n-1}x^{n-1}\)
设\(c_k = F(w_n^{-k})\)
则\(c_k = \sum_{i=0}^{n-1}d_iw_n^{-ki}\)
带入\(d_i\)，
于是
\(c_k ~= \sum_{i=0}^{n-1}[\sum _{j=0}^{n-1}a_j(w_n^i)^j]w_n^{-ki}\)
\(~~~~~~=\sum _{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(w_n^i)^{j-k}\)
令\(S(j,k)=\sum_{i=0}^{n-1}(w_n^i)^{j-k}\),
设\(j-k = \delta\)
则\(S(j,k)=w_n^0+w^\delta_n+...+w^{(n-1)\delta}_n\)。
公比为\(w^\delta_n\)。
当\(w^\delta_n=1\)，即\(\delta=0\)时
\(S(j,k)=n\)，此时\(\delta=0\)，可得\(j=k\)。
\(w^\delta_n\neq1\)，即\(\delta\neq0\)时,\(S(j,k)=\frac{w_n^0[(w_n^\delta)^n-1]}{w_n^\delta-1} = \frac{1[(w_n^n)^\delta-1]}{w_n^\delta-1}=\frac{1[(1)^\delta-1]}{w_n^\delta-1}=0\),此时\(j\neq k\)
所以\(j=k\)时\(S(j,k)=n\)，\(j\neq k\)时\(S(j,k)=0\)
带入原式：
\(c_k = \sum_{j=0}^{n-1}a_jS(j,k) = a_kn\)
故\(a_k=\frac{c_k}{n}\)。
现在我们对照\(DFT\)就可以得到\(IDFT\)的过程了，在\(DFT\)中我们对一个系数多项式代入\(w_n^k\)的到一个点值向量，而\(IDFT\)就是对一个以点值为系数的系数多项式带入\(w_n^{-k}\)得到的点值向量并将每一项除\(n\)，得\((\frac{c_0}{n},\frac{c_1}{n},...,\frac{c_{n-1}}{n})\)。
注意\(w_n^{-k}\)是\(w_n^k\)的共轭复数。
五 、蝴蝶变换
相对于迭代的算法，递归需要调用内存栈从而浪费多余的时间和空间。
而\(FFT\)分治过程中存在一个规律从而可以让我们把一个递归的算法使用迭代的方式写出来，我们称之为蝴蝶变换。
对于一个八次的多项式\({\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}}\)
在分治过程中是这样(：
原序列：0 1 2 3 4 5 6 7
现序列：0 4 2 6 1 5 3 7
观察二进制表示，发现现序列相当于原序列的二进制翻转。
于是我们可以自底向上的迭代解决这个问题。
六、代码
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <stack>
using namespace std;
#define ll long long
const double pi = acos(-1);
const ll N = 5e6+50;
ll n,m;
ll limit = 1;//二进制位数
ll L,R[N];//二进制位数、二进制翻转数组
struct Complex{
double x,y;
Complex(double x = 0,double y = 0):x(x),y(y){}
}a[N],b[N];
Complex operator*(Complex p,Complex q){
return Complex(p.x*q.x-p.y*q.y,p.x*q.y+p.y*q.x);
}
Complex operator+(Complex p,Complex q){
return Complex(p.x+q.x,p.y+q.y);
}
Complex operator-(Complex p,Complex q){
return Complex(p.x-q.x,p.y-q.y);
}
void fft(Complex *A,ll type){//把一个多项式转换成点值表示法
for(ll i = 0;i < limit;i++)if(i < R[i])swap(A[i],A[R[i]]);
for(ll mid = 1;mid < limit;mid<<=1){
Complex wn(cos(pi/mid),type*sin(pi/mid));//单位根
for(ll len = mid<<1,pos = 0;pos < limit;pos += len){
Complex w(1,0);
for(ll k = 0;k < mid;++k,w = w*wn){
Complex x = A[pos+k];
Complex y = w*A[pos+mid+k];
A[pos+k] = x+y;
A[pos+mid+k] = x-y;
}
}
}
if(type == 1)return ;
for(ll i = 0;i <= limit;i++){
a[i].x /= limit;
}
}
string p,q;
stack<ll> st;
int main() {
cin>>p>>q;n = p.length()-1,m = q.length()-1;
for(ll i = n;i >= 0;i--)a[i].x = p[i]-'0';
for(ll i = m;i >= 0;i--)b[i].x = q[i]-'0';
while(limit <= n+m)limit<<=1,L++;//长度
for(ll i = 0;i < limit;i++){
R[i] = (R[i>>1]>>1) | ((i&1)<<(L-1));
//在原序列中i与i/2的关系是:i是i/2的左移
//那么反转之后就需要右移，同时处理尾数
}
fft(a,1);
fft(b,1);
for(ll i = 0;i <= limit;i++)a[i] = a[i]*b[i];
fft(a,-1);//逆变换
for(ll i = n+m;i > 0;i--){
ll tmp = (ll)(a[i].x+0.5);
st.push(tmp%10);
a[i-1].x+=tmp/10;
}
cout<<(ll)(a[0].x+0.5);
while(!st.empty()){
cout<<st.top();
st.pop();
}
return 0;
}
浙公网安备 33010602011771号