📖 附录 G：数学算法基础
⏱ 预计阅读时间：50 分钟 | 难度：🟡 中等
前置条件
在学习本附录之前，请确保你已掌握：
- 基本循环与数组（第 2.2 章）
- 模运算的基本性质（第 3.1 章）
🎯 学习目标
学完本附录后，你将能够：
- 用埃氏筛和欧拉线性筛高效枚举质数
- 用快速幂在 O(log N) 内计算大幂次取模
- 用扩展欧几里得求逆元
- 理解区间 DP 的核心框架并解决石子合并类问题
- 运用矩阵快速幂加速线性递推
G.1 质数筛法
为什么需要筛法？
判断一个数是否是质数，朴素做法是枚举 [2, √N]，时间 O(√N)。
但若要一次性求出 [2, N] 内所有质数，朴素做法是 O(N√N)，当 N = 10^7 时太慢。
筛法 利用「质数的倍数是合数」这一事实，批量标记合数。
G.1.1 埃拉托斯特尼筛（埃氏筛）
核心思想：
从 2 开始，每发现一个质数，就把它的所有倍数标记为合数。
📄 从 2 开始，每发现一个质数，就把它的所有倍数标记为合数。
筛选 N = 20 的过程：
初始：2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
标记 2 的倍数（从 2*2=4 开始）：
删除 4, 6, 8, 10, 12, 14, 16, 18, 20
标记 3 的倍数（从 3*3=9 开始）：
删除 9, 15 （12, 18 已经被 2 删除）
4 已被删除，跳过
标记 5 的倍数（从 5*5=25 > 20，停止）
最终质数：2, 3, 5, 7, 11, 13, 17, 19
从 ii 而不是 2i 开始的原因：
2×i, 3×i, ..., (i-1)×i 这些倍数，已经被更小的质数筛过了（因为它们有更小的质因子）。
📄 2×i, 3×i, ..., (i-1)×i 这些倍数，已经被更小的质数筛过了（因为它们有更小的质因子）。
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e7 + 5;
bool is_prime[MAXN]; // is_prime[i] = true 表示 i 是质数
vector<int> primes; // 所有质数列表
// 埃氏筛：筛出 [2, n] 内的所有质数
// 时间复杂度：O(N log log N)
void sieve_eratosthenes(int n) {
fill(is_prime + 2, is_prime + n + 1, true); // 初始全是质数
for (int i = 2; (long long)i * i <= n; i++) {
if (is_prime[i]) {
// 从 i*i 开始标记（更小的倍数已经被筛过）
for (int j = i * i; j <= n; j += i)
is_prime[j] = false;
}
}
// 收集质数
for (int i = 2; i <= n; i++)
if (is_prime[i]) primes.push_back(i);
}
int main() {
sieve_eratosthenes(100);
cout << "100 以内的质数：";
for (int p : primes) cout << p << " ";
cout << endl;
cout << "质数个数：" << primes.size() << endl;
// 输出：2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
// 个数：25
return 0;
}
G.1.2 欧拉线性筛（最优筛法）
埃氏筛的问题：一个合数可能被多个质数重复标记（如 12 被 2 和 3 各标记一次）。
欧拉线性筛 保证每个合数只被其最小质因子筛一次，时间复杂度严格 O(N)。
核心规则： 对于每个 i，只用「i 乘以不超过 i 的最小质因子的那些质数」来筛合数。
📄 C++ 完整代码
const int MAXN = 1e7 + 5;
int min_prime[MAXN]; // min_prime[i] = i 的最小质因子
vector<int> primes;
// 欧拉线性筛：O(N)
void sieve_linear(int n) {
fill(min_prime, min_prime + n + 1, 0); // 0 表示还没被筛
for (int i = 2; i <= n; i++) {
if (!min_prime[i]) {
// i 没有被任何更小的数筛过 → i 是质数
min_prime[i] = i; // 质数的最小质因子是自身
primes.push_back(i);
}
// 用 i 去筛 i 的倍数
for (int p : primes) {
if ((long long)i * p > n) break;
min_prime[i * p] = p; // p 是 i*p 的最小质因子
if (i % p == 0) {
// p 是 i 的最小质因子
// 若继续用更大的质数 q 筛 i*q，
// i*q 的最小质因子仍是 p（因为 p | i | i*q）
// 所以 i*q 会被 (i/p)*p*q 中较小的部分筛到
// 此处终止，保证每个合数只被筛一次
break;
}
}
}
}
int main() {
sieve_linear(30);
cout << "30 以内的质数：";
for (int p : primes) cout << p << " ";
cout << endl;
// 利用最小质因子做质因数分解
auto factorize = [&](int x) {
cout << x << " = ";
while (x > 1) {
int p = min_prime[x];
int cnt = 0;
while (x % p == 0) { x /= p; cnt++; }
cout << p << "^" << cnt << " ";
}
cout << endl;
};
factorize(360); // 360 = 2^3 3^2 5^1
return 0;
}
线性筛的额外能力： 同时计算积性函数（欧拉函数 φ、莫比乌斯函数 μ、约数个数等），这在数论题中非常有用。
G.1.3 筛法对比
| 方法 | 时间复杂度 | 空间 | 优势 |
|---|---|---|---|
| 朴素判质数 | O(√N) 每个 | O(1) | 单个数判断 |
| 埃氏筛 | O(N log log N) | O(N) | 代码简单，实际快 |
| 欧拉线性筛 | O(N) | O(N) | 同时记录最小质因子 |
| Bitset 埃氏筛 | O(N log log N / 64) | O(N/8) | 超大 N（≥ 10^8）时更快 |
G.2 快速幂
问题
计算 $a^n \mod p$，其中 n 可能高达 $10^{18}$。
朴素方法（循环乘 n 次）需要 O(N) 次运算，根本不可行。
核心思想：二进制分解指数
将 n 写成二进制形式，利用递推：
$$a^n = \begin{cases} (a^{n/2})^2 & n \text{ 为偶数} \ (a^{(n-1)/2})^2 \times a & n \text{ 为奇数} \end{cases}$$
示例： 计算 $3^{13}$
📄 Code 完整代码
13 = 1101₂ = 8 + 4 + 1
3^13 = 3^8 × 3^4 × 3^1
计算过程（从低位到高位）：
base = 3, exp = 13, result = 1
exp & 1 = 1 → result = 1 × 3 = 3
base = 3 × 3 = 9，exp = 6
exp & 1 = 0 → result 不变（= 3）
base = 9 × 9 = 81，exp = 3
exp & 1 = 1 → result = 3 × 81 = 243
base = 81 × 81 = 6561，exp = 1
exp & 1 = 1 → result = 243 × 6561 = 1594323
base = ..., exp = 0（结束）
验证：3^13 = 1594323 ✓
完整实现
📄 查看代码：完整实现
// 快速幂：计算 base^exp % mod
// 时间复杂度：O(log exp)
// 空间复杂度：O(1)
long long fast_pow(long long base, long long exp, long long mod) {
long long result = 1;
base %= mod; // 先对底数取模
while (exp > 0) {
// 若当前位（最低位）为 1，乘入结果
if (exp & 1)
result = result * base % mod;
// 底数平方，指数右移一位
base = base * base % mod;
exp >>= 1;
}
return result;
}
// 使用示例
int main() {
long long MOD = 1e9 + 7;
cout << fast_pow(2, 10, MOD) << endl; // 1024
cout << fast_pow(3, 100, MOD) << endl; // 981350898
cout << fast_pow(2, 1e18, MOD) << endl; // 通过快速幂计算，不会超时
return 0;
}
模逆元（求 a 关于 mod 的逆元）
当 mod 是质数时，由费马小定理 $a^{p-1} \equiv 1 \pmod{p}$，得 $a^{-1} \equiv a^{p-2} \pmod{p}$。
📄 C++ 完整代码
// 求 a 关于质数 mod 的模逆元
long long mod_inv(long long a, long long mod) {
return fast_pow(a, mod - 2, mod);
}
// 应用：计算组合数 C(n, k) mod p
long long C(int n, int k, long long mod) {
if (k > n || k < 0) return 0;
// 预处理阶乘和逆阶乘
vector<long long> fact(n + 1), inv_fact(n + 1);
fact[0] = 1;
for (int i = 1; i <= n; i++) fact[i] = fact[i-1] * i % mod;
inv_fact[n] = mod_inv(fact[n], mod);
for (int i = n - 1; i >= 0; i--) inv_fact[i] = inv_fact[i+1] * (i+1) % mod;
return fact[n] % mod * inv_fact[k] % mod * inv_fact[n-k] % mod;
}
G.3 区间动态规划（区间 DP）
什么是区间 DP？
区间 DP 是一类特殊的 DP，状态定义在区间上：
$$dp[i][j] = \text{将区间 } [i, j] \text{ 处理后的最优值}$$
关键特征：
- 大区间的最优解由小区间合并得到
- 必须先算小区间，再算大区间
- 枚举分割点 k，将 [i,j] 分为 [i,k] 和 [k+1,j]
状态转移框架
dp[i][j] = min/max over all k in [i, j-1]:
dp[i][k] + dp[k+1][j] + cost(i, j)
遍历顺序：按区间长度
📄 查看代码：遍历顺序：按区间长度
// 区间 DP 标准模板
int n;
int dp[305][305]; // dp[i][j] = 区间 [i,j] 的最优值
// 初始化：单个元素
for (int i = 1; i <= n; i++)
dp[i][i] = 基础情况;
// 按区间长度从小到大
for (int len = 2; len <= n; len++) { // 区间长度
for (int i = 1; i + len - 1 <= n; i++) { // 左端点
int j = i + len - 1; // 右端点
dp[i][j] = INF; // 初始为极值
// 枚举分割点 k
for (int k = i; k < j; k++) {
dp[i][j] = min(dp[i][j], dp[i][k] + dp[k+1][j] + cost(i, j));
}
}
}
// 答案
return dp[1][n];
G.3.1 经典例题：石子合并
N 堆石子排成一列，每次合并相邻两堆，代价为合并后石子总数，求合并所有石子的最小（或最大）代价。
状态定义： dp[i][j] = 合并 [i, j] 堆石子的最优代价
转移： 枚举最后一次合并的分割点 k，即 [i,k] 先合并完，[k+1,j] 先合并完，再把这两堆合并
$$dp[i][j] = \min_{i \le k < j} {dp[i][k] + dp[k+1][j] + sum[j] - sum[i-1]}$$
其中 sum[j] - sum[i-1] 是区间 [i,j] 的石子总数（也是最后合并这步的代价）。
📄 其中 `sum[j] - sum[i-1]` 是区间 [i,j] 的石子总数（也是最后合并这步的代价）。
#include <bits/stdc++.h>
using namespace std;
int n;
int a[305];
int sum[305]; // 前缀和
int dp[305][305]; // dp[i][j] = 合并 [i,j] 的最小代价
int main() {
cin >> n;
for (int i = 1; i <= n; i++) {
cin >> a[i];
sum[i] = sum[i-1] + a[i];
}
// 单堆：代价 = 0（已经"合并完了"）
// dp[i][i] = 0（默认初始化）
// 区间 DP
for (int len = 2; len <= n; len++) {
for (int i = 1; i + len - 1 <= n; i++) {
int j = i + len - 1;
dp[i][j] = INT_MAX;
for (int k = i; k < j; k++) {
int cost = sum[j] - sum[i-1]; // 最后合并这步的代价
dp[i][j] = min(dp[i][j], dp[i][k] + dp[k+1][j] + cost);
}
}
}
cout << dp[1][n] << endl;
return 0;
}
追踪示例（a = [3, 5, 2, 1]）：
📄 Code 完整代码
前缀和：sum = [0, 3, 8, 10, 11]
初始：dp[1][1]=0, dp[2][2]=0, dp[3][3]=0, dp[4][4]=0
len=2：
dp[1][2] = dp[1][1] + dp[2][2] + sum[2]-sum[0] = 0+0+8 = 8
dp[2][3] = 0+0+7 = 7
dp[3][4] = 0+0+3 = 3
len=3：
dp[1][3]:
k=1: dp[1][1]+dp[2][3]+(sum[3]-sum[0]) = 0+7+10 = 17
k=2: dp[1][2]+dp[3][3]+(sum[3]-sum[0]) = 8+0+10 = 18
dp[1][3] = min(17,18) = 17
dp[2][4]:
k=2: dp[2][2]+dp[3][4]+(sum[4]-sum[1]) = 0+3+8 = 11
k=3: dp[2][3]+dp[4][4]+(sum[4]-sum[1]) = 7+0+8 = 15
dp[2][4] = 11
len=4：
dp[1][4]:
k=1: dp[1][1]+dp[2][4]+11 = 0+11+11 = 22
k=2: dp[1][2]+dp[3][4]+11 = 8+3+11 = 22
k=3: dp[1][3]+dp[4][4]+11 = 17+0+11 = 28
dp[1][4] = 22
答案：22
G.3.2 变形：括号匹配
给定一个括号序列，求添加最少多少个括号使其合法。
状态定义： dp[i][j] = 使 s[i..j] 合法需要添加的最少括号数
转移：
- 若 s[i] 和 s[j] 匹配（一对括号）：
dp[i][j] = dp[i+1][j-1] - 枚举分割点：
dp[i][j] = min(dp[i][k] + dp[k+1][j]) - 单个字符：
dp[i][i] = 1（必须添加一个配对括号）
📄 3. 单个字符：`dp[i][i] = 1`（必须添加一个配对括号）
#include <bits/stdc++.h>
using namespace std;
bool match(char a, char b) {
return (a == '(' && b == ')') ||
(a == '[' && b == ']') ||
(a == '{' && b == '}');
}
int main() {
string s;
cin >> s;
int n = s.size();
vector<vector<int>> dp(n, vector<int>(n, 0));
// 单个字符需要 1 个括号配对
for (int i = 0; i < n; i++) dp[i][i] = 1;
for (int len = 2; len <= n; len++) {
for (int i = 0; i + len - 1 < n; i++) {
int j = i + len - 1;
dp[i][j] = INT_MAX;
// 两端匹配
if (len >= 2 && match(s[i], s[j]))
dp[i][j] = min(dp[i][j], (len == 2 ? 0 : dp[i+1][j-1]));
// 枚举分割点
for (int k = i; k < j; k++)
dp[i][j] = min(dp[i][j], dp[i][k] + dp[k+1][j]);
}
}
cout << dp[0][n-1] << endl;
return 0;
}
G.4 矩阵快速幂（加速线性递推）
适用场景
如果一个递推关系的形式是：
$$\begin{pmatrix} f(n) \ f(n-1) \end{pmatrix} = M \times \begin{pmatrix} f(n-1) \ f(n-2) \end{pmatrix}$$
则可以用矩阵快速幂，在 O(K³ log N) 内计算第 N 项（K 为状态向量大小）。
以 Fibonacci 数列为例
$$f(n) = f(n-1) + f(n-2)$$
转化为矩阵形式：
$$\begin{pmatrix} f(n) \ f(n-1) \end{pmatrix} = \begin{pmatrix} 1 & 1 \ 1 & 0 \end{pmatrix} \times \begin{pmatrix} f(n-1) \ f(n-2) \end{pmatrix}$$
因此：
$$\begin{pmatrix} f(n) \ f(n-1) \end{pmatrix} = \begin{pmatrix} 1 & 1 \ 1 & 0 \end{pmatrix}^{n-1} \times \begin{pmatrix} f(1) \ f(0) \end{pmatrix}$$
📄 C++ 完整代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> Matrix;
const ll MOD = 1e9 + 7;
const int K = 2; // 状态维度
// 矩阵乘法
Matrix mat_mul(const Matrix& A, const Matrix& B) {
int n = A.size();
Matrix C(n, vector<ll>(n, 0));
for (int i = 0; i < n; i++)
for (int k = 0; k < n; k++)
for (int j = 0; j < n; j++)
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
return C;
}
// 矩阵快速幂：M^n
Matrix mat_pow(Matrix M, ll n) {
int sz = M.size();
// 初始化为单位矩阵
Matrix result(sz, vector<ll>(sz, 0));
for (int i = 0; i < sz; i++) result[i][i] = 1;
while (n > 0) {
if (n & 1) result = mat_mul(result, M);
M = mat_mul(M, M);
n >>= 1;
}
return result;
}
// 计算第 n 个 Fibonacci 数 mod MOD
ll fibonacci(ll n) {
if (n <= 1) return n;
// 转移矩阵
Matrix trans = {{1, 1}, {1, 0}};
Matrix result = mat_pow(trans, n - 1);
// 初始状态 [f(1), f(0)] = [1, 0]
// 结果 = result * [1, 0]^T 的第一个元素
return result[0][0]; // result[0][0]*f(1) + result[0][1]*f(0) = result[0][0]
}
int main() {
for (int i = 0; i <= 10; i++)
cout << "f(" << i << ") = " << fibonacci(i) << endl;
// 大数：f(10^18) mod (10^9+7)
cout << "f(10^18) = " << fibonacci(1e18) << endl;
return 0;
}
⚠️ 常见错误
| 错误 | 原因 | 修复方案 |
|---|---|---|
埃氏筛 i*i 溢出 | i 较大时 i*i 超 int 范围 | 用 (long long)i*i <= n 或 i <= n/i |
| 快速幂底数未取模 | 第一行漏了 base %= mod | 总是先 base %= mod |
| 区间 DP 初始化错误 | 忘记 dp[i][i] = 基础情况 | 单元素的基础情况必须手动设置 |
| 区间 DP 遍历顺序错误 | 没有按区间长度从小到大 | 最外层循环必须是 len |
| 矩阵乘法模运算位置 | 累加时溢出再取模 | 每次乘加后立即 % MOD |
💪 练习题（共 8 道，全部含完整解答）
🟢 基础练习（1~3）
题目 1：质因数分解
利用欧拉线性筛的 min_prime 数组，在 O(log N) 时间内对任意正整数做质因数分解，输出格式为 2^3 * 3^2 * 5。
✅ 完整解答
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e6 + 5;
int min_prime[MAXN];
vector<int> primes;
void sieve(int n) {
for (int i = 2; i <= n; i++) {
if (!min_prime[i]) { min_prime[i] = i; primes.push_back(i); }
for (int p : primes) {
if ((long long)i * p > n) break;
min_prime[i * p] = p;
if (i % p == 0) break;
}
}
}
void factorize(int x) {
cout << x << " = ";
bool first = true;
while (x > 1) {
int p = min_prime[x], cnt = 0;
while (x % p == 0) { x /= p; cnt++; }
if (!first) cout << " * ";
cout << p << "^" << cnt;
first = false;
}
cout << "\n";
}
int main() {
sieve(1e6);
factorize(360); // 2^3 * 3^2 * 5^1
factorize(1000000); // 2^6 * 5^6
factorize(97); // 97^1（质数）
return 0;
}
关键： min_prime[x] 是 x 的最小质因子，每次除尽后继续处理 x/p^k，循环 O(log x) 次。
题目 2：组合数取模
计算 C(N, K) mod (10^9+7)，N 和 K 可高达 10^6。
要求预处理阶乘和逆元，查询 O(1)。
✅ 完整解答
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 1e9 + 7;
const int MAXN = 1e6 + 5;
ll fact[MAXN], inv_fact[MAXN];
ll fast_pow(ll base, ll exp, ll mod) {
ll result = 1; base %= mod;
while (exp > 0) {
if (exp & 1) result = result * base % mod;
base = base * base % mod;
exp >>= 1;
}
return result;
}
void preprocess(int n) {
fact[0] = 1;
for (int i = 1; i <= n; i++) fact[i] = fact[i-1] * i % MOD;
inv_fact[n] = fast_pow(fact[n], MOD - 2, MOD); // 费马小定理
for (int i = n - 1; i >= 0; i--) inv_fact[i] = inv_fact[i+1] * (i+1) % MOD;
}
ll C(int n, int k) {
if (k < 0 || k > n) return 0;
return fact[n] % MOD * inv_fact[k] % MOD * inv_fact[n-k] % MOD;
}
int main() {
preprocess(1e6);
cout << C(10, 3) << "\n"; // 120
cout << C(1000000, 500000) << "\n"; // 大数取模
cout << C(5, 0) << "\n"; // 1
return 0;
}
预处理流程：
1. 正向递推 fact[0..N]：O(N)
2. 用快速幂求 fact[N] 的逆元：O(log MOD)
3. 反向递推所有逆阶乘：inv_fact[i] = inv_fact[i+1] * (i+1)，O(N)
4. 查询 C(n,k)：O(1)
题目 3：石子合并（区间 DP 基础）
N 堆石子排成一排，每次合并相邻两堆，代价为合并后的总石子数。求总代价最小值。
✅ 完整解答
#include <bits/stdc++.h>
using namespace std;
int main() {
int n; cin >> n;
vector<int> a(n + 1);
vector<int> sum(n + 1, 0);
for (int i = 1; i <= n; i++) { cin >> a[i]; sum[i] = sum[i-1] + a[i]; }
vector<vector<int>> dp(n + 1, vector<int>(n + 1, 0));
// len = 区间长度，从 2 开始（长度 1 代价为 0）
for (int len = 2; len <= n; len++) {
for (int i = 1; i + len - 1 <= n; i++) {
int j = i + len - 1;
dp[i][j] = INT_MAX;
for (int k = i; k < j; k++) {
int cost = sum[j] - sum[i-1]; // 合并 [i..j] 的代价
dp[i][j] = min(dp[i][j], dp[i][k] + dp[k+1][j] + cost);
}
}
}
cout << dp[1][n] << "\n";
return 0;
}
输入示例：
4
3 5 2 1
输出：22
追踪： 见 G.3.1 节的详细分步追踪。
🟡 进阶练习（4~6）
题目 4：矩阵链乘（经典区间 DP）
给 N 个矩阵，第 i 个的维度为 p[i-1] × p[i]（共 N+1 个维度值）。
求计算连乘 M1 × M2 × ... × MN 的最少乘法次数（通过改变括号顺序）。
提示： dp[i][j] = 计算第 i 到第 j 个矩阵乘积的最少乘法次数。
转移：枚举最后一次分割点 k，代价为 dp[i][k] + dp[k+1][j] + p[i-1]*p[k]*p[j]。
✅ 完整解答
#include <bits/stdc++.h>
using namespace std;
int main() {
int n; cin >> n;
vector<int> p(n + 1);
for (int& x : p) cin >> x;
// dp[i][j] = 计算矩阵 i..j 的最少乘法次数（1-indexed）
vector<vector<long long>> dp(n + 1, vector<long long>(n + 1, 0));
for (int len = 2; len <= n; len++) {
for (int i = 1; i + len - 1 <= n; i++) {
int j = i + len - 1;
dp[i][j] = LLONG_MAX;
for (int k = i; k < j; k++) {
long long cost = (long long)p[i-1] * p[k] * p[j];
dp[i][j] = min(dp[i][j], dp[i][k] + dp[k+1][j] + cost);
}
}
}
cout << dp[1][n] << "\n";
return 0;
}
示例（N=3，矩阵维度 10×30，30×5，5×60）：
p = [10, 30, 5, 60]
dp[1][2] = 10×30×5 = 1500
dp[2][3] = 30×5×60 = 9000
dp[1][3]:
k=1: dp[1][1] + dp[2][3] + 10×30×60 = 0 + 9000 + 18000 = 27000
k=2: dp[1][2] + dp[3][3] + 10×5×60 = 1500 + 0 + 3000 = 4500
dp[1][3] = 4500
输出：4500（先算 M2×M3，再算 M1×(M2M3)）
题目 5：Fibonacci 第 N 项（矩阵快速幂）
求斐波那契数列第 N 项 f(N) mod (10^9+7)，N 可高达 10^18。
✅ 完整解答
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> Matrix;
const ll MOD = 1e9 + 7;
Matrix mat_mul(const Matrix& A, const Matrix& B) {
int n = A.size();
Matrix C(n, vector<ll>(n, 0));
for (int i = 0; i < n; i++)
for (int k = 0; k < n; k++)
for (int j = 0; j < n; j++)
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
return C;
}
Matrix mat_pow(Matrix M, ll n) {
int sz = M.size();
Matrix result(sz, vector<ll>(sz, 0));
for (int i = 0; i < sz; i++) result[i][i] = 1; // 单位矩阵
while (n > 0) {
if (n & 1) result = mat_mul(result, M);
M = mat_mul(M, M);
n >>= 1;
}
return result;
}
ll fibonacci(ll n) {
if (n <= 1) return n;
// [f(n), f(n-1)] = [[1,1],[1,0]]^(n-1) * [f(1), f(0)]
Matrix trans = {{1, 1}, {1, 0}};
auto R = mat_pow(trans, n - 1);
// R * [1, 0]^T = R[0][0]*1 + R[0][1]*0 = R[0][0]
return R[0][0];
}
int main() {
cout << fibonacci(1) << "\n"; // 1
cout << fibonacci(10) << "\n"; // 55
cout << fibonacci(50) << "\n"; // 586268941
ll n = (ll)1e18;
cout << fibonacci(n) << "\n"; // 某个大数 mod 10^9+7
return 0;
}
为什么递推不行？ f(10^18) 需要 10^18 步，即便每步只是加法也需要 ~32 年（按 10^9 步/秒）。矩阵快速幂只需 ~120 次矩阵乘法。
题目 6：牛买卖（区间 DP 变形）
N 个人依次报价（整数），你可以买入或卖出，每次操作只涉及相邻时刻。
合并相邻两个操作（买-卖一对），代价为两个价格之差的绝对值。求将所有操作合并成一笔交易的最小总代价。
（即：括号匹配最少插入次数的数值版本）
提示： 这与石子合并类似，用区间 DP，dp[i][j] = 将第 i 到第 j 笔操作合并完毕的最小代价。
✅ 完整解答
本质是：对区间 [i,j] 枚举分割点 k，左右子区间各自合并后，再将两段合并，追加代价 |price[i] - price[j]|（合并时左端与右端的价差）。
#include <bits/stdc++.h>
using namespace std;
int main() {
int n; cin >> n;
vector<int> p(n + 1);
for (int i = 1; i <= n; i++) cin >> p[i];
vector<vector<int>> dp(n + 1, vector<int>(n + 1, 0));
for (int len = 2; len <= n; len++) {
for (int i = 1; i + len - 1 <= n; i++) {
int j = i + len - 1;
dp[i][j] = INT_MAX;
for (int k = i; k < j; k++) {
// 左边合并后最终价格是 p[k]（向左看），右边是 p[k+1]（向右看）
// 合并两段的额外代价：两段之间的「接口代价」= |p[i] - p[j]|
// （用石子合并的模型，cost = 区间端点差值）
dp[i][j] = min(dp[i][j],
dp[i][k] + dp[k+1][j] + abs(p[i] - p[j]));
}
}
}
cout << dp[1][n] << "\n";
return 0;
}
此题展示了区间 DP 的灵活性：只需修改「合并代价」的计算方式，即可处理不同的区间合并问题。
🔴 挑战练习（7~8）
题目 7：1000 以内所有质数之和
用埃氏筛求出 1000 以内的所有质数，输出它们的个数和总和。
再用线性筛改写，比较结果是否一致。
✅ 完整解答
#include <bits/stdc++.h>
using namespace std;
// 埃氏筛版本
void eratosthenes(int n) {
vector<bool> is_prime(n + 1, true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; (long long)i * i <= n; i++)
if (is_prime[i])
for (int j = i * i; j <= n; j += i)
is_prime[j] = false;
int cnt = 0; long long sum = 0;
for (int i = 2; i <= n; i++)
if (is_prime[i]) { cnt++; sum += i; }
cout << "埃氏筛：" << cnt << " 个质数，总和 = " << sum << "\n";
}
// 欧拉线性筛版本
void euler(int n) {
vector<int> min_p(n + 1, 0);
vector<int> primes;
for (int i = 2; i <= n; i++) {
if (!min_p[i]) { min_p[i] = i; primes.push_back(i); }
for (int p : primes) {
if ((long long)i * p > n) break;
min_p[i * p] = p;
if (i % p == 0) break;
}
}
long long sum = 0;
for (int p : primes) sum += p;
cout << "线性筛：" << primes.size() << " 个质数，总和 = " << sum << "\n";
}
int main() {
eratosthenes(1000);
euler(1000);
// 两者输出应完全一致：168 个质数，总和 = 76127
return 0;
}
题目 8：字符串加密（快速幂综合应用）
RSA 加密中，加密公式为：密文 = 明文^e mod n。
给定明文 M（一个整数）、指数 e、模数 n，计算密文。
其中 M, e, n 可高达 10^18。
✅ 完整解答
直接套用快速幂模板。注意当 M 和 n 都接近 10^18 时，M * M 会溢出 long long，需要用 __int128 或「龟速乘（二进制分解乘法）」。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 lll;
// 快速幂（支持 base 和 mod 都高达 10^18）
ll fast_pow(ll base, ll exp, ll mod) {
ll result = 1;
base %= mod;
while (exp > 0) {
if (exp & 1)
result = (lll)result * base % mod; // __int128 避免溢出
base = (lll)base * base % mod;
exp >>= 1;
}
return result;
}
// 龟速乘（不用 __int128 的替代方案）
ll safe_mul(ll a, ll b, ll mod) {
ll result = 0;
a %= mod;
while (b > 0) {
if (b & 1) result = (result + a) % mod;
a = (a + a) % mod;
b >>= 1;
}
return result;
}
ll fast_pow_safe(ll base, ll exp, ll mod) {
ll result = 1; base %= mod;
while (exp > 0) {
if (exp & 1) result = safe_mul(result, base, mod);
base = safe_mul(base, base, mod);
exp >>= 1;
}
return result;
}
int main() {
ll M, e, n;
cin >> M >> e >> n;
cout << fast_pow(M, e, n) << "\n"; // 密文
// 测试：M=2, e=10, n=1000 → 2^10=1024 mod 1000 = 24
cout << fast_pow(2, 10, 1000) << "\n"; // 24
return 0;
}
为什么用 __int128？
long long 最大约 9.2×10^18，但 base * base 可达 (10^18)^2 = 10^36，远超范围。
__int128 支持到约 1.7×10^38，足以处理中间计算。
💡 章节联系： 质数筛是数论题的基础工具；快速幂是「取模运算」场景（组合数、大幂次、RSA）的标配；区间 DP 是 USACO Gold 的高频题型，掌握了石子合并后还可以扩展到括号匹配、矩阵链乘、凸包剖分等。