多项式
- Q: 什么是 \(FFT\)?
A: \(Fast Fourier Transform\) 快速傅里叶变换。 - Q: 什么是 \(DFT\)?
A: \(Discrete Fourier Transform\) 离散傅里叶变换。 - Q: 什么是 \(IDFT\)?
A: \(Inverse Discrete Fourier Transform\) 逆·离散傅里叶变换。 - Q: 什么是 \(NTT\)?
A: \(number theory Transform\) 数论变换。 - Q: 什么是 \(MTT\)?
A: 任意模数数论变换。
\(FFT\)
或许有很多人会在这里强调好些引理,但我觉得在推导到那步时能够继续推下去就可以了。
没有必要过分强调。。。
这里将多项式和函数直接对应
其中 \(f_{x_i}=f(x_i)\)
该方法用于在\(O(n\log n)\)求离散卷积。
什么是离散卷积?
给出定义,令\(y=f\times g\)
或者这么表示。
此时 所构成的多项式\({y}\) \(=\) \(f\times g\)
所以譬如我们都知道的十进制竖式乘法就是一种卷积。
引入概念:点值表达式
一个多项式可以表达成 \(f=\sum\limits_{i=0}^{n}a_i\times x^i\),我们称为其为系数表达式
而此时的多项式我们可以理解成一个 \(n\) 次的函数
而在初中二年级我们就知道,任意 \(n+1\) 个点都可以确定一个 \(n\) 次的函数。
这就是为什么小学一些让你找规律的题其实并没有准确答案。
记点集为\(\begin{Bmatrix}(x_0,y_0),(x_1,y_1)\cdots(x_n,y_n)\end{Bmatrix}\)
则\(\forall i\in[0,n] f(x_i)=y_i\)
由此我们将系数表达 \(\begin{Bmatrix} a_0,a_1\cdots a_n\end{Bmatrix}\) 转换到点值表达 \(\begin{Bmatrix}(x_0,y_0),(x_1,y_1)\cdots(x_n,y_n)\end{Bmatrix}\) 的操作称为一次傅里叶变换
而从点值重新转回系数表达的插值操作称为一次逆傅里叶变换
点值表达式的好处就在于,当两个点值表达式相乘时:
即\(\ \ \ \ \ h(x_i)=y_i \times z_i\),\(h\) 的点值表达则为 \(\begin{Bmatrix}(x_0,y_0\times z_0),(x_1,y_1\times z_1)\cdots(x_n,y_n\times z_n)\end{Bmatrix}\)
我们实现了 \(O(n)\) 的离散卷积乘法。
但是我们还不知道如何将系数表达式转换为点值表达式
一般的我们可以,随便找 \(n+1\) 个数然后代入计算。
但 \(O(2n)\) 的复杂度我们并不允许。
我们突然意识到,随便取 \(n+1\) 个数是不可能的,这辈子都不可能。
这里我们引入概念,\(n\) 次单位根的概念:
满足 \(x^n=1\) 的所有的复数 \(x\)。
引入复数概念
没啥特别的,高中基础吧。
其中 记 \(Re(z) = a\) 表示复数 \(z\)的实部。
记 \(Im(z)=b\)表示复数 \(z\) 的虚部。
记 \(\mid z \mid=\sqrt{a^2+b^2}\) 表示复数 \(z\) 的模长。
复数加法法则:实部相加,虚部相加。
而加法的实质意义是平行四边形法则。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1000.jpg %}
复数乘法法则:相当于拆括号。
而乘法的实质意义是旋转相似。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1001.jpg %}
每一个复数 \(z\) 将 \((Re(z),Im(z))\) 在平面上表示出。
发现对于任意的一个复数都可以在平面上表示。
于是乎,我们将这个平面称为复平面。
复数相乘几何意义:模长相乘,副角相加。
引入单位根概念
2次单位根:\(1\),\(-1\)。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0800.jpg %}
3次单位根:\(1\),\(\omega\),\(\omega ^{2}\)。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0801.jpg %}
4次单位根:\(1\),\(i\),\(-1\),\(,-i\)。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0802.jpg %}
我们发现 \(2^k\)这种单位根有很好的轴对称以及中心对称性,为了简化问题,下文所说的 \(n\) 均默认 \(n=2^k,k\in Z\)。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0803.jpg %}
我们沿逆时针一圈把单位根标序号,对于第 \(i\) 个我们记录为 \(\omega_{n}^{i-1}\).
于是我们生成了
\(\omega_{n}^{0}\),\(\omega_{n}^{1}\cdots \omega_{n}^{n-1}\)
\(\omega_{n}^{k}\) 性质何在?
首先我们这么理解为 \(\omega_{n}\)的 \(k\) 次幂
由欧拉定理:
其中 \(\varTheta\) 为复数的辐角。
根据其定义,有以下性质。
-
\[\omega_{2n}^{2k}=\omega_{n}^{k} \]
-
\[\omega_{n}^{k}=-\omega_{n}^{k+\frac{n}{2}} \]
我们已经确定了选那些数了,剩下的只剩带入求出该点的函数值!
快速傅里叶变换
我们如何快速求出他们 \(\omega_{n}^{0}\) , \(\omega_{n}^{1}\cdots \omega_{n}^{n-1}\) 的函数值?
按次数奇偶划分。
划分成两个多项式。
发现:
发现什么了吗? 在 \(\omega_{n}^{k}\) 和 \(\omega_{n}^{k+\frac{n}{2}}\) 的函数值,只差个负号。
而我们每次这么操作时所需求解的区间就会折半,分治的思想。
所以快速傅里叶复杂度保证在 \(O(n\log n)\) 了。
逆·快速傅里叶变换
我们引进范德蒙德矩阵
这里的 \((x_i,y_i)\) 即为点值表达。
\(a_i\) 为系数表达。
我们记录
在得知 \((x_i,y_i)\) 后我们只要算出 \(V\) 的逆矩阵 \(V^{-1}\)
好了,现在大力求逆矩阵就完了,所以我们开始高斯消元。
ちょっとまって
这尼玛不还是 \(O(n^3)\)?
哦,我们突然发现(真的是发现)好像\(V^{-1}_{i\ ,\ j}=\dfrac{\omega _{n}^{-ij}}{n}\)
你丫不是扯淡?
证明
发现\(\omega_{n}^{k(i-j)}\)为等比数列
当 \(i\neq j\) 时
当 \(i=j\) 时,原式 \(=\dfrac{n\times \omega_{n}^{0}}{n}=1\)
综上:得出 \(P\) 为 \(n\) 阶单位阵。
根据逆矩阵定义
这时我们只要在 \(FFT\) 的过程中记录一个 \(flag\) 在系数上\(×-1\) 即可
void fft(Complex* a,int len,int f)
{
if(len==1) return;
Complex a0[len/2],a1[len/2];
for(int i=0;i<len;i++,i++)
{
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}
fft(a0,len/2,f); fft(a1,len/2,f);
Complex wn(cos(f*2.0*pie/len),sin(f*2.0*pie/len));
Complex w(1,0);
for(int i=0;i<(len/2);i++)
{
a[i]=a0[i]+w*a1[i];
a[i+len/2]=a0[i]-w*a1[i];
w=w*wn;
}
}
还没完!
\(FFT\)还没完,我们还要减小常数
1.小trick
w*a1[i]
被算了两遍,没意义,所以
for(int i=0;i<(len/2);i++)
{
Complex t=w*a1[i];
a[i]=a0[i]+t;
a[i+len/2]=a0[i]-t;
w=w*wn;
}
说实话,仅仅这么改屁用没有。
2.我不用递归啦
观察最后得到的数列:
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1004.jpg %}
观察化成二进制后正好反转了。
最后的数列我们得到了。
得到递推式:
for(int i=0;i<=lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
然后从底层向上迭代。
具体如下:
const double p=acos(-1.0);
inline void fft(Complex *a,int f)
{
for(int i=0;i<lim;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid=mid<<1)
{
Complex wn(cos(1.0*p/mid),f*sin(1.0*p/mid));
for(int r=mid<<1,j=0;j<lim;j+=r)
{
Complex w(1.0,0.0);
for(int k=0;k<mid;k++,w=w*wn)
{
Complex x=a[k+j],y=w*a[k+j+mid];
a[k+j]=x+y;
a[k+j+mid]=x-y;
}
}
}
}
\(FFT\) 到此结束了。
\(NTT\)
由于\(FFT\) 炸精还慢,在模数特别时,\(NTT\)诞生了。
阶
在\(\bmod\ m\) 意义下,当 \((m,a)=1\) 时,使得 \(a^r\equiv 1 \pmod m\) 的最小的 \(r\) ,叫做 \(a\) 关于 \(m\) 的阶。记为 \(\delta_{n}(a)=r\) 。
性质:
- 若 \((a,m)=1\),且 \(a^n=1\pmod m\),则 \(\delta_{m}(a)\mid n\)。
- \(\delta_{m}(a)\mid \varphi(m)\)。
原根
若 \(\delta_{m}(a)=\varphi(m)\) 则称 \(a\) 是 \(\bmod\ m\)的一个原根。
性质:
- 原根存在 \(\Leftrightarrow\) \(m=2,4,p^e,2p^e\)。
- \(m\) 有 \(\varphi(\varphi(m))\) 个原根。
- 设 \(g\) 为 \(m\) 的一个原根,那么所有的原根为\(g,g^2,g^3,\ \cdots \ ,g^{\varphi(i)}\)。
发现有些性质和单位根相同,所以我们就使用原根代替单位根,实现快速数论卷积。
实现与 \(FFT\) 类似。
\(MTT\)
我们的 \(NTT\) 十分依赖模数,这使 \(NTT\) 很鸡肋。
理论基础(雾(
一般认为,两个FFT跑得比三个NTT稍微快一点。
- 有一种方法是将非 \(NTT\) 模数拆成 \(3\) 个原根模数,最后 \(Crt\) 合并。但是这是九次 \(NTT\) 慢死你。。。
- 根据command_block大佬的博客,这里提供一种五次 \(FFT\) 的思路。
首先是拆系数我们将 \(f\) 拆成 \(f=a_1\times d + b_1,d=2^{15}\);\(g=a_2\times d + b_2,d=2^{15}\)
\(f\times g=a_1\times a_2\times d^2+(a_1\times b_2+a_2\times b_1)\times d+b_1\times b_2\)
emmm 四次 \(DFT\),三次 \(IDFT\),相当于 \(7\)次\(FFT\) 慢死了。。
ちょっとまって
我们看:这种结构让我们想到了虚数乘法时的样子。
我们设 \(A_1\) 代表 \(a_1\),\(B_1\) 代表 \(b_1\),\(A_2\) 代表 \(a_2\),\(B_2\) 代表 \(b_2\)。
设复多项式 \(F=A_1+iA_2\) , \(G=B_1+iB_2\)
\(P_1=F\times G=(A_1\times B_1-A_2\times B_2)+i(A_1\times B_2+A_2\times B_1)\)
再设 \(H=A_1-iA_2\)
\(P_2=H\times G=(A_1\times B_1+A_2\times B_2)+i(A_1\times B_2-A_2\times B_1)\)
那么 \(P_1+P_2=2(A_1\times B_1+iA_1\times B_2)\)
\(A_1\times B_1=Re(P_1+P_2)\) , \(A_1\times B_2=Im(P_1+P_2)\)
所以现在可以将\(A_1\times A_2,A_1\times B_2,A_2\times B_1,B_1\times B_2\)
所以将 \(F,G,H\) 转为点值表达花费 \(3\) 次 \(FFT\)。
\(P_1,P_2\) 转回系数表达花费 \(2\) 次 \(FFT\)。
总记 \(5\) 次 \(FFT\)
\(5FFT<9NTT-FFT\) 优化胜利啦
\(FWT\)
待更。。。