多项式学习笔记

哎，果然，不写还是一点不会。只讲解原理，代码实现上的优化不讲。目录写完再补。

一、若干定义与基础

1.多项式和多项式的表达

多项式是若干个单项式的和，每个单项式称作多项式的项，定义多项式最高次项的次数是多项式的度，记作 \(\deg A\)，一个多项式一定可以表示成形如 \(A(x)=\sum\limits_{i=0}^{\deg A}a_i\times x^i\) 的形式，这也被叫做多项式的系数表达。

对与度为 \(n\) 的多项式，将其抽象为二维平面上的函数，其点值表达即取函数上 \(n\) 个不同的点值，写作 \(A(x)=\{(x_0,A(x_0)),(x_1,A(x_1))\dots(x_{n-1},A(x_{n-1}))\}\)。

系数表达和点值表达都可以唯一确定一个多项式，系数表达显然。对于点值表达，我们构造如下矩阵：

\[\begin{bmatrix}x_0^0&x_0^1&\cdots&x_0^n\\x_1^0&x_1^1&\cdots&x_1^n\\\cdots&\cdots&\cdots&\cdots\\x_{n-1}^0&x_{n-1}^1&\cdots&x_{n-1}^n\end{bmatrix}\times\begin{bmatrix}a_0\\a_1\\\cdots\\a_n\end{bmatrix}=\begin{bmatrix}A(x_0)\\A(x_1)\\\cdots\\A(x_{n-1})\end{bmatrix} \]

其中乘法的左矩阵与乘法的结果已知。于是原问题等价于求做矩阵的逆。

发现这样的矩阵非常特殊，事实上，我们称这样的矩阵为范德蒙德矩阵，而范德蒙德矩阵一定有逆。

因此我们知道，多项式的两种表达都可以唯一确定一个多项式，且两种表达方式可以互相转换。我们称系数到点值的过程为多点求值，点值到系数为插值。

2.向量，复数，原根

以上是多项式算法的数学基础，请熟练掌握他们的一些性质，不多赘述，可以自己去学。
复数相关，原根相关。

具体地，需要你知道：

\[\omega_n^k=\frac{2\pi}{n}\times(\cos k+i\sin k) \]

\[\omega_{2n}^{2k}=\omega_n^k,w_n^{\frac{n}{2}+k}=-\omega_n^{\frac{k}{n}} \]

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

二、多项式卷积

1.多项式乘法

对于两个多项式 \(A,B\) 定义他们的乘法运算 \(A\times B=\sum\limits_{i=0}^{\deg A}\sum\limits_{j=0}^{\deg B}A_iB_jx^{i+j}\)

直接按式子计算，显然是 \(O(n^2)\)，我们希望有更优秀的算法。于是有了 FFT。

发现点值表达的乘法是 \(O(n)\) 的，只需要对应项对应相乘即可。考虑是不是能有一种方法，把相乘的两个多项式变成点值表达，相乘得到结果多项式的点值表达后再插值回去，通过这种方式来避免直接相乘？答案是可以的。

\(\textbf{FFT, FAST FOURIER TRANSFORM}\)

先看多点求值这一步，我们来推一推式子，有很妙的一点：

\[A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \]

我们把次数为单数和双数的项分别提出来,在这里不妨设 \(n\) 是双数。

\[A1(x)=a_0+a_2x^2+\cdots+a_{n-2}x^{n-2} \]

\[A2(x)=a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1} \]

发现 \(A2\) 可以改写为：

\[A2(x)=x\times(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}) \]

于是不妨令 \(A2\) 为 \(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}\)

很容易发现 \(A1\) 和 \(A2\) 的各项次数完全相同，且每一项的次数都是偶数，因此可以进一步把他们改写成关于 \(x^2\) 的多项式，于是我们有。

\[A(x)=A1(x^2)+xA2(x^2) \]

于是我们的式子有：

\[A(x)=A1(\omega_n^{2k})+w_n^kA2(\omega_n^{2k}) \]

\[=A1(\omega_{\frac{n}{2}}^k)+\omega_n^kA2(\omega_{\frac{n}{2}}^k) \]

