前置知识

复数

\(i^2=-1\)
复数为 \(z=a+bi\)
\(a+bi,a-bi\) 互为共轭复数
一个复数对应复平面上的任意一个点，一个复数与复平面上 \((a,b)\) 一一对应。复数与向量之间的运算有相似的地方，假设复数到复平面原点的距离为 \(r\) ，那么复数也可以用 \((r\cos \alpha,r\sin \alpha)\) 来表示。
加法为向量加法：\(z_1+z_2=(a+c)+(b+d)i\)
乘法：模长相乘幅角相加。\((ac-bd)+(ad+bc)i\)

也可以表示成 \((r_1r_2\cos(\theta_1+\theta_2),r_1r_2\sin(\theta_1+\theta_2))\)。实质是向量旋转。

单位根

\(n\) 次单位根为 \(x^n=1\) 的 \(n\) 个根。

\(\omega_n=\cos\frac{2\pi}{n}+\sin\frac{2\pi}{n}i\)

所谓的 \(n\) 个根为 \(\omega_n^1,\omega_n^2...\omega_n^n\)。

注意到：

\[\begin{aligned} \omega_n^2&=\cos\frac{2\pi}{n}\times \cos\frac{2\pi}{n}+\sin\frac{2\pi}{n}\times \sin\frac{2\pi}{n}+(\sin\frac{2\pi}{n}\times \cos\frac{2\pi}{n}+\cos\frac{2\pi}{n}\times \sin\frac{2\pi}{n})i\\ &=\cos \frac{2\pi\times2}{n}+\sin\frac{2\pi\times 2}{n}i \end{aligned} \]

推广后得到：

\[\omega_{n}^k=\cos\frac{2k\pi}{n}+\sin\frac{2k\pi}{n}i \]

单位根满足以下性质：

\(\omega_n^k=\omega_{2n}^{2k}\)
证明：\(\omega_n^k=(\cos\frac{2k\pi}{n},\sin\frac{2k\pi}{n})=(\cos\frac{4k\pi}{2n},\sin\frac{4k\pi}{2n})=\omega_{2n}^{2k}\)
\(\omega_n^{k+\frac n2}=-\omega_{n}^k\)
\(\omega_n^0=\omega_n^n=1\)

有关于上面若干性质的证明都比较显然，这里不再赘述。

多项式

原始表示法：\(f(x)=\sum\limits_{i=0}^na_ix^i\)
原始的多项式乘法：\(\sum\limits_{i=0}^{n_1}\sum\limits_{j=0}^{n_2}a_ib_jx^{i+j}\) 。
点值表示法：\(n+1\) 个点确定一个 \(n\) 次多项式，所以我们可以用 \((x_1,f(x_1)),(x_2,f(x_2)),...(x_{n+1},f(x_{n+1}))\) ，来表示一个 \(n\) 次多项式。
点值表示法多项式相乘：

不难发现，多项式相乘后的点值为 \((x_1,g(x_1)f(x_1)),(x_2,g(x_2)f(x_2)...x_{n+1}g(x_{n+1})f(x_{n+1}))\)

多项式乘法

容易发现原式的多项式乘法复杂度是 \(O(n^2)\) 而点值表示法师 \(O(n)\)，所以我们考虑原始表示法和点值表示法的相互转化，我们称 \(DFT(f)\) 为原始表示法 \(f\) 化为点值表示法，\(IDFT(f)\) 表示点值表示法 \(f\) 化为原始表示法。那么我们进行多项式乘法需要完成的步骤有三个：

\(DFT(f),DFT(g)\)
\(DFT(f\times g)=DFT(f)\times DFT(g)\)
\(IDFT(DFT(f\times g))\)

其中，第二步可以 \(O(n)\)。现在只需要解决第一个步骤和第三个步骤。

DFT

以下我们假设多项式项数 \(n\) 为 \(2^k\) 我们可以强行补到 \(2^k\) 次方的形式。令多项式 \(A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\)

