多项式与生成函数

多项式与生成函数

多项式

有关多项式

定义

在数学中，几个单项式的和，叫做多项式 。多项式中的每个单项式叫做多项式的项，这些单项式中的最高项次数，就是这个多项式的次数。其中多项式中不含字母的项叫做常数项。——baidu

假设有 $n$ 次多项式 $A$ ，其可以表示为：

$$A = \sum_{i = 0} ^ n a_i x^i$$

其中 $a_i$ 表示第 $i$ 此项的系数。

运算

多项式也有其四则运算，有 $n$ 次多项式 $A$ 和 $m$ 次多项式 $B$，

加法

$$A + B = \sum_{i = 0} ^ n a_i x^i + \sum_{i = 0}^m b_i x^i = \sum_{i = 0} ^ {\max(n, m)} (a_i +b_i)x^i$$

新的多项式次数为 $max(n,m)$。

减法

和加法相同，

$$A - B = \sum_{i = 0} ^ n a_i x^i - \sum_{i = 0}^m b_i x^i = \sum_{i = 0} ^ {\max(n, m)} (a_i -b_i)x^i$$

乘法

$$A \ast B = \sum_{i = 0} ^ {n + m}(\sum_{j + k = i} a_j b_k) x ^i$$

新的多项式次数为 $n + m$，直接展开复杂度为 $O(n^2)$，有更优的计算方式。

除法

详见多项式除法。

拉格朗日插值定理

有$n + 1$个点值即可确定一个$n$次多项式。

假设我们有 $n + 1$ 个点值为 $(x_i, y_i)$。我们可以构造出对应的多项式为，

$$f(x) = \sum_{i = 0} ^ n y_i \frac{\prod _{j \not = i}(x - x_j)}{\prod_{j \not = i}(x_i - x_j)}$$

复杂度为 $O(n^2)$，当点值连续时可以通过预处理阶乘优化到 $O(n)$。

此时，

$$f(x) = \sum_{i = 0} ^ n y_i \frac{\prod_{j = 0}^n(x - j)}{(x-i)i!(n-i)!(-1)^{n-i}}$$

不连续时有更快的多项式快速插值算法。

秦九昭算法

秦九昭算法可以用来进行单次 $O(n)$的多项式求值，对于 $n$ 次多项式 $A = \sum_{i= 0} ^ n a_i x^i$。

将次式一步步展开为，

$$\begin{aligned} A &= a_n x^n +a_{n- 1} x ^ {n - 1} + \cdots a_2 x^2 +a_1 x + a_0 \\ &= x (a_n x^{n - 1} + \cdots a_2 x + a_1) + a_0 \\ &= \cdots \end{aligned}$$

最后拆成了多个一次式，直接计算即可。但是对于多点求值的时候需要更优的算法。

导数

导数(Derivative)，也称为导函数值，又名微商。在数学中发挥着很重要的作用。

导数是函数的局部性质。一个函数在某一点的导数描述了这个函数在这一点附近的变化率。如果函数的自变量和取值都是实数的话，函数在某一点的导数就是该函数所代表的曲线在这一点上的切线斜率。导数的本质是通过极限的概念对函数进行局部的线性逼近。 ——––baidu

其实通俗的来讲导数就是一个函数的变化率，比较经典的就是二次函数求切线。

假设我们要求某一光滑曲线在 $M$ 点的切线，先在函数图像上取一点 $N$，假设直线 $MT$ 为 $M$ 点的切线，当 $\angle {NMT}$ 无限小时此时的直线 $MN$ 就已经成为切线，设直线 $MN$ 的斜率为 $\tan \phi$，有：

$$\tan \phi = \frac{f(x) - f(x_0)}{x-x_0}$$

当$N$逼近与$M$时，

$$\tan \phi = \lim\limits_{x \to x_0} \frac{f(x) - f(x_0)}{x-x_0}$$

令 $\Delta x = x-x_0$，所以

$$\tan \phi = \lim\limits_{x \to 0} \frac{f(x_0+\Delta x) - f(x_0)}{\Delta x}$$

此处 $\tan \phi$ 即为函数在 $x_0$ 处导数的取值，现在引出导数真正的定义 。

假设函数 $y = f(x)$ 在点 $x_0$ 处的邻域有定义，当自变量 $x$ 在 $x_0$ 处取得增量 $\Delta x$(仍然在邻域内)，相应的函数取得增量 $\Delta y = f(x_0+\Delta x) - f(x_0)$，如果 $\frac{\Delta y}{\Delta x}$ 在 $\Delta x \to 0$ 时的极限存在，称为函数 $y = f(x)$ 在$x_0$处可导，它的导数可以写成

$$f'(x_0) = \lim\limits_{x \to 0} \frac{f(x_0+\Delta x) - f(x_0)}{\Delta x}$$

$f'(x_0)$ 也可以记成 $\frac{\mathrm{d}y}{\mathrm{d}x}$ 或者 $\frac{\mathrm{d}f(x)}{\mathrm{d}x}$,另外函数在 $x_0$ 处可导的充分必要条件是，函数在 $x_0$ 处的左右两侧的导数都 必须存在，并且相等 。

常见函数的导数：

1. $f(x) = C$，$C$为常数，$f'(x) = 0$
2. $f(x) = x^n, f'(x) = n x ^ {n - 1}$
3. $f(x) = \sin x, f'(x) = \cos x$
4. $f(x) = \cos x, f'(x) = -\sin x$
5. $f(x) = a^x , f'(x) = a^x \ln a$
6. $f(x) = \log_ax , f'(x) = \frac{1}{x \ln a}$
7. $f(x) = \ln x, f'(x) = \frac{1}{x}$
8. $f[g(x)]'= f'[g(x)]g'(x)$
9. $[\frac{f(x)}{g(x)}]' = \frac{f'(x)g(x) - f(x)g'(x)}{g(x) ^ 2}$
10. $(f(x) \pm g(x))' = f'(x) \pm g'(x)$
11. $(f(x)g(x))' = f'(x)g(x) + f(x)g'(x)$

微分

首先思考一个简答的问题，一个正方形金属薄片边长为 $x$，受温度变化影响边长增加了 $\Delta x$，求其面积变化了多少。

设正方形面积为 $A$ ,这样我们就得到了一个函数 $A = x^2$，计算 $\Delta A$ 为:

$$\Delta A = (x_0+ \Delta x)^2 - x_0^2 = 2x_0\Delta x + (\Delta x)^2$$

首先第一部分 $2 x_0 \Delta x$ 为一个线性函数，第二部分中当 $\Delta x \to 0$ 时，$(\Delta x)^2$是比
$\Delta x$ 高阶的无穷小，即 $(\Delta x)^2 =\omicron(\Delta x)$

好了增量 $\Delta y$ 现在已经可以表示为：

$$\Delta y = A \Delta x + \omicron(\Delta x )$$

其中 $A$ 是不依赖 $\Delta x$ 的常数，并且 $\Delta y$ 与 $A\Delta x$ 的差为：

$$\Delta y - A\Delta x =\omicron(\Delta x)$$

是比 $\Delta x$ 更高阶的无穷小，这样当 $A \ne 0$，且 $|\Delta x |$ 无穷小时，$\Delta y$ 就可以近似的表示为 $\Delta y = A \Delta x$。

接下来给出定义 ，设函数 $y = f(x)$ 在某区间内有定义，$x_0$ 及 $x_0+ \Delta x$ 在这段区间内，如果函数的增量

$$\Delta y = f(x_0 + \Delta x) - f(x_0)$$

可以表示为

$$\Delta y = A \Delta x + \omicron(\Delta x)$$

其中 $A$ 是不依赖于 $\Delta x$ 常数，那么称函数 $y = f(x)$ 在 $x_0$ 处是可微的，而 $A\Delta x$ 叫做函数 $y = f(x)$ 在 $x_0$ 处相当于自变量增量 $\Delta x$ 的 微分(Differential) 记作 $dy$ ，即

$$\mathrm{d}y = A \Delta x$$

积分

积分(integral) 是微积分学与数学分析的一个核心概念。通常分为定积分和不定积分两种。（主要区别就是定积分得到的结果是一个数，不定积分得到的是函数）

tips：不定积分是求导运算的逆运算。 这一结论被称为微积分基本定理 (fundamental theorem of calculus)。

直观的说，对于一个给定的正实数值函数，在一个实数区间上的定积分可以理解为在坐标平面上，由曲线，直线以及轴围围成的 曲边梯形的面积值（一种确定的实数值）