这个时候，如果我们取另一组 \(x=\omega_n^{k+\frac{n}{2}}\) 再带到式子里

\[A(x)=A1((\omega_n^{k+\frac{n}{2}})^2)+\omega_n^{k+\frac{n}{2}}A2((\omega_n^{k+\frac{n}{2}})^2) \]

\[=A1((-\omega_n^k)^2)-\omega_n^kA2((-\omega_n^k)^2) \]

\[=A1(\omega_n^{2k})-w_n^kA2(\omega_n^{2k}) \]

\[=A1(\omega_{\frac{n}{2}}^k)-\omega_n^kA2(\omega_{\frac{n}{2}}^k) \]

欸，你发现什么？是不是上面的式子和下面的式子就差了一个正负号？

这说明什么？是不是如果我们让 \(k\) 从 \(1\) 遍历到 \(\frac{n}{2}\)，\(k+\frac{n}{2}\) 也同样经过了 \(\frac{n}{2}\) 到 \(n\)，\(\omega_n^k\) 是可以 \(O(1)\) 求得。也就是说，我们只用求解 \(1\) 到 \(\frac{n}{2}\) 这部分的 \(A1\) 和 \(A2\) 的值，就可以得到另一半的所有结果。

然后你发现我们本质是要求 \(A1,A2\) 这两个度为 \(n/2\) 的多项式的点值，问题是等价的，但是很明显规模缩小了一半，可以递归做。

很明显递归最多 \(\log n\) 层，原多项式的每一项系数都被分到了一个要求解的多项式里，所以每层是 \(O(n)\)，我们做一遍多点求值是 \(O(n\log n)\) 的。

具体实现上，由于我们每次会将规模除二，所以如果 \(n\) 是 \(2\) 的次幂是很好算的，一般取大于等于 \(n\) 的最小幂次。

代码给一下：

inline void fft(std::vector<cp>a, int lm){
    if(lm == 1) return; //最底层时一定只有常数项
    std::vector<cp> a1, a2; //cp:复数
    for(int i = 0; i <= lm; i += 2)
        a1.pb(a[i]), a2.pb(a[i + 1]);
    fft(a1, lm >> 1), fft(a2, lm >> 1);
    cp wn = cp(cos(2.*pi / lm), sin(2.*pi) / lm), w = cp(1, 0); //定义
    for(int i = 0; i < lm >> 1; i++, w = w * wn){
        int f = a1[i], g = w * a2[i];
        a[i] = f + g, a[i + (lm >> 1)] = f - g; //直接计算
    }
}

当然有优化，但是那不是原理上的，我也不想讲，只给代码了。

int r[N];
inline void fft(std::vector<cp> a, int lm){
    for(int i = 1; i <= n; i++)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1));
    for(int i = 1; i <= n; i++)
        if(i < r[i]) std::swap(a[i], a[r[i]]);
    for(int md = 1; md < lm; md <<= 1){
        cp wn = (cos(pi/md), sin(pi/md));
        for(int i = 0; i < md; i += md << 1){
            cp w = cp(1, 0);
            for(int j = 0; j < md; j++, w = w * wn){
                cp f = a[i + j], g = w * a[i + j + md];
                a[i + j] = f + g, a[i + j + md] = f - g;
            }
        }
    }
}

好了，现在我们已经会把一个多项式 \(n\log n\) 地变成点值表达，把两个要相乘的多项式乘起来，得到了结果的点值表达，可喜可贺！

问题来了，我们要如何把点值表达再变回系数表达？

不如回顾之前对多项式基础的概述。

\[\begin{bmatrix}x_0^0&x_0^1&\cdots&x_0^n\\x_1^0&x_1^1&\cdots&x_1^n\\\cdots&\cdots&\cdots&\cdots\\x_{n-1}^0&x_{n-1}^1&\cdots&x_{n-1}^n\end{bmatrix}\times\begin{bmatrix}a_0\\a_1\\\cdots\\a_n\end{bmatrix}=\begin{bmatrix}A(x_0)\\A(x_1)\\\cdots\\A(x_{n-1})\end{bmatrix} \]

这样的矩阵，我们知道等号左边的其中两项，我们想要求出剩下的。

