多项式合集

多项式合集

拉格朗日插值

问题背景

给出 $n$ 个点 $(x_i,y_i)$，令这 $n$ 个点确定的多项式为 $L(x)$，求 $L(k)\mathrm{m}\mathrm{o}\mathrm{d}998244353$ 的值。

结论

$$L(x)=\sum _{i=1}^n{y}_i{l}_i(x)$$

其中每个 $l_i(x)$ 为拉格朗日基本多项式，表达式为

$$l_i(x)=\prod _{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}$$

其特点是 $l_i(x_i)=1$，$\forall j\ne i$ 有 $l_i(x_j)=0$

推导

抛开拉插，这道题明显可以列方程组然后使用高斯消元求解，但是复杂度为 $O(n^3)$ 且精度问题明显，所以拉格朗日是这样考虑的：

对于每个点 $P_i(x_i,y_i)$，构造一个 $n-1$ 次多项式 $l_i(x)$ 使其在 $x_i$ 上取值为 $1$，在其余 $x_j$ 上为 $0$。构造的结果就是上面的结论：

$$l_i(x)=\prod _{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}$$

这个多项式的正确性还是很显然的。然后我们也知道这个多项式它就是唯一的。

然后考虑构造答案：很显然对于点 $P_i(x_i,y_i)$，只有 $l_i(x_i)$ 的取值为 $1$，其他的都为 $0$。所以答案的正确性也是比较显然的：对于 $x_i$，只有 $y_i{l}_i(x_i)$ 产生了贡献，其余的都是 $0$。故这个多项式是正确的。

所以回到一开始，我们需要的就是

$$f(k)=\sum _{i=1}^n{y}_i\prod _{j=1,j\ne i}^n\frac{k-x_j}{x_i-x_j}$$

由于模数是质数，所以使用费马小定理求逆元，跑得飞快。

复杂度 $O(n^2)$，求逆元就是个很小的常数

#include <cstdio>
#include <cctype>
#define il inline

typedef long long ll;

inline ll read()
{
    char c = getchar();
    ll s = 0;
    bool x = 0;
    while (!isdigit(c))
        x = x | (c == '-'), c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return x ? -s : s;
}

const ll maxn = 2e3 + 5, mod = 998244353;
ll x[maxn], y[maxn];

ll pow(ll base, ll p)
{
    ll ans = 1;
    base = (base + mod) % mod;
    for (; p; p >>= 1)
    {
        if (p & 1)
            ans = ans * base % mod;
        base = base * base % mod;
    }
    return ans;
}

il ll inv(ll n)
{
    return pow(n, mod - 2);
}

int main()
{
    ll n = read(), k = read();
    for (int i = 1 ; i <= n; ++i)
        x[i] = read(), y[i] = read();
    ll ans = 0;
    for (int i = 1; i <= n; ++i)
    {
        ll prod1 = 1, prod2 = 1;
        for (int j = 1; j <= n; ++j)
        {
            if (i == j)
                continue;
            prod1 = prod1 * (k - x[j]) % mod;
            prod2 = prod2 * (x[i] - x[j]) % mod;
        }
        ans = (ans + prod1 * y[i] % mod * inv(prod2) % mod + mod) % mod;
    }
    printf("%lld\n", ans);
    return 0;
}

拉格朗日插值与范德蒙矩阵

可以考虑将这 $n+1$ 个点值表示为如下形式：

$$\left[\begin{array}{ccccc}{x}_0^0& x_0^1& x_0^2& \cdots & x_0^n\\ x_1^0& x_1^1& x_1^2& \cdots & x_1^n\\ ⋮& ⋮& ⋮& & ⋮\\ x_n^0& x_n^1& x_n^2& \cdots & x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]$$

左边这个矩阵就是所谓的范德蒙德矩阵，记作 $\mathbf{V}$，系数列向量记作 $\mathbf{A}$，右边的记作 $\mathbf{B}$，则很明显：

$$\mathbf{V}\mathbf{A}=\mathbf{B}$$

打开来看清楚些实际就是多项式 $f$ 在每个点处的值：

$$y_j=f(x_j)=\sum _{i=0}^n{a}_i{x}_j^i$$

我们把两边都乘上 ${\mathbf{V}}^{-1}$：

$$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]={\mathbf{V}}^{-1}\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]$$

就得到了 $a_i$ 一定可以表示为某种形如

$$a_k=\sum \left[\begin{array}{c}⋮\end{array}\right]y_k$$

的形式，即 $a_k$ 只与 $x_i$ 与 $y_k$ 有关。

所以不难发现对于一个要求的 $f(x_{\ominus })$，都可以被表示为如下形式

$$f(x_{\ominus })=\sum {\delta }_k(x_{\ominus })y_k$$

${\delta }_k(x)$ 构造的过程即需要考虑 $x=x_k$ 时 ${\delta }_j(x)=0\wedge {\delta }_k(x)=1$，其中 $k\ne j$。${\delta }_j(x_k)=0$ 说明每一个 ${\delta }_j$ 都要有 $(x-x_k)$ 这个因式，然后又因为 ${\delta }_k(x_k)=1$，所以最终构造出来就是上面的结果：

$$f(x)=\sum _{i=1}^n{y}_i\prod _{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}$$

我们其实也可以利用拉格朗日插值来求范德蒙矩阵的逆阵，复杂度 $O(n^2)$

开始全家桶之前

形式化定义

约定：$f_i$ 表示 $f(x)$ 在 $x^i$ 处的系数，即一个多项式可以表示为 $\sum _{i=0}{f}_i{x}^i$ 的形式。

两个多项式的加减法定义为

$$f(x)\pm g(x)=\sum _{i=0}(f_i\pm g_i)x^i$$

复杂度 $O(n)$

两个多项式的乘法（加法卷积）定义为：

$$f(x)\ast g(x)=\sum _{i=0}{x}^i\sum _{j=0}{f}_j{g}_{i-j}$$

不难发现其正确性。可以手动模拟一下多项式的乘法看看是不是这样子的。其本质也就是卷完之后合并同类项。朴素的做的话复杂度是 $O(n^2)$ 的，下面要讲的 FFT/NTT 可以加速到 $O(n\log n)$

界