我们选取 \(n\) 个点来使上面的多项式变为点值表示法，而我们选择的 \(n\) 个点（此时多项式次数为 \(n-1\) ）为 \(\omega_n^0,\omega_n^1...\omega_n^{n-1}\)。

我们令：

\[\begin{aligned} A_1(x)=a_0+a_2x+...a_{n-2}x^{\frac n 2-1}\\ A_2(x)=a_1+a_3x+...a_{n-1}x^{\frac n2-1} \end{aligned} \]

不难发现，这时有：

\[A(x)=A_1(x^2)+xA_2(x^2) \]

我们令 \(x=\omega_n^k\)

发现 \(A(\omega_{n}^k)=A_1(\omega_n^{2k})+\omega_{n}^kA_2(\omega_{n}^{2k})=A_1(\omega_{\frac{n}{2}}^k)+\omega_{n}^kA_2(\omega_{\frac{n}{2}}^k)\)

我们还有 \(A(\omega_n^{k+\frac n2})=A_1(\omega_{\frac n2}^k)-\omega_n^kA_2(\omega_{\frac n2}^k)\)

所以复杂度应该为 \(T(n)=2T(\frac{n}{2})+O(n)=O(n\log n)\)。

IDFT

我们令 \(B=DFT(A)\) 即 \(b_i=A(\omega_n^i)\) ，我们设 \(C(x)=\sum\limits_{i=0}^{n-1} b_ix^i\)，然后让 \(x=\omega_n^{-k}\)

我们就有了下面这个式子：

\[\begin{aligned} C(\omega_n^{-k})&=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(\omega_{n}^i)^j(\omega_n^{-k})^i\\ &=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_{n}^{j-k})^i \end{aligned} \]

令 \(S(i)=\sum\limits_{i=0}^{n-1}(\omega_{n}^{j-k})^i\)，然后我们讨论一下会发现：

\(j=k\Rightarrow S(i)=n\)
\(j\not =k\Rightarrow S(i)=0\)

\(\text{Update on 2022.6.18}\)

这个性质的证明已经忘了好几次了，简单证明一下：
我们可以用等比数列拆开，可以得到：\(\frac{\omega_{n}^{(j-k)n}-1}{\omega_n^{(j-k)}-1}\)，分子是 \(0\)。\(\text{Q.E.D}\)

所以有：

\[C(\omega_n^{-k})=a_k\times n \]

所以我们直接按照 DFT 那样去做 IDFT 即可。

NTT 快速数论变换

NTT 适用于在模意义下做多项式乘法。

在 FFT 中我们用到了的性质：

\(\omega_n^n=1\)
\(\omega_{n}^1,\omega_n^2...\omega_n^n\) 互不相等
\(\omega_n^k=\omega_{2n}^{2k},\omega_{n}^{k+\frac{n}{2}}=-\omega_n^k\)
\(\sum\limits_{i=0}^{n-1}(\omega_n^k)^i=n\times [k=0]\)

假设 NTT 是在 \(\bmod\ p\) 意义下的，其中 \(p\) 为一个质数。设 \(g\) 为 \(p\) 的原根。

设 \(p-1=q\times 2^r=q\times n\) 我们令 \(n\) 为多项式项数。

设 \(\omega_n=g^q\)，我们来证明这个东西仍然满足上面的性质。

\(\omega_{n}^n=g^{qn}=g^{p-1}\equiv 1(\bmod p)\)
假设有两个相等，则有 \(g^{qi}\equiv g^{qj}(\bmod p)\Rightarrow g^{qi-qj}\equiv 1(\bmod p)\)，显然有 \(i-j\not = n\)。所以如果成立，与原根定义不符合，所以成立。
前半段：\(g^{qk}=g^{\frac{q}{2}\times 2k}\) 。所以前半段成立。