你可以高斯消元，然后你发现你的复杂度变成了 \(O(n^3)\)。有一种科技叫拉格朗日插值，可以 \(n\log n\) 但是我不会，以后要是会一下的话会写。

但是不用这么麻烦! 还记得我们说的吗？左边这个矩阵是一个非常特殊的矩阵，他的名字叫范德蒙德矩阵，而范德蒙德矩阵一定有逆！具体地，我们可以构造：

\[\frac{1}{n}\times\begin{bmatrix}x_0^0&x_0^{-1}&\cdots&x_0^{-n}\\x_1^0&x_1^{-1}&\cdots&x_1^{-n}\\\cdots&\cdots&\cdots&\cdots\\x_{n-1}^0&x_{n-1}^{-1}&\cdots&x_{n-1}^{-n}\end{bmatrix} \]

易于验证他的正确性。那不就好办了吗？想要求出系数表达，等于给等号左右两边同时左乘以上的逆矩阵，发现这一步相当于把我们求出来的点值当作多项式的系数，做一遍 FFT 。在这个过程中，我们的 \(\omega_n^k\) 变成了 \(-\omega_n^k\)。但是不妨看我们以上的推导。

\[A(x)=A1(\omega_n^{2k})+w_n^kA2(\omega_n^{2k}) \]

发现 \(A1,A2\) 是两个关于 \((\omega_n^k)^2\) 的多项式，夏然我们改变 \(\omega\) 的正负不影响这一步，唯一有变化的是 \(\omega_n^k\) 事实上我们在计算的时候把它变成 \(-\omega_n^k\) 就行。

所以我们得到点值表达后，可以对之再做一遍 FFT 过程中只改变单位复根，在最后乘上 \(\frac{1}{n}\) 就可以完成一趟多项式乘法。

代码不是问题，懒得写了。

\(\textbf{NTT, NUMBER THEORETIC TRANSFORM}\)

NTT 是 FFT 的数论实现，因此在学习 NTT 前请保证你会 FFT。

你发现，欸 FFT 用了复数，所以要用浮点数，但是浮点数很慢啊，而且有误差。于是我们使用另一种东西来替代复根这个东西，于是我们选择了原根。

原根的一些概念自己学去，不 dfs 式讲课。

对于原根，有以下我们要用的性质。

令模数 \(p=qn+1\) 其中 \(n\) 是 \(2\) 的次幂。

\[g^{qn}=q^{p-1}=g^{\varphi(p)}\equiv1(\bmod\ p) \]

我们令 \(g_n=g^q(\bmod\ p)\) 等价为复数中的 \(\omega_n^1\)，那么

\[g_n^i\times g_n^j=g^{iq}+g^{jq}=g^{(i+j)q}=g_n^{i+j} \]

\(g_n^1,g_n^2\cdots g_n^n\) 互不相同，构成长为 \(n\) 的循环。

\[g_n^n\equiv 1,g_n^{\frac{n}{2}}\equiv-1 \]

我们发现，在 FFT 推导中用到的一切复数性质，在原根上都能找到很好的对应。

所以我们考虑把 FFT 中的复数换成原根。在这个过程中，式子的推导不成问题，问题是如何求解 \(g_n\)

有 \(\omega_n=\omega_{p-1}^{\frac{p-1}{n}}=\omega_{\varphi(p)}^{\frac{\varphi(p)}{n}}\)（由定义）。显然这个式子的循环节是 \(\varphi(p)\) 故有 \(g=\omega_n\)

所以在 NTT 的时候，\(g_n=\omega_n=\omega_{\varphi(p)}^{\frac{\varphi(p)}{n}}=g^{\frac{p-1}{n}}\)

需要注意，以上式子有意义当且仅当 \(n\) 整除 \(p-1\)。故我们一般选择 \(998244353=119\times2^{23}+1\) 这种质数，因为它可以承载 \(2^{23}\) 这种规模的运算。它的一个原根是 \(3\)，想整活的老哥可以试试 \(114514\) 事实上它是 \(998244353\) 的原根。

代码吗，把 FFT 里的 \(\omega\) 换成我们这里的原根就行。

NTT 显然会更快，而且不存在精度问题。