定理:
把多项式 \(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;
}