后半段：只需要证明：\(\omega_n^{\frac n2}=-1\) 。而 \(\omega_n^n=1,(\omega_n^{\frac{n}{2}})^2=1,\omega_n^{\frac{n}{2}}\not = 1\) 。所以成立。
这个性质依赖于前面的性质，显然成立。

所以我们找到了一个单位根的替代品。主要，在实践中，我们通常取模数为 \(998244353=2^{23}\times 7\times 17+1\)。这样，只要这个多项式长度不超过 \(2^{23}\) ，都可以用 NTT 而不是 FFT。

另一个常用的模数为 \(1004535809\)，里面有一个 \(2^20\) 的因子，其原根也为 \(3\)。

从线性变换角度理解 FFT

我们可以用矩阵来表示从系数表示法到点值表示法的转换，设 \(a_i\) 为 \(x^i\) 的系数，则经过上面的讨论，变换矩阵应该为：

\[\begin{pmatrix} 1&1&1&\cdots&1\\ \omega_{n}^1&\omega_n^2&\omega_n^3&\cdots&\omega_n^{n-1}\\ \omega_{n}^{2}&\omega_n^{4}&\omega_n^6&\cdots&\omega_n^{(n-1)2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \omega_n^{n-1}&\omega_n^{2(n-1)}&\omega_n^{3(n-1)}&\cdots&\omega_{n}^{(n-1)(n-1)} \end{pmatrix} \]

通过右乘系数矩阵 \((a_0,a_1,\cdots,a_{n-1})^T\) 来得到点值矩阵。

蝴蝶变换

在具体实现中，递归版本的 FFT 和 NTT 的常数过大，蝴蝶变换可以帮助我们进行非递归版的 FFT 或 NTT。

蝴蝶变换的定义很简单，就是把一个二进制数 \(x\) 在给定二进制位数 \(K\) 的前提下进行翻转得到的结果，可以简单的通过 \(O(n)\) 求出。

事实上，FFT 做完之后原来的下标到现在的位置就是原位置经过蝴蝶变换得到的结果。

应该说，这个不难证明。只要各位读者愿意把多项式长度 \(n=8\) 时候的 FFT 过程画出来就不难验证结论的正确性，至于证明，读者可以细心观察二进制位最后为 \(0\) 和为 \(1\) 时的变换结果。

昨晚蝴蝶变换之后就便于我们做 FFT 而不用复制数组。代码如下：

inline void NTT(int *f,int len,int flag){
    for(int i=0;i<len;i++) if(i<tr[i]) swap(f[i],f[tr[i]]);
    for(int i=2;i<=len;i<<=1){
        int tg=ksm(g,(mod-1)/i,mod),l=i>>1;
        if(flag==-1) tg=ksm(tg,mod-2,mod);
        for(int j=0;j<len;j+=i){
            int buf=1;
            for(int k=j;k<j+l;k++){
                int tt=1ll*buf*f[k+l]%mod;
                f[k+l]=((f[k]-tt)%mod+mod)%mod;
                f[k]=(f[k]+tt)%mod;buf=1ll*buf*tg%mod;
            }
        }
    }
    if(flag==1) return;
    int inv=ksm(len,mod-2,mod);
    for(int i=0;i<len;i++) f[i]=1ll*f[i]*inv%mod;
}

变换完之后我们就可以再也不用关注下标的问题了，根据 DFT 的式子，我们只需要知道 \(A_1,A_2\) 的对应点值即可。

多项式算法简单扩展

多项式求逆

我们考虑给定 \(A(x)\) 要求一个 \(B(x)\)，使得 \(A(x)B(x)\equiv 1 (\bmod x^n)\)

如果 \(n=1\) 显然这个问题就是一个求逆元的事情。这启发我们尝试把答案逐渐缩小。

这里假设我们已经找到了一个 \(C(x)\) 使得 \(A(x)C(x)\equiv 1\bmod x^{\lceil\frac x2\rceil}\)

我们尝试用 \(A(x),C(x)\) 去表示 \(B(x)\) 。

不难发现，有：

\[\begin{aligned} &A(x)B(x)\equiv 1\bmod x^{\lceil\frac x2\rceil}\Rightarrow A(x)(B(x)-C(x))\equiv 0 \bmod x^{\lceil\frac x2\rceil}\\ &\Rightarrow B(x)\equiv C(x)\bmod x^{\lceil\frac x2\rceil}\\ &\Rightarrow x^{\lceil\frac x2\rceil}|B(x)-C(x)\\ &\Rightarrow x^n|(B(x)-C(x))^2\Rightarrow B^2(x)-2B(x)C(x)+C^2(x)\equiv 0\bmod x^n\\ &\Rightarrow A(x)B^2(x)-2A(x)B(x)C(x)+A(x)C^2(x)\equiv 0\bmod x^n\\ &\Rightarrow B(x)-2C(x)+A(x)C^2(x)\equiv 0 \bmod x^n\\ &\Rightarrow B(x)\equiv 2C(x)-A(x)C^2(x)\bmod x^n \end{aligned} \]

所以我们就可以用多项式乘法来做了。

时间复杂度上限为 \(T(n)=T(\frac{n}{2})+n\log n=O(n\log n)\)。

多项式开根

给定 \(A(x)\) 找到一个 \(B(x)\) 使得 \(B^2(x)\equiv A(x)\bmod x^n\)

我们同样考虑缩小这个问题，假设我们已经求到了一个 \(C(x),s.t.C^2(x)\equiv A(x)\bmod x^{\lceil\frac x2\rceil}\)

那么我们有：

\[\begin{aligned} &(C^2(x)-A(x))^2\equiv 0\bmod x^{n}\\ &\Rightarrow C^4(x)-2C^2(x)A(x)+A^2(x)\equiv 0\bmod x^{n}\\ &\Rightarrow C^4(x)+2C^2(x)A(x)+A^2(x)\equiv 4C^2(x)A(x)\bmod x^n\\ &\Rightarrow (\frac{C^2(x)+A(x)}{2C(x)})^2\equiv A(x)\bmod x^n\\ &\Rightarrow B(x)=\frac{C^2(x)+A(x)}{2C(x)} \end{aligned} \]

运用多项式求逆就可以做了。

多项式除法

给定 \(A(x),B(x)\) 要求找到 \(C(x),D(x)\) 使得 \(A(x)=C(x)B(x)+D(x)\)。其中需要满足 \(D(x)\) 的次数小于 \(B(x)\) 的次数。

我们设 \(A\) 的次数为 \(n\) ，\(B\) 的次数为 \(m\) ，那么 \(C\) 的次数为 \(n-m\) ，可以认为 \(D\) 的次数为 \(m-1\)。

我们设 \(inv(A)=x^n\times A(\frac 1x)\) ，那么显然有 \(A=inv(inv(A))\) 。通俗来讲，\(inv\) 就是把系数倒着写一遍。

那么有：

\[\begin{aligned} &inv(A)=inv(C\times B+D)\\ &\Rightarrow x^n\times A(\frac 1x)=x^{n-m}\times C(\frac 1x)\times x^m\times B(\frac 1x)+x^{n-m+1}\times x^{m-1}\times D(\frac 1x)\\ &\Rightarrow inv(A)=inv(C)\times inv(B)+x^{n-m+1}inv(D)\\ &\Rightarrow inv(A)\equiv inv(C)\times inv(B)\bmod x^{n-m+1}\\ \end{aligned} \]

我们可以利用多项式求逆来求解一个 \(C\) ，然后代入原式求 \(D\) 即可。

多项式求对数

给定 \(f(x)\) 要求找到 \(g(x),s.t.g(x)\equiv \ln f(x)\bmod x^n\)

我们首先关注这里的对数的定义是什么，我们对两边求导可以得到：

\[g'(x)\equiv \frac{f'(x)}{f(x)}\bmod x^n \]

