OI 中的多项式(或许完工?)
OI 中的多项式
目前已有知识点:
拉格朗日插值
模板
给定 \(n\) 个点 \((x_i,y_i)\),要求求出一个多项式函数 \(f(x)\),满足 \(f(x)\) 的图像经过这 \(n\) 个点,并且 \(f(x)\) 的度为 \(n-1\)。给出一个 \(k\),求 \(f(k)(\bmod 998\ 244\ 353)\)。(from Luogu P4781)
一个很朴素的办法就是直接列出方程组,然后用高斯消元求解。这样的时间复杂度是 \(\mathcal O(n^3)\) 的,并不算一个高效的算法。而使用拉格朗日插值可以在 \(\mathcal O(n^2)\) 的时间内解决这个问题,是很高效的一种算法。
设拉格朗日的基本多项式为
那么会发现对于一个 \(l_j\),会有 \(l_j(x_j)=1\) 并且 \(l_j(x_i)=0,\forall i \neq j\)。
有了这个基本多项式,那么可以就构造出函数 \(f(x)\):
容易发现,这个函数一定是满足题目中给出的 \(n\) 个点的。将两个式子合并到一起,就可以得到答案:
由于题目让我们求 \(f(k)\) 的值,所以可以直接将上面式子中的 \(x\) 替换成 \(k\) 进行计算。
快速傅里叶变换 FFT
前置知识
复数
复数 \(z\) 是形如 \(a+bi\) 的数,其中 \(i^2=-1\),\(a\) 称为实部,\(b\) 称为虚部。
复平面是一个笛卡尔平面,有两条坐标轴,横轴是实部,纵轴是虚部。那么根据复数的形式可知,复数可以用复平面上的一个向量 \((a,b)\) 表示。以横轴正方向为始边,向量为终边,逆时针转过的角度叫做复数 \(z\) 的辐角 \(\theta\),并且 \(\theta\) 必须为正数。
复数的模 \(|z|=\sqrt{a^2+b^2}\),与向量的模相同。
复数运算
加减法
加减法很好理解,实部和虚部对应相加即可。
乘法
类似二项式乘法的计算方法。
看一下在复平面上的意义。发现 \(\theta_0=\theta_1+\theta_2\),并且 \(|z_0|=|z_1|\times|z_2|\)。
欧拉公式
单位根
复数域下,满足 \(x^n=1\) 的 \(x\) 被称作 \(n\) 次单位根。
显然这样的 \(x\) 会有 \(n\) 个。按照辐角大小排序后第 \(k\) 大的根是 \(e^{i\frac{2k\pi}{n}}\)。
可以发现,这 \(n\) 个单位根在复平面上刚好 \(n\) 等分单位圆。
本原单位根
代数上的不好理解,就不提了。复平面上辐角最小的单位根就叫做本原单位根。
将 \(n\) 次本原单位根记作 \(\omega_n\),那么有 \(\omega_n=e^{i\frac{2k\pi}n}\)。根据虚数乘法在复平面上的表示方法,可以发现第 \(k\) 大的 \(n\) 次单位根可以表示成 \(\omega_n^k\)。
关于本原单位根,有一个很重要的性质:
设 \(n = 2m\),则有:
关于这两个式子的证明,可以画出复平面过后然后用旋转的想法去证明,这里不再给出(其实是懒)。
离散傅里叶变换 DFT
考虑两个多项式做乘法,如果直接用系数表示法的话显然时间复杂度是 \(\mathcal O(n^2)\) 的,并且很难优化。而使用点值表示法的话,计算乘法的时间复杂度仅是 \(\mathcal O(n)\),但是将系数表示法转化为点值表示法的时间复杂度是 \(\mathcal O(n^2)\),并且将点值表示法转化称为系数表示法需要高斯消元,时间复杂度是 \(\mathcal O(n^3)\) 的,比直接系数表示法算慢多了。那么如果可以找到一种能够快速将系数表示法转化成点值表示法并且同样快速的转化回来的算法,那么用点值表示法计算将会快很多。
将系数表示法转化为点值表示法的变换叫做 离散傅里叶变换(DFT),将点值表示法转化为系数表示法的变换叫做 离散傅里叶逆变换(IDFT)。
将一个 \(n-1\) 次多项式表示成为 \(A(x)\),方便下面的书写。
考虑构造出一个长度为 \(n\) 的数列 \(\{b_i\}\),其中 \(b_k=A(\omega_n^k)\),即(将 \(\sum\) 内的变量写作 \(j\) 的原因是避免与虚数的 \(i\) 混淆):
这里的 \(b_k\) 就表示了多项式在 \(\omega_n^k\) 处的点值。如果直接求这一过程,时间复杂度是 \(\mathcal O(n^2)\) 的,如果使用 快速傅里叶变换(FFT) 就可以在 \(\mathcal O(n\log n)\) 的时间内求出这一点值。
快速傅里叶变换 FFT
将 \(A(x)\) 按照指数的奇偶性分成两组(\(n=2m\)):
将后面那个求和式中的最后一个 \(x\) 拿出来:
假设:
那么有:
如果知道了 \(A_0\) 和 \(A_1\) 的点值,那么就可以 \(\mathcal O(n)\) 去计算 \(A(x)\) 在各点的点值了,并且会发现求 \(A_0,A_1\) 的过程与求解 \(A\) 的方式是一样的,所以计算 \(A_0,A_1\) 可以直接递归处理。但是计算一下时间复杂度,会发现这样仍然是 \(\mathcal O(n^2)\) 的,因此考虑去优化这个做法。
在前置知识中提到了关于本原单位根的两个公式:
那么根据第二个公式,启发我们只计算 \(A(x)\) 前 \(m\) 个位置的点值,然后尝试发掘前 \(m\) 个位置和后 \(m\) 个位置的关系。
设 \(0\le k<m\),那么上面推出来的式子可以化成:
那么根据第一个式子可以将这个式子化成:
然后对于后 \(m\) 个位置,可以写成:
用第二个式子可以化成:
然后效仿前一个的推法:
会发现求 \(A(x)\) 的时候只需要求到 \(A_0,A_1\) 前 \(m\) 项的值即可,相比于之前,每次递归计算的值少了一半,所以时间复杂度变成了 \(\mathcal O(n\log n)\)。
离散傅里叶逆变换 IDFT
目前已知了一个多项式经过 DFT 后变成的点值表达式,即上面所说的 \(\{b_k\}\) 数列,也就是:
要求原系数表达式数列 \(\{a_k\}\)。
这里直接给出结论(不会证):
如果直接求解的话仍然是 \(\mathcal O(n^2)\) 的,所以考虑用 快速傅里叶逆变换(IFFT) 来进行这一操作。
快速傅里叶逆变换 IFFT
会发现求解的这个式子和 FFT 是高度一致的。因为 \(\omega_n^{-k} = \omega_n^{n+k}\),所以可以直接用 FFT,然后将计算出来的 \(\{b_k\}\) 翻转过来就是 IFFT 得到的式子。
这样,我们得到了一个可以 \(\mathcal O(n\log n)\) 计算多项式乘法的算法。用洛谷的模板题给出代码。
蝶形变换
会发现分治的过程可以用迭代的方式代替,这样就省去了动态开空间的时间。
具体做法就是先算出每个位置的值在分治到叶子节点的时候的位置,然后换到对应位置后直接从分治树的底端向上迭代进行计算。
for (int i = 1; i < lim; ++i)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (tim - 1));
代码结合例题给出。
例题
模板
给定多项式 \(F(x),G(x)\),求 \(F(x)*G(x)\)。(from Luogu P3803)
值得一说的是,这道题在之前递归做法因为常数过大是被卡了的,但是现在可以比较轻松跑过去了。
模板(高精乘)
给定 \(n,m\),求 \(n\times m\)(\(n,m\le 10^{10^6}\))。(from Luogu P1919)
显然正常的高精度乘法是 \(\mathcal O(n^2)\) 的,会直接 T 掉,发现一个数可以被化成多项式的形式,然后直接用 FFT 做就可以了。
需要注意最后可能某一位的系数会大于 \(9\),处理一下进位就行了。
快速数论变换 NTT
前置知识
原根
如果一个数 \(m\) 存在一个 \(g\),满足 \(\gcd(g,m)=1\),且 \(g\) 在模 \(m\) 下的阶等于 \(\varphi(m)\),那么称 \(g\) 是 \(m\) 的原根。
这个定义并不好理解,需要用到一些群论的知识。但是这个原根具有很好的性质:\(0\le i < p\) 的所有 \(g^i\) 在模 \(p\) 意义下两两不同。
关于原根的定义了解到这里就足够理解 NTT 了。
NTT
FFT 中的所有系数都需要用 double
来表示,因此计算速度会很慢,并且还会丢失精度。所以就有一种叫做 NTT 的算法可以更高效地解决多项式乘法的问题。
考虑将答案对一个数 \(p\) 取模,那么会发现在 FFT 中我们用到的单位根的性质与模 \(p\) 意义下的原根高度相似。
也就是说可以将 \(\omega_n^k\) 等价成为 \((g^{\frac{p-1}{n}})^k\)。
单位根的两个重要性质:\(\omega_n^{2k}=\omega_m^k\) 和 \(\omega_n^k = -\omega_n^{k + m}\) 在原根上也具有。
因此只需要将 FFT 的过程中使用到的单位根换成原根即可。
例题
给定多项式 \(F(x),G(x)\),求 \(F(x)* G(x)\)。(from Luogu P3803)
多项式ln
前置知识:
这里给出多项式求导和积分的式子:
\(C\) 可以为任意常数。
\(\ln\) 的求导法则:\((\ln x)'=\dfrac 1 x\)。
复合函数求导法则:\(f(g(x))'=f'(g(x))g'(x)\)
多项式 ln
给定多项式 \(A(x)\),求 \(B(x)\) 满足 \(B(x)\equiv \ln A(x)\pmod{x^n}\)。(from Luogu P4725)
对两侧求导,得到 \(B'(x)\equiv (\ln A(x))'\pmod{x^n}\)。展开得:\((\ln A(x))'=\dfrac {A'(x)} {A(x)}\)。因此得到 \(B'(x)\equiv \dfrac{A'(x)}{A(x)}\pmod{x^n}\)。所以使用多项式求逆即可算出 \(B'(x)\),然后对 \(B'(x)\) 积分一次就可以得到 \(B(x)\)。
多项式牛顿迭代
有多项式 \(f(x)\),且有 \(f(g(x))\equiv 0\pmod{x^n}\),求 \(\bmod\ x^n\) 意义下的 \(g(x)\)。
可以采用倍增的方式,假设已经求得了 \(\bmod\ x^{\lceil n/2\rceil}\) 意义下的 \(g_0(x)\),那么将 \(f(g(x))\) 在 \(g_0(x)\) 处泰勒展开得到:
容易发现有 \(g(x)\equiv g_0(x)\pmod{x^{\lceil n/2\rceil}}\),那么 \(g(x)-g_0(x)\) 的最低非零系数项的次数为 \(\lceil\dfrac n 2\rceil\),所以发现,当上述式子的 \(i\ge 2\) 时,\((g(x)-g_0(x))^i\) 在模 \(x^n\) 意义下恒为 \(0\)。因此这个无穷项求和的式子只有前两项非 \(0\),因此有:
推出:
多项式求逆
给定多项式 \(A(x)\),求 \(B(x)\) 满足 \(A(x)*B(x)\equiv 1\pmod{x^n}\) 。(from Luogu P4238)
转化一下就是求 \(\dfrac 1 {B(x)} - A(x)\equiv 0\pmod{x^n}\)。
设 \(\dfrac 1{B_0(x)}-A(x)\equiv 0\pmod{x^{\lceil n/2\rceil}}\)。尝试套用牛顿迭代的式子。令 \(F(B_0(x))\equiv\dfrac 1 {B_0(x)}-A(x)\pmod{x^n}\),那么有:
当 \(n=1\) 时,显然 \([x^0]B(x)=\dfrac1{[x^0]A(x)}\)。
对于时间复杂度,有递推式 \(\mathcal T(n)=\mathcal T(\dfrac n 2)+\mathcal O(n\log n)=\mathcal O(n\log n)\)。
多项式开根
给定多项式 \(A(x)\),求 \(B(x)\) 满足 \(B^2(x)\equiv A(x)\pmod{x^n}\)。(from Luogu P5205)
即 \(B^2(x)-A(x)\equiv 0\pmod{x^n}\)。
与多项式求逆做法相同,设 \(B_0^2(x)-A(x)\equiv 0\pmod{x^{\lceil x/2\rceil}}\)。那么有:
除法使用多项式求逆即可。
对于 \(n=1\) 时,有 \([x^0]B(x)=\sqrt{[x^0]A(x)}\),当不保证 \([x^0]A(x)=1\) 时需要用到二次剩余。
时间复杂度同样是 \(\mathcal O(n\log n)\),但是常数很大。
多项式 exp
给定多项式 \(A(x)\),求 \(B(x)\) 满足 \(B(x)\equiv e^{A(x)}\pmod{x^n}\)。(from Luogu P4726)
对两侧取 \(\ln\) 得到 \(\ln B(x)-A(x)\equiv 0\pmod{x^n}\)。
还是设 \(\ln B_0(x)-A(x)\equiv 0\pmod{x^{\lceil n/2\rceil}}\),有:
需要使用多项式 ln。
对于 \(n=1\),当 \([x^0]A(x)=0\) 时,有 \([x^0]B(x)=1\),否则 \([x^0]B(x)\) 模意义下不存在。
时间复杂度仍然是 \(\mathcal O(n\log n)\),常数也很大。
生成函数
前置知识: 数列
普通生成函数
定义
定义一个序列 \(a\) 的普通生成函数 OGF 为:
\(a\) 既可以是有穷序列,也可以是无穷序列。
如果 \(a\) 是有穷序列,比如 \(a=\langle1,2,3,4\rangle\)(下标统一从 \(0\) 开始),那么序列 \(a\) 的普通生成函数是 \(F(x)=1+2x+3x^2+4x^3\)。
如果 \(a\) 是无穷序列:
- \(a=\langle1,1,1,\cdots\rangle\) 的普通生成函数是 \(F(x)=\displaystyle\sum\limits_{n=0}^{\infty}x^n\)。
- \(a=\langle1,2,3,\cdots\rangle\) 的普通生成函数是 \(F(x)=\displaystyle\sum\limits_{n=0}^{\infty}(n+1)x^n\)。
- \(a=\langle1,2,4,8,\cdots\rangle\) 的普通生成函数是 \(F(x)=\displaystyle\sum\limits_{n=0}^{\infty}2^nx^n\)。
可以理解成为 \(x^n\) 前的那个式子就是这个序列的通项公式。
运算
加减法:
乘法(卷积):
其实卷积就是多项式乘法。
封闭形式
如果直接带着这么多求和号进行运算会非常麻烦,因此可以用一种特殊的办法来表示普通生成函数。
以 \(F(x)=\displaystyle\sum\limits_{n=0}^{\infty}x^n\) 为例,有方程:
解得:
这就称为 \(F(x)\) 的封闭形式。
等比数列 \(\langle1,k,k^2,k^3,\cdots\rangle\) 的普通生成函数是 \(F(x)=\displaystyle\sum\limits_{n=0}^{\infty}k^nx^n\),有:
指数生成函数
定义
序列 \(a\) 的指数生成函数 EGF 定义为 \(F(x)=\displaystyle\sum\limits_{n=0}^{\infty}a_n\dfrac{x^n}{n!}\)
运算
加减法:
乘法(卷积):
封闭形式
序列 \(a=\langle1,1,1,\cdots\rangle\) 的指数生成函数是:
序列 \(a=\langle1,k,k^2,k^3,\cdots\rangle\) 的指数生成函数是:
例题
有 \(n\) 堆糖,每一堆糖有 \(a_i\) 个,你至少需要拿走 \(l\) 个糖,并且最多拿走 \(r\) 个糖,求你有多少种拿法(所有数据小于 \(1000\))。
如果将每堆糖的方案表示成一个序列(下标表示取的个数,对应位置上的数表示这样取的方案数),那么可以得知第 \(i\) 堆糖的方案的普通生成函数为 \(F(x)=\displaystyle\sum\limits_{n=0}^{a_i}x_n\)。
由卷积的定义可以得知最后的方案数就是将每一堆糖的普通生成函数卷积起来就行了,即 \(A(x)=\displaystyle\prod\limits_{i=1}^{n}F_i(x)\)。
由于数据范围很小,因此可以直接分治乘起来,或者用 NTT 或 FFT 计算。
有 \(1,3,5,7,9\) 五个数,要求组成一个 \(n\) 位数,每位上的数字可以相等,并且要求 \(3,7\) 的个数都是偶数个,问方案数。
延续上面的思路,用指数生成函数表示每一个数的方案:
- \(1,5,9\) 为 \(F(x)=\displaystyle\sum\limits_{i=0}^{\infty}\dfrac{x^i}{i!}=e^x\)
- \(3,7\) 为 \(F(x)=\displaystyle\sum\limits_{i=0}^{\infty}\dfrac{x^{2i}}{{2i!}}=\dfrac{e^x+e^{-x}}{2}\)
那么总的方案的指数生成函数就是这五个值卷起来的结果:\(A(x)=(e^x)^3\cdot(\dfrac{e^x+e^{-x}}2)^2=\dfrac{e^{5x}+2e^{3x}+e^x}4\)。
将这个式子从封闭形式转回来得到通项公式:\(A_n=\dfrac{5^n+2\times3^n+1}4\)。