模板索引：数学相关
基本数论概念
- 整除：对两个正整数 \(a, b \ (b \le a)\)，如果存在一个整数 \(k\)，使得 \(a = bk\)，则称 \(b\) 整除 \(a\)，记作 \(b | a\)。
- 算术基本定理：对任何一个大于 \(1\) 的整数 \(x\)，\(x\) 均能表示成若干个质因子的乘积。且这个表示法是唯一的。
- 因数：对一个正整数 \(x\)，所有能整除 \(x\) 的正整数均为 \(x\) 的因
数。 - 质数：因数只有 \(1\) 和自身的正整数叫质数；\(1\) 不是质数。
- 互质的定义：\(gcd(a, b) = 1\) .
- 因数分解：\(O(\sqrt{x})\) 对任意正整数 x 分解因子。
- 质因子分解：\(O(\sqrt{x})\) 对任意正整数 x 分解质因子。
- 欧几里得算法：\(gcd(a, b) = gcd(b, a \mod b)\)
线性筛
点击查看代码
void prime(){
for(int i = 2; i <= n; i++){
if(!vis[i]) p[++cnt] = i;
for(int j = 1; j <= cnt && 1ll * i * p[j] <= n; j++){
vis[i * p[j]] = 1;
if(i % p[j] == 0) break;
}
}
}
此外，如果预处理每次记下 \(n\) 的最小质因子，那么可以 \(O(log n)\) 完成一次质因子分解。
快速幂
ll work(ll x, ll y, ll m) {
if (y == 1) return x;
if (y == 0) return 1;
ll tmp = work(x, y / 2, m) % m;
if (y % 2 == 1) return tmp * tmp % m * x % m;
else return tmp * tmp % m;
}
二元一次不定方程
方程 \(ax + by = c\) 有解当且仅当 \(gcd(a, b) | c\).
利用拓展欧几里得算法可解 ax + by = gcd(a, b).
解得后两边同乘使右边变为c即可.
点击查看代码
long long exgcd(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
long long d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
P1516 青蛙的约会
大致题意：给一个长 \(L\) 的环，原点为 \(0\)，两只青蛙出发点坐标分别是 \(x，y\)，每次分别能跳 \(m，n\) 米，跳一次所花时间相同，都往正方向跳，问你最少跳几次能会面.
对于 \(100\%\) 的数据，\(1 \le x, y, m, n \le 2 \times 10^9\)，\(x \ne y\)，\(1 \le L \le 2.1 \times 10^9\)。
问题等价解非齐次方程 \(k(m - n) + zL = y - x\).
点击查看代码
ll x, y, m, n, L;
cin >> x >> y >> m >> n >> L;
ll k, z;
if(m < n){
swap(m, n);
swap(x, y);
}
ll gcd = exgcd(m - n, L, k, z); // (m - n)k + Lz = gcd(m - n, L);
if ((y - x) % gcd != 0)
{
cout << "Impossible" << endl;
return 0;
}
k *= (y - x) / gcd; // k' * (m - n) + L * z' = gcd(m - n, L) * (y - x) / gcd(m - n, L) = y - x
k = (k % (L / gcd) + L / gcd) % (L / gcd); // 找到最小的k
// k，z 的通解都是一个等差数列；
// k 的通解公差是 L / gcd(m - n, L)，所以 k_0 + t * (L / gcd(m - n, L)) 都是解
// 推导可以看oi-wiki
// 于是在[0，L / gcd)之间一定存在一个解
cout << k << endl;
欧拉函数
- 定义欧拉函数 \(φ(n)\) 表示小于 \(n\) 的数中与 \(n\) 互质的数的个数。（定义为小于等于不会影响，因为gcd(n, n） = n; 显然不影响计数。
- 特别的，\(φ(1) = 1\)
欧拉函数的性质
- 积性：如果 \(\gcd(a,b)=1\)，则 \(\varphi(a\times b)=\varphi(a)\times \varphi(b)\)。
- 欧拉反演：\(\sum_{d\mid n}\varphi(d)=n\)。
- 性质三：对任意质数 \(p\)，\(\varphi(p^k)=p^k-p^{k-1}\)。（证明：显然只有 \(p,2p,3p,\ldots,p^{k-1}p\) 与 \(p^k\) 不互质，共 \(p^{k-1}\) 个。）
计算欧拉函数
-
单个欧拉函数：设 \(n\) 的唯一分解式是 \(n=p_1^{k_1}p_2^{k_2}\cdots p_s^{k_s}\)。根据积性，\(\varphi(n)=\prod_{i=1}^{s}\varphi\left(p_i^{k_i}\right)\)。后者根据性质三可计算。
-
线性筛欧拉函数：在线性筛的同时可以筛出欧拉函数，设 \(p\) 是 \(i\) 的最小质因子，分三种情况讨论：
- \(i\) 为质数：\(\varphi(i)=i-1\). （i 自己不互质，其他都互质）
- \(p\) 和 \(\dfrac{i}{p}\) 互质：\(\varphi(i)=\varphi\left(\dfrac{i}{p}\right)\varphi(p)\) （利用积性）
- \(p\) 是 \(\dfrac{i}{p}\) 的质因子：\(\varphi(i)=\varphi\left(\dfrac{i}{p}\right)\times p\).
-
对第三条的证明：设 \(i=p^k q\)，其中 \(q\) 不含质因子 \(p\)。则
\[\varphi(i)=\varphi(p^k)\varphi(q)=(p^k-p^{k-1})\times \varphi(q)=p(p^{k-1}-p^{k-2})\times \varphi(q). \]逆用性质三，注意到
\[(p^{k-1}-p^{k-2})\times \varphi(q)=\varphi(p^{k-1})\varphi(q)=\varphi\left(\dfrac{i}{p}\right) \]故 $$ \varphi(i)=p \times \varphi\left(\dfrac{i}{p}\right) $$
线性筛欧拉函数
点击查看代码
phi[1] = 1;
for (int i = 2; i < N; i++) {
if (!vis[i]) {
p[++cnt] = i;
phi[i] = i - 1; // 情况一
}
for (int j = 1; j <= cnt && 1LL * i * p[j] < N; j++) {
int t = i * p[j];
vis[t] = 1;
phi[t] = phi[i] * (p[j] - 1); // 情况二
if (i % p[j] == 0) {
phi[t] = phi[i] * p[j]; // 情况三
break;
}
}
}
应用例题：P2158 [SDOI2008]仪仗队
先特判坐标轴与（2， 2），由于对称，只考虑上/下三角区域。
容易发现能被看到的人坐标应该满足 \(gcd(x, y) = 1\)。
暴力枚举每个点跑 \(gcd\) 复杂度都平方了，对于固定 \(x\)，发现满足条件的 \(y\) 的个数和意义都等于其欧拉函数值。
亦发现坐标轴也满足这个结论
逆元
- 定义
在模 \(p\) 同余下，对每个 \(x\)，能否找到一个整数 \(x^{-1}\in[1,p)\)，使任何数除以 \(x\) 等价于乘 \(x^{-1}\)？
这样的 \(x^{-1}\) 称为 \(x\) 在模 \(p\) 同余下的 逆元。
显然，逆元 \(x\) 应满足 \(x x^{-1}\equiv 1\ (\mathrm{mod}\ p)\)。
逆元的性质
-
逆元存在性定理
定理：\(x\) 在模 \(p\) 同余下存在逆元当且仅当 \(gcd(x, p) = 1\)。 -
逆元的唯一性
模 p 同余下，一个整数 x 的逆元若存在，则唯一。 -
逆元的单射性
在模质数 p 同余下，[1, p − 1] 内所有整数的逆元互不相同。 -
一个数的逆元的逆元等于它自身
利用拓展欧几里得算法求逆元
模板题：P1082同余方程
求关于 $ x$ 的同余方程 $ a x \equiv 1 \pmod {b}$ 的最小正整数解。
好处： \(p\) 可以不是质数
对任意给定的 \(x\) 和 \(p\)，显然 \(x \times x^{-1} \equiv 1 \pmod p\)。可以利用这个式子计算 \(x^{-1}\)。
上式等价于 \(xx^{-1}=kp+1 \Longleftrightarrow xx^{-1}+pk=1\)。
其中 \(x\) 和 \(p\) 是已知量，所以上式是一个不定方程，\(x^{-1}\) 和 \(k\) 是变量。用 exgcd 求解可以得到 \(x^{-1}\)。
根据不定方程有解的条件，\(\gcd(x,p)\mid 1 \Longleftrightarrow \gcd(x,p)=1\)。
这和逆元存在性定理相互印证。
点击查看代码
ll exgcd(ll a, ll b, ll &x, ll &y){ // ax + by = gcd(a, b)
if(b == 0){
x = 1; y = 0;
return a;
}
ll d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int main(){
ll a, b;
cin >> a >> b;
ll x, y;
ll d = exgcd(a, b, x, y);
cout << (x % b + b) % b << endl;
// 注意x可能为负数，输出时需要调整为正数
return 0;
}
欧拉定理及其退化（费马小定理）
暂略。
利用费马小定理求逆元
p必须是质数
根据费马小定理，\(a^{p-1}\equiv 1 \pmod p\)。
两侧同乘 \(a^{-1}\)，得到 \(a^{-1}\equiv a^{p-2}\pmod p\)。
于是快速幂求出 \(a^{p-2}\pmod p\) 即可。
P2613 【模板】有理数取余
给出一个有理数 \(c=\frac{a}{b}\)，求 \(c \bmod 19260817\) 的值。
这个值被定义为 \(bx\equiv a\pmod{19260817}\) 的解。
对于所有数据，保证 \(0\leq a \leq 10^{10001}\)，\(1 \leq b \leq 10^{10001}\)，且 \(a, b\) 不同时是 \(19260817\) 的倍数。
线性预处理逆元
- 方法一：顺推法，得到 \(i\) 在模 \(p\) 意义下的乘法逆元
\(p\) 应当为质数使得对所有 \(1 \le i \le n\) 都有逆元，即 \(gcd(i, p) = 1\).
设 \(inv_i\) 表示 \(i\) 的逆元，有递推式：
\(inv_i \equiv -\left\lfloor \dfrac{p}{i} \right\rfloor \times inv_{p \bmod i}\pmod p\)
边界为 \(inv_1=1\)。
点击查看代码
inv[1] = 1;
for (int i = 2; i <= n; i++)
inv[i] = p - (p / i) * inv[p % i] % p;
- 方法二： 逆推法，得到 \(fac[i]\) 的逆元
点击查看代码
fac[0] = 1;
for(int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i % p; // 先算模意义下阶乘
inv[n] = qpow(fac[n], p - 2, p);
for(int i = n - 1; i >= 1; i--) inv[i] = inv[i + 1] * (i + 1) % p; // 阶乘的逆元
for(int i = 1; i <= n; i++) cout << fac[i - 1] * inv[i] % p << endl;
中国剩余定理 CRT
要求 \(m_i\) 两两互质。
记 \(M=\prod m_i\)，\(M_i=M/m_i\)，\(t_i\) 是 \(M_i\) 在模 \(m_i\) 下的逆元。
则方程组的唯一通解是：
简单证明： 容易验证 \(x=\sum_{i=1}^{n} a_i t_i M_i\) 是方程的一组特解。
假设 \(x_1,x_2\) 均是方程组特解，显然 \(x_1-x_2 \equiv 0 \pmod{m_i}\)，
于是 \(M=\prod m_i \mid (x_1-x_2)\)。即两个解至少相差 \(M\)。
此外 \(M\equiv 0 \pmod{m_i}\)，于是对任何特解 \(x\)，\(x+kM \equiv x \equiv a_i \pmod{m_i}\)。
于是 \(x+kM\) 也是一个特解。由此，通解唯一性得证。
P1495 【模板】中国剩余定理（CRT)
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n;
ll m[20], M[20], MM = 1, a[20];
ll t[20];
ll inv(ll a, ll p){
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
// 规范在模p的界限内
}
int main(){
cin >> n;
for(int i = 1; i <= n; i++){
cin >> m[i] >> a[i];
MM *= m[i];
}
for(int i = 1; i <= n; i++){
M[i] = MM/m[i];
t[i] = inv(M[i], m[i]);
}
ll x = 0;
for(int i = 1; i <= n; i++){
x += (__int128)a[i] * t[i] % MM * M[i] % MM;
x %= MM;
}
cout << x << endl;
return 0;
}
附：防止爆 longlong：
- 方法一：利用long double
- 一般win环境，模数大于 \(1e9\)，小于等于 \(1e12\)左右。
- 常见的 GCC/Clang + x86 的 80-bit 扩展精度，long double 往往是 64 位有效精度（LDBL_MANT_DIG = 64）。
这种情况下这个技巧通常能撑到 p 接近 1e18 量级
点击查看代码
ll mul(ll a, ll b, ll p){
ll d = (ll)(a * (long double) b / p + 0.5);
ll c = a * b - d * p;
return c < 0 ? c + p : c;
}
线性代数相关
高斯消元（带增广矩阵）
int gaosi(int n, int m){
int curi = 1; // 表示枚举哪一行，也是下一个主元要放的行号，等于主元个数+1
for(int j = 1; j <= m && curi <= n; j++){ // 按列枚举
// 1. 选主元：找到当前列j中，从行curi到n绝对值最大的行t
int t = curi;
for(int i = curi + 1; i <= n; i++)
if(fabs(a[i][j]) > fabs(a[t][j])) // 找到这一列的最大非0元素，减少浮点误差
t = i;
// 2. 如果当前列全为0，跳过这一列
if(fabs(a[t][j]) < eps) continue; // 这一列没找到非零主元
// 3. 交换行t和行curi（从列j开始交换，前面的列已处理）
for (int k = j; k <= m + 1; k++) // 把非0元素所在行交换到当前行
swap(a[t][k], a[curi][k]);
// 4. 主元归一：将行curi的主元位置变为1（倒序避免主元被提前修改）
for(int k = m + 1; k >= j; k--) // 主元归一，其他的也要相对应除以主元位
a[curi][k] /= a[curi][j]; // 注意要倒着写，比较巧
// 5. 消去其他所有行的第j列（包括上面的行，直接得到行最简形）
for(int i = 1; i <= n; i++) // 用当前主元行 curi，把其他所有行 i 的第 j 列消成 0。
if(i != curi && fabs(a[i][j]) > eps ) // 只消去非零行
for(int k = m + 1; k >= j; k--) // 注意要倒着写，比较巧
a[i][k] -= a[curi][k] * a[i][j];
curi++; // 主元个数+1，下一个主元放在下一行
}
// 6. 检查是否存在矛盾行：0x + 0y + ... = 非0
for (int i = curi; i <= n; i++) {
if (fabs(a[i][m + 1]) > eps) {
return 0; // 无解
}
}
// 7. 根据主元个数判断解的情况
curi--; // ccuri 等于主元个数加1，因为它维护的是当前在哪行找新的主元
if (curi == m) {
return 1; // 主元个数=未知数个数，唯一解0
} else {
return 2; // 主元个数<未知数个数，无穷多解
}
}
容斥原理
3 元容斥公式
对于三个集合 \(A_1, A_2, A_3\)，有：
\(n\) 元容斥公式
对于 \(n\) 个集合 \(A_1, A_2, \dots, A_n\)，有：
也可以写成更紧凑的形式：
四元容斥的简要证明
如果把图拆成：\(E_i\)：恰好落在 i 个集合里（\(1 \le i \le 4\)）
那么：
但容斥各项对应的是：
代进去：
化简正好就是：
矩阵快速幂
警示
填初状态向量的时候不要填反了。
比如【矩阵快速幂优化线性递推式】，列向量的第一个应该是 \(f_{n-1}\)，最后一个才是 \(f_1\)。
矩阵快速幂优化线性递推式
对于由前面的项的线性组合得到的后面的项的公式就叫线性组合公式，比如斐波那契数列。
但其通项公式并不容易求出，或者比较复杂，此时如果我们需要一个相当大的项时，\(O(N)\) 的复杂度也无法求出。
我们引入矩阵乘计算线性递推公式
考虑对于以下线性递推公式构建系数矩阵：
第一行用于表示 \(a_n\) ，后面的行分别用于保留 \(a_i\)。
注意最后一行的倒数第二个位置才是 \(1\)，因为 $ a_{n - r} $ 会自然消失掉。
其他位置都是 \(0\).
例题P1962斐波那契数列
求斐波那契数列第 \(N\) 项。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7;
void mmul(ll a[2][2], ll b[2][2], ll c[2][2], ll x, ll y, ll z){
// a, b是系数矩阵，c是结果矩阵
// a的尺寸是x*y，b的尺寸是y*z，c的尺寸是x*z
ll d[2][2];
//一定要把结果装到临时变量里，再复制到c
// 不能直接赋值给c
// 否则当a和c是同一个数组时就错了
for(int i = 0; i < x; i++)
for(int j = 0; j < z; j++){ // 枚举c_{i, j}
d[i][j] = 0;
for(int k = 0; k < y; k++)
d[i][j] = (d[i][j] + a[i][k] * b[k][j]) % mod;
}
for(int i = 0; i < x; i++)
for(int j = 0; j < z; j++) // 注意是z
c[i][j] = d[i][j];
}
void mqpow(ll X[2][2], ll y){
ll Tmp[2][2] = {
{1, 0},
{0, 1}, // 单位矩阵
};
while(y){// 二进制拆分y次幂
if(y & 1) mmul(X, Tmp, Tmp, 2, 2, 2); // Tmp = X*Tmp
mmul(X, X, X, 2, 2, 2); // 注意这里是X自己迭代
y >>= 1;
}
for(int i = 0; i < 2; i++)
for(int j = 0; j < 2; j++)
X[i][j] = Tmp[i][j];
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
ll n;
cin >> n;
if(n <= 2) {
cout << 1 << endl;
return 0;
}
ll K[2][2] = { // 构造的系数矩阵
{1, 1},
{1, 0},
};
ll A[2][2]={
{1}, // 等价于{1, 0}，但A我们当作 2 * 1 矩阵来使用
{1}
};
mqpow(K, n - 2);
mmul(K, A, A, 2, 2, 1);
cout << A[0][0]<< endl;
return 0;
}
重载乘法写法（更推荐）
不重载的写法对乘法顺序，初始化下标问题，快速幂易用错对象等问题有极大优化和避免，符合人类习惯。
注意，数组传递会传自身指针，故可以修改实参，但我们自己写的结构体是值传递，所以要返回好结果数组。
点击查看代码
struct matrix{
ll a[2][2];
matrix() {memset(a, 0, sizeof(a));}
matrix(int x){ // 单位矩阵初始化
memset(a, 0, sizeof(a));
for(int i = 0; i < 2; i++) a[i][i] = 1;
}
matrix operator*(const matrix &b) const{
matrix res;
for(int i = 0; i < 2; i++)
for(int j = 0; j < 2; j++)
for(int k = 0; k < 2; k++)
(res.a[i][j] += a[i][k] * b.a[k][j]) %= mod;
return res;
}
} K, A;
matrix mqpow(matrix x, ll y){
matrix tmp = matrix(1);
while(y){
if(y & 1) tmp = x * tmp;
x = x * x;
y >>= 1;
}
return tmp;
}
矩阵快速幂+kmp 求构造不含模式串的文本串方案数
非常非常经典的一类题
题目描述
阿申准备报名参加 GT 考试，准考证号为 \(N\) 位数\(X_1,X_2…X_N\ (0\le X_i\le 9)\)，他不希望准考证号上出现不吉利的数字。
他的不吉利数字\(A_1,A_2,\cdots, A_M\ (0\le A_i\le 9)\) 有 \(M\) 位，不出现是指 \(X_1,X_2\cdots X_N\) 中没有一段恰好等于 \(A_1,A_2,\cdots ,A_M\)，\(A_1\) 和 \(X_1\) 可以为 \(0\)。输入格式
第一行输入 \(N,M,K\) 接下来一行输入 \(M\) 位的数。
输出格式
阿申想知道不出现不吉利数字的号码有多少种，输出模 \(K\) 取余的结果。
对于全部数据，\(N\leq10^9\)，\(M\leq 20\)，\(K\leq10000\)。
题目大意：给定模式串，求构造长为 \(n\) 的文本串不含模式串的方案数。
考虑 \(DP\)，设 \(f[i][j]\) 表示已经构造了长为 \(i\) 的序列，且当前后缀与模式串的最长匹配前缀长度恰好为 \(j\) 的方案数，于是有转移方程：
其中 \(g[k][j]\) 当前已匹配了不吉利数字的前 \(k\) 位时，新填入一个数字后，恰好变成匹配了前 \(j\) 位的数字有多少种（\(0...9\) 共 \(10\) 个数字中有几个满足条件）。
这里的 \(j\) 很可能小于 \(k\)，如果新填入一个字符匹配上了模式串的下一位，则 \(g[k][k + 1] + 1\)，如不然，则跳失配数组找最长匹配。
这里自然可以用 \(kmp\) 来解决！
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, m, mod;
struct matrix{
ll a[30][30];
matrix(){memset(a, 0, sizeof(a));}
matrix(int x){ // 单位矩阵初始化
memset(a, 0, sizeof(a));
for(int i = 0; i <= m - 1; i++) a[i][i] = 1;
}
matrix operator*(const matrix &b) const{
matrix res;
for(int i = 0; i <= m - 1; i++)
for(int j = 0; j <= m - 1; j++)
for(int k = 0; k <= m - 1; k++)
(res.a[i][j] += a[i][k] * b.a[k][j]) %= mod;
return res;
}
} G, K;
void pre(string s){
int n = s.length();
vector<int> pi = prefix(s); // 求kmp的pi数组，函数代码略
for(int i = 0; i <= m - 1; i++)
for(char num = '0'; num <= '9'; num++){
int j = i;
while(j && num != s[j]) j = pi[j - 1];
if(s[j] == num) j++;
if(j == m) continue; // 完整匹配了，不合法，跳过
G.a[i][j]++; // 从i跳到j的方案数 + 1
}
}
int main(){
cin >> n >> m >> mod;
string str;
cin >> str;
pre(str);
// 注意K只使用第一行，表示一个阶段（即文本串当前长度）多个状态的向量组合。
// 9 是指匹配长为 0 的有 9 种方案；1 是指匹配长为 1 的有一种方案。这就是初阶段，即 n = 1 的状态。
K.a[0][0] = 9; K.a[0][1] = 1;
G = mqpow(G, n - 1);
K = K * G;
ll ans = 0;
for(int i = 0; i <= m - 1; i++) (ans += K.a[0][i]) %= mod;
cout << ans << endl;
return 0;
}
矩阵快速幂优化图上递推
矩阵乘法的思想和 \(floyd\) 这类会有一定的重合。
P3758 可乐
有一个机器人，初始在 \(1\) 号城市上。这个机器人每一秒都会随机触发以下三种行为之一： 停在原地，去下一个相邻的城市，自爆。在第 \(0\) 秒时可乐机器人在 \(1\) 号城市，问在 \(N\) 个城市，\(M\) 条道路的图，经过 \(t\) 秒，可乐机器人的行为方案数是多少？
输入格式
第一行输入两个正整数 \(N\)，\(M\)。\(N\) 表示城市个数，\(M\) 表示道路个数。
接下来 \(M\) 行每行两个整数 \(u\)，\(v\)，表示 \(u\)，\(v\) 之间有一条道路。保证两座城市之间只有一条路相连，且没有任何一条道路连接两个相同的城市。
最后一行是一个整数 \(t\)，表示经过的时间。输入样例
3 2 1 2 2 3 2输出样例
8样例解释
- \(1\) ->爆炸。
- \(1\) -> \(1\) ->爆炸。
- \(1\) -> \(2\) ->爆炸。
- \(1\) -> \(1\) -> \(1\)。
- \(1\) -> \(1\) -> \(2\)。
- \(1\) -> \(2\) -> \(1\)。
- \(1\) -> \(2\) -> \(2\)。
- \(1\) -> \(2\) -> \(3\)。
共 \(8\) 种。数据范围与约定
- 对于 \(20\%\) 的数据，保证 \(t \leq 1000\)。
- 对于\(100\%\)的数据，保证 \(1 < t \leq 10^6\)，\(1 \leq N \leq30\)，\(0 < M < 100\)，\(1 \leq u, v \leq N\)。
第一步：从最基础的事实出发——邻接矩阵的幂就是路径计数
先暂时忘掉爆炸，考虑一个更纯粹的问题：给一张图，从 \(i\) 出发恰好走 \(k\) 步到 \(j\) 的走法有多少种？
先考虑朴素动态规划，令 \(dp_{i, j}\) 表示第 \(i\) 秒，走到第 \(j\) 个城市的方案数。
转移方程为 \(dp[s + 1][v] = \sum_{u = 0} ^ n {dp[s][u] * A[u][v]}\)。
定义“阶段”表示第 \(i\) 秒，“状态”表示现在在第 \(j\) 个城市。
把一个阶段 \(i\) 的所有状态写在一行成为一个行向量 \(DP_i\)，便容易知道 \(DP_t = DP_0 \times A^t\)，这就是以下矩阵快速幂思路的底层逻辑。
定义一个矩阵 \(A\)，其中 \(A[i][j]\) 表示“从 \(i\) 一步走到 \(j\) 的方法数”。注意这里是方法数，不是 0/1 可达性。
断言： \((A^k)[i][j]\) 就是从 \(i\) 出发恰好 \(k\) 步到 \(j\) 的走法数。
用归纳法证。当 \(k=1\) 时，
这就是定义。假设 \(k\) 步时成立，看 \(k+1\) 步：
这个求和公式在组合意义上是什么？一条从 \(i\) 走到 \(j\) 的 \(k+1\) 步路径，必然可以拆成“前 \(k\) 步从 \(i\) 走到某个中间点 \(m\)”加上“最后一步从 \(m\) 走到 \(j\)”。枚举这个中间点 \(m\)，每种选择下的方案数就是两段方案数的乘积（乘法原理），再把所有可能的 \(m\) 加起来（加法原理）——正好就是上面的式子。
所以矩阵乘法的定义公式
本身就是“在中间点处拆路径”这个组合原理的数学化身。 不是巧合，是设计如此。这也解释了为什么矩阵乘法在那么多计数问题里都能派上用场——它就是“乘法原理 + 加法原理”的代数封装。
第二步：如何把“三种行为”构造成转移矩阵
现在回到这道题。机器人在每一秒有三种选择，我们要把它们都写成“从某个状态转移到某个状态”。
- 停在原地： 相当于从 \(i\) 转移到 \(i\)。给矩阵加一个自环，即所有 \(A[i][i]=1\)。
- 走到相邻城市： 如果 \(i\) 和 \(j\) 有边，就 \(A[i][j]=1\)。
- 自爆： 这是最需要技巧的一种。爆炸后机器人不在任何城市了，怎么办？
引入一个虚拟的“已爆炸”状态，记为 \(0\) 号点。这样整个系统有 \(N+1\) 个状态：\(0\) 号是“已爆炸”，\(1\) 到 \(N\) 是正常的城市。
然后在矩阵里做两件事：
一件是所有正常城市都能一步转移到 \(0\) 号点，即 \(A[i][0]=1\) 对 \(i=1,\dots,N\) 成立——这代表“在这一秒选择了自爆”。
另一件更关键：给 \(0\) 号点也加一个自环，\(A[0][0]=1\)。
根据初始邻接关系和其他边造出的转移矩阵，就可以理解为 \(1\) 秒内从任意城市到任意城市的转移方案数。
第三步：答案统计
通过矩阵快速幂很容易得到 \(A^t\)，而答案就是 \(\sum_{v = 0}^{n}{A[1][v]}\)。
原因是，第 \(0\) 秒时动态规划的初状态可以写为 \(f_0 = [0 1 0 0 0 0 ... 0]\) （注意这里第下标 \(0\) 为 \(0\) 是我们建立的虚点），表示初始在节点 \(1\)，我们将初状态向量放在左边，左乘转移矩阵，得到末状态再求和就是答案了，而显然求和结果等于 \(A^t\)的第下标 \(1\) 行求和结果，那就不必要显示再计算了。
为什么这里要用行向量作状态左乘转移矩阵？
因为如果使用列向量作为状态并右乘转移矩阵，则构建转移矩阵的时候，加边 \(u -> v\) 的时候应该令 \(A[v][u] = 1\) 而不是 \(A[u][v] = 1\)，即对矩阵 \(A\) 做转置。
首先用行向量肯定左乘，列向量肯定右乘，不然不满足矩阵乘法维度要求。
我们单看一个状态的转移过程：\(dp_{s+1}[v] = \sum_{u = 0} ^ n {dp_s[u] * A[u][v]}\)，即转移到第 \(v\) 个城市，需【枚举起点 \(u\)】 随后乘 【\(u\) 到 \(v\) 的方案数】再求和。
如果你使用列向量，并且不转置转移矩阵，当列向量右乘转移矩阵的时候，转移矩阵第 \(i\) 行表示的是第 \(i\) 个城市到每一个城市的转移方案，左矩阵行乘右矩阵列时，相当于固定来源城市，得到的是其抵达每个城市的方案数量。
但是我们需要的是固定终点城市，枚举来源城市，恰好反了！
对小边权拆点实现0/1权重图
一定要注意拆点之后要把数组同倍数开大！
其实这个技术应该放在图论部分，但是我刚学到矩阵快速幂用这个方法，索性先写这里了。
P4159 迷路
如果边权只有 \(0, 1\)，那就是上一道题了。
但是这个题边权是 \(0-9\)。
我们拆点，把点 \(i\) 拆成 \(pos(i, 0-9)\) 十个点，其中 \(pos(i, 0)\) 是每个点的真实点，其他都是辅助点，有点类似分层图。
其他比赛时非系统学习的方法总结
值域前缀差
原出处：2025辽宁省赛
给定一个长度为n的正整数组,请您求出有多少个非空子数组*满足:该非空子数组的正整数之和能被出现在该非空子数组中的最大数字整除。
一个数组a是一个数组b的非空子数组,当且仅当a可以通过从b的开头删除零个或者多个数以及从结尾删除零个或者多个数而得到,并且a含有至少一个数。
数字(digit)是指构成数(number)的0,1,2,3,4,5,6,7,8,9。例如，出现在数组 [213] 中的数字有1,2,3,最大数字是3;出现在数组 [2025,11,15] 中的数字有0,1,2,5,最大数字是5。
原来我考虑单调栈，给每个区间找一个具体的代表最大值，当最大值出现多次时，区间可能被左边那个最大值“认领”，也可能被右边那个“认领”，就容易重复或遗漏。
而 GPT 给出的 F(x,k) 方法，完全放弃了“找一个代表”这件事，直接换了一个更抽象的角度来看区间：
令 F(x, k) 统计的是：最大值 ≤ k 且和能被 x 整除的区间个数。
那么 “最大值恰好等于 x” 的区间，不就是 “≤ x 但 不是 ≤ x−1” 吗，即直接 \(F(x,x) - F(x,x-1)\) 就是答案，从而绕过原方法要一直讨论区间，产生复杂难处理的容斥。
此外，假设固定模 \(x\) 时找答案，原来考虑设置余数为 \(1-9\) 的9个等价类，后来发现其实只需要记录以前区间余数为 \(x\) 的个数，这样以 y 为余数且能被 x 整除的区间，就是减掉前缀中余数也是 y 的，剩下的区间不就是整除的，记录个数即可，无需 DP 转移。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int t;
int n;
const int maxn = 1e5 + 10;
int a[maxn], b[maxn];
int turn(int x);
ll f(int x, int k){
int cnt[10];
memset(cnt, 0, sizeof(cnt));
int pre = 0;
ll ans = 0;
cnt[0] = 1; // 注意前缀可以和自己匹配的情况
for(int i = 1; i <= n; i++){
if(b[i] > k){
pre = 0;
memset(cnt, 0, sizeof(cnt));
cnt[0] = 1;
continue;
}
pre = (pre + a[i]) % x;
ans += cnt[pre];
cnt[pre]++;
}
return ans;
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> t;
while(t--){
cin >> n;
for(int i = 1; i <= n; i++) {
cin >> a[i];
b[i] = turn(a[i]);
}
ll ans = 0;
for(int k = 1; k <= 9; k++){
// 枚举以模 k 为 0 且最大值是 k 的区间
// 模 k 情况下，最大值为 k 的区间数减掉为k - 1的区间数，自然就是k的，绕过容斥
ans += f(k, k);
ans -= f(k, k - 1);
}
cout << ans << endl;
}
return 0;
}
// 找一个数的最大数字
int turn(int x){
int tmp = 0;
while(x){
tmp = max(tmp, x % 10);
x /= 10;
}
return tmp;
}
总结
核心思想： 统计“区间（连续子数组）中某个整体度量恰好为 \(x\)”的个数时，如果可以将 “度量 ≤ k” 轻松转化为 “禁止某些元素出现”，就能通过 值域前缀差 \(F(x) - F(x-1)\) 直接得到答案，而不需要为每个区间指定一个“代表元素”。
迁移示例
示例 1：统计“最大值恰好为 \(m\)”的子数组个数
原问题：给定数组 \(a\)，问有多少个子数组的最大值恰好等于 \(m\)（\(m\) 是整个数组的最大值之一）。经典做法用单调栈，但用前缀差可以这样做：
定义 \(f([l,r]) = \max(a_l,\dots,a_r)\)。
定义 \(G(k) =\) 最大值 \(\le k\) 的子数组个数。
那么最大值恰好为 \(m\) 的数量 = \(G(m) - G(m-1)\)。
计算 \(G(k)\)：把所有 \(a_i > k\) 的位置视为障碍，数组被切割成多个连续段，每段内所有元素都 \(\le k\)。在每一段内，任何子数组都满足最大值 \(\le k\)，因此该段贡献的子数组数就是 \(\binom{len}{2} + len\)。
累加所有段的子数组数即得 \(G(k)\)。
要点：完全不需要单调栈找边界，只需逐段计数。
示例 2：统计“不同整数的个数恰好为 \(K\)”的子数组
问题：给定数组，求有多少个子数组，其中不同数字的个数恰好为 \(K\)。
设 \(f([l,r]) =\) 区间内不同数字个数。
定义 \(H(k) =\) 不同数字个数 \(\le k\) 的子数组个数。
明显，\(H\) 可以通过滑动窗口在 \(O(nk)\) 内求出（但不可行）。这里更适合用另一种转换：恰好 \(K\) = 至多 \(K\) 个不同 - 至多 \(K-1\) 个不同。然后分别用双指针计算至多 \(k\) 个不同的子数组个数（经典滑动窗口，对于每个左端点找到最远的右端点）。
准确答案 = \(\text{at most } K - \text{at most } (K-1)\)。
这个思路在很多“恰好 \(K\) 个不同”问题里是标准解法，本质正是我们的“前缀差”。
示例 3：统计“最大公因数恰好为 \(d\)”的子序列个数（数论类）
问题：给定一个整数集合，问有多少个非空子序列，其最大公因数（GCD）恰好为 \(d\)。
设 \(f(S) = \gcd(S)\)。
定义 \(U(k) =\) 所有元素都是 \(k\) 的倍数的子序列个数。注意：如果子序列里所有数都是 \(k\) 的倍数，那么 \(\gcd\) 一定是 \(k\) 的某个倍数，即 \(\gcd\) 是 \(k\) 的倍数。这不是 \(\le k\)，而是“是 \(k\) 的倍数”。
我们转而统计 \(C(m) =\) 所有元素都是 \(m\) 的倍数的子序列个数 = \(2^{\text{cnt}[m]} - 1\)，其中 \(\text{cnt}[m]\) 是数组中能被 \(m\) 整除的数的个数。
那么 \(\gcd\) 恰好为 \(d\) 的个数，可以通过 倍数容斥/莫比乌斯反演 从 \(C\) 求得，这本质上也是基于“至少/倍数”的前缀和思想（但方向在倍数偏序上）。稍作变体：如果定义 \(G(k) =\) 所有元素都 \(\le k\) 且 gcd 是... 这里不太直接。不过经典做法是 \(ans[d] = C[d] - \sum_{j\ge 2} ans[d\cdot j]\)，这正是对倍数维度做差分。
如果你希望完全类比“不超过”的形式，可以设值域上偏序为整除，定义 \(F(x) =\) 所有元素都是 \(x\) 的约数的子序列个数（比较别扭）。但更常见的是用“倍数容斥”，它本质上和前缀差是一样的思想：将“精确 \(d\)”转化为“\(d\) 的倍数”与“更严格倍数”的差。
示例 4：统计“恰好包含某个特定子串”的字符串个数
问题：求长度为 \(n\) 的字符串中，恰好出现 \(k\) 次模式串 \(P\) 的个数。经典用 DP 或自动机，但也可用前缀差转化：
设 \(F(t) =\) 出现次数 \(\le t\) 的字符串个数（用 DP 加一维记录出现次数上限），则恰好 \(k\) = \(F(k) - F(k-1)\)。
这里虽然直接 DP 也可以求恰好，但用前缀差可以把“不超过”的 DP 状态简化（可以累加，不需要精确匹配次数的转移限制）。
浙公网安备 33010602011771号