有些时候，题目只对多项式的前若干项感兴趣，所以我们给运算设定一个上界，即 $(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。意思就是只考虑这个多项式的前 $n$ 项，从 $x^n$ 开始以后的全部舍掉。

不难发现由加法和乘法是从低位到高位贡献的，所以

$$\begin{array}{rl}(f(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n\pm g(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n& =(f(x)\pm g(x))\mathrm{m}\mathrm{o}\mathrm{d}{x}^n\\ (f(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\ast (g(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n& =(f(x)\ast g(x))\mathrm{m}\mathrm{o}\mathrm{d}{x}^n\end{array}$$

下面我们就开始学习多项式的各种操作吧

快速傅里叶变换（FFT)

FFT 可以加速卷积，让时间复杂度从 $O(n^2)$ 降到 $O(n\log n)$，学习 FFT 的基础操作前，需要先了解复数，因为 FFT 就是基于单位复数根的良好性质实现的。

复数基础

（数学选修 2-2 内容）

定义虚数单位 ${\mathrm{i}}^2=\sqrt{-1}$，把形如 $a+b\mathrm{i}(a,b\in \mathbb{R})$ 的数称为复数，所有复数的集合称为复数集 $\mathbb{C}$。

复数一般使用 $z$ 表示，表示为 $z=a+b\mathrm{i}$，这种形式称为复数的代数形式。$a$ 被称为复数的实部，$b$ 称为复数的虚部，未加说明的情况下一般认为 $a,b\in \mathbb{R}$。很明显地，当 $a=0\wedge b\ne 0$ 时，这个复数为纯虚数，当 $b=0$ 时，这个复数为实数。

每个复数 $a+b\mathrm{i}$ 都能对应平面直角坐标系里面的一个点 $(a,b)$，同样的也可以对应一个向量 $(a,b)$。故定义复数的模为 $\sqrt{a^2+b^2}$。

定义复数的加法与乘法：
$$\begin{array}{rl}& (a+b\mathrm{i})+(c+d\mathrm{i})\\ =& (a+c)+(b+d)\mathrm{i}\end{array}$$

$$\begin{array}{rl}& (a+b\mathrm{i})(c+d\mathrm{i})\\ =& ac+ad\mathrm{i}+cb\mathrm{i}+bd{\mathrm{i}}^2\\ =& (ac-bd)+(ad+bc)\mathrm{i}\end{array}$$

这都是比较显然的。

容易看出复数满足很多实数的运算律。

定义复数 $z=a+b\mathrm{i}$ 的共轭复数为 $\overline{z}=a-b\mathrm{i}$，不难发现 $z$ 与 $\overline{z}$ 关于实轴对称。
$$z\overline{z}=(a+b\mathrm{i})(a-b\mathrm{i})=a^2+b^2=|z{|}^2$$
复数既然可以对应平面直角坐标系中的向量，不难发现其可以使用其模长与辐角来表示：
$$z=a+b\mathrm{i}⟺z=r(\cos \theta +\mathrm{i}\sin \theta )$$
其中 $r$ 为 $z$ 的模长，$\theta$ 为其辐角。即我们可以把一个复数表示成二元组 $(r,\theta )$ 的形式。

现在考虑两个复数 $(r_1,{\theta }_1)$ 与 $(r_2,{\theta }_2)$ 相乘得到的结果：
$$\begin{array}{rl}(r_1,{\theta }_1)\times (r_2,{\theta }_2)& =r_1(\cos {\theta }_1+\mathrm{i}\sin {\theta }_1)r_2(\cos {\theta }_2+\mathrm{i}\sin {\theta }_2)\\ & =r_1{r}_2(\cos {\theta }_1\cos {\theta }_2-\sin {\theta }_1\sin {\theta }_2+\mathrm{i}\sin {\theta }_1\cos {\theta }_2+\mathrm{i}\sin {\theta }_2\cos {\theta }_1)\\ & =r_1{r}_2\left(\cos ({\theta }_1+{\theta }_2)+\mathrm{i}\sin ({\theta }_1+{\theta }_2)\right)\\ & =(r_1{r}_2,{\theta }_1+{\theta }_2)\end{array}$$
于是我们可以概括复数乘法的法则：模长相乘，辐角相加。（上述推导需要掌握基本的三角恒等变换）

从欧拉公式到单位圆

给出复数指数幂的定义：
$${\mathrm{e}}^{x+y\mathrm{i}}=e^x(\cos y+\mathrm{i}\sin y)$$
这个公式是由我也不会证明的泰勒展开推导出来的：
$$\begin{array}{rl}\sin (x)& =x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\frac{x^9}{9!}+\cdots =\sum _{k=1}^{\infty }\frac{(-1{)}^{k-1}{x}^{2k-1}}{(2k-1)!}\\ \cos (x)& =1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\frac{x^8}{8!}+\cdots =\sum _{k=0}^{\infty }\frac{(-1{)}^k{x}^{2k}}{(2k)!}\\ {\mathrm{e}}^x& =1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots =\sum _{k=0}^{\infty }\frac{x^k}{k!}\end{array}$$
将 $x+y\mathrm{i}$ 代入进去即可推导。

如果 $x=0$，我们就得到大名鼎鼎的欧拉公式：
$${\mathrm{e}}^{x\mathrm{i}}=\cos x+\mathrm{i}\sin x$$
更特殊地，如果 $x=\pi$，得到的就是下面这个神奇的式子：
$${\mathrm{e}}^{\pi \mathrm{i}}=-1$$
复平面上我们可以定义类似于平面直角坐标系上的单位圆，单位圆上的所有复数构成集合 { z : ∣ z ∣ = 1 } \{z: |z| = 1\} {z:∣z∣=1}。这些复数都可以表示为 $\cos \theta +\mathrm{i}\sin \theta$ 或 $e^{\theta \mathrm{i}}$ 的形式。

多项式的表示法

系数表示法：顾名思义
f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n    ⟺    f ( x ) = { a 0 , a 1 , a 2 , ⋯   , a n } = [ x 0 x 1 x 2 ⋯ x n ] [ a 0 a 1 a 2 ⋮ a n ] f(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n\iff f(x) = \{a_0,a_1,a_2,\cdots,a_n\} = \begin{bmatrix} x^0 & x^1 & x^2 &\cdots & x^n \end{bmatrix} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_n \end{bmatrix} f(x)=a0​+a1​x+a2​x2+⋯+an​xn⟺f(x)={a0​,a1​,a2​,⋯,an​}=[x0​x1​x2​⋯​xn​]⎣⎢⎢⎢⎢⎢⎡​a0​a1​a2​⋮an​​⎦⎥⎥⎥⎥⎥⎤​
点值表示法：

我们知道由一个多项式在 $n+1$ 个点上的取值是可以唯一确定一个多项式的，其本质也就是线性方程组的解。所以一个 $n$ 次多项式可以用 $n+1$ 个点表示：

f ( x ) = { ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n , y n ) } f(x) = \{(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\} f(x)={(x0​,y0​),(x1​,y1​),⋯,(xn​,yn​)}

或者：

$$\left[\begin{array}{ccccc}{x}_0^0& x_0^1& x_0^2& \cdots & x_0^n\\ x_1^0& x_1^1& x_1^2& \cdots & x_1^n\\ ⋮& ⋮& ⋮& & ⋮\\ x_n^0& x_n^1& x_n^2& \cdots & x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]$$

通过下面的这个形式我们看得出来其就是一个典型的线性方程组的形式，不难证明其解的唯一性。

并且我们发现点值表示法有一个很明显的优势：可以在 $O(n)$ 的时间内将两个多项式乘起来，只需把对应点的 $y$ 乘起来即可。

通俗的来说，FFT 实现的就是快速求多项式乘法的过程：先把系数表示法转成点值表示法（DFT，离散傅里叶变换），乘完之后再把点值还原为插值（IDFT，离散傅里叶逆变换）。可是朴素的 DFT 需要的时间复杂度为 $O(n^2)$，IDFT 还回其系数需要高斯消元是 $O(n^3)$ 的。而 FFT 利用了一些很特殊很特殊的值加速了 DFT 和 IDFT 的过程，使得总时间复杂度降低到了 $O(n\log n)$。

单位复数根

解这个方程：
$$x^n=1$$
我们会发现这个方程在实数范围内只有 $1$ 或者 $2$ 个解。然而代数基本定理告诉我们这样的方程有 $n$ 个复数域上的解。由模长相乘辐角相加我们知道因为最终 $x^n=1$，所以这些满足条件的 $x$ 的模长必定也是 $1$。然后需要满足他们的辐角的 $n$ 倍能被 $2\pi$ 整除。

不难发现其就是 $n$ 等分单位圆：

我们记 $n$ 次单位根的第 $k$ 个记为 ${\omega }_n^k$，不难发现 ${\omega }_k^n={\mathrm{e}}^{\frac{2k\pi i}{n}}$。由此可见，单位复数根具有一些非常好的性质比如：
$$\begin{array}{rl}{\omega }_n^0={\omega }_n^n& =1\\ {\omega }_n^k& ={\omega }_{2n}^{2k}\\ {\omega }_{2n}^{k+n}& =-{\omega }_{2n}^k\\ {\left({\omega }_{2n}^{k+n}\right)}^2& ={\omega }_n^k\end{array}$$
利用这些性质，我们可以加速 DFT 的过程。FFT 就是利用分治思想加速求每个 $f({\omega }_n^k)$ 的值

DFT

此时 DFT 的分治思想就是分开考虑奇次项和偶次项：

考虑
$$f(x)=a_0{x}^0+a_1{x}^1+a_2{x}^2+\cdots$$
将其分为两个多项式
$$\begin{array}{rl}f(x)& =a_0{x}^0+a_2{x}^2+a_4{x}^4+a_6{x}^6+a_8{x}^8+\cdots +a_1{x}^1+a_3{x}^3+a_5{x}^5+a_7{x}^7+a_9{x}^9+\cdots \\ & =a_0{x}^0+a_2{x}^2+a_4{x}^4+a_6{x}^6+a_8{x}^8+\cdots +x(a_1{x}^0+a_3{x}^2+a_5{x}^4+a_7{x}^6+\cdots )\end{array}$$
考虑两个新多项式：
$$\begin{array}{rl}{f}_0(x)& =a_0{x}^0+a_2{x}^1+a_4{x}^2+a_6{x}^3+\cdots \\ f_1(x)& =a_1{x}^0+a_3{x}^1+a_5{x}^2+a_7{x}^3+\cdots \end{array}$$
不难发现
$$f(x)=f_0(x^2)+x{f}_1(x^2)$$
利用单位复数根的性质：
$$\begin{array}{rl}\mathrm{D}\mathrm{F}\mathrm{T}(f({\omega }_n^k))& =\mathrm{D}\mathrm{F}\mathrm{T}(f_0({\omega }_n^{2k}))+{\omega }_n^k\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_n^{2k}))\\ & =\mathrm{D}\mathrm{F}\mathrm{T}(f_0({\omega }_{\frac{n}{2}}^k))+{\omega }_n^k\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_{\frac{n}{2}}^k))\end{array}$$

$$\begin{array}{rl}\mathrm{D}\mathrm{F}\mathrm{T}(f({\omega }_n^{k+\frac{n}{2}}))& =\mathrm{D}\mathrm{F}\mathrm{T}(f_0({\omega }_n^{2k+n}))+{\omega }_n^{k+\frac{n}{2}}\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_n^{2k+n}))\\ & =\mathrm{D}\mathrm{F}\mathrm{T}(f_0({\omega }_n^n{\omega }_n^{2k}))-{\omega }_n^k\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_n^n{\omega }_n^{2k}))\\ & =\mathrm{D}\mathrm{F}\mathrm{T}(f_0({\omega }_{\frac{n}{2}}^k))-{\omega }_n^k\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_{\frac{n}{2}}^k))\end{array}$$