如果一个函数的积分存在，并且有限，就说这个函数是可积的。一般来说，被积函数不一定只有一个变量，积分域也可以是不同维度的空间，甚至是没有直观几何意义的抽象空间。对于只有一个变量 $x$ 的实值函数 $f$ 在闭区间 $[a,b]$ 上的积分记作：

$$\int _{a}^{b} f(x) \mathrm{d}x$$

其中的 $\mathrm{d}x$ 就是积分变量。

黎曼积分(Riemann Integral)，就是所说的正常积分，定积分。其求积分值的核心思想就是通过无限逼近来确定这个积分值，需要注意的是，如果 $f(x)$ 去负值，则对应的面积值也为负值。

对于一个闭区间 $[a,b]$ 的分割$P$就是指在这一个区间中取一个有限的点列 $a = x_0 < x_1 < x_2 < \cdots < x_n = b$ 每个子区间长度的最大值定义为 $\lambda = \max(x_{i+1}-x_i)$。

定义取样分割为在进行分割 $P$ 后在每一个子区间 $[x_{i},x_{i+1}]$ 中取出一点 $x_i \le t_i \le x_{i+1}$,其区间长度最大值仍用 $\lambda$ 来表示。

对于一个在区间 $[a,b]$ 有定义的实值函数 $f$， $f$ 关于取样分割 $x_0 \cdots x_n , t_0 \cdots t_{n-1}$ 的黎曼和表示为：

$$\sum_{i=0}^{n-1}f(t_i)(x_{i+1}-x_i)$$

上述式子中的每一项是子区间长度 $x_{i+1}- x_i$ 与在 $t_i$ 处的 $f(t_i)$ 的乘积，其实直观的来说就是将一个曲线梯形分割成无限个小的图形将其面积相加。

看来需要一个更严格的定义，我们需要把 $\lambda$ 趋近于0函数值才能更精确。

设 $S$ 是函数 $f$ 在闭区间 $[a,b]$ 上的黎曼积分，当且仅当 $\forall \epsilon > 0, \exists \delta > 0$，使得 $\forall x_0 \cdots x_n , t_0 \cdots t_{n-1}$，只要它的子区间长度最大值 $\lambda \le \delta$ 就有：

$$\bigg| \sum_{i = 0}^{n-1} f(t_i)(x_{i+1} -x_i)-S \bigg| < \epsilon$$

也就是说，对于一个函数 $f$，如果在闭区间 $[a,b]$ 上无论如何取样分割，只要它的子区间长度最大值足够小，函数 $f$ 的黎曼和都会趋向一个确定的值，那么在闭区间 $[a,b]$ 上的黎曼积分存在，并且定义为黎曼和的极限，这时候我们称函数 $f$ 为黎曼可积的。

还有一个非常重要的定理， 若$F(x)$为$f(x)$的不定积分 ，则$f(x)$为$F(x)$ 的导数；反之亦然。

根据以上定理，我们有：

$$\int_a^b f(x) \mathrm{d}x = F(b)-F(a)$$

常见积分公式：

1. $\int k \mathrm{d}x = kx+c$，$k$为常数
2. $\int x^a \mathrm{d} x = \frac{1}{a + 1} x ^ { a + 1} + c (a \not = -1)$
3. $\int \frac{1}{x} \mathrm{d}x = \ln |x| + c$
4. $\int a^x \mathrm{d}x = \frac{a^x}{\ln a} + c (a > 0, a \not = 1)$
5. $\int \log_a x \mathrm{d} x = \frac{1}{\ln_a} (x \ln x - x) + c$

牛顿迭代

牛顿迭代可以用来求方程的近似的一个或者多个解，首先设方程函数 $f(x) = m$ 可以将方程改为 $g(x) = f(x) - m = 0$ 我们只需要求出 $g(x) = 0$ 的解，就可以求出 $f(x) = m$ 的解。

对于求解 $f(x) = 0$， 设 $r$ 是所求零点， 选取 $x_0$ 作为 $r$ 的初始近似值，则我们可以过点 $(x_0, f(x_0))$ 做曲线 $y = f(x)$ 的切线 $L$， 我们知道切线与 $x$ 轴有交点， 我们已知切线 $L$ 的方程为 $y =f(x_0) + f'(x_0)(x - x_0)$ 此时$L$ 与 $x$ 轴的交点为 $x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$ 然后再以 $(x_1, f(x_1))$ 点作切线， 求出与 $x$ 轴的交点， 重复以上过程直到 $f(x_n)$ 无限接近与 $0$， 即可，其中$n$ 次迭代公式为

$$x_{n + 1} = x_n - \frac{f(x_n)}{f'(x_n)}$$

泰勒展开

泰勒公式，也称泰勒展开式。是用一个函数在某点的信息，描述其附近取值的公式。如果函数足够平滑，在已知函数在某一点的各阶导数值的情况下，泰勒公式可以利用这些导数值来做系数，构建一个多项式近似函数，求得在这一点的邻域中的值。

就是用一个多项式函数去逼近一个给定的函数，但是必须逼近的时候从某个函数图像上的某个点展开。

设 $n$ 是一个正整数。如果定一个包含 $a$ 的区间上的函数 $f$ 在 $a$ 处 $n + 1$ 次可导， 那么对于这个区间上的任意 $x$ 都有，

$$f(x) = \sum_{i = 0} ^ n \frac{f^{(i)}(a)}{i!}(x - a) ^ i + R_i(x)$$

其中的多项式称为函数在 $a$ 处的泰勒展开式， $R_i(x)$ 是泰勒公式的余项是 $(x - a) ^ i$ 的高阶无穷小。

若 $a = 0$ 即可得到 $f(x) = f(x) = \sum_{i = 0} ^ n \frac{f^{(i)}(0)}{i!}x ^ i + R_i(x)$

多项式相关算法

快速傅里叶变换(FFT)

快速傅里叶变换，主要用于加速多项式乘法，对于两个多项式 $A$ 和 $B$，FFT可以将 $O(n^2)$ 的朴素算法优化为 $O(n \log n)$。

单位根

以单位圆点为起点，单位园的 $n$ 等分点为终点，在单位圆上可以得到 $n$ 个复数，设幅角为正且最小的复数为 $\omega_n$，称为 $n$ 次单位根，即

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

有欧拉公式

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

特别的有 $\omega_n ^ 0 = \omega_n ^ n = 1$。

另外还有以下性质

1. $\omega _ n ^ k = e ^{\frac{2 \pi i k}{n}}$
2. $\omega_{dn} ^ {dk} = \omega_n ^k$
3. $\omega_n ^k = a + bi, \omega _ n ^ {-k} = a - bi$
4. $\omega_n ^{k +\frac{n}{2}} = - \omega_n ^ k$

离散傅里叶变换(DFT)

离散傅里叶变换(DFT) 主要是利用分支思想根据一个 $n$ 次的多项式可以由 $n + 1$ 个点唯一确定

对于一个 $n$ 次多项式，

$$A = \sum_{i = 0} ^ n a_i x ^ i$$

其系数进行奇偶分类，得到

$$A_0 = a_0 +a_2 x + a_4 x^2 + \cdots \\ A_1 = a_1 +a_3 x + a_5 x^2 + \cdots \\$$

原来的式子变成了

$$A(x) = A_0(x ^ 2) + x \cdot A_1(x ^ 2)$$

然后将 $\omega_ n ^k$ 和 $\omega_n ^ {k + \frac{n}{2}}$ 带入得

$$A(\omega_n ^ k) = A_0(\omega_n ^ {2k}) + \omega_n ^ k \cdot A_1(\omega_n ^ {2k}) \\ A(\omega_n ^ {k + \frac{n}{2}}) = A_0(\omega_n ^ {2k}) - \omega_n ^ {k} \cdot A_1(\omega_n ^ {2k})$$

可以发现两个式子只有常数不一样，也就是说只要知道了，$[0, k - 1]$ 内的点值就可以得到 $[\frac{n}{2}, k +\frac{n}{2} - 1]$ 之间的点值。

时间复杂度为 $O(n \log n)$。

离散傅里叶逆变换(IDFT)

离散傅里叶变换(IDFT) ，可以将点值快速转化为系数，从而得出系数多项式。

需要用到单位根反演：

$$\frac{1}{n} \sum_{i = 0} ^ {n - 1} \omega_n ^ {ix} = [x \bmod n = 0]$$

证明：

由于 $\omega _ n ^ {ix} = \omega_n ^ {(i-1)x} \times \omega_n ^ x$，

所以 $\omega_n ^ {ix}$ 为等比数列，

$$\frac{1}{n} \sum_{i = 0} ^ {n - 1} \omega_n ^ {ix} = \left\{ \begin{aligned} &\frac{1}{n} \sum_{i = 0} ^ {n - 1} 1 ^ i = 1 &x \bmod n = 0 \\ &\frac{1}{n} \times \frac{1 - \omega_n ^ {nx}}{1 - \omega_n ^ x} = 0 & x \bmod n \not = 0 \end{aligned} \right.$$

接下来就是IDFT的过程：将对应点值共轭带入求出后除以 $len$ 即可。

递归实现

由于需要用到复数需要自己封装。

void FFT(Complex *a, int len, int type)
{
    if(len == 1)return;
    Complex a1[(len >> 1) + 10], a2[(len >> 1) + 10];
    for(int i = 0; i <= len; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    FFT(a1, len >> 1, type);
    FFT(a2, len >> 1, type);
    Complex x = Complex(cos(2.0 * pi / len), type * sin(2.0 * pi / len));
    Complex w = Complex(1, 0);
    for(int i = 0; i < (len >> 1); i++, w = w * x)
    {
        a[i] = a1[i] + w * a2[i];
        a[i + (len >> 1)] = a1[i] - w * a2[i];
    }
}

最后需要除以长度，但是递归实现。

蝴蝶优化

观察分治的过程最后的结果其实就是对应二进制的反转，可以对此对分治进行优化。

预处理一个’rev’数组，

for(int i = 0; i <= len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));

然后FFT的过程就变成了

void FFT(Complex *a, int len, int type)
{
    for(int i = 0; i < len; i++)
        if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int k = 1; k < len; k <<= 1)
    {
        Complex x(cos(pi / k), type * sin(pi / k));
        for(int i = 0; i < len; i += k << 1)
        {
            Complex w(1, 0);
            for(int j = 0; j < k; j++)
            {
                Complex y = a[i + j];
                Complex z = w * a[i + j + k];
                a[i + j] = y + z;
                a[i + j + k] = y - z;
                w = w * x;
            }
        }
    }
    if(type == -1)
        for(int i = 0; i < len; i++)
            a[i].a /= len;
}

快速数论变换(NTT)

快速数论变换(NTT) 相比于FFT虽然时间复杂度均为 $O(n \log n)$，但是FFT会出现精度问题，不好解决取模，NTT可以很好的解决这个问题，但是常数会大。

原根

原根定义为：设 $m$ 为正整数，$a$ 是整数，若 $a \bmod m$ 的阶等于 $\varphi (m)$，则称 $a$ 为 $\bmod m$ 的一个原根。

原根有一个很重要的性质可以支持像FFT中单位根一样的运算，即：若 $P$ 为素数， 假设一个数 $g$ 是 $P$ 的原根， 那么 $g^i \bmod P$ 的结果两两不同。

可以得到：

$$\omega_n \equiv g^{\frac{p - 1}{n}} \pmod p$$

然后我们就可以将FFT中的 $\omega _n$ 替换为 $g^{\frac{p - 1}{n}}$

但是注意的是NTT对模数有要求，其模数必须要满足原根的定义，否则是不能使用NTT的，比如 $998244353$ 就为NTT模数， 其原根为 $3$。

void NTT(int *a, int len, int type)
{
    for(int i = 0; i < len; i++)
        if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int k = 1; k < len; k <<= 1)
    {
        int x = qpow(type == 1 ? g : gi, (mod - 1) / (k << 1));
        for(int i = 1; i < len; i += k << 1)
        {
            int w = 1;
            for(int j = 0; j < k; j++)
            {
                int y = a[i + j];
                int z = w * a[i + j + k] % mod;
                a[i + j] = (y + z) % mod;
                a[i + j + k] = (y - z + mod) % mod;
                w = w * x % mod;
            }
        }
    }
    if(type == -1)
    {
        int iv = qpow(len, mod - 2);
        for(int i = 0; i < len; i++)
            a[i] = a[i] * iv % mod;
    }
}

