fft

· · 个人记录

  1. Q: 什么是 FFT? A: Fast Fourier Transform 快速傅里叶变换。
  2. Q: 什么是 DFT? A: Discrete Fourier Transform 离散傅里叶变换。
  3. Q: 什么是 IDFT? A: Inverse Discrete Fourier Transform 逆·离散傅里叶变换。
  4. Q: 什么是 NTT? A: number theory Transform 数论变换。
  5. Q: 什么是 MTT? A: 任意模数数论变换。

FFT

或许有很多人会在这里强调好些引理，但我觉得在推导到那步时能够继续推下去就可以了。

没有必要过分强调。。。

这里将多项式和函数直接对应

其中 f_{x_i}=f(x_i)

该方法用于在O(n\log n)求离散卷积。

什么是离散卷积？

给出定义，令y=f\times g

y_n=\sum _{i=0}^{n}f_i\times g_{n-i}

或者这么表示。

y_n=\sum _{i=- \infty}^{\infty}f_i\times g_{n-i}

此时 所构成的多项式{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} 的操作称为一次傅里叶变换。

而从点值重新转回系数表达的插值操作称为一次逆傅里叶变换

点值表达式的好处就在于，当两个点值表达式相乘时：

f=\begin{Bmatrix}(x_0,y_0),(x_1,y_1)\cdots(x_n,y_n)\end{Bmatrix} g=\begin{Bmatrix}(x_0,z_0),(x_1,z_1)\cdots(x_n,z_n)\end{Bmatrix} h=f\times g h=\begin{Bmatrix}(x_0,y_0),(x_1,y_1)\cdots(x_n,y_n)\end{Bmatrix}\times \begin{Bmatrix}(x_0,z_0),(x_1,z_1)\cdots(x_n,z_n)\end{Bmatrix} \forall i\in[0,n],\ \ \ \ \ h(x_i)=f(x_i)\times g(x_i)

即\ \ \ \ \ 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(n^2) 的复杂度我们并不满意。

我们突然意识到，随便取 n+1 个数是不可能的，这辈子都不可能。

这里我们引入概念，n 次单位根的概念： 满足 x^n=1 的所有的复数 x。

引入复数概念

没啥特别的，高中基础吧。

z=a+bi\ ,\ \ (a\in \Re,\ b\in \Re,\ i^2=-1\ )

其中 记 Re(z) = a 表示复数 z的实部。

记 Im(z)=b表示复数 z 的虚部。

记 \mid z \mid=\sqrt{a^2+b^2} 表示复数 z 的模长。

复数加法法则：实部相加，虚部相加。

Eg:(1+3i)+(2-i)=(3+2i)

而加法的实质意义是平行四边形法则。

复数乘法法则：相当于拆括号。

Eg:(1+3i)\times(2+2i)=1\times 2+3i\times 2+1\times 2i+3i\times 2i=2-6+6i+2i=-4+8i

而乘法的实质意义是旋转相似。

每一个复数 z 将 (Re(z),Im(z)) 在平面上表示出。

发现对于任意的一个复数都可以在平面上表示。

于是乎，我们将这个平面称为复平面。

复数相乘几何意义：模长相乘，副角相加。

引入单位根概念

2次单位根：1,-1。

3次单位根：1,\omega,\omega ^{2}。

4次单位根：1,i,-1,,-i。

我们发现 2^k这种单位根有很好的轴对称以及中心对称性，为了简化问题，下文所说的 n 均默认 n=2^k,k\in Z。

我们沿逆时针一圈把单位根标序号，对于第 i 个我们记录为 \omega_{n}^{i-1}.

于是我们生成了