其中 $k<\frac{n}{2}$。不难发现只要我们求得出 $\mathrm{D}\mathrm{F}\mathrm{T}(f_0({\omega }_{\frac{n}{2}}^k))$ 与 $\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_{\frac{n}{2}}^k))$ 的话，就可以同时求出 $\mathrm{D}\mathrm{F}\mathrm{T}(f({\omega }_n^k))$ 和 $\mathrm{D}\mathrm{F}\mathrm{T}(f({\omega }_n^{k+\frac{n}{2}}))$。接下来再对 $f_0$ 与 $f_1$ 递归 DFT 即可。其时间复杂度函数是形如下面这样的：
$$T(n)=T(n/2)+O(n)$$
所以总复杂度为 $\Theta (n\log n)$

实际实现的时候一定要注意传进去的系数一定要是 $2^m$ 个的，不然分治的过程中左右不一样会出问题。第一次传进去的时候就高位补 $0$，补成最高项次数为 $2^{m-1}$ 的多项式。

void dft(int lim, complex *a)
{
    if (lim == 1) return;//常数项直接返回
    complex a1[lim >> 1], a2[lim >> 1];
    for (int i = 0; i < lim; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];//把系数按照奇偶分开
    dft(lim >> 1, a1, type);//求 DFT(f_0())
    dft(lim >> 1, a2, type);//求 DFT(f_1())
    complex Wn = complex(cos(2.0 * pi / lim), sin(2.0 * pi / lim)), w = complex(1, 0);
    for (int i = 0; i < (lim >> 1); ++i, w = w * Wn)
    {
        a[i] = a1[i] + w * a2[i];//求 DFT(f(\omega_n^k))
        a[i + (lim >> 1)] = a1[i] - w * a2[i];//求 DFT(f(\omega_n^{k+\fracn2}))
    }
    return;
}

IDFT

好了现在假装我们已经求出了两个多项式的点值表达并已经将他们乘起来，但是我们最终还是要把他还原回去到系数表示的。这个过程就叫做 IDFT。

其实就是我们需要求解下面关于 $a$ 的线性方程组：

$$\left[\begin{array}{ccccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& ({\omega }_n^0{)}^2& \cdots & ({\omega }_n^0{)}^n\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& ({\omega }_n^1{)}^2& \cdots & ({\omega }_n^1{)}^n\\ ⋮& ⋮& ⋮& & ⋮\\ ({\omega }_n^n{)}^0& ({\omega }_n^n{)}^1& ({\omega }_n^n{)}^2& \cdots & ({\omega }_n^n{)}^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]$$

我们将其乘上左边矩阵的逆：

$$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]={\left[\begin{array}{ccccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& ({\omega }_n^0{)}^2& \cdots & ({\omega }_n^0{)}^n\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& ({\omega }_n^1{)}^2& \cdots & ({\omega }_n^1{)}^n\\ ⋮& ⋮& ⋮& & ⋮\\ ({\omega }_n^n{)}^0& ({\omega }_n^n{)}^1& ({\omega }_n^n{)}^2& \cdots & ({\omega }_n^n{)}^n\end{array}\right]}^{-1}\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]$$

模相同的正交列向量构成的矩阵的逆是转置的模分之一倍，所以：