快速沃尔什变换(FWT)

给定两个长度为 $2 ^ n$ 的两个序列 $A,B$， 求序列 $C$

$$C_i = \sum_{j \oplus k = i} A_j \times B_k$$

其中 $\oplus$ 表示位运算与，或，异或。

或运算

首先求序列 $FWT[A] = \sum_{i = i | j} A _ j$，来求出满足条件的 $i$ 的子集，显然会有

$$C_i = \sum_{i = j | k} A_j \times B_k \Rightarrow FWT[C] = FWT[A] \times FWT[B]$$

接下来就是考虑如何进行 $FWT$ 运算， 有

$$FWT[A] = merge(FWT[A_0], FWT[A_0] + FWT[A_1]) \\ IFWT[A] = merge(IFWT[A_0], IFWT[A_1] - FWT[A_0])$$

void FWT_or(int *a, int len, int type)
{
    for(int k = 1; k < len; k <<= 1)
    {
        for(int i = 0; i < len; i += (k << 1))
        {
            for(int j = 0; j < k; j++)
            {
                a[i + j + k] = ((a[i + j + k] + a[i + j] * type + mod) % mod + mod) % mod;
            }
        }
    }
}

与运算

同或运算，有

$$FWT[A] = merge(FWT[A_0] + FWT[A_1], FWT[A_1])\\ IFWT[A] = merge(IFWT[A_0] - IFWT[A_1], IFWT[A_1])$$

void FWT_and(int *a, int len, int type)
{
    for(int k = 1; k < len; k <<= 1)
    {
        for(int i = 0; i < len; i += (k << 1))
        {
            for(int j = 0; j < k; j++)
            {
                a[i + j] = ((a[i + j] + a[i + j + k] * type + mod) % mod + mod) % mod;
            }
        }
    }
}

异或运算

推导得

$$FWT[A] = merge(FWT[A_0] + FWT[A_1], FWT[A_0] - FWT[A_1]) \\ IFWT[A] = merge(\frac{IFWT[A_0] + IFWT[A_1]}{2}, \frac{IFWT[A_0] - IFWT[A_1]}{2}) \\$$

void FWT_xor(int *a, int len, int type)
{
    for(int k = 1; k < len; k <<= 1)
    {
        for(int i = 0; i < len; i += (k << 1))
        {
            for(int j = 0; j < k; j++)
            {
                int x = a[i + j], y = a[i + j + k];
                a[i + j] = ((x + y) % mod * type + mod) % mod;
                a[i + j + k] = ((x - y + mod) % mod * type + mod) % mod;
            }
        }
    }
}

任意模数快速数论变换

普通的NTT对模数是有要求的其必须满足原根的相关定义，模数必须可以写成 $a \cdot 2 ^ k + 1$ 的形式。

$$469762049 = 7 \times 2 ^ {26} + 1 (g = 3) \\ 998244353 = 119 \times 2 ^{23} + 1(g = 3) \\ 1004535809 = 479 \times 2 ^ {21} + 1(g = 3)$$

如果题目中模数为 $1e9 + 7$，那么NTT就会受到限制，然后就可以使用任意模数NTT，（也可以称为三模数NTT），

计算时可以先找三个大质数， 分别计算结果，然后用中国剩余定理CRT合并即可。

首先记三次NTT的结果为：

$$ans \equiv a_1 \pmod {p_1} \\ ans \equiv a_2 \pmod {p_2} \\ ans \equiv a_3 \pmod {p_3}$$

先合并前两个得到:

$$ans \equiv a_4 \pmod {p_1 p_2}$$

将其转化为等式为：

$$ans = k p_1 p_2 + a_4$$

接着求 $k$：

$$k = (a_3 - a_4)p_1^{-1} p_2 ^ {-1} \pmod {p_3}$$

所以：

$$ans \equiv kp_1 p_2 + a_4 \pmod {p_1p_2p_3}$$

拆系数FFT

拆系数FFT也可以用来解决任意模数多项式乘法，由于普通的FFT精度误差大，需要进行优化，考虑将一个较大的数拆成两个较小的数，有 $n$ 次多项式 $A$ 和 $m$ 次多项式 $B$,

$$A = \sum_{i = 0} ^ n a_i x ^ i \\ B = \sum_{i = 0} ^ m b_i x ^ i \\$$

假设有一实数 $P$，将多项式的系数拆为

$$A = \sum_{i = 0} ^ n (x_1 \sqrt p + y_1) x ^ i \\ B = \sum_{i = 0} ^ m (x_2 \sqrt p + y_2) x ^ i \\$$

再次拆解成两个多项式，

$$A = \sum_{i = 0} ^ n x_1 \sqrt p x ^ i + \sum_{i = 0} ^ n y_1 x ^ i = \sqrt p C_1 + D_1 \\ B = \sum_{i = 0} ^ m x_2 \sqrt p x ^ i + \sum_{i = 0} ^ m y_2 x ^ i = \sqrt p C_2 + D_2 \\$$

