笔记：快速傅里叶/数论变换

本文主要讲解 \(O(nlogn)\) 的多项式乘法。

参考了 OI Wiki。

系数表达和点值表达

系数表示法，即由 \(F(x) = \sum_{i=0} a_i x\) 确定的多项式。

由 \(n+1\) 个点可以唯一确定一个 \(n\) 次多项式，即由 \((x_0,F(x_0)),(x_1,F(x_1)),\cdots(x_n,F(x_n))\) 也能确定一个多项式。

两个多项式相乘，若使用系数表示，复杂度通常为 \(O(n^2)\) ，若使用点值表示，只需要把值乘起来，复杂度为 \(O(n)\) 。

快速傅里叶变换

思想

点值表示可以很快的进行多项式乘法，通过选取合适的自变量，可以快速将系数表示转化成点值表示，乘法后再快速转回系数表示。

FFT 基于分治思想，通过选取有对称关系的点，可以在 \(O(nlogn)\) 的复杂度进行转换。

单位复根可以完成这个过程。

过程

对于多项式 \(F(x)=a_0+a_1x+\cdots a_n x^n\)，做如下转换:

\[F(x)=a_0+a_1x+\cdots a_n x^n\\= (a_0+a_2x^2+\cdots a_{n-1}x^{n-1})+x(a_1+a_3x^2+\cdots a_n x^{n-1})\\= G(x^2)+xH(x^2)\]

单位复根有许多优美的性质，设 \(\omega_n^i\) 为 \(n\) 阶第 \(i\) 个单位根（\(n\) 为偶数）：

\[F(\omega_n^i)=G(\omega_n^{2i})+\omega_n^iH(\omega_n^{2i})=\\ G(\omega_{n/2}^{i})+\omega_n^iH(\omega_{n/2}^i)\]

由 \(\omega_n^{i+n/2}=-\omega_n^i\) 得到：

\[F(\omega_n^{i+n/2})=G(\omega_{n/2}^i)-\omega_n^iH(\omega_{n/2}^i) \]

我们可以用 \(G(\omega_{n/2}^{i})\) 和 \(H(\omega_{n/2}^i)\) 表达 \(F(\omega_n^i)\) 和 \(F(\omega_n^{i+n/2})\)。

逆变换

从系数表示到点值表示是一个线性变换的过程:

\[\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{bmatrix}=A \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix}\]

\(A\) 是由单位复根和 \(1\) 组成的矩阵，求解 \(A\) 的逆，为每一项取倒数并处以 \(n\)，而 \(\frac{1}{\omega_k}=e^{-\frac{2\pi i}{k}}=\cos(\frac{2 \pi}{k})+i\sin(-\frac{2 \pi}{k})\)，仅相差一个负号，将 \(\frac{1}{\omega_k}\) 视为单位复根，执行相同的步骤，并在最后除以长度 \(n\) 即可。

优化

上述递归过程可以使用非递归实现。

位逆序变换

对 \(G,H\) 进行分割时，我们会把偶数项移到前面，技术项移到后面，可以证明，分割到最后，\(x_i\) 的位置会从 \(i\) 变成 \(i\) 的二进制翻转的位置。

例如 \(x_{0\cdots 7}\)，转换后变成 \(\{x_0,x_4,x_2,x_6,x_1,x_5,x_3,x_7\}\)。

设 \(f_s\) 表示 \(s\) 的二进制翻转后的值，转移是方便的。

蝴蝶变换

在从下往上递推的过程中，由这两个式子：

\[F(\omega_n^{i})=G(\omega_{n/2}^i)+\omega_n^iH(\omega_{n/2}^i)\\F(\omega_n^{i+n/2})=G(\omega_{n/2}^i)-\omega_n^iH(\omega_{n/2}^i) \]

变换数组中，\(F_i\) 存储了 \(G(\omega_{n/2}^i)\) ，\(F_{i+n/2}\) 存储了 \(H(\omega_{n/2}^i)\)，我们直接计算 \(F(\omega_n^{i})\) 和 \(F(\omega_n^{i+n/2})\) 并覆盖即可。

通过两个变换，我们不需要申请额外的空间，同时使用了非递归实现，速度更快。

快速数论变换

FFT 存在精度丢失的问题，且不能取模。快速数论变换 NTT 可以解决模意义下的变换。

阶和原根

满足 \(a^n \equiv 1 ~(mod ~p)\) 存在的最小整数 \(n\) 称为 \(a\) 模 \(m\) 的阶，\(n=ord_m(a)\)。

若 \(\gcd(g,m)=1\)，且 \(ord_m(g)=\varphi(m)\)，则 \(g\) 为模 \(m\) 的原根。

\(g,g^2, \cdots ,g^{\varphi(m)}\) 构成 \(m\) 的简化剩余系，当 \(m\) 是质数时，这些数模 \(m\) 两两不同，且根据欧拉定理，\(g^{a+\varphi(m)} \equiv g^a ~(mod ~m)\)。

NTT 数论变换

现在假定 \(m\) 是质数 \(p\)，如果 \(\varphi(p)\) 包含较多的 \(2\) 因子，可以使用 \(p\) 的原根代替复数单位根，用 \(g^\frac{p-1}{n}\) 代替 \(\omega_n^1\)，设为 \(G_n\)，有 \(G_n^{a} \equiv -G_n^{a+n/2} ~(mod~p)\)。

常见的质数有 \(998244353=7 \times 17 \times 2^{23}+1,g=3\)，能处理大约 \(8 \times 10^6\) 次的多项式。

求逆变换时，相似的用 \(\frac{1}{g}^\frac{p-1}{n}\) 当作单位根。

//inv=qpow(Bint(3),p-2);
int rev[N];
void change(Bint f[],int len)
{
    for(int i=0;i<len;i++)  
    {
        rev[i]=rev[i>>1]>>1;
        if(i&1) rev[i]|=len>>1;
    }
    for(int i=0;i<len;i++)
        if(i<rev[i]) swap(f[i],f[rev[i]]);
}

void ntt(Bint f[],int len,int op)
{
    change(f,len);
    for(int h=2;h<=len;h<<=1)
    {
        Bint step=qpow((op==1)?Bint(3):inv,(p-1)/h);
        for(int j=0;j<len;j+=h)
        {
            Bint w(1);
            for(int k=j;k<j+h/2;k++)
            {
                Bint u=f[k],v=w*f[k+h/2];
                f[k]=u+v; f[k+h/2]=u-v;
                w=w*step;
            }
        }
    }
    if(op==-1)
    {
        Bint iv=qpow(Bint(len),p-2);
        for(int i=0;i<len;i++) f[i]=f[i]*iv;
    }
}