很明显等式右边是可以算的，那么我们算出后用不定积分算一下就可以得到 \(g\)。

一般来说，题目保证常数项等于 \(1\)。

多项式 exp

给定 \(f(x)\) 要求求解 \(g(x),s.t.g(x)\equiv e^{f(x)}\bmod x^n\)

首先引入 \(f(x)\) 的泰勒级数：

\[f(x)=\sum\limits_{i=0}^{\infty}\frac{f^{(i)}(x_0)}{i!}\times (x-x_0)^i \]

我们还是首先要考虑这里 exp 的定义是什么，我们对两边求对可以得到：

\[\ln g(x)\equiv f(x)\bmod x^n\Rightarrow \ln g(x)-f(x)\equiv 0\bmod x^n \]

我们令 \(h(g(x))=\ln g(x)-f(x)\) ，注意这里函数 \(h\) 的自变量是多项式 \(g(x)\) 而不是 \(x\) ，也就是说，这并不是一个复合函数，这里的 \(f(x)\) 相当于一个常量。而我们期望得到一个 \(h(g(x))=0\) 的一个解。

我们对 \(h\) 多次求导可以得到：

\[\begin{aligned} h^{(1)}(g(x))&=\frac{1}{g(x)}\\ h^{(2)}(g(x))&=-\frac{1}{g^2(x)}\\ h^{(3)}(g(x))&=2\frac{1}{g^{3}(x)}\\ h^{(n)}(g(x))&=(n-1)!\times (-1)^{n-1}\times \frac{1}{g^n(x)} \end{aligned} \]