进行多项式乘法可得,

$$A \ast B = (\sqrt p C_1 + D_1) \ast (\sqrt p C_2 + D_2) = pC_1C_2 + \sqrt p (C_1D_2 + C_2D_1) + D_1D_2$$

由于多次进行FFT常数极大，不过有些式子可以合并，降低FFT次数。

多项式乘法逆

定义多项式 $F ^ {-1}$ 为多项式 $F$ 的乘法逆元，满足

$$F \ast F ^ {-1} \equiv 1 \pmod {x^n}$$

假设我们已得知 $F \ast G' \equiv 1 \pmod {x ^ \frac{n}{2}}, F \ast G \equiv 1 \pmod {x ^ n}$

可以推出，

$$F \ast G \equiv 1 \pmod {x ^ n} \\ G - G' \equiv 0 \pmod {x ^ n} \\ (G - G') ^ 2 \equiv 0 \pmod {x^n}$$

进一步展开得，

$$G'^ 2 - 2GG' + G ^ 2 \equiv 0 \pmod {x ^ n}$$

然后两边同时卷上 $F$，

$$FG'^2 - 2G' + G \equiv 0 \pmod {x ^ n}$$

得出

$$G \equiv 2G' - FG'^ 2 \pmod {x ^ n}$$

可以递归实现，时间复杂度 $O(n \log n)$。

多项式对数函数

定义多项式对数函数为

$$G = \ln (F) \pmod {x ^ n}$$

假设有多项式 $F$ 和多项式 $G$， 记 $G = \ln F \pmod {x ^ n}$