\omega_{n}^{0}$,$\omega_{n}^{1}\cdots \omega_{n}^{n-1} 首先我们这么理解为 $\omega_{n}$的 $k$ 次幂 由欧拉定理： $$e^{i\varTheta}=\cos {\varTheta}\ \ + i\sin{\varTheta}$$ 其中 $\varTheta$ 为复数的辐角。 根据其定义，有以下性质。 $$\omega_{n}^{k}=e^{ i\frac{2\pi k}{n}}=\cos {\dfrac{2 \pi k}{n}}+i\sin{\dfrac{2 \pi k}{n}}$$ 1. $$\omega_{2n}^{2k}=\omega_{n}^{k}$$ 3. $$\omega_{n}^{k}=-\omega_{n}^{k+\frac{n}{2}}$$ 我们已经确定了选那些数了，剩下的只剩带入求出该点的函数值！ --- ### 快速傅里叶变换 我们如何快速求出他们 $\omega_{n}^{0}$ , $\omega_{n}^{1}\cdots \omega_{n}^{n-1}$ 的函数值？ $$f(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1},\ \ (x=\omega_{n}^{k})$$ 按次数奇偶划分。 $$f(x)=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1}),\ \ (x=\omega_{n}^{k})$$ 划分成两个多项式。 $$f(x)=f_1(x^2)+xf_2(x^2),\ \ (x=\omega_{n}^{k})$$ 发现： $$f(\omega_{n}^{k})=f_1(\omega_{n}^{2k})+\omega_{n}^{k}f_2(\omega_{n}^{2k})$$ $$f(\omega_{n}^{k+\frac{n}{2}})=f_1((-\omega_{n}^{k})^2)-\omega_{n}^{k}f_2((-\omega_{n}^{k})^2)$$ $$\Leftrightarrow$$ $$f(\omega_{n}^{k+\frac{n}{2}})=f_1(\omega_{n}^{2k})-\omega_{n}^{k}f_2(\omega_{n}^{2k})$$ 发现什么了吗？ 在 $\omega_{n}^{k}$ 和 $\omega_{n}^{k+\frac{n}{2}}$ 的函数值，只差个负号。 所以我们可以只求出 $[0,\dfrac{n}{2}] $ 的函数值就可以了。 而我们每次这么操作时所需求解的区间就会折半，分治的思想。 所以快速傅里叶复杂度保证在 $O(n\log n)$ 了。 ### 逆·快速傅里叶变换 我们引进范德蒙德矩阵 $$ \left [ \begin{matrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x^{n-1}_{n-1} \end{matrix} \right] \left [ \begin{matrix} a_0 \\ a_1 \\ \vdots\\ a_{n-1} \end{matrix} \right] = \left [ \begin{matrix} y_0 \\ y_1 \\ \vdots\\ y_{n-1} \end{matrix} \right] $$ 这里的 $(x_i,y_i)$ 即为点值表达。 $a_i$ 为系数表达。 我们记录 $$ V=\left [ \begin{matrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x^{n-1}_{n-1} \end{matrix} \right]$$ 在得知 $(x_i,y_i)$ 后我们只要算出 $V$ 的逆矩阵 $V^{-1} \begin{matrix} a_0 \\ a_1 \\ \vdots\\ a_{n-1} \end{matrix} \right] = V^{-1} \left [ \begin{matrix} y_0 \\ y_1 \\ \vdots\\ y_{n-1} \end{matrix} \right]

好了，现在大力求逆矩阵就完了，所以我们开始高斯消元。

ちょっとまって

这尼玛不还是 O(n^3)?

哦，我们突然发现（真的是发现）好像 V_{i\ ,\ j}^{-1}=\dfrac{\omega_{n}^{-ij} }{n}

你丫不是扯淡？

证明

P=V\times V^{-1}\ , P_{\ i,\ j}=\sum _{k=0}^{n-1}\frac{\omega_{n}^{ki}\times\omega _{n}^{-jk}}{n}=\sum _{k=0}^{n-1}\frac{\omega_{n}^{k(i-j)}}{n}

发现\omega_{n}^{k(i-j)}为等比数列

当 i\neq j 时

\dfrac{1-w_{n}^{(i-j)n}}{1-\omega_{n}^{i-j}}

根据定义,

\omega_{n}^{(i-j)n}=(\omega_{n}^{n})^{(i-j)}=1

原式=0

当 i=j 时，原式 =\dfrac{n\times \omega_{n}^{0}}{n}=1 综上：得出 P 为 n 阶单位阵。 根据逆矩阵定义

V^{-1}_{i\ ,\ j}=\dfrac{\omega _{n}^{-ij}}{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;
    }
}

还没完！

#### 1.小trick `w*a1[i]` 被算了两遍，没意义，所以 ``` cpp 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.我不用递归啦 观察最后得到的数列： ![](https://cdn.jsdelivr.net/gh/protons-z/[email protected]/images/Sol/2020-04-1003.jpg) ![](https://cdn.jsdelivr.net/gh/protons-z/[email protected]/images/Sol/2020-04-1004.jpg) 观察化成二进制后正好反转了。 最后的数列我们得到了。 得到递推式： ``` cpp for(int i=0;i<=lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1)); ``` 然后从底层向上迭代。 具体如下： ``` cpp 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 。 性质：

  1. 若 (a,m)=1,且 a^n=1\pmod m,则 \delta_{m}(a)\mid n。

原根

若 \delta_{m}(a)=\varphi(m) 则称 a 是 \bmod\ m的一个原根。 性质：

  1. 原根存在 \Leftrightarrow m=2,4,p^e,2p^e。
  3. 设 g 为 m 的一个原根，那么所有的原根为g,g^2,g^3,\ \cdots \ ,g^{\varphi(i)}。

发现有些性质和单位根相同，所以我们就使用原根代替单位根，实现快速数论卷积。 实现与 FFT 类似。

MTT

我们的 NTT 十分依赖模数，这使 NTT 很鸡肋。

理论基础（雾（

一般认为,两个FFT跑得比三个NTT稍微快一点。

  1. 有一种方法是将非 NTT 模数拆成 3 个原根模数，最后 Crt 合并。但是这是九次 NTT 慢死你。。。
  2. 根据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。

    总记 $5$ 次 $FFT

FWT

待更。。。