多项式
i_love_xqh
·
·
算法·理论
拉格朗日插值
基本流程
对于一个 n 次函数 f(x)，给出 n+1 个点值 (x_i,y_i)，保证 x_i 互不相同，求 f(k)。
考虑 crt 的思想，对于第 i 个点值 (x_i,y_i)，构造函数 f_i(x) 满足 f_i(x_i)=y_i，且对于 j\ne i，f_i(x_j)=0。
由
f_i(x)=\lambda\prod_{j=0,j\ne i}^n(x-x_j)\\
f_i(x_i)=y_i
可以推出 f_i(x) 的表达式
f_i(x)=y_i\prod_{j=0,j\ne i}^{n}\frac{x-x_j}{x_i-x_j}
所以 f(x) 的表达式即为所有 f_i(x) 的和
f(x)=\sum_{i=0}^ny_i\prod_{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}
直接计算时间复杂度为 O(n^2)，用分治 NTT 可以优化到 O(n\log^2 n)。
特殊情况
如果 x_i=i，则 f(x) 可以变为
\begin{aligned}
f(x)&=\sum_{i=0}^{n}y_i\frac{x(x-1)\cdots(x-i+1)(x-i-1)(x-i-2)\cdots(x-n)}{i(i-1)\cdots1(-1)(-2)\cdots(i-n)}\\
&=\sum_{i=0}^{n}y_i\frac{p_is_i}{i!(n-i)!(-1)^{n-i}}
\end{aligned}
其中 p_i=\prod\limits_{j=0}^{i-1}(x-j)，s_i=\prod\limits_{j=i+1}^{n}(x-j)，通过预处理可以 O(n) 得到，再预处理出阶乘及逆元，整体时间复杂度就是 O(n)。
动态加点
因为
f(x)=\sum_{i=0}^ny_i\prod_{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}
令 g(x)=\prod\limits_{i=0}^n(x-x_i)，h(i)=\prod\limits_{j=0,j\ne i}^n\frac{y_i}{x_i-x_j}，所以有
f(x)=g(x)\sum_{i=0}^{n}h(i)
每次加点 (x_{n+1},y_{n+1}) 时对于 0\le i\le n，h(i) 都会乘以 \frac{1}{x_i-x_{n+1}}，然后 h(n+1) 直接算就行了，所以查询和修改操作都是 O(n) 的。
例题
[NOIP2020] 微信步数
先把无解判掉，即 n 步后每一维位移量为 0 且偏差量不超过 w_i。
转换问题，即枚举走的步数，统计有多少种方案不会走出去。进一步发现对于每一维，可供选择的位置是一段区间 [l_i,r_i]，所以方案数即 \prod(r_i-l_i+1)。可以通过统计位移量 s_i 的极值来计算。时间复杂度最大是 O(nkV)。
考虑优化，令 l_{i,j} 为第一轮经过 i 步第 j 维位移量的最小值，r_{i,j} 同理，d_j 表示第 j 维的 n 步后的位移量。相当对于 [w_j-r_{i,j},w_j] 和 [1,-l_{i,j}] 这两段区间都不能选。对于第二轮，那么 l'_{i,j} 要么就是 l_{i,j}，要么就是 l_{i,j}+d_j，r_{i,j} 同理。发现每一轮位移量是相同的，那么对于经过 i 步，每一轮减少的方案数也是相同的，设 p_{i,j}=r'_{i,j}-r_{i,j}+l_{i,j}-l'_{i,j} 即经过 i 步第 j 维减少的方案数，g_j 表示第一轮后第 j 维的方案数。所以我们可以将统计答案改写成以下式子：
\sum_{i=1}^{n}\sum_{x=0}^{\infin}\prod_{j=1}^{k}(g_j-p_{i,j}-p_{n,j}x)
因为第 x+2 轮走了 i 步，第 2\sim x+1 轮每一轮第 j 维会减少 p_{n,j} 种方案，所以就是第 x+2 轮经过 i 步后的方案数。容易发现第二维上限是可以算出的，后面又是关于 x 的 k 次多项式，所以可以在 k^2 的时间复杂度求出多项式系数，然后交换求和顺序变成求自然数幂和，即
\sum_{i=1}^{n}\sum_{j=0}^{k}h_j\sum_{x=0}^{\text{lim}}x^j
补充知识
复数
定义复数域 \mathbb{C}，对于 z\in \mathbb{C}，z 可以唯一表示为 a+b\mathrm{i}，其中 a,b\in \R，定义 \mathrm{i}^2=-1。
基本运算法则：
-
-
z_1z_2=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)\mathrm{i}
- 如果 z_1=a+b\mathrm{i}，z_2=a-b\mathrm{i}，则 z_1 和 z_2 互为共轭。
建立复平面，实轴为 x，虚轴为 y，则任意 z 在复平面上都可以表示成一个点坐标 Z(a,b)。
考虑用向量表示 z，即 (\theta,r)，其中 r=\lvert z\rvert=\sqrt{a^2+b^2}，又称为 z 的模长。\theta 为 \overrightarrow{OZ} 与 x 正半轴的夹角。
简单推导可以得到 z=(\theta,r)=r(\cos\theta+\mathrm{i}\sin\theta)，由欧拉公式 e^{\mathrm{i}\theta}=\cos \theta+\mathrm{i}\sin \theta，可以得到 z=r\cdot e^{\mathrm{i}\theta}。
于是可得：
单位根
定义 \omega_n 为 n 次单位根，即 \omega_n^{n}=1。
所以 \omega_n^n=1=e^{\mathrm{i}\cdot 0}，于是 \omega_n=e^{\frac{\mathrm{i}2k\pi}{n}}，为了方便，令 \omega_{n}=e^{\frac{\mathrm{i}2\pi}{n}}=\cos(\frac{2\pi}{n})+\mathrm{i}\sin(\frac{2\pi}{n})。
那么 n 个根可以表示为 \omega^{0}_{n},\omega^{1}_{n},\cdots,\omega^{n-1}_{n}。
对数
a^c=b\Longleftrightarrow \log_ab=c
\ln a=\log_ea
\log_{a}(bc)=\log_ab+\log_ac
\log_{a}(\frac{b}{c})=\log_ab-\log_ac
\log_{a}(b^n)=n\log_ab
\log_{a}b=\frac{\log_cb}{\log_ca}
极限
定义
设 f(x) 在 x_0 的邻域有定义，若对于 \forall \varepsilon >0，\exist \delta>0，使得当 0<|x-x_0|<\delta 时，有 |f(x)-a|<\varepsilon，则记 a=\lim\limits_{x\to x_0}f(x)。
运算法则
若 \lim\limits_{x\to x_0}f(x)=a，\lim\limits_{x\to x_0}g(x)=b，则
-
\lim\limits_{x\to x_0}(f+g)(x)=a+b
-
\lim\limits_{x\to x_0}(fg)(x)=ab
- 如果 b\ne 0，\lim\limits_{x\to x_0}(\frac{f}{g})=\frac{a}{b}
如果 g(x) 在 u 处有极限 a，f(x) 在 a 处有极限 b，且存在 (u-\delta,u)\cup (u,u+\delta) 使 g(x)\ne a。则
\lim\limits_{x\to u}(f\circ g)(x)=b
-
\lim[c\cdot f(x)]=c\cdot \lim f(x)
-
\lim[f(x)]^n=[\lim f(x)]^n
夹逼定理
设 g(x)\le f(x)\le h(x)，且 \lim\limits_{x\to x_0}g(x)=\lim\limits_{x\to x_0}h(x)=a，则 \lim\limits_{x\to x_0}f(x)=a。
可证得重要结论
\lim\limits_{x\to 0}\frac{\sin x}{x}=1
单调有界准则
设函数 f(x) 在 x_0 的某左邻域单调且有界，则其左极限 \lim_{x\to x_0^{-}}f(x) 必定存在。
类似还有右极限和趋于无穷的单调有界准则。
可证得重要结论
e=\lim\limits_{x\to \infty}\left(1+\frac{1}{x}\right)^{x}
存在。
导数
定义
设函数 f(x) 在 x_0 的某邻域 U(x_0) 内有定义，若极限
\lim_{x\to x_0}\frac{f(x)-f(x_0)}{x-x_0}
存在，则称 f(x) 在 x_0 处可导，该极限为 f(x) 在 x_0 处的导数，记为 f'(x_0)。也可改写为
f'(x_0)=\lim_{\Delta x\to 0}\frac{f(x_0+\Delta x)-f(x_0)}{\Delta x}
如果 f(x) 在区间 I 上处处可导，定义 f(x) 在区间上的导函数 f'(x) 为
f'(x)=\lim_{\Delta x\to 0}\frac{f(x+\Delta x)-f(x)}{\Delta x},x\in I
如果 y=f(x)，那么 f'(x)=y'=\frac{dy}{dx}=\frac{d}{dx}f(x)。
### 运算法则
若函数 $f(x)$ 和 $g(x)$ 可导，则
- $(f+g)'=f'+g'
-
(fg)'=f'g+fg'
- 如果在 g\ne 0 处，(\frac{f}{g})'=\frac{f'g-fg'}{g^2}
链式法则：若 u=g(x) 可导，v=f(u) 可导，则复合函数 f\circ g 可导，且 (f\circ g)'(x)=f'(u)g'(x)=f'(g(x))g'(x)。
指数对数化：指数和底数都含 x 的函数叫幂指函数，可通过指数对数化，f(x)^{g(x)}=e^{\ln f(x)^{g(x)}}=e^{g(x)\ln f(x)}，然后通过链式法则和乘法法则求导。
常用结论
-
-
-
- $y=\log_a x$，则 $a^{y}=x$。同时对 $x$ 求导得 $a^{y}\ln a\cdot y'=1$，移项即得 $y'=\frac{1}{x\ln a}$。
微分
定义
对于 y=f(x)，如果当 \Delta x\to 0 时，\Delta y 可表示为 \Delta y=k\Delta x+o(\Delta x)，其中 k 是常量，则称 f(x) 在 x_0 处可微，dy=k\Delta x，记为微分符号。
### 定理
罗尔定理、拉格朗日中值定理、柯西中值定理。
### 洛必达法则
若
- $\lim\limits_{x\to a}f(x)=\lim\limits_{x\to a}F(x)=0
-
f(x)$ 与 $F(x)$ 在 $a$ 的某去心邻域可导且 $F'(x)\ne 0
-
\lim\limits_{x\to a}\frac{f'(x)}{F'(x)}=A
则 \lim\limits_{x\to a}\frac{f(x)}{F(x)}=\lim\limits_{x\to a}\frac{f'(x)}{F'(x)}=A。
以上是 0 比 0 型，也有 \infty 比 \infty 型。
泰勒中值定理
若函数 f(x) 在包含 x_0 在内的开区间 (a,b) 内 (n+1) 阶可导，则对 \forall x\in(a,b) 有
f(x)=f(x_0)+\sum_{k=1}^{n}\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k+R_n(x)
其中（\xi 为 x 与 x_0 之间的某个值）
R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-x_0)^{n+1}
上式称为 f(x) 在 x_0 处展开的 n 阶泰勒公式，其中 R_n(x) 称为拉格朗日余项。
特别地，如果 f(x) 无穷阶可导，则它的无穷阶泰勒公式
\sum_{k=0}^{\infty}\frac{f^{(n)}(x_0)}{k!}(x-x_0)^k
称为它的泰勒级数。
在泰勒公式中，取 x_0=0，得到麦克劳林公式
f(x)=f(0)+\sum_{k=1}^{n}\frac{f^{(k)}(0)}{k!}x^k+o(x^n)
常见函数的麦克劳林展开：
-
e^x=1+x+\frac{x^2}{2!}+\cdots+\frac{x^n}{n!}+\cdots
-
\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}\cdots+(-1)^{n-1}\frac{x^{2n-1}}{(2n-1)!}+\cdots
-
\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}+\cdots +(-1)^n\frac{x^{2n}}{(2n)!}$+$\cdots
积分
本质是求导的逆运算。
不定积分
若 F'(x)=f(x) 或 d(F(x))=f(x)dx，则 F(x) 为 f(x) 的一个原函数。将 f(x) 的原函数集合（原函数加常数仍是原函数）称为它的不定积分，记
\int f(x)dx=F(x)+C
不定积分表（部分）
\int x^{a}dx=\frac{x^{a+1}}{a+1}+C(a\ne -1)\\
\int \frac{dx}{x}=\ln |x|+C
第一类换元积分法（凑微分）
\begin{aligned}
\int f(x)dx&=\int g(\varphi(x))\varphi'(x)dx\\
&=\int g(u)\frac{du}{dx}dx\Big|_{u=\varphi(x)}\\
&=\int g(u)du\\
&=G(u)+C
\end{aligned}
关键在于把 f(x) 凑成 g(u)u'。
如
\begin{aligned}
\int 2xe^{x^2}dx&=\int u'e^{u}dx\Big|_{u=x^2}\\
&=\int e^udu\\
&=e^u+C\\
&=e^{x^{2}}+C
\end{aligned}
第二类换元积分法
第一类反过来：
\begin{aligned}
\int g(u)du&=\int g(\varphi(x))\varphi'(x)dx\Big|_{\varphi (x)=u}\\
&=\int f(x)dx\\
&=F(x)+C\\
&=F(\varphi^{-1}(u))+C
\end{aligned}
例：
\begin{aligned}
\int \sqrt{a^2-x^2}dx&=\int \sqrt{a^2-(a\sin t)^2}(a\sin t)'dt\Big|_{a\sin t=x}\\
&=\int \sqrt{a^2(1-\sin^2 t)}a\cos tdt\\
&=a^2\int \cos^2tdt\\
&=\frac{a^2}{2}\arcsin \frac{x}{a}+\frac{x}{2}\sqrt{a^2-x^2}+C
\end{aligned}
还有三角代换、倒代换。
分部积分法
由 (uv)'=u'v+uv' 得
\int uv'dx=uv-\int u'vdx
可简写为
\int udv=uv-\int vdu
### 定积分
$$
\int_a^bf(x)dx=\lim_{\Delta x_i\to 0}\sum_{i=1}^{n}f(\xi_i)\Delta x_i
$$
可以看作 $[a,b]$ 区间围成图形的面积。
可以求弧长、体积、侧面积，在 OI 中可用来估计时间复杂度。
#### 牛顿-莱布尼茨公式
微积分基本定理。
$$
\int_a^b f(x)dx=F(x)\Big|_{b}^{a}=F(b)-F(a)
$$
其中 $F(x)$ 为 $f(x)$ 在 $[a,b]$ 上的原函数。
#### 斯特林公式
$n!\sim \sqrt{2\pi n}(\frac{n}{e})^n$，可以用来估计时间复杂度。
#### 数值积分
用于求定积分的近似值。
辛普森公式：
$$
\int_a^bf(x)dx\approx \frac{b-a}{b}\left(f(a)+4f(\frac{a+b}{2})+f(b)\right)
$$
以及自适应辛普森算法。
# 快速傅里叶变换
## FFT
对于多项式 $F(x)=\sum\limits_{i=0}^{n-1}a_ix^{i}$，$G(x)=\sum\limits_{i=0}^{m-1}b_ix^{i}$，求 $H=F*G$。
设 $H=\sum\limits_{i=0}^{n+m-2}c_ix^i$，由于这是关于 $x$ 的 $n+m-2$ 次函数。朴素的话有 $c_i=\sum\limits_{j=0}^{i}a_jb_{i-j}$。
考虑快速求出 $H$ 的 $n+m-1$ 个点值，再快速求出 $H$。
### DFT
因为 $H(x)=F(x)\cdot G(x)$，考虑选出 $N=n+m-1$ 个点 $x_0,\cdots,x_{N-1}$，考虑 $F$，有
$$
\begin{bmatrix}
x_0^{0} & x_0^{1} & \cdots & x_{0}^{N-1} \\
x_1^{0} & x_1^{1} & \cdots & x_{1}^{N-1} \\
\vdots & \vdots & \ddots & \vdots \\
x_{N-1}^{0} & x_{N-1}^{1} & \cdots & x_{N-1}^{N-1}
\end{bmatrix}
\begin{bmatrix}
a_0 \\
a_1 \\
\vdots \\
a_{N-1}
\end{bmatrix}=
\begin{bmatrix}
F(x_0) \\
F(x_1) \\
\vdots \\
F(x_{N-1})
\end{bmatrix}
$$
直接算时间复杂度是 $O(N^2)$ 的，考虑递归的思想。
通过神奇的观察可以发现，当 $x_i=\omega_{N}^i$ 时，左边矩阵不同的数只有 $N$ 个，接下来就好办了。
令 $x=\omega_n^{i}$。定义 $F_0(x)$ 和 $F_1(x)$。
$$
F(x)=a_0x^0+a_1x^1+a_2x^2+a_3x^3+\cdots+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\\
F_0(x)=a_0x^0+a_2x^1+\cdots+a_{n-2}x^{\frac{n-2}{2}}\\
F_1(x)=a_1x^0+a_3x^1+\cdots+a_{n-1}x^{\frac{n-2}{2}}
$$
于是有 $F(x)=F_0(x^2)+xF_1(x^2)$。令 $m=\frac{n}{2}$，容易发现 $(\omega_{n}^{i})^2=\omega_{m}^{i}$，所以发现变成子任务，然后分类讨论：
- $0\le i<m$，$F(\omega_{n}^{i})=F_0(\omega_{m}^i)+\omega_{n}^{i}F_1(\omega_{m}^{i})$。
- $0\le i<m$，$F(\omega_{n}^{i+m})=F(-\omega_{n}^{i})=F_0(\omega_{m}^{i})-\omega_{n}^{i}F_1(w_{m}^{i})
### IDFT
已知左边矩阵和答案矩阵求系数矩阵，容易发现即求左边矩阵的逆矩阵。
经过观察可以发现，令左边矩阵为 $X$，$X_{i,j}=\omega_{N}^{ij}$，则 $X^{-1}_{i,j}=\frac{1}{N}\omega_{N}^{-ij}$。通过直接计算可以证明。
### 实现细节
- 可以手写复数类，容易发现 $\omega_{N}^{i}$ 与 $\omega_{N}^{-i}$ 互为共轭，所以 DFT 和 IDFT 可以调用一个函数，传一个参进去即可。
- 可以去掉递归，只需要求出每个位置 $i$ 在底层时对应哪个位置 $i'$ 即可，通过观察可以发现 $i'$ 事实上是 $i$ 的二进制翻转。然后模拟递归过程，从底层开始，每次合并 $2^{k}$ 个数。
- 对于 $F_0$ 和 $F_1$，将 $F_0(\omega_{m}^{i})$ 的值存在 $i$ 的位置，$F_1(\omega_{m}^{i})$ 的值存在 $m+i$ 的位置，这样就不用新开空间了。
- 求单位根。如果已知 $0\sim m-1$ 的 $\omega_{m}^{i}$，那么可以直接算出 $\omega_{n}^{i}$，对于 $i\equiv 0\pmod 2$，$\omega_{n}^{i}=\omega_{m}^{i/2}$，否则只需要先通过欧拉公式求出 $\omega_{n}$，则 $\omega_{n}^{i}=\omega_{n}^{i-1}\omega_n$。
```cpp
void fft(comp a[],int n,int ty){
for(int i=0,j=0;i<n;i++){
if(i<j)swap(a[i],a[j]);
for(int k=n>>1;(j^=k)<k;k>>=1);
}
w[0]=comp{1,0};
for(int i=1;i<n;i<<=1){
int t=i*2;
comp x=comp{cos(2*pi/t),ty*sin(2*pi/t)};
for(int j=t-2;j>=0;j-=2)w[j]=w[j>>1],w[j+1]=w[j]*x;
for(int j=0;j<n;j+=t){
for(int k=j;k<j+i;k++){
comp t0=a[k],t1=a[k+i]*w[k-j];
a[k]=t0+t1;
a[k+i]=t0-t1;
}
}
}
if(ty==-1){
for(int i=0;i<n;i++)a[i].a/=n,a[i].b/=n;
}
}
```
## 例题
### 01序列问题
给出两个长度分别为 $n$、$m$ 的 01 序列 $A$、$B$，有 $Q$ 次询问，每次询问当把 $A$ 的第 $i$ 位和 $B$ 的第 $j$ 位对齐时，$A$,$B$ 公共部分有多少对对齐且相同的数。
---
不如将 $B$ 序列翻转得到 $B'$，则 $A\times B'$ 的 $i$ 次项就是 $\sum\limits_{j=0}^{i}A_{j}B'_{i-j}$，即将 $A$ 与 $B$ 的第 $m-i$ 位对齐都是 $1$ 的位置数。
然后分别将 01 取反再跑一次就得到都是 $0$ 的位置数。
### [[HNOI2017] 礼物](https://www.luogu.com.cn/problem/P3723)
如果将 $A$ 序列加 $x$（这里允许 $x$ 为负数，相当于 $B$ 序列加 $-x$），则答案为
$$
\sum_{i=0}^{n-1}(A_i+x-B_i)^2=S^2(A)+S^2(B)+2x(S(A)-S(B))+nx^2-2\sum_{i=0}^{n-1}A_iB_i
$$
容易发现前 $2$ 项为定值，中间 $2$ 项可以通过枚举 $x$ 算，只需要考虑求最后一项的最大值。
将 $B$ 翻转并复制一份即可。
# 快速数论变换
## NTT
对于 FFT，只能处理 $10^4\times 10^4\times 10^6$ 级别问题，其中 $10^4$ 是值域，$10^6$ 是长度。否则会爆 double。
如果题目中有模数且可以用 NTT，则能处理 $p\times p\times 10^6$ 级别问题。并且将复数运算代替为取模运算，常数更加优秀。
---
考虑在模意义下的多项式乘法，如果模数 $p$ 有原根，那么本质上原根和单位根是相同的。$g^{0,1,\cdots,p-2}$ 相当于 $\omega_{p-1}^{0,1,\cdots,p-2}$，都是互不相同且能达到效果的。并且易得 $\omega_{n}=g^{\frac{p-1}{n}}$，因为 $n$ 次方均为 $1$。也就是说当 $n\mid p-1$ 时，单位根才存在。
而在 DFT 和 IDFT 的过程中，要求长度是 $2$ 的幂次，也就是说要满足 $2^k\mid p-1$，这也就限制了一些模数不可取。
常见的模数有
- $998244353=7\times 17\times 2^{23}+1$，原根为 $3$，能处理长度 $\le 2^{23}$。
- $1004535809=479\times 2^{21}+1$，原根为 $3$。
- $469762049 = 7\times 2^{26}+1$，原根为 $3$。
在 IDFT 时，要求 $\omega_n^{-i}$，本质就是 $g^{-\frac{p-1}{n}}=(g^{-1})^{\frac{p-1}{n}}$。
```cpp
void ntt(int a[],int n,int ty){
for(int i=0,j=0;i<n;i++){
if(i<j)swap(a[i],a[j]);
for(int k=n>>1;(j^=k)<k;k>>=1);
}
g[0]=1;
for(int i=1;i<n;i<<=1){
int t=i*2;
int x=ksm(ty==1?3:(P+1)/3,(P-1)/t);
for(int j=t-2;j>=0;j-=2)g[j]=g[j>>1],g[j+1]=(ll)g[j]*x%P;
for(int j=0;j<n;j+=t){
for(int k=j;k<j+i;k++){
int t0=a[k],t1=(ll)a[k+i]*g[k-j]%P;
a[k]=mod(t0+t1);
a[k+i]=mod(t0+P-t1);
}
}
}
if(ty==-1){
int x=ksm(n,P-2);
for(int i=0;i<n;i++)a[i]=(ll)a[i]*x%P;
}
}
```
## 差卷积
在一些题目当中会出现 $c_i=\sum\limits_{j=i}^{n}a_{j}b_{j-i}$ 的形式，因为 $j-(j-i)=i$，所以叫做差卷积。
只需
$$
\begin{aligned}
c_i&=\sum_{j=i}^{n}a_jb_{j-i}\\
&=\sum_{j=0}^{n-i}a_{i+j}b_{j}\\
&=\sum_{j=0}^{n-i}a'_{n-i-j}b_{j}
\end{aligned}
$$
其中 $a'$ 为 $a$ 的翻转。发现相当于求 $a'$ 与 $b$ 的和卷积的第 $n-i$ 项。
## 任意模数卷积
值域 $\le 10^9$，$n\le 10^5$。
### 9 NTT
因为不保证模数是 NTT 模数，所以不能直接 NTT。
考虑到最后答案 $\le 10^{23}$，所以可以三模 NTT。即取三个 NTT 模数分别得出结果，然后 crt 即可求出真实答案，再进行对模数取模。
### 8 FFT
因为答案太大，会爆 double，考虑拆位。
将 $x$ 拆成 $\lfloor\frac{x}{s}\rfloor s+x\bmod s$，其中 $s=\sqrt{10^9}$。然后 $A_1$ 存前半部分，$A_2$ 存后半部分，那么 $A*B=(A_1s+A_2)*(B_1s+B_2)=A_1*B_1s^2+(A_1*B_2+A_2*B_1)s+A_2*B_2$。
这样卷积答案就只有 $10^{14}$ 量级，但是总共要跑 $8$ 次 FFT。
### 5 FFT
考虑构造复数多项式 $F=A_1+A_2\mathrm{i}$，$G=A_1-A_2\mathrm{i}$ 和 $H=B_1+B_2\mathrm{i}$。发现通过算出 $F*H$ 和 $G*H$ 可以解出 $8$ FFT 需要的值。总共只有 $5$ 次 FFT。
$O(\text{5 FFT})<O(\text{9 NTT})$。
4 FFT 去哪了呢？
# 快速沃尔什变换
解决 $c_i=\sum\limits_{j\oplus k=i}a_jb_k$，其中 $\oplus$ 是一个位运算，如 $\mathbf{xor}$、$\mathbf{or}$ 和 $\mathbf{and}$。
## 思想
### 异或卷积
考虑构造线性变换 $\mathcal{T}$ 使得 $\mathcal{T}(F)_i\cdot \mathcal{T}(G)_i=\mathcal{T}(H)_i$。相当于对于一个向量乘以矩阵 $M$ 得到一个新的向量。可以这样写
$$
\left(\sum_{j=0}^{n-1}a_jM_{j,i}\right)\left(\sum_{j=0}^{n-1}b_{j}M_{j,i}\right)=\sum_{j=0}^{n-1}c_jM_{j,i}
$$
可以变形得到
$$
\begin{aligned}
\sum_{j=0}^{n-1}\sum_{k=0}^{n-1}a_jb_kM_{j,i}M_{k,i}&=\sum_{j=0}^{n-1}\sum_{k\ \mathbf{xor}\ l=j}a_{k}b_{l}M_{j,i}\\
&=\sum_{j=0}^{n-1}\sum_{k=0}^{n-1}a_jb_kM_{j\ \mathbf{xor}\ k,i}
\end{aligned}
$$
于是可以得到 $M_{j,i}M_{k,i}=M_{j\ \mathbf{xor}\ k,i}$。然后经过一顿推导，发现
$$
F_i(k)=F_{i-1}(k)+F_{i-1}(k+2^{i-1})\\
F_i(k+2^{i-1})=F_{i-1}(k)-F_{i-1}(k+2^{i-1})
$$
然后就可以 $O(n\log n)$ 算了。逆变换相当于解二元一次方程，最后除以 $n$ 即可。
```cpp
void fwt(int a[],int n,int ty){
for(int i=1;i<n;i<<=1){
int t=i*2;
for(int j=0;j<n;j+=t){
for(int k=j;k<j+i;k++){
int t0=a[k],t1=a[k+i];
a[k]=mod(t0+t1);
a[k+i]=mod(t0+p-t1);
}
}
}
if(ty==-1){
int x=ksm(n,p-2);
for(int i=0;i<n;i++)a[i]=(ll)a[i]*x%p;
}
}
```
### 或卷积
同理可以得到 $M_{j,i}M_{k,i}=M_{j\ \mathbf{and}\ k,i}$，构造 $M_{x,y}=[x\subseteq y]$ 满足条件。
那么 $\mathcal{T}(F)_i=\sum\limits_{j\subseteq i}a_j$，相当于做高维前缀和，有
$$
F_i(k)=F_{i-1}(k)\\
F_i(k+2^{i-1})=F_{i-1}(k)+F_{i-1}(k+2^{i-1})
$$
逆变换易得。时间复杂度 $O(n\log n)$。
```cpp
void fwt(int a[],int n,int ty){
for(int i=1;i<n;i<<=1){
int t=i*2;
for(int j=0;j<n;j+=t){
for(int k=j;k<j+i;k++){
a[k+i]=mod(a[k]*ty+a[k+i]);
}
}
}
}
```
### 与卷积
对于与卷积，与或卷积类似，构造 $M_{x,y}=[y\subseteq x]$ 满足条件。然后即做高维后缀和，有
$$
F_i(k)=F_{i-1}(k)+F_{i-1}(k+2^{i-1})\\
F_i(k+2^{i-1})=F_{i-1}(k+2^{i-1})
$$
代码与或卷积类似。
## 子集卷积
求
$$
h(i)=\sum_{j\subseteq i}f(j)g(i\ \mathbf{xor}\ j)
$$
等价于
$$
h(i)=\sum_{j\ \mathbf{or}\ k=i,|j|+|k|=|i|}f(j)g(k)
$$
其中 $|x|$ 表示 $\text{popcount}(x)$。
把长度看成 $2^n$。如果只有第一个条件我们可以 $O(2^n n)$ 求出，考虑增加第二个条件。
发现有点类似和卷积的形式，新开一维状态表示 $1$ 的个数，有
$$
h_w(i)=\sum_{j\ \mathbf{or}\ k=i}\sum_{t=0}^{n}f_t(j)g_{w-t}(k)
$$
于是我们便可以枚举 $w=0\sim n$ 求出 $\mathcal{T}(F_{w})$ 和 $\mathcal{T}(G_{w})$，然后用 $O(2^{n}n^2)$ 的时间复杂度合并，然后再每一位逆变换回去即可。总体时间复杂度 $O(2^nn^2)
常见用于优化状压 dp。
多项式全家桶
https://www.luogu.com.cn/paste/1kz4jufv
多项式求逆
已知 A(x)，求 B(x) 使 A(x)B(x)\equiv 1\pmod {x^n}。
取模是因为 B 可能有无限多项，故只取前 n 项。O(n^2) 暴力是容易的，直接递推即可；考虑用倍增优化。
假设已经求出 B 的前 m 项，记为 B_m，考虑求出 B_{2m}。我们可以得到到两个式子：
A(x)B_{2m}(x)\equiv 1\pmod {x^{2m}}\\
B_{2m}(x)-B_{m}(x)\equiv 0\pmod {x^{m}}
将第二个式子平方可得
B_{2m}(x)^2-2B_{2m}(x)B_m(x)+B_m(x)^2\equiv 0\pmod {x^{2m}}
再同时乘上 A(x)，由第一个式子得
B_{2m}(x)-2B_{m}(x)+B_m(x)^2A(x)\equiv 0\pmod{x^{2m}}
通过移项可得
B_{2m}(x)\equiv 2B_m(x)-A(x)B_m(x)^2\pmod{x^{2m}}
每次计算时 A 只取前 2m 项，那么时间复杂度为 \sum\limits_{i=1}^{\log_2 n}i2^i=n\log n。
多项式除法
已知 A(x) 和 B(x)，求 D(x) 和 R(x) 使得
A(x)=B(x)D(x)+R(x)
假设 A(x) 为 n 次，B(x) 为 m 次，则要求 R(x) 低于 m 次。D(x) 显然是 n-m 次，令 R(x) 为 m-1 次（高位如果不存在用 0 补位）。
考虑将系数翻转，即 A'(x)=A(x^{-1})x^n，因为
A(x^{-1})=B(x^{-1})D(x^{-1})+R(x^{-1})
所以
A(x^{-1})x^n=B(x^{-1})x^{m}D(x^{-1})x^{n-m}+R(x^{-1})x^n\\
A'(x)=B'(x)D'(x)+R'(x)x^{n-m+1}
可以得到
A'(x)\equiv B'(x)D'(x)\pmod {x^{n-m+1}}
于是可以算出 D'(x)，且是 n-m 次的，符合要求。
然后 R(x)=A(x)-B(x)D(x)。总时间复杂度 O(n\log n)。
线性递推
已知 f_0\sim f_{m-1}，d_{0}\sim d_{m-1}，且
f_{n}=\sum_{i=0}^{m-1}f_{n-i-1}d_i
求 f_{n}。
通过矩阵快速幂可以做到 O(m^3\log n)，考虑优化。
有特征方程 x^{m}=\sum\limits_{i=0}^{m-1}x^{m-i-1}d_i。任何 f_{n} 都可以用 f_{0\sim m-1} 表示，所以只需求出 x^{n}。
考虑倍增，已知 x^{k}，如何求出 x^{2k}。
x^{k}=\sum_{i=0}^{m-1}a_ix^{i}
将两边同时平方，需要将 x^{m\sim 2m-2} 的项进行降幂处理。事实上，这个操作等价于将它对 (x^{m}-\sum\limits_{i=0}^{m-1}x^{m-i-1}d_i) 取模，因为该式结果为 0，而降幂本质上是不断减去这个多项式乘以一个整式，而结果不变，所以是等价的。
取模用多项式除法即可，于是总时间复杂度为 O(m\log m\log n)。
多项式多点求值
已知 F(x)，求 F(x_{1\sim m})。
先考虑在 m 较小的情况下，想到构造 G(x) 使得 G(x_i)=F(x_i)，且 G(x) 的次数不超过 m 次。
仿照线性递推降幂的思路，由于
A(x)=\prod_{i=1}^{m}(x-x_i)
对于 x=x_i 结果为 0，且次数为 m，所以将 F(x) 对 A(x) 取模后便能得到 G(x)。时间复杂度瓶颈在带值 O(m^2)。
考虑进行分治。令 A_{l,r}(x) 为 \prod\limits_{i=l}^{r}(x-x_i)，那么将 F(x) 从上至下一直取模到 A_{i,i}(x)，容易发现此时的 F(x) 就是 F(x_i)，因为 F(x) 只剩常数项，且取模不改变 F(x_i) 的值。时间复杂度为 O(n\log^2 n)。
转置原理优化
https://www.luogu.com.cn/article/5l5gpirn
多项式开方
已知 A(x)，求 B(x)^2\equiv A(x)\pmod {x^{n}}。
暴力可以 O(n^2) 做，仿照求逆的思路，考虑倍增。
假设已知 B_m(x)，求 B_{2m}(x)。可以得到
B_{2m}(x)-B_{m}(x)\equiv 0\pmod {x^{m}}\\
B_{2m}(x)^2-2B_{2m}(x)B_{m}(x)+B_m(x)^2\equiv 0\pmod {x^{2m}}\\
因为
B_{2m}(x)^2\equiv A(x)\pmod{x^{2m}}
所以
B_{2m}(x)\equiv\frac{A(x)+B_{m}(x)^2}{2B_{m}(x)}\pmod {x^{2m}}
多项式求导/积分
设 F(x)=\sum\limits_{i=0}^{n}a_ix^i，根据求导和积分法则可以得到
F'(x)=\frac{d}{dx}F(x)=\sum_{i=0}^{n-1}a_{i+1}(i+1)x^i\\
\int F(x)dx=\sum_{i=1}^{n+1}\frac{a_{i-1}}{i}x^i
时间复杂度显然为 O(n)。
多项式 ln（exp）
已知 A(x)，求 B(x)=\ln A(x)，即 e^{B(x)}=A(x)。
考虑先对 $\ln A(x)$ 求导，然后再积分回去
$$
\frac{d}{dx}\ln A(x)=\frac{d}{dA(x)}\ln A(x)\frac{d}{dx}A(x)=\frac{A'(x)}{A(x)}\\
\ln A(x)=\int\frac{A'(x)}{A(x)}dx
$$
总时间复杂度为 $O(n\log n)$。
---
已知 $A(x)$，求 $B(x)=e^{A(x)}$。
$B(x)$ 存在的条件为 $[x^0]A(x)=0$。
求解可以通过牛顿迭代。
>牛顿迭代指的是，已知 $G(x)$，求 $F(x)$ 使得 $G(F(x))\equiv 0\pmod {x^n}$。
>
>比如说想求 $B(x)\equiv \sqrt {A(x)}\pmod {x^n}$，则构造 $G(B(x))=A(x)-B^2(x)$，考虑如何求解。
>
>考虑倍增，已知 $B_m(x)$，如何求 $B_{2m}(x)$。
>
>我们让 $G(B_{2m}(x))$ 在 $B_{m}(x)$ 处泰勒展开
>$$
>G(B_{2m}(x))=\sum_{i=0}^{\infty}\frac{G^{(i)}(B_{m}(x))}{i!}(B_{2m}(x)-B_m(x))^i\equiv 0\pmod {x^{2m}}
>$$
>
>容易发现 $B_{2m}(x)-B_m(x)$ 的最低次至少 $x^{m}$，所以当 $i\ge 2$ 时在 $\bmod x^{2m}$ 后都变成 $0$ 了，所以又可以写成
>$$
>G(B_{2m}(x))\equiv G(B_{m}(x))+G'(B_{m}(x))(B_{2m}(x)-B_{m}(x))\pmod {x^{2m}}
>$$
>
>然后便得到
>$$
>B_{2m}(x)\equiv B_{m}(x)-\frac{G(B_m(x))}{G'(B_m(x))}\pmod {x^{2m}}
>$$
回归本题，只需构造 $G(B(x))=A(x)-\ln B(x)$ 即可。注意 $G$ 求导是关于 $B_m(x)$，所以 $A(x)$ 是常量。
$$
B_{2m}(x)\equiv B_{m}(x)-\frac{A(x)-\ln B_{m}(x)}{-\frac{1}{B_m(x)}}\equiv B_m(x)(1-\ln B_m(x) +A(x))\pmod {x^{2m}}
$$
时间复杂度 $O(n\log n)$。
## 多项式快速幂
已知 $A(x)$，求 $B(x)=A^k(x)\bmod x^{n}$。
考虑两边取对数，得 $\ln B(x)=k\ln A(x)$，算出 $\ln B(x)$ 后再 $\exp$ 回去，接下来讨论两个问题：
- $k$ 应该对 $p$ 取模。考虑 $e^{x}$ 的泰勒展开，将 $x$ 替换为 $k\ln A(x)$，发现 $x$ 并不在指数位，所以可以直接取模。
- $[x^0]A(x)\ne 1$。显然不能直接求 $\ln$，找到 $A(x)$ 的最低次有效项 $a_ix^i$，所以
$$
\begin{aligned}
B(x)&=\left(a_ix^i\frac{A(x)}{a_ix^i}\right)^k\\
&=a_i^kx^{ik}\exp\left(k\ln \frac{A(x)}{a_ix^i}\right)
\end{aligned}
$$
因为对 $x^{n}$ 取模，所以若 $ik$ 太大则结果为全 $0$ 多项式。
## 快速拉格朗日插值
我们要求
$$
f(x)=\sum_{i=0}^ny_i\prod_{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}
$$
观察到分母和 $i,j$ 都有关，尝试拿出来单独求
$$
g_i(x)=\prod_{j\ne i}(x-x_j)
$$
容易发现对于 $k\ne i$，$g_k(x_i)=0$，所以可以把所有的 $g_i(x)$ 加起来，不影响 $g_i(x_i)$ 的值。
$$
G(x)=\sum_{i=1}^{n}\prod_{j\ne i}\left(x-x_j\right)
$$
这里已经可以直接分治 NTT 求了，但是常数较大，注意到
$$
H(x)=\prod_{i=1}^{n}(x-x_i)\\
\begin{aligned}
H'(x)&=\sum_{i=1}^{n}(x-x_i)'\prod_{j\ne i}(x-x_j)\\
&=\sum_{i=1}^{n}\prod_{j\ne i}(x-x_j)\\
&=G(x)
\end{aligned}
$$
所以只需要分治 NTT 求出 $H(x)$ 后求导即得 $G(x)$，然后用多项式多点求值求出 $G(x_{1\sim n})$。
现在即求
$$
f(x)=\sum_{i=0}^{n}\frac{y_i}{G(x_i)}\prod_{j\ne i}(x-x_j)
$$
直接上分治 NTT，有
$$
f_{i,i}(x)=\frac{y_i}{G(x_i)}\\
f_{l,r}(x)=f_{l,mid}(x)H_{mid+1,r}+f_{mid+1,r}(x)H_{l,mid}
$$
总时间复杂度 $O(n\log ^2 n)$。
# 生成函数
解决无标号计数问题常用普通生成函数，有标号计数问题则是指数生成函数。
## 普通生成函数
对于数列 $a$，定义其普通生成函数（OGF）为
$$
F(x)=\sum_{n\ge 0}a_nx^{n}
$$
容易发现 $F(x)G(x)$ 为数列 $\sum_{i=0}^{n}a_ib_{n-i}$ 的 OGF。
常见封闭形式有：
- 等比数列
$$
\sum_{n\ge 0}p^nx^n=\frac{1}{1-px}
$$
- 组合数列
$$
\sum_{n\ge 0}\binom{m+n-1}{n-1}=\frac{1}{(1-x)^m}=\left(\sum_{n\ge 0}x^n\right)^{m}
$$
可用隔板法证明。
- 等差幂数列
$$
\sum_{n\ge 0}x^{nk}=\frac{1}{1-x^k}
$$
将形式幂级数转换为封闭形式的作用是方便乘除。
## 指数生成函数
对于数列 $a$，定义其指数生成函数（EGF）为
$$
\hat{F}(x)=\sum_{n\ge 0}a_n\frac{x^n}{n!}
$$
容易发现 $\hat{F}(x)\hat{G}(x)$ 为数列 $\sum_{i=0}^{n}\binom{n}{i}a_ib_{n-i}$ 的 EGF。有组合意义。
注意到
$$
\hat{F}(x)=\sum_{n\ge 0}\frac{x^n}{n!}=e^x
$$
因为 $e^x$ 的麦克劳林展开为它。
## EGF 与 EXP 的组合意义
设 $\hat{F}(x)$ 是 $f(x)$ 的 EGF。
如果 $\exp \hat{F}(x)=\hat{G}(x)$，则 $\hat{G}(x)$ 是
$$
g_n=\sum_{[n]=\{S_1,S_2,\cdots,S_k\}}\prod_{i=1}^{k}f_{|S_i|}
$$
或者说
$$
g_n=\sum_{i=1}^{n}\binom{n-1}{i-1}f_ig_{n-i}
$$
的 EGF。
---
代数角度推导：
$$
g_n=\sum_{i=1}^{n}\binom{n-1}{i-1}f_ig_{n-i}\\
\frac{g_n}{(n-1)!}=\sum_{i=1}^{n}\frac{f_i}{(i-1)!}\frac{g_{n-i}}{(n-i)!}\\
\hat{G}'(x)x=\hat{F}'(x)\hat{G}(x)x\\
\hat{F}(x)=\int \frac{\hat{G}'(x)}{\hat{G}(x)}dx=\ln \hat{G}(x)
$$
---
求 $n$ 个点的简单有标号无向连通图数量。以下均是有标号的前提下。
设 $f(i)$ 为 $i$ 个点的简单无向连通图数量，$g(i)$ 为 $i$ 个点的简单无向图数量。所以 $g(i)=2^{\binom{i}{2}}$，即看每条边选不选。令 $\hat{F}(x)$ 和 $\hat{G}(x)$ 分别为它们的指数生成函数。
看上面两个式子，$g(i)$ 显然可以拆分成若干个连通块，所以符合第一个式子；然后考虑暴力，枚举一个连通块，其他的点随便连。但是为了避免算重，固定枚举 $1$ 号点所在连通块的大小，并在剩下 $n-1$ 个点中选出在连通块内的点，于是便能得到第二个式子。通过第二个式子便能在 $O(n^2)$ 的时间复杂度内求得 $f(n)$。
如果用上面的结论，就有 $\exp \hat{F}(x)=\hat{G}(x)$，如果是对的就能直接通过取 $\ln$ 算出 $F(x)$，考虑一下它的**组合意义**：
$$
\exp \hat{F}(x)=\sum_{i=0}^{\infty}\frac{\hat{F}(x)^i}{i!}
$$
考虑 $[x^k]$。分子表示的是 $k$ 个点包含 $i$ 个有序连通块的方案数，那么除以 $i!$ 即表示 $k$ 个点包含 $i$ 个无序连通块的方案数，那么相加起来即表示 $k$ 个点包含若干无序连通块的方案数，也就正好是 $[x^k]\hat{G}(x)$。所以就能对应上。
## 应用
有 $m$ 种物品，体积为 $a_i$，数量为 $b_i$，对 $1\sim n$ 求填满容积为 $x$ 背包的方案数。
设 $f(i)$ 为容积为 $i$ 的方案数，$F(x)$ 为它的普通生成函数。容易发现
$$
\begin{aligned}
F(x)&=\prod_{i=1}^{m}\left(1+x^{a_i}+x^{2a_i}+\cdots+x^{a_ib_i}\right)\\
&=\prod_{i=1}^{m}\frac{1-x^{a_i(b_i+1)}}{1-x^{a_i}}
\end{aligned}
$$
要求 $F(x)$，考虑两边取 $\ln$，
$$
\ln F(x)=\sum_{i=1}^{m}\ln(1-x^{a_i(b_i+1)})-\ln(1-x^{a_i})
$$
考虑如何求 $\ln (1-x^{k})$，可以求导后积分回去，即
$$
\begin{aligned}
\int \frac{-kx^{k-1}}{1-x^{k}} dx&=\int (-kx^{k-1})(1+x^k+x^{2k}+\cdots)dx\\
&=\int -k(x^{k-1}+x^{2k-1}+\cdots)dx\\
&=-\sum_{j=1}^{\infty}\frac{x^{kj}}{j}
\end{aligned}
$$
所以
$$
\ln F(x)=\sum_{i=1}^{m}\sum_{j=1}^{\infty}\frac{-x^{a_i(b_i+1)j}+x^{a_ij}}{j}
$$
我们只关心 $[x^k]F(x)(k\le n)$，所以如果枚举 $i$ ，枚举量是 $O(n\log n)$ 级别的，求出系数后再 $\exp$ 回去即可。总时间复杂度 $O(n\log n)$。
---
有集合 $S$，已知 $f(1\sim n)$ 为 $i$ 表示成 $S$ 中元素之和的方案数，求 $S$。
仿照上题，$s_k=[k\in S]$，有
$$
\ln F(x)=-\sum_{k=1}^{\infty}s_k\ln (1-x^k)
$$
写成
$$
\sum_{i=1}^{n}f(i)x^i=\sum_{k=1}^{\infty}\sum_{i=1}^{\infty}\frac{x^{ik}}{i}s_k
$$
每次比较 $x_i$ 的系数然后更新即可。
---
已知 $a_{1\sim m}$，$f(k)=\sum a_i^k$，求 $f(1\sim n)$。
设 $s_k=\sum_{i=1}^{m}[a_i=k]$。
$$
\begin{aligned}
F(x)&=\sum_{n=0}^{\infty}\sum_{i=1}^{m}(a_i x)^n\\
&=\sum_{i=1}^{m}\frac{1}{1-a_ix}
\end{aligned}
$$
分式相加不太好处理，考虑同时积分
$$
\begin{aligned}
\int F(x)dx&=\sum_{i=1}^{m}\int \frac{1}{1-a_ix}dx\\
&=\sum_{i=1}^{m}-\frac{1}{a_i}\ln(1-a_ix)
\end{aligned}
$$
有个 $\frac{1}{a_i}$ 不太好搞，考虑把 $F(x)$ 写为
$$
\sum_{i=1}^{m}\left(\frac{a_ix}{1-a_ix}+1\right)=m+x\sum_{i=1}^{m}\frac{a_i}{1-a_ix}
$$
然后就变成
$$
\sum_{i=1}^{m}-\ln(1-a_ix)=-\ln\left(\prod_{i=1}^{m}(1-a_ix)\right)
$$
然后分治 NTT 即可。