$$G' \equiv (\ln F)' \pmod {x ^ n} \\ G' \equiv (\ln' F) F' \pmod {x ^ n} \\ G' \equiv \frac{F'}{F} \pmod {x ^ n}$$

需要注意的是这种方法需要保证 $F_0 = 1$。

多项式指数函数

定义多项式指数函数为

$$G(x) = e ^ {F(x)} \pmod {x ^ n}$$

牛顿迭代

牛顿迭代用于求函数零点，通过不断地切线逼近所求值，但最终也只是近似值，迭代的次数越多，精确度越高，误差越小。

假如我们要对一个非常大的数 $a$ 开方，手算，利用牛顿法来解决这个问题，其实本质上是求得 $f(x) = x ^2 - a$ 精确到整数得零点，假设我们已经求得了一个近似值 $x_0$，那么我们只需要过 $(x_0, f(x_0))$ 这个点， 作这个函数图像的切线，取切线与 $x$ 轴的交点作为新的 $x_0$。

假设我们要求一个函数 $f(x)$ 的零点， 初始近似值是 $x_0$，则切线方程为

$$y = f'(x_0)(x - x_0) + f(x_0)$$

令 $y = 0$，得到 $x = x_0 - \frac{f(x_0)}{f'(x_0)}$。

假设我们现在要求 $F(G(x)) \equiv 0$，然后利用上面的式子每一次令

$$G(x) = G_0(x) - \frac{F(G_0(x))}{F'(G_0(x))}$$

然后就可以很快的逼近真实值。

接下来推一下多项式exp

$$B(x) \equiv e ^ {A(x)} \pmod {x ^ n} \\ \ln B(x) - A(x) \equiv 0 \pmod {x^ n}$$

现在问题变为了使得 $F(G(x)) = \ln G(x) - A(x) \equiv 0$。

然后求导，

$$F'(G_0(x)) = \frac{1}{G_0(x)}$$

然后接着带入上面牛顿迭代的式子，

$$G(x) = {G_0(x)(1 - \ln G_0(x) + A(x))}$$

每次迭代，使用多项式求 $\ln$，然后再做一遍多项式乘法，然后就可以得到答案，时间复杂度 $O(n \log n)$。
需要保证$F_0 = 0$。

多项式开根

多项式开根用来解决

$$G^2(x) \equiv F(x) \pmod {x^n}$$

假设我们有 $G'^2(x) \equiv F(x) \pmod {x ^ {\frac{n}{2}}}, H(G(x)) = G^2(x) - F$，求 $G^2(x) \equiv F(x) \pmod {x ^ n}$，

$$G'^2 (x) \equiv F(x) \mod x ^ {\frac{n}{2}} , G^2(x) \equiv F(x) \pmod {x ^ {\frac{n}{2}}} \\ G^2(x) - F \equiv 0 \pmod {x ^ {\frac{n}{2}}} \\ H(G) \equiv 0 \pmod {x ^ {\frac{n}{2}}} \\ G \equiv G' - \frac{H(G')}{H'(G')} \pmod {x ^ n} \\ G \equiv \frac{G'^2 + F} {2G'} \pmod {x ^ n}$$

需要保证 $F_0 = 1$。

多项式幂函数

多项式幂函数是用来解决

$$G(x) \equiv (F(x)) ^ k \mod x ^ n$$

先求一遍 $\ln$ 然后乘以 $k$ 再使用 $\exp$，就好啦。

需保证$F_0 = 1$。

多项式的一些普通情况

多项式求ln

不保证 $F_0 = 1$。不存在，有定理：

在模意义下当且仅当 $F_0 = 1$， $F(x)$ 有对数多项式问题。

多项式求exp

不保证 $F_0 = 0$ 。同多项式求 $\ln$。

多项式开根

不保证 $F_0 = 1$，但保证 $F_0$ 是 $\bmod 998244353$ 下的二次剩余。

边界求一遍二次剩余即可。

多项式幂函数

不保证 $F_0 = 1$。可以先找到系数不为$0$的一项，然后让式子除以这一项最后再乘回来就好了

$$F(x)^k = \bigg(\frac{F(x)}{x ^ t} \bigg) ^ k x ^ {tk}$$

分治FFT/NTT

给定序列 $g$ 和 $f$，其中

$$f_i = \sum_{j = 1} ^ i f_{i - j}g_j$$

可以采用如同CDQ分治的方法统计前面的对后面的贡献，时间复杂度 $O(n \log ^ 2 n)$。

void Divide(int l, int r)
{
    if(l == r)return;
    int mid = (l + r) >> 1;
    Divide(l, mid);    
    int n = mid - l + 1;
    int m = r - l + 1;
    int len; init(n, m, len);
    for(int i = 0; i <= len; i++)
        A[i] = B[i] = 0;
    for(int i = l; i <= mid; i++)
        A[i - l] = F[i];
    for(int i = 1; i <= r - l; i++)
        B[i] = G[i];
    NTT(A, len, 1); NTT(B, len, 1);
    for(int i = 0; i < len; i++)
        A[i] = A[i] * B[i] % mod;
    NTT(A, len, - 1);
    for(int i = mid + 1; i <= r; i++)
        F[i] = (F[i] + A[i - l]) % mod; 
    Divide(mid + 1, r);
}

求 $f$， 这里还有一个多项式求逆的方法

设 $F(x) = \sum_{i = 0} ^ {\infty} f_i x ^ i , G(x) = \sum_{i = 0} ^ {\infty}g_i x ^ i$，且 $g_0 = 0$，

所以有

$$F(x) G(x) = \sum_{i = 0} ^ {\infty} \sum_{j + k = i}f_jg_k = F(x) - f_0 x ^ 0 \\ F(x)G(x) \equiv F(x) - f_0 \pmod {x ^ n} \\ F(x) \equiv (1 - G(x)) ^ {-1} \pmod {x ^ n}$$

多项式除法

给定一个 $n$ 次多项式 $F(x)$ 和一个 $m$ 次多项式 $G(x)$， 求多项式 $A(x),B(x)$，满足：

• $A(x)$ 次数为 $n - m$， $B(x)$ 次数小于 $m$
• $F(x) = A(x) \ast G(x) + B(x)$

所有运算在模 $998244353$ 下进行。
定义一种让多项式反转的操作为 $A'(x) = x^n A(\frac{1}{x})$，然后化简式子

$$F(x) = A(x) \ast G(x) + B(x) \\ x^n F(\frac{1}{x}) = x^{n - m} A(\frac{1}{x}) \ast x^m G(\frac{1}{x}) + x^{n - m + 1}\ast x ^{m - 1} B(\frac{1}{x}) \\ F'(x) = A'(x) \ast G'(x) + x^{n - m + 1} \ast B'(x) \\ F'(x) \equiv A'(x) \ast G'(x) \pmod {x^{n - m + 1}} \\ A'(x) \equiv F'(x) \ast G'^{-1}(x) \pmod {x^{n - m + 1}}$$

先进行多项式求逆，然后再推 $B(x) = F(x) - A(x) \ast G(x)$。

多项式多点求值

多项式朴素算法求值，单次是 $O(n)$ 的，多次 $O(n^2)$，对于多点求值时，可以用类似分治的算法优化到 $O(n \log ^ 2 n)$ 。

对于多项式 $F$，和一个点 $x_0$，我们可以构造多项式 $F(x) = G(x) \ast Q(x) + R(x)$， 若 $G(x) = x - x_0$， 就可以得到 $R(x_0)$ 即为对应的点值， 可以用分治的思想，进行多项式取模， 分别对 $\prod_{i = l} ^ mid (x - x_i)$ 和 $\prod _ {i = mid + 1} ^ r (x - x_i)$ 取模后递归计算。

多项式快速插值

对于 $n$ 个已知的点值，直接使用拉格朗日插值，时间复杂度是 $O(n ^ 2)$ 的，使用分治可以将其优化到 $O(n \log ^ 2 n)$ 。

首先对于拉格朗日插值公式，

$$f(x) = \sum_{i = 0} ^ n y_i \prod_{i \not = j} \frac{x - x_j}{x_i - x_j}$$

转化成，

$$f(x) = \sum_{i = 0} ^ n \frac{y_i}{\prod_{i \not = j} (x_i - x_j)} \prod_{i \not = j} (x - x_j)$$

设 $g(x) = \prod_{i = 1} ^ n (x - x_i)$ 然后相当于求 $\frac{g(x)}{x - x_i}$，根据洛必达法则，若

$$\lim_{x \to a} f(x) = 0, \lim_{x \to a} g(x) = 0$$

有

$$\lim_{x \to a} \frac{f(x)}{g(x)} = \lim_{x \to a} \frac{f'(x)}{g'(x)}$$

然后代入后可以得出，所求即为 $g'(x_i)$，这样就可以先算出 $g$ 然后求导进行一次多点求值，然后考虑分治，

设 $f_{l,r}$ 表示 $l$ 到 $r$ 之间的点插值得出的答案，则有

$$\begin{aligned} f_{l, r} &= \sum_{i = l} ^ r \frac{y_i}{g'(x_i)} \prod_{j = l, i \not = j} ^ r (x - x_j) \\ &= \prod_{j = mid + 1} ^ r (x - x_j) \sum_{i = l} ^ {mid} \frac{y_i}{g'(x_i)} \prod _ {j = l, i \not = j} ^ {mid} (x - x_j) + \prod_{j = l} ^ {mid} (x - x_j) \sum_{i = mid + 1} ^ r \frac{y_i}{g'(x_i)} \prod_{j = mid + 1, i \not = j} ^ r (x - x_j) \\ &= \prod_{j = mid + 1} ^ r (x - x_j) f_{l, mid} + \prod_{j = l} ^ {mid} (x - x_j) f_{mid + 1, r} \end{aligned}$$

直接分治计算即可。

多项式复合函数

有 $F(x), G(x)$，求

$$H(x) \equiv F(G(x)) \pmod {x^{n + 1}}$$

即：

$$H(x) \equiv \sum_{i = 0} ^n[x^i] F(x) \times G(x)^i \pmod {x^{n + 1}}$$

对 $998244353$ 取模。

设 $m = \sqrt n$，则

$$\sum_{i = 0} ^ n [x^i] F(x)G(x) ^ i = \sum_{i = 0} ^ {m - 1} \sum_{j = 0} ^ {m - 1}[x^{im + j}]F(x)G(x)^{im + j} = \sum_{i = 0}^{m - 1}G(x)^{im}\sum_{j = 0} ^ {m - 1}[x^{im+j}]F(x)G(x)^j$$

然后预处理 $G(x)^{im}$ 和 $G(x)^j$，其它的直接暴力计算， 时间复杂度 $O(n^2 + n \sqrt n \log n)$。

多项式复合函数逆

有 $F(x)$，求

$$G(F(x)) \equiv x \pmod {x ^ n}$$

对 $998244353$ 取模。

有拉格朗日反演公式：

$$[x^n]G(x) = \frac{1}{n}[x^{n - 1}](\frac{x}{F(x)})^n$$

这个是求单项系数的， 需要求的是所有系数。

利用解决多项式复合函数的方法， 令 $m = \sqrt n$， 则

$$G(x) = \sum_{i = 1}^{n}\bigg (\frac{1}{i}[x^{i - 1}] (\frac{x}{F(x)})^i \bigg) = \sum_{i = 0}^{m - 1}\sum_{j = 1} ^{m }\bigg(\frac{1}{im + j}[x^{im + j - 1}](\frac{x}{F(x)})^{im + j} \bigg)x^{im + j} = \sum_{i = 0}^{m - 1}\sum_{j = 1} ^ m \bigg( \frac{1}{im + j} [x^{im + j - 1}](\frac{x}{F(x)})^im (\frac{x}{F(x)})^j \bigg)x^{im + j}$$

然后就是和多项式复合函数的一样的处理方法，时间复杂度 $O(n^2 + \sqrt n \log n)$。

下降幂多项式乘法

给定两个下降幂多项式求 $F \ast G$。

先构造出点值的EGF，$F' = \sum_{i = 0} ^ {\infty} \frac{F(i)}{i!} x ^ i$,

$$F' = \sum_{i = 0} ^ {\infty} \frac{x ^ i}{i!} \sum_{j = 0} ^ n \frac{i!}{(i - j)!} \times F_j \\ F' = \sum_{i = 0} ^ {\infty} x ^ i \sum_{j = 0} ^ n \frac{1}{(i - j)!} \times F_j\\$$

换一下枚举顺序，

$$F' = \sum_{i = 0} ^ n F_i \sum_{j = i} ^ {\infty} \frac{1}{(j - i)!} x ^ j \\ F' = \sum_{i = 0} ^ n F_i x ^ i \sum_{j = 0} ^ {\infty} \frac{1}{j!} x ^ j \\$$

最后可以得出，

$$F' = e ^ x \sum_{i = 0} ^ n F_i x ^ i$$

也就是说直接卷上 $\sum _ {i = 0} ^ {\infty} \frac{1}{i!} x ^ i$ 即可得到下降幂多项式点值的EGF。

得到点值EGF之后，直接相乘，然后再卷上 $e ^ {-x} = \sum_{i = 0} ^ {\infty} \frac{(-1) ^ i}{i!} x ^ i$ 就可以得到卷积结果。

普通多项式转下降幂多项式

根据下降幂多项式和其点值EGF的关系，我们可以先对普通多项式进行多点求值，然后求出EGF，在卷上$e ^ {-x}$ 即可得出。 时间复杂度 $O(n \log ^ 2 n)$。

还可以用分治NTT，同样复杂度，但常数更小。

下降幂多项式转普通多项式

先卷上$e ^ x$然后就可求出点值的EGF， 然后由于点值是连续的，可以直接阶乘优化插值，时间复杂度$O(n log ^ 2 n)$

同样也有分治NTT的做法。

多项式三角函数

可以使用幂函数解决

$$\cos F(x) = \frac{e ^ {iF(x)} + e ^ {-iF(x)}}{2} \\ \sin F(x) = \frac{e ^ {iF(x)} - e ^ {-iF(x)}}{2i}$$

多项式反三角函数

见洛谷P5265

多项式相关

模板

预处理和IO

#define ll long long
#define poly vector <ll>
const int N = 4e5 + 10; // size
const int mod = 998244353;
const int g = 3;
const int gi = 332748118;
const int inv2 = (mod + 1) >> 1;

namespace MOD
{
    ll add(ll x, ll y){return x + y >= mod ? x + y - mod : x + y;}
    ll dec(ll x, ll y){return x - y < 0 ? x - y + mod : x - y;}
}
using namespace MOD;

// 快模

namespace IO
{
    inline ll read()
    {
        ll x = 0, f = 1;char ch = getchar();
        while(!isdigit(ch)){if(ch == '-')f = -1;ch = getchar();}
        while(isdigit(ch)){x =((x << 3ll) + (x << 1ll)) + (ch ^ 48ll);ch = getchar();}
        return x * f;
    }
    inline void write(ll x)
    {
        char ch[100];int len = 0;
        if(x == 0)ch[++len] = '0';
        if(x < 0)putchar('-'), x = -x;
        while(x != 0){ch[++len] = x % 10ll + '0'; x /= 10ll;}
        while(len)putchar(ch[len--]);
    }
}

using namespace IO;

// 快读

inline void read(poly &A, int n)
{
    A.resize(n + 1);
    for(int i = 0; i <= n; i++)
        A[i] = read();
}

inline void print(poly B)
{
    for(auto x : B)write(x), putchar(' ');
}

// 读入和输出多项式

inline ll qpow(ll a, ll b)
{
    ll t = 1;
    while(b != 0)
    {
        if(b & 1)t = t * a % mod;
        a = a * a % mod; b >>= 1;
    }
    return t;
}

inline ll inv(ll x)
{
    return qpow(x, mod - 2);
}

// 快速幂和逆元

ll fv[N];

inline void initfv()
{
    fv[1] = 1;
    for(ll i = 2; i <= (N - 10) / 4; i++)
        fv[i] = fv[mod % i] * (mod - mod / i) % mod;
}

//预处理逆元， 积分用

ll p[3][50];

inline void prework()
{
    for(int k = 1, i = 1; k < N; k <<= 1, i++)
    {
        p[0][i] = qpow(g, (mod - 1) / (k << 1));
        p[1][i] = qpow(gi, (mod - 1) / (k << 1));
        p[2][i] = inv(k << 1);
    }
}

// 预处理快速幂

多项式乘法

主要使用NTT

inline void init(int n, int m, int &len)
{
    len = 1; int cnt = 0;
    while(len <= (n + m))len <<= 1, cnt++;
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
}

// 预处理rev数组

inline void NTT(poly &a, int len, int type)
{
    for(int i = 0; i < len; i++)
        if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int k = 1, o = 1; k < len; k <<= 1, o++)
    {
        ll x = type == 1 ? p[0][o] : p[1][o];
        for(int i = 0; i < len; i += k << 1)
        {
            ll w = 1;
            for(int j = 0; j < k; j++)
            {
                ll y = a[i + j];
                ll z = w * a[i + j + k] % mod;
                a[i + j] = add(y, z);
                a[i + j + k] = dec(y, z);
                w = w * x % mod;
            }
        }
    }
    if(type == -1)
    {
        ll iv = p[2][(int)log2(len)];
        for(int i = 0; i < len; i++)
            a[i] = a[i] * iv % mod;
    }
}

//NTT

多项式基本运算

inline poly operator + (poly A, poly B)
{
    int n = max(A.size(), B.size());
    A.resize(n); B.resize(n);
    for(int i = 0; i < n; i++)
        A[i] = add(A[i], B[i]);
    return A;
}
inline poly operator - (poly A, poly B)
{
    int n = max(A.size(), B.size());
    A.resize(n); B.resize(n);
    for(int i = 0; i < n; i++)
        A[i] = dec(A[i], B[i]);
    return A;
}
inline poly operator * (poly A, poly B)
{
    int len, n = A.size() - 1, m = B.size() - 1;
    init(n, m, len);
    A.resize(len); B.resize(len);
    NTT(A, len, 1); NTT(B, len, 1);
    for(int i = 0; i < len; i++)
        A[i] = A[i] * B[i] % mod;
    NTT(A, len, -1); A.resize(n + m + 1);
    return A;
}

分治NTT

inline void Divide(int l, int r, poly &F, poly &G)
{
    if(l == r)return;
    int mid = (l + r) >> 1;
    Divide(l, mid, F, G);
    int n = mid - l + 1;
    int m = r - l + 1;
    poly A; A.resize(n);
    poly B; B.resize(m);
    for(int i = l; i <= mid; i++)
        A[i - l] = F[i];
    for(int i = 1; i <= r - l; i++)
        B[i] = G[i];
    A = A * B;
    for(int i = mid + 1; i <= r; i++)
        F[i] = add(F[i], A[i - l]);
    Divide(mid + 1, r, F, G);
}

多项式求逆

inline void Inverse(poly &A, poly &B, int n)
{
    if(n == 1)return void(B[0] = inv(A[0]));
    Inverse(A, B, (n + 1) >> 1);
    int len; init(n, n, len);
    poly C; C.resize(len); B.resize(len);
    for(int i = 0; i < n; i++)C[i] = A[i];
    NTT(C, len, 1); NTT(B, len, 1);
    for(int i = 0; i < len; i++)
        B[i] = dec(2, B[i] * C[i] % mod) * B[i] % mod;
    NTT(B, len, -1); B.resize(n);
}

多项式求导

inline void Diff(poly &A, poly &B)
{
    int n = A.size() - 1;
    B.resize(n);
    for(ll i = 1; i <= n; i++)
        B[i - 1] = i * A[i] % mod;
}

多项式积分

inline void Integ(poly &A, poly &B)
{
    int n = A.size();
    B.resize(n + 1);
    for(int i = 1; i <= n; i++)
        B[i] = A[i - 1] * fv[i] % mod;
    B[0] = 0;
}

多项式对数函数

inline void Ln(poly &A, poly &B)
{
    int n = A.size() - 1;
    poly F; poly G; G.resize(1);
    Diff(A, F); Inverse(A, G, n + 1);
    F = F * G;
    Integ(F, B); B.resize(n + 1);
}

多项式指数函数

inline void Exp(poly &A, poly &B, int n)
{
    if(n == 1)return void(B[0] = 1);
    Exp(A, B, (n + 1) >> 1); 
    poly G; B.resize(n);
    Ln(B, G); G[0] -= 1; B = (A - G) * B;
    B.resize(n);
}

多项式开根

inline void Sqrt(poly &A, poly &B, int n)
{
    if(n == 1)return void(B[0] = 1);
    Sqrt(A, B, (n + 1) >> 1);
    int len; init(n, n, len); B.resize(n);
    poly G, H; H.resize(1); Inverse(B, H, n);
    G.resize(len); H.resize(len); B.resize(len);
    for(int i = 0; i < n; i++)G[i] = A[i];
    NTT(H, len, 1); NTT(B, len, 1); NTT(G, len, 1);
    for(int i = 0; i < len; i++)
        B[i] = add(B[i], G[i] * H[i] % mod) * inv2 % mod;
    NTT(B, len, -1); B.resize(n);
}

多项式幂函数

inline void Pow(poly &A, poly &B, ll k)
{
    int n = A.size() - 1;
    poly C; Ln(A, C);
    for(int i = 0; i <= n; i++)
        C[i] = C[i] * k % mod;
    Exp(C, B, n + 1);
}

多项式带余除法

inline void Mod(poly &F, poly &G, poly &A, poly &B)
{
    poly f = F, g = G, h; 
    int n = F.size() - 1, m = G.size() - 1;
    reverse(f.begin(), f.end());
    reverse(g.begin(), g.end());
    g.resize(n - m + 1);
    h.resize(1); Inverse(g, h, n - m + 1);
    A = f * h; A.resize(n - m + 1);
    reverse(A.begin(), A.end());
    B = F - (A * G); B.resize(m);
}

多项式多点求值

inline void getV(int l, int r, int p)
{
    if(l == r)
    {
        v[p].resize(2);
        v[p][0] = mod - a[l];
        v[p][1] = 1;
        return;
    }
    int mid = (l + r) >> 1;
    getV(l, mid, p << 1);
    getV(mid + 1, r, p << 1 | 1);
    v[p] = v[p << 1] * v[p << 1 | 1];
}

inline void solve(int l, int r, int p, poly &f)
{
    if(l == r)return write(f[0]), putchar('\n'), void();
    if(r - l + 1 <= 1000)
    {
        for(int i = l; i <= r; i++)
        {
            ll x = 1, sum = 0;
            for(int j = 0; j < f.size(); j++)
            {
                sum = (sum + x * f[j] % mod) % mod;
                x = x * a[i] % mod;
            }
            write(sum); putchar('\n');
        }
        return;
    }
    int mid = (l + r) >> 1;poly L, R, A;
    Mod(f, v[p << 1], A, L);
    solve(l, mid, p << 1, L);
    Mod(f, v[p << 1 | 1], A, R);
    solve(mid + 1, r, p << 1 | 1, R);
}

多项式快速插值

inline void Larange(int l, int r, int p, poly &f)
{
    if(l == r)
    {
        f.resize(1);
        f[0] = Y[l] * inv(val[l]) % mod;
        return;
    }
    int mid = (l + r) >> 1; poly L, R;
    Larange(l, mid, p << 1, L);
    Larange(mid + 1, r, p << 1 | 1, R);
    f = (L * v[p << 1 | 1]) + (R * v[p << 1]); 
}

多项式模板全套

namespace Poly
{
    #define ll long long
    #define poly vector <ll>
    const int N = 4e5 + 10;
    const int mod = 998244353;
    const int g = 3;
    const int gi = 332748118;
    const int inv2 = (mod + 1) >> 1;

    namespace MOD
    {
        ll add(ll x, ll y){return x + y >= mod ? x + y - mod : x + y;}
        ll dec(ll x, ll y){return x - y < 0 ? x - y + mod : x - y;}
    }
    using namespace MOD;

    namespace IO
    {
        inline ll read()
        {
            ll x = 0, f = 1;char ch = getchar();
            while(!isdigit(ch)){if(ch == '-')f = -1;ch = getchar();}
            while(isdigit(ch)){x =((x << 3ll) + (x << 1ll)) + (ch ^ 48ll);ch = getchar();}
            return x * f;
        }
        inline void write(ll x)
        {
            char ch[100];int len = 0;
            if(x == 0)ch[++len] = '0';
            if(x < 0)putchar('-'), x = -x;
            while(x != 0){ch[++len] = x % 10ll + '0'; x /= 10ll;}
            while(len)putchar(ch[len--]);
        }
    }

    using namespace IO;

    int rev[N];

    ll p[3][50], fv[N];

    inline void read(poly &A, int n)
    {
        A.resize(n + 1);
        for(int i = 0; i <= n; i++)
            A[i] = read();
    }

    inline void print(poly B)
    {
        for(auto x : B)write(x), putchar(' ');
    }

    inline ll qpow(ll a, ll b)
    {
        ll t = 1;
        while(b != 0)
        {
            if(b & 1)t = t * a % mod;
            a = a * a % mod; b >>= 1;
        }
        return t;
    }

    inline ll inv(ll x)
    {
        return qpow(x, mod - 2);
    }

    inline void initfv()
    {
        fv[1] = 1;
        for(ll i = 2; i <=(N - 10) / 4; i++)
            fv[i] = fv[mod % i] * (mod - mod / i) % mod;
    }

    inline void prework()
    {
        for(int k = 1, i = 1; k < N; k <<= 1, i++)
        {
            p[0][i] = qpow(g, (mod - 1) / (k << 1));
            p[1][i] = qpow(gi, (mod - 1) / (k << 1));
            p[2][i] = inv(k << 1);
        }
    }

    inline void init(int n, int m, int &len)
    {
        len = 1; int cnt = 0;
        while(len <= (n + m))len <<= 1, cnt++;
        for(int i = 0; i < len; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
    }

    inline void NTT(poly &a, int len, int type)
    {
        for(int i = 0; i < len; i++)
            if(i < rev[i])swap(a[i], a[rev[i]]);
        for(int k = 1, o = 1; k < len; k <<= 1, o++)
        {
            ll x = type == 1 ? p[0][o] : p[1][o];
            for(int i = 0; i < len; i += k << 1)
            {
                ll w = 1;
                for(int j = 0; j < k; j++)
                {
                    ll y = a[i + j];
                    ll z = w * a[i + j + k] % mod;
                    a[i + j] = add(y, z);
                    a[i + j + k] = dec(y, z);
                    w = w * x % mod;
                }
            }
        }
        if(type == -1)
        {
            ll iv = p[2][(int)log2(len)];
            for(int i = 0; i < len; i++)
                a[i] = a[i] * iv % mod;
        }
    }

    inline poly operator + (poly A, poly B)
    {
        int n = max(A.size(), B.size());
        A.resize(n); B.resize(n);
        for(int i = 0; i < n; i++)
            A[i] = add(A[i], B[i]);
        return A;
    }
    inline poly operator - (poly A, poly B)
    {
        int n = max(A.size(), B.size());
        A.resize(n); B.resize(n);
        for(int i = 0; i < n; i++)
            A[i] = dec(A[i], B[i]);
        return A;
    }
    inline poly operator * (poly A, poly B)
    {
        int len, n = A.size() - 1, m = B.size() - 1;
        init(n, m, len);
        A.resize(len); B.resize(len);
        NTT(A, len, 1); NTT(B, len, 1);
        for(int i = 0; i < len; i++)
            A[i] = A[i] * B[i] % mod;
        NTT(A, len, -1); A.resize(n + m + 1);
        return A;
    }

    inline void Divide(int l, int r, poly &F, poly &G)
    {
        if(l == r)return;
        int mid = (l + r) >> 1;
        Divide(l, mid, F, G);
        int n = mid - l + 1;
        int m = r - l + 1;
        poly A; A.resize(n);
        poly B; B.resize(m);
        for(int i = l; i <= mid; i++)
            A[i - l] = F[i];
        for(int i = 1; i <= r - l; i++)
            B[i] = G[i];
        A = A * B;
        for(int i = mid + 1; i <= r; i++)
            F[i] = add(F[i], A[i - l]);
        Divide(mid + 1, r, F, G);
    }

    inline void Inverse(poly &A, poly &B, int n)
    {
        if(n == 1)return void(B[0] = inv(A[0]));
        Inverse(A, B, (n + 1) >> 1);
        int len; init(n, n, len);
        poly C; C.resize(len); B.resize(len);
        for(int i = 0; i < n; i++)C[i] = A[i];
        NTT(C, len, 1); NTT(B, len, 1);
        for(int i = 0; i < len; i++)
            B[i] = dec(2, B[i] * C[i] % mod) * B[i] % mod;
        NTT(B, len, -1); B.resize(n);
    }

    inline void Diff(poly &A, poly &B)
    {
        int n = A.size() - 1;
        B.resize(n);
        for(ll i = 1; i <= n; i++)
            B[i - 1] = i * A[i] % mod;
    }

    inline void Integ(poly &A, poly &B)
    {
        int n = A.size();
        B.resize(n + 1);
        for(int i = 1; i <= n; i++)
            B[i] = A[i - 1] * fv[i] % mod;
        B[0] = 0;
    }

    inline void Ln(poly &A, poly &B)
    {
        int n = A.size() - 1;
        poly F; poly G; G.resize(1);
        Diff(A, F); Inverse(A, G, n + 1);
        F = F * G;
        Integ(F, B); B.resize(n + 1);
    }

    inline void Exp(poly &A, poly &B, int n)
    {
        if(n == 1)return void(B[0] = 1);
        Exp(A, B, (n + 1) >> 1); 
        poly G; B.resize(n);
        Ln(B, G); G[0] -= 1; B = (A - G) * B;
        B.resize(n);
    }

    inline void Sqrt(poly &A, poly &B, int n)
    {
        if(n == 1)return void(B[0] = 1);
        Sqrt(A, B, (n + 1) >> 1);
        int len; init(n, n, len); B.resize(n);
        poly G, H; H.resize(1); Inverse(B, H, n);
        G.resize(len); H.resize(len); B.resize(len);
        for(int i = 0; i < n; i++)G[i] = A[i];
        NTT(H, len, 1); NTT(B, len, 1); NTT(G, len, 1);
        for(int i = 0; i < len; i++)
            B[i] = add(B[i], G[i] * H[i] % mod) * inv2 % mod;
        NTT(B, len, -1); B.resize(n);
    }

    inline void Pow(poly &A, poly &B, ll k)
    {
        int n = A.size() - 1;
        poly C; Ln(A, C);
        for(int i = 0; i <= n; i++)
            C[i] = C[i] * k % mod;
        Exp(C, B, n + 1);
    }

    inline void Mod(poly &F, poly &G, poly &A, poly &B)
    {
        poly f = F, g = G, h; 
        int n = F.size() - 1, m = G.size() - 1;
        if(n < m)return void(B = F);
        reverse(f.begin(), f.end());
        reverse(g.begin(), g.end());
        g.resize(n - m + 1);
        h.resize(1); Inverse(g, h, n - m + 1);
        A = f * h; A.resize(n - m + 1);
        reverse(A.begin(), A.end());
        B = F - (A * G); B.resize(m);
    }

    inline poly operator % (poly A, poly B)
    {
        poly F, G;
        Mod(A, B, F, G);
        return G;
    }
}

using namespace Poly;

多项式模板纯享版

namespace Poly
{
    #define ll long long  
    #define poly vector <ll>
    const int N = 4e5 + 10;const int mod = 998244353;const int g = 3;const int gi = 332748118;const int inv2 = (mod + 1) >> 1;
    namespace MOD{ll add(ll x, ll y){return x + y >= mod ? x + y - mod : x + y;}ll dec(ll x, ll y){return x - y < 0 ? x - y + mod : x - y;}}using namespace MOD;
    namespace IO{inline ll read(){ll x = 0, f = 1;char ch = getchar();while(!isdigit(ch)){if(ch == '-')f = -1;ch = getchar();}while(isdigit(ch)){x =((x << 3ll) + (x << 1ll)) + (ch ^ 48ll);ch = getchar();}return x * f;}inline void write(ll x){char ch[100];int len = 0;if(x == 0)ch[++len] = '0';if(x < 0)putchar('-'), x = -x;while(x != 0){ch[++len] = x % 10ll + '0'; x /= 10ll;}while(len)putchar(ch[len--]);}}using namespace IO;int rev[N];ll p[3][50], fv[N];
    inline void read(poly &A, int n){A.resize(n + 1);for(int i = 0; i <= n; i++)A[i] = read();}
    inline void print(poly B){for(auto x : B)write(x), putchar(' ');}
    inline ll qpow(ll a, ll b){ll t = 1;while(b != 0){if(b & 1)t = t * a % mod;a = a * a % mod; b >>= 1;}return t;}
    inline ll inv(ll x){return qpow(x, mod - 2);}
    inline void initfv(){fv[1] = 1;for(ll i = 2; i <=(N - 10) / 4; i++)fv[i] = fv[mod % i] * (mod - mod / i) % mod; }
    inline void prework(){for(int k = 1, i = 1; k < N; k <<= 1, i++){p[0][i] = qpow(g, (mod - 1) / (k << 1));p[1][i] = qpow(gi, (mod - 1) / (k << 1));p[2][i] = inv(k << 1);}}
    inline void init(int n, int m, int &len){len = 1; int cnt = 0;while(len <= (n + m))len <<= 1, cnt++;for(int i = 0; i < len; i++)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));}
    inline void NTT(poly &a, int len, int type){for(int i = 0; i < len; i++)if(i < rev[i])swap(a[i], a[rev[i]]);for(int k = 1, o = 1; k < len; k <<= 1, o++){ll x = type == 1 ? p[0][o] : p[1][o];for(int i = 0; i < len; i += k << 1){ll w = 1;for(int j = 0; j < k; j++){ll y = a[i + j];ll z = w * a[i + j + k] % mod;a[i + j] = add(y, z);a[i + j + k] = dec(y, z);w = w * x % mod;}}}if(type == -1){ll iv = p[2][(int)log2(len)];for(int i = 0; i < len; i++)a[i] = a[i] * iv % mod;}}
    inline poly operator + (poly A, poly B){int n = max(A.size(), B.size());A.resize(n); B.resize(n);for(int i = 0; i < n; i++)A[i] = add(A[i], B[i]);return A;}
    inline poly operator - (poly A, poly B){int n = max(A.size(), B.size());A.resize(n); B.resize(n);for(int i = 0; i < n; i++)A[i] = dec(A[i], B[i]);return A;}
    inline poly operator * (poly A, poly B){int len, n = A.size() - 1, m = B.size() - 1;init(n, m, len);A.resize(len); B.resize(len);NTT(A, len, 1); NTT(B, len, 1);for(int i = 0; i < len; i++)A[i] = A[i] * B[i] % mod;NTT(A, len, -1); A.resize(n + m + 1);return A;}
    inline void Divide(int l, int r, poly &F, poly &G){if(l == r)return;int mid = (l + r) >> 1;Divide(l, mid, F, G);int n = mid - l + 1;int m = r - l + 1;poly A; A.resize(n);poly B; B.resize(m);for(int i = l; i <= mid; i++)A[i - l] = F[i];for(int i = 1; i <= r - l; i++)B[i] = G[i];A = A * B;for(int i = mid + 1; i <= r; i++)F[i] = add(F[i], A[i - l]);Divide(mid + 1, r, F, G);}
    inline void Inverse(poly &A, poly &B, int n){if(n == 1)return void(B[0] = inv(A[0]));Inverse(A, B, (n + 1) >> 1);int len; init(n, n, len);poly C; C.resize(len); B.resize(len);for(int i = 0; i < n; i++)C[i] = A[i];NTT(C, len, 1); NTT(B, len, 1);for(int i = 0; i < len; i++)B[i] = dec(2, B[i] * C[i] % mod) * B[i] % mod;NTT(B, len, -1); B.resize(n);}
    inline void Diff(poly &A, poly &B){int n = A.size() - 1;B.resize(n);for(ll i = 1; i <= n; i++)B[i - 1] = i * A[i] % mod;}
    inline void Integ(poly &A, poly &B){int n = A.size();B.resize(n + 1);for(int i = 1; i <= n; i++)B[i] = A[i - 1] * fv[i] % mod;B[0] = 0;}
    inline void Ln(poly &A, poly &B){int n = A.size() - 1;poly F; poly G; G.resize(1);Diff(A, F); Inverse(A, G, n + 1);F = F * G;Integ(F, B); B.resize(n + 1);}
    inline void Exp(poly &A, poly &B, int n){if(n == 1)return void(B[0] = 1);Exp(A, B, (n + 1) >> 1); poly G; B.resize(n);Ln(B, G); G[0] -= 1; B = (A - G) * B;B.resize(n);}
    inline void Sqrt(poly &A, poly &B, int n){if(n == 1)return void(B[0] = 1);Sqrt(A, B, (n + 1) >> 1);int len; init(n, n, len); B.resize(n);poly G, H; H.resize(1); Inverse(B, H, n);G.resize(len); H.resize(len); B.resize(len);for(int i = 0; i < n; i++)G[i] = A[i];NTT(H, len, 1); NTT(B, len, 1); NTT(G, len, 1);for(int i = 0; i < len; i++)B[i] = add(B[i], G[i] * H[i] % mod) * inv2 % mod;NTT(B, len, -1); B.resize(n);}
    inline void Pow(poly &A, poly &B, ll k){int n = A.size() - 1;poly C; Ln(A, C);for(int i = 0; i <= n; i++)C[i] = C[i] * k % mod;Exp(C, B, n + 1);}
    inline void Mod(poly &F, poly &G, poly &A, poly &B){poly f = F, g = G, h;int n = F.size() - 1, m = G.size() - 1;if(n < m)return void(B = F); reverse(f.begin(), f.end());reverse(g.begin(), g.end());g.resize(n - m + 1);h.resize(1); Inverse(g, h, n - m + 1);A = f * h; A.resize(n - m + 1);reverse(A.begin(), A.end());B = F - (A * G); B.resize(m);}
    inline poly operator % (poly A, poly B){poly F, G;Mod(A, B, F, G);return G;}
}
using namespace Poly;

生成函数

普通生成函数(OGF)

定义序列 $a$ 的生成函数为

$$F(x) = \sum_{i = 0} ^ {\infty} a_ix^i$$

运算

其加法运算为

$$F(x) \pm G(x) = \sum_{i = 0} ^ {\infty} (a_i \pm b_i) x^ i$$

乘法运算，也就是卷积为

$$F(x) \ast G(x) = \sum_{i = 0} ^ {\infty} \bigg( \sum_{j = 0} ^ i a_jb_{i - j} \bigg)x^i$$

其为序列 $\sum_ {i = 0}^n a_i b_{n - i}$ 的生成函数。

封闭形式

根据等比数列求和公式可得

$$F(x) = \sum_{i = 0} ^ {\infty} a_i x^i = a_0 \times \frac{1 - x ^ {\infty}}{1 - x} = \frac{a_0}{1 - x}$$

牛顿二项式定理

定义组合数为

$$\dbinom{n}{m} = \frac{n ^ {\underline{m}}}{m!}, n \in \mathbb{C}, m \in \mathbb{N}$$

指数型生成函数(EGF)

定义序列 $a$ 的指数型生成函数为

$$F(x) = \sum_{i = 0} ^ {\infty} a_i \frac{x^i}{i !}$$

对于两个序列 $a$ 和 序列 $b$， 设它们的指数生成函数分别为 $F(x)$ 和 $G(x)$， 有

$$\begin{aligned} F \ast G &= \sum_{i = 0} ^ {\infty} a_i \frac{x ^ i}{i!} \sum_{j = 0} ^ {\infty} b_j \frac{x ^ j}{j!} \\ &= \sum_{i = 0} ^ {\infty} x ^ i \sum_{j = 0} ^ i a_jb_{i - j} \frac{1}{j!(i - j)!} \\ &= \sum_{i = 0} ^ {\infty} \frac{x ^ i}{i!} \sum_{j = 0} ^ i \dbinom{i}{j} a_j b_{i - j} \end{aligned}$$

两个函数的卷积对应的是序列 $\sum_{i = 0} ^ n \dbinom{n}{i} a_i b_{n - i}$ 的指数型生成函数。

显然有

$$e^x = \sum_{i = 1}^{\infty} \frac{e^0 x^i}{i!} = \sum_{i = 1}^{\infty}\frac{x^i}{i!} \\ \frac{e^x + e^{-x}}{2} = \sum_{i = 0}^{\infty} \frac{x^{2i}}{(2i)!} \\ \frac{e^x - e ^ {-x}}{2} =\sum_{i = 0}^{\infty} \frac{x ^ {2i + 1}}{(2i+1)!}$$

可以应用到组合数学中。