$${\left[\begin{array}{ccccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& ({\omega }_n^0{)}^2& \cdots & ({\omega }_n^0{)}^n\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& ({\omega }_n^1{)}^2& \cdots & ({\omega }_n^1{)}^n\\ ⋮& ⋮& ⋮& & ⋮\\ ({\omega }_n^n{)}^0& ({\omega }_n^n{)}^1& ({\omega }_n^n{)}^2& \cdots & ({\omega }_n^n{)}^n\end{array}\right]}^{-1}=\frac{1}{n+1}\left[\begin{array}{ccccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& ({\omega }_n^{-0}{)}^2& \cdots & ({\omega }_n^{-0}{)}^n\\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& ({\omega }_n^{-1}{)}^2& \cdots & ({\omega }_n^{-1}{)}^n\\ ⋮& ⋮& ⋮& & ⋮\\ ({\omega }_n^{-n}{)}^0& ({\omega }_n^{-n}{)}^1& ({\omega }_n^{-n}{)}^2& \cdots & ({\omega }_n^{-n}{)}^n\end{array}\right]$$

所以不难发现，IDFT 其实就是再做了一遍 DFT，只不过是反起来的。只是算出来最后的系数结果都要除以点值的个数，反应在代码里面就是那个 lim 变量。

不难发现 ${\omega }_n^k$ 的共轭就是虚部取反，所以可以在 DFT 函数里面传一个参数表示是否为 IDFT。

这样子一个递归版的 FFT 就写完了，总体的代码如下：

#include <cstdio>
#include <cctype>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

const int maxn = 2e6 + 5;
const double pi = acos(-1.0);

inline int read()
{
    char c = getchar();
    int s = 0;
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return s;
}

struct complex
{
    double x, y;
    complex(double xx = 0, double yy = 0)
    {
        x = xx, y = yy;
    }
} a[maxn], b[maxn];

complex operator+(const complex &a, const complex &b) {return complex(a.x + b.x, a.y + b.y);}
complex operator-(const complex &a, const complex &b) {return complex(a.x - b.x, a.y - b.y);}
complex operator*(const complex &a, const complex &b) {return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void dft(int lim, complex *a, int type)//type = 1 DFT；type = -1 IDFT
{
    if (lim == 1) return;//返回常数项
    complex a1[lim >> 1], a2[lim >> 1];
    for (int i = 0; i < lim; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    dft(lim >> 1, a1, type);
    dft(lim >> 1, a2, type);
    complex Wn = complex(cos(2.0 * pi / lim), type * sin(2.0 * pi / lim)), w = complex(1, 0);
    for (int i = 0; i < (lim >> 1); ++i, w = w * Wn)
    {
        a[i] = a1[i] + w * a2[i];
        a[i + (lim >> 1)] = a1[i] - w * a2[i];
    }
    return;
}

int main()
{
    int n = read(), m = read();
    FOR(i, 0, n) a[i].x = read();
    FOR(i, 0, m) b[i].x = read();
    int lim = 1;
    while (lim <= n + m) lim <<= 1;//lim一定要大于 n + m
    dft(lim, a, 1);
    dft(lim, b, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i];//点值乘起来
    dft(lim, a, -1);//IDFT还回去
    FOR(i, 0, n + m)
        printf("%d ", (int)(a[i].x / lim + 0.5));//最后要除那个数然后还原回去，四舍五入
    return 0;
}

位逆序置换

然而，上面的代码连模板都跑不过去……

考虑继续优化 DFT 的过程。递归的过程中开了大量的空间并且常数巨大，考虑非递归写法。

只考虑我们对 $0$ 到 $7$ 操作：

递归的过程：

original		0	1	2	3	4	5	6	7
recursion#1		0	2	4	6	1	3	5	7
recursion#2		0	4	2	6	1	5	3	7
recursion#3		0	4	2	6	1	5	3	7
original bin	000	001	010	011	100	101	110	111
now bin			000	100	010	110	001	101	011	111

可见递归到最后的结果无非就是一个二进制反转。

所以我们可以考虑非递归，一开始就先把所有的数放到最后的位置，然后迭代的时候一步步还回去即可。这个过程就是位逆序置换（蝴蝶变换）

考虑处理出 $x$ 二进制位翻转之后的数 $R(x)$。易知 $R(0)=0$。我们可以从小到大求 $R(x)$。很明显，$\lfloor x/2\rfloor$ 的二进制位是 $x$ 右移一位，那么如果知道了 $R(\lfloor x/2\rfloor )$ 就可以很容易的求出 $R(x)$，再分 $x$ 的奇偶性判断就可以了。
$$R(x)=\lfloor \frac{R(\lfloor x/2\rfloor )}{2}\rfloor +(x\mathrm{m}\mathrm{o}\mathrm{d}2)\times \frac{len}{2}$$
举个例子：翻转 $(10101110{)}_2$，首先我们知道它的二分之一倍为 $(01010111{)}_2$，其翻转结果为 $(11101010{)}_2$，除以二变为 $(01110101{)}_2$，由于它是偶数所以前面不用补 $1$。不难发现其就是一开始要求的翻转结果。

预处理翻转结果的代码：

while (lim <= n + m) lim <<= 1;
FOR(i, 0, lim - 1)
    rev[i] = ((rev[i >> 1] >> 1) | (((i & 1) ? (lim >> 1) : 0)));

然后在处理翻转的时候只需要下面几行：

FOR(i, 0, lim - 1)
    if (i < rev[i])
        myswap(a[i], a[rev[i]]);

不难验证其正确性。

而且观察我们在求 $\mathrm{D}\mathrm{F}\mathrm{T}(f({\omega }_n^k))$ 时我们需要算两遍 ${\omega }_n^k\mathrm{D}\mathrm{F}\mathrm{T}(f_1({\omega }_{\frac{n}{2}}^k))$，复数的乘法常数很大，考虑使用临时变量记录以降低常数。

这样子的话迭代版的 DFT 过程就很好写了：

void DFT(int lim, complex *a, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(a[i], a[rev[i]]);//先预处理翻转完了的结果
    for (int p = 2; p <= lim; p <<= 1)//模拟合并答案的过程，即为所谓的 n
    {
        int len = p >> 1;//即上面的 n / 2
        complex Wp = complex(cos(2 * pi / p), type * sin(2 * pi / p));//处理出 p 次单位根
        for (int k = 0; k < lim; k += p)//对每一个进行合并
        {
            complex w = complex(1, 0);//处理 \omega_p^0
            for (int l = k; l < k + len; ++l, w = w * Wp)//开始合并
            {
                //此时的 a[l] 就是之前的 a1[i]，a[len + l] 就是之前的 a2[i]
                complex tmp = w * a[len + l];
                a[len + l] = a[l] - tmp;//相当于上面的 a[i + (lim >> 1)] = a1[i] - w * a2[i]
                a[l] = a[l] + tmp;//相当于上面的 a[i] = a1[i] + w * a2[i]
            }
        }
    }
}

多项式乘法的实现

总的一个非递归版 FFT 的实现如下（洛谷 P3803）：

#include <cstdio>
#include <cctype>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

const int maxn = 3e6 + 5;
const double pi = acos(-1.0);

inline int read()
{
    char c = getchar();
    int s = 0;
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return s;
}

template<typename T> inline void myswap(T &a, T &b)
{
    T t = a;
    a = b;
    b = t;
    return;
}

struct complex
{
    double x, y;
    complex(double xx = 0, double yy = 0)
    {
        x = xx, y = yy;
    }
} a[maxn], b[maxn];

int rev[maxn];

complex operator+(const complex &a, const complex &b) {return complex(a.x + b.x, a.y + b.y);}
complex operator-(const complex &a, const complex &b) {return complex(a.x - b.x, a.y - b.y);}
complex operator*(const complex &a, const complex &b) {return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void DFT(int lim, complex *a, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(a[i], a[rev[i]]);//先预处理翻转完了的结果
    for (int p = 2; p <= lim; p <<= 1)//模拟合并答案的过程，即为所谓的 n
    {
        int len = p >> 1;//即上面的 n / 2
        complex Wp = complex(cos(2 * pi / p), type * sin(2 * pi / p));//处理出 p 次单位根
        for (int k = 0; k < lim; k += p)//对每一个进行合并
        {
            complex w = complex(1, 0);//处理 \omega_p^0
            for (int l = k; l < k + len; ++l, w = w * Wp)//开始合并
            {
                //此时的 a[l] 就是之前的 a1[i]，a[len + l] 就是之前的 a2[i]
                complex tmp = w * a[len + l];
                a[len + l] = a[l] - tmp;//相当于上面的 a[i + (lim >> 1)] = a1[i] - w * a2[i]
                a[l] = a[l] + tmp;//相当于上面的 a[i] = a1[i] + w * a2[i]
            }
        }
    }
}

int main()
{
    int n = read(), m = read();
    FOR(i, 0, n) a[i].x = read();
    FOR(i, 0, m) b[i].x = read();
    int lim = 1;
    while (lim <= n + m) lim <<= 1;//补齐高位
    FOR(i, 0, lim - 1)
        rev[i] = ((rev[i >> 1] >> 1) | (((i & 1) ? (lim >> 1) : 0)));//先处理翻转完的结果
    DFT(lim, a, 1);//DFT
    DFT(lim, b, 1);//DFT
    FOR(i, 0, lim)
        a[i] = a[i] * b[i];//对处理出来的点值进行乘法
    DFT(lim, a, -1);//IDFT
    FOR(i, 0, n + m)
        printf("%d ", (int)(a[i].x / lim + 0.5));
    return 0;
}

使用 FFT 来求高精度整数乘法的实现（洛谷 P1919）：

#include <cstdio>
#include <cstring>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)
#define DEC(i, a, b) for (int i = a; i >= b; --i)

template<typename T> inline void myswap(T &a, T &b) {T t = a; a = b; b = t; return;}

typedef double db;

const int maxn = 3000000 + 5;
const db pi = acos(-1.0);

struct cmplx
{
    db x, y;
    cmplx(db xx = 0, db yy = 0) {x = xx, y = yy;}
} a[maxn], b[maxn];

cmplx operator+(const cmplx &a, const cmplx &b) {return cmplx(a.x + b.x, a.y + b.y);}
cmplx operator-(const cmplx &a, const cmplx &b) {return cmplx(a.x - b.x, a.y - b.y);}
cmplx operator*(const cmplx &a, const cmplx &b) {return cmplx(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

char s1[maxn], s2[maxn];
int rev[maxn], ans[maxn];

void DFT(cmplx *f, int lim, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(f[i], f[rev[i]]);
    for (int p = 2; p <= lim; p <<= 1)
    {
        int len = p >> 1;
        cmplx Wp(cos(2.0 * pi / p), type * sin(2.0 * pi / p));
        for (int k = 0; k < lim; k += p)
        {
            cmplx w(1, 0);
            for (int l = k; l < k + len; ++l, w = w * Wp)
            {
                cmplx tmp = w * f[l + len];
                f[l + len] = f[l] - tmp;
                f[l] = f[l] + tmp;
            }
        }
    }
    return;
}

int main()
{
    scanf("%s\n%s", s1, s2);
    int n1 = -1, n2 = -1;
    DEC(i, strlen(s1) - 1, 0) a[++n1].x = s1[i] - '0';
    DEC(i, strlen(s2) - 1, 0) b[++n2].x = s2[i] - '0';
    int lim = 1;
    while (lim <= n1 + n2) lim <<= 1;
    FOR(i, 0, lim - 1)
        rev[i] = ((rev[i >> 1] >> 1) | (((i & 1) ? (lim >> 1) : 0)));
    DFT(a, lim, 1);
    DFT(b, lim, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i];
    DFT(a, lim, -1);
    FOR(i, 0, lim)
        ans[i] = (int)(a[i].x / lim + 0.5);
    FOR(i, 0, lim)
        if (ans[i] >= 10) ans[i + 1] += ans[i] / 10, ans[i] %= 10, lim += (i == lim);
    while (!ans[lim] && lim > -1) --lim;
    if (lim == -1) puts("0");
    else DEC(i, lim, 0) printf("%d", ans[i]);
    return 0;
}

当然，千万要记得 IDFT 还回去的时候要除以 lim，实在怕记不住就在 DFT 函数里面加几句话直接处理好

if (type == -1)
    FOR(i, 0, lim - 1)
        f[i].x /= lim;

针对多项式乘法：三次变两次优化

我们发现我们在做多项式乘法的时候，需要先 DFT $A(x)$，$B(x)$，乘在一起之后再 IDFT 还回来 $C(x)$，一共进行了三次这样的操作。考虑如何减少我们调用 DFT 的次数。

可以把 $B(x)$ 的系数放到 $A(x)$ 系数的虚部上面，即 $a+b\mathrm{i}$，然后 DFT 一下 $A(x)$ 再求个平方，得到 $A^2(x)$，再 IDFT 回去。我们可以发现得到的系数都是 $(a+b\mathrm{i}{)}^2=a^2-b^2+2ab\mathrm{i}$ 的形式的，所以只需要取出虚部再除以二就得到答案了。

这样的写法可以减小常数，跑的比 NTT 还快。

#include <cstdio>
#include <cctype>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

typedef double db;

const int maxn = 3e6 + 5;
const db pi = acos(-1.0);

inline int read()
{
    char c = getchar();
    int s = 0;
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return s;
}

template<typename T> inline void myswap(T &a, T &b)
{
    T t = a;
    a = b;
    b = t;
    return;
}

struct cmplx
{
    db x, y;
    cmplx(db xx = 0, db yy = 0)
    {
        x = xx, y = yy;
    }
} a[maxn];

int rev[maxn];

cmplx operator+(const cmplx &a, const cmplx &b) {return cmplx(a.x + b.x, a.y + b.y);}
cmplx operator-(const cmplx &a, const cmplx &b) {return cmplx(a.x - b.x, a.y - b.y);}
cmplx operator*(const cmplx &a, const cmplx &b) {return cmplx(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void DFT(cmplx *f, int lim, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(f[i], f[rev[i]]);
    for (int p = 2; p <= lim; p <<= 1)
    {
        int len = p >> 1;
        cmplx Wp(cos(2 * pi / p), type * sin(2 * pi / p));
        for (int k = 0; k < lim; k += p)
        {
            cmplx w(1, 0);
            for (int l = k; l < k + len; ++l, w = w * Wp)
            {
                cmplx tmp = w * f[len + l];
                f[len + l] = f[l] - tmp;
                f[l] = f[l] + tmp;
            }
        }
    }
}

int main()
{
    int n = read(), m = read();
    FOR(i, 0, n) a[i].x = read();
    FOR(i, 0, m) a[i].y = read();
    int lim = 1;
    while (lim <= n + m) lim <<= 1;
    FOR(i, 0, lim - 1)
        rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0));
    DFT(a, lim, 1);
    FOR(i, 0, lim - 1)
        a[i] = a[i] * a[i];
    DFT(a, lim, -1);
    FOR(i, 0, n + m)
        printf("%d ", (int)(a[i].y / lim / 2.0 + 0.5));
    return 0;
}

快速数论变换（NTT）

有了 FFT，我们已经有能力在 $O(n\log n)$ 的时间内求出两个多项式的卷积了。但是 FFT 也有它的缺点：复数采用的浮点运算不仅造成精度的问题，还会增大常数。遗憾的是数学家们已经证明了 $\mathbb{C}$ 中只有单位复数根满足 FFT 的要求。

考虑到利用多项式的计数题很多都是模意义下的，所以自然希望为单位复数根找一个模意义下的替代品。此时就进入下面的前置知识：原根。

原根

设整数 $r,n$ 满足 $r\perp n\wedge r\ne 0\wedge n>0$，使得 $r^x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$ 的最小正整数 $x$ 称为 $r$ 模 $n$ 的阶，记为 ${\mathrm{o}\mathrm{r}\mathrm{d}}_{n}r$ 或 ${\delta }_n(r)$。

若 $r,n\in {\mathbb{N}}^{+}\wedge r\perp n$，当 ${\mathrm{ord}}_{n}r=\varphi (n)$ 时，称 $r$ 是模 $n$ 的原根或者 $n$ 的原根。

NTT

对于质数 $p=qn+1(n=2^m)$，原根 $g$ 满足 $g^{qn}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，将 $g_n=g^q(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 看作 ${\omega }_n$ 的等价，其满足相似的性质，比如 $g_n^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p),g_n^{n/2}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

常见的质数
$$\begin{array}{rlr}p& =998244353=7\times 17\times 2^{23}+1,& g=3\\ p& =1004535809=479\times 2^{21}+1,& g=3\end{array}$$
迭代到长度为 $l$ 时，$g_l=g^{\frac{p-1}{l}}$

直接看代码：

#include <cstdio>
#include <cctype>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

typedef long long ll;

const ll G = 3;
const ll mod = 998244353;
const int maxn = 3e6 + 5;

inline int read()
{
    char c = getchar();
    int s = 0;
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return s;
}

template<typename T> inline void myswap(T &a, T &b)
{
    T t = a;
    a = b;
    b = t;
    return;
}

ll pow(ll base, ll p = mod - 2)
{
    ll ret = 1;
    for (; p; p >>= 1)
    {
        if (p & 1)
            ret = ret * base % mod;
        base = base * base % mod;
    }
    return ret;
}

int rev[maxn];
ll f[maxn], g[maxn];
const ll invG = pow(G);

void NTT(ll *f, int lim, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(f[i], f[rev[i]]);
    for (int p = 2; p <= lim; p <<= 1)
    {
        int len = p >> 1;
        ll tG = pow(type ? G : invG, (mod - 1) / p);
        for (int k = 0; k < lim; k += p)
        {
            ll buf = 1;
            for (int l = k; l < k + len; ++l, buf = buf * tG % mod)
            {
                ll tmp = buf * f[len + l] % mod;
                f[len + l] = f[l] - tmp;
                if (f[len + l] < 0) f[len + l] += mod;//及时取模
                f[l] = f[l] + tmp;
                if (f[l] > mod) f[l] -= mod;//及时取模
            }
        }
    }
    ll invlim = pow(lim);//最后还回去，除以lim相当于乘上lim的逆元
    if (!type)
        FOR(i, 0, lim - 1)
            f[i] = (f[i] * invlim % mod);
    return;
}

int main()
{
    int n = read(), m = read();
    FOR(i, 0, n) f[i] = read();
    FOR(i, 0, m) g[i] = read();
    int lim = 1;
    while (lim <= n + m) lim <<= 1;
    FOR(i, 0, lim - 1)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
    NTT(f, lim, 1), NTT(g, lim, 1);
    FOR(i, 0, lim - 1)
        f[i] = f[i] * g[i] % mod;
    NTT(f, lim, 0);
    FOR(i, 0, n + m)
        printf("%d ", (int)f[i]);
    return 0;
}

FFT/NTT 优化卷积的一些例子

在继续之前，我们先来看看 FFT/NTT 的一些应用。（高精度乘法就不说了，记得最后进位就可以了）

• 优化一般的卷积
• 和生成函数一起食用
• 字符串匹配（你没看错）

洛谷 P3338 [ZJOI2014]力

题意：给定 { q } \{q\} {q}，定义
$$F_i=\sum _{j=1}^{i-1}\frac{q_i{q}_j}{(i-j{)}^2}-\sum _{j=i+1}^n\frac{q_i{q}_j}{(i-j{)}^2}$$
求
$$E_i=\frac{F_i}{q_i}$$
考虑暴力的话，这道题是 $O(n^2)$ 的，过不去，考虑转化式子：
$$\begin{array}{rl}{E}_i& =\frac{F_i}{q_i}\\ & =\sum _{j=1}^{i-1}\frac{q_j}{(i-j{)}^2}-\sum _{j=i+1}^n\frac{q_j}{(i-j{)}^2}\end{array}$$
我们尝试将其化为卷积的形式，令 $f_i=q_i$，且 $f_0=0$；$g_i=\frac{1}{i^2}$，且 $g_0=0$，回代：
$$\begin{array}{rl}{E}_i& =\sum _{j=0}^i{f}_j{g}_{i-j}-\sum _{j=i}^n{f}_j{g}_{j-i}\end{array}$$
左边的部分已经是一个卷积的形式了，考虑继续化简右边。此时我们可以使用一个翻转的技巧，令 $f_i^{\prime }=f_{n-i}$，$t=n-i$，则右半边的式子可以化为 $\sum _{j=0}^t{f}_{t-j}^{\prime }{g}_j$。现在两边都化为卷积的形式了，可以愉快的使用 FFT 加速了。

即我们设多项式 $A(x)=\sum _{i=0}^n{f}_i{x}^n$，$B(x)=\sum _{i=0}^n{g}_i{x}^n$，$C(x)=\sum _{i=0}^n{f}_i^{\prime }$。再令 $L(x)=A(x)B(x)$，$R(x)=B(x)C(x)$，不难发现答案 $E_i=l_i-r_{n-i}$，其中 $l_i$ 和 $r_i$ 分别为 $L(x)$ 和 $R(x)$ 中 $x^i$ 的系数。

int main()
{
    int n; scanf("%d", &n);
    FOR(i, 1, n)
    {
        scanf("%lf", &a[i].x);
        b[i].x = (1.0 / i / i);
        c[n - i].x = a[i].x;
    }
    int lim = 1;
    while (lim <= (n << 1)) lim <<= 1;
    FOR(i, 0, lim)
        rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0));
    DFT(a, lim, 1), DFT(b, lim, 1), DFT(c, lim, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i], c[i] = b[i] * c[i];
    DFT(a, lim, -1), DFT(c, lim, -1);
    FOR(i, 1, n)
        printf("%.3lf\n", a[i].x - c[n - i].x);
    return 0;
}

洛谷 P3723 [AH2017/HNOI2017]礼物

题意：给定两个序列 { x } \{x\} {x} 和 { y } \{y\} {y}，可以整体平移序列或者整体加/减某个数，求最终序列

$$\sum _{i=1}^n(x_i-y_i{)}^2$$

的最小值。

分析：设整体加减的数为 $c$（$c$ 可正可负），我们需要最小化的就是下面这个式子：

$$\sum _{i=1}^n(x_i-y_i+c{)}^2$$

展开上面的式子，由 $(x_i-y_i+c{)}^2=x_i^2+y_i^2+c^2+2x_{i}c-2y_{i}c-2x_i{y}_i$ 可以得到原式可化简为

$$\sum x_i^2+\sum y_i^2+n{c}^2+2c\sum x_i-2c\sum y_i-2\sum x_i{y}_i$$

（下标省略）

不难发现我们只需要最大化 $\sum x_i{y}_i$ 就可以啦。

设 $f_k$ 为旋转了 $k$ 个单位后 $\sum x_i{y}_i$ 的取值，先把 { x } \{x\} {x} 倍长一波，则

$$f_k=\sum _{i=1}^n{x}_{i+k}{y}_k$$

翻转 $y$：

$$f_k=\sum _{i=1}^n{x}_{i+k}{y}_{n-i+1}^{\prime }$$

考虑多项式 $f(t)={\sum }_{i=1}^n{x}_i{t}^i$，$g(t)={\sum }_{i=1}^n{y}_i{t}^i$，令 $h(t)=f(t)\ast g(t)$，不难发现其 $t^{n+k+1}$ 的系数即为 $f_k$。因此可以使用 FFT/NTT 将倍长过的 { x } \{x\} {x} 与翻转过的 { y } \{y\} {y} 卷起来，然后把结果从第 $n+1$ 到第 $2n$ 处找最值就可以了

NTT 的实现：

ll a[maxn], b[maxn];
ll suma, sumb, suma2, sumb2, n, m;
int rev[maxn];

void NTT(ll *f, int lim, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            swap(f[i], f[rev[i]]);
    for (int p = 2; p <= lim; p <<= 1)
    {
        int len = p >> 1;
        ll Gp = pow(type ? G : invG, (mod - 1) / p);
        for (int k = 0; k < lim; k += p)
        {
            ll buf = 1;
            for (int l = k; l < k + len; ++l, buf = buf * Gp % mod)
            {
                ll tmp = buf * f[l + len] % mod;
                f[l + len] = f[l] - tmp;
                if (f[l + len] < 0) f[l + len] += mod;
                f[l] = f[l] + tmp;
                if (f[l] > mod) f[l] -= mod;
            }
        }
    }
    ll invlim = pow(lim);
    if (!type)
        FOR(i, 0, lim - 1)
            f[i] = f[i] * invlim % mod;
    return;
}

int main()
{
    n = read(), m = read();
    FOR(i, 1, n)
        a[i] = a[i + n] = read(), suma += a[i], suma2 += a[i] * a[i];
    FOR(i, 1, n)
        b[n - i + 1] = read(), sumb += b[n - i + 1], sumb2 += b[n - i + 1] * b[n - i + 1];
    int lim = 1;
    while (lim <= 3 * n) lim <<= 1;
    FOR(i, 0, lim - 1)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
    NTT(a, lim, 1), NTT(b, lim, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i] % mod;//千万不要忘记取模
    NTT(a, lim, 0);
    ll ans = 1e18;
    FOR(i, 1, n)
        FOR(j, -m, m)
            ans = min(ans, suma2 + sumb2 + j * j * n + 2 * j * (suma - sumb) - 2 * a[i + n]);
    printf("%lld\n", ans);
    return 0;
}

BZOJ3771 Triple

题意：有 $n$ 把价值分别为 $a_i$ 的斧子，河神可能拿走 1 - 3 把，问每种可能的损失价值及其对应方案数。（不计顺序）

思路：这是一道生成函数的入门题。

考虑设出多项式 $A(x)$，其系数有 $A[a_i]=1$，代表选一把的。则你可能会觉得答案为 $A(x)+A^2(x)+A^3(x)$。但是这样是显然不对的。为什么？

因为这样的话同一个元素可能被选两次或三次，对于这种情况定义 $B(x)$ 和 $C(x)$ 满足 $B[2a_i]=1$ 和 $C[3a_i]=1$，代表同时选两次/三次的，减掉这些方案数就可以了。然后需要注意顺序问题：

选一把的答案为 $A(x)$，不难发现选两种的即为 $\frac{A^2(x)-B(x)}{2}$，选三种的比较麻烦：不能同时选两种一样的，即减去 $3A(x)B(x)$，，但是选三种同样的又会被多减两次，最后除以 $3!$ 去掉顺序问题，所以最终答案为：

$$A(x)+\frac{A^2(x)-B(x)}{2}+\frac{A^3(x)-3A(x)B(x)+2C(x)}{6}$$

生成函数的卷积使用 NTT 或 FFT 优化即可。注意此时 NTT 模数要取一个更大的质数。不知道为什么生成函数能这样对应的可以意会一下多项式卷积的定义式以及这些系数的组合意义。

ll f1[maxn], f2[maxn], f3[maxn], ans[maxn];
ll g[maxn], t[maxn];

int main()
{
    int n = read();
    while (n--)
    {
        int tmp = read();
        ++f1[tmp], ++g[tmp], ++ans[tmp];
        ++f2[tmp << 1], ++f3[tmp * 3];
    }
    int lim = 1;
    while (lim <= (40000 * 3 + 5)) lim <<= 1;
    FOR(i, 0, lim - 1)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
    NTT(f1, lim, 1), NTT(g, lim, 1);
    FOR(i, 0, lim - 1)
        g[i] = f1[i] * g[i] % mod;
    NTT(g, lim, 0);
    FOR(i, 0, lim - 1)
        ans[i] += (g[i] - f2[i]) / 2;
    NTT(g, lim, 1);
    FOR(i, 0, lim - 1)
        g[i] = f1[i] * g[i] % mod;
    NTT(g, lim, 0);
    NTT(f2, lim, 1);
    FOR(i, 0, lim - 1)
        f2[i] = f2[i] * f1[i] % mod;
    NTT(f2, lim, 0);
    FOR(i, 0, lim - 1)
    {
        ans[i] += (g[i] - 3 * f2[i] + 2 * f3[i]) / 6;
        if (ans[i]) printf("%d %lld\n", i, ans[i]);
    }
    return 0;
}

FFT/NTT 与字符串匹配

字符串下标从 $1$ 开始

最一般的情况

考虑文本串 $S$ 和模式串 $T$，串长 $n=|S|$，$m=|T|$，保证 $n\ge m$，现在需要找出 $T$ 在 $S$ 中出现的每个位置。直接跑 KMP 就可以了，但是这个不是要提的重点。考虑串 $S$ 的第 $i$ 个字符为 $S[i]$，那么匹配就可以写成 $S[i]-T[j]=0$，这个应该是比较好想的。

假设 $T$ 在 $S$ 的第 $i$ 位开始成功匹配，则我们有

$$\sum _{j=1}^m(S[i+j-1]-T[j]{)}^2=0$$

为了防止正负号相互抵消，所以需要平方。由于 $i+j-1+j$ 不是定值，不符合我们需要的卷积的形式，所以翻转一下 $T$ 让其变为 $T^{\prime }$：

$$\sum _{j=1}^m(S[i+j-1]-T[m-j+1]{)}^2=0$$

打开来我们就会发现

$$\sum _{j=1}^m(S[i+j-1{]}^2+T[j{]}^2-2S[i+j-1]T[m-j+1])=0$$

$i+j-1+m-j+1=i+m$，为定值。

所以上面的式子就可以化成

$$\sum _{j=1}^{m}S[i+j-1{]}^2+\sum _{j=1}^{m}T[j{]}^2-2\sum _{x+y=i+m}S[x]T[y]$$

第一项直接前缀和就可以解决，第二项常数，第三项使用 FFT/NTT。

更加好理解地，设 $f(x)={\sum }_{i+j=x+m}S[i]T[j]$，我们只需要求出这个 $f(x)$ 就可以了。

问题来了，这样难写复杂度高常数大全方位被 KMP 吊打的算法有什么存在的意义吗？对不起还真的有：

带通配符的字符串匹配

请看例题 洛谷 P4173 残缺的字符串。仍然是字符串匹配，但是每个串都有通配符，这个时候 KMP 就显得无能为力了。怎么办呢？好好思考一下两个字符如何才能匹配：

• 两个字符完全一样
• 其中至少一个为通配符

两者是逻辑或的关系，我们魔改一下上面的式子，不难发现我们只需要把通配符的值设为 $0$ 就可以解决了：

定义匹配函数 $F(x)$ 表示 $S$ 的第 $x$ 位开始和 $T$ 是否匹配，匹配的话 $F(x)=0$。

$$F(x)=\sum _{j=1}^m(S[i+j-1]-T[j]{)}^{2}S[i+j-1]T[j]$$

化简：

$$\begin{array}{rl}F(x)& =\sum _{j=1}^m(S[x+j-1]-T[j]{)}^{2}S[x+j-1]T[j]\\ & =\sum _{j=1}^m(S[x+j-1]-T^{\prime }[m-j+1]{)}^{2}S[x+j-1]T^{\prime }[m-j+1]\\ & =\sum _{j=1}^m(S[x+j-1{]}^2+T^{\prime }[m-j+1{]}^2-2S[x+j-1]T^{\prime }[m-j+1])S[x+j-1]T^{\prime }[m-j+1]\\ & =\sum _{j=1}^{m}S[x+j-1{]}^3{T}^{\prime }[m-j+1]+\sum _{j=1}^{m}S[x+j-1]T^{\prime }[m-j+1{]}^3-2\sum _{j=1}^{m}S[x+j-1{]}^2{T}^{\prime }[m-j+1{]}^2\\ & =\sum _{i+j=x+m}S[i{]}^3{T}^{\prime }[j]+\sum _{i+j=x+m}S[i]T^{\prime }[j{]}^3-2\sum _{i+j=x+m}S[i{]}^{2}T[j{]}^2\end{array}$$

于是问题就解决了，只需要用 NTT/FFT 计算出上面三项恶心的东西出来就 OK 了。一共进行 $7$ 次 NTT 即可。

需要注意的是最后枚举答案的时候只能枚举到 $n-m+1$ 处，否则只有 $35$ 分。

int f[maxn << 1], g[maxn << 1], f2[maxn << 1], g2[maxn << 1], f3[maxn << 1], g3[maxn << 1];
int ans[maxn << 1], vec[maxn << 1], tot;

int n, m;
char a[maxn], b[maxn];

int main()
{
    m = readInt(), n = readInt();
    scanf("%s", a + 1);
    scanf("%s", b + 1);
    FOR(i, 1, m)
    {
        g[i] = (a[m - i + 1] == '*') ? 0 : a[m - i + 1] - 'a' + 1;
        g2[i] = g[i] * g[i], g3[i] = g2[i] * g[i];
    }
    FOR(i, 1, n)
    {
        f[i] = (b[i] == '*') ? 0 : b[i] - 'a' + 1;
        f2[i] = f[i] * f[i], f3[i] = f2[i] * f[i];
    }
    int lim = 1;
    while (lim <= n + m) lim <<= 1;
    NTT(f, lim, 1), NTT(f2, lim, 1), NTT(f3, lim, 1);
    NTT(g, lim, 1), NTT(g2, lim, 1), NTT(g3, lim, 1);
    FOR(i, 0, lim - 1)
        ans[i] = (1ll * f3[i] * g[i] % mod + 1ll * f[i] * g3[i] % mod - 2ll * f2[i] * g2[i] % mod) % mod;
    NTT(ans, lim, 0);
    int cnt = 0;
    FOR(i, 1, n - m + 1)
        if (ans[i + m] == 0)
            ++cnt, vec[++tot] = i;
    printf("%d\n", cnt);
    FOR(i, 1, tot) printf("%d ", vec[i]);
    return 0;
}

另外一道例题

请看 CF528D Fuzzy Search。题意：字符串匹配，$1\le |T|\le |S|\le 2\times 10^5$，字符集只有 $\text{ATCG}$。$T$ 在 $S$ 中的第 $i$ 个位置出现当且仅当 $\forall j\in [1,T]$，$\exists p$ 使得 $|i+j-1-p|\le k\wedge S[p]=T[j]$。即偏移量不能超过 $k$。

注意到字符集很小，只有 $4$ 个字符，所以我们可以把字符串 01 化，分开考虑每个字母。比如 $S=\text{ATCGAA}$，$T=\text{ACAA}$，现在只考虑字母 $\text{A}$，把 $\text{A}$ 化成 $1$ 而其他的化为 $0$，则 $S=\text{100011}$，$T=\text{1011}$。然后假设 $k=1$，把能扩展的都往两边扩展，则 $S=\text{110111}$。

这个时候，我们就可以设匹配 $F(x,c)$ 表示 $T$ 从 $S$ 的第 $x$ 位开始字符 $c$ 能匹配的数量，最终答案为 A ( x ) = ∑ c ∈ { A,T,C,G } F ( x , c ) A(x) = \sum_{c\in\lbrace\texttt{A,T,C,G}\rbrace}F(x,c) A(x)=∑c∈{A,T,C,G}​F(x,c)，$T$ 匹配成功当且仅当 $A(x)=m$，即所有字符都匹配到了，否则失败。

不难发现 $F(x,c)={\sum }_{j=1}^{m}S[x+j-1]T[j]$，老套路翻转一下变为 $F(x,c)={\sum }_{j=1}^{m}S[x+j-1]T[m-j+1]={\sum }_{i+j=x+m}S[i]T[j]$。这就是喜闻乐见的卷积形式了。NTT 直接上就完了。

int f[maxn << 1], g[maxn << 1];
int ans[maxn << 1];

int n, m, k, lim = 1;
char s[maxn << 1], t[maxn << 1];

void proc(char c)
{
    clr(f, lim), clr(g, lim);
    for (int i = 1, lst = -1e9; i <= n; ++i)
    {
        if (s[i] == c)
            lst = i;
        if (i - lst <= k)
            f[i] = 1;
    }
    for (int i = n, lst = 1e9; i; --i)
    {
        if (s[i] == c)
            lst = i;
        if (lst - i <= k)
            f[i] = 1;
    }
    FOR(i, 1, m)
        g[i] = (t[m - i + 1] == c);
    NTT(f, lim, 1), NTT(g, lim, 1);
    FOR(i, 0, lim - 1)
        f[i] = 1ll * f[i] * g[i] % mod;
    NTT(f, lim, 0);
    FOR(i, 1, n)
        ans[i] += f[i + m];
}

int main()
{
    n = readInt(), m = readInt(), k = readInt();
    scanf("%s", s + 1);
    scanf("%s", t + 1);
    while (lim <= n + m) lim <<= 1;
    FOR(i, 0, 3) proc("ATCG"[i]);
    int cnt = 0;
    FOR(i, 1, n)
        cnt += (ans[i] == m);
    printf("%d\n", cnt);
    return 0;
}

关于封装

以后的全家桶会大量使用 NTT 等基础操作，考虑实现一个常数较小的封装：

首先是各种 #define，由于我们在进行多项式运算的时候需要考虑界的问题，因此一定要把超过界了的给清零，不然可能出现各种奇奇怪怪的问题比如多卷了之类的：

#define ll long long
#define ull unsigned ll
#define FOR(i, a, b) for (int i = a; i <= b; ++i)
#define clr(f, n) memset(f, 0, (sizeof(int)) * (n))
#define cpy(f, g, n) memcpy(f, g, (sizeof(int)) * (n))

memset 和 memcpy 的用法建议自己去查。

接下来是一些基本的东西：

const ll G = 3, mod = 998244353;
const int maxn = ((1 << 21) + 500);

ll qpow(ll base, ll p = mod - 2)
{
    ll ret = 1;
    for (; p; p >>= 1)
    {
        if (p & 1)
            ret = ret * base % mod;
        base = base * base % mod;
    }
    return ret;
}

const ll invG = qpow(G);

没什么说的

NTT 时需要用到的位逆序置换：

int rev[maxn << 1], revlim;

void get_rev(int lim)
{
    if (lim == revlim) return;
    revlim = lim;
    FOR(i, 0, lim - 1)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
    return;
}

这样子可以在需要进行很多次 NTT 的时候智能的求出对应需要的 rev。

NTT 和乘法的封装：关于 static 关键字相关的建议自己查一下。

void NTT(int *g, int n, int type)
{
    get_rev(n);
    static ull f[maxn << 1], w[maxn];
    w[0] = 1;
    FOR(i, 0, n - 1)
        f[i] = (((long long)mod << 5ll) + g[rev[i]]) % mod;//防止负数带来影响
    for (int l = 1; l < n; l <<= 1)
    {
        ull tmp = qpow(type ? G : invG, (mod - 1) / (l << 1));
        FOR(i, 1, l - 1) w[i] = w[i - 1] * tmp % mod;//预处理“单位根”
        for (int i = 0; i < n; i += (l << 1))
        {
            for (int j = 0; j < l; ++j)
            {
                ll tt = w[j] * f[i + j + l] % mod;
                f[i + j + l] = f[i + j] + mod - tt;
                f[i + j] += tt;
            }
        }
        if (l == (1 << 10))
            FOR(i, 0, n - 1) f[i] %= mod;
    }
    if (!type)
    {
        ull invn  = qpow(n);
        FOR(i, 0, n - 1)
            g[i] = f[i] % mod * invn % mod;
    }
    else FOR(i, 0, n - 1)
        g[i] = f[i] % mod;
    return;
}

void times(int *f, int *g, int len, int lim)//len 表示两个多项式的最高次数，lim 为最终需要的项数
{
    static int sav[maxn << 1];//临时变量
    int n = 1;
    while (n < (len << 1)) n <<= 1;
    clr(sav, n), cpy(sav, g, n);
    NTT(f, n, 1); NTT(sav, n, 1);
    FOR(i, 0, n - 1)
        f[i] = 1ll * f[i] * sav[i] % mod;
    NTT(f, n, 0);
    clr(f + lim, n - lim), clr(sav, n);//把界以上的部分清掉，把 sav 清干净
    return;
}

有了如上的封装，我们在写 P3803 时主函数里面就可以帅气的写：

int main()
{
    n = readInt(), m = readInt();
    FOR(i, 0, n) f[i] = readInt();
    FOR(i, 0, m) g[i] = readInt();
    times(f, g, max(m, n), m + n + 1);
    FOR(i, 0, m + n) printf("%d ", f[i]);
    return 0;
}

就做完了。

多项式乘法逆

定义

当两个多项式 $F(x)$ 与 $G(x)$ 在每一项系数模 $p$ 时有 $F(x)\ast G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$ 时，称 $F(x)$ 与 $G(x)$ 互为乘法逆元。此处 $(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$ 代表次数高于 $n$ 的项都不考虑。

需要的前置知识：NTT

求法

考虑倍增。假设我们要求满足 $F(x)\ast G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^k)$ 的 $G(x)$，并且已经求出了满足 $F(x)\ast G^{\prime }(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{x}{2}\rceil })$ 的 $G‘(x)$。那么我们由
$$F(x)\ast G^{\prime }(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{x}{2}\rceil })$$
必然可以推出
$$F(x)\ast G(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{x}{2}\rceil })$$
那么
$$G(x)\equiv G^{\prime }(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{x}{2}\rceil })$$
作差，
$$G(x)-G^{\prime }(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{x}{2}\rceil })$$
将两边同时平方，不难发现界会从 $\lceil \frac{n}{2}\rceil$ 变为 $n$。
$$G^2(x)-2G(x)G^{\prime }(x)+G^{\prime 2}(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
现在要求的是 $G(x)$，而我们发现 $G^2(x)$ 不好处理，怎么办？同时乘以 $F(x)$ 就可以消掉一个 $G(x)$：
$$G(x)-2G^{\prime }(x)+G^{\prime 2}(x)F(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
所以我们得到了
$$G(x)\equiv 2G^{\prime }(x)-G^{\prime 2}(x)F(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
根据这个，我们就可以从 $G^{\prime }(x)$ 推出 $G(x)$ 的值出来了。从上往下递归求解，到常数项的时候直接费马小定理求逆元然后一步步回溯上去。不难发现复杂度为 $T(n)=T(n/2)+O(n\log n)$，由主定理知总复杂度为 $O(n\log n)$。

实现

递归版多项式求逆：

int tmp[maxn << 1];//这是临时数组

void invpoly(int *f, int *ans, int m)
{
    if (m == 1)
        return ans[0] = qpow(f[0]), void();//常数项就直接返回
    invpoly(f, ans, (m + 1) >> 1);//先递归求出 m / 2 向上取整的情况
    int n = 1;
    while (n < (m << 1)) n <<= 1;
    cpy(tmp, f, m), clr(tmp + m, n - m);//把 f 数组的前 n 项都复制进 tmp 里面并把 tmp 高于 m 的地方全部清零
    NTT(ans, n, 1), NTT(tmp, n, 1);
    FOR(i, 0, n - 1)
        ans[i] = 1ll * (2ll - 1ll * ans[i] * tmp[i] % mod + mod) % mod * ans[i] % mod;//根据公式计算
    NTT(ans, n, 0);
    clr(ans + m, n - m);//高于 m 的舍弃
    return;
}

int f[maxn << 1], ans[maxn << 1];

int n, m;

int main()
{
    n = readInt();
    FOR(i, 0, n - 1) f[i] = readInt();
    invpoly(f, ans, n);
    FOR(i, 0, n - 1) printf("%d ", ans[i]);
    return 0;
}

当然我们也可以不递归，考虑递推实现（其实快不了多少的）。递推无非就是从 $1$ 开始倍增向上走。

void invpoly(int *f, int m)
{
    int n;
    for (n = 1; n < m; n <<= 1);
    static int w[maxn << 1], r[maxn << 1], sav[maxn << 1];
    w[0] = qpow(f[0]);
    for (int len = 2; len <= n; len <<= 1)//len 代表当前的界
    {
        FOR(i, 0, (len >> 1) - 1)
            r[i] = (w[i] << 1) % mod;//处理 2G'(x)
        cpy(sav, f, len);
        NTT(w, len << 1, 1);
        FOR(i, 0, (len << 1) - 1)
            w[i] = 1ll * w[i] * w[i] % mod;
        NTT(sav, len << 1, 1);
        FOR(i, 0, (len << 1) - 1)
            w[i] = 1ll * w[i] * sav[i] % mod;
        NTT(w, len << 1, 0);
        clr(w + len, len);
        FOR(i, 0, len - 1)
            w[i] = (r[i] - w[i] + mod) % mod;
    }
    cpy(f, w, m);//把答案还回到 f 数组里面
    clr(sav, n << 1);//清零
    clr(w, n << 1);//清零
    clr(r, n << 1);//清零
    return;
}

请注意在递推进行乘法的时候一定要把空间开成 len << 1，为什么呢？计算 $G^{\prime 2}(x)F(x)$ 的时候，$G^{\prime 2}(x)$ 含 len 项，$F(x)$ 也含 len 项，卷起来就是 len << 1 项了，如果不开够的话 NTT 把点值还原回去的时候是会出问题的。

多项式的导数/积分

回顾一下一些基础的导数公式：
$$\begin{array}{rl}f(x)=e^x& ⟹f^{\prime }(x)=e^x\\ f(x)=\ln x& ⟹f^{\prime }(x)=\frac{1}{x}\\ f(x)=a{x}^k& ⟹f^{\prime }(x)=ak{x}^{k-1}\\ (f(x)\pm g(x){)}^{\prime }& =f^{\prime }(x)\pm g^{\prime }(x)\\ (f(g(x)){)}^{\prime }& =f^{\prime }(g(x))\times g^{\prime }(x)\end{array}$$
我们在这里定义一下多项式的求导：
$$f(x)=\sum _{i=0}^n{a}_i{x}^i⟹f^{\prime }(x)=\sum _{i=0}^{n-1}(i+1)a_{i+1}{x}^i$$
积分为求导的逆运算：
$$f(x)=\sum _{i=0}^n{a}_i{x}_i⟹\int f(x)\mathrm{d}x=\sum _{i=1}^{n+1}\frac{a_{i-1}{x}^i}{i}$$
所以多项式求导和求积分的代码就很容易写出来了，当然需要一开始线性预处理一下逆元。

void derivate(int *f, int m)
{
    FOR(i, 1, m - 1)
        f[i - 1] = 1ll * f[i] * i % mod;
    f[m - 1] = 0;
    return;
}

int inv[maxn];

void initinv(int lim)
{
    inv[1] = 1;
    FOR(i, 2, lim)
        inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
    return;
}

void intergrate(int *f, int m)
{
    DEC(i, m, 1)
        f[i] = 1ll * f[i - 1] * inv[i] % mod;
    f[0] = 0;
}