我们还是和以前一样，设 \(\ln g_0(x)\equiv f(x)\bmod x^{\lceil\frac n2\rceil}\)，假设这个 \(g_0(x)\) 是已知的。我们现在要求 \(g_1(x),s.t.\ln g_1(x)\equiv f(x)\bmod x^n\)。我们写出 \(h(g_1(x))\) 当 \(x_0=g_0(x)\) 的泰勒级数。有：

\[h(g_1(x))=\sum\limits_{i=0}^{\infty}\frac{h^{(i)}(g_0(x))}{i!}\times (g_1(x)-g_0(x))^i \]

两边都对 \(x^n\) 取模可以得到：

\[h(g_1(x))\equiv \sum\limits_{i=0}^{\infty}\frac{h^{(i)}(g_0(x))}{i!}\times (g_1(x)-g_0(x))^i\bmod x^n \]

注意到实际上有：

\[g_0(x)\equiv g_1(x)\bmod x^{\lceil\frac n2\rceil}\Rightarrow (g_0(x)-g_1(x))^2\equiv 0\bmod x^n\Rightarrow (g_0(x)-g_1(x))^i\equiv 0\bmod x^n(i\ge 2,i\in \mathbb{N}) \]

所以有：

\[\begin{aligned} h(g_1(x))&\equiv \sum\limits_{i=0}^{\infty}\frac{h^{(i)}(g_0(x))}{i!}\times (g_1(x)-g_0(x))^i\bmod x^n\\ &\equiv \sum\limits_{i=0}^{1}\frac{h^{(i)}(g_0(x))}{i!}\times (g_1(x)-g_0(x))^i\bmod x^n\\ &\equiv h(g_0(x))+h'(g_0(x))\times (g_1(x)-g_0(x))\bmod x^n \end{aligned} \]

而 \(h(g_1(x))\equiv 0\) ：

\[\begin{aligned} &h(g_0(x))+ h'(g_0(x))\times (g_1(x)-g_0(x))\equiv 0\bmod x^n\\ &\Rightarrow \ln g_0(x)-f(x)+\frac{1}{g_0(x)}\times (g_1(x)-g_0(x))\equiv 0 \bmod x^n\\ &\Rightarrow g_1(x)\equiv g_0(x)(f(x)-\ln g_0(x))+g_0(x)\equiv 0 \bmod x^n \end{aligned} \]

所以 \(g_1(x)\) 就可以算了。